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

Merge pull request #185 from broadinstitute/gda_poolcaller_fix_47921867 

Corner case fix to General Ploidy SNP likelihood model.
This commit is contained in:
Mark DePristo 2013-04-24 04:14:05 -07:00
parent e83d9bef59 2ab270cf3f
commit df6ba74395
1 changed files with 23 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidySNPGenotypeLikelihoods.java
View File

@ -227,7 +227,7 @@ public class GeneralPloidySNPGenotypeLikelihoods extends GeneralPloidyGenotypeLi
* @param capBaseQualsAtMappingQual Cap base at mapping qual * @param capBaseQualsAtMappingQual Cap base at mapping qual
* @param minBaseQual Minimum base quality to consider * @param minBaseQual Minimum base quality to consider
* @param errorModel Site error model * @param errorModel Site error model
* @return Number of bases added * @return Number of bases added - only good bases actually added to GLs are counted.
*/ */
private int add(ReadBackedPileup pileup, boolean ignoreBadBases, boolean capBaseQualsAtMappingQual, int minBaseQual, ErrorModel errorModel) { private int add(ReadBackedPileup pileup, boolean ignoreBadBases, boolean capBaseQualsAtMappingQual, int minBaseQual, ErrorModel errorModel) {
// Number of [A C G T]'s in pileup, in that order // Number of [A C G T]'s in pileup, in that order
@ -235,28 +235,29 @@ public class GeneralPloidySNPGenotypeLikelihoods extends GeneralPloidyGenotypeLi
for (byte b: BaseUtils.BASES) for (byte b: BaseUtils.BASES)
numSeenBases.add(0); numSeenBases.add(0);
if (hasReferenceSampleData) { int nGoodBases = 0;
// count number of elements in pileup // count number of elements in pileup
for (PileupElement elt : pileup) { for (PileupElement elt : pileup) {
byte obsBase = elt.getBase(); byte obsBase = elt.getBase();
byte qual = qualToUse(elt, ignoreBadBases, capBaseQualsAtMappingQual, minBaseQual); byte qual = qualToUse(elt, ignoreBadBases, capBaseQualsAtMappingQual, minBaseQual);
if ( qual == 0 ) if ( qual == 0 )
continue; continue;
int idx = 0; int idx = 0;
for (byte base:BaseUtils.BASES) { for (byte base:BaseUtils.BASES) {
int cnt = numSeenBases.get(idx); int cnt = numSeenBases.get(idx);
numSeenBases.set(idx++,cnt + (base == obsBase?1:0)); numSeenBases.set(idx++,cnt + (base == obsBase?1:0));
}
} }
if (VERBOSE) nGoodBases++;
System.out.format("numSeenBases: %d %d %d %d\n",numSeenBases.get(0),numSeenBases.get(1),numSeenBases.get(2),numSeenBases.get(3));
} }
if (VERBOSE)
System.out.format("numSeenBases: %d %d %d %d\n",numSeenBases.get(0),numSeenBases.get(1),numSeenBases.get(2),numSeenBases.get(3));
computeLikelihoods(errorModel, myAlleles, numSeenBases, pileup); computeLikelihoods(errorModel, myAlleles, numSeenBases, pileup);
return pileup.getNumberOfElements(); return nGoodBases;
} }
/** /**
@ -281,7 +282,8 @@ public class GeneralPloidySNPGenotypeLikelihoods extends GeneralPloidyGenotypeLi
double p1 = 0.0; double p1 = 0.0;
if (!hasReferenceSampleData) { if (!hasReferenceSampleData) {
// no error model: loop throught pileup to compute likalihoods just on base qualities // no error model: loop through pileup to compute likelihoods just on base qualities
// In this case, vector numObservations is not used directly for GL computation
for (final PileupElement elt : pileup) { for (final PileupElement elt : pileup) {
final byte obsBase = elt.getBase(); final byte obsBase = elt.getBase();
final byte qual = qualToUse(elt, true, true, mbq); final byte qual = qualToUse(elt, true, true, mbq);