src/c++/lib/genotyping/BreakpointGenotyper.cpp
/**
 * return genotype likelihood given genotype (using Poisson model with internal parameters)
 * @param lambda Poisson distribution parameter
 * @param gv Genotype vector
 * @param read_counts Read count vector for each allele
 */
return genotype likelihood given genotype (using Poisson model with internal parameters)
@param lambda Poisson distribution parameter
@param gv Genotype vector
@param read_counts Read count vector for each allele
return genotype likelihood given genotype ( using Poisson model with internal parameters ) @param lambda Poisson distribution parameter @param gv Genotype vector @param read_counts Read count vector for each allele
double BreakpointGenotyper::genotypeLikelihood(
    double lambda, const GenotypeVector& gv, const vector<int32_t>& read_counts) const
{
    double log_phi;
    auto it = genotype_prior_.find(gv);
    if (it == genotype_prior_.end())
    {
        log_phi = 0;
    }
    else
    {
        log_phi = it->second;
    }
    vector<int> allele_ploidy;
    allele_ploidy.resize(n_alleles_, 0);
    for (unsigned int al = 0; al < n_alleles_; ++al)
    {
        for (const auto g : gv)
        {
            if (al == g)
            {
                ++allele_ploidy[al];
            }
        }
    }
    double gl = log_phi;
    for (unsigned int al = 0; al < n_alleles_; ++al)
    {
        if (allele_ploidy[al] == 0)
        {
            const double eps = (allele_error_rate_.size() == 1 ? allele_error_rate_[0] : allele_error_rate_[al]);
            const poisson_distribution<> error_distribution(lambda * eps);
            gl += log(pdf(error_distribution, read_counts[al]));
        }
        else
        {
            const double mu
                = (haplotype_read_fraction_.size() == 1 ? haplotype_read_fraction_[0] : haplotype_read_fraction_[al]);
            const double mu_with_ploidy = allele_ploidy[al] * mu;
            const poisson_distribution<> allele_count_distribution(lambda * mu_with_ploidy);
            gl += log(pdf(allele_count_distribution, read_counts[al]));
        }
        if (std::isinf(gl))
        {
            return -std::numeric_limits<double>::max();
        }
    }
    return gl;
}
double BreakpointGenotyper :: genotypeLikelihood ( double lambda , const GenotypeVector & gv , const vector < int32_t > & read_counts ) const { double log_phi ; auto it = genotype_prior_ . find ( gv ) ; if ( it == genotype_prior_ . end ( ) ) { log_phi = 0 ; } else { log_phi = it -> second ; } vector < int > allele_ploidy ; allele_ploidy . resize ( n_alleles_ , 0 ) ; for ( unsigned int al = 0 ; al < n_alleles_ ; ++ al ) { for ( const auto g : gv ) { if ( al == g ) { ++ allele_ploidy [ al ] ; } } } double gl = log_phi ; for ( unsigned int al = 0 ; al < n_alleles_ ; ++ al ) { if ( allele_ploidy [ al ] == 0 ) { const double eps = ( allele_error_rate_ . size ( ) == 1 ? allele_error_rate_ [ 0 ] : allele_error_rate_ [ al ] ) ; const poisson_distribution < > error_distribution ( lambda * eps ) ; gl += log ( pdf ( error_distribution , read_counts [ al ] ) ) ; } else { const double mu = ( haplotype_read_fraction_ . size ( ) == 1 ? haplotype_read_fraction_ [ 0 ] : haplotype_read_fraction_ [ al ] ) ; const double mu_with_ploidy = allele_ploidy [ al ] * mu ; const poisson_distribution < > allele_count_distribution ( lambda * mu_with_ploidy ) ; gl += log ( pdf ( allele_count_distribution , read_counts [ al ] ) ) ; } if ( std :: isinf ( gl ) ) { return - std :: numeric_limits < double > :: max ( ) ; } } return gl ; }
return genotype likelihood given genotype (using Poisson model with internal parameters)
 @param lambda Poisson distribution parameter
 @param gv Genotype vector
 @param read_counts Read count vector for each allele
return genotype likelihood given genotype ( using Poisson model with internal parameters ) @param lambda Poisson distribution parameter @param gv Genotype vector @param read_counts Read count vector for each allele
// compute how many copies of each allele we have in this GT // compute GL by summing all allele contributions // no copies -> all reads supporting this allele will be errors