@Override protected void finish() { collector.finish(); final MetricsFile<RnaSeqMetrics, Integer> file = getMetricsFile(); collector.addAllLevelsToFile(file); file.write(OUTPUT); boolean atLeastOneHistogram = false; for (final Histogram<Integer> histo : file.getAllHistograms()) { atLeastOneHistogram = atLeastOneHistogram || !histo.isEmpty(); } // Generate the coverage by position plot if (CHART_OUTPUT != null && atLeastOneHistogram) { final int rResult = RExecutor.executeFromClasspath("picard/analysis/rnaSeqCoverage.R", OUTPUT.getAbsolutePath(), CHART_OUTPUT.getAbsolutePath(), INPUT.getName(), this.plotSubtitle); if (rResult != 0) { throw new PicardException("Problem invoking R to generate plot."); } } }
@Test(dataProvider = "TheoreticalSensitivityConstantDepthDataProvider") public void testSensitivityAtConstantDepth(final double expected, final File metricsFile, final double alleleFraction, final int depth, final int sampleSize, final double tolerance) throws Exception { // This tests Theoretical Sensitivity assuming a uniform depth with a distribution of base quality scores. // Because this only tests sensitivity at a constant depth, we use this for testing the code at high depths. final MetricsFile<?, Integer> metrics = new MetricsFile<>(); try (final FileReader metricsFileReader = new FileReader(metricsFile)) { metrics.read(metricsFileReader); } final List<Histogram<Integer>> histograms = metrics.getAllHistograms(); final Histogram<Integer> qualityHistogram = histograms.get(1); // We ensure that even using different random seeds we converge to roughly the same value. for (int i = 0; i < 3; i++) { double result = TheoreticalSensitivity.sensitivityAtConstantDepth(depth, qualityHistogram, 3, sampleSize, alleleFraction, i); Assert.assertEquals(result, expected, tolerance); } }
@Override protected void finish() { collector.finish(); final MetricsFile<RnaSeqMetrics, Integer> file = getMetricsFile(); collector.addAllLevelsToFile(file); file.write(OUTPUT); boolean atLeastOneHistogram = false; for (final Histogram<Integer> histo : file.getAllHistograms()) { atLeastOneHistogram = atLeastOneHistogram || !histo.isEmpty(); } // Generate the coverage by position plot if (CHART_OUTPUT != null && atLeastOneHistogram) { final int rResult = RExecutor.executeFromClasspath("picard/analysis/rnaSeqCoverage.R", OUTPUT.getAbsolutePath(), CHART_OUTPUT.getAbsolutePath(), INPUT.getName(), this.plotSubtitle); if (rResult != 0) { throw new PicardException("Problem invoking R to generate plot."); } } }
@Test(dataProvider = "TheoreticalSensitivityDataProvider") public void testSensitivity(final double expected, final File metricsFile, final double alleleFraction, final int sampleSize) throws Exception { // This tests Theoretical Sensitivity using distributions on both base quality scores // and the depth histogram. // We use a pretty forgiving tolerance here because for these tests // we are not using large enough sample sizes to converge. final double tolerance = 0.02; final MetricsFile<?, Integer> metrics = new MetricsFile<>(); try (final FileReader metricsFileReader = new FileReader(metricsFile)) { metrics.read(metricsFileReader); } final List<Histogram<Integer>> histograms = metrics.getAllHistograms(); final Histogram<Integer> depthHistogram = histograms.get(0); final Histogram<Integer> qualityHistogram = histograms.get(1); final double result = TheoreticalSensitivity.theoreticalSensitivity(depthHistogram, qualityHistogram, sampleSize, 3, alleleFraction); Assert.assertEquals(result, expected, tolerance); }
/** * Histogram Width was being incorrectly set causing histograms to be trimmed an removed inappropriately. * See https://github.com/broadinstitute/picard/issues/253 * Test to be sure that the right number of histograms are being output. */ @Test public void testHistogramWidthIsSetProperly() throws IOException { final File input = new File(TEST_DATA_DIR, "insert_size_metrics_test.sam"); final File outfile = File.createTempFile("test", ".insert_size_metrics"); final File pdf = File.createTempFile("test", ".pdf"); outfile.deleteOnExit(); pdf.deleteOnExit(); final String[] args = new String[]{ "INPUT=" + input.getAbsolutePath(), "OUTPUT=" + outfile.getAbsolutePath(), "Histogram_FILE=" + pdf.getAbsolutePath(), "LEVEL=null", "LEVEL=READ_GROUP" }; Assert.assertEquals(runPicardCommandLine(args), 0); final MetricsFile<InsertSizeMetrics, Comparable<?>> output = new MetricsFile<InsertSizeMetrics, Comparable<?>>(); output.read(new FileReader(outfile)); Assert.assertEquals(output.getAllHistograms().size(), 5); }
@Override public void test() throws IOException { final MetricsFile<DuplicationMetrics, Double> metricsOutput = testMetrics(); // Check contents of set size bin against expected values if (!expectedSetSizeMap.isEmpty()) { boolean checked = false; for (final Histogram<Double> histo : metricsOutput.getAllHistograms()) { final String label = histo.getValueLabel(); for (final Double bin : histo.keySet()) { final String binStr = String.valueOf(bin); final List<String> labelBinStr = Arrays.asList(label, binStr); if (expectedSetSizeMap.containsKey(labelBinStr)) { checked = true; Histogram.Bin<Double> binValue = histo.get(bin); final double actual = binValue.getValue(); final double expected = expectedSetSizeMap.get(labelBinStr); Assert.assertEquals(actual, expected); } } } if (!checked) { Assert.fail("Could not not find matching entry for expectedSetSizeMap in metrics."); } } }
@Test(dataProvider = "equivalanceHetVsArbitrary") public void testHetVsArbitrary(final File metricsFile, final double tolerance, final int sampleSize) throws Exception { // This test compares Theoretical Sensitivity for arbitrary allele fractions with the theoretical het sensitivity // model. Since allele fraction of 0.5 is equivalent to a het, these should provide the same answer. final MetricsFile<?, Integer> metrics = new MetricsFile<>(); try (final FileReader metricsFileReader = new FileReader(metricsFile)) { metrics.read(metricsFileReader); } final List<Histogram<Integer>> histograms = metrics.getAllHistograms(); final Histogram<Integer> depthHistogram = histograms.get(0); final Histogram<Integer> qualityHistogram = histograms.get(1); final double[] qualityDistribution = TheoreticalSensitivity.normalizeHistogram(qualityHistogram); final double[] depthDistribution = TheoreticalSensitivity.normalizeHistogram(depthHistogram); final double resultFromTS = TheoreticalSensitivity.theoreticalSensitivity(depthHistogram, qualityHistogram, sampleSize, 3, 0.5); final double resultFromTHS = TheoreticalSensitivity.hetSNPSensitivity(depthDistribution, qualityDistribution, sampleSize, 3); Assert.assertEquals(resultFromTS, resultFromTHS, tolerance); }
@Test(dataProvider = "sumOfGaussiansDataProvider") public void testDrawSumOfQScores(final File metricsFile, final int altDepth, final double tolerance) throws Exception { final MetricsFile<TheoreticalSensitivityMetrics, Integer> metrics = new MetricsFile<>(); try (final FileReader metricsFileReader = new FileReader(metricsFile)) { metrics.read(metricsFileReader); } final List<Histogram<Integer>> histograms = metrics.getAllHistograms(); final Histogram<Integer> qualityHistogram = histograms.get(1); final TheoreticalSensitivity.RouletteWheel qualityRW = new TheoreticalSensitivity.RouletteWheel(TheoreticalSensitivity.trimDistribution(TheoreticalSensitivity.normalizeHistogram(qualityHistogram))); final Random randomNumberGenerator = new Random(51); // Calculate mean and deviation of quality score distribution to enable Gaussian sampling below final double averageQuality = qualityHistogram.getMean(); final double standardDeviationQuality = qualityHistogram.getStandardDeviation(); for (int k = 0; k < 1; k++) { int sumOfQualitiesFull = IntStream.range(0, altDepth).map(n -> qualityRW.draw()).sum(); int sumOfQualities = TheoreticalSensitivity.drawSumOfQScores(altDepth, averageQuality, standardDeviationQuality, randomNumberGenerator.nextGaussian()); Assert.assertEquals(sumOfQualitiesFull, sumOfQualities, sumOfQualitiesFull * tolerance); } }
for (final Histogram<Integer> histogram : output.getAllHistograms()) { if (histogram.getValueLabel().equals("count_WHOLE_GENOME")) { Assert.assertEquals(histogram.get(0).getValue(), 1080d);
final Histogram<Integer> coverageHistogram = output.getAllHistograms().get(0); Assert.assertEquals(coverageHistogram.get(0).getValue(), 10.0); Assert.assertEquals(coverageHistogram.get(1).getValue(), 10.0);
@Test() public void testRawBqDistributionWithSoftClips() throws IOException { final String input = TEST_DATA_DIR + "chrMReadsWithClips.sam"; final File outFile = File.createTempFile("test", ".TargetedMetrics_Coverage"); outFile.deleteOnExit(); final String[] args = new String[] { "TARGET_INTERVALS=" + singleIntervals, "INPUT=" + input, "OUTPUT=" + outFile.getAbsolutePath(), "REFERENCE_SEQUENCE=" + referenceFile, "LEVEL=ALL_READS", "AMPLICON_INTERVALS=" + singleIntervals, "SAMPLE_SIZE=" + 0 }; Assert.assertEquals(runPicardCommandLine(args), 0); final MetricsFile<TargetedPcrMetrics, Comparable<Integer>> output = new MetricsFile<>(); output.read(new FileReader(outFile)); Assert.assertEquals(output.getMetrics().size(), 1); for (final TargetedPcrMetrics metrics : output.getMetrics()) { Assert.assertEquals(metrics.TOTAL_READS, 2); } Assert.assertEquals(output.getNumHistograms(), 2); final Histogram<Comparable<Integer>> histogram = output.getAllHistograms().get(1); Assert.assertTrue(TestNGUtil.compareDoubleWithAccuracy(histogram.getSumOfValues(), 62,0.01)); Assert.assertTrue(TestNGUtil.compareDoubleWithAccuracy(histogram.get(32).getValue(), 52D, 0.01)); Assert.assertTrue(TestNGUtil.compareDoubleWithAccuracy(histogram.get(33).getValue(), 10D, 0.01)); }
for (final Histogram<Integer> histogram : output.getAllHistograms()) { if (histogram.getValueLabel().equals("count_WHOLE_GENOME")) { Assert.assertEquals(histogram.get(0).getValue(), 274d);
final CollectWgsMetrics.WgsMetrics metrics = output.getMetrics().get(0); final Histogram<Integer> highQualityDepthHistogram = output.getAllHistograms().get(0); final Histogram<Integer> baseQHistogram = output.getAllHistograms().get(1);