PathSeq Illumina adapter trimming and simple repeat masking #3354

mwalker174 · 2017-07-24T23:01:52Z

Adds additional filtering steps to the PathSeq filter to 1) trim adapter sequences and 2) mimic a simple filter used in RepeatMasker that masks windows with excessive A/T or G/C content.

codecov-io · 2017-07-24T23:41:43Z

Codecov Report

Merging #3354 into master will increase coverage by 0.029%.
The diff coverage is 94.667%.

@@               Coverage Diff               @@
##              master     #3354       +/-   ##
===============================================
+ Coverage     80.472%   80.501%   +0.029%     
- Complexity     17507     17535       +28     
===============================================
  Files           1173      1175        +2     
  Lines          63376     63451       +75     
  Branches        9878      9894       +16     
===============================================
+ Hits           51000     51079       +79     
+ Misses          8428      8423        -5     
- Partials        3948      3949        +1

Impacted Files	Coverage Δ	Complexity Δ
...bender/tools/spark/pathseq/PathSeqFilterSpark.java	`70.968% <ø> (ø)`	`7 <0> (ø)`	⬇️
...itute/hellbender/tools/spark/pathseq/PSFilter.java	`92.617% <100%> (+0.531%)`	`33 <1> (+1)`	⬆️
...ools/spark/pathseq/PSFilterArgumentCollection.java	`80% <100%> (+1.429%)`	`2 <0> (ø)`	⬇️
...ellbender/transformers/AdapterTrimTransformer.java	`92.857% <92.857%> (ø)`	`12 <12> (?)`
...nder/transformers/SimpleRepeatMaskTransformer.java	`94.286% <94.286%> (ø)`	`11 <11> (?)`
...nstitute/hellbender/utils/clipping/ClippingOp.java	`84.365% <0%> (+1.629%)`	`91% <0%> (+2%)`	⬆️
...stitute/hellbender/utils/clipping/ReadClipper.java	`71.038% <0%> (+1.639%)`	`75% <0%> (+2%)`	⬆️

TedBrookings

I don't see any major problems so I'm pre-approving. The biggest issue is the comment I made in IlluminaAdapterTrimTransformer.java, which could be a logic error or could be down to lack of domain knowledge on my part.

TedBrookings · 2017-07-27T19:31:01Z

src/main/java/org/broadinstitute/hellbender/transformers/IlluminaAdapterTrimTransformer.java

+                        if (++numMismatches > maxMismatches) break;
+                    }
+                }
+                if (numMismatches <= maxMismatches) {


I am not familiar with adaptor sequence matching, but this seems like a potential problem. Someone would naively think that increasing maxMismatches will allow more poor matches through without affecting the already passing reads, but in principle this could cause a poorer match to be selected if a "good enough" match is found before the best one. Maybe better to keep a bestNumMismatches variable and a bestAlignmentStart, then update read on bestAlignmentStart if bestNumMismatches < maxMismatches

I think what you are saying is that if some adapter ends with the sequence of another adapter, or if a particular adapter ends with itself e.g.:

AGCAT

GAGCAT

Both end in AGCAT but if we had the read ending in GAGCAT and (2) happens to be tested before (1), it would only trim 4 bases instead of 5. Also if you had an adapter AAAAA and minimum clip length of 4, it would only clip 4 bases.

This is a very good point I didn't consider. I don't believe this is the case for any standard adapter sequences (Picard's MarkIlluminaAdapters also doesn't seem to account for it), but since this class exposes the use of custom adapters we should check for it.

In light of this, I'm also renaming the class to AdapterTrimTransformer

TedBrookings · 2017-07-27T20:09:06Z

src/main/java/org/broadinstitute/hellbender/transformers/SimpleRepeatMaskTransformer.java

+    final int windowSize, threshAT, threshGC;
+
+    public SimpleRepeatMaskTransformer(final int threshAT, final int threshGC, final int windowSize) {
+        this.threshAT = threshAT;


Maybe good to put in an assert that threshAT and threshGC are less than windowSize (if not, it indicates a typo or logic error).

TedBrookings · 2017-07-27T20:10:10Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PSFilter.java

-import org.broadinstitute.hellbender.transformers.BaseQualityClipReadTransformer;
-import org.broadinstitute.hellbender.transformers.BaseQualityReadTransformer;
-import org.broadinstitute.hellbender.transformers.DUSTReadTransformer;
+import org.broadinstitute.hellbender.transformers.*;


I don't feel strongly about this, but maybe better to import the classes individually so that dependencies can more easily be managed.

Intellij has this option of "optimizing the imports", under Menu section "Code". Try that?

This is what IntelliJ does

TedBrookings · 2017-07-27T20:12:53Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PSFilter.java

+
+            //Apply simple repeat masking
+            //See "Low-complexity DNA and simple repeats" at http://www.repeatmasker.org/webrepeatmaskerhelp.html
+            reads = reads.map(new ReadTransformerSparkifier(new SimpleRepeatMaskTransformer(29, 29, 30)));


If these numbers may need to be tweaked, it would be good to add them to filterArgs. If they're final numbers, it may still be good to name them as constants (e.g. shortWindowLen, maxShortWindowGCCounts or similar, so that it's easier to figure out what they do).

I agree with @TedBrookings . Generally we prefer to have no magic numbers baked in the code. And we do a favor for our future selfs not to have a magic number hiding here and there.

Added these as constants to the PSFilter class.

SHuang-Broad

A improvement comment, and a question on CIGAR in IlluminaAdapterTrimTransformer that could potentially be a more serious issue.

SHuang-Broad · 2017-07-28T16:13:13Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PSFilter.java

-import org.broadinstitute.hellbender.transformers.BaseQualityClipReadTransformer;
-import org.broadinstitute.hellbender.transformers.BaseQualityReadTransformer;
-import org.broadinstitute.hellbender.transformers.DUSTReadTransformer;
+import org.broadinstitute.hellbender.transformers.*;


Intellij has this option of "optimizing the imports", under Menu section "Code". Try that?

SHuang-Broad · 2017-07-28T16:19:16Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PSFilter.java

+
+            //Apply simple repeat masking
+            //See "Low-complexity DNA and simple repeats" at http://www.repeatmasker.org/webrepeatmaskerhelp.html
+            reads = reads.map(new ReadTransformerSparkifier(new SimpleRepeatMaskTransformer(29, 29, 30)));


I agree with @TedBrookings . Generally we prefer to have no magic numbers baked in the code. And we do a favor for our future selfs not to have a magic number hiding here and there.

SHuang-Broad · 2017-07-28T16:22:41Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PathSeqFilterSpark.java

@@ -47,8 +49,9 @@
 *
 */
 @CommandLineProgramProperties(summary = "First, low-quality and repetitive sequences reads are filtered: \n\t(1) Remove secondary and " +
-        "supplementary reads; \n\t(2) mask repetitive sequences using the sDUST algorithm; \n\t(3) hard clip according to masked " +
-        "and low-quality bases; \n\t(4) remove reads below minimum length; \n\t(5) mask low-quality bases; \n\t(6) filter reads with" +
+        "supplementary reads; \n\t(2) trim adapter sequences; \n\t(3) mask sequences with excessive A/T or G/C content; " +


a technically-more correct spelling is "secondary and supplementary reads" -> "secondary and supplementary alignments", but may not make more sense in your application.

Fair point. Done

SHuang-Broad · 2017-07-28T16:24:53Z

src/main/java/org/broadinstitute/hellbender/tools/spark/pathseq/PathSeqFilterSpark.java

@@ -47,8 +49,9 @@
 *
 */
 @CommandLineProgramProperties(summary = "First, low-quality and repetitive sequences reads are filtered: \n\t(1) Remove secondary and " +
-        "supplementary reads; \n\t(2) mask repetitive sequences using the sDUST algorithm; \n\t(3) hard clip according to masked " +
-        "and low-quality bases; \n\t(4) remove reads below minimum length; \n\t(5) mask low-quality bases; \n\t(6) filter reads with" +
+        "supplementary reads; \n\t(2) trim adapter sequences; \n\t(3) mask sequences with excessive A/T or G/C content; " +


IMHO, this summary is a bit too detailed. It exposes implementation detail.

I think I'm going to leave it as is for now until we get official documentation up, otherwise it's impossible to tell what the parameters do.

SHuang-Broad · 2017-07-28T16:26:46Z

src/main/java/org/broadinstitute/hellbender/transformers/IlluminaAdapterTrimTransformer.java

+ * Trims adapter sequences from read ends. If the adapter sequence ends before the end of the read, trims the remaining
+ * bases as well.
+ */
+public class IlluminaAdapterTrimTransformer implements ReadTransformer {


Do you intend this class to be final?

Yes, will also make SimpleRepeatMaskTransformer final.

SHuang-Broad · 2017-07-28T16:30:48Z

src/main/java/org/broadinstitute/hellbender/transformers/IlluminaAdapterTrimTransformer.java

+ * bases as well.
+ */
+public class IlluminaAdapterTrimTransformer implements ReadTransformer {
+    public static final long serialVersionUID = 1L;


why public?

Whoops - done.

SHuang-Broad · 2017-07-28T16:34:55Z

src/main/java/org/broadinstitute/hellbender/transformers/SimpleRepeatMaskTransformer.java

+    public static final long serialVersionUID = 1L;
+    private static final byte MASKED_BASE_QUAL = 2;
+
+    final int windowSize, threshAT, threshGC;


Maybe make them private? And, I am personally OK with declaring several primitives on the same line, but it most likely goes against our coding style. So your call.

Made private and put each variable on a separate line

SHuang-Broad · 2017-07-28T16:35:13Z

src/main/java/org/broadinstitute/hellbender/transformers/SimpleRepeatMaskTransformer.java

+ * Masks read bases with a supra-threshold number of A/T's or G/C's within a given window size.
+ */
+public class SimpleRepeatMaskTransformer implements ReadTransformer {
+    public static final long serialVersionUID = 1L;


why public?

SHuang-Broad · 2017-07-28T17:19:07Z

src/main/java/org/broadinstitute/hellbender/transformers/SimpleRepeatMaskTransformer.java

+    }
+
+    @Override
+    public GATKRead apply(final GATKRead read) {


A suggestion for how to implement it.

I'll start with what my understanding is:
This function, rolls along the provided sequence with a fixed-size window in a base--by-base fashion, calculates the AT/GC (with N) content within the window and mask out the bases and changes the qual accordingly if GC content is high or low.

Now my suggestion:
Since there are many reads to be checked against, I would imagine this could be a hot spot. So making it fast is probably a worthy thing to do:

instead of summing up base counts for each window, why not store a previousATNcount and previousGCNcount, considering that each time we step into a new window, these two variables are going to change only +/-1. This would get rid of the basesCount array;

as you have already did in the current implementation, the changes in bases count would be +/-1, why not, instead of setting bases and quals in each window, simply mark if the current window defined by start and end should be masked, and finally maximally extend the start and end, then go change those bases and quals

in maskWindow, instead of closed interval, why not make it semi-open [start, end)? this would also get rid of some of the "-1"'s in the apply() function.

instead of summing up base counts for each window, why not store a previousATNcount and previousGCNcount, considering that each time we step into a new window, these two variables are going to change only +/-1. This would get rid of the basesCount array;

We could use two variables instead of an array, but I'm not sure if it's faster. I prefer using the array because it avoids a couple of extra if statements.

as you have already did in the current implementation, the changes in bases count would be +/-1, why not, instead of setting bases and quals in each window, simply mark if the current window defined by start and end should be masked, and finally maximally extend the start and end, then go change those bases and quals

This could save a lot of time. Done.

in maskWindow, instead of closed interval, why not make it semi-open [start, end)? this would also get rid of some of the "-1"'s in the apply() function.

Done.

SHuang-Broad · 2017-07-28T17:32:12Z

src/main/java/org/broadinstitute/hellbender/transformers/IlluminaAdapterTrimTransformer.java

+                    //We have a match. Clip from the beginning of the adapter to the end of the read.
+                    final String basesString = read.getBasesString();
+                    final String qualString = new String(read.getBaseQualities());
+                    read.setBases(basesString.substring(0, alignmentStart).getBytes());


This part confuses me: by "clipping", do you mean set the bases to 'N' and qual to 2? In addition, a read should have a CIGAR, and clipping operations should update the CIGAR (either soft or hard), otherwise the object is in a incoherent state.

Hard clip. Good point about cigars. I found a handy ReadClipper class to use that will take care of this.

SHuang-Broad · 2017-08-02T16:54:01Z

👍

mwalker174 requested review from TedBrookings and SHuang-Broad July 27, 2017 15:47

TedBrookings approved these changes Jul 27, 2017

View reviewed changes

SHuang-Broad requested changes Jul 28, 2017

View reviewed changes

SHuang-Broad assigned mwalker174 Jul 28, 2017

mwalker174 added 3 commits July 31, 2017 12:53

PathSeq Illumina adapter trimming and simple repeat masking

c761b90

Update PathSeqFilterSpark documentation

be1a9b3

Addressed reviewer comments

09d2929

mwalker174 force-pushed the mw_adapter_repeat branch from cafde17 to 09d2929 Compare July 31, 2017 16:54

mwalker174 mentioned this pull request Jul 31, 2017

Use ReadClipper in BaseQualityClipReadTransformer #3388

Merged

mwalker174 closed this Jul 31, 2017

mwalker174 reopened this Jul 31, 2017

Removed unused variable MASKED_BASE_QUAL in SimpleRepeatMaskTransformer

5ac37a9

SHuang-Broad approved these changes Aug 2, 2017

View reviewed changes

mwalker174 merged commit c5b3c9a into master Aug 2, 2017

mwalker174 deleted the mw_adapter_repeat branch August 2, 2017 21:03

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

PathSeq Illumina adapter trimming and simple repeat masking #3354

PathSeq Illumina adapter trimming and simple repeat masking #3354

mwalker174 commented Jul 24, 2017

codecov-io commented Jul 24, 2017 •

edited

Loading

TedBrookings left a comment

TedBrookings Jul 27, 2017

mwalker174 Jul 31, 2017

mwalker174 Jul 31, 2017

TedBrookings Jul 27, 2017

mwalker174 Jul 31, 2017

TedBrookings Jul 27, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

TedBrookings Jul 27, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad left a comment

SHuang-Broad Jul 28, 2017

SHuang-Broad Jul 28, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017 •

edited

Loading

SHuang-Broad Jul 28, 2017

mwalker174 Jul 31, 2017

SHuang-Broad commented Aug 2, 2017

PathSeq Illumina adapter trimming and simple repeat masking #3354

PathSeq Illumina adapter trimming and simple repeat masking #3354

Conversation

mwalker174 commented Jul 24, 2017

codecov-io commented Jul 24, 2017 • edited Loading

Codecov Report

TedBrookings 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

SHuang-Broad 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mwalker174 Jul 31, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SHuang-Broad commented Aug 2, 2017

codecov-io commented Jul 24, 2017 •

edited

Loading

mwalker174 Jul 31, 2017 •

edited

Loading