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

Joint Estimation model now emits a reasonable slod 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1969 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2009-11-03 21:12:42 +00:00
parent 11d950abe0
commit a545859c62
1 changed files with 39 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java
View File

@ -300,31 +300,7 @@ public class JointEstimateGenotypeCalculationModel extends GenotypeCalculationMo
double bestAFguess = findMaxEntry(alleleFrequencyPosteriors[indexOfMax]).second / (double)(frequencyEstimationPoints-1); double bestAFguess = findMaxEntry(alleleFrequencyPosteriors[indexOfMax]).second / (double)(frequencyEstimationPoints-1);
ArrayList<Genotype> calls = new ArrayList<Genotype>(); ArrayList<Genotype> calls = new ArrayList<Genotype>();
GenotypeLocusData locusdata = GenotypeWriterFactory.createSupportedGenotypeLocusData(OUTPUT_FORMAT, ref, loc); // first, populate the sample-specific data
if ( locusdata != null ) {
if ( locusdata instanceof ConfidenceBacked ) {
((ConfidenceBacked)locusdata).setConfidence(phredScaledConfidence);
}
if ( locusdata instanceof SLODBacked ) {
// TODO -- generate strand score
double strandScore = 0.0;
// calculate strand score
//EMOutput forward = runEM(ref, contexts, priors, StratifiedContext.FORWARD);
//EMOutput reverse = runEM(ref, contexts, priors, StratifiedContext.REVERSE);
//double forwardLod = (forward.getPofD() + reverse.getPofNull()) - overall.getPofNull();
//double reverseLod = (reverse.getPofD() + forward.getPofNull()) - overall.getPofNull();
//logger.debug("forward lod=" + forwardLod + ", reverse lod=" + reverseLod);
//double strandScore = Math.max(forwardLod - lod, reverseLod - lod);
//logger.debug(String.format("SLOD=%f", lod, strandScore));
((SLODBacked)locusdata).setSLOD(strandScore);
}
if ( locusdata instanceof AlleleFrequencyBacked ) {
((AlleleFrequencyBacked)locusdata).setAlleleFrequency(bestAFguess);
}
}
for ( String sample : GLs.keySet() ) { for ( String sample : GLs.keySet() ) {
// create the call // create the call
@ -347,6 +323,44 @@ public class JointEstimateGenotypeCalculationModel extends GenotypeCalculationMo
calls.add(call); calls.add(call);
} }
// next, the general locus data
// note that calculating strand bias involves overwriting data structures, so we do that last
GenotypeLocusData locusdata = GenotypeWriterFactory.createSupportedGenotypeLocusData(OUTPUT_FORMAT, ref, loc);
if ( locusdata != null ) {
if ( locusdata instanceof ConfidenceBacked ) {
((ConfidenceBacked)locusdata).setConfidence(phredScaledConfidence);
}
if ( locusdata instanceof AlleleFrequencyBacked ) {
((AlleleFrequencyBacked)locusdata).setAlleleFrequency(bestAFguess);
}
if ( locusdata instanceof SLODBacked ) {
// the overall lod
double overallLog10PofNull = Math.log10(alleleFrequencyPosteriors[indexOfMax][0]);
double overallLog10PofF = Math.log10(PofFs[indexOfMax]);
double lod = overallLog10PofF - overallLog10PofNull;
// the forward lod
initializeGenotypeLikelihoods(ref, contexts, StratifiedContext.FORWARD);
calculateAlleleFrequencyPosteriors(ref, loc);
double forwardLog10PofNull = Math.log10(alleleFrequencyPosteriors[indexOfMax][0]);
double forwardLog10PofF = Math.log10(PofFs[indexOfMax]);
// the reverse lod
initializeGenotypeLikelihoods(ref, contexts, StratifiedContext.REVERSE);
calculateAlleleFrequencyPosteriors(ref, loc);
double reverseLog10PofNull = Math.log10(alleleFrequencyPosteriors[indexOfMax][0]);
double reverseLog10PofF = Math.log10(PofFs[indexOfMax]);
double forwardLod = forwardLog10PofF + reverseLog10PofNull - overallLog10PofNull;
double reverseLod = reverseLog10PofF + forwardLog10PofNull - overallLog10PofNull;
//logger.debug("forward lod=" + forwardLod + ", reverse lod=" + reverseLod);
double strandScore = Math.max(forwardLod - lod, reverseLod - lod);
//logger.debug(String.format("SLOD=%f", strandScore));
((SLODBacked)locusdata).setSLOD(strandScore);
}
}
return new Pair<List<Genotype>, GenotypeLocusData>(calls, locusdata); return new Pair<List<Genotype>, GenotypeLocusData>(calls, locusdata);
} }
} }