From 62cee2fb5bf4e4c11e9317c3e8bc86525793b30d Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Mon, 18 Jun 2012 21:36:27 -0400 Subject: [PATCH] Feature request from Tim that could be useful to all: there's now an --interval_padding argument that specifies how many basepairs to add to each of the intervals provided with -L (on both ends). This is particularly useful when trying to run over the exome plus flanks and don't want to have to pre-compute the flanks (just use e.g. --interval_padding 50). Added integration test to cover this feature. --- .../sting/commandline/IntervalBinding.java | 17 ++++++++------- .../sting/gatk/GenomeAnalysisEngine.java | 21 ++++++++++++++++--- .../arguments/GATKArgumentCollection.java | 6 ++++++ .../sting/utils/GenomeLocParser.java | 21 +++++++++++++++++++ .../sting/utils/interval/IntervalUtils.java | 20 ++++++++++++++++++ .../interval/IntervalIntegrationTest.java | 19 +++++++++++++++-- 6 files changed, 91 insertions(+), 13 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/commandline/IntervalBinding.java b/public/java/src/org/broadinstitute/sting/commandline/IntervalBinding.java index 1f873ffbd..3ad84c9a5 100644 --- a/public/java/src/org/broadinstitute/sting/commandline/IntervalBinding.java +++ b/public/java/src/org/broadinstitute/sting/commandline/IntervalBinding.java @@ -29,17 +29,14 @@ import org.broad.tribble.AbstractFeatureReader; import org.broad.tribble.Feature; import org.broad.tribble.FeatureCodec; import org.broad.tribble.FeatureReader; -import org.broad.tribble.readers.AsciiLineReader; import org.broadinstitute.sting.gatk.GenomeAnalysisEngine; import org.broadinstitute.sting.gatk.refdata.ReferenceDependentFeatureCodec; import org.broadinstitute.sting.gatk.refdata.tracks.FeatureManager; import org.broadinstitute.sting.utils.GenomeLoc; +import org.broadinstitute.sting.utils.GenomeLocParser; import org.broadinstitute.sting.utils.exceptions.UserException; import org.broadinstitute.sting.utils.interval.IntervalUtils; -import java.io.File; -import java.io.FileInputStream; -import java.io.IOException; import java.util.*; /** @@ -73,7 +70,11 @@ public final class IntervalBinding { return stringIntervals; } - public List getIntervals(GenomeAnalysisEngine toolkit) { + public List getIntervals(final GenomeAnalysisEngine toolkit) { + return getIntervals(toolkit.getGenomeLocParser()); + } + + public List getIntervals(final GenomeLocParser genomeLocParser) { List intervals; if ( featureIntervals != null ) { @@ -83,17 +84,17 @@ public final class IntervalBinding { final FeatureCodec codec = new FeatureManager().getByName(featureIntervals.getTribbleType()).getCodec(); if ( codec instanceof ReferenceDependentFeatureCodec ) - ((ReferenceDependentFeatureCodec)codec).setGenomeLocParser(toolkit.getGenomeLocParser()); + ((ReferenceDependentFeatureCodec)codec).setGenomeLocParser(genomeLocParser); try { FeatureReader reader = AbstractFeatureReader.getFeatureReader(featureIntervals.getSource(), codec, false); for ( Feature feature : reader.iterator() ) - intervals.add(toolkit.getGenomeLocParser().createGenomeLoc(feature)); + intervals.add(genomeLocParser.createGenomeLoc(feature)); } catch (Exception e) { throw new UserException.MalformedFile(featureIntervals.getSource(), "Problem reading the interval file", e); } } else { - intervals = IntervalUtils.parseIntervalArguments(toolkit.getGenomeLocParser(), stringIntervals); + intervals = IntervalUtils.parseIntervalArguments(genomeLocParser, stringIntervals); } return intervals; diff --git a/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java b/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java index d77f266bd..6fa70f437 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java @@ -592,7 +592,7 @@ public class GenomeAnalysisEngine { // if include argument isn't given, create new set of all possible intervals GenomeLocSortedSet includeSortedSet = (argCollection.intervals == null ? GenomeLocSortedSet.createSetFromSequenceDictionary(this.referenceDataSource.getReference().getSequenceDictionary()) : - loadIntervals(argCollection.intervals, argCollection.intervalSetRule)); + loadIntervals(argCollection.intervals, argCollection.intervalSetRule, argCollection.intervalPadding)); // if no exclude arguments, can return parseIntervalArguments directly if ( argCollection.excludeIntervals == null ) @@ -619,16 +619,31 @@ public class GenomeAnalysisEngine { * @param rule interval merging rule * @return A sorted, merged list of all intervals specified in this arg list. */ - protected GenomeLocSortedSet loadIntervals( List> argList, IntervalSetRule rule ) { + protected GenomeLocSortedSet loadIntervals( final List> argList, final IntervalSetRule rule ) { + return loadIntervals(argList, rule, 0); + } + + /** + * Loads the intervals relevant to the current execution + * @param argList argument bindings; might include filenames, intervals in samtools notation, or a combination of the above + * @param rule interval merging rule + * @param padding how much to pad the intervals + * @return A sorted, merged list of all intervals specified in this arg list. + */ + protected GenomeLocSortedSet loadIntervals( final List> argList, final IntervalSetRule rule, final int padding ) { List allIntervals = new ArrayList(); for ( IntervalBinding intervalBinding : argList ) { - List intervals = intervalBinding.getIntervals(this); + List intervals = intervalBinding.getIntervals(this.getGenomeLocParser()); if ( intervals.isEmpty() ) { logger.warn("The interval file " + intervalBinding.getSource() + " contains no intervals that could be parsed."); } + if ( padding > 0 ) { + intervals = IntervalUtils.getIntervalsWithFlanks(this.getGenomeLocParser(), intervals, padding); + } + allIntervals = IntervalUtils.mergeListsBySetOperator(intervals, allIntervals, rule); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java b/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java index bd830a920..13c737a2e 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java +++ b/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java @@ -104,6 +104,12 @@ public class GATKArgumentCollection { @Argument(fullName = "interval_merging", shortName = "im", doc = "Indicates the interval merging rule we should use for abutting intervals", required = false) public IntervalMergingRule intervalMerging = IntervalMergingRule.ALL; + /** + * For example, '-L chr1:100' with a padding value of 20 would turn into '-L chr1:80-120'. + */ + @Argument(fullName = "interval_padding", shortName = "ip", doc = "Indicates how many basepairs of padding to include around each of the intervals specified with the -L/--intervals argument", required = false) + public int intervalPadding = 0; + @Input(fullName = "reference_sequence", shortName = "R", doc = "Reference sequence file", required = false) public File referenceFile = null; diff --git a/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java b/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java index e33ff6f30..b3f4542af 100644 --- a/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java +++ b/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java @@ -587,6 +587,27 @@ public final class GenomeLocParser { return createGenomeLoc(contigName, contigIndex, start, stop, true); } + /** + * Creates a loc padded in both directions by maxBasePairs size (if possible). + * @param loc The original loc + * @param padding The number of base pairs to pad on either end + * @return The contiguous loc of length up to the original length + 2*padding (depending on the start/end of the contig). + */ + @Requires({"loc != null", "padding > 0"}) + public GenomeLoc createPaddedGenomeLoc(final GenomeLoc loc, final int padding) { + if (GenomeLoc.isUnmapped(loc)) + return loc; + final String contigName = loc.getContig(); + final SAMSequenceRecord contig = contigInfo.getSequence(contigName); + final int contigIndex = contig.getSequenceIndex(); + final int contigLength = contig.getSequenceLength(); + + final int start = Math.max(1, loc.getStart() - padding); + final int stop = Math.min(contigLength, loc.getStop() + padding); + + return createGenomeLoc(contigName, contigIndex, start, stop, true); + } + /** * Creates a loc to the right (starting at the loc stop + 1) of maxBasePairs size. * @param loc The original loc diff --git a/public/java/src/org/broadinstitute/sting/utils/interval/IntervalUtils.java b/public/java/src/org/broadinstitute/sting/utils/interval/IntervalUtils.java index ea1eaeb51..c96226405 100644 --- a/public/java/src/org/broadinstitute/sting/utils/interval/IntervalUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/interval/IntervalUtils.java @@ -648,6 +648,26 @@ public class IntervalUtils { return expanded; } + /** + * Returns a list of intervals between the passed int locs. Does not extend UNMAPPED locs. + * @param parser A genome loc parser for creating the new intervals + * @param locs Original genome locs + * @param basePairs Number of base pairs on each side of loc + * @return The list of intervals between the locs + */ + public static List getIntervalsWithFlanks(final GenomeLocParser parser, final List locs, final int basePairs) { + + if (locs.size() == 0) + return Collections.emptyList(); + + final List expanded = new ArrayList(); + for ( final GenomeLoc loc : locs ) { + expanded.add(parser.createPaddedGenomeLoc(loc, basePairs)); + } + + return sortAndMergeIntervals(parser, expanded, IntervalMergingRule.ALL).toList(); + } + private static LinkedHashMap> splitByContig(List sorted) { LinkedHashMap> splits = new LinkedHashMap>(); GenomeLoc last = null; diff --git a/public/java/test/org/broadinstitute/sting/utils/interval/IntervalIntegrationTest.java b/public/java/test/org/broadinstitute/sting/utils/interval/IntervalIntegrationTest.java index 4d8e096c1..875765ca8 100644 --- a/public/java/test/org/broadinstitute/sting/utils/interval/IntervalIntegrationTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/interval/IntervalIntegrationTest.java @@ -207,11 +207,26 @@ public class IntervalIntegrationTest extends WalkerTest { " -R " + hg18Reference + " -o %s" + " -L " + validationDataLocation + "intervalTest.3.vcf", - 1, // just one output file - Arrays.asList(md5)); + 1, // just one output file + Arrays.asList(md5)); executeTest("testComplexVCF", spec); } + @Test(enabled = true) + public void testComplexVCFWithPadding() { + String md5 = "649ee93d50739c656e94ec88a32c7ffe"; + WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec( + "-T CountLoci" + + " --interval_padding 2" + + " -I " + validationDataLocation + "OV-0930.normal.chunk.bam" + + " -R " + hg18Reference + + " -o %s" + + " -L " + validationDataLocation + "intervalTest.3.vcf", + 1, // just one output file + Arrays.asList(md5)); + executeTest("testComplexVCFWithPadding", spec); + } + @Test(enabled = true) public void testMergingWithComplexVCF() { String md5 = "6d7fce9fee471194aa8b5b6e47267f03";