Fixing `MergeBamAlignment creates invalid BAMs in the presence of secondary alignments` #932

DenisVerkhoturov · 2017-09-19T12:48:58Z

Description

The problem originates in HitsForInserts#coordinateByHitIndex. This method makes decision whether read is paired or not on assumption that all the records have HI-tag but ignores other criteria such as RNEXT (Ref. name of the mate/next read) and PNEXT (Position of the mate/next read) SAM fields. So paired reads are determined as two unrelated ones in the case when HI-tag is absent.

Implementation

Added new method (coorditanteByMate) for correlating pair-aware alignments. It does not require presence of HI-tags and performs work based on the data contained in fields: RNAME, RNEXT, POS, PNEXT.
Added MergeBamAlignmentTest#testPairedReadsAreProcessedProperlyInAbsenceOfHitIndexTags running MergeBamAlignment with SAM file contained pair-aware alignments without HI-tag.
Added testCoordinateByMateMethodWorksProperly for the new method.
Changed BestMapqPrimaryAlignmentSelectionStrategy. It uses implemented method coordinateByMate instead of coordinateByHitIndex now. We kept the old method due to backward compatibility.

Performance

All the benchmarks were measured by JMH (Mode - thrpt, Cnt - 100, Units - ops/s)

Benchmark		100 paired alignments		500 paired alignments		10_000 paired alignments
Method	Case	Score	Error	Score	Error	Score	Error
coordinateByHitIndex	HI-tag and mate is absent	272_334.682	± 3852.273	27_056.185	± 543.235	95.384	± 0.354
	HI-tag is absent	282_456.059	± 5726.669	35_603.756	± 954.916	696.495	± 7.419
	HI-tag is present	280_571.088	± 1900.595	27_734.340	± 914.012	95.101	± 0.935
coordinateByMate	HI-tag and mate is absent	180_912.877	± 1983.565	17_965.070	± 242.398	55.060	± 2.967
	HI-tag is absent	174_695.474	± 1786.275	17_139.180	± 128.948	51.079	± 0.542
	HI-tag is present	185_531.494	± 1975.027	18_120.429	± 157.350	57.640	± 3.658

New method is 1.7 times slower but this caused by using required fields that allow rearranging alignments in any situations. And the only situation `coordinateByIndex` is much faster is when we have no **Hi-tag** and pair-aware alignments are being rearranged incorrectly (case described in issue).

Resolves

MergeBamAlignment creates invalid BAMs in the presence of secondary alignments

Checklist (never delete this)

Never delete this, it is our record that procedure was followed. If you find that for whatever reason one of the checklist points doesn't apply to your PR, you can leave it unchecked but please add an explanation below.

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

coveralls · 2017-09-19T13:08:20Z

Coverage decreased (-0.03%) to 76.734% when pulling 31d042e on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into 9a7ba13 on broadinstitute:master.

pshapiro4broad

Just a few comments about the Java code, not about the code's meaning.

pshapiro4broad · 2017-09-19T13:48:47Z

src/main/java/picard/sam/HitsForInsert.java

+            final Integer hiTagValueOfFirst = firstOfPairOrFragment.get(i).getIntegerAttribute(SAMTag.HI.name());
+            final Integer hiTagValueOfSecond = secondOfPair.get(i).getIntegerAttribute(SAMTag.HI.name());
+
+            if (nonNull(hiTagValueOfFirst) && nonNull(hiTagValueOfSecond)) {


What is the rationale for using these functions instead of comparing against null directly? My understanding is that these methods are provided so we have convenient lambda forms for streaming operations, etc. Using them inline like this doesn't look like an improvement to me.

I prefer to use functional style instead of procedural even if it could decrease performance a little bit because it usually increases code readability. Especially for functions of Java standard library. I could revert all cases when I use static import if it is the concern. So is it critical?

Not sure what you mean about functional vs procedural... one form is an expression consisting of a function call, the other is a binary expression using the!= operator. Both are equally functional and procedural, as I understand these terms.

Readability is a subjective argument; to me, using != null is more readable. Given that it requires a static import to avoid making the code much wordier than the original, I don't think it's worth the change.

Speaking about operators vs functions, I have always thought that since Java doesn't allow to make partial application of an operator it is convenient to use functions in such cases. Let's imagine we have equals function that takes two objects and checks their equality. So isNull function is just a partial application of equals function with fixed null as one of the arguments. Then nonNull is negated isNull. So that way, for the sake of consistency it is easier to use functions instead of operators in Java not only in streams but everywhere in mine opinion. Is it a little bit excessively? 😄

It sounds like you have a more enlightened view of functional programming in Java than I do. :)

Thanks for the explanation, I now understand the distinction you're making between expressions and functions, and why it's more consistent to use functions everywhere. I can't say I'm convinced but I'll definitely be giving it more thought.

pshapiro4broad · 2017-09-19T13:50:21Z

src/main/java/picard/sam/HitsForInsert.java

 import java.util.Comparator;
 import java.util.List;
+import java.util.Objects;
+
+import static java.util.Objects.isNull;


In general static imports should be avoided as they make code harder to read.

Replaced by operators as was discussed here.

pshapiro4broad · 2017-09-19T13:51:32Z

src/main/java/picard/sam/HitsForInsert.java

+            final Integer hiTagValueOfSecond = secondOfPair.get(i).getIntegerAttribute(SAMTag.HI.name());
+
+            if (nonNull(hiTagValueOfFirst) && nonNull(hiTagValueOfSecond)) {
+                if (hiTagValueOfFirst < hiTagValueOfSecond) secondOfPair.add(i, null);


The original code didn't do this, but the new code should use { } around all if/else body statements, even if they're a single line.

pshapiro4broad · 2017-09-19T13:55:43Z

src/main/java/picard/sam/HitsForInsert.java

+        return nonNull(first) && nonNull(second)
+                && first.getMateAlignmentStart() == second.getAlignmentStart()
+                && first.getAlignmentStart() == second.getMateAlignmentStart()
+                && Objects.equals(first.getReferenceName(), second.getMateReferenceName())


Why use Object.equals() here? is it possible that one of these is null?

Since both SamRecod#mReferenceName and SamRecord#mMateReferenceName are nullable we need to ensure that they are non-null or simply delegate it to Objects.equals to avoid NullPointerException. This is a theoretical case that may hardly ever occur in runtime. So, do we need to prevent this case or can I unwrap this check?

It's OK, I was wondering if this was added to fix a known problem or was just precautionary.

Looking at SamRecord I see that it's common to check these fields for null so it seems wise to do it here too.

It would be great if picard had nullable constraints in the code base! That would be a huge PR though.

Optional can be rather helpful here and it would be just great if we could use Try and Lazy types from vavr to reduce service-code and to clarify Picard logic!=D

pshapiro4broad · 2017-09-19T13:56:24Z

src/main/java/picard/sam/HitsForInsert.java

+    /**
+     * @return true if records belong to the same pairwise alignment
+     */
+    private static boolean arePair(SAMRecord first, SAMRecord second) {


This looks like it should be a non-static method on SAMRecord, as it only deals with fields in that object and takes that object as its argument.

SAMRecord is part of htsjdk and this was the only constraint I haven't done it. Could I just make this function a SamRecord#arePair method and make associated pull-request to htsjdk?

Could you please review associated request to htsjdk?

coveralls · 2017-09-20T11:09:26Z

Coverage decreased (-0.04%) to 76.726% when pulling e418168 on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into 9a7ba13 on broadinstitute:master.

coveralls · 2017-09-20T11:58:15Z

Coverage decreased (-0.04%) to 76.726% when pulling 11a9c1f on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into 9a7ba13 on broadinstitute:master.

coveralls · 2017-09-20T11:58:15Z

Coverage decreased (-0.04%) to 76.726% when pulling 11a9c1f on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into 9a7ba13 on broadinstitute:master.

coveralls · 2017-09-20T16:26:25Z

Coverage decreased (-0.03%) to 76.73% when pulling 96bc12a on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into 9a7ba13 on broadinstitute:master.

yfarjoun · 2017-09-22T07:11:10Z

@DenisVerkhoturov there are many "esthetic" changes mixed in with substantive changes. Could you please rebase and separate into two commits, one that is purely esthetic and one that is substantive? it will make it much easier to understand what you are doing.

DenisVerkhoturov · 2017-09-22T10:27:17Z

@yfarjoun Sure, done=)

coveralls · 2017-09-22T10:46:14Z

Coverage decreased (-0.03%) to 77.699% when pulling e4c8c0d on EpamLifeSciencesTeam:epam-ls_MergeBamAlignment into b63d73e on broadinstitute:master.

tbl3rd · 2017-12-13T16:58:09Z

src/main/java/picard/sam/HitsForInsert.java

+    }
+
+    /**
+     * Identifies weather records present pairwise alignment or not.


weather -=> whether

Thank you, fixed here.

tbl3rd · 2017-12-13T17:01:35Z

src/main/java/picard/sam/HitsForInsert.java

+    void coordinateByMate() {
+        final List<SAMRecord> newFirstOfPairOrFragment = new ArrayList<>();
+        final List<SAMRecord> newSecondOfPair = new ArrayList<>();
+


Looks like this file preferred ++i to i++.

As you might see (here for instance), there was a usage of a postfix increment in this file before. So I just followed this example.
All increment operators in this file changed to prefix form.

tbl3rd · 2017-12-13T17:37:36Z

src/test/java/picard/sam/MergeBamAlignmentTest.java

+
+            hitsForInsert.addFirstOfPairOrFragment(second);
+        }
+


Is 3 some function of HitsForInsert::numHits?

Thank you, fixed here.

tbl3rd · 2017-12-13T17:44:34Z

src/main/java/picard/sam/HitsForInsert.java


-        // Now renumber any correlated alignments, and remove hit index if no correlated read.


Thank you for breaking this out and cleaning up hi, and so on.

Thank you for your feedback=)

…ary alignments

DenisVerkhoturov changed the title ~~fix: MergeBamAlignment creates invalid BAMs in the presence of second…~~ Fixing MergeBamAlignment creates invalid BAMs in the presence of secondary alignments Sep 19, 2017

pshapiro4broad reviewed Sep 19, 2017

View reviewed changes

DenisVerkhoturov force-pushed the epam-ls_MergeBamAlignment branch from e418168 to 11a9c1f Compare September 20, 2017 11:36

DenisVerkhoturov mentioned this pull request Sep 21, 2017

Sam record pair identification depends on required SAM fields samtools/htsjdk#994

Closed

5 tasks

DenisVerkhoturov force-pushed the epam-ls_MergeBamAlignment branch from 96bc12a to 58082fa Compare September 22, 2017 10:25

tbl3rd reviewed Dec 13, 2017

View reviewed changes

DenisVerkhoturov added 3 commits January 25, 2018 13:32

fix: MergeBamAlignment creates invalid BAMs in the presence of second…

d9efb8e

…ary alignments

refactor: Fix javadoc, diamond operator, java 8 comparator

340248d

fix: Fixes concerning @tbl3rd's review

e4c8c0d

DenisVerkhoturov force-pushed the epam-ls_MergeBamAlignment branch from 58082fa to e4c8c0d Compare January 25, 2018 13:20

gbggrant requested a review from yfarjoun February 13, 2018 14:32

yfarjoun self-assigned this Jun 4, 2018

Merge branch 'master' into epam-ls_MergeBamAlignment

90bce0b

jacarey mentioned this pull request Jul 29, 2020

Expose INCLUDE_SECONDARY_ALIGNMENTS for MergeBamAlignment task. fulcrumgenomics/dagr#381

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixing `MergeBamAlignment creates invalid BAMs in the presence of secondary alignments` #932

Fixing `MergeBamAlignment creates invalid BAMs in the presence of secondary alignments` #932

DenisVerkhoturov commented Sep 19, 2017 •

edited

Loading

coveralls commented Sep 19, 2017 •

edited

Loading

pshapiro4broad left a comment

pshapiro4broad Sep 19, 2017

DenisVerkhoturov Sep 20, 2017

pshapiro4broad Sep 20, 2017

DenisVerkhoturov Sep 20, 2017

DenisVerkhoturov Sep 21, 2017

pshapiro4broad Sep 21, 2017

pshapiro4broad Sep 19, 2017

DenisVerkhoturov Sep 20, 2017

pshapiro4broad Sep 19, 2017

DenisVerkhoturov Sep 20, 2017 •

edited

Loading

pshapiro4broad Sep 19, 2017

DenisVerkhoturov Sep 20, 2017

pshapiro4broad Sep 20, 2017

DenisVerkhoturov Sep 20, 2017

pshapiro4broad Sep 19, 2017

DenisVerkhoturov Sep 20, 2017

DenisVerkhoturov Sep 21, 2017

coveralls commented Sep 20, 2017 •

edited

Loading

coveralls commented Sep 20, 2017

coveralls commented Sep 20, 2017

coveralls commented Sep 20, 2017 •

edited

Loading

yfarjoun commented Sep 22, 2017

DenisVerkhoturov commented Sep 22, 2017

coveralls commented Sep 22, 2017 •

edited

Loading

tbl3rd Dec 13, 2017

DenisVerkhoturov Jan 25, 2018

tbl3rd Dec 13, 2017

DenisVerkhoturov Jan 25, 2018

tbl3rd Dec 13, 2017

DenisVerkhoturov Jan 25, 2018

tbl3rd Dec 13, 2017

DenisVerkhoturov Jan 25, 2018


		// Now renumber any correlated alignments, and remove hit index if no correlated read.

Fixing MergeBamAlignment creates invalid BAMs in the presence of secondary alignments #932

Are you sure you want to change the base?

Fixing MergeBamAlignment creates invalid BAMs in the presence of secondary alignments #932

Conversation

DenisVerkhoturov commented Sep 19, 2017 • edited Loading

Description

Implementation

Performance

Resolves

Checklist (never delete this)

Content

Review

coveralls commented Sep 19, 2017 • edited Loading

pshapiro4broad 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

DenisVerkhoturov Sep 20, 2017 • edited Loading

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

coveralls commented Sep 20, 2017 • edited Loading

coveralls commented Sep 20, 2017

coveralls commented Sep 20, 2017

coveralls commented Sep 20, 2017 • edited Loading

yfarjoun commented Sep 22, 2017

DenisVerkhoturov commented Sep 22, 2017

coveralls commented Sep 22, 2017 • edited Loading

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

Fixing `MergeBamAlignment creates invalid BAMs in the presence of secondary alignments` #932

Fixing `MergeBamAlignment creates invalid BAMs in the presence of secondary alignments` #932

DenisVerkhoturov commented Sep 19, 2017 •

edited

Loading

coveralls commented Sep 19, 2017 •

edited

Loading

DenisVerkhoturov Sep 20, 2017 •

edited

Loading

coveralls commented Sep 20, 2017 •

edited

Loading

coveralls commented Sep 20, 2017 •

edited

Loading

coveralls commented Sep 22, 2017 •

edited

Loading