max-intronlen controls "concordant pair" alignment flags - that should not be the case

[therealgenna]

I use hisat2 version 2.2.1.

I've noticed that when I set small --max-intronlen some of my longer read pairs (fragments) show as discordantly mapped. They have flags 81 and 161, instead of 83 and 163 [missing the read mapped in proper pair (0x2) flag], as well as YT:Z field is set to DP instead of CP.

If I increase --max-intronlen sufficiently, the alignment is identical, with the correct flags that the pair is properly mapped.

I think this is a bug, because this is just based on the unsequenced gap between the two reads; this gap does not have to be an intron! It could (and likely is) a genomic sequence, which would have been sequenced had I sequenced more bases in each read.

I could not find this issue documented, so I report it here.

Here's an example:
(Human genome hg38)

R1.fa:

@FRAG 1:N:0:0
GTGGCGGCTCCATAGCTCAGGGGTTAGAGCACTGGTCTTGTAAACCAGGGGTC
+
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCC

R2.fa:

@FRAG 2:N:0:0
AACAGTCTAGTTTCTTCCTCGCTTACTGCAGCGTACAAACACCTATTTTCCCA
+
CCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCCCCCCCCCC;CCCCCCCCCCCC

Issuing

hisat2 -x /path/to/human_gencode-v28 -1 R1.fa -2 R2.fa --no-unal --no-softclip --score-min L,0,-0.6 --min-intronlen 20 --max-intronlen 189 -S e.sam > e.log 2>&1; grep -v '^@' e.sam

yields

FRAG        81      chr14   20631185        60      53M     =       20630943        -295    GACCCCTGGTTTACAAGACCAGTGCTCTAACCCCTGAGCTATGGAGCCGCCAC   CCCCCCC;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   AS:i:0  XN:i:0  XM:i:0XO:i:0  XG:i:0  NM:i:0  MD:Z:53 YS:i:0  YT:Z:DP NH:i:1
FRAG        161     chr14   20630943        60      53M     =       20631185        295     AACAGTCTAGTTTCTTCCTCGCTTACTGCAGCGTACAAACACCTATTTTCCCA   CCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCCCCCCCCCC;CCCCCCCCCCCC   AS:i:0  XN:i:0  XM:i:0XO:i:0  XG:i:0  NM:i:0  MD:Z:53 YS:i:0  YT:Z:DP NH:i:1

while issuing (only changing --max-intronlen from 189 to 190)

hisat2 -x /path/to/human_gencode-v28 -1 R1.fa -2 R2.fa --no-unal --no-softclip --score-min L,0,-0.6 --min-intronlen 20 --max-intronlen 190 -S e.sam
 > e.log 2>&1; grep -v '^@' e.sam

yields identical alignment but also correctly setting "properly paired" flags:

FRAG        83      chr14   20631185        60      53M     =       20630943        -295    GACCCCTGGTTTACAAGACCAGTGCTCTAACCCCTGAGCTATGGAGCCGCCAC   CCCCCCC;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   AS:i:0  XN:i:0  XM:i:0XO:i:0  XG:i:0  NM:i:0  MD:Z:53 YS:i:0  YT:Z:CP NH:i:1
FRAG        163     chr14   20630943        60      53M     =       20631185        295     AACAGTCTAGTTTCTTCCTCGCTTACTGCAGCGTACAAACACCTATTTTCCCA   CCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCCCCCCCCCC;CCCCCCCCCCCC   AS:i:0  XN:i:0  XM:i:0XO:i:0  XG:i:0  NM:i:0  MD:Z:53 YS:i:0  YT:Z:CP NH:i:1

My reads are both 53nts, and CIGARs are all 53M. The insert size reported by IGV is 295 (20631185 + 53 - 1 - 20630943 + 1 = 295.)
The unsequenced gap between the two aligned reads is thus 295 - 53*2 = 189.
So hisat2 uses max intron length parameter to control allowable fragment size, even though the gap between R1 and R2 does not have to be an intron.