Search in sources :

Example 1 with Interval

use of htsjdk.samtools.util.Interval in project gatk by broadinstitute.

the class SplitIntervals method onTraversalStart.

@Override
public void onTraversalStart() {
    ParamUtils.isPositive(scatterCount, "scatter count must be > 0.");
    if (!outputDir.exists() && !outputDir.mkdir()) {
        throw new RuntimeIOException("Unable to create directory: " + outputDir.getAbsolutePath());
    }
    // in general dictionary will be from the reference, but using -I reads.bam or -F variants.vcf
    // to use the sequence dict from a bam or vcf is also supported
    final SAMSequenceDictionary sequenceDictionary = getBestAvailableSequenceDictionary();
    final List<SimpleInterval> intervals = hasIntervals() ? intervalArgumentCollection.getIntervals(sequenceDictionary) : IntervalUtils.getAllIntervalsForReference(sequenceDictionary);
    final IntervalList intervalList = new IntervalList(sequenceDictionary);
    intervals.stream().map(si -> new Interval(si.getContig(), si.getStart(), si.getEnd())).forEach(intervalList::add);
    final IntervalListScatterer scatterer = new IntervalListScatterer(subdivisionMode);
    final List<IntervalList> scattered = scatterer.scatter(intervalList, scatterCount, false);
    final DecimalFormat formatter = new DecimalFormat("0000");
    IntStream.range(0, scattered.size()).forEach(n -> scattered.get(n).write(new File(outputDir, formatter.format(n) + "-scattered.intervals")));
}
Also used : IntStream(java.util.stream.IntStream) CommandLineProgramProperties(org.broadinstitute.barclay.argparser.CommandLineProgramProperties) IntervalListScatterer(org.broadinstitute.hellbender.tools.picard.interval.IntervalListScatterer) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Argument(org.broadinstitute.barclay.argparser.Argument) DecimalFormat(java.text.DecimalFormat) IntervalList(htsjdk.samtools.util.IntervalList) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) VariantProgramGroup(org.broadinstitute.hellbender.cmdline.programgroups.VariantProgramGroup) IntervalArgumentCollection(org.broadinstitute.hellbender.cmdline.argumentcollections.IntervalArgumentCollection) ParamUtils(org.broadinstitute.hellbender.utils.param.ParamUtils) File(java.io.File) GATKTool(org.broadinstitute.hellbender.engine.GATKTool) Interval(htsjdk.samtools.util.Interval) List(java.util.List) IntervalUtils(org.broadinstitute.hellbender.utils.IntervalUtils) RuntimeIOException(htsjdk.samtools.util.RuntimeIOException) RuntimeIOException(htsjdk.samtools.util.RuntimeIOException) IntervalListScatterer(org.broadinstitute.hellbender.tools.picard.interval.IntervalListScatterer) IntervalList(htsjdk.samtools.util.IntervalList) DecimalFormat(java.text.DecimalFormat) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) File(java.io.File) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Interval(htsjdk.samtools.util.Interval)

Example 2 with Interval

use of htsjdk.samtools.util.Interval in project gatk by broadinstitute.

the class CollectRnaSeqMetricsTest method testMultiLevel.

@Test
public void testMultiLevel() throws Exception {
    final String sequence = "chr1";
    final String ignoredSequence = "chrM";
    // Create some alignments that hit the ribosomal sequence, various parts of the gene, and intergenic.
    final SAMRecordSetBuilder builder = new SAMRecordSetBuilder(true, SAMFileHeader.SortOrder.coordinate, false);
    // Set seed so that strandedness is consistent among runs.
    builder.setRandomSeed(0);
    final int sequenceIndex = builder.getHeader().getSequenceIndex(sequence);
    final SAMReadGroupRecord rg1 = new SAMReadGroupRecord("2");
    rg1.setSample("Sample");
    rg1.setLibrary("foo");
    builder.setReadGroup(rg1);
    builder.addPair("pair1", sequenceIndex, 45, 475);
    builder.addPair("pair2", sequenceIndex, 90, 225);
    builder.addFrag("frag1", sequenceIndex, 150, true);
    builder.addFrag("frag2", sequenceIndex, 450, true);
    final SAMReadGroupRecord rg2 = new SAMReadGroupRecord("3");
    rg2.setSample("Sample");
    rg2.setLibrary("bar");
    builder.setReadGroup(rg2);
    builder.addPair("pair3", sequenceIndex, 120, 600);
    builder.addFrag("frag3", sequenceIndex, 225, false);
    builder.addPair("rrnaPair", sequenceIndex, 400, 500);
    builder.addFrag("ignoredFrag", builder.getHeader().getSequenceIndex(ignoredSequence), 1, false);
    final File samFile = BaseTest.createTempFile("tmp.collectRnaSeqMetrics.", ".sam");
    try (final SAMFileWriter samWriter = new SAMFileWriterFactory().makeSAMWriter(builder.getHeader(), false, samFile)) {
        for (final SAMRecord rec : builder.getRecords()) samWriter.addAlignment(rec);
    }
    // Create an interval list with one ribosomal interval.
    final Interval rRnaInterval = new Interval(sequence, 300, 520, true, "rRNA");
    final IntervalList rRnaIntervalList = new IntervalList(builder.getHeader());
    rRnaIntervalList.add(rRnaInterval);
    final File rRnaIntervalsFile = BaseTest.createTempFile("tmp.rRna.", ".interval_list");
    rRnaIntervalList.write(rRnaIntervalsFile);
    // Generate the metrics.
    final File metricsFile = BaseTest.createTempFile("tmp.", ".rna_metrics");
    final String[] args = new String[] { "--input", samFile.getAbsolutePath(), "--output", metricsFile.getAbsolutePath(), "--REF_FLAT", getRefFlatFile(sequence).getAbsolutePath(), "--RIBOSOMAL_INTERVALS", rRnaIntervalsFile.getAbsolutePath(), "--STRAND_SPECIFICITY", "SECOND_READ_TRANSCRIPTION_STRAND", "--IGNORE_SEQUENCE", ignoredSequence, "--LEVEL", "SAMPLE", "--LEVEL", "LIBRARY" };
    runCommandLine(args);
    final MetricsFile<RnaSeqMetrics, Comparable<?>> output = new MetricsFile<>();
    output.read(new FileReader(metricsFile));
    for (final RnaSeqMetrics metrics : output.getMetrics()) {
        if (metrics.LIBRARY == null) {
            Assert.assertEquals(metrics.PF_ALIGNED_BASES, 396);
            Assert.assertEquals(metrics.PF_BASES, 432);
            Assert.assertEquals(metrics.RIBOSOMAL_BASES.longValue(), 108L);
            Assert.assertEquals(metrics.CODING_BASES, 136);
            Assert.assertEquals(metrics.UTR_BASES, 51);
            Assert.assertEquals(metrics.INTRONIC_BASES, 50);
            Assert.assertEquals(metrics.INTERGENIC_BASES, 51);
            Assert.assertEquals(metrics.CORRECT_STRAND_READS, 3);
            Assert.assertEquals(metrics.INCORRECT_STRAND_READS, 4);
            Assert.assertEquals(metrics.IGNORED_READS, 1);
        } else if (metrics.LIBRARY.equals("foo")) {
            Assert.assertEquals(metrics.PF_ALIGNED_BASES, 216);
            Assert.assertEquals(metrics.PF_BASES, 216);
            Assert.assertEquals(metrics.RIBOSOMAL_BASES.longValue(), 36L);
            Assert.assertEquals(metrics.CODING_BASES, 89);
            Assert.assertEquals(metrics.UTR_BASES, 51);
            Assert.assertEquals(metrics.INTRONIC_BASES, 25);
            Assert.assertEquals(metrics.INTERGENIC_BASES, 15);
            Assert.assertEquals(metrics.CORRECT_STRAND_READS, 3);
            Assert.assertEquals(metrics.INCORRECT_STRAND_READS, 2);
            Assert.assertEquals(metrics.IGNORED_READS, 0);
        } else if (metrics.LIBRARY.equals("bar")) {
            Assert.assertEquals(metrics.PF_ALIGNED_BASES, 180);
            Assert.assertEquals(metrics.PF_BASES, 216);
            Assert.assertEquals(metrics.RIBOSOMAL_BASES.longValue(), 72L);
            Assert.assertEquals(metrics.CODING_BASES, 47);
            Assert.assertEquals(metrics.UTR_BASES, 0);
            Assert.assertEquals(metrics.INTRONIC_BASES, 25);
            Assert.assertEquals(metrics.INTERGENIC_BASES, 36);
            Assert.assertEquals(metrics.CORRECT_STRAND_READS, 0);
            Assert.assertEquals(metrics.INCORRECT_STRAND_READS, 2);
            Assert.assertEquals(metrics.IGNORED_READS, 1);
        }
    }
}
Also used : MetricsFile(htsjdk.samtools.metrics.MetricsFile) SAMFileWriter(htsjdk.samtools.SAMFileWriter) SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord) SAMFileWriterFactory(htsjdk.samtools.SAMFileWriterFactory) SAMRecordSetBuilder(htsjdk.samtools.SAMRecordSetBuilder) IntervalList(htsjdk.samtools.util.IntervalList) SAMRecord(htsjdk.samtools.SAMRecord) FileReader(java.io.FileReader) MetricsFile(htsjdk.samtools.metrics.MetricsFile) File(java.io.File) Interval(htsjdk.samtools.util.Interval) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test) CommandLineProgramTest(org.broadinstitute.hellbender.CommandLineProgramTest)

Example 3 with Interval

use of htsjdk.samtools.util.Interval in project gatk by broadinstitute.

the class HetPulldownCalculatorUnitTest method inputGetTumorHetPulldown.

@DataProvider(name = "inputGetTumorHetPulldown")
public Object[][] inputGetTumorHetPulldown() {
    final Pulldown tumorHetPulldown = new Pulldown(normalHeader);
    tumorHetPulldown.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4));
    tumorHetPulldown.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6));
    tumorHetPulldown.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8));
    tumorHetPulldown.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9));
    tumorHetPulldown.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5));
    final IntervalList normalHetIntervals = new IntervalList(tumorHeader);
    normalHetIntervals.add(new Interval("1", 11522, 11522));
    normalHetIntervals.add(new Interval("1", 12098, 12098));
    normalHetIntervals.add(new Interval("1", 14630, 14630));
    normalHetIntervals.add(new Interval("2", 14689, 14689));
    normalHetIntervals.add(new Interval("2", 14982, 14982));
    return new Object[][] { { normalHetIntervals, tumorHetPulldown } };
}
Also used : IntervalList(htsjdk.samtools.util.IntervalList) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) AllelicCount(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Interval(htsjdk.samtools.util.Interval) DataProvider(org.testng.annotations.DataProvider)

Example 4 with Interval

use of htsjdk.samtools.util.Interval in project gatk by broadinstitute.

the class HetPulldownCalculatorUnitTest method inputGetPileupBaseCount.

@DataProvider(name = "inputGetPileupBaseCount")
public Object[][] inputGetPileupBaseCount() throws IOException {
    try (final SamReader bamReader = SamReaderFactory.makeDefault().open(NORMAL_BAM_FILE)) {
        final IntervalList intervals = new IntervalList(bamReader.getFileHeader());
        intervals.add(new Interval("1", 100, 100));
        intervals.add(new Interval("1", 11000, 11000));
        intervals.add(new Interval("1", 14000, 14000));
        intervals.add(new Interval("1", 14630, 14630));
        final SamLocusIterator locusIterator = new SamLocusIterator(bamReader, intervals);
        final Nucleotide.Counter baseCounts1 = makeBaseCounts(0, 0, 0, 0);
        final Nucleotide.Counter baseCounts2 = makeBaseCounts(0, 9, 0, 0);
        final Nucleotide.Counter baseCounts3 = makeBaseCounts(12, 0, 0, 0);
        final Nucleotide.Counter baseCounts4 = makeBaseCounts(0, 0, 8, 9);
        if (!locusIterator.hasNext()) {
            throw new SAMException("Can't get locus to start iteration. Check that " + NORMAL_BAM_FILE.toString() + " contains 1:0-16000.");
        }
        final SamLocusIterator.LocusInfo locus1 = locusIterator.next();
        final SamLocusIterator.LocusInfo locus2 = locusIterator.next();
        final SamLocusIterator.LocusInfo locus3 = locusIterator.next();
        final SamLocusIterator.LocusInfo locus4 = locusIterator.next();
        locusIterator.close();
        return new Object[][] { { locus1, baseCounts1 }, { locus2, baseCounts2 }, { locus3, baseCounts3 }, { locus4, baseCounts4 } };
    }
}
Also used : SamLocusIterator(htsjdk.samtools.util.SamLocusIterator) IntervalList(htsjdk.samtools.util.IntervalList) Nucleotide(org.broadinstitute.hellbender.utils.Nucleotide) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Interval(htsjdk.samtools.util.Interval) DataProvider(org.testng.annotations.DataProvider)

Example 5 with Interval

use of htsjdk.samtools.util.Interval in project gatk by broadinstitute.

the class Pulldown method getIntervals.

/** Returns a new instance of an IntervalList, constructed from the intervals of the internally held
     * AllelicCounts.  This IntervalList is modifiable and does not change with the state of the Pulldown.   */
public IntervalList getIntervals() {
    final IntervalList intervals = new IntervalList(header);
    intervals.addall(getCounts().stream().map(AllelicCount::getInterval).map(si -> new Interval(si.getContig(), si.getStart(), si.getEnd())).collect(Collectors.toList()));
    return intervals;
}
Also used : IntervalList(htsjdk.samtools.util.IntervalList) Interval(htsjdk.samtools.util.Interval)

Aggregations

Interval (htsjdk.samtools.util.Interval)24 IntervalList (htsjdk.samtools.util.IntervalList)23 File (java.io.File)10 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)8 SAMFileHeader (htsjdk.samtools.SAMFileHeader)6 DataProvider (org.testng.annotations.DataProvider)6 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)5 Test (org.testng.annotations.Test)5 QueryInterval (htsjdk.samtools.QueryInterval)4 AllelicCount (org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount)4 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)3 IndexedFastaSequenceFile (htsjdk.samtools.reference.IndexedFastaSequenceFile)3 SAMFileWriter (htsjdk.samtools.SAMFileWriter)2 SAMFileWriterFactory (htsjdk.samtools.SAMFileWriterFactory)2 SAMRecord (htsjdk.samtools.SAMRecord)2 SAMRecordSetBuilder (htsjdk.samtools.SAMRecordSetBuilder)2 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)2 MetricsFile (htsjdk.samtools.metrics.MetricsFile)2 RuntimeIOException (htsjdk.samtools.util.RuntimeIOException)2 SamLocusIterator (htsjdk.samtools.util.SamLocusIterator)2