From 78ce822b6f56fc5c6cc43be77f0faa47fbeabba6 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Mon, 19 Nov 2012 09:07:04 -0500 Subject: [PATCH 1/9] Protect against NPE when using non-GATK reports for inputs expecting valid GATK reports --- .../broadinstitute/sting/gatk/report/GATKReportVersion.java | 3 +++ 1 file changed, 3 insertions(+) diff --git a/public/java/src/org/broadinstitute/sting/gatk/report/GATKReportVersion.java b/public/java/src/org/broadinstitute/sting/gatk/report/GATKReportVersion.java index b51fb17f0..1079d9b91 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/report/GATKReportVersion.java +++ b/public/java/src/org/broadinstitute/sting/gatk/report/GATKReportVersion.java @@ -80,6 +80,9 @@ public enum GATKReportVersion { * @return The version as an enum. */ public static GATKReportVersion fromHeader(String header) { + if ( header == null ) + throw new UserException.BadInput("The GATK report has no version specified in the header"); + if (header.startsWith("##:GATKReport.v0.1 ")) return GATKReportVersion.V0_1; From ff180a8e02eaeffbc42226c789c6c6946affae68 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Mon, 19 Nov 2012 09:09:57 -0500 Subject: [PATCH 2/9] Significant refactoring of the Haplotype Caller to handle problems with GGA. The main fix is that we now maintain a mapping from 'original' allele to 'Smith-Waterman-based' allele so that we no longer need to do a (buggy) matching throughout the calling process. --- .../haplotypecaller/GenotypingEngine.java | 253 ++++++------------ .../haplotypecaller/HaplotypeCaller.java | 6 +- .../LikelihoodCalculationEngine.java | 78 +++--- .../LikelihoodCalculationEngineUnitTest.java | 5 +- .../broadinstitute/sting/utils/Haplotype.java | 16 +- 5 files changed, 144 insertions(+), 214 deletions(-) diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java index d91df82e2..9fc636efe 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java @@ -31,7 +31,6 @@ import net.sf.samtools.Cigar; import net.sf.samtools.CigarElement; import org.apache.commons.lang.ArrayUtils; import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine; -import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext; import org.broadinstitute.sting.utils.*; import org.broadinstitute.sting.utils.collections.Pair; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; @@ -52,153 +51,17 @@ public class GenotypingEngine { noCall.add(Allele.NO_CALL); } - // WARN - // This function is the streamlined approach, currently not being used by default - // WARN - // WARN: This function is currently only being used by Menachem. Slated for removal/merging with the rest of the code. - // WARN - @Requires({"refLoc.containsP(activeRegionWindow)", "haplotypes.size() > 0"}) - public List>>> assignGenotypeLikelihoodsAndCallHaplotypeEvents( final UnifiedGenotyperEngine UG_engine, - final ArrayList haplotypes, - final byte[] ref, - final GenomeLoc refLoc, - final GenomeLoc activeRegionWindow, - final GenomeLocParser genomeLocParser ) { - // Prepare the list of haplotype indices to genotype - final ArrayList allelesToGenotype = new ArrayList(); - - for( final Haplotype h : haplotypes ) { - allelesToGenotype.add( Allele.create(h.getBases(), h.isReference()) ); - } - final int numHaplotypes = haplotypes.size(); - - // 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 double[] genotypeLikelihoods = new double[numHaplotypes * (numHaplotypes+1) / 2]; - final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(haplotypes, sample); - int glIndex = 0; - for( int iii = 0; iii < numHaplotypes; iii++ ) { - for( int jjj = 0; jjj <= iii; jjj++ ) { - genotypeLikelihoods[glIndex++] = haplotypeLikelihoodMatrix[iii][jjj]; // for example: AA,AB,BB,AC,BC,CC - } - } - genotypes.add(new GenotypeBuilder(sample, noCall).PL(genotypeLikelihoods).make()); - } - final VariantCallContext call = UG_engine.calculateGenotypes(new VariantContextBuilder().loc(activeRegionWindow).alleles(allelesToGenotype).genotypes(genotypes).make(), UG_engine.getUAC().GLmodel); - if( call == null ) { return Collections.emptyList(); } // exact model says that the call confidence is below the specified confidence threshold so nothing to do here - - // Prepare the list of haplotypes that need to be run through Smith-Waterman for output to VCF - final ArrayList haplotypesToRemove = new ArrayList(); - for( final Haplotype h : haplotypes ) { - if( call.getAllele(h.getBases()) == null ) { // exact model removed this allele from the list so no need to run SW and output to VCF - haplotypesToRemove.add(h); - } - } - haplotypes.removeAll(haplotypesToRemove); - - if( OUTPUT_FULL_HAPLOTYPE_SEQUENCE ) { - final List>>> returnVCs = new ArrayList>>>(); - // set up the default 1-to-1 haplotype mapping object - final HashMap> haplotypeMapping = new HashMap>(); - for( final Haplotype h : haplotypes ) { - final ArrayList list = new ArrayList(); - list.add(h); - haplotypeMapping.put(call.getAllele(h.getBases()), list); - } - returnVCs.add( new Pair>>(call,haplotypeMapping) ); - return returnVCs; - } - - final ArrayList>>> returnCalls = new ArrayList>>>(); - - // Using the cigar from each called haplotype figure out what events need to be written out in a VCF file - final TreeSet startPosKeySet = new TreeSet(); - int count = 0; - if( DEBUG ) { System.out.println("=== Best Haplotypes ==="); } - for( final Haplotype h : haplotypes ) { - if( DEBUG ) { - System.out.println( h.toString() ); - System.out.println( "> Cigar = " + h.getCigar() ); - } - // Walk along the alignment and turn any difference from the reference into an event - h.setEventMap( generateVCsFromAlignment( h, h.getAlignmentStartHapwrtRef(), h.getCigar(), ref, h.getBases(), refLoc, "HC" + count++ ) ); - startPosKeySet.addAll(h.getEventMap().keySet()); - } - - // Create the VC merge priority list - final ArrayList priorityList = new ArrayList(); - for( int iii = 0; iii < haplotypes.size(); iii++ ) { - priorityList.add("HC" + iii); - } - - // 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() ) { - final ArrayList eventsAtThisLoc = new ArrayList(); - for( final Haplotype h : haplotypes ) { - final HashMap eventMap = h.getEventMap(); - final VariantContext vc = eventMap.get(loc); - if( vc != null && !containsVCWithMatchingAlleles(eventsAtThisLoc, vc) ) { - eventsAtThisLoc.add(vc); - } - } - - // Create the allele mapping object which maps the original haplotype alleles to the alleles present in just this event - final ArrayList> alleleMapper = createAlleleMapper( loc, eventsAtThisLoc, haplotypes ); - - // Merge the event to find a common reference representation - final VariantContext mergedVC = VariantContextUtils.simpleMerge(genomeLocParser, eventsAtThisLoc, priorityList, VariantContextUtils.FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED, VariantContextUtils.GenotypeMergeType.PRIORITIZE, false, false, null, false, false); - - final HashMap> alleleHashMap = new HashMap>(); - int aCount = 0; - for( final Allele a : mergedVC.getAlleles() ) { - alleleHashMap.put(a, alleleMapper.get(aCount++)); // BUGBUG: needs to be cleaned up and merged with alleleMapper - } - - if( DEBUG ) { - System.out.println("Genotyping event at " + loc + " with alleles = " + mergedVC.getAlleles()); - //System.out.println("Event/haplotype allele mapping = " + alleleMapper); - } - - // Grab the genotype likelihoods from the appropriate places in the haplotype likelihood matrix -- calculation performed independently per sample - final GenotypesContext myGenotypes = GenotypesContext.create(haplotypes.get(0).getSampleKeySet().size()); - for( final String sample : haplotypes.get(0).getSampleKeySet() ) { // BUGBUG: assume all haplotypes saw the same samples - final int myNumHaplotypes = alleleMapper.size(); - final double[] genotypeLikelihoods = new double[myNumHaplotypes * (myNumHaplotypes+1) / 2]; - final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sample, alleleMapper); - int glIndex = 0; - for( int iii = 0; iii < myNumHaplotypes; iii++ ) { - for( int jjj = 0; jjj <= iii; jjj++ ) { - genotypeLikelihoods[glIndex++] = haplotypeLikelihoodMatrix[iii][jjj]; // for example: AA,AB,BB,AC,BC,CC - } - } - - // using the allele mapping object translate the haplotype allele into the event allele - final Genotype g = new GenotypeBuilder(sample) - .alleles(findEventAllelesInSample(mergedVC.getAlleles(), call.getAlleles(), call.getGenotype(sample).getAlleles(), alleleMapper, haplotypes)) - .phased(loc != startPosKeySet.first()) - .PL(genotypeLikelihoods).make(); - myGenotypes.add(g); - } - returnCalls.add( new Pair>>( - new VariantContextBuilder(mergedVC).log10PError(call.getLog10PError()).genotypes(myGenotypes).make(), alleleHashMap) ); - } - } - return returnCalls; - } - // BUGBUG: Create a class to hold this complicated return type @Requires({"refLoc.containsP(activeRegionWindow)", "haplotypes.size() > 0"}) - public List>>> assignGenotypeLikelihoodsAndCallIndependentEvents( final UnifiedGenotyperEngine UG_engine, - final ArrayList haplotypes, - final byte[] ref, - final GenomeLoc refLoc, - final GenomeLoc activeRegionWindow, - final GenomeLocParser genomeLocParser, - final ArrayList activeAllelesToGenotype ) { + public List>>> assignGenotypeLikelihoodsAndCallIndependentEvents( final UnifiedGenotyperEngine UG_engine, + final List haplotypes, + final byte[] ref, + final GenomeLoc refLoc, + final GenomeLoc activeRegionWindow, + final GenomeLocParser genomeLocParser, + final List activeAllelesToGenotype ) { - final ArrayList>>> returnCalls = new ArrayList>>>(); + final ArrayList>>> returnCalls = new ArrayList>>>(); // Using the cigar from each called haplotype figure out what events need to be written out in a VCF file final TreeSet startPosKeySet = new TreeSet(); @@ -261,7 +124,15 @@ public class GenotypingEngine { if( eventsAtThisLoc.isEmpty() ) { continue; } // Create the allele mapping object which maps the original haplotype alleles to the alleles present in just this event - final ArrayList> alleleMapper = createAlleleMapper( loc, eventsAtThisLoc, haplotypes ); + Map> alleleMapper = createAlleleMapper( loc, eventsAtThisLoc, haplotypes ); + + final Allele refAllele = eventsAtThisLoc.get(0).getReference(); + final ArrayList alleleOrdering = new ArrayList(alleleMapper.size()); + alleleOrdering.add(refAllele); + for ( final Allele allele : alleleMapper.keySet() ) { + if ( !refAllele.equals(allele) ) + alleleOrdering.add(allele); + } // Sanity check the priority list for( final VariantContext vc : eventsAtThisLoc ) { @@ -283,12 +154,6 @@ public class GenotypingEngine { final VariantContext mergedVC = VariantContextUtils.simpleMerge(genomeLocParser, eventsAtThisLoc, priorityList, VariantContextUtils.FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED, VariantContextUtils.GenotypeMergeType.PRIORITIZE, false, false, null, false, false); if( mergedVC == null ) { continue; } - HashMap> alleleHashMap = new HashMap>(); - int aCount = 0; - for( final Allele a : mergedVC.getAlleles() ) { - alleleHashMap.put(a, alleleMapper.get(aCount++)); // BUGBUG: needs to be cleaned up and merged with alleleMapper - } - if( DEBUG ) { System.out.println("Genotyping event at " + loc + " with alleles = " + mergedVC.getAlleles()); //System.out.println("Event/haplotype allele mapping = " + alleleMapper); @@ -299,7 +164,7 @@ public class GenotypingEngine { for( final String sample : haplotypes.get(0).getSampleKeySet() ) { // BUGBUG: assume all haplotypes saw the same samples final int numHaplotypes = alleleMapper.size(); final double[] genotypeLikelihoods = new double[numHaplotypes * (numHaplotypes+1) / 2]; - final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sample, alleleMapper); + final double[][] haplotypeLikelihoodMatrix = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sample, alleleMapper, alleleOrdering); int glIndex = 0; for( int iii = 0; iii < numHaplotypes; iii++ ) { for( int jjj = 0; jjj <= iii; jjj++ ) { @@ -313,23 +178,23 @@ public class GenotypingEngine { 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> alleleHashMapTrim = new HashMap>(); + final HashMap> alleleHashMapTrim = new HashMap>(); 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), alleleHashMap.get(call.getAlleles().get(iii))); + alleleHashMapTrim.put(vcCallTrim.getAlleles().get(iii), alleleMapper.get(call.getAlleles().get(iii))); } call = vcCallTrim; - alleleHashMap = alleleHashMapTrim; + alleleMapper = alleleHashMapTrim; } - returnCalls.add( new Pair>>(call, alleleHashMap) ); + returnCalls.add( new Pair>>(call, alleleMapper) ); } } } return returnCalls; } - protected static void cleanUpSymbolicUnassembledEvents( final ArrayList haplotypes ) { + protected static void cleanUpSymbolicUnassembledEvents( final List haplotypes ) { final ArrayList haplotypesToRemove = new ArrayList(); for( final Haplotype h : haplotypes ) { for( final VariantContext vc : h.getEventMap().values() ) { @@ -348,7 +213,7 @@ public class GenotypingEngine { haplotypes.removeAll(haplotypesToRemove); } - protected void mergeConsecutiveEventsBasedOnLD( final ArrayList haplotypes, final TreeSet startPosKeySet, final byte[] ref, final GenomeLoc refLoc ) { + protected void mergeConsecutiveEventsBasedOnLD( final List haplotypes, final TreeSet 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; } @@ -395,7 +260,9 @@ public class GenotypingEngine { final ArrayList haplotypeList = new ArrayList(); haplotypeList.add(h); for( final String sample : haplotypes.get(0).getSampleKeySet() ) { - final double haplotypeLikelihood = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods( haplotypeList, sample )[0][0]; + final HashSet sampleSet = new HashSet(1); + sampleSet.add(sample); + final double haplotypeLikelihood = LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods( sampleSet, haplotypeList )[0][0]; if( thisHapVC == null ) { if( nextHapVC == null ) { x11 = MathUtils.approximateLog10SumLog10(x11, haplotypeLikelihood); } else { x12 = MathUtils.approximateLog10SumLog10(x12, haplotypeLikelihood); } @@ -489,37 +356,71 @@ public class GenotypingEngine { @Requires({"haplotypes.size() >= eventsAtThisLoc.size() + 1"}) @Ensures({"result.size() == eventsAtThisLoc.size() + 1"}) - protected static ArrayList> createAlleleMapper( final int loc, final ArrayList eventsAtThisLoc, final ArrayList haplotypes ) { - final ArrayList> alleleMapper = new ArrayList>(); - final ArrayList refList = new ArrayList(); + protected static Map> createAlleleMapper( final int loc, final List eventsAtThisLoc, final List haplotypes ) { + + final Allele refAllele = eventsAtThisLoc.get(0).getReference(); + + final Map> alleleMapper = new HashMap>(eventsAtThisLoc.size()+1); for( final Haplotype h : haplotypes ) { if( h.getEventMap().get(loc) == null ) { // no event at this location so this is a reference-supporting haplotype - refList.add(h); + if ( !alleleMapper.containsKey(refAllele) ) + alleleMapper.put(refAllele, new ArrayList()); + alleleMapper.get(refAllele).add(h); + } else if ( h.isArtificialHaplotype() ) { + if ( !alleleMapper.containsKey(h.getArtificialAllele()) ) + alleleMapper.put(h.getArtificialAllele(), new ArrayList()); + alleleMapper.get(h.getArtificialAllele()).add(h); } else { - boolean foundInEventList = false; for( final VariantContext vcAtThisLoc : eventsAtThisLoc ) { if( h.getEventMap().get(loc).hasSameAllelesAs(vcAtThisLoc) ) { - foundInEventList = true; + final Allele altAllele = vcAtThisLoc.getAlternateAllele(0); + if ( !alleleMapper.containsKey(altAllele) ) + alleleMapper.put(altAllele, new ArrayList()); + alleleMapper.get(altAllele).add(h); + break; } } - if( !foundInEventList ) { // event at this location isn't one of the genotype-able options (during GGA) so this is a reference-supporting haplotype - refList.add(h); - } } } - alleleMapper.add(refList); - for( final VariantContext vcAtThisLoc : eventsAtThisLoc ) { - final ArrayList list = new ArrayList(); - for( final Haplotype h : haplotypes ) { - if( h.getEventMap().get(loc) != null && h.getEventMap().get(loc).hasSameAllelesAs(vcAtThisLoc) ) { - list.add(h); + + for( final Haplotype h : haplotypes ) { + if ( h.getEventMap().get(loc) == null || h.isArtificialHaplotype() ) + continue; + + Allele matchingAllele = null; + for ( final Map.Entry> alleleToTest : alleleMapper.entrySet() ) { + if ( alleleToTest.getKey().equals(refAllele) ) + continue; + + final Haplotype artificialHaplotype = alleleToTest.getValue().get(0); + if ( isSubSetOf(artificialHaplotype.getEventMap(), h.getEventMap()) ) { + matchingAllele = alleleToTest.getKey(); + break; } } - alleleMapper.add(list); + + if ( matchingAllele == null ) + matchingAllele = refAllele; + alleleMapper.get(matchingAllele).add(h); } + return alleleMapper; } + protected static boolean isSubSetOf(final Map subset, final Map superset) { + + for ( final Map.Entry fromSubset : subset.entrySet() ) { + final VariantContext fromSuperset = superset.get(fromSubset.getKey()); + if ( fromSuperset == null ) + return false; + + if ( !fromSuperset.hasAlternateAllele(fromSubset.getValue().getAlternateAllele(0)) ) + return false; + } + + return true; + } + @Ensures({"result.size() == haplotypeAllelesForSample.size()"}) protected static List findEventAllelesInSample( final List eventAlleles, final List haplotypeAlleles, final List haplotypeAllelesForSample, final ArrayList> alleleMapper, final ArrayList haplotypes ) { if( haplotypeAllelesForSample.contains(Allele.NO_CALL) ) { return noCall; } diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java index a185ba6af..2b739a321 100755 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java @@ -421,10 +421,8 @@ public class HaplotypeCaller extends ActiveRegionWalker implem // subset down to only the best haplotypes to be genotyped in all samples ( in GGA mode use all discovered haplotypes ) final ArrayList bestHaplotypes = ( UG_engine.getUAC().GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ? likelihoodCalculationEngine.selectBestHaplotypes( haplotypes ) : haplotypes ); - for( final Pair>> callResult : - ( GENOTYPE_FULL_ACTIVE_REGION && UG_engine.getUAC().GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES - ? genotypingEngine.assignGenotypeLikelihoodsAndCallHaplotypeEvents( UG_engine, bestHaplotypes, fullReferenceWithPadding, getPaddedLoc(activeRegion), activeRegion.getExtendedLoc(), getToolkit().getGenomeLocParser() ) - : genotypingEngine.assignGenotypeLikelihoodsAndCallIndependentEvents( UG_engine, bestHaplotypes, fullReferenceWithPadding, getPaddedLoc(activeRegion), activeRegion.getLocation(), getToolkit().getGenomeLocParser(), activeAllelesToGenotype ) ) ) { + for( final Pair>> callResult : + genotypingEngine.assignGenotypeLikelihoodsAndCallIndependentEvents( UG_engine, bestHaplotypes, fullReferenceWithPadding, getPaddedLoc(activeRegion), activeRegion.getLocation(), getToolkit().getGenomeLocParser(), activeAllelesToGenotype ) ) { if( DEBUG ) { System.out.println(callResult.getFirst().toStringWithoutGenotypes()); } final Map stratifiedReadMap = LikelihoodCalculationEngine.partitionReadsBasedOnLikelihoods( getToolkit().getGenomeLocParser(), perSampleReadList, perSampleFilteredReadList, callResult, UG_engine.getUAC().CONTAMINATION_FRACTION, UG_engine.getUAC().contaminationLog ); diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java index a0924623b..543987e74 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java @@ -148,34 +148,21 @@ public class LikelihoodCalculationEngine { return Math.min(b1.length, b2.length); } - @Requires({"haplotypes.size() > 0"}) - @Ensures({"result.length == result[0].length", "result.length == haplotypes.size()"}) - public static double[][] computeDiploidHaplotypeLikelihoods( final ArrayList haplotypes, final String sample ) { - // set up the default 1-to-1 haplotype mapping object, BUGBUG: target for future optimization? - final ArrayList> haplotypeMapping = new ArrayList>(); - for( final Haplotype h : haplotypes ) { - final ArrayList list = new ArrayList(); - list.add(h); - haplotypeMapping.add(list); - } - return computeDiploidHaplotypeLikelihoods( sample, haplotypeMapping ); - } - // 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 ArrayList> haplotypeMapping ) { + public static double[][] computeDiploidHaplotypeLikelihoods( final String sample, final Map> haplotypeMapping, final List alleleOrdering ) { final TreeSet sampleSet = new TreeSet(); sampleSet.add(sample); - return computeDiploidHaplotypeLikelihoods(sampleSet, haplotypeMapping); + return computeDiploidHaplotypeLikelihoods(sampleSet, haplotypeMapping, 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 samples, final ArrayList> haplotypeMapping ) { + public static double[][] computeDiploidHaplotypeLikelihoods( final Set samples, final Map> haplotypeMapping, final List alleleOrdering ) { - final int numHaplotypes = haplotypeMapping.size(); + final int numHaplotypes = alleleOrdering.size(); final double[][] haplotypeLikelihoodMatrix = new double[numHaplotypes][numHaplotypes]; for( int iii = 0; iii < numHaplotypes; iii++ ) { Arrays.fill(haplotypeLikelihoodMatrix[iii], Double.NEGATIVE_INFINITY); @@ -184,9 +171,9 @@ public class LikelihoodCalculationEngine { // 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++ ) { - for( int jjj = 0; jjj <= iii; jjj++ ) { - for( final Haplotype iii_mapped : haplotypeMapping.get(iii) ) { - for( final Haplotype jjj_mapped : haplotypeMapping.get(jjj) ) { + 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); @@ -200,12 +187,48 @@ public class LikelihoodCalculationEngine { } haplotypeLikelihoodMatrix[iii][jjj] = Math.max(haplotypeLikelihoodMatrix[iii][jjj], haplotypeLikelihood); } - } + } } } // normalize the diploid likelihoods matrix - return normalizeDiploidLikelihoodMatrixFromLog10( haplotypeLikelihoodMatrix ); + return normalizeDiploidLikelihoodMatrixFromLog10( haplotypeLikelihoodMatrix ); + } + + // 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 samples, final List 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); + 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++ ) { + // 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 ); } @Requires({"likelihoodMatrix.length == likelihoodMatrix[0].length"}) @@ -296,14 +319,7 @@ public class LikelihoodCalculationEngine { final Set sampleKeySet = haplotypes.get(0).getSampleKeySet(); // BUGBUG: assume all haplotypes saw the same samples final ArrayList bestHaplotypesIndexList = new ArrayList(); bestHaplotypesIndexList.add( findReferenceIndex(haplotypes) ); // always start with the reference haplotype - // set up the default 1-to-1 haplotype mapping object - final ArrayList> haplotypeMapping = new ArrayList>(); - for( final Haplotype h : haplotypes ) { - final ArrayList list = new ArrayList(); - list.add(h); - haplotypeMapping.add(list); - } - final double[][] haplotypeLikelihoodMatrix = computeDiploidHaplotypeLikelihoods( sampleKeySet, haplotypeMapping ); // all samples pooled together + final double[][] haplotypeLikelihoodMatrix = computeDiploidHaplotypeLikelihoods( sampleKeySet, haplotypes ); // all samples pooled together int hap1 = 0; int hap2 = 0; @@ -347,7 +363,7 @@ public class LikelihoodCalculationEngine { public static Map partitionReadsBasedOnLikelihoods( final GenomeLocParser parser, final HashMap> perSampleReadList, final HashMap> perSampleFilteredReadList, - final Pair>> call, + final Pair>> call, final double downsamplingFraction, final PrintStream downsamplingLog ) { final Map returnMap = new HashMap(); diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngineUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngineUnitTest.java index e82946690..19ced9f42 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngineUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngineUnitTest.java @@ -10,7 +10,6 @@ import org.broadinstitute.sting.BaseTest; import org.broadinstitute.sting.utils.Haplotype; import org.broadinstitute.sting.utils.MathUtils; import org.testng.Assert; -import org.testng.annotations.BeforeClass; import org.testng.annotations.DataProvider; import org.testng.annotations.Test; @@ -102,7 +101,9 @@ public class LikelihoodCalculationEngineUnitTest extends BaseTest { haplotypes.add(haplotype); } } - return LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(haplotypes, "myTestSample"); + final HashSet sampleSet = new HashSet(1); + sampleSet.add("myTestSample"); + return LikelihoodCalculationEngine.computeDiploidHaplotypeLikelihoods(sampleSet, haplotypes); } } diff --git a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java index b30d47074..6de15e18b 100755 --- a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java +++ b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java @@ -49,6 +49,7 @@ public class Haplotype { private int alignmentStartHapwrtRef; public int leftBreakPoint = 0; public int rightBreakPoint = 0; + private Allele artificialAllele = null; /** * Create a simple consensus sequence with provided bases and a uniform quality over all bases of qual @@ -71,6 +72,11 @@ public class Haplotype { this(bases, 0); } + public Haplotype( final byte[] bases, final Allele artificialAllele ) { + this(bases, 0); + this.artificialAllele = artificialAllele; + } + public Haplotype( final byte[] bases, final GenomeLoc loc ) { this(bases); this.genomeLocation = loc; @@ -171,6 +177,14 @@ public class Haplotype { this.cigar = cigar; } + public boolean isArtificialHaplotype() { + return artificialAllele != null; + } + + public Allele getArtificialAllele() { + return artificialAllele; + } + @Requires({"refInsertLocation >= 0"}) public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation ) { // refInsertLocation is in ref haplotype offset coordinates NOT genomic coordinates @@ -182,7 +196,7 @@ public class Haplotype { newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, 0, haplotypeInsertLocation)); // bases before the variant newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, altAllele.getBases()); // the alt allele of the variant newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, haplotypeInsertLocation + refAllele.length(), bases.length)); // bases after the variant - return new Haplotype(newHaplotypeBases); + return new Haplotype(newHaplotypeBases, altAllele); } public static class HaplotypeBaseComparator implements Comparator, Serializable { From f0b8a0228fef45f23478c1a12be1cd58633c0873 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Mon, 19 Nov 2012 09:57:55 -0500 Subject: [PATCH 3/9] Quick fix for HC refactoring: when copying over Haplotype objects, make sure to copy over the artificial allele used to create it too. --- .../gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java | 2 ++ public/java/src/org/broadinstitute/sting/utils/Haplotype.java | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java index fd46e4e69..33fa49543 100755 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java @@ -369,6 +369,8 @@ public class SimpleDeBruijnAssembler extends LocalAssemblyEngine { h.setAlignmentStartHapwrtRef( swConsensus2.getAlignmentStart2wrt1() ); h.setCigar( AlignmentUtils.leftAlignIndel(swConsensus2.getCigar(), ref, h.getBases(), swConsensus2.getAlignmentStart2wrt1(), 0) ); + if ( haplotype.isArtificialHaplotype() ) + h.setArtificialAllele(haplotype.getArtificialAllele()); h.leftBreakPoint = leftBreakPoint; h.rightBreakPoint = rightBreakPoint; if( swConsensus2.getCigar().toString().contains("S") || swConsensus2.getCigar().getReferenceLength() != activeRegionStop - activeRegionStart ) { // protect against SW failures diff --git a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java index 6de15e18b..af4e31698 100755 --- a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java +++ b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java @@ -185,6 +185,10 @@ public class Haplotype { return artificialAllele; } + public void setArtificialAllele(final Allele artificialAllele) { + this.artificialAllele = artificialAllele; + } + @Requires({"refInsertLocation >= 0"}) public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation ) { // refInsertLocation is in ref haplotype offset coordinates NOT genomic coordinates From 937ac7290f7caa7f0cae996608f5a68358f10c09 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 20 Nov 2012 16:13:29 -0500 Subject: [PATCH 4/9] Lots more GGA fixes for the HC now that I understand what's going on internally. Integration tests pass except for the GGA test which I believe now produces better results. --- .../haplotypecaller/GenotypingEngine.java | 82 ++++++++++++------- .../LikelihoodCalculationEngine.java | 4 +- .../SimpleDeBruijnAssembler.java | 22 +++-- .../HaplotypeCallerIntegrationTest.java | 3 +- .../broadinstitute/sting/utils/Haplotype.java | 17 ++-- .../sting/utils/HaplotypeUnitTest.java | 2 +- 6 files changed, 87 insertions(+), 43 deletions(-) diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java index 9fc636efe..beec8a92e 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java @@ -62,6 +62,7 @@ public class GenotypingEngine { final List activeAllelesToGenotype ) { final ArrayList>>> returnCalls = new ArrayList>>>(); + 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 final TreeSet startPosKeySet = new TreeSet(); @@ -70,7 +71,7 @@ public class GenotypingEngine { for( final Haplotype h : haplotypes ) { // Walk along the alignment and turn any difference from the reference into an event h.setEventMap( generateVCsFromAlignment( h, h.getAlignmentStartHapwrtRef(), h.getCigar(), ref, h.getBases(), refLoc, "HC" + count++ ) ); - if( activeAllelesToGenotype.isEmpty() ) { startPosKeySet.addAll(h.getEventMap().keySet()); } + if( !in_GGA_mode ) { startPosKeySet.addAll(h.getEventMap().keySet()); } if( DEBUG ) { System.out.println( h.toString() ); System.out.println( "> Cigar = " + h.getCigar() ); @@ -80,10 +81,10 @@ public class GenotypingEngine { } cleanUpSymbolicUnassembledEvents( haplotypes ); - if( activeAllelesToGenotype.isEmpty() && 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 + 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( !activeAllelesToGenotype.isEmpty() ) { // we are in GGA mode! + if( in_GGA_mode ) { for( final VariantContext compVC : activeAllelesToGenotype ) { startPosKeySet.add( compVC.getStart() ); } @@ -95,7 +96,7 @@ public class GenotypingEngine { final ArrayList eventsAtThisLoc = new ArrayList(); // the overlapping events to merge into a common reference view final ArrayList priorityList = new ArrayList(); // used to merge overlapping events into common reference view - if( activeAllelesToGenotype.isEmpty() ) { + if( !in_GGA_mode ) { for( final Haplotype h : haplotypes ) { final HashMap eventMap = h.getEventMap(); final VariantContext vc = eventMap.get(loc); @@ -129,9 +130,8 @@ public class GenotypingEngine { final Allele refAllele = eventsAtThisLoc.get(0).getReference(); final ArrayList alleleOrdering = new ArrayList(alleleMapper.size()); alleleOrdering.add(refAllele); - for ( final Allele allele : alleleMapper.keySet() ) { - if ( !refAllele.equals(allele) ) - alleleOrdering.add(allele); + for( final VariantContext vc : eventsAtThisLoc ) { + alleleOrdering.add(vc.getAlternateAllele(0)); } // Sanity check the priority list @@ -154,6 +154,16 @@ public class GenotypingEngine { final VariantContext mergedVC = VariantContextUtils.simpleMerge(genomeLocParser, eventsAtThisLoc, priorityList, VariantContextUtils.FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED, VariantContextUtils.GenotypeMergeType.PRIORITIZE, false, false, null, false, false); 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> updatedAlleleMapper = new HashMap>(alleleMapper.size()); + for ( int i = 0; i < mergedVC.getNAlleles(); i++ ) { + final Allele oldAllele = alleleOrdering.get(i); + final Allele newAllele = mergedVC.getAlleles().get(i); + updatedAlleleMapper.put(newAllele, alleleMapper.get(oldAllele)); + alleleOrdering.set(i, newAllele); + } + alleleMapper = updatedAlleleMapper; + if( DEBUG ) { System.out.println("Genotyping event at " + loc + " with alleles = " + mergedVC.getAlleles()); //System.out.println("Event/haplotype allele mapping = " + alleleMapper); @@ -358,48 +368,48 @@ public class GenotypingEngine { @Ensures({"result.size() == eventsAtThisLoc.size() + 1"}) protected static Map> createAlleleMapper( final int loc, final List eventsAtThisLoc, final List haplotypes ) { - final Allele refAllele = eventsAtThisLoc.get(0).getReference(); - final Map> alleleMapper = new HashMap>(eventsAtThisLoc.size()+1); + final Allele refAllele = eventsAtThisLoc.get(0).getReference(); + alleleMapper.put(refAllele, new ArrayList()); + for( final VariantContext vc : eventsAtThisLoc ) + alleleMapper.put(vc.getAlternateAllele(0), new ArrayList()); + + final ArrayList undeterminedHaplotypes = new ArrayList(haplotypes.size()); for( final Haplotype h : haplotypes ) { if( h.getEventMap().get(loc) == null ) { // no event at this location so this is a reference-supporting haplotype - if ( !alleleMapper.containsKey(refAllele) ) - alleleMapper.put(refAllele, new ArrayList()); alleleMapper.get(refAllele).add(h); - } else if ( h.isArtificialHaplotype() ) { - if ( !alleleMapper.containsKey(h.getArtificialAllele()) ) - alleleMapper.put(h.getArtificialAllele(), new ArrayList()); + } else if( h.isArtificialHaplotype() && loc == h.getArtificialAllelePosition() && alleleMapper.containsKey(h.getArtificialAllele()) ) { alleleMapper.get(h.getArtificialAllele()).add(h); } else { + boolean haplotypeIsDetermined = false; for( final VariantContext vcAtThisLoc : eventsAtThisLoc ) { if( h.getEventMap().get(loc).hasSameAllelesAs(vcAtThisLoc) ) { - final Allele altAllele = vcAtThisLoc.getAlternateAllele(0); - if ( !alleleMapper.containsKey(altAllele) ) - alleleMapper.put(altAllele, new ArrayList()); - alleleMapper.get(altAllele).add(h); + alleleMapper.get(vcAtThisLoc.getAlternateAllele(0)).add(h); + haplotypeIsDetermined = true; break; } } + + if( !haplotypeIsDetermined ) + undeterminedHaplotypes.add(h); } } - for( final Haplotype h : haplotypes ) { - if ( h.getEventMap().get(loc) == null || h.isArtificialHaplotype() ) - continue; - + for( final Haplotype h : undeterminedHaplotypes ) { Allele matchingAllele = null; - for ( final Map.Entry> alleleToTest : alleleMapper.entrySet() ) { - if ( alleleToTest.getKey().equals(refAllele) ) + for( final Map.Entry> alleleToTest : alleleMapper.entrySet() ) { + // don't test against the reference allele + if( alleleToTest.getKey().equals(refAllele) ) continue; final Haplotype artificialHaplotype = alleleToTest.getValue().get(0); - if ( isSubSetOf(artificialHaplotype.getEventMap(), h.getEventMap()) ) { + if( isSubSetOf(artificialHaplotype.getEventMap(), h.getEventMap(), true) ) { matchingAllele = alleleToTest.getKey(); break; } } - if ( matchingAllele == null ) + if( matchingAllele == null ) matchingAllele = refAllele; alleleMapper.get(matchingAllele).add(h); } @@ -407,20 +417,36 @@ public class GenotypingEngine { return alleleMapper; } - protected static boolean isSubSetOf(final Map subset, final Map superset) { + protected static boolean isSubSetOf(final Map subset, final Map superset, final boolean resolveSupersetToSubset) { for ( final Map.Entry fromSubset : subset.entrySet() ) { final VariantContext fromSuperset = superset.get(fromSubset.getKey()); if ( fromSuperset == null ) return false; - if ( !fromSuperset.hasAlternateAllele(fromSubset.getValue().getAlternateAllele(0)) ) + List supersetAlleles = fromSuperset.getAlternateAlleles(); + if ( resolveSupersetToSubset ) + supersetAlleles = resolveAlternateAlleles(fromSubset.getValue().getReference(), fromSuperset.getReference(), supersetAlleles); + + if ( !supersetAlleles.contains(fromSubset.getValue().getAlternateAllele(0)) ) return false; } return true; } + private static List resolveAlternateAlleles(final Allele targetReference, final Allele actualReference, final List currentAlleles) { + if ( targetReference.length() <= actualReference.length() ) + return currentAlleles; + + final List newAlleles = new ArrayList(currentAlleles.size()); + final byte[] extraBases = Arrays.copyOfRange(targetReference.getBases(), actualReference.length(), targetReference.length()); + for ( final Allele a : currentAlleles ) { + newAlleles.add(Allele.extend(a, extraBases)); + } + return newAlleles; + } + @Ensures({"result.size() == haplotypeAllelesForSample.size()"}) protected static List findEventAllelesInSample( final List eventAlleles, final List haplotypeAlleles, final List haplotypeAllelesForSample, final ArrayList> alleleMapper, final ArrayList haplotypes ) { if( haplotypeAllelesForSample.contains(Allele.NO_CALL) ) { return noCall; } diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java index 543987e74..304f8d5cb 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java @@ -196,8 +196,8 @@ public class LikelihoodCalculationEngine { } // 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()"}) + @Requires({"haplotypeList.size() > 0"}) + @Ensures({"result.length == result[0].length", "result.length == haplotypeList.size()"}) public static double[][] computeDiploidHaplotypeLikelihoods( final Set samples, final List haplotypeList ) { final int numHaplotypes = haplotypeList.size(); diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java index 33fa49543..4f072d720 100755 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java @@ -278,9 +278,10 @@ public class SimpleDeBruijnAssembler extends LocalAssemblyEngine { final int activeRegionStart = refHaplotype.getAlignmentStartHapwrtRef(); final int activeRegionStop = refHaplotype.getAlignmentStartHapwrtRef() + refHaplotype.getCigar().getReferenceLength(); - for( final VariantContext compVC : activeAllelesToGenotype ) { // for GGA mode, add the desired allele into the haplotype + // for GGA mode, add the desired allele into the haplotype + for( final VariantContext compVC : activeAllelesToGenotype ) { for( final Allele compAltAllele : compVC.getAlternateAlleles() ) { - final Haplotype insertedRefHaplotype = refHaplotype.insertAllele(compVC.getReference(), compAltAllele, activeRegionStart + compVC.getStart() - activeRegionWindow.getStart()); + final Haplotype insertedRefHaplotype = refHaplotype.insertAllele(compVC.getReference(), compAltAllele, activeRegionStart + compVC.getStart() - activeRegionWindow.getStart(), compVC.getStart()); if( !addHaplotype( insertedRefHaplotype, fullReferenceWithPadding, returnHaplotypes, activeRegionStart, activeRegionStop ) ) { return returnHaplotypes; //throw new ReviewedStingException("Unable to add reference+allele haplotype during GGA-enabled assembly: " + insertedRefHaplotype); @@ -290,15 +291,24 @@ public class SimpleDeBruijnAssembler extends LocalAssemblyEngine { for( final DefaultDirectedGraph graph : graphs ) { for ( final KBestPaths.Path path : KBestPaths.getKBestPaths(graph, NUM_BEST_PATHS_PER_KMER_GRAPH) ) { + final Haplotype h = new Haplotype( path.getBases( graph ), path.getScore() ); if( addHaplotype( h, fullReferenceWithPadding, returnHaplotypes, activeRegionStart, activeRegionStop ) ) { - if( !activeAllelesToGenotype.isEmpty() ) { // for GGA mode, add the desired allele into the haplotype if it isn't already present + + // for GGA mode, add the desired allele into the haplotype if it isn't already present + if( !activeAllelesToGenotype.isEmpty() ) { final HashMap eventMap = GenotypingEngine.generateVCsFromAlignment( h, h.getAlignmentStartHapwrtRef(), h.getCigar(), fullReferenceWithPadding, h.getBases(), refLoc, "HCassembly" ); // BUGBUG: need to put this function in a shared place for( final VariantContext compVC : activeAllelesToGenotype ) { // for GGA mode, add the desired allele into the haplotype if it isn't already present final VariantContext vcOnHaplotype = eventMap.get(compVC.getStart()); - if( vcOnHaplotype == null || !vcOnHaplotype.hasSameAllelesAs(compVC) ) { + + // This if statement used to additionally have: + // "|| !vcOnHaplotype.hasSameAllelesAs(compVC)" + // but that can lead to problems downstream when e.g. you are injecting a 1bp deletion onto + // a haplotype that already contains a 1bp insertion (so practically it is reference but + // falls into the bin for the 1bp deletion because we keep track of the artificial alleles). + if( vcOnHaplotype == null ) { for( final Allele compAltAllele : compVC.getAlternateAlleles() ) { - addHaplotype( h.insertAllele(compVC.getReference(), compAltAllele, activeRegionStart + compVC.getStart() - activeRegionWindow.getStart()), fullReferenceWithPadding, returnHaplotypes, activeRegionStart, activeRegionStop ); + addHaplotype( h.insertAllele(compVC.getReference(), compAltAllele, activeRegionStart + compVC.getStart() - activeRegionWindow.getStart(), compVC.getStart()), fullReferenceWithPadding, returnHaplotypes, activeRegionStart, activeRegionStop ); } } } @@ -370,7 +380,7 @@ public class SimpleDeBruijnAssembler extends LocalAssemblyEngine { h.setAlignmentStartHapwrtRef( swConsensus2.getAlignmentStart2wrt1() ); h.setCigar( AlignmentUtils.leftAlignIndel(swConsensus2.getCigar(), ref, h.getBases(), swConsensus2.getAlignmentStart2wrt1(), 0) ); if ( haplotype.isArtificialHaplotype() ) - h.setArtificialAllele(haplotype.getArtificialAllele()); + h.setArtificialAllele(haplotype.getArtificialAllele(), haplotype.getArtificialAllelePosition()); h.leftBreakPoint = leftBreakPoint; h.rightBreakPoint = rightBreakPoint; if( swConsensus2.getCigar().toString().contains("S") || swConsensus2.getCigar().getReferenceLength() != activeRegionStop - activeRegionStart ) { // protect against SW failures diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java index 6828dbcb5..a57462d1d 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java @@ -29,9 +29,10 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest { HCTest(NA12878_BAM, "", "baabae06c85d416920be434939124d7f"); } + // 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", "39da622b309597d7a0b082c8aa1748c9"); + HCTest(CEUTRIO_BAM, "--max_alternate_alleles 3 -gt_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf", "f2d0309fdf50d5827e9c60ed0dd07e3f"); } private void HCTestComplexVariants(String bam, String args, String md5) { diff --git a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java index af4e31698..30fdce75d 100755 --- a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java +++ b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java @@ -50,7 +50,8 @@ public class Haplotype { public int leftBreakPoint = 0; public int rightBreakPoint = 0; private Allele artificialAllele = null; - + private int artificialAllelePosition = -1; + /** * Create a simple consensus sequence with provided bases and a uniform quality over all bases of qual * @@ -72,9 +73,10 @@ public class Haplotype { this(bases, 0); } - public Haplotype( final byte[] bases, final Allele artificialAllele ) { + protected Haplotype( final byte[] bases, final Allele artificialAllele, final int artificialAllelePosition ) { this(bases, 0); this.artificialAllele = artificialAllele; + this.artificialAllelePosition = artificialAllelePosition; } public Haplotype( final byte[] bases, final GenomeLoc loc ) { @@ -185,12 +187,17 @@ public class Haplotype { return artificialAllele; } - public void setArtificialAllele(final Allele artificialAllele) { + public int getArtificialAllelePosition() { + return artificialAllelePosition; + } + + public void setArtificialAllele(final Allele artificialAllele, final int artificialAllelePosition) { this.artificialAllele = artificialAllele; + this.artificialAllelePosition = artificialAllelePosition; } @Requires({"refInsertLocation >= 0"}) - public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation ) { + public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation, final int genomicInsertLocation ) { // refInsertLocation is in ref haplotype offset coordinates NOT genomic coordinates final int haplotypeInsertLocation = ReadUtils.getReadCoordinateForReferenceCoordinate(alignmentStartHapwrtRef, cigar, refInsertLocation, ReadUtils.ClippingTail.RIGHT_TAIL, true); if( haplotypeInsertLocation == -1 || haplotypeInsertLocation + refAllele.length() >= bases.length ) { // desired change falls inside deletion so don't bother creating a new haplotype @@ -200,7 +207,7 @@ public class Haplotype { newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, 0, haplotypeInsertLocation)); // bases before the variant newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, altAllele.getBases()); // the alt allele of the variant newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, haplotypeInsertLocation + refAllele.length(), bases.length)); // bases after the variant - return new Haplotype(newHaplotypeBases, altAllele); + return new Haplotype(newHaplotypeBases, altAllele, genomicInsertLocation); } public static class HaplotypeBaseComparator implements Comparator, Serializable { diff --git a/public/java/test/org/broadinstitute/sting/utils/HaplotypeUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/HaplotypeUnitTest.java index ddffb6e4c..13db1d39e 100644 --- a/public/java/test/org/broadinstitute/sting/utils/HaplotypeUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/HaplotypeUnitTest.java @@ -159,7 +159,7 @@ public class HaplotypeUnitTest extends BaseTest { final VariantContext vc = new VariantContextBuilder().alleles(alleles).loc("1", loc, loc + h1refAllele.getBases().length - 1).make(); h.setAlignmentStartHapwrtRef(0); h.setCigar(cigar); - final Haplotype h1 = h.insertAllele(vc.getReference(), vc.getAlternateAllele(0), loc); + final Haplotype h1 = h.insertAllele(vc.getReference(), vc.getAlternateAllele(0), loc, vc.getStart()); final Haplotype h1expected = new Haplotype(newHap.getBytes()); Assert.assertEquals(h1, h1expected); } From cc7680e6010ab184752c9f4c20b96e055178478c Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Sun, 4 Nov 2012 14:40:17 -0800 Subject: [PATCH 5/9] NA12878 knowledge base backed by MongoDB -- Idea is simply to create a persistent database of all TP/FP sites on chr20 in NA12878. Individual callsets can be imported, and a consensus algorithm is run over all callsets in the database to create a consensus collection, which can be used to assess NA12878 callsets for GATK and methods development -- Framework for representing simple VariantContexts and Genotypes in MongoDB, querying for records, and iterating over them in the GATK -- Not hooked up to Tribble, but could be done reasonably easily now (future TODO) -- Tools to import callsets, create consensus callsets, import and export reviews -- Scripts to reset the knowledge base and repopulate it with the standard data files (Eric will expand) -- Actually scales to all of chr20, includes AssessNA12878 that reads a VCF and itemizes it against the truth data set -- ImportCallset can load OMNI, HM3, CEU best practices, mills/devine sites and genotypes, properly marking sites as poly/mono/unk as well as TP/FP/UNK based on command line parameters -- Added shell scripts that start up a local mongo db, that connect to a local or BI hosted mongo for NA12878.db for debugging, and a setupNA12878db script that can load OMNI, HM3, CEU best practices, Mills/Devine into the db and then update the consensus. -- Reviewed sites can be exported to a VCF, and imported again, as a mechanism to safely store the only non-recoverable data from the Mongo DB. -- Created a NA12878DBWalker that manages the outer DB interaction, and that all MongoDB interacting walkers inherit from. Added a NA12878DBArgumentCollection.java consolating all of the common command line arguments (though strictly not necessary as all of this occurs in the root walker) UnitTests -- Can connect to a test knowledge base for development and unit testing -- PolymorphicStatus, TruthStatus, SiteIterator -- NA12878KBUnitTestBase provides simple utilities for connecting to the test mongo db, getting calls, etc -- MongoVariantContext tests creation, matching, and encoding -> writing -> read -> decoding from the mongodb AssessNA12878 -- Generic tool for comparing a NA12878 callset against the knowledge base. See http://gatkforums.broadinstitute.org/discussion/1848/using-the-na12878-knowledge-base for detailed documentation -- Performs trivial filtering on FS, MQ, QD for SNPs and non-SNPs to separate out variants likely to be filtered from those that are honest-to-goodness FPs Misc -- Ability to provide Description for Simplified GATK report --- .../sting/gatk/report/GATKReport.java | 13 ++++++- .../broadinstitute/sting/utils/GenomeLoc.java | 14 ++++++++ .../org/broadinstitute/sting/utils/Utils.java | 4 +++ .../sting/utils/codecs/vcf/VCFUtils.java | 34 +++++++++++++++++++ 4 files changed, 64 insertions(+), 1 deletion(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/report/GATKReport.java b/public/java/src/org/broadinstitute/sting/gatk/report/GATKReport.java index 47bc48f81..6685ee12a 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/report/GATKReport.java +++ b/public/java/src/org/broadinstitute/sting/gatk/report/GATKReport.java @@ -271,7 +271,18 @@ public class GATKReport { * @return a simplified GATK report */ public static GATKReport newSimpleReport(final String tableName, final String... columns) { - GATKReportTable table = new GATKReportTable(tableName, "A simplified GATK table report", columns.length); + return newSimpleReportWithDescription(tableName, "A simplified GATK table report", columns); + } + + /** + * @see #newSimpleReport(String, String...) but with a customized description + * @param tableName + * @param desc + * @param columns + * @return + */ + public static GATKReport newSimpleReportWithDescription(final String tableName, final String desc, final String... columns) { + GATKReportTable table = new GATKReportTable(tableName, desc, columns.length); for (String column : columns) { table.addColumn(column, ""); diff --git a/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java b/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java index 6df9c9f1d..4d2c26a79 100644 --- a/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java +++ b/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java @@ -315,6 +315,20 @@ public class GenomeLoc implements Comparable, Serializable, HasGenome return ( comparison == -1 || ( comparison == 0 && this.getStop() < that.getStart() )); } + /** + * Tests whether this genome loc starts at the same position as that. + * + * i.e., do this and that have the same contig and the same start position + * + * @param that genome loc to compare to + * @return true if this and that have the same contig and the same start position + */ + @Requires("that != null") + public final boolean startsAt( GenomeLoc that ) { + int comparison = this.compareContigs(that); + return comparison == 0 && this.getStart() == that.getStart(); + } + /** * Tests whether any portion of this contig is before that contig. * @param that Other contig to test. diff --git a/public/java/src/org/broadinstitute/sting/utils/Utils.java b/public/java/src/org/broadinstitute/sting/utils/Utils.java index b780d0966..1d12d6f8b 100755 --- a/public/java/src/org/broadinstitute/sting/utils/Utils.java +++ b/public/java/src/org/broadinstitute/sting/utils/Utils.java @@ -293,6 +293,10 @@ public class Utils { } } + public static String join(final String separator, final T ... objects) { + return join(separator, Arrays.asList(objects)); + } + public static String dupString(char c, int nCopies) { char[] chars = new char[nCopies]; Arrays.fill(chars, c); diff --git a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java index be87e7306..a8aefb703 100755 --- a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java @@ -30,12 +30,17 @@ import net.sf.samtools.SAMSequenceRecord; import org.apache.commons.io.FilenameUtils; import org.apache.log4j.Logger; import org.broad.tribble.Feature; +import org.broad.tribble.FeatureCodecHeader; +import org.broad.tribble.readers.PositionalBufferedStream; import org.broadinstitute.sting.commandline.RodBinding; import org.broadinstitute.sting.gatk.GenomeAnalysisEngine; import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource; +import org.broadinstitute.sting.utils.collections.Pair; import org.broadinstitute.sting.utils.variantcontext.VariantContext; import java.io.File; +import java.io.FileInputStream; +import java.io.IOException; import java.util.*; /** @@ -317,4 +322,33 @@ public class VCFUtils { assembly = "hg19"; return assembly; } + + /** + * Read all of the VCF records from source into memory, returning the header and the VariantContexts + * + * @param source the file to read, must be in VCF4 format + * @return + * @throws IOException + */ + public static Pair> readVCF(final File source) throws IOException { + // read in the features + final List vcs = new ArrayList(); + final VCFCodec codec = new VCFCodec(); + PositionalBufferedStream pbs = new PositionalBufferedStream(new FileInputStream(source)); + FeatureCodecHeader header = codec.readHeader(pbs); + pbs.close(); + + pbs = new PositionalBufferedStream(new FileInputStream(source)); + pbs.skip(header.getHeaderEnd()); + + final VCFHeader vcfHeader = (VCFHeader)header.getHeaderValue(); + + while ( ! pbs.isDone() ) { + final VariantContext vc = codec.decode(pbs); + if ( vc != null ) + vcs.add(vc); + } + + return new Pair>(vcfHeader, vcs); + } } \ No newline at end of file From ff87642a91c569982fed062d34cf1bf63caab63b Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 20 Nov 2012 22:29:56 -0500 Subject: [PATCH 6/9] Enable cycle covariate unit tests --- .../sting/utils/recalibration/CycleCovariateUnitTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java index c3d93b2cb..deb0931d6 100644 --- a/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java @@ -24,7 +24,7 @@ public class CycleCovariateUnitTest { covariate.initialize(RAC); } - @Test(enabled = false) + @Test(enabled = true) public void testSimpleCycles() { short readLength = 10; GATKSAMRecord read = ReadUtils.createRandomRead(readLength); From 72e2d569c540e758284d9304caf44bf2b2e6ca35 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 20 Nov 2012 22:41:57 -0500 Subject: [PATCH 7/9] The user can now set the maximum allowable cycle on the command-line with --maximum_cycle_value. This value is (now) enforced in the Cycle covariate and a User Error is thrown if the maximum value is passed (with a helpful error message). Added unit tests to cover this new functionality. --- .../bqsr/RecalibrationArgumentCollection.java | 18 ++++++++++--- .../covariates/CycleCovariate.java | 7 +++++- .../recalibration/CycleCovariateUnitTest.java | 25 ++++++++++++++++++- 3 files changed, 44 insertions(+), 6 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/RecalibrationArgumentCollection.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/RecalibrationArgumentCollection.java index e5704a1e2..c64482151 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/RecalibrationArgumentCollection.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/RecalibrationArgumentCollection.java @@ -102,13 +102,10 @@ public class RecalibrationArgumentCollection { @Argument(fullName = "no_standard_covs", shortName = "noStandard", doc = "Do not use the standard set of covariates, but rather just the ones listed using the -cov argument", required = false) public boolean DO_NOT_USE_STANDARD_COVARIATES = false; - ///////////////////////////// - // Debugging-only Arguments - ///////////////////////////// /** * This calculation is critically dependent on being able to skip over known polymorphic sites. Please be sure that you know what you are doing if you use this option. */ - @Hidden + @Advanced @Argument(fullName = "run_without_dbsnp_potentially_ruining_quality", shortName = "run_without_dbsnp_potentially_ruining_quality", required = false, doc = "If specified, allows the recalibrator to be used without a dbsnp rod. Very unsafe and for expert users only.") public boolean RUN_WITHOUT_DBSNP = false; @@ -139,6 +136,13 @@ public class RecalibrationArgumentCollection { @Argument(fullName = "indels_context_size", shortName = "ics", doc = "size of the k-mer context to be used for base insertions and deletions", required = false) public int INDELS_CONTEXT_SIZE = 3; + /** + * The cycle covariate will generate an error if it encounters a cycle greater than this value. + * This argument is ignored if the Cycle covariate is not used. + */ + @Argument(fullName = "maximum_cycle_value", shortName = "maxCycle", doc = "the maximum cycle value permitted for the Cycle covariate", required = false) + public int MAXIMUM_CYCLE_VALUE = 500; + /** * A default base qualities to use as a prior (reported quality) in the mismatch covariate model. This value will replace all base qualities in the read for this default value. Negative value turns it off (default is off) */ @@ -176,9 +180,15 @@ public class RecalibrationArgumentCollection { @Argument(fullName = "binary_tag_name", shortName = "bintag", required = false, doc = "the binary tag covariate name if using it") public String BINARY_TAG_NAME = null; + + ///////////////////////////// + // Debugging-only Arguments + ///////////////////////////// + @Hidden @Argument(fullName = "default_platform", shortName = "dP", required = false, doc = "If a read has no platform then default to the provided String. Valid options are illumina, 454, and solid.") public String DEFAULT_PLATFORM = null; + @Hidden @Argument(fullName = "force_platform", shortName = "fP", required = false, doc = "If provided, the platform of EVERY read will be forced to be the provided String. Valid options are illumina, 454, and solid.") public String FORCE_PLATFORM = null; diff --git a/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/CycleCovariate.java b/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/CycleCovariate.java index 5d0d94b69..a9b6c7152 100755 --- a/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/CycleCovariate.java +++ b/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/CycleCovariate.java @@ -49,7 +49,7 @@ import java.util.EnumSet; public class CycleCovariate implements StandardCovariate { - private static final int MAXIMUM_CYCLE_VALUE = 1000; + private int MAXIMUM_CYCLE_VALUE; private static final int CUSHION_FOR_INDELS = 4; private String default_platform = null; @@ -59,6 +59,8 @@ public class CycleCovariate implements StandardCovariate { // Initialize any member variables using the command-line arguments passed to the walkers @Override public void initialize(final RecalibrationArgumentCollection RAC) { + this.MAXIMUM_CYCLE_VALUE = RAC.MAXIMUM_CYCLE_VALUE; + if (RAC.DEFAULT_PLATFORM != null && !NGSPlatform.isKnown(RAC.DEFAULT_PLATFORM)) throw new UserException.CommandLineException("The requested default platform (" + RAC.DEFAULT_PLATFORM + ") is not a recognized platform."); @@ -88,6 +90,9 @@ public class CycleCovariate implements StandardCovariate { final int MAX_CYCLE_FOR_INDELS = readLength - CUSHION_FOR_INDELS - 1; for (int i = 0; i < readLength; i++) { + if ( cycle > MAXIMUM_CYCLE_VALUE ) + throw new UserException("The maximum allowed value for the cycle is " + MAXIMUM_CYCLE_VALUE + ", but a larger cycle was detected in read " + read.getReadName() + ". Please use the --maximum_cycle_value argument to increase this value (at the expense of requiring more memory to run)"); + final int substitutionKey = keyFromCycle(cycle); final int indelKey = (i < CUSHION_FOR_INDELS || i > MAX_CYCLE_FOR_INDELS) ? -1 : substitutionKey; values.addCovariate(substitutionKey, indelKey, indelKey, i); diff --git a/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java index deb0931d6..b73b1a311 100644 --- a/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/recalibration/CycleCovariateUnitTest.java @@ -1,6 +1,7 @@ package org.broadinstitute.sting.utils.recalibration; import org.broadinstitute.sting.gatk.walkers.bqsr.RecalibrationArgumentCollection; +import org.broadinstitute.sting.utils.exceptions.UserException; import org.broadinstitute.sting.utils.recalibration.covariates.CycleCovariate; import org.broadinstitute.sting.utils.sam.GATKSAMReadGroupRecord; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; @@ -53,9 +54,31 @@ public class CycleCovariateUnitTest { for (short i = 0; i < values.length; i++) { short actual = Short.decode(covariate.formatKey(values[i][0])); int expected = init + (increment * i); - // System.out.println(String.format("%d: %d, %d", i, actual, expected)); Assert.assertEquals(actual, expected); } } + @Test(enabled = true, expectedExceptions={UserException.class}) + public void testMoreThanMaxCycleFails() { + int readLength = RAC.MAXIMUM_CYCLE_VALUE + 1; + GATKSAMRecord read = ReadUtils.createRandomRead(readLength); + read.setReadPairedFlag(true); + read.setReadGroup(new GATKSAMReadGroupRecord("MY.ID")); + read.getReadGroup().setPlatform("illumina"); + + ReadCovariates readCovariates = new ReadCovariates(read.getReadLength(), 1); + covariate.recordValues(read, readCovariates); + } + + @Test(enabled = true) + public void testMaxCyclePasses() { + int readLength = RAC.MAXIMUM_CYCLE_VALUE; + GATKSAMRecord read = ReadUtils.createRandomRead(readLength); + read.setReadPairedFlag(true); + read.setReadGroup(new GATKSAMReadGroupRecord("MY.ID")); + read.getReadGroup().setPlatform("illumina"); + + ReadCovariates readCovariates = new ReadCovariates(read.getReadLength(), 1); + covariate.recordValues(read, readCovariates); + } } From c54fc94505e6d1b913081b8854d5063d5ed11593 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 20 Nov 2012 23:19:59 -0500 Subject: [PATCH 8/9] Protect against features that start off the end of the read (otherwise, Arrays.fill fails) --- .../gatk/walkers/bqsr/BaseRecalibrator.java | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/BaseRecalibrator.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/BaseRecalibrator.java index 9506510a9..b415bb1f5 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/BaseRecalibrator.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/bqsr/BaseRecalibrator.java @@ -268,16 +268,25 @@ public class BaseRecalibrator extends ReadWalker implements NanoSche } protected boolean[] calculateKnownSites( final GATKSAMRecord read, final List features ) { - final int BUFFER_SIZE = 0; final int readLength = read.getReadBases().length; final boolean[] knownSites = new boolean[readLength]; Arrays.fill(knownSites, false); for( final Feature f : features ) { int featureStartOnRead = ReadUtils.getReadCoordinateForReferenceCoordinate(read.getSoftStart(), read.getCigar(), f.getStart(), ReadUtils.ClippingTail.LEFT_TAIL, true); // BUGBUG: should I use LEFT_TAIL here? - if( featureStartOnRead == ReadUtils.CLIPPING_GOAL_NOT_REACHED ) { featureStartOnRead = 0; } + if( featureStartOnRead == ReadUtils.CLIPPING_GOAL_NOT_REACHED ) { + featureStartOnRead = 0; + } + int featureEndOnRead = ReadUtils.getReadCoordinateForReferenceCoordinate(read.getSoftStart(), read.getCigar(), f.getEnd(), ReadUtils.ClippingTail.LEFT_TAIL, true); - if( featureEndOnRead == ReadUtils.CLIPPING_GOAL_NOT_REACHED ) { featureEndOnRead = readLength; } - Arrays.fill(knownSites, Math.max(0, featureStartOnRead - BUFFER_SIZE), Math.min(readLength, featureEndOnRead + 1 + BUFFER_SIZE), true); + if( featureEndOnRead == ReadUtils.CLIPPING_GOAL_NOT_REACHED ) { + featureEndOnRead = readLength; + } + + if( featureStartOnRead > readLength ) { + featureStartOnRead = featureEndOnRead = readLength; + } + + Arrays.fill(knownSites, Math.max(0, featureStartOnRead), Math.min(readLength, featureEndOnRead + 1), true); } return knownSites; }