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

For Mauricio: now, even in GENOTYPE_GIVEN_ALLELES mode, the VariantCallContext (which now inherits directly from VC) will report reference calls as confidently called if they pass the threshold even if the QUAL of the record itself is low because we were forced to have an ALT allele. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5485 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2011-03-22 02:42:28 +00:00
parent ab6a815184
commit 0ee687e49d
3 changed files with 25 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
View File

@ -202,14 +202,14 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
// the base was confidently callable
sum.nBasesCalledConfidently += value.confidentlyCalled ? 1 : 0;
// can't make a confident variant call here
if ( value.vc == null )
// can't make a call here
if ( !value.shouldEmit )
return sum;
try {
// we are actually making a call
sum.nCallsMade++;
writer.add(value.vc, value.refBase);
writer.add(value, value.refBase);
} catch (IllegalArgumentException e) {
throw new IllegalArgumentException(e.getMessage() + "; this is often caused by using the --assume_single_sample_reads argument with the wrong sample name");
}

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

@ -297,7 +297,7 @@ public class UnifiedGenotyperEngine {
// estimate our confidence in a reference call and return
if ( vc.getNSamples() == 0 )
return (UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES ?
estimateReferenceConfidence(stratifiedContexts, genotypePriors.getHeterozygosity(), false, 1.0) :
estimateReferenceConfidence(vc, stratifiedContexts, genotypePriors.getHeterozygosity(), false, 1.0) :
new VariantCallContext(generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext), refContext.getBase(), false));
// 'zero' out the AFs (so that we don't have to worry if not all samples have reads at this position)
@ -336,7 +336,7 @@ public class UnifiedGenotyperEngine {
if ( UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES && !passesEmitThreshold(phredScaledConfidence, bestAFguess) ) {
// technically, at this point our confidence in a reference call isn't accurately estimated
// because it didn't take into account samples with no data, so let's get a better estimate
return estimateReferenceConfidence(stratifiedContexts, genotypePriors.getHeterozygosity(), true, 1.0 - PofF);
return estimateReferenceConfidence(vc, stratifiedContexts, genotypePriors.getHeterozygosity(), true, 1.0 - PofF);
}
// create the genotypes
@ -423,7 +423,7 @@ public class UnifiedGenotyperEngine {
vcCall = variantContexts.iterator().next(); // we know the collection will always have exactly 1 element.
}
VariantCallContext call = new VariantCallContext(vcCall, passesCallThreshold(phredScaledConfidence));
VariantCallContext call = new VariantCallContext(vcCall, confidentlyCalled(phredScaledConfidence, PofF));
call.setRefBase(refContext.getBase());
return call;
}
@ -491,7 +491,7 @@ public class UnifiedGenotyperEngine {
AFs[i] = AlleleFrequencyCalculationModel.VALUE_NOT_CALCULATED;
}
private VariantCallContext estimateReferenceConfidence(Map<String, StratifiedAlignmentContext> contexts, double theta, boolean ignoreCoveredSamples, double initialPofRef) {
private VariantCallContext estimateReferenceConfidence(VariantContext vc, Map<String, StratifiedAlignmentContext> contexts, double theta, boolean ignoreCoveredSamples, double initialPofRef) {
if ( contexts == null )
return null;
@ -518,7 +518,7 @@ public class UnifiedGenotyperEngine {
P_of_ref *= 1.0 - (theta / 2.0) * MathUtils.binomialProbability(0, depth, 0.5);
}
return new VariantCallContext(QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING);
return new VariantCallContext(vc, QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING, false);
}
protected void printVerboseData(String pos, VariantContext vc, double PofF, double phredScaledConfidence, double[] normalizedPosteriors) {
@ -651,6 +651,11 @@ public class UnifiedGenotyperEngine {
return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING;
}
protected boolean confidentlyCalled(double conf, double PofF) {
return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING ||
(UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES && QualityUtils.phredScaleErrorRate(PofF) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING);
}
protected void computeAlleleFrequencyPriors(int N) {
// calculate the allele frequency priors for 1-N
double sum = 0.0;

21
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java
View File

@ -35,29 +35,32 @@ import org.broad.tribble.util.variantcontext.VariantContext;
*
* Useful helper class to communicate the results of calculateGenotype to framework
*/
public class VariantCallContext {
public VariantContext vc = null;
public class VariantCallContext extends VariantContext {
public byte refBase;
// Was the site called confidently, either reference or variant?
public boolean confidentlyCalled = false;
// Should this site be emitted?
public boolean shouldEmit = true;
VariantCallContext(VariantContext vc, boolean confidentlyCalledP) {
this.vc = vc;
super(vc);
this.confidentlyCalled = confidentlyCalledP;
}
VariantCallContext(VariantContext vc, boolean confidentlyCalledP, boolean shouldEmit) {
super(vc);
this.confidentlyCalled = confidentlyCalledP;
this.shouldEmit = shouldEmit;
}
VariantCallContext(VariantContext vc, byte ref, boolean confidentlyCalledP) {
this.vc = vc;
super(vc);
this.refBase = ref;
this.confidentlyCalled = confidentlyCalledP;
}
// blank variant context => we're a ref site
VariantCallContext(boolean confidentlyCalledP) {
this.confidentlyCalled = confidentlyCalledP;
}
public void setRefBase(byte ref) {
this.refBase = ref;
}