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

Checking in the tiger team changes. LOD calculation modified. -qScale is back in case people need it. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3990 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
rpoplin 2010-08-09 20:41:03 +00:00
parent 0eeb659aa3
commit 3eee3183fd
3 changed files with 7 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantGaussianMixtureModel.java
View File

@ -545,9 +545,6 @@ public final class VariantGaussianMixtureModel extends VariantOptimizationModel
sum += pVarInCluster[kkk]; // * clusterTruePositiveRate[kkk]; sum += pVarInCluster[kkk]; // * clusterTruePositiveRate[kkk];
} }
// if ( sum > 1 )
// System.out.printf("Bad pVar, bad!");
return sum; return sum;
} }
@ -837,7 +834,6 @@ public final class VariantGaussianMixtureModel extends VariantOptimizationModel
private void evaluateGaussiansForSingleVariant( final double[] annotations, final double[] pVarInCluster ) { private void evaluateGaussiansForSingleVariant( final double[] annotations, final double[] pVarInCluster ) {
final int numAnnotations = annotations.length; final int numAnnotations = annotations.length;
//final double mult[] = new double[numAnnotations];
final double evalGaussianPDFLog10[] = new double[maxGaussians]; final double evalGaussianPDFLog10[] = new double[maxGaussians];
for( int kkk = 0; kkk < maxGaussians; kkk++ ) { for( int kkk = 0; kkk < maxGaussians; kkk++ ) {

20
java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java
View File

@ -92,6 +92,8 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
private double QUAL_STEP = 0.1; private double QUAL_STEP = 0.1;
@Argument(fullName="singleton_fp_rate", shortName="fp_rate", doc="Prior expectation that a singleton call would be a FP", required=false) @Argument(fullName="singleton_fp_rate", shortName="fp_rate", doc="Prior expectation that a singleton call would be a FP", required=false)
private double SINGLETON_FP_RATE = 0.5; private double SINGLETON_FP_RATE = 0.5;
@Argument(fullName="quality_scale_factor", shortName="qScale", doc="Multiply all final quality scores by this value. Needed to normalize the quality scores.", required=false)
private double QUALITY_SCALE_FACTOR = 20.0;
///////////////////////////// /////////////////////////////
// Private Member Variables // Private Member Variables
@ -99,7 +101,6 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
private VariantGaussianMixtureModel theModel = null; private VariantGaussianMixtureModel theModel = null;
private VCFWriter vcfWriter; private VCFWriter vcfWriter;
private Set<String> ignoreInputFilterSet = null; private Set<String> ignoreInputFilterSet = null;
private int numUnstable = 0;
//--------------------------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------------------------
// //
@ -195,18 +196,13 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
final double totalPrior = 1.0 - ((1.0 - acPrior) * (1.0 - knownPrior)); final double totalPrior = 1.0 - ((1.0 - acPrior) * (1.0 - knownPrior));
if( MathUtils.compareDoubles(totalPrior, 1.0, 1E-8) == 0 || MathUtils.compareDoubles(totalPrior, 0.0, 1E-8) == 0 ) { if( MathUtils.compareDoubles(totalPrior, 1.0, 1E-8) == 0 || MathUtils.compareDoubles(totalPrior, 0.0, 1E-8) == 0 ) {
throw new StingException("Some is wrong with the prior that was entered by the user: Prior = " + totalPrior); // TODO - fix this up later throw new StingException("Something is wrong with the prior that was entered by the user: Prior = " + totalPrior); // TODO - fix this up later
} }
double pVar = theModel.evaluateVariant( vc ); final double pVar = theModel.evaluateVariant( vc );
if( pVar > 1.0 ) { final double lod = (Math.log10(totalPrior) + Math.log10(pVar)) - ((Math.log10(1.0 - totalPrior)) + Math.log10(1.0));
pVar = 0.99; variantDatum.qual = QUALITY_SCALE_FACTOR * QualityUtils.lodToPhredScaleErrorRate(lod);
numUnstable++;
}
final double lod = (Math.log10(totalPrior) + Math.log10(pVar)) - ((Math.log10(1.0 - totalPrior)) + Math.log10(1.0 - pVar));
variantDatum.qual = 10.0 * QualityUtils.lodToPhredScaleErrorRate(lod);
mapList.add( variantDatum ); mapList.add( variantDatum );
@ -247,10 +243,6 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
vcfWriter.close(); vcfWriter.close();
if( numUnstable > 0 ) {
logger.warn("WARNING: Found " + numUnstable + " variant(s) with pVar > 1, Most likely numerical instability during clustering !!!!!");
}
final VariantDataManager dataManager = new VariantDataManager( reduceSum, theModel.dataManager.annotationKeys ); final VariantDataManager dataManager = new VariantDataManager( reduceSum, theModel.dataManager.annotationKeys );
reduceSum.clear(); // Don't need this ever again, clean up some memory reduceSum.clear(); // Don't need this ever again, clean up some memory

2
java/test/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrationWalkersIntegrationTest.java
View File

@ -37,7 +37,7 @@ public class VariantRecalibrationWalkersIntegrationTest extends WalkerTest {
@Test @Test
public void testVariantRecalibrator() { public void testVariantRecalibrator() {
HashMap<String, String> e = new HashMap<String, String>(); HashMap<String, String> e = new HashMap<String, String>();
e.put( validationDataLocation + "yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf", "4570c93c96de61342d75c1c658e2d5c1" ); e.put( validationDataLocation + "yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf", "df4e6f55c714b68b13e47b61d4bd0cd5" );
for ( Map.Entry<String, String> entry : e.entrySet() ) { for ( Map.Entry<String, String> entry : e.entrySet() ) {
String vcf = entry.getKey(); String vcf = entry.getKey();