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

Fixed GQ estimate when chosen genotype isn't the most likely according to the GLs. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3362 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2010-05-14 19:17:46 +00:00
parent 1538dc0144
commit 32389dc0a9
2 changed files with 19 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java
View File

@ -289,12 +289,23 @@ public class DiploidGenotypeCalculationModel extends JointEstimateGenotypeCalcul
int genotype = indexes[index];
double score;
if ( genotype == GenotypeType.REF.ordinal() )
score = matrix[index][GenotypeType.REF.ordinal()] - Math.max(matrix[index][GenotypeType.HET.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]);
else if ( genotype == GenotypeType.HET.ordinal() )
score = matrix[index][GenotypeType.HET.ordinal()] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]);
else // ( genotype == GenotypeType.HOM.ordinal() )
score = matrix[index][GenotypeType.HOM.ordinal()] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HET.ordinal()]);
int maxEntry = MathUtils.maxElementIndex(matrix[index]);
// if the max value is for the most likely genotype, we can compute next vs. next best
if ( genotype == maxEntry ) {
if ( genotype == GenotypeType.REF.ordinal() )
score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.HET.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]);
else if ( genotype == GenotypeType.HET.ordinal() )
score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]);
else // ( genotype == GenotypeType.HOM.ordinal() )
score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HET.ordinal()]);
}
// otherwise, we need to calculate the probability of the genotype
else {
double[] normalized = MathUtils.normalizeFromLog10(matrix[index]);
double chosenGenotype = normalized[genotype];
score = -1.0 * Math.log10(1.0 - chosenGenotype);
}
samplesToGenotypes.put(sample, new Pair<Integer, Double>(genotype, Math.abs(score)));
index++;

4
java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
View File

@ -77,8 +77,8 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
@Test
public void testParameter() {
HashMap<String, String> e = new HashMap<String, String>();
e.put( "-genotype", "310d3f8ac2a1100d82fe6241e201e1d8" );
e.put( "-all_bases", "470196b3f9a2ed622383e25552972a3e" );
e.put( "-genotype", "03e228232eeb0a32ea3791daec29248c" );
e.put( "-all_bases", "fe894f335e048d89fccd7a7b39a105a1" );
e.put( "--min_base_quality_score 26", "a1142c69f1951a05d5a33215c69133af" );
e.put( "--min_mapping_quality_score 26", "796d1ea790481a1a257609279f8076bb" );
e.put( "--max_mismatches_in_40bp_window 5", "379bd56d4402e974127e0677ff51d040" );