population/function/functions.hpp

Go to the documentation of this file.

#ifndef GENESIS_POPULATION_FUNCTION_FUNCTIONS_H_
#define GENESIS_POPULATION_FUNCTION_FUNCTIONS_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/sample_counts_filter.hpp"
#include "genesis/population/sample_counts.hpp"
#include "genesis/population/variant.hpp"
#include "genesis/sequence/reference_genome.hpp"
#include "genesis/utils/text/char.hpp"
 
#include <array>
#include <iosfwd>
#include <utility>
#include <string>
#include <vector>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Bases and Counts
// =================================================================================================
 
inline constexpr bool is_valid_base( char c )
{
    // Can't use a function call here, to comply with C++11 constexpr rules.
    // c = utils::to_upper( c );
    return (
        c == 'A' || c == 'a' ||
        c == 'C' || c == 'c' ||
        c == 'G' || c == 'g' ||
        c == 'T' || c == 't'
    );
}
 
inline constexpr bool is_valid_base_or_n( char c )
{
    // Can't use a function call here, to comply with C++11 constexpr rules.
    // c = utils::to_upper( c );
    return (
        c == 'A' || c == 'a' ||
        c == 'C' || c == 'c' ||
        c == 'G' || c == 'g' ||
        c == 'T' || c == 't' ||
        c == 'N' || c == 'n'
    );
}
 
SampleCounts::size_type get_base_count( SampleCounts const& sample, char base );
 
void set_base_count( SampleCounts& sample, char base, SampleCounts::size_type value );
 
// =================================================================================================
//     Sorting
// =================================================================================================
 
template<typename T>
std::array<size_t, 4> nucleotide_sorting_order( std::array<T, 4> const& values )
{
    // Sort quickly via sorting network, putting large values first.
    // See https://stackoverflow.com/a/25070688/4184258
    auto indices = std::array<size_t, 4>{{ 0, 1, 2, 3 }};
    if( values[indices[0]] < values[indices[1]] ) {
        std::swap( indices[0], indices[1] );
    }
    if( values[indices[2]] < values[indices[3]] ) {
        std::swap( indices[2], indices[3] );
    }
    if( values[indices[0]] < values[indices[2]] ) {
        std::swap( indices[0], indices[2] );
    }
    if( values[indices[1]] < values[indices[3]] ) {
        std::swap( indices[1], indices[3] );
    }
    if( values[indices[1]] < values[indices[2]] ) {
        std::swap( indices[1], indices[2] );
    }
 
    // Now they are sorted, largest ones first.
    assert( values[indices[0]] >= values[indices[1]] );
    assert( values[indices[1]] >= values[indices[2]] );
    assert( values[indices[2]] >= values[indices[3]] );
 
    return indices;
}
 
template<typename T>
std::array<size_t, 6> sample_counts_sorting_order( std::array<T, 6> const& v )
{
    // Implementation inspired by https://stackoverflow.com/a/2792216/4184258
    // We did not test if this is faster than a sorting network here. Fast enough for now anyway.
 
    // Obtain the index at which each value ends up in the sorting order, largest one first.
    int i0 = (v[0] <  v[1]) + (v[0] <  v[2]) + (v[0] <  v[3]) + (v[0] <  v[4]) + (v[0] <  v[5]);
    int i1 = (v[1] <= v[0]) + (v[1] <  v[2]) + (v[1] <  v[3]) + (v[1] <  v[4]) + (v[1] <  v[5]);
    int i2 = (v[2] <= v[0]) + (v[2] <= v[1]) + (v[2] <  v[3]) + (v[2] <  v[4]) + (v[2] <  v[5]);
    int i3 = (v[3] <= v[0]) + (v[3] <= v[1]) + (v[3] <= v[2]) + (v[3] <  v[4]) + (v[3] <  v[5]);
    int i4 = (v[4] <= v[0]) + (v[4] <= v[1]) + (v[4] <= v[2]) + (v[4] <= v[3]) + (v[4] <  v[5]);
    int i5 = (v[5] <= v[0]) + (v[5] <= v[1]) + (v[5] <= v[2]) + (v[5] <= v[3]) + (v[5] <= v[4]);
    assert( i0 + i1 + i2 + i3 + i4 + i5 == 15 );
 
    // At those indices in the result, set the position that they need to point to.
    std::array<size_t, 6> order;
    order[i0] = 0;
    order[i1] = 1;
    order[i2] = 2;
    order[i3] = 3;
    order[i4] = 4;
    order[i5] = 5;
 
    // Now everything is sorted.
    assert( v[order[0]] >= v[order[1]] );
    assert( v[order[1]] >= v[order[2]] );
    assert( v[order[2]] >= v[order[3]] );
    assert( v[order[3]] >= v[order[4]] );
    assert( v[order[4]] >= v[order[5]] );
    return order;
}
 
SortedSampleCounts sorted_sample_counts( SampleCounts const& sample );
 
std::pair<SortedSampleCounts, SortedSampleCounts> sorted_average_sample_counts(
    SampleCounts const& sample_a,
    SampleCounts const& sample_b
);
 
SortedSampleCounts sorted_sample_counts(
    Variant const& variant, bool reference_first, SampleCountsFilterPolicy filter_policy
);
 
// =================================================================================================
//     Allele Count
// =================================================================================================
 
size_t allele_count( SampleCounts const& sample );
 
size_t allele_count( SampleCounts const& sample, size_t min_count );
 
size_t allele_count( SampleCounts const& sample, size_t min_count, size_t max_count );
 
// =================================================================================================
//     Merging
// =================================================================================================
 
void merge_inplace( SampleCounts& p1, SampleCounts const& p2 );
 
SampleCounts merge( SampleCounts const& p1, SampleCounts const& p2 );
 
SampleCounts merge( std::vector<SampleCounts> const& p, SampleCountsFilterPolicy filter_policy );
 
inline SampleCounts merge_sample_counts( Variant const& v, SampleCountsFilterPolicy filter_policy )
{
    return merge( v.samples, filter_policy );
}
 
inline constexpr size_t nucleotide_sum( SampleCounts const& sample )
{
    return sample.a_count + sample.c_count + sample.g_count + sample.t_count;
}
 
inline size_t total_nucleotide_sum( Variant const& variant, SampleCountsFilterPolicy filter_policy )
{
    return nucleotide_sum( merge_sample_counts( variant, filter_policy ));
}
 
inline constexpr size_t sample_counts_sum( SampleCounts const& sample )
{
    return
        sample.a_count +
        sample.c_count +
        sample.g_count +
        sample.t_count +
        sample.n_count +
        sample.d_count
    ;
}
 
inline size_t total_sample_counts_sum( Variant const& variant, SampleCountsFilterPolicy filter_policy )
{
    return sample_counts_sum( merge_sample_counts( variant, filter_policy ));
}
 
// =================================================================================================
//     Consensus
// =================================================================================================
 
std::pair<char, double> consensus( SampleCounts const& sample );
 
std::pair<char, double> consensus( SampleCounts const& sample, bool is_covered );
 
char guess_reference_base(
    Variant const& variant,
    bool force = false,
    SampleCountsFilterPolicy filter_policy = SampleCountsFilterPolicy::kOnlyPassing
);
 
char guess_alternative_base(
    Variant const& variant,
    bool force = false,
    SampleCountsFilterPolicy filter_policy = SampleCountsFilterPolicy::kOnlyPassing
);
 
void guess_and_set_ref_and_alt_bases(
    Variant& variant,
    bool force = false,
    SampleCountsFilterPolicy filter_policy = SampleCountsFilterPolicy::kOnlyPassing
);
 
void guess_and_set_ref_and_alt_bases(
    Variant& variant,
    char ref_base,
    bool force = false,
    SampleCountsFilterPolicy filter_policy = SampleCountsFilterPolicy::kOnlyPassing
);
 
void guess_and_set_ref_and_alt_bases(
    Variant& variant,
    genesis::sequence::ReferenceGenome const& ref_genome,
    bool force = false,
    SampleCountsFilterPolicy filter_policy = SampleCountsFilterPolicy::kOnlyPassing
);
 
// =================================================================================================
//     Output
// =================================================================================================
 
std::ostream& operator<<( std::ostream& os, SampleCounts const& bs );
 
} // namespace population
} // namespace genesis
 
#endif // include guard