From 2ed4968ed07cf30fad3af622867ee50a62b15d73 Mon Sep 17 00:00:00 2001 From: Ron Levine Date: Wed, 1 Feb 2017 17:45:34 -0500 Subject: [PATCH] Must use a RIBOSOMAL_INTERVALS file if RRNA_FRAGMENT_PERCENTAGE = 0.0 --- .../java/picard/analysis/CollectRnaSeqMetrics.java | 11 +++- .../analysis/directed/RnaSeqMetricsCollector.java | 8 ++- .../picard/analysis/CollectRnaSeqMetricsTest.java | 62 ++++++++++++++++++++++ 3 files changed, 78 insertions(+), 3 deletions(-) diff --git a/src/main/java/picard/analysis/CollectRnaSeqMetrics.java b/src/main/java/picard/analysis/CollectRnaSeqMetrics.java index dfeb67111..aed80f0bf 100644 --- a/src/main/java/picard/analysis/CollectRnaSeqMetrics.java +++ b/src/main/java/picard/analysis/CollectRnaSeqMetrics.java @@ -135,6 +135,15 @@ public static void main(final String[] argv) { } @Override + protected String[] customCommandLineValidation() { + // No ribosomal intervals file and rRNA fragment percentage = 0 + if ( RIBOSOMAL_INTERVALS == null && RRNA_FRAGMENT_PERCENTAGE == 0 ) { + throw new PicardException("Must use a RIBOSOMAL_INTERVALS file if RRNA_FRAGMENT_PERCENTAGE = 0.0"); + } + return super.customCommandLineValidation(); + } + + @Override protected void setup(final SAMFileHeader header, final File samFile) { if (CHART_OUTPUT != null) IOUtil.assertFileIsWritable(CHART_OUTPUT); @@ -143,7 +152,7 @@ protected void setup(final SAMFileHeader header, final File samFile) { LOG.info("Loaded " + geneOverlapDetector.getAll().size() + " genes."); final Long ribosomalBasesInitialValue = RIBOSOMAL_INTERVALS != null ? 0L : null; - final OverlapDetector ribosomalSequenceOverlapDetector = RnaSeqMetricsCollector.makeOverlapDetector(samFile, header, RIBOSOMAL_INTERVALS); + final OverlapDetector ribosomalSequenceOverlapDetector = RnaSeqMetricsCollector.makeOverlapDetector(samFile, header, RIBOSOMAL_INTERVALS, LOG); final HashSet ignoredSequenceIndices = RnaSeqMetricsCollector.makeIgnoredSequenceIndicesSet(header, IGNORE_SEQUENCE); diff --git a/src/main/java/picard/analysis/directed/RnaSeqMetricsCollector.java b/src/main/java/picard/analysis/directed/RnaSeqMetricsCollector.java index 8bfb3300a..36083c5c5 100644 --- a/src/main/java/picard/analysis/directed/RnaSeqMetricsCollector.java +++ b/src/main/java/picard/analysis/directed/RnaSeqMetricsCollector.java @@ -6,6 +6,7 @@ import htsjdk.samtools.util.Histogram; import htsjdk.samtools.util.Interval; import htsjdk.samtools.util.IntervalList; +import htsjdk.samtools.util.Log; import htsjdk.samtools.util.OverlapDetector; import htsjdk.samtools.util.SequenceUtil; import picard.PicardException; @@ -60,12 +61,15 @@ public RnaSeqMetricsCollector(final Set accumulationLev return new PerUnitRnaSeqMetricsCollector(sample, library, readGroup, ribosomalInitialValue); } - public static OverlapDetector makeOverlapDetector(final File samFile, final SAMFileHeader header, final File ribosomalIntervalsFile) { + public static OverlapDetector makeOverlapDetector(final File samFile, final SAMFileHeader header, final File ribosomalIntervalsFile, final Log log) { - OverlapDetector ribosomalSequenceOverlapDetector = new OverlapDetector(0, 0); + final OverlapDetector ribosomalSequenceOverlapDetector = new OverlapDetector(0, 0); if (ribosomalIntervalsFile != null) { final IntervalList ribosomalIntervals = IntervalList.fromFile(ribosomalIntervalsFile); + if (ribosomalIntervals.size() == 0) { + log.warn("The RIBOSOMAL_INTERVALS file, " + ribosomalIntervalsFile.getAbsolutePath() + " does not contain intervals"); + } try { SequenceUtil.assertSequenceDictionariesEqual(header.getSequenceDictionary(), ribosomalIntervals.getHeader().getSequenceDictionary()); } catch (SequenceUtil.SequenceListsDifferException e) { diff --git a/src/test/java/picard/analysis/CollectRnaSeqMetricsTest.java b/src/test/java/picard/analysis/CollectRnaSeqMetricsTest.java index fd25a8503..84da166c8 100644 --- a/src/test/java/picard/analysis/CollectRnaSeqMetricsTest.java +++ b/src/test/java/picard/analysis/CollectRnaSeqMetricsTest.java @@ -29,12 +29,15 @@ import htsjdk.samtools.util.IntervalList; import htsjdk.samtools.util.StringUtil; import org.testng.Assert; +import org.testng.annotations.DataProvider; import org.testng.annotations.Test; +import picard.PicardException; import picard.cmdline.CommandLineProgramTest; import picard.annotation.RefFlatReader.RefFlatColumns; import java.io.File; import java.io.FileReader; +import java.io.IOException; import java.io.PrintStream; public class CollectRnaSeqMetricsTest extends CommandLineProgramTest { @@ -112,6 +115,65 @@ public void testBasic() throws Exception { Assert.assertEquals(metrics.PCT_R2_TRANSCRIPT_STRAND_READS, 0.666667); } + @DataProvider(name = "rRnaIntervalsFiles") + public static Object[][] rRnaIntervalsFiles() throws IOException { + return new Object[][] { + {null}, + {File.createTempFile("tmp.rRna.", ".interval_list")} + }; + } + + @Test(dataProvider = "rRnaIntervalsFiles") + public void testNoIntevalsNoFragPercentage(final File rRnaIntervalsFile) throws Exception { + final SAMRecordSetBuilder builder = new SAMRecordSetBuilder(true, SAMFileHeader.SortOrder.coordinate); + + // Add a header but no intervals + if ( rRnaIntervalsFile != null ) { + final IntervalList rRnaIntervalList = new IntervalList(builder.getHeader()); + rRnaIntervalList.write(rRnaIntervalsFile); + rRnaIntervalsFile.deleteOnExit(); + } + + // Create some alignments + final String sequence = "chr1"; + final String ignoredSequence = "chrM"; + + // Set seed so that strandedness is consistent among runs. + builder.setRandomSeed(0); + final int sequenceIndex = builder.getHeader().getSequenceIndex(sequence); + builder.addPair("pair1", sequenceIndex, 45, 475); + builder.addFrag("ignoredFrag", builder.getHeader().getSequenceIndex(ignoredSequence), 1, false); + + final File samFile = File.createTempFile("tmp.collectRnaSeqMetrics.", ".sam"); + samFile.deleteOnExit(); + + final SAMFileWriter samWriter = new SAMFileWriterFactory().makeSAMWriter(builder.getHeader(), false, samFile); + for (final SAMRecord rec : builder.getRecords()) samWriter.addAlignment(rec); + samWriter.close(); + + // Generate the metrics. + final File metricsFile = File.createTempFile("tmp.", ".rna_metrics"); + metricsFile.deleteOnExit(); + + final String rRnaIntervalsPath = rRnaIntervalsFile != null ? rRnaIntervalsFile.getAbsolutePath() : null; + final String[] args = new String[]{ + "INPUT=" + samFile.getAbsolutePath(), + "OUTPUT=" + metricsFile.getAbsolutePath(), + "REF_FLAT=" + getRefFlatFile(sequence).getAbsolutePath(), + "RIBOSOMAL_INTERVALS=" + rRnaIntervalsPath, + "RRNA_FRAGMENT_PERCENTAGE=" + 0.0, + "STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND", + "IGNORE_SEQUENCE=" + ignoredSequence + }; + try { + Assert.assertEquals(runPicardCommandLine(args), 0); + } catch(final Exception e) { + if (rRnaIntervalsFile != null) { + Assert.assertEquals(e.getCause().getClass(), PicardException.class); + } + } + } + @Test public void testMultiLevel() throws Exception { final String sequence = "chr1";