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

Merge branch 'master' of ssh://copper.broadinstitute.org/humgen/gsa-scr1/gsa-engineering/git/unstable

This commit is contained in:
Menachem Fromer 2012-04-18 12:42:28 -04:00
parent 5d3c764ed2 8a84456626
commit 70777c7af8
6 changed files with 139 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

129
public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyRecalibration.java
View File

@ -35,16 +35,14 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.PartitionBy;
import org.broadinstitute.sting.gatk.walkers.PartitionType;
import org.broadinstitute.sting.gatk.walkers.RodWalker;
import org.broadinstitute.sting.gatk.walkers.TreeReducible;
import org.broadinstitute.sting.utils.SampleUtils;
import org.broadinstitute.sting.utils.codecs.vcf.*;
import org.broadinstitute.sting.utils.collections.NestedHashMap;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.text.XReadLines;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContextBuilder;
import java.io.File;
import java.io.FileNotFoundException;
import java.util.*;
/**
@ -86,8 +84,8 @@ import java.util.*;
*
*/
@PartitionBy(PartitionType.NONE)
public class ApplyRecalibration extends RodWalker<Integer, Integer> {
@PartitionBy(PartitionType.LOCUS)
public class ApplyRecalibration extends RodWalker<Integer, Integer> implements TreeReducible<Integer> {
/////////////////////////////
// Inputs
@ -98,9 +96,9 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
@Input(fullName="input", shortName = "input", doc="The raw input variants to be recalibrated", required=true)
public List<RodBinding<VariantContext>> input;
@Input(fullName="recal_file", shortName="recalFile", doc="The input recal file used by ApplyRecalibration", required=true)
private File RECAL_FILE;
protected RodBinding<VariantContext> recal;
@Input(fullName="tranches_file", shortName="tranchesFile", doc="The input tranches file describing where to cut the data", required=true)
private File TRANCHES_FILE;
protected File TRANCHES_FILE;
/////////////////////////////
// Outputs
@ -112,7 +110,7 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
// Command Line Arguments
/////////////////////////////
@Argument(fullName="ts_filter_level", shortName="ts_filter_level", doc="The truth sensitivity level at which to start filtering", required=false)
private double TS_FILTER_LEVEL = 99.0;
protected double TS_FILTER_LEVEL = 99.0;
@Argument(fullName="ignore_filter", shortName="ignoreFilter", doc="If specified the variant recalibrator will use variants even if the specified filter name is marked in the input VCF file", required=false)
private String[] IGNORE_INPUT_FILTERS = null;
@Argument(fullName = "mode", shortName = "mode", doc = "Recalibration mode to employ: 1.) SNP for recalibrating only SNPs (emitting indels untouched in the output VCF); 2.) INDEL for indels; and 3.) BOTH for recalibrating both SNPs and indels simultaneously.", required = false)
@ -123,8 +121,6 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
/////////////////////////////
final private List<Tranche> tranches = new ArrayList<Tranche>();
final private Set<String> inputNames = new HashSet<String>();
final private NestedHashMap lodMap = new NestedHashMap();
final private NestedHashMap annotationMap = new NestedHashMap();
final private Set<String> ignoreInputFilterSet = new TreeSet<String>();
//---------------------------------------------------------------------------------------------------------------
@ -174,20 +170,6 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
final VCFHeader vcfHeader = new VCFHeader(hInfo, samples);
vcfWriter.writeHeader(vcfHeader);
try {
logger.info("Reading in recalibration table...");
for ( final String line : new XReadLines( RECAL_FILE ) ) {
final String[] vals = line.split(",");
lodMap.put( Double.parseDouble(vals[3]), vals[0], Integer.parseInt(vals[1]), Integer.parseInt(vals[2]) ); // value comes before the keys
annotationMap.put( vals[4], vals[0], Integer.parseInt(vals[1]), Integer.parseInt(vals[2]) ); // value comes before the keys
}
} catch ( FileNotFoundException e ) {
throw new UserException.CouldNotReadInputFile(RECAL_FILE, e);
} catch ( Exception e ) {
throw new UserException.MalformedFile(RECAL_FILE, "Could not parse LOD and annotation information in input recal file. File is somehow malformed.");
}
}
//---------------------------------------------------------------------------------------------------------------
@ -202,52 +184,75 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
return 1;
}
for( VariantContext vc : tracker.getValues(input, context.getLocation()) ) {
if( vc != null ) {
if( VariantRecalibrator.checkRecalibrationMode( vc, MODE ) && (vc.isNotFiltered() || ignoreInputFilterSet.containsAll(vc.getFilters())) ) {
VariantContextBuilder builder = new VariantContextBuilder(vc);
String filterString = null;
final List<VariantContext> VCs = tracker.getValues(input, context.getLocation());
final List<VariantContext> recals = tracker.getValues(recal, context.getLocation());
final Double lod = (Double) lodMap.get( vc.getChr(), vc.getStart(), vc.getEnd() );
final String worstAnnotation = (String) annotationMap.get( vc.getChr(), vc.getStart(), vc.getEnd() );
if( lod == null ) {
throw new UserException("Encountered input variant which isn't found in the input recal file. Please make sure VariantRecalibrator and ApplyRecalibration were run on the same set of input variants. First seen at: " + vc );
}
for( final VariantContext vc : VCs ) {
// Annotate the new record with its VQSLOD and the worst performing annotation
builder.attribute(VariantRecalibrator.VQS_LOD_KEY, String.format("%.4f", lod));
builder.attribute(VariantRecalibrator.CULPRIT_KEY, worstAnnotation);
if( VariantRecalibrator.checkRecalibrationMode( vc, MODE ) && (vc.isNotFiltered() || ignoreInputFilterSet.containsAll(vc.getFilters())) ) {
for( int i = tranches.size() - 1; i >= 0; i-- ) {
final Tranche tranche = tranches.get(i);
if( lod >= tranche.minVQSLod ) {
if( i == tranches.size() - 1 ) {
filterString = VCFConstants.PASSES_FILTERS_v4;
} else {
filterString = tranche.name;
}
break;
}
}
if( filterString == null ) {
filterString = tranches.get(0).name+"+";
}
if( !filterString.equals(VCFConstants.PASSES_FILTERS_v4) ) {
builder.filters(filterString);
}
vcfWriter.add( builder.make() );
} else { // valid VC but not compatible with this mode, so just emit the variant untouched
vcfWriter.add( vc );
final VariantContext recalDatum = getMatchingRecalVC(vc, recals);
if( recalDatum == null ) {
throw new UserException("Encountered input variant which isn't found in the input recal file. Please make sure VariantRecalibrator and ApplyRecalibration were run on the same set of input variants. First seen at: " + vc );
}
final String lodString = recalDatum.getAttributeAsString(VariantRecalibrator.VQS_LOD_KEY, null);
if( lodString == null ) {
throw new UserException("Encountered a malformed record in the input recal file. There is no lod for the record at: " + vc );
}
final double lod;
try {
lod = Double.valueOf(lodString);
} catch (NumberFormatException e) {
throw new UserException("Encountered a malformed record in the input recal file. The lod is unreadable for the record at: " + vc );
}
VariantContextBuilder builder = new VariantContextBuilder(vc);
String filterString = null;
// Annotate the new record with its VQSLOD and the worst performing annotation
builder.attribute(VariantRecalibrator.VQS_LOD_KEY, lodString); // use the String representation so that we don't lose precision on output
builder.attribute(VariantRecalibrator.CULPRIT_KEY, recalDatum.getAttribute(VariantRecalibrator.CULPRIT_KEY));
for( int i = tranches.size() - 1; i >= 0; i-- ) {
final Tranche tranche = tranches.get(i);
if( lod >= tranche.minVQSLod ) {
if( i == tranches.size() - 1 ) {
filterString = VCFConstants.PASSES_FILTERS_v4;
} else {
filterString = tranche.name;
}
break;
}
}
if( filterString == null ) {
filterString = tranches.get(0).name+"+";
}
if( !filterString.equals(VCFConstants.PASSES_FILTERS_v4) ) {
builder.filters(filterString);
}
vcfWriter.add( builder.make() );
} else { // valid VC but not compatible with this mode, so just emit the variant untouched
vcfWriter.add( vc );
}
}
return 1; // This value isn't used for anything
}
private static VariantContext getMatchingRecalVC(final VariantContext target, final List<VariantContext> recalVCs) {
for( final VariantContext recalVC : recalVCs ) {
if ( target.getEnd() == recalVC.getEnd() ) {
return recalVC;
}
}
return null;
}
//---------------------------------------------------------------------------------------------------------------
//
// reduce
@ -262,6 +267,10 @@ public class ApplyRecalibration extends RodWalker<Integer, Integer> {
return 1; // This value isn't used for anything
}
public Integer treeReduce( final Integer lhs, final Integer rhs ) {
return 1; // This value isn't used for anything
}
public void onTraversalDone( final Integer reduceSum ) {
}
}

35
public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java
View File

@ -30,14 +30,16 @@ import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeader;
import org.broadinstitute.sting.utils.codecs.vcf.VCFWriter;
import org.broadinstitute.sting.utils.collections.ExpandingArrayList;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContextBuilder;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.*;
/**
* Created by IntelliJ IDEA.
@ -285,11 +287,28 @@ public class VariantDataManager {
(TRUST_ALL_POLYMORPHIC || !trainVC.hasGenotypes() || trainVC.isPolymorphicInSamples());
}
public void writeOutRecalibrationTable( final PrintStream RECAL_FILE ) {
public void writeOutRecalibrationTable( final VCFWriter recalWriter ) {
// we need to sort in coordinate order in order to produce a valid VCF
Collections.sort( data, new Comparator<VariantDatum>() {
public int compare(VariantDatum vd1, VariantDatum vd2) {
return vd1.loc.compareTo(vd2.loc);
}} );
// create dummy alleles to be used
final List<Allele> alleles = new ArrayList<Allele>(2);
alleles.add(Allele.create("N", true));
alleles.add(Allele.create("<VQSR>", false));
// to be used for the important INFO tags
final HashMap<String, Object> attributes = new HashMap<String, Object>(3);
for( final VariantDatum datum : data ) {
RECAL_FILE.println(String.format("%s,%d,%d,%.4f,%s",
datum.contig, datum.start, datum.stop, datum.lod,
(datum.worstAnnotation != -1 ? annotationKeys.get(datum.worstAnnotation) : "NULL")));
attributes.put(VCFConstants.END_KEY, datum.loc.getStop());
attributes.put(VariantRecalibrator.VQS_LOD_KEY, String.format("%.4f", datum.lod));
attributes.put(VariantRecalibrator.CULPRIT_KEY, (datum.worstAnnotation != -1 ? annotationKeys.get(datum.worstAnnotation) : "NULL"));
VariantContextBuilder builder = new VariantContextBuilder("VQSR", datum.loc.getContig(), datum.loc.getStart(), datum.loc.getStart(), alleles).attributes(attributes);
recalWriter.add(builder.make());
}
}
}

6
public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDatum.java
View File

@ -25,6 +25,8 @@
package org.broadinstitute.sting.gatk.walkers.variantrecalibration;
import org.broadinstitute.sting.utils.GenomeLoc;
/**
* Created by IntelliJ IDEA.
* User: rpoplin
@ -46,9 +48,7 @@ public class VariantDatum implements Comparable<VariantDatum> {
public double originalQual;
public double prior;
public int consensusCount;
public String contig;
public int start;
public int stop;
public GenomeLoc loc;
public int worstAnnotation;
public MultivariateGaussian assignment; // used in K-means implementation

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

@ -37,6 +37,8 @@ import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.R.RScriptExecutor;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.codecs.vcf.StandardVCFWriter;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeader;
import org.broadinstitute.sting.utils.collections.ExpandingArrayList;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.io.Resource;
@ -136,9 +138,11 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
// Outputs
/////////////////////////////
@Output(fullName="recal_file", shortName="recalFile", doc="The output recal file used by ApplyRecalibration", required=true)
private PrintStream RECAL_FILE;
protected File recalFile = null;
protected StandardVCFWriter recalWriter = null;
@Output(fullName="tranches_file", shortName="tranchesFile", doc="The output tranches file used by ApplyRecalibration", required=true)
private File TRANCHES_FILE;
protected File TRANCHES_FILE;
/////////////////////////////
// Additional Command Line Arguments
@ -150,7 +154,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
* that this parameter is used for display purposes only and isn't used anywhere in the algorithm!
*/
@Argument(fullName="target_titv", shortName="titv", doc="The expected novel Ti/Tv ratio to use when calculating FDR tranches and for display on the optimization curve output figures. (approx 2.15 for whole genome experiments). ONLY USED FOR PLOTTING PURPOSES!", required=false)
private double TARGET_TITV = 2.15;
protected double TARGET_TITV = 2.15;
/**
* See the input VCF file's INFO field for a list of all available annotations.
@ -170,7 +174,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
@Output(fullName="rscript_file", shortName="rscriptFile", doc="The output rscript file generated by the VQSR to aid in visualization of the input data and learned model", required=false)
private File RSCRIPT_FILE = null;
@Argument(fullName="ts_filter_level", shortName="ts_filter_level", doc="The truth sensitivity level at which to start filtering, used here to indicate filtered variants in the model reporting plots", required=false)
private double TS_FILTER_LEVEL = 99.0;
protected double TS_FILTER_LEVEL = 99.0;
/////////////////////////////
// Debug Arguments
@ -224,6 +228,10 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
if( !dataManager.checkHasTruthSet() ) {
throw new UserException.CommandLineException( "No truth set found! Please provide sets of known polymorphic loci marked with the truth=true ROD binding tag. For example, -B:hapmap,VCF,known=false,training=true,truth=true,prior=12.0 hapmapFile.vcf" );
}
final VCFHeader vcfHeader = new VCFHeader( null, Collections.<String>emptySet() );
recalWriter = new StandardVCFWriter(recalFile, getMasterSequenceDictionary(), false);
recalWriter.writeHeader(vcfHeader);
}
//---------------------------------------------------------------------------------------------------------------
@ -246,9 +254,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
// Populate the datum with lots of fields from the VariantContext, unfortunately the VC is too big so we just pull in only the things we absolutely need.
dataManager.decodeAnnotations( datum, vc, true ); //BUGBUG: when run with HierarchicalMicroScheduler this is non-deterministic because order of calls depends on load of machine
datum.contig = vc.getChr();
datum.start = vc.getStart();
datum.stop = vc.getEnd();
datum.loc = getToolkit().getGenomeLocParser().createGenomeLoc(vc);
datum.originalQual = vc.getPhredScaledQual();
datum.isSNP = vc.isSNP() && vc.isBiallelic();
datum.isTransition = datum.isSNP && VariantContextUtils.isTransition(vc);
@ -345,7 +351,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
}
logger.info( "Writing out recalibration table..." );
dataManager.writeOutRecalibrationTable( RECAL_FILE );
dataManager.writeOutRecalibrationTable( recalWriter );
if( RSCRIPT_FILE != null ) {
logger.info( "Writing out visualization Rscript file...");
createVisualizationScript( dataManager.getRandomDataForPlotting( 6000 ), goodModel, badModel, lodCutoff );

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

@ -66,6 +66,14 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
executeTest("test Multiple SNP alleles", spec);
}
@Test
public void testBadRead() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " -nosl -NO_HEADER -glm BOTH -I " + validationDataLocation + "badRead.test.bam -o %s -L 1:22753424-22753464", 1,
Arrays.asList("7678827a2ee21870a41c09d28d26b996"));
executeTest("test bad read", spec);
}
// --------------------------------------------------------------------------------------------------------------
//
// testing compressed output

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

@ -27,7 +27,7 @@ public class VariantRecalibrationWalkersIntegrationTest extends WalkerTest {
VRTest lowPass = new VRTest("phase1.projectConsensus.chr20.raw.snps.vcf",
"0ddd1e0e483d2eaf56004615cea23ec7", // tranches
"58780f63182e139fdbe17f6c18b5b774", // recal file
"f8e21a1987960b950db1f0d98be45352", // recal file
"f67d844b6252a55452cf4167b77530b1"); // cut VCF
@DataProvider(name = "VRTest")
@ -74,7 +74,7 @@ public class VariantRecalibrationWalkersIntegrationTest extends WalkerTest {
VRTest indel = new VRTest("combined.phase1.chr20.raw.indels.sites.vcf",
"6d7ee4cb651c8b666e4a4523363caaff", // tranches
"4759b111a5aa53975d46e0f22c7983bf", // recal file
"ee5b408c8434a594496118875690c438", // recal file
"5d7e07d8813db96ba3f3dfe4737f83d1"); // cut VCF
@DataProvider(name = "VRIndelTest")
@ -118,5 +118,21 @@ public class VariantRecalibrationWalkersIntegrationTest extends WalkerTest {
Arrays.asList(params.cutVCFMD5));
executeTest("testApplyRecalibrationIndel-"+params.inVCF, spec);
}
@Test
public void testApplyRecalibrationSnpAndIndelTogether() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-R " + b37KGReference +
" -T ApplyRecalibration" +
" -L 20:1000100-1000500" +
" -mode BOTH" +
" -NO_HEADER" +
" -input " + validationDataLocation + "VQSR.mixedTest.input" +
" -o %s" +
" -tranchesFile " + validationDataLocation + "VQSR.mixedTest.tranches" +
" -recalFile " + validationDataLocation + "VQSR.mixedTest.recal",
Arrays.asList("08060b7f5c9cf3bb1692b50c58fd5a4b"));
executeTest("testApplyRecalibrationSnpAndIndelTogether", spec);
}
}