afs_estimate.cpp

Go to the documentation of this file.

/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2021 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 <lczech@carnegiescience.edu>
    Department of Plant Biology, Carnegie Institution For Science
    260 Panama Street, Stanford, CA 94305, USA
*/
 
#include "genesis/population/functions/afs_estimate.hpp"
 
#include "genesis/sequence/functions/quality.hpp"
 
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <stdexcept>
 
namespace genesis {
namespace population {
 
// =================================================================================================
//     Allele Frequency Spectrum Estimate
// =================================================================================================
 
AfsPileupRecord convert_to_afs_pileup_record( SimplePileupReader::Record const& record )
{
    // Prepare return value.
    AfsPileupRecord result;
 
    // We expect the pileup file to contain exactly one sample. No mpileup here.
    if( record.samples.size() != 1 ) {
        throw std::runtime_error(
            "Invalid pileup record line that does not contain exactly one sample."
        );
    }
 
    // Get the sample stored in the record and check it.
    auto const& smp = record.samples[0];
    if( smp.read_bases.size() != smp.phred_scores.size() ) {
        throw std::runtime_error(
            "Invalid pileup record line where number of read bases "
            "differs from number of phred scores."
        );
    }
 
    // We might reserve too much here, if not all bases are either ancestral or derived.
    // That should not matter all too much in terms of memory - but is better for speed,
    // so let's do it this way.
    auto const base_cnt = smp.read_bases.size();
    result.alleles.reserve( base_cnt );
    result.phred_scores.reserve( base_cnt );
 
    // Helper function to check that the ancestral/reference base is valid.
    auto check_ancestral_base_char_ = []( char r ){
        if( r != 'A' && r != 'C' && r != 'G' && r != 'T' && r != 'N' ) {
            throw std::runtime_error(
                "Invalid ancestral base in pileup record line that is not in [ACGTN]."
            );
        }
    };
 
    // Get the derived base, depending on what our input is.
    char derived_base;
    if( smp.ancestral_base != '\0' ) {
        check_ancestral_base_char_( smp.ancestral_base );
        derived_base = smp.ancestral_base;
    } else if( record.reference_base != '\0' ) {
        check_ancestral_base_char_( record.reference_base );
        derived_base = record.reference_base;
 
    // TODO: else folded: simply use the most common as anc, and second most common (if present)
    // as derived allele. could use the base counts code here again, if we made AfsPileupRecord
    // derived from BaseCounts.
 
    } else {
        throw std::runtime_error(
            "Cannot convert pileup record line for allele frequency estimation, "
            "as it misses the ancestral allele information."
        );
    }
 
    // Convert bases to just bool values of whether they are ancestral or derived,
    // and skip all others.
    for( size_t i = 0; i < base_cnt; ++i ) {
        if( smp.read_bases[i] == record.reference_base ) {
            result.alleles.emplace_back( false );
            result.phred_scores.emplace_back( smp.phred_scores[i] );
        } else if( smp.read_bases[i] == derived_base ) {
            result.alleles.emplace_back( true );
            result.phred_scores.emplace_back( smp.phred_scores[i] );
        }
    }
    assert( result.alleles.size() == result.phred_scores.size() );
 
    return result;
}
 
std::vector<double> prob_cond_true_freq(
    size_t n, // n: number of haplotypes in the pool
    std::vector<bool> const& alleles, // vo: vector of alleles (1 if derived, 0 otherwise), size r
    std::vector<unsigned char> const& phred_scores, // SE: vector of error probabilities, size r
    bool unfolded // unfolded: 1 if unfolded, 0 else
) {
    // result p: conditional prob of vo given z,r,SE,unfolded. size n+1*1
    if( unfolded ) {
        return prob_cond_true_freq_unfolded( n, alleles, phred_scores );
    } else {
        // Compute for alleles and their inverse (ancestral and derived flipped), then average.
        auto const prob1 = prob_cond_true_freq_unfolded( n, alleles, phred_scores, false );
        auto const prob2 = prob_cond_true_freq_unfolded( n, alleles, phred_scores, true );
        assert( prob1.size() == n + 1 );
        assert( prob2.size() == n + 1 );
 
        auto result = std::vector<double>( n + 1 );
        for( size_t i = 0; i < n + 1; ++i ) {
            result[i] = ( prob1[i] + prob2[i] ) / 2.0;
        }
        return result;
    }
}
 
std::vector<double> prob_cond_true_freq_unfolded(
    size_t n, // n: number of haplotypes in the pool
    std::vector<bool> const& alleles, // vo: vector of alleles (1 if derived, 0 otherwise), size r
    std::vector<unsigned char> const& phred_scores, // SE: vector of error probabilities, size r
    bool invert_alleles
) {
    // Prepare result p: conditional prob of vo given z,r,SE,unfolded. size n+1*1
    // auto result = std::vector<double>( n + 1, 0.0 );
    auto result = std::vector<double>( n + 1, 1.0 );
 
    // r: number of reads at the position
    // nd: n converted to double
    auto const r = alleles.size();
    auto const nd = static_cast<double>( n );
    if( alleles.size() != phred_scores.size() ) {
        throw std::invalid_argument(
            "Invalid alleles and their error probabilities with different size."
        );
    }
 
    // Process all reads/alleles with their respective error probs
    for( size_t i = 0; i < r; ++i ) {
        // We only need the actual error probability here, so instead of allocating a vector for
        // all of them, we just convert from phread to err prob here on the fly.
        auto const err_prob = sequence::phred_score_to_error_probability( phred_scores[i] );
        assert( 0.0 < err_prob && err_prob <= 1.0 );
 
        // Equations (1) and (2) from Boitard et al 2012, 10.1093/molbev/mss090,
        // where we sum over all possible states for n chromosomes, with different
        // sums for whether we have the derived or ancestral allele.
        if( alleles[i] ^ invert_alleles ) {
            for( size_t j = 0; j < n + 1; ++j ) {
                auto const jn = static_cast<double>( j ) / nd;
                auto const term1 = ( 1.0 - err_prob ) * jn;
                auto const term2 = err_prob * ( 1.0 - jn );
                // result[j] += std::log( term1 + term2 );
                result[j] *= ( term1 + term2 );
            }
        } else {
            for( size_t j = 0; j < n + 1; ++j ) {
                auto const jn = static_cast<double>( j ) / nd;
                auto const term1 = ( 1.0 - err_prob ) * ( 1.0 - jn );
                auto const term2 = err_prob * jn;
                // result[j] += std::log( term1 + term2 );
                result[j] *= ( term1 + term2 );
            }
        }
    }
 
    // for( size_t j = 0; j < n + 1; ++j ) {
    //     result[j] = std::exp( result[j] );
    // }
    return result;
}
 
} // namespace population
} // namespace genesis