Thanks to Guillermo I was able to isolate an instance of where the MLEAC > AN. It turns out that this is valid, e.g. when PLs are all 0s for a sample we no-call it but it's allowed to factor into the MLE (since that's the contract with the exact model). Removing the check in UG and instead protecting for it in the AlleleCount stratification.

This commit is contained in:
Eric Banks 2012-08-20 01:16:23 -04:00
parent 7fa76f719b
commit 97b191f578
2 changed files with 5 additions and 10 deletions

View File

@ -444,16 +444,10 @@ public class UnifiedGenotyperEngine {
if ( alleleCountsofMLE.size() > 0 ) {
attributes.put(VCFConstants.MLE_ALLELE_COUNT_KEY, alleleCountsofMLE);
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());
// the MLEAC is allowed to be larger than the AN (e.g. in the case of all PLs being 0, the GT is ./. but the exact model may arbitrarily choose an AC>1)
for ( int AC : alleleCountsofMLE )
MLEfrequencies.add((double)AC / (double)AN);
MLEfrequencies.add(Math.min(1.0, (double)AC / (double)AN));
attributes.put(VCFConstants.MLE_ALLELE_FREQUENCY_KEY, MLEfrequencies);
}

View File

@ -46,7 +46,8 @@ public class AlleleCount extends VariantStratifier {
int AC = 0; // by default, the site is considered monomorphic
if ( eval.hasAttribute(VCFConstants.MLE_ALLELE_COUNT_KEY) && eval.isBiallelic() ) {
AC = eval.getAttributeAsInt(VCFConstants.MLE_ALLELE_COUNT_KEY, 0);
// the MLEAC is allowed to be larger than the AN (e.g. in the case of all PLs being 0, the GT is ./. but the exact model may arbitrarily choose an AC>1)
AC = Math.min(eval.getAttributeAsInt(VCFConstants.MLE_ALLELE_COUNT_KEY, 0), nchrom);
} else if ( eval.hasAttribute(VCFConstants.ALLELE_COUNT_KEY) && eval.isBiallelic() ) {
AC = eval.getAttributeAsInt(VCFConstants.ALLELE_COUNT_KEY, 0);
} else if ( eval.isVariant() ) {
@ -56,7 +57,7 @@ public class AlleleCount extends VariantStratifier {
// make sure that the AC isn't invalid
if ( AC > nchrom )
throw new UserException.MalformedVCF(String.format("The AC or MLEAC value (%d) at position %s:%d " +
throw new UserException.MalformedVCF(String.format("The AC value (%d) at position %s:%d " +
"is larger than the number of chromosomes over all samples (%d)", AC,
eval.getChr(), eval.getStart(), nchrom));