window_average.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FUNCTION_WINDOW_AVERAGE_H_
#define GENESIS_POPULATION_FUNCTION_WINDOW_AVERAGE_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2024 Lucas Czech
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
    Contact:
    Lucas Czech <lucas.czech@sund.ku.dk>
    University of Copenhagen, Globe Institute, Section for GeoGenetics
    Oster Voldgade 5-7, 1350 Copenhagen K, Denmark
*/
 
#include "genesis/population/filter/filter_stats.hpp"
#include "genesis/population/filter/filter_status.hpp"
#include "genesis/population/filter/sample_counts_filter.hpp"
#include "genesis/population/filter/variant_filter.hpp"
#include "genesis/population/genome_locus_set.hpp"
#include "genesis/population/sample_counts.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/population/window/base_window.hpp"
#include "genesis/population/window/window_view.hpp"
 
#include <cassert>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Window Averaging
// =================================================================================================
 
enum class WindowAveragePolicy
{
    kWindowLength,
 
    kAvailableLoci,
 
    kValidLoci,
 
    kValidSnps,
 
    kSum,
 
    kProvidedLoci
};
 
template<class D>
inline size_t get_window_length( BaseWindow<D> const& window )
{
    // If the window is of type GenomeWindowStream, its total length is the sum of all lengths
    // of the chromosomes that the genome has covered.
    if( window.is_whole_genome() ) {
        size_t window_length = 0;
        for( auto const& chr_len : window.chromosomes() ) {
            window_length += chr_len.second;
        }
        return window_length;
    }
 
    // In all other cases, we simply use the window width function.
    return window.width();
}
 
template<class D>
inline size_t get_window_provided_loci_count(
    BaseWindow<D> const& window,
    std::shared_ptr<GenomeLocusSet> provided_loci
) {
    // We need a provided loci mask for this function.
    if( !provided_loci ) {
        throw std::invalid_argument(
            "Cannot comput window average denominator from provided loci mask, "
            "as no such mask was provided."
        );
    }
 
    // Helper function that takes a chromosome and positions on it, and returns the number
    // of set loci on the provided mask in that window.
    auto get_chr_loci_count_ = [&provided_loci]( std::string const& chr, size_t first, size_t last )
    {
        // Position checks. Should not happen if our internal usage is correct.
        if( first == 0 || last == 0 || first > last ) {
            throw std::invalid_argument(
                "Invalid positions first=" + std::to_string( first ) + " last=" +
                std::to_string(last) + " on chromosome \"" + chr +
                "\" for computing provided loci mask window denominator."
            );
        }
 
        // Get the chromosome. This can fail if the user did not provide a fitting mask.
        if( ! provided_loci->has_chromosome( chr )) {
            throw std::runtime_error(
                "Cannot compute provided loci on chromosome \"" +  chr +
                "\", as the provided loci mask does not contain the chromosome."
            );
        }
        auto const& bv = provided_loci->chromosome_positions( chr );
 
        // Mask check. In our internal usage, this should not fail, but we check anyway,
        // in case this function is called with a mask that is not meant for the given purpose.
        if( bv.get(0) ) {
            throw std::invalid_argument(
                "Invvalid provided loci mask with bit 0 set."
            );
        }
 
        // Another check based on user data. Can fail if the user did not provide a fitting mask.
        if( last >= bv.size() ) {
            throw std::runtime_error(
                "Cannot compute provided loci on chromosome \"" +  chr +
                "\", as the provided loci mask for the chromosome has length " +
                std::to_string( bv.size() - 1 ) + ", but the window covers positions " +
                std::to_string( first ) + "-" + std::to_string( last )
            );
        }
 
        // Finally, we have checked everything. Our first and last position are both inclusive,
        // while the bitvector count uses past-the-end, so we need to add one here for the last.
        return bv.count( first, last + 1 );
    };
 
    // If the window is a WindowStream over a whole genome, we use all its chromosomes.
    // This might not cover all chromosomes that the provided loci have data for,
    // in case a region filter was applied, so we want to account for that.
    // We are also rather strict in the process, to avoid accidental mismatches on the user side.
    // We count all positions between the beginning of the chromosome, and the end as either
    // specified in the data, or, if a sequence dict was given to the genome window stream, from that.
    if( window.is_whole_genome() ) {
        size_t loci_count = 0;
        for( auto const& chr_len : window.chromosomes() ) {
            loci_count += get_chr_loci_count_( chr_len.first, 1, chr_len.second );
        }
        return loci_count;
    }
 
    // Here we are in the normal case for all other window types. For those, we just simply use
    // their first and last positions. For ChromosomeWindowStream, this is similar as for the above
    // whole genome WindowStream, and is also be based on the sequence dict, if given, or on the data.
    return get_chr_loci_count_(
        window.chromosome(), window.first_position(), window.last_position()
    );
}
 
template<class D>
inline double window_average_denominator(
    WindowAveragePolicy policy,
    BaseWindow<D> const& window,
    std::shared_ptr<GenomeLocusSet> provided_loci,
    VariantFilterStats const& variant_filter_stats,
    SampleCountsFilterStats const& sample_counts_filter_stats
) {
    // We cannot have processed more samples than there were passing variants,
    // as we only should have processed a sample once its variant is found to be passing.
    // In case of the FST processor, we make sure that only one of the samples gets recorded
    // in the stats, so this works there as well. We do not simply assert this here,
    // as a misuse of the filters could result in an error here, which would be on the user
    // side, and so is an exception.
    // We skip this test when using the sum anyway, as in those cases, we might not have
    // the correct filter stats available in the first place.
    if(
        policy != WindowAveragePolicy::kSum &&
        sample_counts_filter_stats.sum() > variant_filter_stats[VariantFilterTag::kPassed]
    ) {
        throw std::invalid_argument(
            "Invalid stat counts with sample_counts_filter_stats.sum() > "
            "variant_filter_stats[VariantFilterTag::kPassed]"
        );
    }
 
    // Now select which value we want to return.
    switch( policy ) {
        case WindowAveragePolicy::kWindowLength: {
            return get_window_length( window );
        }
        case WindowAveragePolicy::kAvailableLoci: {
            auto const missing = variant_filter_stats_category_counts(
                variant_filter_stats, VariantFilterTagCategory::kMissingInvalid
            );
            return variant_filter_stats.sum() - missing;
        }
        case WindowAveragePolicy::kValidLoci: {
            // Here, we use the number of positions that passed all total variant filters
            // except for being a SNP, as well as the per-sample count of postions that
            // furthermore passed all sample positions.
            auto const valid_non_snps = variant_filter_stats_category_counts(
                variant_filter_stats, VariantFilterTagCategory::kInvariant
            );
            auto const valid_snps = sample_counts_filter_stats_category_counts(
                sample_counts_filter_stats, SampleCountsFilterTagCategory::kPassed
            );
            return valid_non_snps + valid_snps;
        }
        case WindowAveragePolicy::kValidSnps: {
            return sample_counts_filter_stats[SampleCountsFilterTag::kPassed];
        }
        case WindowAveragePolicy::kSum: {
            return 1.0;
        }
        case WindowAveragePolicy::kProvidedLoci: {
            return get_window_provided_loci_count( window, provided_loci );
        }
        default: {
            throw std::invalid_argument( "Invalid WindowAveragePolicy" );
        }
    }
    return std::numeric_limits<double>::quiet_NaN();
}
 
} // namespace population
} // namespace genesis
 
#endif // include guard