//if ( stratifiedContexts.size() == 0 )
// return null;
if ( !vc.isBiallelic() )
return null;
final GenotypesContext genotypes = vc.getGenotypes();
if ( !vc.hasGenotypes() )
return null;
double ratioHom = 0.0;
double ratioHet = 0.0;
double weightHom = 0.0;
double weightHet = 0.0;
double overallNonDiploid = 0.0;
for ( Genotype genotype : genotypes ) {
if ( vc.isSNP() ) {
final int[] counts = getCounts(genotype, stratifiedContexts, vc);
// If AD was not calculated, we can't continue
if(counts == null)
continue;
final int n_allele = counts.length;
int count_sum = 0;
for(int i=0; i<n_allele; i++){
count_sum += counts[i];
}
double pTrue = 1.0 - Math.pow(10.0,-genotype.getGQ() / (double) 10 );
if ( genotype.isHet() ) {
final int otherCount = count_sum - (counts[0] + counts[1]);
// sanity check
if ( counts[0] + counts[1] == 0 )
continue;
// weight the allele balance by genotype quality so that e.g. mis-called homs don't affect the ratio too much
ratioHet += pTrue * ((double)counts[0] / (double)(counts[0] + counts[1]));
weightHet += pTrue;
overallNonDiploid += ( (double) otherCount )/((double) count_sum*genotypes.size());
} else if ( genotype.isHom() ) {
final int alleleIdx = genotype.isHomRef() ? 0 : 1 ;
final int alleleCount = counts[alleleIdx];
int bestOtherCount = 0;
for(int i=0; i<n_allele; i++){
if( i == alleleIdx )
continue;
if( counts[i] > bestOtherCount )
bestOtherCount = counts[i];
}
final int otherCount = count_sum - alleleCount;
ratioHom += pTrue*( (double) alleleCount)/((double) (alleleCount+bestOtherCount));
weightHom += pTrue;
overallNonDiploid += ((double ) otherCount)/((double) count_sum*genotypes.size());
}
// Allele Balance for indels was not being computed correctly (since there was no allele matching). Instead of
// prolonging the life of imperfect code, I've decided to delete it. If someone else wants to try again from
// scratch, be my guest - but make sure it's done correctly! [EB]
}