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

Merge branch 'master' of github.com:broadinstitute/gsa-unstable

This commit is contained in:
Ami Levy-Moonshine 2012-12-05 15:07:12 -05:00
parent d0b8cc7773 d0cab795b7
commit 5d78a61f7a
37 changed files with 569 additions and 430 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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 " + privateTestDir + "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

4
public/java/src/org/broadinstitute/sting/gatk/executive/MicroScheduler.java
View File

@ -43,7 +43,6 @@ import org.broadinstitute.sting.utils.AutoFormattingTime;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler;
import org.broadinstitute.sting.utils.progressmeter.ProgressMeter;
import org.broadinstitute.sting.utils.threading.ThreadEfficiencyMonitor;
@ -346,9 +345,6 @@ public abstract class MicroScheduler implements MicroSchedulerMBean {
for ( final TraversalEngine te : allCreatedTraversalEngines)
te.shutdown();
// horrible hack to print nano scheduling information across all nano schedulers, if any were used
NanoScheduler.printCombinedRuntimeProfile();
allCreatedTraversalEngines.clear();
availableTraversalEngines.clear();
}

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);

10
public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java
View File

@ -495,4 +495,14 @@ public class GenomeLoc implements Comparable<GenomeLoc>, Serializable, HasGenome
public long sizeOfOverlap( final GenomeLoc that ) {
return ( this.overlapsP(that) ? Math.min( getStop(), that.getStop() ) - Math.max( getStart(), that.getStart() ) + 1L : 0L );
}
/**
* Returns the maximum GenomeLoc of this and other
* @param other another non-null genome loc
* @return the max of this and other
*/
public GenomeLoc max(final GenomeLoc other) {
final int cmp = this.compareTo(other);
return cmp == -1 ? other : this;
}
}

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 {

21
public/java/src/org/broadinstitute/sting/utils/nanoScheduler/InputProducer.java
View File

@ -2,7 +2,6 @@ package org.broadinstitute.sting.utils.nanoScheduler;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.MultiThreadedErrorTracker;
import org.broadinstitute.sting.utils.SimpleTimer;
import java.util.Iterator;
import java.util.concurrent.BlockingQueue;
@ -19,11 +18,6 @@ class InputProducer<InputType> implements Runnable {
*/
final Iterator<InputType> inputReader;
/**
* Our timer (may be null) that we use to track our input costs
*/
final SimpleTimer inputTimer;
/**
* Where we put our input values for consumption
*/
@ -51,16 +45,13 @@ class InputProducer<InputType> implements Runnable {
public InputProducer(final Iterator<InputType> inputReader,
final MultiThreadedErrorTracker errorTracker,
final SimpleTimer inputTimer,
final BlockingQueue<InputValue> outputQueue) {
if ( inputReader == null ) throw new IllegalArgumentException("inputReader cannot be null");
if ( errorTracker == null ) throw new IllegalArgumentException("errorTracker cannot be null");
if ( inputTimer == null ) throw new IllegalArgumentException("inputTimer cannot be null");
if ( outputQueue == null ) throw new IllegalArgumentException("OutputQueue cannot be null");
this.inputReader = inputReader;
this.errorTracker = errorTracker;
this.inputTimer = inputTimer;
this.outputQueue = outputQueue;
}
@ -94,16 +85,15 @@ class InputProducer<InputType> implements Runnable {
* @throws InterruptedException
*/
private synchronized InputType readNextItem() throws InterruptedException {
inputTimer.restart();
if ( ! inputReader.hasNext() ) {
// we are done, mark ourselves as such and return null
readLastValue = true;
inputTimer.stop();
return null;
} else {
// get the next value, and return it
final InputType input = inputReader.next();
inputTimer.stop();
if ( input == null )
throw new IllegalStateException("inputReader.next() returned a null value, breaking our contract");
nRead++;
return input;
}
@ -121,6 +111,9 @@ class InputProducer<InputType> implements Runnable {
final InputType value = readNextItem();
if ( value == null ) {
if ( ! readLastValue )
throw new IllegalStateException("value == null but readLastValue is false!");
// add the EOF object so our consumer knows we are done in all inputs
// note that we do not increase inputID here, so that variable indicates the ID
// of the last real value read from the queue
@ -133,8 +126,10 @@ class InputProducer<InputType> implements Runnable {
}
latch.countDown();
} catch (Exception ex) {
} catch (Throwable ex) {
errorTracker.notifyOfError(ex);
} finally {
// logger.info("Exiting input thread readLastValue = " + readLastValue);
}
}

67
public/java/src/org/broadinstitute/sting/utils/nanoScheduler/NSRuntimeProfile.java
View File

@ -1,67 +0,0 @@
package org.broadinstitute.sting.utils.nanoScheduler;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.AutoFormattingTime;
import org.broadinstitute.sting.utils.SimpleTimer;
/**
* Holds runtime profile (input, read, map) times as tracked by NanoScheduler
*
* User: depristo
* Date: 9/10/12
* Time: 8:31 PM
*/
public class NSRuntimeProfile {
final SimpleTimer outsideSchedulerTimer = new SimpleTimer("outside");
final SimpleTimer inputTimer = new SimpleTimer("input");
final SimpleTimer mapTimer = new SimpleTimer("map");
final SimpleTimer reduceTimer = new SimpleTimer("reduce");
/**
* Combine the elapsed time information from other with this profile
*
* @param other a non-null profile
*/
public void combine(final NSRuntimeProfile other) {
outsideSchedulerTimer.addElapsed(other.outsideSchedulerTimer);
inputTimer.addElapsed(other.inputTimer);
mapTimer.addElapsed(other.mapTimer);
reduceTimer.addElapsed(other.reduceTimer);
}
/**
* Print the runtime profiling to logger
*
* @param logger
*/
public void log(final Logger logger) {
log1(logger, "Input time", inputTimer);
log1(logger, "Map time", mapTimer);
log1(logger, "Reduce time", reduceTimer);
log1(logger, "Outside time", outsideSchedulerTimer);
}
/**
* @return the total runtime for all functions of this nano scheduler
*/
//@Ensures("result >= 0.0")
public double totalRuntimeInSeconds() {
return inputTimer.getElapsedTime()
+ mapTimer.getElapsedTime()
+ reduceTimer.getElapsedTime()
+ outsideSchedulerTimer.getElapsedTime();
}
/**
* Print to logger.info timing information from timer, with name label
*
* @param label the name of the timer to display. Should be human readable
* @param timer the timer whose elapsed time we will display
*/
//@Requires({"label != null", "timer != null"})
private void log1(final Logger logger, final String label, final SimpleTimer timer) {
final double myTimeInSec = timer.getElapsedTime();
final double myTimePercent = myTimeInSec / totalRuntimeInSeconds() * 100;
logger.info(String.format("%s: %s (%5.2f%%)", label, new AutoFormattingTime(myTimeInSec), myTimePercent));
}
}

82
public/java/src/org/broadinstitute/sting/utils/nanoScheduler/NanoScheduler.java
View File

@ -57,16 +57,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
boolean debug = false;
private NSProgressFunction<InputType> progressFunction = null;
/**
* Tracks the combined runtime profiles across all created nano schedulers
*/
final static private NSRuntimeProfile combinedNSRuntimeProfiler = new NSRuntimeProfile();
/**
* The profile specific to this nano scheduler
*/
final private NSRuntimeProfile myNSRuntimeProfile = new NSRuntimeProfile();
/**
* Create a new nanoscheduler with the desire characteristics requested by the argument
*
@ -94,9 +84,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
this.inputExecutor = Executors.newSingleThreadExecutor(new NamedThreadFactory("NS-input-thread-%d"));
this.masterExecutor = Executors.newSingleThreadExecutor(new NamedThreadFactory("NS-master-thread-%d"));
}
// start timing the time spent outside of the nanoScheduler
myNSRuntimeProfile.outsideSchedulerTimer.start();
}
/**
@ -123,11 +110,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
* After this call, execute cannot be invoked without throwing an error
*/
public void shutdown() {
myNSRuntimeProfile.outsideSchedulerTimer.stop();
// add my timing information to the combined NS runtime profile
combinedNSRuntimeProfiler.combine(myNSRuntimeProfile);
if ( nThreads > 1 ) {
shutdownExecutor("inputExecutor", inputExecutor);
shutdownExecutor("mapExecutor", mapExecutor);
@ -137,19 +119,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
shutdown = true;
}
public void printRuntimeProfile() {
myNSRuntimeProfile.log(logger);
}
public static void printCombinedRuntimeProfile() {
if ( combinedNSRuntimeProfiler.totalRuntimeInSeconds() > 0.1 )
combinedNSRuntimeProfiler.log(logger);
}
protected double getTotalRuntime() {
return myNSRuntimeProfile.totalRuntimeInSeconds();
}
/**
* Helper function to cleanly shutdown an execution service, checking that the execution
* state is clean when it's done.
@ -245,8 +214,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
if ( map == null ) throw new IllegalArgumentException("map function cannot be null");
if ( reduce == null ) throw new IllegalArgumentException("reduce function cannot be null");
myNSRuntimeProfile.outsideSchedulerTimer.stop();
ReduceType result;
if ( ALLOW_SINGLE_THREAD_FASTPATH && getnThreads() == 1 ) {
result = executeSingleThreaded(inputReader, map, initialValue, reduce);
@ -254,7 +221,6 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
result = executeMultiThreaded(inputReader, map, initialValue, reduce);
}
myNSRuntimeProfile.outsideSchedulerTimer.restart();
return result;
}
@ -273,28 +239,19 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
while ( true ) {
// start timer to ensure that both hasNext and next are caught by the timer
myNSRuntimeProfile.inputTimer.restart();
if ( ! inputReader.hasNext() ) {
myNSRuntimeProfile.inputTimer.stop();
break;
} else {
final InputType input = inputReader.next();
myNSRuntimeProfile.inputTimer.stop();
// map
myNSRuntimeProfile.mapTimer.restart();
final long preMapTime = LOG_MAP_TIMES ? 0 : myNSRuntimeProfile.mapTimer.currentTimeNano();
final MapType mapValue = map.apply(input);
if ( LOG_MAP_TIMES ) logger.info("MAP TIME " + (myNSRuntimeProfile.mapTimer.currentTimeNano() - preMapTime));
myNSRuntimeProfile.mapTimer.stop();
if ( i++ % this.bufferSize == 0 && progressFunction != null )
if ( progressFunction != null )
progressFunction.progress(input);
// reduce
myNSRuntimeProfile.reduceTimer.restart();
sum = reduce.apply(mapValue, sum);
myNSRuntimeProfile.reduceTimer.stop();
}
}
@ -320,6 +277,7 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
while ( true ) {
// check that no errors occurred while we were waiting
handleErrors();
// checkForDeadlocks();
try {
final ReduceType result = reduceResult.get(100, TimeUnit.MILLISECONDS);
@ -341,6 +299,26 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
}
}
// private void checkForDeadlocks() {
// if ( deadLockCheckCounter++ % 100 == 0 ) {
// logger.info("Checking for deadlocks...");
// final ThreadMXBean bean = ManagementFactory.getThreadMXBean();
// final long[] threadIds = bean.findDeadlockedThreads(); // Returns null if no threads are deadlocked.
//
// if (threadIds != null) {
// final ThreadInfo[] infos = bean.getThreadInfo(threadIds);
//
// logger.error("!!! Deadlock detected !!!!");
// for (final ThreadInfo info : infos) {
// logger.error("Thread " + info);
// for ( final StackTraceElement elt : info.getStackTrace() ) {
// logger.error("\t" + elt.toString());
// }
// }
// }
// }
// }
private void handleErrors() {
if ( errorTracker.hasAnErrorOccurred() ) {
masterExecutor.shutdownNow();
@ -380,7 +358,7 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
// Create the input producer and start it running
final InputProducer<InputType> inputProducer =
new InputProducer<InputType>(inputReader, errorTracker, myNSRuntimeProfile.inputTimer, inputQueue);
new InputProducer<InputType>(inputReader, errorTracker, inputQueue);
inputExecutor.submit(inputProducer);
// a priority queue that stores up to bufferSize elements
@ -389,7 +367,7 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
new PriorityBlockingQueue<MapResult<MapType>>();
final Reducer<MapType, ReduceType> reducer
= new Reducer<MapType, ReduceType>(reduce, errorTracker, myNSRuntimeProfile.reduceTimer, initialValue);
= new Reducer<MapType, ReduceType>(reduce, errorTracker, initialValue);
try {
int nSubmittedJobs = 0;
@ -408,7 +386,8 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
// wait for all of the input and map threads to finish
return waitForCompletion(inputProducer, reducer);
} catch (Exception ex) {
} catch (Throwable ex) {
// logger.warn("Reduce job got exception " + ex);
errorTracker.notifyOfError(ex);
return initialValue;
}
@ -486,16 +465,12 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
final InputType input = inputWrapper.getValue();
// map
myNSRuntimeProfile.mapTimer.restart();
final long preMapTime = LOG_MAP_TIMES ? 0 : myNSRuntimeProfile.mapTimer.currentTimeNano();
final MapType mapValue = map.apply(input);
if ( LOG_MAP_TIMES ) logger.info("MAP TIME " + (myNSRuntimeProfile.mapTimer.currentTimeNano() - preMapTime));
myNSRuntimeProfile.mapTimer.stop();
// enqueue the result into the mapResultQueue
result = new MapResult<MapType>(mapValue, jobID);
if ( jobID % bufferSize == 0 && progressFunction != null )
if ( progressFunction != null )
progressFunction.progress(input);
} else {
// push back the EOF marker so other waiting threads can read it
@ -508,7 +483,8 @@ public class NanoScheduler<InputType, MapType, ReduceType> {
mapResultQueue.put(result);
final int nReduced = reducer.reduceAsMuchAsPossible(mapResultQueue);
} catch (Exception ex) {
} catch (Throwable ex) {
// logger.warn("Map job got exception " + ex);
errorTracker.notifyOfError(ex);
} finally {
// we finished a map job, release the job queue semaphore

9
public/java/src/org/broadinstitute/sting/utils/nanoScheduler/Reducer.java
View File

@ -4,7 +4,6 @@ import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.MultiThreadedErrorTracker;
import org.broadinstitute.sting.utils.SimpleTimer;
import java.util.concurrent.CountDownLatch;
import java.util.concurrent.PriorityBlockingQueue;
@ -34,7 +33,6 @@ class Reducer<MapType, ReduceType> {
final CountDownLatch countDownLatch = new CountDownLatch(1);
final NSReduceFunction<MapType, ReduceType> reduce;
final SimpleTimer reduceTimer;
final MultiThreadedErrorTracker errorTracker;
/**
@ -61,20 +59,16 @@ class Reducer<MapType, ReduceType> {
* reduceTimer
*
* @param reduce the reduce function to apply
* @param reduceTimer the timer to time the reduce function call
* @param initialSum the initial reduce sum
*/
public Reducer(final NSReduceFunction<MapType, ReduceType> reduce,
final MultiThreadedErrorTracker errorTracker,
final SimpleTimer reduceTimer,
final ReduceType initialSum) {
if ( errorTracker == null ) throw new IllegalArgumentException("Error tracker cannot be null");
if ( reduce == null ) throw new IllegalArgumentException("Reduce function cannot be null");
if ( reduceTimer == null ) throw new IllegalArgumentException("reduceTimer cannot be null");
this.errorTracker = errorTracker;
this.reduce = reduce;
this.reduceTimer = reduceTimer;
this.sum = initialSum;
}
@ -125,10 +119,7 @@ class Reducer<MapType, ReduceType> {
nReducesNow++;
// apply reduce, keeping track of sum
reduceTimer.restart();
sum = reduce.apply(result.getValue(), sum);
reduceTimer.stop();
}
numJobsReduced++;

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

@ -1054,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();
}

65
public/java/src/org/broadinstitute/sting/utils/progressmeter/ProgressMeter.java
View File

@ -26,6 +26,7 @@ package org.broadinstitute.sting.utils.progressmeter;
import com.google.java.contract.Ensures;
import com.google.java.contract.Invariant;
import com.google.java.contract.Requires;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.exceptions.UserException;
@ -143,6 +144,12 @@ public class ProgressMeter {
/** We use the SimpleTimer to time our run */
private final SimpleTimer timer = new SimpleTimer();
private GenomeLoc maxGenomeLoc = null;
private String positionMessage = "starting";
private long nTotalRecordsProcessed = 0;
final ProgressMeterDaemon progressMeterDaemon;
/**
* Create a new ProgressMeter
*
@ -177,21 +184,15 @@ public class ProgressMeter {
targetSizeInBP = processingIntervals.coveredSize();
// start up the timer
progressMeterDaemon = new ProgressMeterDaemon(this);
start();
}
/**
* Forward request to notifyOfProgress
*
* Assumes that one cycle has been completed
*
* @param loc our current location. Null means "in unmapped reads"
* @param nTotalRecordsProcessed the total number of records we've processed
* Start up the progress meter, printing initialization message and starting up the
* daemon thread for periodic printing.
*/
public void notifyOfProgress(final GenomeLoc loc, final long nTotalRecordsProcessed) {
notifyOfProgress(loc, false, nTotalRecordsProcessed);
}
@Requires("progressMeterDaemon != null")
private synchronized void start() {
timer.start();
lastProgressPrintTime = timer.currentTime();
@ -199,6 +200,8 @@ public class ProgressMeter {
logger.info("[INITIALIZATION COMPLETE; STARTING PROCESSING]");
logger.info(String.format("%15s processed.%s runtime per.1M.%s completed total.runtime remaining",
"Location", processingUnitName, processingUnitName));
progressMeterDaemon.start();
}
/**
@ -216,19 +219,41 @@ public class ProgressMeter {
* Synchronized to ensure that even with multiple threads calling notifyOfProgress we still
* get one clean stream of meter logs.
*
* Note this thread doesn't actually print progress, unless must print is true, but just registers
* the progress itself. A separate printing daemon periodically polls the meter to print out
* progress
*
* @param loc Current location, can be null if you are at the end of the processing unit
* @param mustPrint If true, will print out info, regardless of time interval
* @param nTotalRecordsProcessed the total number of records we've processed
*/
private synchronized void notifyOfProgress(final GenomeLoc loc, boolean mustPrint, final long nTotalRecordsProcessed) {
public synchronized void notifyOfProgress(final GenomeLoc loc, final long nTotalRecordsProcessed) {
if ( nTotalRecordsProcessed < 0 ) throw new IllegalArgumentException("nTotalRecordsProcessed must be >= 0");
// weird comparison to ensure that loc == null (in unmapped reads) is keep before maxGenomeLoc == null (on startup)
this.maxGenomeLoc = loc == null ? loc : (maxGenomeLoc == null ? loc : loc.max(maxGenomeLoc));
this.nTotalRecordsProcessed = Math.max(this.nTotalRecordsProcessed, nTotalRecordsProcessed);
// a pretty name for our position
this.positionMessage = maxGenomeLoc == null
? "unmapped reads"
: String.format("%s:%d", maxGenomeLoc.getContig(), maxGenomeLoc.getStart());
}
/**
* Actually try to print out progress
*
* This function may print out if the progress print is due, but if not enough time has elapsed
* since the last print we will not print out information.
*
* @param mustPrint if true, progress will be printed regardless of the last time we printed progress
*/
protected synchronized void printProgress(final boolean mustPrint) {
final long curTime = timer.currentTime();
final boolean printProgress = mustPrint || maxElapsedIntervalForPrinting(curTime, lastProgressPrintTime, progressPrintFrequency);
final boolean printLog = performanceLog != null && maxElapsedIntervalForPrinting(curTime, lastPerformanceLogPrintTime, PERFORMANCE_LOG_PRINT_FREQUENCY);
if ( printProgress || printLog ) {
final ProgressMeterData progressData = takeProgressSnapshot(loc, nTotalRecordsProcessed);
final ProgressMeterData progressData = takeProgressSnapshot(maxGenomeLoc, nTotalRecordsProcessed);
final AutoFormattingTime elapsed = new AutoFormattingTime(progressData.getElapsedSeconds());
final AutoFormattingTime bpRate = new AutoFormattingTime(progressData.secondsPerMillionBP());
@ -241,13 +266,8 @@ public class ProgressMeter {
lastProgressPrintTime = curTime;
updateLoggerPrintFrequency(estTotalRuntime.getTimeInSeconds());
// a pretty name for our position
final String posName = loc == null
? (mustPrint ? "done" : "unmapped reads")
: String.format("%s:%d", loc.getContig(), loc.getStart());
logger.info(String.format("%15s %5.2e %s %s %5.1f%% %s %s",
posName, progressData.getUnitsProcessed()*1.0, elapsed, unitRate,
positionMessage, progressData.getUnitsProcessed()*1.0, elapsed, unitRate,
100*fractionGenomeTargetCompleted, estTotalRuntime, timeToCompletion));
}
@ -296,13 +316,18 @@ public class ProgressMeter {
*/
public void notifyDone(final long nTotalRecordsProcessed) {
// print out the progress meter
notifyOfProgress(null, true, nTotalRecordsProcessed);
this.nTotalRecordsProcessed = nTotalRecordsProcessed;
this.positionMessage = "done";
printProgress(true);
logger.info(String.format("Total runtime %.2f secs, %.2f min, %.2f hours",
timer.getElapsedTime(), timer.getElapsedTime() / 60, timer.getElapsedTime() / 3600));
if ( performanceLog != null )
performanceLog.close();
// shutdown our daemon thread
progressMeterDaemon.done();
}
/**

60
public/java/src/org/broadinstitute/sting/utils/progressmeter/ProgressMeterDaemon.java 100644
View File

@ -0,0 +1,60 @@
package org.broadinstitute.sting.utils.progressmeter;
/**
* Daemon thread that periodically prints the progress of the progress meter
*
* User: depristo
* Date: 12/4/12
* Time: 9:16 PM
*/
public final class ProgressMeterDaemon extends Thread {
/**
* How frequently should we poll and print progress?
*/
private final static long POLL_FREQUENCY_MILLISECONDS = 10 * 1000;
/**
* Are we to continue periodically printing status, or should we shut down?
*/
boolean done = false;
/**
* The meter we will call print on
*/
final ProgressMeter meter;
/**
* Create a new ProgressMeterDaemon printing progress for meter
* @param meter the progress meter to print progress of
*/
public ProgressMeterDaemon(final ProgressMeter meter) {
if ( meter == null ) throw new IllegalArgumentException("meter cannot be null");
this.meter = meter;
setDaemon(true);
setName("ProgressMeterDaemon");
}
/**
* Tells this daemon thread to shutdown at the next opportunity, as the progress
* metering is complete.
*/
public final void done() {
this.done = true;
}
/**
* Start up the ProgressMeterDaemon, polling every tens of seconds to print, if
* necessary, the provided progress meter. Never exits until the JVM is complete,
* or done() is called, as the thread is a daemon thread
*/
public void run() {
while (! done) {
meter.printProgress(false);
try {
Thread.sleep(POLL_FREQUENCY_MILLISECONDS);
} catch (InterruptedException e) {
throw new RuntimeException(e);
}
}
}
}

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;

5
public/java/test/org/broadinstitute/sting/utils/nanoScheduler/InputProducerUnitTest.java
View File

@ -2,7 +2,6 @@ package org.broadinstitute.sting.utils.nanoScheduler;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.MultiThreadedErrorTracker;
import org.broadinstitute.sting.utils.SimpleTimer;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
@ -46,7 +45,7 @@ public class InputProducerUnitTest extends BaseTest {
final LinkedBlockingDeque<InputProducer<Integer>.InputValue> readQueue =
new LinkedBlockingDeque<InputProducer<Integer>.InputValue>(queueSize);
final InputProducer<Integer> ip = new InputProducer<Integer>(elements.iterator(), new MultiThreadedErrorTracker(), new SimpleTimer(), readQueue);
final InputProducer<Integer> ip = new InputProducer<Integer>(elements.iterator(), new MultiThreadedErrorTracker(), readQueue);
final ExecutorService es = Executors.newSingleThreadExecutor();
@ -94,7 +93,7 @@ public class InputProducerUnitTest extends BaseTest {
final LinkedBlockingDeque<InputProducer<Integer>.InputValue> readQueue =
new LinkedBlockingDeque<InputProducer<Integer>.InputValue>();
final InputProducer<Integer> ip = new InputProducer<Integer>(elements.iterator(), new MultiThreadedErrorTracker(), new SimpleTimer(), readQueue);
final InputProducer<Integer> ip = new InputProducer<Integer>(elements.iterator(), new MultiThreadedErrorTracker(), readQueue);
final ExecutorService es = Executors.newSingleThreadExecutor();
es.submit(ip);

37
public/java/test/org/broadinstitute/sting/utils/nanoScheduler/NanoSchedulerUnitTest.java
View File

@ -188,17 +188,6 @@ public class NanoSchedulerUnitTest extends BaseTest {
Assert.assertTrue(callback.callBacks >= test.nExpectedCallbacks(), "Not enough callbacks detected. Expected at least " + test.nExpectedCallbacks() + " but saw only " + callback.callBacks);
nanoScheduler.shutdown();
// TODO -- need to enable only in the case where there's serious time spend in
// TODO -- read /map / reduce, otherwise the "outside" timer doesn't add up
final double myTimeEstimate = timer.getElapsedTime();
final double tolerance = 0.1;
if ( false && myTimeEstimate > 0.1 ) {
Assert.assertTrue(nanoScheduler.getTotalRuntime() > myTimeEstimate * tolerance,
"NanoScheduler said that the total runtime was " + nanoScheduler.getTotalRuntime()
+ " but the overall test time was " + myTimeEstimate + ", beyond our tolerance factor of "
+ tolerance);
}
}
@Test(enabled = true && ! DEBUG, dataProvider = "NanoSchedulerBasicTest", dependsOnMethods = "testMultiThreadedNanoScheduler", timeOut = NANO_SCHEDULE_MAX_RUNTIME)
@ -243,7 +232,7 @@ public class NanoSchedulerUnitTest extends BaseTest {
for ( final int nThreads : Arrays.asList(8) ) {
for ( final boolean addDelays : Arrays.asList(true, false) ) {
final NanoSchedulerBasicTest test = new NanoSchedulerBasicTest(bufSize, nThreads, 1, 1000000, false);
final int maxN = addDelays ? 10000 : 100000;
final int maxN = addDelays ? 1000 : 10000;
for ( int nElementsBeforeError = 0; nElementsBeforeError < maxN; nElementsBeforeError += Math.max(nElementsBeforeError / 10, 1) ) {
tests.add(new Object[]{nElementsBeforeError, test, addDelays});
}
@ -259,17 +248,22 @@ public class NanoSchedulerUnitTest extends BaseTest {
executeTestErrorThrowingInput(10, new NullPointerException(), exampleTest, false);
}
@Test(enabled = true, expectedExceptions = ReviewedStingException.class, timeOut = 10000)
@Test(enabled = true, expectedExceptions = ReviewedStingException.class, timeOut = 1000)
public void testInputErrorIsThrown_RSE() throws InterruptedException {
executeTestErrorThrowingInput(10, new ReviewedStingException("test"), exampleTest, false);
}
@Test(enabled = true, expectedExceptions = NullPointerException.class, dataProvider = "NanoSchedulerInputExceptionTest", timeOut = 10000, invocationCount = 1)
public void testInputErrorDoesntDeadlock(final int nElementsBeforeError, final NanoSchedulerBasicTest test, final boolean addDelays ) throws InterruptedException {
@Test(enabled = true, expectedExceptions = NullPointerException.class, dataProvider = "NanoSchedulerInputExceptionTest", timeOut = 1000, invocationCount = 1)
public void testInputRuntimeExceptionDoesntDeadlock(final int nElementsBeforeError, final NanoSchedulerBasicTest test, final boolean addDelays ) throws InterruptedException {
executeTestErrorThrowingInput(nElementsBeforeError, new NullPointerException(), test, addDelays);
}
private void executeTestErrorThrowingInput(final int nElementsBeforeError, final RuntimeException ex, final NanoSchedulerBasicTest test, final boolean addDelays) {
@Test(enabled = true, expectedExceptions = ReviewedStingException.class, dataProvider = "NanoSchedulerInputExceptionTest", timeOut = 1000, invocationCount = 1)
public void testInputErrorDoesntDeadlock(final int nElementsBeforeError, final NanoSchedulerBasicTest test, final boolean addDelays ) throws InterruptedException {
executeTestErrorThrowingInput(nElementsBeforeError, new Error(), test, addDelays);
}
private void executeTestErrorThrowingInput(final int nElementsBeforeError, final Throwable ex, final NanoSchedulerBasicTest test, final boolean addDelays) {
logger.warn("executeTestErrorThrowingInput " + nElementsBeforeError + " ex=" + ex + " test=" + test + " addInputDelays=" + addDelays);
final NanoScheduler<Integer, Integer, Integer> nanoScheduler = test.makeScheduler();
nanoScheduler.execute(new ErrorThrowingIterator(nElementsBeforeError, ex, addDelays), test.makeMap(), test.initReduce(), test.makeReduce());
@ -279,9 +273,9 @@ public class NanoSchedulerUnitTest extends BaseTest {
final int nElementsBeforeError;
final boolean addDelays;
int i = 0;
final RuntimeException ex;
final Throwable ex;
private ErrorThrowingIterator(final int nElementsBeforeError, RuntimeException ex, boolean addDelays) {
private ErrorThrowingIterator(final int nElementsBeforeError, Throwable ex, boolean addDelays) {
this.nElementsBeforeError = nElementsBeforeError;
this.ex = ex;
this.addDelays = addDelays;
@ -290,7 +284,12 @@ public class NanoSchedulerUnitTest extends BaseTest {
@Override public boolean hasNext() { return true; }
@Override public Integer next() {
if ( i++ > nElementsBeforeError ) {
throw ex;
if ( ex instanceof Error )
throw (Error)ex;
else if ( ex instanceof RuntimeException )
throw (RuntimeException)ex;
else
throw new RuntimeException("Bad exception " + ex);
} else if ( addDelays ) {
maybeDelayMe(i);
return i;

5
public/java/test/org/broadinstitute/sting/utils/nanoScheduler/ReducerUnitTest.java
View File

@ -2,7 +2,6 @@ package org.broadinstitute.sting.utils.nanoScheduler;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.MultiThreadedErrorTracker;
import org.broadinstitute.sting.utils.SimpleTimer;
import org.broadinstitute.sting.utils.Utils;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
@ -93,7 +92,7 @@ public class ReducerUnitTest extends BaseTest {
final List<List<MapResult<Integer>>> jobGroups = Utils.groupList(allJobs, groupSize);
final ReduceSumTest reduce = new ReduceSumTest();
final Reducer<Integer, Integer> reducer = new Reducer<Integer, Integer>(reduce, new MultiThreadedErrorTracker(), new SimpleTimer(), 0);
final Reducer<Integer, Integer> reducer = new Reducer<Integer, Integer>(reduce, new MultiThreadedErrorTracker(), 0);
final TestWaitingForFinalReduce waitingThread = new TestWaitingForFinalReduce(reducer, expectedSum(allJobs));
final ExecutorService es = Executors.newSingleThreadExecutor();
@ -155,7 +154,7 @@ public class ReducerUnitTest extends BaseTest {
private void runSettingJobIDTwice() throws Exception {
final PriorityBlockingQueue<MapResult<Integer>> mapResultsQueue = new PriorityBlockingQueue<MapResult<Integer>>();
final Reducer<Integer, Integer> reducer = new Reducer<Integer, Integer>(new ReduceSumTest(), new MultiThreadedErrorTracker(), new SimpleTimer(), 0);
final Reducer<Integer, Integer> reducer = new Reducer<Integer, Integer>(new ReduceSumTest(), new MultiThreadedErrorTracker(), 0);
reducer.setTotalJobCount(10);
reducer.setTotalJobCount(15);

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);
}
}