Gap penalties check, options.c line 246

[lmjakt]

Summary:

Suggest changing:

if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)

to

if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e >= mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)

to avoid triggering error when using option -E2 or -E2,2.

And to consider modifying the error message for this operation.

Details

Use of alignment option -E2,2 or -E2 results in:

[ERROR] dual gap penalties violating E1>E2 and O1+E1<O2+E2

This is triggered by the following conditional on line 246 of options.c:

if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2))

Where the mo members q, q2, e, e2 correspond to O1, O2, E1 and E2 respectively (if I understand correctly).

(mo->q != mo->q2 || mo->e != mo->e2)

is TRUE since the default O1 and O2 values are 4 and 24.

!(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)

is also TRUE since E1 == E2, which means that the second half will not be evaluated and hence doesn't matter in this scenario.

If the conditional was mo->e >= mo->e2, then mo->q + mo->e < mo->q2 + mo->e2 would be evaluated, and be TRUE since (2 + 4 < 2 + 24). Negating the second half will then give an overall FALSE value and the error should not be triggered.

The error message seems misleading to me. The message states that E1 > E2, which in this case is not TRUE for -E2,2. In the default scenario opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e1 = 1, E1 is larger than E2 but this does not trigger the error. The error message also states that O1+E1<O2+E2, but it seems to me that the opposite condition is checked for due to the negation of the AND result. And in any case, under the default scenario this is TRUE, (2 + 4 < 1 + 24) and this is clearly not an error.

I'm not completely sure as to the intent of the check, but regardless, it seems that the error message should use 'or' rather than 'and', since the negation of an AND combination will be TRUE if either conditional is FALSE. Changing the error message to more clearly reflect the conditional:

[ERROR] dual gap penalties violating !(E1>E2 and O1+E1<O2+E2)

would seem more accurate; but perhaps a bit cryptic to the normal user.

My reason for trying -E2,2 is that I'm looking at pore-C data and finding many alignments with long indels. These indels may be caused by religation of nearby sites and I wish to avoid this.

Somewhat unrelated; I very much appreciate minimap2 (and htslib). Not only because it works very well, but also because the code for it is nicely structured and clearly written making it easy to reason about things. And that is not that common.

Thanks,

Martin

[lh3]

The test is intentional. When you use -O4,4 -E2,2 or simply -O4 -E2, minimap2 will use the standard affine gap penalty which is faster than dual-affine. As we are talking about source code, see these lines:

minimap2/align.c

Lines 356 to 360 in 6d49eb6

	} else if (opt->q == opt->q2 && opt->e == opt->e2) { // affine gap
	ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, ksw_flag, ez);
	} else { // dual affine gap
	ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, ksw_flag, ez);
	}

Nonetheless, I agree that the error message can be improved. I will leave this issue open and address when I make major changes to minimap2 next time.

[lh3]

BTW, thanks for the caution on pore-C. That may indeed happen and I haven't thought about it before. Thinking more about pore-C, I guess you probably want to disable long-range chaining as well with -r500,500.

[lmjakt]

Thanks for the quick response and the suggestion about the long-range chaining.

I'm not sure whether it would cause any other problems, or if it would be advisable, but I think my suggested change would work if combined with another one character change on line 356. The only real advantage would be to allow specifying -Ex without also giving -Oy and to make it somewhat easier to express the requirements in the manual (see details below).

The test is intentional. When you use -O4,4 -E2,2 or simply -O4 -E2, minimap2 will use the standard affine gap penalty which is faster than dual-affine. As we are talking about source code, see these lines:

minimap2/align.c

Lines 356 to 360 in 6d49eb6

	} else if (opt->q == opt->q2 && opt->e == opt->e2) { // affine gap
	ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, ksw_flag, ez);
	} else { // dual affine gap
	ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, ksw_flag, ez);
	}

If line 356 is changed to:

} else if (opt->q <= opt->q2 && opt->e == opt->e2) {

Then it will still do the correct thing even if opt->q < opt->q2 as the condition guarantees that q and e will always give a smaller (or equal) penalty than q2 and e2 (and opt->q should never be larger than opt->q2).

This would allow the use of -Ex (or -Ex,x) without also specifying -Oy. You would still not be able to use just -Oy without specifying the extension penalties, but that is implicit from the requirement:

(e >= e2) && (e + q <= e2 + q2)

And that is easy to express in the manual:

E INT1[,INT2]
Gap extension penalty [2,1]. A gap of length k costs min{O1+kE1,O2+kE2}. In the splice mode, the second gap penalties are not used.
Note that $E1 \geq E2$ AND $O1+E1 \leq O2+E2$ must be TRUE.

This suggests (to me) that:

minimap2/options.c

Line 246 in 6d49eb6

if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {

could be simplified to:

if (mo->e < mo->e2 || mo->q + mo->e > mo->q2 + mo->e2) {

if line 356 in align.c is changed (as suggested above) to:

} else if (opt->q <= opt->q2 && opt->e == opt->e2) {

As long as there are no other consequences of allowing e == e2 && q < q2.

However, I think the only advantage to making these changes is that it makes it easier to express the requirements in the manual (to avoid, "if E1 and E2 are the same then O1 must equal O2" in addition to specifying the required inequalities). I find combinatorial inequalities hard to reason about, and easy to get wrong, so changing expressions in code that works is possibly not a good idea.

Nonetheless, I agree that the error message can be improved. I will leave this issue open and address when I make major changes to minimap2 next time.

Thanks. It might also be an idea to mention the requirement in the manual (which I think would be simpler given the suggested changes). I can, if it would be useful, make the changes (warnings, code or manual) and send a pull request.