From 3308a110bdacb79043ba2e07ea472f2a1c8932dd Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 26 Apr 2017 12:28:51 -0700 Subject: [PATCH] Adding MAX_TARGET_COVERAGE to targeted metrics Also adds the metric to HsMetrics and TargetedPcrMetrics. --- src/main/java/picard/analysis/directed/HsMetrics.java | 3 +++ .../picard/analysis/directed/TargetMetricsCollector.java | 8 ++++++++ .../java/picard/analysis/directed/TargetedPcrMetrics.java | 3 +++ .../java/picard/analysis/directed/CollectHsMetricsTest.java | 12 +++++++----- 4 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/main/java/picard/analysis/directed/HsMetrics.java b/src/main/java/picard/analysis/directed/HsMetrics.java index 1e6abaa92..f0e9fe23b 100644 --- a/src/main/java/picard/analysis/directed/HsMetrics.java +++ b/src/main/java/picard/analysis/directed/HsMetrics.java @@ -121,6 +121,9 @@ /** The median coverage of a target region. */ public double MEDIAN_TARGET_COVERAGE; + /** The maximum coverage of reads that mapped to target regions of an experiment. */ + public long MAX_TARGET_COVERAGE; + /** The number of aligned, de-duped, on-bait bases out of the PF bases available. */ public double PCT_USABLE_BASES_ON_BAIT; diff --git a/src/main/java/picard/analysis/directed/TargetMetricsCollector.java b/src/main/java/picard/analysis/directed/TargetMetricsCollector.java index e11a6f234..650f4b287 100644 --- a/src/main/java/picard/analysis/directed/TargetMetricsCollector.java +++ b/src/main/java/picard/analysis/directed/TargetMetricsCollector.java @@ -637,6 +637,9 @@ private void calculateTargetCoverageMetrics() { // the number of bases we counted towards the depth histogram plus those that got thrown out by the coverage cap long totalCoverage = 0; + // the maximum depth at any target base + long maxDepth = 0; + // The "how many target bases at at-least X" calculations. // downstream code relies on this array being sorted in ascending order final int[] targetBasesDepth = {0, 1, 2, 10, 20, 30, 40, 50, 100}; @@ -657,6 +660,7 @@ private void calculateTargetCoverageMetrics() { for (final int depth : c.getDepths()) { totalCoverage += depth; highQualityCoverageHistogramArray[Math.min(depth, coverageCap)]++; + maxDepth = Math.max(maxDepth, depth); // Add to the "how many target bases at at-least X" calculations. for (int i = 0; i < targetBasesDepth.length; i++) { @@ -677,6 +681,7 @@ private void calculateTargetCoverageMetrics() { // we do this instead of highQualityDepthHistogram.getMean() because the histogram imposes a coverage cap metrics.MEAN_TARGET_COVERAGE = (double) totalCoverage / metrics.TARGET_TERRITORY; metrics.MEDIAN_TARGET_COVERAGE = highQualityDepthHistogram.getMedian(); + metrics.MAX_TARGET_COVERAGE = maxDepth; // compute the coverage value such that 80% of target bases have better coverage than it i.e. 20th percentile // this roughly measures how much we must sequence extra such that 80% of target bases have coverage at least as deep as the current mean coverage @@ -1008,6 +1013,9 @@ public String toString() { /** The median coverage of targets. */ public double MEDIAN_TARGET_COVERAGE; + /** The maximum coverage of reads that mapped to target regions of an experiment. */ + public long MAX_TARGET_COVERAGE; + /** The fraction of targets that did not reach coverage=1 over any base. */ public double ZERO_CVG_TARGETS_PCT; diff --git a/src/main/java/picard/analysis/directed/TargetedPcrMetrics.java b/src/main/java/picard/analysis/directed/TargetedPcrMetrics.java index 96b580a8c..fc2207b9d 100644 --- a/src/main/java/picard/analysis/directed/TargetedPcrMetrics.java +++ b/src/main/java/picard/analysis/directed/TargetedPcrMetrics.java @@ -93,6 +93,9 @@ /** The median coverage of reads that mapped to target regions of an experiment. */ public double MEDIAN_TARGET_COVERAGE; + /** The maximum coverage of reads that mapped to target regions of an experiment. */ + public long MAX_TARGET_COVERAGE; + /** The fold by which the amplicon region has been amplified above genomic background. */ public double FOLD_ENRICHMENT; diff --git a/src/test/java/picard/analysis/directed/CollectHsMetricsTest.java b/src/test/java/picard/analysis/directed/CollectHsMetricsTest.java index 287d59cfc..6c72ef768 100644 --- a/src/test/java/picard/analysis/directed/CollectHsMetricsTest.java +++ b/src/test/java/picard/analysis/directed/CollectHsMetricsTest.java @@ -28,15 +28,15 @@ public String getCommandLineProgramName() { return new Object[][] { // two reads, each has 100 bases. bases in one read are medium quality (20), in the other read poor quality (10). // test that we exclude half of the bases - {TEST_DIR + "/lowbaseq.sam", intervals, 1, 10, true, 2, 200, 0.5, 0.0, 0.50, 0.0, 1000}, + {TEST_DIR + "/lowbaseq.sam", intervals, 1, 10, true, 2, 200, 0.5, 0.0, 0.50, 0.0, 1, 1000}, // test that read 2 (with mapping quality 1) is filtered out with minimum mapping quality 2 - {TEST_DIR + "/lowmapq.sam", intervals, 2, 0, true, 2, 202, 0, 0.0, 0.505, 0.0, 1000}, + {TEST_DIR + "/lowmapq.sam", intervals, 2, 0, true, 2, 202, 0, 0.0, 0.505, 0.0, 1, 1000}, // test that we clip overlapping bases - {TEST_DIR + "/overlapping.sam", intervals, 0, 0, true, 2, 202, 0, 0.5, 0.505, 0, 1000}, + {TEST_DIR + "/overlapping.sam", intervals, 0, 0, true, 2, 202, 0, 0.5, 0.505, 0, 1, 1000}, // test that we do not clip overlapping bases - {TEST_DIR + "/overlapping.sam", intervals, 0, 0, false, 2, 202, 0, 0.0, 0.505, 0.505, 1000}, + {TEST_DIR + "/overlapping.sam", intervals, 0, 0, false, 2, 202, 0, 0.0, 0.505, 0.505, 2, 1000}, // A read 10 base pairs long. two intervals: one maps identically to the read, other does not overlap at all - {TEST_DIR + "/single-short-read.sam", twoSmallIntervals, 20, 20, true, 1, 10, 0.0, 0.0, 0.5, 0.0, 1000 } + {TEST_DIR + "/single-short-read.sam", twoSmallIntervals, 20, 20, true, 1, 10, 0.0, 0.0, 0.5, 0.0, 1, 1000 } }; } @@ -53,6 +53,7 @@ public void runCollectTargetedMetricsTest(final String input, final double pctExcOverlap, final double pctTargetBases1x, final double pctTargetBases2x, + final long maxTargetCoverage, final int sampleSize) throws IOException { final File outfile = File.createTempFile("CollectHsMetrics", ".hs_metrics", TEST_DIR); @@ -82,6 +83,7 @@ public void runCollectTargetedMetricsTest(final String input, Assert.assertEquals(metrics.PCT_EXC_OVERLAP, pctExcOverlap); Assert.assertEquals(metrics.PCT_TARGET_BASES_1X, pctTargetBases1x); Assert.assertEquals(metrics.PCT_TARGET_BASES_2X, pctTargetBases2x); + Assert.assertEquals(metrics.MAX_TARGET_COVERAGE, maxTargetCoverage); } }