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

Merge branch 'master' of gsa2:/humgen/gsa-scr1/chartl/dev/unstable

This commit is contained in:
Chris Hartl 2012-12-05 11:20:28 -05:00
parent c30dc3c38f 0c925856cb
commit 430d6a07f2
36 changed files with 518 additions and 302 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

97
build.xml
View File

@ -185,10 +185,7 @@
<include name="**/*.class"/>
</fileset>
<patternset id="dependency.mask" includes="*.jar">
<exclude name="testng*.jar" />
<exclude name="bcel*.jar" />
</patternset>
<patternset id="dependency.mask" includes="*.jar" />
<path id="external.dependencies">
<fileset dir="${lib.dir}" erroronmissingdir="false">
@ -205,6 +202,16 @@
<pathelement location="${scala.classes}" />
</path>
<path id="build.results">
<!-- Ensure that GenomeAnalysisTK.jar comes first in the path, as it contains overrides for certain classes in our dependencies -->
<pathelement location="${dist.dir}/GenomeAnalysisTK.jar" />
<!-- After GenomeAnalysisTK.jar we include all of the other jars in the dist directory -->
<fileset dir="${dist.dir}" erroronmissingdir="false">
<patternset refid="dependency.mask" />
<exclude name="GenomeAnalysisTK.jar" />
</fileset>
</path>
<fileset id="external.source.files" dir="${external.dir}" erroronmissingdir="false">
<include name="**/*.java" />
</fileset>
@ -226,20 +233,20 @@
<!-- the path for resources that need to go into the GATK jar;
any additional resources should go into this set. -->
<path id="gatk.resources">
<fileset dir="${basedir}">
<include name="${java.public.source.dir}/**/templates/*" />
<include name="${java.private.source.dir}/**/templates/*" if="include.private" />
<include name="${java.protected.source.dir}/**/templates/*" if="include.protected" />
<fileset dir="${java.public.source.dir}">
<include name="**/resources/*" />
<include name="**/templates/*" />
</fileset>
<fileset dir="${java.private.source.dir}" erroronmissingdir="false">
<include name="**/resources/*" if="include.private" />
<include name="**/templates/*" if="include.private" />
</fileset>
<fileset dir="${java.protected.source.dir}" erroronmissingdir="false">
<include name="**/resources/*" if="include.protected" />
<include name="**/templates/*" if="include.protected" />
</fileset>
</path>
<path id="build.results">
<fileset dir="${dist.dir}">
<patternset refid="dependency.mask" />
</fileset>
</path>
<!-- ******************************************************************************** -->
<!-- Ivy Retrieve -->
<!-- ******************************************************************************** -->
@ -672,6 +679,24 @@
</jar>
</target>
<target name="na12878kb.jar" depends="gatk.compile,init.jar">
<jar jarfile="${dist.dir}/na12878kb.jar">
<fileset dir="${java.classes}">
<include name="org/broadinstitute/sting/utils/GenomeLocParser*.class"/>
<include name="org/broadinstitute/sting/utils/GenomeLoc.class"/>
<include name="org/broadinstitute/sting/utils/HasGenomeLocation.class"/>
<include name="org/broadinstitute/sting/utils/BaseUtils.class"/>
<include name="org/broadinstitute/sting/utils/Utils.class"/>
<include name="org/broadinstitute/sting/utils/exceptions/**/*.class"/>
<include name="org/broadinstitute/sting/gatk/walkers/na12878kb/core/**/*.class"/>
<include name="net/sf/picard/reference/FastaSequenceFile.class"/>
</fileset>
<fileset dir="${java.private.source.dir}">
<include name="org/broadinstitute/sting/gatk/walkers/na12878kb/core/resources/**/*"/>
</fileset>
</jar>
</target>
<target name="gatk.jar" depends="gatk.compile, init.jar, R.script.stage" description="generate the GATK distribution">
<jar jarfile="${dist.dir}/GenomeAnalysisTK.jar">
<path refid="gatk.resources"/>
@ -1103,15 +1128,10 @@
</path>
<path id="testng.default.classpath">
<pathelement location="${java.classes}" />
<pathelement location="${scala.classes}" />
<path refid="build.results" />
<pathelement location="${java.contracts.dir}" />
<pathelement location="${java.test.classes}" />
<pathelement location="${scala.test.classes}" />
<pathelement location="${R.tar.dir}" />
<path refid="R.script.source.path" />
<pathelement location="${key.dir}" />
<path refid="external.dependencies" />
</path>
<!-- Test targets -->
@ -1119,9 +1139,6 @@
<target name="test.init.compile">
<mkdir dir="${java.test.classes}"/>
<mkdir dir="${scala.test.classes}"/>
<antcall target="resolve">
<param name="ivy.conf" value="test"/>
</antcall>
</target>
<target name="test.java.internal.compile" depends="dist,test.init.compile">
@ -1129,10 +1146,8 @@
<javac fork="true" memoryMaximumSize="512m" destdir="${java.test.classes}" debug="true" optimize="on" tempdir="${java.io.tmpdir}">
<src refid="java.test.source.path" />
<classpath>
<path refid="external.dependencies" />
<pathelement location="${java.classes}"/>
<path refid="build.results" />
<pathelement location="${java.contracts.dir}"/>
<pathelement location="${testng.jar}"/>
</classpath>
<compilerarg value="-proc:none"/>
</javac>
@ -1143,11 +1158,9 @@
<javac fork="true" memoryMaximumSize="512m" destdir="${java.test.classes}" debug="true" optimize="on" tempdir="${java.io.tmpdir}" srcdir="${external.dir}">
<include name="*/test/**/*.java"/>
<classpath>
<path refid="external.dependencies" />
<path refid="build.results" />
<pathelement location="${java.test.classes}"/>
<pathelement location="${java.classes}"/>
<pathelement location="${java.contracts.dir}"/>
<pathelement location="${testng.jar}"/>
</classpath>
<compilerarg value="-proc:none"/>
</javac>
@ -1160,9 +1173,8 @@
<scalac fork="true" jvmargs="-Xmx512m" destdir="${scala.test.classes}" deprecation="yes" unchecked="yes">
<src refid="scala.test.source.path" />
<classpath>
<path refid="scala.dependencies"/>
<path refid="build.results"/>
<pathelement location="${java.test.classes}"/>
<pathelement location="${testng.jar}"/>
</classpath>
</scalac>
</target>
@ -1369,14 +1381,13 @@
<!-- Fast test target that cuts major corners for speed. Requires that a full build has been done first. Java-only, single test class only -->
<!-- Usage: ant fasttest -Dsingle=TestClass -->
<target name="fasttest" depends="init.javaonly,init,test.init">
<target name="fasttest" depends="init.javaonly,init">
<condition property="not.clean">
<and>
<available file="${build.dir}" />
<available file="${lib.dir}" />
<available file="${dist.dir}" />
<available file="${java.test.classes}" />
<available file="${testng.jar}" />
</and>
</condition>
<fail message="fasttest requires a NON-CLEAN working directory (INCLUDING test classes). Do a full test build using ant test.compile first." unless="not.clean" />
@ -1394,13 +1405,27 @@
<javac fork="true" memoryMaximumSize="512m" destdir="${java.test.classes}" debug="true" optimize="on" tempdir="${java.io.tmpdir}">
<src refid="java.test.source.path" />
<classpath>
<path refid="external.dependencies" />
<pathelement location="${java.classes}"/>
<pathelement location="${testng.jar}"/>
<path refid="external.dependencies" />
</classpath>
<compilerarg value="-proc:none"/>
</javac>
<!-- fasttest uses the unpackaged class files in its test classpath to avoid having to rebuild the jars in dist/ -->
<path id="testng.fasttest.classpath">
<pathelement location="${java.classes}" />
<pathelement location="${scala.classes}" />
<pathelement location="${java.contracts.dir}" />
<pathelement location="${java.test.classes}" />
<pathelement location="${scala.test.classes}" />
<pathelement location="${R.tar.dir}" />
<path refid="R.script.source.path" />
<pathelement location="${key.dir}" />
<path refid="external.dependencies" />
<path refid="java.source.path" /> <!-- Terrible hack to allow fasttest to see resource files stored in the source tree -->
</path>
<property name="testng.classpath" value="testng.fasttest.classpath" />
<run-test testtype="${single}" outputdir="${report}/${single}" runfailed="false"/>
</target>
</project>

11
ivy.xml
View File

@ -24,11 +24,8 @@
<ivy-module version="1.0">
<info organisation="org.broadinstitute" module="Sting"/>
<configurations defaultconfmapping="test->default">
<configurations>
<conf name="default" description="the core dependencies for the GATK"/>
<conf name="test" extends="default" description="external dependencies used for testing and metrics"/>
<conf name="scala" extends="default" description="the dependencies for scala"/>
<conf name="queue" extends="scala" description="the dependencies for Queue"/>
</configurations>
<dependencies defaultconf="default">
<dependency org="net.sf" name="sam" rev="latest.integration"/>
@ -83,9 +80,9 @@
<dependency org="org.scala-lang" name="scala-library" rev="2.9.2"/>
<!-- testing and evaluation dependencies -->
<dependency org="org.testng" name="testng" rev="5.14.1" conf="test"/>
<dependency org="org.uncommons" name="reportng" rev="1.1.2" conf="test"/>
<dependency org="com.google.code.caliper" name="caliper" rev="1.0-SNAPSHOT" conf="test"/>
<dependency org="org.testng" name="testng" rev="5.14.1"/>
<dependency org="org.uncommons" name="reportng" rev="1.1.2"/>
<dependency org="com.google.code.caliper" name="caliper" rev="1.0-SNAPSHOT"/>
<!-- Contracts for Java and dependencies -->
<dependency org="com.google.code.cofoja" name="cofoja" rev="1.0-r139"/>

141
protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java
View File

@ -30,12 +30,16 @@ import com.google.java.contract.Requires;
import net.sf.samtools.Cigar;
import net.sf.samtools.CigarElement;
import org.apache.commons.lang.ArrayUtils;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine;
import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.sting.utils.variantcontext.*;
import java.io.PrintStream;
import java.util.*;
public class GenotypingEngine {
@ -43,23 +47,27 @@ public class GenotypingEngine {
private final boolean DEBUG;
private final static List<Allele> noCall = new ArrayList<Allele>(); // used to noCall all genotypes until the exact model is applied
private final static Allele SYMBOLIC_UNASSEMBLED_EVENT_ALLELE = Allele.create("<UNASSEMBLED_EVENT>", false);
private final VariantAnnotatorEngine annotationEngine;
public GenotypingEngine( final boolean DEBUG ) {
public GenotypingEngine( final boolean DEBUG, final VariantAnnotatorEngine annotationEngine ) {
this.DEBUG = DEBUG;
this.annotationEngine = annotationEngine;
noCall.add(Allele.NO_CALL);
}
// BUGBUG: Create a class to hold this complicated return type
@Requires({"refLoc.containsP(activeRegionWindow)", "haplotypes.size() > 0"})
public List<Pair<VariantContext, Map<Allele, List<Haplotype>>>> assignGenotypeLikelihoodsAndCallIndependentEvents( final UnifiedGenotyperEngine UG_engine,
final List<Haplotype> haplotypes,
final byte[] ref,
final GenomeLoc refLoc,
final GenomeLoc activeRegionWindow,
final GenomeLocParser genomeLocParser,
final List<VariantContext> activeAllelesToGenotype ) {
public List<VariantContext> assignGenotypeLikelihoodsAndCallIndependentEvents( final UnifiedGenotyperEngine UG_engine,
final List<Haplotype> haplotypes,
final List<String> samples,
final Map<String, PerReadAlleleLikelihoodMap> haplotypeReadMap,
final Map<String, ArrayList<GATKSAMRecord>> perSampleFilteredReadList,
final byte[] ref,
final GenomeLoc refLoc,
final GenomeLoc activeRegionWindow,
final GenomeLocParser genomeLocParser,
final List<VariantContext> activeAllelesToGenotype ) {
final ArrayList<Pair<VariantContext, Map<Allele, List<Haplotype>>>> returnCalls = new ArrayList<Pair<VariantContext, Map<Allele, List<Haplotype>>>>();
final List<VariantContext> returnCalls = new ArrayList<VariantContext>();
final boolean in_GGA_mode = !activeAllelesToGenotype.isEmpty();
// Using the cigar from each called haplotype figure out what events need to be written out in a VCF file
@ -79,8 +87,8 @@ public class GenotypingEngine {
}
cleanUpSymbolicUnassembledEvents( haplotypes );
if( !in_GGA_mode && haplotypes.get(0).getSampleKeySet().size() >= 10 ) { // if not in GGA mode and have at least 10 samples try to create MNP and complex events by looking at LD structure
mergeConsecutiveEventsBasedOnLD( haplotypes, startPosKeySet, ref, refLoc );
if( !in_GGA_mode && samples.size() >= 10 ) { // if not in GGA mode and have at least 10 samples try to create MNP and complex events by looking at LD structure
mergeConsecutiveEventsBasedOnLD( haplotypes, samples, haplotypeReadMap, startPosKeySet, ref, refLoc );
}
if( in_GGA_mode ) {
for( final VariantContext compVC : activeAllelesToGenotype ) {
@ -90,7 +98,7 @@ public class GenotypingEngine {
// Walk along each position in the key set and create each event to be outputted
for( final int loc : startPosKeySet ) {
if( loc >= activeRegionWindow.getStart() && loc <= activeRegionWindow.getStop() ) {
if( loc >= activeRegionWindow.getStart() && loc <= activeRegionWindow.getStop() ) { // genotyping an event inside this active region
final ArrayList<VariantContext> eventsAtThisLoc = new ArrayList<VariantContext>(); // the overlapping events to merge into a common reference view
final ArrayList<String> priorityList = new ArrayList<String>(); // used to merge overlapping events into common reference view
@ -153,7 +161,7 @@ public class GenotypingEngine {
if( mergedVC == null ) { continue; }
// let's update the Allele keys in the mapper because they can change after merging when there are complex events
Map<Allele, List<Haplotype>> updatedAlleleMapper = new HashMap<Allele, List<Haplotype>>(alleleMapper.size());
final Map<Allele, List<Haplotype>> updatedAlleleMapper = new HashMap<Allele, List<Haplotype>>(alleleMapper.size());
for ( int i = 0; i < mergedVC.getNAlleles(); i++ ) {
final Allele oldAllele = alleleOrdering.get(i);
final Allele newAllele = mergedVC.getAlleles().get(i);
@ -167,12 +175,14 @@ public class GenotypingEngine {
//System.out.println("Event/haplotype allele mapping = " + alleleMapper);
}
final Map<String, PerReadAlleleLikelihoodMap> alleleReadMap = convertHaplotypeReadMapToAlleleReadMap( haplotypeReadMap, alleleMapper, UG_engine.getUAC().CONTAMINATION_FRACTION, UG_engine.getUAC().contaminationLog );
// Grab the genotype likelihoods from the appropriate places in the haplotype likelihood matrix -- calculation performed independently per sample
final GenotypesContext genotypes = GenotypesContext.create(haplotypes.get(0).getSampleKeySet().size());
for( final String sample : haplotypes.get(0).getSampleKeySet() ) { // BUGBUG: assume all haplotypes saw the same samples
final GenotypesContext genotypes = GenotypesContext.create(samples.size());
for( final String sample : samples ) {
final int numHaplotypes = alleleMapper.size();
final double[] genotypeLikelihoods = new double[numHaplotypes * (numHaplotypes+1) / 2];
final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sample, alleleMapper, alleleOrdering);
final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sample, alleleReadMap, alleleOrdering);
int glIndex = 0;
for( int iii = 0; iii < numHaplotypes; iii++ ) {
for( int jjj = 0; jjj <= iii; jjj++ ) {
@ -181,27 +191,57 @@ public class GenotypingEngine {
}
genotypes.add( new GenotypeBuilder(sample).alleles(noCall).PL(genotypeLikelihoods).make() );
}
VariantContext call = UG_engine.calculateGenotypes(new VariantContextBuilder(mergedVC).genotypes(genotypes).make(), UG_engine.getUAC().GLmodel);
final VariantContext call = UG_engine.calculateGenotypes(new VariantContextBuilder(mergedVC).genotypes(genotypes).make(), UG_engine.getUAC().GLmodel);
if( call != null ) {
if( call.getAlleles().size() != mergedVC.getAlleles().size() ) { // some alleles were removed so reverseTrimming might be necessary!
final VariantContext vcCallTrim = VariantContextUtils.reverseTrimAlleles(call);
// also, need to update the allele -> haplotype mapping
final HashMap<Allele, List<Haplotype>> alleleHashMapTrim = new HashMap<Allele, List<Haplotype>>();
for( int iii = 0; iii < vcCallTrim.getAlleles().size(); iii++ ) { // BUGBUG: this is assuming that the original and trimmed alleles maintain the same ordering in the VC
alleleHashMapTrim.put(vcCallTrim.getAlleles().get(iii), alleleMapper.get(call.getAlleles().get(iii)));
}
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap = filterToOnlyOverlappingReads( genomeLocParser, alleleReadMap, perSampleFilteredReadList, call );
VariantContext annotatedCall = annotationEngine.annotateContext(stratifiedReadMap, call);
call = vcCallTrim;
alleleMapper = alleleHashMapTrim;
if( annotatedCall.getAlleles().size() != mergedVC.getAlleles().size() ) { // some alleles were removed so reverseTrimming might be necessary!
annotatedCall = VariantContextUtils.reverseTrimAlleles(annotatedCall);
}
returnCalls.add( new Pair<VariantContext, Map<Allele,List<Haplotype>>>(call, alleleMapper) );
returnCalls.add( annotatedCall );
}
}
}
return returnCalls;
}
private static Map<String, PerReadAlleleLikelihoodMap> filterToOnlyOverlappingReads( final GenomeLocParser parser,
final Map<String, PerReadAlleleLikelihoodMap> perSampleReadMap,
final Map<String, ArrayList<GATKSAMRecord>> perSampleFilteredReadList,
final VariantContext call ) {
final Map<String, PerReadAlleleLikelihoodMap> returnMap = new HashMap<String, PerReadAlleleLikelihoodMap>();
final GenomeLoc callLoc = parser.createGenomeLoc(call);
for( final Map.Entry<String, PerReadAlleleLikelihoodMap> sample : perSampleReadMap.entrySet() ) {
final PerReadAlleleLikelihoodMap likelihoodMap = PerReadAlleleLikelihoodMap.getBestAvailablePerReadAlleleLikelihoodMap();
for( final Map.Entry<GATKSAMRecord,Map<Allele,Double>> mapEntry : sample.getValue().getLikelihoodReadMap().entrySet() ) {
// only count the read if it overlaps the event, otherwise it is not added to the output read list at all
if( callLoc.overlapsP(parser.createGenomeLoc(mapEntry.getKey())) ) { // BUGBUG: This uses alignment start and stop, NOT soft start and soft end...
for( final Map.Entry<Allele,Double> alleleDoubleEntry : mapEntry.getValue().entrySet() ) {
likelihoodMap.add(mapEntry.getKey(), alleleDoubleEntry.getKey(), alleleDoubleEntry.getValue());
}
}
}
// add all filtered reads to the NO_CALL list because they weren't given any likelihoods
for( final GATKSAMRecord read : perSampleFilteredReadList.get(sample.getKey()) ) {
// only count the read if it overlaps the event, otherwise it is not added to the output read list at all
if( callLoc.overlapsP(parser.createGenomeLoc(read)) ) {
for( final Allele allele : call.getAlleles() ) {
likelihoodMap.add(read, allele, 0.0);
}
}
}
returnMap.put(sample.getKey(), likelihoodMap);
}
return returnMap;
}
protected static void cleanUpSymbolicUnassembledEvents( final List<Haplotype> haplotypes ) {
final ArrayList<Haplotype> haplotypesToRemove = new ArrayList<Haplotype>();
for( final Haplotype h : haplotypes ) {
@ -221,7 +261,41 @@ public class GenotypingEngine {
haplotypes.removeAll(haplotypesToRemove);
}
protected void mergeConsecutiveEventsBasedOnLD( final List<Haplotype> haplotypes, final TreeSet<Integer> startPosKeySet, final byte[] ref, final GenomeLoc refLoc ) {
// BUGBUG: ugh, too complicated
protected Map<String, PerReadAlleleLikelihoodMap> convertHaplotypeReadMapToAlleleReadMap( final Map<String, PerReadAlleleLikelihoodMap> haplotypeReadMap,
final Map<Allele, List<Haplotype>> alleleMapper,
final double downsamplingFraction,
final PrintStream downsamplingLog ) {
final Map<String, PerReadAlleleLikelihoodMap> alleleReadMap = new HashMap<String, PerReadAlleleLikelihoodMap>();
for( final Map.Entry<String, PerReadAlleleLikelihoodMap> haplotypeReadMapEntry : haplotypeReadMap.entrySet() ) { // for each sample
final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap = PerReadAlleleLikelihoodMap.getBestAvailablePerReadAlleleLikelihoodMap();
for( final Map.Entry<Allele, List<Haplotype>> alleleMapperEntry : alleleMapper.entrySet() ) { // for each output allele
final List<Haplotype> mappedHaplotypes = alleleMapperEntry.getValue();
for( final Map.Entry<GATKSAMRecord, Map<Allele,Double>> readEntry : haplotypeReadMapEntry.getValue().getLikelihoodReadMap().entrySet() ) { // for each read
double maxLikelihood = Double.NEGATIVE_INFINITY;
for( final Map.Entry<Allele,Double> alleleDoubleEntry : readEntry.getValue().entrySet() ) { // for each input allele
if( mappedHaplotypes.contains( new Haplotype(alleleDoubleEntry.getKey().getBases())) ) { // exact match of haplotype base string
maxLikelihood = Math.max( maxLikelihood, alleleDoubleEntry.getValue() );
}
}
perReadAlleleLikelihoodMap.add(readEntry.getKey(), alleleMapperEntry.getKey(), maxLikelihood);
}
}
perReadAlleleLikelihoodMap.performPerAlleleDownsampling(downsamplingFraction, downsamplingLog); // perform contamination downsampling
alleleReadMap.put(haplotypeReadMapEntry.getKey(), perReadAlleleLikelihoodMap);
}
return alleleReadMap;
}
protected void mergeConsecutiveEventsBasedOnLD( final List<Haplotype> haplotypes,
final List<String> samples,
final Map<String, PerReadAlleleLikelihoodMap> haplotypeReadMap,
final TreeSet<Integer> startPosKeySet,
final byte[] ref,
final GenomeLoc refLoc ) {
final int MAX_SIZE_TO_COMBINE = 15;
final double MERGE_EVENTS_R2_THRESHOLD = 0.95;
if( startPosKeySet.size() <= 1 ) { return; }
@ -265,12 +339,13 @@ public class GenotypingEngine {
}
}
// count up the co-occurrences of the events for the R^2 calculation
final ArrayList<Haplotype> haplotypeList = new ArrayList<Haplotype>();
haplotypeList.add(h);
for( final String sample : haplotypes.get(0).getSampleKeySet() ) {
for( final String sample : samples ) {
final HashSet<String> sampleSet = new HashSet<String>(1);
sampleSet.add(sample);
final double haplotypeLikelihood = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods( sampleSet, haplotypeList )[0][0];
final List<Allele> alleleList = new ArrayList<Allele>();
alleleList.add(Allele.create(h.getBases()));
final double haplotypeLikelihood = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods( sampleSet, haplotypeReadMap, alleleList )[0][0];
if( thisHapVC == null ) {
if( nextHapVC == null ) { x11 = MathUtils.approximateLog10SumLog10(x11, haplotypeLikelihood); }
else { x12 = MathUtils.approximateLog10SumLog10(x12, haplotypeLikelihood); }

35
protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
View File

@ -203,9 +203,6 @@ public class HaplotypeCaller extends ActiveRegionWalker<Integer, Integer> implem
// the genotyping engine
private GenotypingEngine genotypingEngine = null;
// the annotation engine
private VariantAnnotatorEngine annotationEngine;
// fasta reference reader to supplement the edges of the reference sequence
private CachingIndexedFastaSequenceFile referenceReader;
@ -250,7 +247,7 @@ public class HaplotypeCaller extends ActiveRegionWalker<Integer, Integer> implem
UG_engine_simple_genotyper = new UnifiedGenotyperEngine(getToolkit(), simpleUAC, logger, null, null, samples, VariantContextUtils.DEFAULT_PLOIDY);
// initialize the output VCF header
annotationEngine = new VariantAnnotatorEngine(Arrays.asList(annotationClassesToUse), annotationsToUse, annotationsToExclude, this, getToolkit());
final VariantAnnotatorEngine annotationEngine = new VariantAnnotatorEngine(Arrays.asList(annotationClassesToUse), annotationsToUse, annotationsToExclude, this, getToolkit());
Set<VCFHeaderLine> headerInfo = new HashSet<VCFHeaderLine>();
@ -283,7 +280,7 @@ public class HaplotypeCaller extends ActiveRegionWalker<Integer, Integer> implem
assemblyEngine = new SimpleDeBruijnAssembler( DEBUG, graphWriter );
likelihoodCalculationEngine = new LikelihoodCalculationEngine( (byte)gcpHMM, DEBUG, pairHMM );
genotypingEngine = new GenotypingEngine( DEBUG );
genotypingEngine = new GenotypingEngine( DEBUG, annotationEngine );
}
//---------------------------------------------------------------------------------------------------------------
@ -405,21 +402,23 @@ public class HaplotypeCaller extends ActiveRegionWalker<Integer, Integer> implem
Collections.sort( haplotypes, new Haplotype.HaplotypeBaseComparator() );
// evaluate each sample's reads against all haplotypes
final HashMap<String, ArrayList<GATKSAMRecord>> perSampleReadList = splitReadsBySample( activeRegion.getReads() );
final HashMap<String, ArrayList<GATKSAMRecord>> perSampleFilteredReadList = splitReadsBySample( filteredReads );
likelihoodCalculationEngine.computeReadLikelihoods( haplotypes, perSampleReadList );
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap = likelihoodCalculationEngine.computeReadLikelihoods( haplotypes, splitReadsBySample( activeRegion.getReads() ) );
final Map<String, ArrayList<GATKSAMRecord>> perSampleFilteredReadList = splitReadsBySample( filteredReads );
// subset down to only the best haplotypes to be genotyped in all samples ( in GGA mode use all discovered haplotypes )
final ArrayList<Haplotype> bestHaplotypes = ( UG_engine.getUAC().GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ? likelihoodCalculationEngine.selectBestHaplotypes( haplotypes ) : haplotypes );
final ArrayList<Haplotype> bestHaplotypes = ( UG_engine.getUAC().GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ? likelihoodCalculationEngine.selectBestHaplotypes( haplotypes, stratifiedReadMap ) : haplotypes );
for( final Pair<VariantContext, Map<Allele, List<Haplotype>>> callResult :
genotypingEngine.assignGenotypeLikelihoodsAndCallIndependentEvents( UG_engine, bestHaplotypes, fullReferenceWithPadding, getPaddedLoc(activeRegion), activeRegion.getLocation(), getToolkit().getGenomeLocParser(), activeAllelesToGenotype ) ) {
if( DEBUG ) { System.out.println(callResult.getFirst().toStringWithoutGenotypes()); }
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap = LikelihoodCalculationEngine.partitionReadsBasedOnLikelihoods( getToolkit().getGenomeLocParser(), perSampleReadList, perSampleFilteredReadList, callResult, UG_engine.getUAC().CONTAMINATION_FRACTION, UG_engine.getUAC().contaminationLog );
final VariantContext annotatedCall = annotationEngine.annotateContext(stratifiedReadMap, callResult.getFirst());
final Map<String, Object> myAttributes = new LinkedHashMap<String, Object>(annotatedCall.getAttributes());
vcfWriter.add( new VariantContextBuilder(annotatedCall).attributes(myAttributes).make() );
for( final VariantContext call : genotypingEngine.assignGenotypeLikelihoodsAndCallIndependentEvents( UG_engine,
bestHaplotypes,
samplesList,
stratifiedReadMap,
perSampleFilteredReadList,
fullReferenceWithPadding,
getPaddedLoc(activeRegion),
activeRegion.getLocation(),
getToolkit().getGenomeLocParser(),
activeAllelesToGenotype ) ) {
vcfWriter.add( call );
}
if( DEBUG ) { System.out.println("----------------------------------------------------------------------------------"); }
@ -474,7 +473,7 @@ public class HaplotypeCaller extends ActiveRegionWalker<Integer, Integer> implem
if( postAdapterRead != null && !postAdapterRead.isEmpty() && postAdapterRead.getCigar().getReadLength() > 0 ) {
final GATKSAMRecord clippedRead = ReadClipper.hardClipLowQualEnds( postAdapterRead, MIN_TAIL_QUALITY );
// protect against INTERVALS with abnormally high coverage
// BUGBUG: remove when positinal downsampler is hooked up to ART/HC
// BUGBUG: remove when positional downsampler is hooked up to ART/HC
if( clippedRead.getReadLength() > 0 && activeRegion.size() < samplesList.size() * DOWNSAMPLE_PER_SAMPLE_PER_REGION ) {
activeRegion.add(clippedRead);
}

156
protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java
View File

@ -71,8 +71,9 @@ public class LikelihoodCalculationEngine {
DEBUG = debug;
}
public void computeReadLikelihoods( final ArrayList<Haplotype> haplotypes, final HashMap<String, ArrayList<GATKSAMRecord>> perSampleReadList ) {
public Map<String, PerReadAlleleLikelihoodMap> computeReadLikelihoods( final ArrayList<Haplotype> haplotypes, final HashMap<String, ArrayList<GATKSAMRecord>> perSampleReadList ) {
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap = new HashMap<String, PerReadAlleleLikelihoodMap>();
int X_METRIC_LENGTH = 0;
for( final Map.Entry<String, ArrayList<GATKSAMRecord>> sample : perSampleReadList.entrySet() ) {
for( final GATKSAMRecord read : sample.getValue() ) {
@ -97,20 +98,16 @@ public class LikelihoodCalculationEngine {
for( final Map.Entry<String, ArrayList<GATKSAMRecord>> sampleEntry : perSampleReadList.entrySet() ) {
//if( DEBUG ) { System.out.println("Evaluating sample " + sample + " with " + perSampleReadList.get( sample ).size() + " passing reads"); }
// evaluate the likelihood of the reads given those haplotypes
computeReadLikelihoods( haplotypes, sampleEntry.getValue(), sampleEntry.getKey() );
stratifiedReadMap.put(sampleEntry.getKey(), computeReadLikelihoods(haplotypes, sampleEntry.getValue(), sampleEntry.getKey()));
}
return stratifiedReadMap;
}
private void computeReadLikelihoods( final ArrayList<Haplotype> haplotypes, final ArrayList<GATKSAMRecord> reads, final String sample ) {
private PerReadAlleleLikelihoodMap computeReadLikelihoods( final ArrayList<Haplotype> haplotypes, final ArrayList<GATKSAMRecord> reads, final String sample ) {
final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap = PerReadAlleleLikelihoodMap.getBestAvailablePerReadAlleleLikelihoodMap();
final int numHaplotypes = haplotypes.size();
final int numReads = reads.size();
final double[][] readLikelihoods = new double[numHaplotypes][numReads];
final int[][] readCounts = new int[numHaplotypes][numReads];
for( int iii = 0; iii < numReads; iii++ ) {
final GATKSAMRecord read = reads.get(iii);
final int readCount = ReadUtils.getMeanRepresentativeReadCount(read);
for( final GATKSAMRecord read : reads ) {
final byte[] overallGCP = new byte[read.getReadLength()];
Arrays.fill( overallGCP, constantGCP ); // Is there a way to derive empirical estimates for this from the data?
Haplotype previousHaplotypeSeen = null;
@ -129,14 +126,12 @@ public class LikelihoodCalculationEngine {
final int haplotypeStart = ( previousHaplotypeSeen == null ? 0 : computeFirstDifferingPosition(haplotype.getBases(), previousHaplotypeSeen.getBases()) );
previousHaplotypeSeen = haplotype;
readLikelihoods[jjj][iii] = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), read.getReadBases(),
readQuals, readInsQuals, readDelQuals, overallGCP, haplotypeStart, jjj == 0);
readCounts[jjj][iii] = readCount;
perReadAlleleLikelihoodMap.add(read, Allele.create(haplotype.getBases()),
pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), read.getReadBases(),
readQuals, readInsQuals, readDelQuals, overallGCP, haplotypeStart, jjj == 0));
}
}
for( int jjj = 0; jjj < numHaplotypes; jjj++ ) {
haplotypes.get(jjj).addReadLikelihoods( sample, readLikelihoods[jjj], readCounts[jjj] );
}
return perReadAlleleLikelihoodMap;
}
private static int computeFirstDifferingPosition( final byte[] b1, final byte[] b2 ) {
@ -148,19 +143,21 @@ public class LikelihoodCalculationEngine {
return Math.min(b1.length, b2.length);
}
// This function takes just a single sample and a haplotypeMapping
@Requires({"haplotypeMapping.size() > 0"})
@Ensures({"result.length == result[0].length", "result.length == haplotypeMapping.size()"})
public static double[][] computeDiploidHaplotypeLikelihoods( final String sample, final Map<Allele, List<Haplotype>> haplotypeMapping, final List<Allele> alleleOrdering ) {
@Requires({"alleleOrdering.size() > 0"})
@Ensures({"result.length == result[0].length", "result.length == alleleOrdering.size()"})
public static double[][] computeDiploidHaplotypeLikelihoods( final String sample,
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap,
final List<Allele> alleleOrdering ) {
final TreeSet<String> sampleSet = new TreeSet<String>();
sampleSet.add(sample);
return computeDiploidHaplotypeLikelihoods(sampleSet, haplotypeMapping, alleleOrdering);
return computeDiploidHaplotypeLikelihoods(sampleSet, stratifiedReadMap, alleleOrdering);
}
// This function takes a set of samples to pool over and a haplotypeMapping
@Requires({"haplotypeMapping.size() > 0"})
@Ensures({"result.length == result[0].length", "result.length == haplotypeMapping.size()"})
public static double[][] computeDiploidHaplotypeLikelihoods( final Set<String> samples, final Map<Allele, List<Haplotype>> haplotypeMapping, final List<Allele> alleleOrdering ) {
@Requires({"alleleOrdering.size() > 0"})
@Ensures({"result.length == result[0].length", "result.length == alleleOrdering.size()"})
public static double[][] computeDiploidHaplotypeLikelihoods( final Set<String> samples,
final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap,
final List<Allele> alleleOrdering ) {
final int numHaplotypes = alleleOrdering.size();
final double[][] haplotypeLikelihoodMatrix = new double[numHaplotypes][numHaplotypes];
@ -170,59 +167,19 @@ public class LikelihoodCalculationEngine {
// compute the diploid haplotype likelihoods
for( int iii = 0; iii < numHaplotypes; iii++ ) {
final Allele iii_allele = alleleOrdering.get(iii);
for( int jjj = 0; jjj <= iii; jjj++ ) {
for( final Haplotype iii_mapped : haplotypeMapping.get(alleleOrdering.get(iii)) ) {
for( final Haplotype jjj_mapped : haplotypeMapping.get(alleleOrdering.get(jjj)) ) {
double haplotypeLikelihood = 0.0;
for( final String sample : samples ) {
final double[] readLikelihoods_iii = iii_mapped.getReadLikelihoods(sample);
final int[] readCounts_iii = iii_mapped.getReadCounts(sample);
final double[] readLikelihoods_jjj = jjj_mapped.getReadLikelihoods(sample);
for( int kkk = 0; kkk < readLikelihoods_iii.length; kkk++ ) {
// Compute log10(10^x1/2 + 10^x2/2) = log10(10^x1+10^x2)-log10(2)
// First term is approximated by Jacobian log with table lookup.
haplotypeLikelihood += readCounts_iii[kkk] * ( MathUtils.approximateLog10SumLog10(readLikelihoods_iii[kkk], readLikelihoods_jjj[kkk]) + LOG_ONE_HALF );
}
}
haplotypeLikelihoodMatrix[iii][jjj] = Math.max(haplotypeLikelihoodMatrix[iii][jjj], haplotypeLikelihood);
}
}
}
}
// normalize the diploid likelihoods matrix
return normalizeDiploidLikelihoodMatrixFromLog10( haplotypeLikelihoodMatrix );
}
// This function takes a set of samples to pool over and a haplotypeMapping
@Requires({"haplotypeList.size() > 0"})
@Ensures({"result.length == result[0].length", "result.length == haplotypeList.size()"})
public static double[][] computeDiploidHaplotypeLikelihoods( final Set<String> samples, final List<Haplotype> haplotypeList ) {
final int numHaplotypes = haplotypeList.size();
final double[][] haplotypeLikelihoodMatrix = new double[numHaplotypes][numHaplotypes];
for( int iii = 0; iii < numHaplotypes; iii++ ) {
Arrays.fill(haplotypeLikelihoodMatrix[iii], Double.NEGATIVE_INFINITY);
}
// compute the diploid haplotype likelihoods
// todo - needs to be generalized to arbitrary ploidy, cleaned and merged with PairHMMIndelErrorModel code
for( int iii = 0; iii < numHaplotypes; iii++ ) {
final Haplotype iii_haplotype = haplotypeList.get(iii);
for( int jjj = 0; jjj <= iii; jjj++ ) {
final Haplotype jjj_haplotype = haplotypeList.get(jjj);
final Allele jjj_allele = alleleOrdering.get(jjj);
double haplotypeLikelihood = 0.0;
for( final String sample : samples ) {
final double[] readLikelihoods_iii = iii_haplotype.getReadLikelihoods(sample);
final int[] readCounts_iii = iii_haplotype.getReadCounts(sample);
final double[] readLikelihoods_jjj = jjj_haplotype.getReadLikelihoods(sample);
for( int kkk = 0; kkk < readLikelihoods_iii.length; kkk++ ) {
for( final Map.Entry<GATKSAMRecord, Map<Allele,Double>> entry : stratifiedReadMap.get(sample).getLikelihoodReadMap().entrySet() ) {
// Compute log10(10^x1/2 + 10^x2/2) = log10(10^x1+10^x2)-log10(2)
// First term is approximated by Jacobian log with table lookup.
haplotypeLikelihood += readCounts_iii[kkk] * ( MathUtils.approximateLog10SumLog10(readLikelihoods_iii[kkk], readLikelihoods_jjj[kkk]) + LOG_ONE_HALF );
haplotypeLikelihood += ReadUtils.getMeanRepresentativeReadCount( entry.getKey() ) *
( MathUtils.approximateLog10SumLog10(entry.getValue().get(iii_allele), entry.getValue().get(jjj_allele)) + LOG_ONE_HALF );
}
}
haplotypeLikelihoodMatrix[iii][jjj] = Math.max(haplotypeLikelihoodMatrix[iii][jjj], haplotypeLikelihood);
haplotypeLikelihoodMatrix[iii][jjj] = haplotypeLikelihood;
}
}
@ -312,13 +269,16 @@ public class LikelihoodCalculationEngine {
@Requires({"haplotypes.size() > 0"})
@Ensures({"result.size() <= haplotypes.size()"})
public ArrayList<Haplotype> selectBestHaplotypes( final ArrayList<Haplotype> haplotypes ) {
public ArrayList<Haplotype> selectBestHaplotypes( final ArrayList<Haplotype> haplotypes, final Map<String, PerReadAlleleLikelihoodMap> stratifiedReadMap ) {
final int numHaplotypes = haplotypes.size();
final Set<String> sampleKeySet = haplotypes.get(0).getSampleKeySet(); // BUGBUG: assume all haplotypes saw the same samples
final Set<String> sampleKeySet = stratifiedReadMap.keySet();
final ArrayList<Integer> bestHaplotypesIndexList = new ArrayList<Integer>();
bestHaplotypesIndexList.add( findReferenceIndex(haplotypes) ); // always start with the reference haplotype
final double[][] haplotypeLikelihoodMatrix = computeDiploidHaplotypeLikelihoods( sampleKeySet, haplotypes ); // all samples pooled together
final List<Allele> haplotypesAsAlleles = new ArrayList<Allele>();
for( final Haplotype h : haplotypes ) { haplotypesAsAlleles.add(Allele.create(h.getBases())); }
final double[][] haplotypeLikelihoodMatrix = computeDiploidHaplotypeLikelihoods( sampleKeySet, stratifiedReadMap, haplotypesAsAlleles ); // all samples pooled together
int hap1 = 0;
int hap2 = 0;
@ -358,52 +318,4 @@ public class LikelihoodCalculationEngine {
}
throw new ReviewedStingException( "No reference haplotype found in the list of haplotypes!" );
}
public static Map<String, PerReadAlleleLikelihoodMap> partitionReadsBasedOnLikelihoods( final GenomeLocParser parser,
final HashMap<String, ArrayList<GATKSAMRecord>> perSampleReadList,
final HashMap<String, ArrayList<GATKSAMRecord>> perSampleFilteredReadList,
final Pair<VariantContext, Map<Allele,List<Haplotype>>> call,
final double downsamplingFraction,
final PrintStream downsamplingLog ) {
final Map<String, PerReadAlleleLikelihoodMap> returnMap = new HashMap<String, PerReadAlleleLikelihoodMap>();
final GenomeLoc callLoc = parser.createGenomeLoc(call.getFirst());
for( final Map.Entry<String, ArrayList<GATKSAMRecord>> sample : perSampleReadList.entrySet() ) {
final PerReadAlleleLikelihoodMap likelihoodMap = PerReadAlleleLikelihoodMap.getBestAvailablePerReadAlleleLikelihoodMap();
final ArrayList<GATKSAMRecord> readsForThisSample = sample.getValue();
for( int iii = 0; iii < readsForThisSample.size(); iii++ ) {
final GATKSAMRecord read = readsForThisSample.get(iii); // BUGBUG: assumes read order in this list and haplotype likelihood list are the same!
// only count the read if it overlaps the event, otherwise it is not added to the output read list at all
if( callLoc.overlapsP(parser.createGenomeLoc(read)) ) {
for( final Allele a : call.getFirst().getAlleles() ) {
double maxLikelihood = Double.NEGATIVE_INFINITY;
for( final Haplotype h : call.getSecond().get(a) ) { // use the max likelihood from all the haplotypes which mapped to this allele (achieved via the haplotype mapper object)
final double likelihood = h.getReadLikelihoods(sample.getKey())[iii];
if( likelihood > maxLikelihood ) {
maxLikelihood = likelihood;
}
}
likelihoodMap.add(read, a, maxLikelihood);
}
}
}
// down-sample before adding filtered reads
likelihoodMap.performPerAlleleDownsampling(downsamplingFraction, downsamplingLog);
// add all filtered reads to the NO_CALL list because they weren't given any likelihoods
for( final GATKSAMRecord read : perSampleFilteredReadList.get(sample.getKey()) ) {
// only count the read if it overlaps the event, otherwise it is not added to the output read list at all
if( callLoc.overlapsP(parser.createGenomeLoc(read)) ) {
for( final Allele a : call.getFirst().getAlleles() ) {
likelihoodMap.add(read, a, 0.0);
}
}
}
returnMap.put(sample.getKey(), likelihoodMap);
}
return returnMap;
}
}

2
protected/java/test/org/broadinstitute/sting/gatk/walkers/bqsr/BQSRIntegrationTest.java
View File

@ -37,6 +37,7 @@ public class BQSRIntegrationTest extends WalkerTest {
" -L " + interval +
args +
" -knownSites " + (reference.equals(b36KGReference) ? b36dbSNP129 : hg18dbSNP132) +
" --allow_potentially_misencoded_quality_scores" + // TODO -- remove me when we get new SOLiD bams
" -o %s";
}
@ -112,6 +113,7 @@ public class BQSRIntegrationTest extends WalkerTest {
" -R " + b36KGReference +
" -I " + privateTestDir + "NA19240.chr1.BFAST.SOLID.hasCSNoCall.bam" +
" -L 1:50,000-80,000" +
" --allow_potentially_misencoded_quality_scores" + // TODO -- remove me when we get new SOLiD bams
" -o %s",
1, // just one output file
UserException.class);

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

@ -177,7 +177,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
executeTest("test using comp track", spec);
}
@Test
@Test(enabled = false) // EB: for some reason this test crashes whenever I run it on my local machine
public void testNoCmdLineHeaderStdout() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
baseCommandNoCmdLineHeaderStdout + " -glm INDEL -L 1:67,225,396-67,288,518", 0,
@ -351,7 +351,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + result.get(0).getAbsolutePath() + " -I " + validationDataLocation +
"low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1,
Arrays.asList("c526c234947482d1cd2ffc5102083a08"));
Arrays.asList("1256a7eceff2c2374c231ff981df486d"));
executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2);
}
@ -436,8 +436,8 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
@Test
public void testNsInCigar() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -I " + validationDataLocation + "testWithNs.bam -o %s -L 8:141799600-141814700", 1,
Arrays.asList("d6d40bacd540a41f305420dfea35e04a"));
"-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -I " + validationDataLocation + "testWithNs.bam -o %s -L 8:141813600-141813700 -out_mode EMIT_ALL_SITES", 1,
Arrays.asList("32f18ba50406cd8c8069ba07f2f89558"));
executeTest("test calling on reads with Ns in CIGAR", spec);
}
@ -457,7 +457,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
@Test
public void testReducedBamSNPs() {
testReducedCalling("SNP", "f5ccbc96d0d66832dd9b3c5cb6507db4");
testReducedCalling("SNP", "dee6590e3b7079890bc3a9cb372c297e");
}
@Test

22
protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
View File

@ -21,19 +21,19 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerMultiSample() {
HCTest(CEUTRIO_BAM, "", "2b39732ff8e0de5bc2ae949aaf7a6f21");
HCTest(CEUTRIO_BAM, "", "d602d40852ad6d2d094be07e60cf95bd");
}
@Test
public void testHaplotypeCallerSingleSample() {
HCTest(NA12878_BAM, "", "8b217638ff585effb9cc70e9a9aa544f");
HCTest(NA12878_BAM, "", "70ad9d53dda4d302b879ca2b7dd5b368");
}
// TODO -- add more tests for GGA mode, especially with input alleles that are complex variants and/or not trimmed
@Test
public void testHaplotypeCallerMultiSampleGGA() {
HCTest(CEUTRIO_BAM, "--max_alternate_alleles 3 -gt_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf",
"541aa8291f03ba33bd1ad3d731fd5657");
HCTest(CEUTRIO_BAM, "--max_alternate_alleles 3 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf",
"e2b3bf420c45c677956a2e4a56d75ea2");
}
private void HCTestComplexVariants(String bam, String args, String md5) {
@ -44,18 +44,18 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerMultiSampleComplex() {
HCTestComplexVariants(privateTestDir + "AFR.complex.variants.bam", "", "fd7170cbde7df04d4fbe1da7903c31c6");
HCTestComplexVariants(privateTestDir + "AFR.complex.variants.bam", "", "883871f8bb4099f69fd804f8a6181954");
}
private void HCTestSymbolicVariants(String bam, String args, String md5) {
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, bam) + " -L 20:5947969-5948369 -L 20:61091236-61091636 --no_cmdline_in_header -o %s -minPruning 2";
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, bam) + " -L 20:5947969-5948369 -L 20:61091236-61091636 --no_cmdline_in_header -o %s -minPruning 1";
final WalkerTestSpec spec = new WalkerTestSpec(base + " " + args, Arrays.asList(md5));
executeTest("testHaplotypeCallerSymbolicVariants: args=" + args, spec);
}
@Test
public void testHaplotypeCallerSingleSampleSymbolic() {
HCTestSymbolicVariants(NA12878_CHR20_BAM, "", "99456fc7207c1fe9f367a0d0afae87cd");
HCTestSymbolicVariants(NA12878_CHR20_BAM, "", "338ab3b7dc3d54df8af94c0811028a75");
}
private void HCTestIndelQualityScores(String bam, String args, String md5) {
@ -66,20 +66,20 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerSingleSampleIndelQualityScores() {
HCTestIndelQualityScores(NA12878_RECALIBRATED_BAM, "", "6c1631785b3f832aecab1a99f0454762");
HCTestIndelQualityScores(NA12878_RECALIBRATED_BAM, "", "aff11b014ca42bfa301bcced5f5e54dd");
}
@Test
public void HCTestProblematicReadsModifiedInActiveRegions() {
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, privateTestDir + "haplotype-problem-4.bam") + " --no_cmdline_in_header -o %s -minPruning 3 -L 4:49139026-49139965";
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("ec437d2d9f3ae07d155983be0155c8ed"));
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("2f4ed6dc969bee041215944a9b24328f"));
executeTest("HCTestProblematicReadsModifiedInActiveRegions: ", spec);
}
@Test
public void HCTestStructuralIndels() {
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, privateTestDir + "AFR.structural.indels.bam") + " --no_cmdline_in_header -o %s -minPruning 6 -L 20:8187565-8187800 -L 20:18670537-18670730";
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("237601bbc39694c7413a332cbb656c8e"));
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("d8d6f2ebe79bca81c8a0911daa153b89"));
executeTest("HCTestStructuralIndels: ", spec);
}
@ -93,7 +93,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
public void HCTestReducedBam() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller -R " + b37KGReference + " --no_cmdline_in_header -I " + privateTestDir + "bamExample.ReducedRead.ADAnnotation.bam -o %s -L 1:67,225,396-67,288,518", 1,
Arrays.asList("40bf739fb2b1743642498efe79ea6342"));
Arrays.asList("d01cb5f77ed5aca1d228cfbce9364c21"));
executeTest("HC calling on a ReducedRead BAM", spec);
}
}

7
protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngineUnitTest.java
View File

@ -51,6 +51,8 @@ public class LikelihoodCalculationEngineUnitTest extends BaseTest {
Assert.assertTrue(compareDoubleArrays(LikelihoodCalculationEngine.normalizeDiploidLikelihoodMatrixFromLog10(likelihoodMatrix2), normalizedMatrix2));
}
// BUGBUG: LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods has changed! Need to make new unit tests!
/*
private class BasicLikelihoodTestProvider extends TestDataProvider {
public Double readLikelihoodForHaplotype1;
public Double readLikelihoodForHaplotype2;
@ -152,10 +154,9 @@ public class LikelihoodCalculationEngineUnitTest extends BaseTest {
logger.warn(String.format("Test: %s", cfg.toString()));
Assert.assertTrue(compareDoubleArrays(calculatedMatrix, expectedMatrix));
}
*/
/**
* Private function to compare 2d arrays
*/
//Private function to compare 2d arrays
private boolean compareDoubleArrays(double[][] b1, double[][] b2) {
if( b1.length != b2.length ) {
return false; // sanity check

16
public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java
View File

@ -206,6 +206,22 @@ public class GATKArgumentCollection {
@Argument(fullName = "baqGapOpenPenalty", shortName="baqGOP", doc="BAQ gap open penalty (Phred Scaled). Default value is 40. 30 is perhaps better for whole genome call sets", required = false)
public double BAQGOP = BAQ.DEFAULT_GOP;
// --------------------------------------------------------------------------------------------------------------
//
// quality encoding checking arguments
//
// --------------------------------------------------------------------------------------------------------------
/**
* Q0 == ASCII 33 according to the SAM specification, whereas Illumina encoding starts at Q64. The idea here is
* simple: we just iterate over all reads and subtract 31 from every quality score.
*/
@Argument(fullName = "fix_misencoded_quality_scores", shortName="fixMisencodedQuals", doc="Fix mis-encoded base quality scores", required = false)
public boolean FIX_MISENCODED_QUALS = false;
@Argument(fullName = "allow_potentially_misencoded_quality_scores", shortName="allowPotentiallyMisencodedQuals", doc="Do not fail when encountered base qualities that are too high and seemingly indicate a problem with the base quality encoding of the BAM file", required = false)
public boolean ALLOW_POTENTIALLY_MISENCODED_QUALS = false;
// --------------------------------------------------------------------------------------------------------------
//
// performance log arguments

8
public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java
View File

@ -30,12 +30,10 @@ import net.sf.samtools.*;
import net.sf.samtools.util.CloseableIterator;
import net.sf.samtools.util.RuntimeIOException;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.gatk.downsampling.*;
import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
import org.broadinstitute.sting.gatk.downsampling.DownsamplingMethod;
import org.broadinstitute.sting.gatk.ReadMetrics;
import org.broadinstitute.sting.gatk.ReadProperties;
import org.broadinstitute.sting.gatk.arguments.ValidationExclusion;
import org.broadinstitute.sting.gatk.downsampling.*;
import org.broadinstitute.sting.gatk.filters.CountingFilteringIterator;
import org.broadinstitute.sting.gatk.filters.ReadFilter;
import org.broadinstitute.sting.gatk.iterators.*;
@ -567,7 +565,7 @@ public class SAMDataSource {
*
* @return the start positions of the first chunk of reads for all BAM files
*/
public Map<SAMReaderID, GATKBAMFileSpan> getInitialReaderPositions() {
protected Map<SAMReaderID, GATKBAMFileSpan> getInitialReaderPositions() {
Map<SAMReaderID, GATKBAMFileSpan> initialPositions = new HashMap<SAMReaderID, GATKBAMFileSpan>();
SAMReaders readers = resourcePool.getAvailableReaders();
@ -585,7 +583,7 @@ public class SAMDataSource {
* @param shard The shard specifying the data limits.
* @return An iterator over the selected data.
*/
public StingSAMIterator getIterator( Shard shard ) {
protected StingSAMIterator getIterator( Shard shard ) {
return getIterator(resourcePool.getAvailableReaders(), shard, shard instanceof ReadShard);
}

4
public/java/src/org/broadinstitute/sting/gatk/iterators/ReadTransformer.java
View File

@ -41,7 +41,7 @@ abstract public class ReadTransformer {
protected ReadTransformer() {}
/**
* Master initialization routine. Called to setup a ReadTransform, using it's overloaded initialialSub routine.
* Master initialization routine. Called to setup a ReadTransform, using it's overloaded initializeSub routine.
*
* @param overrideTime if not null, we will run this ReadTransform at the time provided, regardless of the timing of this read transformer itself
* @param engine the engine, for initializing values
@ -59,7 +59,7 @@ abstract public class ReadTransformer {
}
/**
* Subclasses must override this to initialize themeselves
* Subclasses must override this to initialize themselves
*
* @param engine the engine, for initializing values
* @param walker the walker we intend to run

2
public/java/src/org/broadinstitute/sting/gatk/report/GATKReport.java
View File

@ -343,7 +343,7 @@ public class GATKReport {
GATKReportTable table = tables.firstEntry().getValue();
if ( table.getNumColumns() != values.length )
throw new ReviewedStingException("The number of arguments in writeRow() must match the number of columns in the table");
throw new ReviewedStingException("The number of arguments in writeRow() " + values.length + " must match the number of columns in the table" + table.getNumColumns());
final int rowIndex = table.getNumRows();
for ( int i = 0; i < values.length; i++ )

6
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java
View File

@ -276,6 +276,12 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
for (PileupElement p : sample.getValue().getBasePileup()) {
// ignore reduced reads because they are always on the forward strand!
// TODO -- when het compression is enabled in RR, we somehow need to allow those reads through into the Fisher test
if ( p.getRead().isReducedRead() )
continue;
if ( ! RankSumTest.isUsableBase(p, false) ) // ignore deletions
continue;

1
public/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
View File

@ -329,7 +329,6 @@ public class PairHMMIndelErrorModel {
getContextHomopolymerLength(readBases,hrunProfile);
fillGapProbabilities(hrunProfile, contextLogGapOpenProbabilities, contextLogGapContinuationProbabilities);
for (Allele a: haplotypeMap.keySet()) {
Haplotype haplotype = haplotypeMap.get(a);

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

@ -81,7 +81,7 @@ public class VariantDataManager {
final double theSTD = standardDeviation(theMean, iii);
logger.info( annotationKeys.get(iii) + String.format(": \t mean = %.2f\t standard deviation = %.2f", theMean, theSTD) );
if( Double.isNaN(theMean) ) {
throw new UserException.BadInput("Values for " + annotationKeys.get(iii) + " annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations. See " + HelpUtils.GATK_FORUM_URL + "discussion/49/using-variant-annotator");
throw new UserException.BadInput("Values for " + annotationKeys.get(iii) + " annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations. See " + HelpUtils.forumPost("discussion/49/using-variant-annotator"));
}
foundZeroVarianceAnnotation = foundZeroVarianceAnnotation || (theSTD < 1E-6);

18
public/java/src/org/broadinstitute/sting/jna/drmaa/v1_0/JnaSession.java
View File

@ -210,13 +210,23 @@ public class JnaSession implements Session {
}
public static void setAttribute(Pointer jt, String name, String value) throws DrmaaException {
checkError(LibDrmaa.drmaa_set_attribute(jt, name, value, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN));
if (getAttrNames().contains(name)) {
checkError(LibDrmaa.drmaa_set_attribute(jt, name, value, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN));
}
else {
throw new InvalidAttributeValueException("Attribute " + name + " is not supported by this implementation of DRMAA");
}
}
public static String getAttribute(Pointer jt, String name) throws DrmaaException {
Memory attrBuffer = new Memory(LibDrmaa.DRMAA_ATTR_BUFFER);
checkError(LibDrmaa.drmaa_get_attribute(jt, name, attrBuffer, LibDrmaa.DRMAA_ATTR_BUFFER_LEN, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN));
return attrBuffer.getString(0);
if (getAttrNames().contains(name)) {
Memory attrBuffer = new Memory(LibDrmaa.DRMAA_ATTR_BUFFER);
checkError(LibDrmaa.drmaa_get_attribute(jt, name, attrBuffer, LibDrmaa.DRMAA_ATTR_BUFFER_LEN, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN));
return attrBuffer.getString(0);
}
else {
throw new InvalidAttributeValueException("Attribute " + name + " is not supported by this implementation of DRMAA");
}
}
public static void setVectorAttribute(Pointer jt, String name, Collection<String> values) throws DrmaaException {

27
public/java/src/org/broadinstitute/sting/utils/Haplotype.java
View File

@ -41,8 +41,6 @@ public class Haplotype {
protected final byte[] bases;
protected final double[] quals;
private GenomeLoc genomeLocation = null;
private HashMap<String, double[]> readLikelihoodsPerSample = null;
private HashMap<String, int[]> readCountsPerSample = null;
private HashMap<Integer, VariantContext> eventMap = null;
private boolean isRef = false;
private Cigar cigar;
@ -94,31 +92,6 @@ public class Haplotype {
return Arrays.hashCode(bases);
}
public void addReadLikelihoods( final String sample, final double[] readLikelihoods, final int[] readCounts ) {
if( readLikelihoodsPerSample == null ) {
readLikelihoodsPerSample = new HashMap<String, double[]>();
}
readLikelihoodsPerSample.put(sample, readLikelihoods);
if( readCountsPerSample == null ) {
readCountsPerSample = new HashMap<String, int[]>();
}
readCountsPerSample.put(sample, readCounts);
}
@Ensures({"result != null"})
public double[] getReadLikelihoods( final String sample ) {
return readLikelihoodsPerSample.get(sample);
}
@Ensures({"result != null"})
public int[] getReadCounts( final String sample ) {
return readCountsPerSample.get(sample);
}
public Set<String> getSampleKeySet() {
return readLikelihoodsPerSample.keySet();
}
public HashMap<Integer, VariantContext> getEventMap() {
return eventMap;
}

2
public/java/src/org/broadinstitute/sting/utils/QualityUtils.java
View File

@ -14,7 +14,7 @@ public class QualityUtils {
public final static double ERROR_RATE_OF_MAX_QUAL_SCORE = qualToErrorProbRaw(MAX_QUAL_SCORE);
public final static double MIN_REASONABLE_ERROR = 0.0001;
public final static byte MAX_REASONABLE_Q_SCORE = 40;
public final static byte MAX_REASONABLE_Q_SCORE = 60; // quals above this value are extremely suspicious
public final static byte MIN_USABLE_Q_SCORE = 6;
public final static int MAPPING_QUALITY_UNAVAILABLE = 255;

2
public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java
View File

@ -414,7 +414,7 @@ public class BAQ {
throw new ReviewedStingException("BAQ tag calculation error. BAQ value above base quality at " + read);
// the original quality is too high, almost certainly due to using the wrong encoding in the BAM file
if ( tag > Byte.MAX_VALUE )
throw new UserException.MalformedBAM(read, "we encountered an extremely high quality score (" + (bq - 64) + ") with BAQ correction factor of " + baq_i + "; the BAM file appears to be using the wrong encoding for quality scores");
throw new UserException.MisencodedBAM(read, "we encountered an extremely high quality score (" + (int)read.getBaseQualities()[i] + ") with BAQ correction factor of " + baq_i);
bqTag[i] = (byte)tag;
}
return new String(bqTag);

14
public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java
View File

@ -240,6 +240,16 @@ public class UserException extends ReviewedStingException {
}
}
public static class MisencodedBAM extends UserException {
public MisencodedBAM(SAMRecord read, String message) {
this(read.getFileSource() != null ? read.getFileSource().getReader().toString() : "(none)", message);
}
public MisencodedBAM(String source, String message) {
super(String.format("SAM/BAM file %s appears to be using the wrong encoding for quality scores: %s; please see the GATK --help documentation for options related to this error", source, message));
}
}
public static class MalformedVCF extends UserException {
public MalformedVCF(String message, String line) {
super(String.format("The provided VCF file is malformed at line %s: %s", line, message));
@ -268,7 +278,7 @@ public class UserException extends ReviewedStingException {
public static class ReadMissingReadGroup extends MalformedBAM {
public ReadMissingReadGroup(SAMRecord read) {
super(read, String.format("Read %s is either missing the read group or its read group is not defined in the BAM header, both of which are required by the GATK. Please use " + HelpUtils.GATK_FORUM_URL + "discussion/59/companion-utilities-replacereadgroups to fix this problem", read.getReadName()));
super(read, String.format("Read %s is either missing the read group or its read group is not defined in the BAM header, both of which are required by the GATK. Please use " + HelpUtils.forumPost("discussion/59/companion-utilities-replacereadgroups to fix this problem"), read.getReadName()));
}
}
@ -344,7 +354,7 @@ public class UserException extends ReviewedStingException {
super(String.format("Lexicographically sorted human genome sequence detected in %s."
+ "\nFor safety's sake the GATK requires human contigs in karyotypic order: 1, 2, ..., 10, 11, ..., 20, 21, 22, X, Y with M either leading or trailing these contigs."
+ "\nThis is because all distributed GATK resources are sorted in karyotypic order, and your processing will fail when you need to use these files."
+ "\nYou can use the ReorderSam utility to fix this problem: " + HelpUtils.GATK_FORUM_URL + "discussion/58/companion-utilities-reordersam"
+ "\nYou can use the ReorderSam utility to fix this problem: " + HelpUtils.forumPost("discussion/58/companion-utilities-reordersam")
+ "\n %s contigs = %s",
name, name, ReadUtils.prettyPrintSequenceRecords(dict)));
}

8
public/java/src/org/broadinstitute/sting/utils/genotyper/PerReadAlleleLikelihoodMap.java
View File

@ -38,10 +38,10 @@ import java.util.*;
public abstract class PerReadAlleleLikelihoodMap {
public static final double INDEL_LIKELIHOOD_THRESH = 0.1;
public static final double INFORMATIVE_LIKELIHOOD_THRESHOLD = 0.1;
protected List<Allele> alleles;
protected Map<GATKSAMRecord,Map<Allele,Double>> likelihoodReadMap;
protected Map<GATKSAMRecord, Map<Allele, Double>> likelihoodReadMap;
public abstract void performPerAlleleDownsampling(final double downsamplingFraction, final PrintStream log);
public abstract ReadBackedPileup createPerAlleleDownsampledBasePileup(final ReadBackedPileup pileup, final double downsamplingFraction, final PrintStream log);
@ -68,7 +68,7 @@ public abstract class PerReadAlleleLikelihoodMap {
}
public void add(PileupElement p, Allele a, Double likelihood) {
add(p.getRead(),a,likelihood);
add(p.getRead(), a, likelihood);
}
public boolean containsPileupElement(PileupElement p) {
@ -120,7 +120,7 @@ public abstract class PerReadAlleleLikelihoodMap {
prevMaxLike = el.getValue();
}
}
return (maxLike - prevMaxLike > INDEL_LIKELIHOOD_THRESH ? mostLikelyAllele : Allele.NO_CALL );
return (maxLike - prevMaxLike > INFORMATIVE_LIKELIHOOD_THRESHOLD ? mostLikelyAllele : Allele.NO_CALL );
}
public static PerReadAlleleLikelihoodMap getBestAvailablePerReadAlleleLikelihoodMap() {

5
public/java/src/org/broadinstitute/sting/utils/help/HelpUtils.java
View File

@ -38,8 +38,9 @@ public class HelpUtils {
public final static String GATK_FORUM_URL = "http://gatkforums.broadinstitute.org/";
public final static String GATK_FORUM_API_URL = "https://gatkforums.broadinstitute.org/api/v1/";
public static String forumPost(String post) {
return GATK_FORUM_URL + post;
}
protected static boolean assignableToClass(ProgramElementDoc classDoc, Class lhsClass, boolean requireConcrete) {
try {

19
public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java
View File

@ -652,23 +652,19 @@ public abstract class AbstractReadBackedPileup<RBP extends AbstractReadBackedPil
PerSamplePileupElementTracker<PE> tracker = (PerSamplePileupElementTracker<PE>) pileupElementTracker;
PerSamplePileupElementTracker<PE> filteredTracker = new PerSamplePileupElementTracker<PE>();
int current = 0;
for (final String sample : tracker.getSamples()) {
PileupElementTracker<PE> perSampleElements = tracker.getElements(sample);
List<PileupElement> filteredPileup = new ArrayList<PileupElement>();
for (PileupElement p : perSampleElements) {
int current = 0;
UnifiedPileupElementTracker<PE> filteredPileup = new UnifiedPileupElementTracker<PE>();
for (PE p : perSampleElements) {
if (positions.contains(current))
filteredPileup.add(p);
}
current++;
if (!filteredPileup.isEmpty()) {
AbstractReadBackedPileup<RBP, PE> pileup = createNewPileup(loc, perSampleElements);
filteredTracker.addElements(sample, pileup.pileupElementTracker);
}
current++;
filteredTracker.addElements(sample, filteredPileup);
}
return (RBP) createNewPileup(loc, filteredTracker);
@ -1058,6 +1054,11 @@ public abstract class AbstractReadBackedPileup<RBP extends AbstractReadBackedPil
public FragmentCollection<PileupElement> toFragments() {
return FragmentUtils.create(this);
}
@Override
public ReadBackedPileup copy() {
return new ReadBackedPileupImpl(loc, (PileupElementTracker<PileupElement>) pileupElementTracker.copy());
}
}

17
public/java/src/org/broadinstitute/sting/utils/pileup/PileupElementTracker.java
View File

@ -34,11 +34,20 @@ import java.util.*;
*/
abstract class PileupElementTracker<PE extends PileupElement> implements Iterable<PE> {
public abstract int size();
public abstract PileupElementTracker<PE> copy();
}
class UnifiedPileupElementTracker<PE extends PileupElement> extends PileupElementTracker<PE> {
private final List<PE> pileup;
@Override
public UnifiedPileupElementTracker<PE> copy() {
UnifiedPileupElementTracker<PE> result = new UnifiedPileupElementTracker<PE>();
for(PE element : pileup)
result.add(element);
return result;
}
public UnifiedPileupElementTracker() { pileup = new LinkedList<PE>(); }
public UnifiedPileupElementTracker(List<PE> pileup) { this.pileup = pileup; }
@ -65,6 +74,14 @@ class PerSamplePileupElementTracker<PE extends PileupElement> extends PileupElem
pileup = new HashMap<String,PileupElementTracker<PE>>();
}
public PerSamplePileupElementTracker<PE> copy() {
PerSamplePileupElementTracker<PE> result = new PerSamplePileupElementTracker<PE>();
for (Map.Entry<String, PileupElementTracker<PE>> entry : pileup.entrySet())
result.addElements(entry.getKey(), entry.getValue());
return result;
}
/**
* Gets a list of all the samples stored in this pileup.
* @return List of samples in this pileup.

8
public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedPileup.java
View File

@ -283,4 +283,12 @@ public interface ReadBackedPileup extends Iterable<PileupElement>, HasGenomeLoca
* @return
*/
public FragmentCollection<PileupElement> toFragments();
/**
* Creates a full copy (not shallow) of the ReadBacked Pileup
*
* @return
*/
public ReadBackedPileup copy();
}

3
public/java/src/org/broadinstitute/sting/utils/sam/GATKSAMRecord.java
View File

@ -397,6 +397,9 @@ public class GATKSAMRecord extends BAMRecord {
else if (op != CigarOperator.HARD_CLIP)
break;
}
if ( softStart < 1 )
softStart = 1;
}
return softStart;
}

67
public/java/src/org/broadinstitute/sting/utils/sam/MisencodedBaseQualityReadTransformer.java 100644
View File

@ -0,0 +1,67 @@
package org.broadinstitute.sting.utils.sam;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.iterators.ReadTransformer;
import org.broadinstitute.sting.gatk.walkers.Walker;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.exceptions.UserException;
/**
* Checks for and errors out (or fixes if requested) when it detects reads with base qualities that are not encoded with
* phred-scaled quality scores. Q0 == ASCII 33 according to the SAM specification, whereas Illumina encoding starts at
* Q64. The idea here is simple: if we are asked to fix the scores then we just subtract 31 from every quality score.
* Otherwise, we randomly sample reads (for efficiency) and error out if we encounter a qual that's too high.
*/
public class MisencodedBaseQualityReadTransformer extends ReadTransformer {
private static final int samplingFrequency = 1000; // sample 1 read for every 1000 encountered
private static final int encodingFixValue = 31; // Illumina_64 - PHRED_33
private boolean disabled;
private boolean fixQuals;
private static int currentReadCounter = 0;
@Override
public ApplicationTime initializeSub(final GenomeAnalysisEngine engine, final Walker walker) {
fixQuals = engine.getArguments().FIX_MISENCODED_QUALS;
disabled = !fixQuals && engine.getArguments().ALLOW_POTENTIALLY_MISENCODED_QUALS;
return ReadTransformer.ApplicationTime.ON_INPUT;
}
@Override
public boolean enabled() {
return !disabled;
}
@Override
public GATKSAMRecord apply(final GATKSAMRecord read) {
if ( fixQuals )
return fixMisencodedQuals(read);
checkForMisencodedQuals(read);
return read;
}
protected static GATKSAMRecord fixMisencodedQuals(final GATKSAMRecord read) {
final byte[] quals = read.getBaseQualities();
for ( int i = 0; i < quals.length; i++ ) {
quals[i] -= encodingFixValue;
}
read.setBaseQualities(quals);
return read;
}
protected static void checkForMisencodedQuals(final GATKSAMRecord read) {
// sample reads randomly for checking
if ( ++currentReadCounter >= samplingFrequency ) {
currentReadCounter = 0;
final byte[] quals = read.getBaseQualities();
for ( final byte qual : quals ) {
if ( qual > QualityUtils.MAX_REASONABLE_Q_SCORE )
throw new UserException.MisencodedBAM(read, "we encountered an extremely high quality score of " + (int)qual);
}
}
}
}

6
public/java/test/org/broadinstitute/sting/utils/baq/BAQUnitTest.java
View File

@ -1,10 +1,6 @@
// our package
package org.broadinstitute.sting.utils.baq;
// the imports for unit testing.
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.testng.Assert;
import org.testng.annotations.Test;
@ -24,7 +20,7 @@ import net.sf.picard.reference.IndexedFastaSequenceFile;
import net.sf.samtools.*;
/**
* Basic unit test for GenomeLoc
* Basic unit test for BAQ calculation
*/
public class BAQUnitTest extends BaseTest {
private SAMFileHeader header;

66
public/java/test/org/broadinstitute/sting/utils/sam/MisencodedBaseQualityUnitTest.java 100644
View File

@ -0,0 +1,66 @@
package org.broadinstitute.sting.utils.sam;
import net.sf.samtools.SAMFileHeader;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.testng.Assert;
import org.testng.annotations.BeforeMethod;
import org.testng.annotations.Test;
import java.util.ArrayList;
import java.util.List;
/**
* Basic unit test for misencoded quals
*/
public class MisencodedBaseQualityUnitTest extends BaseTest {
private static final String readBases = "AAAAAAAAAA";
private static final byte[] badQuals = { 59, 60, 62, 63, 64, 61, 62, 58, 57, 56 };
private static final byte[] goodQuals = { 60, 60, 60, 60, 60, 60, 60, 60, 60, 60 };
private static final byte[] fixedQuals = { 28, 29, 31, 32, 33, 30, 31, 27, 26, 25 };
private SAMFileHeader header;
@BeforeMethod
public void before() {
header = ArtificialSAMUtils.createArtificialSamHeader(1, 1, 1000);
}
private GATKSAMRecord createRead(final boolean useGoodBases) {
GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "foo", 0, 10, readBases.getBytes(), useGoodBases ? goodQuals : badQuals);
read.setCigarString("10M");
return read;
}
@Test(enabled = true)
public void testGoodQuals() {
final List<GATKSAMRecord> reads = new ArrayList<GATKSAMRecord>(10000);
for ( int i = 0; i < 10000; i++ )
reads.add(createRead(true));
testEncoding(reads);
}
@Test(enabled = true, expectedExceptions = {UserException.class})
public void testBadQualsThrowsError() {
final List<GATKSAMRecord> reads = new ArrayList<GATKSAMRecord>(10000);
for ( int i = 0; i < 10000; i++ )
reads.add(createRead(false));
testEncoding(reads);
}
@Test(enabled = true)
public void testFixBadQuals() {
final GATKSAMRecord read = createRead(false);
final GATKSAMRecord fixedRead = MisencodedBaseQualityReadTransformer.fixMisencodedQuals(read);
for ( int i = 0; i < fixedQuals.length; i++ )
Assert.assertEquals(fixedQuals[i], fixedRead.getBaseQualities()[i]);
}
private void testEncoding(final List<GATKSAMRecord> reads) {
for ( final GATKSAMRecord read : reads )
MisencodedBaseQualityReadTransformer.checkForMisencodedQuals(read);
}
}

3
public/packages/GATKEngine.xml
View File

@ -36,6 +36,9 @@
<dir name="org/broadinstitute/sting/utils/R" includes="**/*.tar.gz" />
<!-- All R scripts in org.broadinstitute.sting -->
<dir name="org/broadinstitute/sting" includes="**/*.R" />
<!-- Resources in org.broadinstitute.sting -->
<dir name="org/broadinstitute/sting" includes="**/resources/*" />
<dir name="org/broadinstitute/sting" includes="**/templates/*" />
<!-- The GATK public key -->
<file path="GATK_public.key" />
</dependencies>

4
public/scala/src/org/broadinstitute/sting/queue/QSettings.scala
View File

@ -72,6 +72,10 @@ class QSettings {
@Argument(fullName="resident_memory_request_parameter", shortName="resMemReqParam", doc="Parameter for resident memory requests. By default not requested.", required=false)
var residentRequestParameter: String = _
@Argument(fullName="job_walltime", shortName="wallTime", doc="Setting the required DRMAA walltime or LSF run limit.", required=false)
@ClassType(classOf[Long])
var jobWalltime: Option[Long] = None
/** The name of the parallel environment (required for SGE, for example) */
@Argument(fullName="job_parallel_env", shortName="jobParaEnv", doc="An SGE style parallel environment to use for jobs requesting more than 1 core. Equivalent to submitting jobs with -pe ARG nt for jobs with nt > 1", required=false)
var parallelEnvironmentName: String = "smp_pe" // Broad default

3
public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala
View File

@ -65,6 +65,9 @@ class DrmaaJobRunner(val session: Session, val function: CommandLineFunction) ex
drmaaJob.setJoinFiles(true)
}
if(!function.wallTime.isEmpty)
drmaaJob.setHardWallclockTimeLimit(function.wallTime.get)
drmaaJob.setNativeSpecification(functionNativeSpec)
// Instead of running the function.commandLine, run "sh <jobScript>"

7
public/scala/src/org/broadinstitute/sting/queue/engine/lsf/Lsf706JobRunner.scala
View File

@ -151,8 +151,11 @@ class Lsf706JobRunner(val function: CommandLineFunction) extends CommandLineJobR
throw new QException("setOption_() returned -1 while setting esub");
}
// LSF specific: get the max runtime for the jobQueue and pass it for this job
request.rLimits(LibLsf.LSF_RLIMIT_RUN) = Lsf706JobRunner.getRlimitRun(function.jobQueue)
if(!function.wallTime.isEmpty)
request.rLimits(LibLsf.LSF_RLIMIT_RUN) = function.wallTime.get.toInt
else
// LSF specific: get the max runtime for the jobQueue and pass it for this job
request.rLimits(LibLsf.LSF_RLIMIT_RUN) = Lsf706JobRunner.getRlimitRun(function.jobQueue)
// Run the command as sh <jobScript>
request.command = "sh " + jobScript

10
public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala
View File

@ -33,6 +33,9 @@ import org.broadinstitute.sting.commandline.Argument
trait CommandLineFunction extends QFunction with Logging {
def commandLine: String
/** Setting the wall time request for DRMAA / run limit for LSF */
var wallTime: Option[Long] = None
/** Upper memory limit */
@Argument(doc="Memory limit", required=false)
var memoryLimit: Option[Double] = None
@ -67,6 +70,9 @@ trait CommandLineFunction extends QFunction with Logging {
super.copySettingsTo(function)
function match {
case commandLineFunction: CommandLineFunction =>
if(commandLineFunction.wallTime.isEmpty)
commandLineFunction.wallTime = this.wallTime
if (commandLineFunction.memoryLimit.isEmpty)
commandLineFunction.memoryLimit = this.memoryLimit
@ -110,6 +116,10 @@ trait CommandLineFunction extends QFunction with Logging {
* Sets all field values.
*/
override def freezeFieldValues() {
if(wallTime.isEmpty)
wallTime = qSettings.jobWalltime
if (jobQueue == null)
jobQueue = qSettings.jobQueue

11
public/scala/test/org/broadinstitute/sting/queue/pipeline/examples/HelloWorldPipelineTest.scala
View File

@ -126,4 +126,15 @@ class HelloWorldPipelineTest {
spec.jobRunners = Seq("GridEngine")
PipelineTest.executeTest(spec)
}
// disabled because our DRMAA implementation doesn't support wallTime
@Test(enabled=false, timeOut=36000000)
def testHelloWorldWithWalltime() {
val spec = new PipelineTestSpec
spec.name = "HelloWorldWithWalltime"
spec.args = "-S public/scala/qscript/org/broadinstitute/sting/queue/qscripts/examples/HelloWorld.scala" +
" -wallTime 100"
spec.jobRunners = PipelineTest.allJobRunners
PipelineTest.executeTest(spec)
}
}