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