zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Put in a sanity check that MLEAC <= AN

This commit is contained in:
Eric Banks 2012-08-17 11:49:53 -04:00
parent 0a706c9105
commit 2676b7fc2e
1 changed files with 10 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
View File

@ -38,6 +38,7 @@ import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.baq.BAQ;
import org.broadinstitute.sting.utils.classloader.PluginManager;
import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
@ -442,10 +443,17 @@ public class UnifiedGenotyperEngine {
// add the MLE AC and AF annotations
if ( alleleCountsofMLE.size() > 0 ) {
attributes.put(VCFConstants.MLE_ALLELE_COUNT_KEY, alleleCountsofMLE);
final double AN = (double)builder.make().getCalledChrCount();
final int AN = builder.make().getCalledChrCount();
// let's sanity check that we don't have an invalid MLE value in there
for ( int MLEAC : alleleCountsofMLE ) {
if ( MLEAC > AN )
throw new ReviewedStingException(String.format("MLEAC value (%d) is larger than AN (%d) at position %s:%d", MLEAC, AN, loc.getContig(), loc.getStart()));
}
final ArrayList<Double> MLEfrequencies = new ArrayList<Double>(alleleCountsofMLE.size());
for ( int AC : alleleCountsofMLE )
MLEfrequencies.add((double)AC / AN);
MLEfrequencies.add((double)AC / (double)AN);
attributes.put(VCFConstants.MLE_ALLELE_FREQUENCY_KEY, MLEfrequencies);
}