New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed a bug in hardClipCigar function that caused incorrect cigar calculation #6280

Merged

yfarjoun merged 9 commits into master from yf_fix_bug_in_hardclip

Dec 2, 2019

Contributor

yfarjoun commented Nov 24, 2019

Problem was ocurring in the presence of insertions and deletions. fixes #6139

Changed ReadClipper unit tests:

The tests in many cases assumed that the unclipped alignment locations do not change when the read is clipped. This is not true: for example if start for cigar
1M1I3M is 100, the unclipped start for 2H3M is 99. All assertUnclipped calls were
removed
Alignment now check that the read length remains to be consistent with the CIGAR, that the aligned bases span are consistent with the CIGAR and that the number of clipped bases from the read is consistent with the requested clipping

Hard clipping in ClippingOps was buggy, thus we introduced new tests for it.
Text was refactored for readability
Clipping in ClippingOps did not treat insertions and deletions in the clipped parts of the CIGAR correctly. This was fixed
Alignment re-calculation after clipping did not work correctly if the initial CIGAR contained insertions and deletions
Hard clipping applied to the hard clipped read did not behave correctly

yfarjoun requested review from davidbenjamin and droazen

November 24, 2019 11:45

gatk-bot commented Nov 24, 2019 •

edited

Loading

Travis reported job failures from build 28153
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28153.13	logs
unit	openjdk8	28153.3	logs

gatk-bot commented Nov 24, 2019 •

edited

Loading

Travis reported job failures from build 28155
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28155.13	logs
unit	openjdk8	28155.3	logs

davidbenjamin requested changes

View reviewed changes

Contributor

davidbenjamin left a comment

Assorted complaints. The fix looks good and the tests are compelling.

src/main/java/org/broadinstitute/hellbender/utils/clipping/ClippingOp.java Outdated

                       readCopied.setBaseQualities(newQuals);
                   }
                   private void applyWriteNs(final GATKRead readCopied) {
                       final byte[] newBases = readCopied.getBases();       //this makes a copy so we can modify in place
-                      overwriteFromStartToStop(newBases, (byte)'N');
+                      overwriteFromStartToStop(newBases, (byte) 'N');

Contributor

davidbenjamin Nov 25, 2019

Personally, I would err on the side of rigor and use Nucleotide.N.encodeAsByte() here.

Contributor Author

yfarjoun Nov 26, 2019

I added a space, and you want me to fix the old code!?? :-)

src/main/java/org/broadinstitute/hellbender/utils/clipping/ClippingOp.java Outdated

                           if (curReadCounter + curReadLength > readBasesBeforeReference) {
                               curReadLength = readBasesBeforeReference - curReadCounter;
                               curRefLength = readBasesBeforeReference - curReadCounter;
                               truncated = true;
+                          } else {

Contributor

davidbenjamin Nov 25, 2019

Make the definition of truncated more explicit:

final boolean truncated = curReadCounter + curReadLength > readBasesBeforeReference;
if (truncated) {
  curReadLength = readBasesBeforeReference - curReadCounter;
  curRefLength = readBasesBeforeReference - curReadCounter;
}

src/main/java/org/broadinstitute/hellbender/utils/clipping/ClippingOp.java Outdated

-                          if (runAsserts) {
-                              assert newCigar.isValid(null, -1) == null;
-                          }
+                          assert !runAsserts || newCigar.isValid(null, -1) == null;

Contributor

davidbenjamin Nov 25, 2019 •

edited

Loading

Asserts in non-test code are not very GATK. Is this really not covered in unit tests? You can defer this since it's tangential to this PR.

Contributor Author

yfarjoun Nov 26, 2019

reverted to keep the history cleaner

src/main/java/org/broadinstitute/hellbender/utils/clipping/ClippingOp.java Outdated

                               } else if (!readHasStarted && cigarElement.getOperator() == CigarOperator.HARD_CLIP) {
                                   totalHardClip += cigarElement.getLength();
                               } else if (!readHasStarted && cigarElement.getOperator() == CigarOperator.DELETION) {
-                                  totalHardClip += cigarElement.getLength();
+                                  //totalHardClip += cigarElement.getLength();

Contributor

davidbenjamin Nov 25, 2019

It seems like you are putting in an effect-less line of code to document that the total is not incremented. Is that because there are a bunch of else ifs and you don't want to give the impression that it's always incremented? How about this:

if (!readHasStarted && (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.CigarOperator.SKIPPED_REGION) {
  totalHardClip += cigarElement.getLength();
}
readHasStarted |= cigarElement.getOperator() != CigarOperator.DELETION &&
                        cigarElement.getOperator() != CigarOperator.SKIPPED_REGION &&
                        cigarElement.getOperator() != CigarOperator.HARD_CLIP;

Also, in the (unmodified) code below it seems that the if (ReadHasStartedBlock) always results in addedHardClips == true. If you agree, could you move the addedHardCips = true statements to the end of that block?

src/main/java/org/broadinstitute/hellbender/utils/clipping/ClippingOp.java Outdated

                           size += elements.get(i).getLength();
                           i++;
                       }
                       return size;
                   }
+                  /**
+                  * Calculates shift of alignment int the newCigar relative to the oldCigar due to

Contributor

davidbenjamin Nov 25, 2019

int -> in

src/main/java/org/broadinstitute/hellbender/utils/read/ArtificialReadUtils.java Outdated

-                      byte [] bases = Utils.dupBytes(base, length);
-                      byte [] quals = Utils.dupBytes(qual, length);
+                      final int length = cigar.getReadLength();
+                      final Random random = new Random(TestUtil.RANDOM_SEED);

Contributor

davidbenjamin Nov 25, 2019

What about Utils.getRandomGenerator() or some such Utils method?

src/test/java/org/broadinstitute/hellbender/utils/clipping/ReadClipperUnitTest.java Outdated

@@ @@ -87,19 +97,24 @@ public void testHardClipByReferenceCoordinates() { @@
                       for (Cigar cigar : cigarList) {
                           GATKRead read = ReadClipperTestUtils.makeReadFromCigar(cigar);
                           int start = getSoftStart(read);
+                          int aln_start = read.getStart();

Contributor

davidbenjamin Nov 25, 2019

camel case

src/test/java/org/broadinstitute/hellbender/utils/clipping/ReadClipperUnitTest.java Outdated

+                  private void assertReadClippingConsistent(final GATKRead original,
+                                                            final GATKRead clipped,
+                                                            final int clipping) {
+                      int clip_diff = original.getLength() - clipped.getLength();

Contributor

davidbenjamin Nov 25, 2019

camel case

src/test/java/org/broadinstitute/hellbender/utils/clipping/ReadClipperUnitTest.java Outdated

+                          Assert.assertTrue(clipped.isUnmapped());
+                          return;
+                      }
+                      int delta_left = clipped.getStart() - original.getStart();

Contributor

davidbenjamin Nov 25, 2019

camels

src/test/java/org/broadinstitute/hellbender/utils/clipping/ReadClipperUnitTest.java Outdated

+                      Assert.assertEquals(delta_left, clipped_left);
+                      int original_end = original.getEnd() ;
+                      int clipped_end = clipped.getEnd() ;
+                      int delta_right = original_end - clipped_end;

Contributor

davidbenjamin Nov 25, 2019

camels

droazen removed their request for review

November 25, 2019 19:12

droazen assigned davidbenjamin

Yossi Farjoun added 3 commits

November 26, 2019 16:39


           Fixed a bug in hardClipCigar function that caused incorrect cigar ca…

b8f6240

…lculation

1. Changed ReadClipper unit tests:
  - The tests in many cases assumed that the unclipped alignment locations do not    change when the read is clipped. This is not true: for example if start for cigar
   1M1I3M is 100, the unclipped start for 2H3M is 99. All assertUnclipped calls were
   removed
   - Alignment now check that the read length remains to be consistent with the CIGAR, that the aligned bases span are consistent with the CIGAR and that the number of clipped bases from the read is consistent with the requested clipping
2. Hard clipping in ClippingOps was buggy, thus we introduced new tests for it.
3. Text was refactored for readability
4. Clipping in ClippingOps did not treat insertions and deletions in the clipped parts of the CIGAR correctly. This was fixed
5. Alignment re-calculation after clipping did not work correctly if the initial CIGAR contained insertions and deletions
6. Hard clipping applied to the hard clipped read did not behave correctly


          - make test bases random

f21512b


          Responding to review comments

3569d07

yfarjoun force-pushed the yf_fix_bug_in_hardclip branch from 780f67b to 3569d07 Compare

November 26, 2019 14:40

gatk-bot commented Nov 26, 2019 •

edited

Loading

Travis reported job failures from build 28178
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28178.13	logs
unit	openjdk8	28178.3	logs

Yossi Farjoun added 5 commits

November 26, 2019 22:26


          - cleanup. and TODO. need fixing before merge.

14d632b


          - opened up convoluted while loop to clarify logic in ClippingOp

e3219e5

- added tests
- fixed bug that was found in the process


          -optimize import

e54e3d8


          - oops

90d7c3d


          - spaces

aad0d75

Contributor Author

yfarjoun commented Nov 27, 2019

builds are now failing due to the changes in the artifical reads generator...I'll make them more flexible so that I can still get more random reads without breaking the other tests....

gatk-bot commented Nov 27, 2019 •

edited

Loading

Travis reported job failures from build 28207
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28207.13	logs
unit	openjdk8	28207.3	logs

gatk-bot commented Nov 27, 2019 •

edited

Loading

Travis reported job failures from build 28205
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28205.13	logs
unit	openjdk8	28205.3	logs

gatk-bot commented Nov 27, 2019 •

edited

Loading

Travis reported job failures from build 28209
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28209.13	logs
unit	openjdk8	28209.3	logs


          - revert random bases in getArtificialReads

d1db393

Contributor Author

yfarjoun commented Dec 1, 2019

back to you @davidbenjamin The tests are also passing this time!

davidbenjamin approved these changes

View reviewed changes

Contributor

davidbenjamin left a comment

It's ready to merge.

yfarjoun merged commit 0a4cca8 into master

yfarjoun deleted the yf_fix_bug_in_hardclip branch

December 2, 2019 07:57

Contributor Author

yfarjoun commented Dec 2, 2019

Thanks for the review @davidbenjamin 🎉

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet