Added handling for indels within CollectHsMetrics and CollectTargetedPcrMetrics. #1273

tfenne · 2019-01-23T22:27:31Z

Description

When the INCLUDE_INDELS parameter is set to true:

Deletions within reads are treated as covering the reference bases if the preceding base is high quality
Inserted bases are counted as "on-target" but don't contribute coverage to any target base(s)

There are two motivations for this change:

I want to be able to use CollectHsMetrics' per-base coverage output to assess how well an assay is covering target bases. If the sample being sequenced has one or more short deletions within the targets this results in a precipitous drop in target coverage. However, I'd argue that those bases in the reference/target are covered, and if covered sufficiently will allow for a deletion to be called in variant calling.
When the het sensitivity code was introduced a bug was introduced that would lead to reads being double-counted in some places if a read overlapped more than one target. E.g. I would expect that, other than for rounding issues, PCT_USABLE_BASES_ON_TARGET + PCT_EXC_* to add up to 1.0. However I would frequently see this adding up to substantially over 1.0.

I've done my best to re-organize the code to make it easier to follow and maintain in future.

Checklist

Content

Added or modified tests to cover changes and any new functionality
Edited the README / documentation (if applicable)
All tests passing on Travis

Review

Final thumbs-up from reviewer
Rebase, squash and reword as applicable

…PcrMetrics. When the INCLUDE_INDELS parameter is set to true: 1. Deletions within reads are treated as covering the reference bases if the preceding base is high quality 2. Inserted bases are counted as "on-target" but don't contribute coverage to any target base(s)

tfenne · 2019-01-23T22:28:24Z

@yfarjoun I'm not really done with this PR, I plan to add some additional tests, but wanted to get this up for you, and anyone else who is interested, to look at and see if you agree with direction. Thanks!

coveralls · 2019-01-23T23:05:09Z

Coverage increased (+0.01%) to 81.493% when pulling 283341e on tfenne:tf_hs_metrics_with_indels into 6ad57dc on broadinstitute:master.

…unPicardCommandLine("I=", input) vs. runPicardCommandLine("I=" + input).

yfarjoun

Looks good overall. thanks for the tests. I have a few comments but nothing major.

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

yfarjoun · 2019-01-28T14:28:06Z

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

- isOnTarget = true;
+ // Then go through the per-target/per-base hq and unfiltered coverage
+ // The cutoff of >= 2 is because even the unfilteredCoverage doesn't want Q1 or Q0 bases!
+ if (qual > 2 && incrementPerTargetCoverage && onTarget) {


Suggested change

if (qual > 2 && incrementPerTargetCoverage && onTarget) {

if (qual >= 2 && incrementPerTargetCoverage && onTarget) {

I'll fix the comment, the code is actually correct here - the het sensitivity code only wants bases > q2, not >= q2.

yfarjoun · 2019-01-28T14:44:36Z

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

+ if (qual > 2 && incrementPerTargetCoverage && onTarget) {
+ for (final Interval target : targets) {
+ if (refPos >= target.getStart() && refPos <= target.getEnd()) {
+ final int targetOffset = refPos - target.getStart();


refactor into overlapsTarget and use below in 619 as well?

if (refPos >= target.getStart() && refPos <= target.getEnd()) {

yfarjoun · 2019-01-28T16:22:28Z

src/test/java/picard/analysis/directed/CollectHsMetricsTest.java

+
+ // Now run CollectHsMetrics four times
+ final File out = Files.createTempFile("hsmetrics.", ".txt").toFile();
+ runPicardCommandLine(Arrays.asList("INCLUDE_INDELS=false", "SAMPLE_SIZE=0", "TI="+ts.getPath(), "BI="+bs.getPath(), "O="+out.getPath(), "I="+withDelBam.getAbsolutePath()));


this might be more readable with a data-provider...no?

I personally find the whole TestNG data provider stuff makes tests harder to read, which is why I avoid it.

yfarjoun · 2019-01-28T16:23:32Z

src/test/java/picard/analysis/directed/CollectHsMetricsTest.java

+ }
+
+ // Write things out to file
+ final File dir = IOUtil.createTempDir("hsmetrics.", ".test");


need to delete the directory yourself.

yfarjoun · 2019-01-28T16:28:05Z

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

@@ -548,63 +550,79 @@ public void acceptRecord(final SAMRecord record) {
 rec = record;
 }

- // Find the target overlaps


Since you are here, in line 331 above, can you replace the use of getNotPrimaryAlignmentFlag() with the non-deprecated isSecondaryAlignment ?

yfarjoun · 2019-01-28T19:07:54Z

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

@@ -548,63 +550,79 @@ public void acceptRecord(final SAMRecord record) {
 rec = record;
 }

- // Find the target overlaps


Also, do you think that the clipping for overlapping bases should be modified in any way if we take indels into account?

I think the clipping is probably fine. The clipping is based on reference positions IIRC, which should be ok.

yfarjoun · 2019-01-31T03:27:53Z

src/main/java/picard/analysis/directed/TargetMetricsCollector.java

+ final int len = cig.getLength();
+
+ for (int i=0; i<len; ++i) {
+ if (op.isAlignment() || (this.includeIndels && op.isIndel())) {


can this line be swapped with the previous line as it seems to encompass all of the for-loop and doesn't depend on i.

The whole if is encompassed, but not the statements on lines 612 and 613 that update readOffset and refOffset. I did at one point try inverting and adding an else that updated the offsets using the length of the element and saw no performance increase, so went back to how it is now.

yfarjoun · 2019-01-31T03:28:51Z

src/test/java/picard/analysis/directed/CollectHsMetricsTest.java

+
+ /** Writes the contents of a SAMRecordSetBuilder out to a file. */
+ File writeBam(final SAMRecordSetBuilder builder, final File f) {
+ final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(builder.getHeader(), false, f);


use try-with-resources?

Co-Authored-By: tfenne <[email protected]>

tfenne · 2019-01-31T20:14:32Z

Thanks for the review @yfarjoun - I've implemented most of your suggestions and commented where I haven't.

yfarjoun

Thanks, I suspect that this can make a big difference for small panels....

tfenne requested a review from yfarjoun January 23, 2019 22:27

tfenne added 3 commits January 23, 2019 16:06

Fix performance regression.

57bc875

Added test for indel handling in hs metrics.

7fd7d8f

Fixed failing test. FWIW the barclayTest doesn't like it if you use r…

f8d0b66

…unPicardCommandLine("I=", input) vs. runPicardCommandLine("I=" + input).

yfarjoun suggested changes Jan 31, 2019

View reviewed changes

Yossi Farjoun and others added 2 commits January 31, 2019 13:12

Apply suggestions from code review

767de36

Co-Authored-By: tfenne <[email protected]>

Changes to address yfarjoun's review comments.

283341e

yfarjoun approved these changes Feb 1, 2019

View reviewed changes

tfenne merged commit 42c8ee0 into broadinstitute:master Feb 1, 2019

tfenne deleted the tf_hs_metrics_with_indels branch February 1, 2019 21:11

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added handling for indels within CollectHsMetrics and CollectTargetedPcrMetrics. #1273

Added handling for indels within CollectHsMetrics and CollectTargetedPcrMetrics. #1273

tfenne commented Jan 23, 2019

tfenne commented Jan 23, 2019

coveralls commented Jan 23, 2019 •

edited

Loading

yfarjoun left a comment

yfarjoun Jan 28, 2019

tfenne Jan 31, 2019

yfarjoun Jan 28, 2019

yfarjoun Jan 28, 2019

tfenne Jan 31, 2019

yfarjoun Jan 28, 2019

yfarjoun Jan 28, 2019

yfarjoun Jan 28, 2019

tfenne Jan 31, 2019

yfarjoun Jan 31, 2019

tfenne Jan 31, 2019

yfarjoun Jan 31, 2019

tfenne commented Jan 31, 2019

yfarjoun left a comment

	if (qual > 2 && incrementPerTargetCoverage && onTarget) {
	if (qual >= 2 && incrementPerTargetCoverage && onTarget) {

Added handling for indels within CollectHsMetrics and CollectTargetedPcrMetrics. #1273

Added handling for indels within CollectHsMetrics and CollectTargetedPcrMetrics. #1273

Conversation

tfenne commented Jan 23, 2019

Description

Checklist

Content

Review

tfenne commented Jan 23, 2019

coveralls commented Jan 23, 2019 • edited Loading

yfarjoun left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tfenne commented Jan 31, 2019

yfarjoun left a comment

Choose a reason for hiding this comment

coveralls commented Jan 23, 2019 •

edited

Loading