SamPairUtil.setMateInformationOnSupplementalAlignment copies TLEN from the mate primary instead of computing it

[nh13]

setMateInformationOnSupplementalAlignment derives the supplementary's TLEN by negating the mate primary's (SamPairUtil.java#L354):

supplemental.setInferredInsertSize(-matePrimary.getInferredInsertSize());

This assumes the supplementary sits where its own primary sits — same reference, same side of the mate. For split alignments, the only records this method is called on, that is false by construction. The sibling path does it correctly: setMateInfo calls computeInsertSize, which guards unmapped (L208) and cross-reference (L211). The supplementary path never calls it.

Reproduced against htsjdk 5.0.0:

--- cross-reference supplementary ---
R1      RNAME=chr17  POS=66798979    RNEXT=chr17  PNEXT=66798853    TLEN=1
R2      RNAME=chr17  POS=66798853    RNEXT=chr17  PNEXT=66798979    TLEN=-1
R1supp  RNAME=chrX   POS=153688036   RNEXT=chr17  PNEXT=66798853    TLEN=1

--- same-reference supplementary beyond its mate ---
R1      RNAME=chr1   POS=1000        RNEXT=chr1   PNEXT=1350        TLEN=450
R2      RNAME=chr1   POS=1350        RNEXT=chr1   PNEXT=1000        TLEN=-450
R1supp  RNAME=chr1   POS=5000        RNEXT=chr1   PNEXT=1350        TLEN=450

Two distinct problems. The first supplementary has RNAME != RNEXT yet a nonzero TLEN. The second spans 5,000–5,099 against a mate at 1,350–1,449, making it the rightmost segment — so the spec's sign rule requires a negative value, and 450 is the primary pair's insert size, describing coordinates the supplementary does not occupy.

Scope

Behaviour unchanged since 9e03608 (2014), which added the method alongside the AbstractAlignmentMerger caller it was written for. Affects Picard MergeBamAlignment and fgbio (Bams.fixMateInfo, ClipBam, TrimPrimers).

Not covered by tests: SamPairUtilTest.testSetMateInfoMateCigarOnSupplementals asserts only the mate CIGAR, never TLEN, and makeSamRecord hardcodes setReferenceIndex(0), so the cross-reference case cannot arise.

Spec status

The SAM spec is silent here rather than violated — hts-specs #522 scoped the TLEN definition to primary reads and left non-primary records undefined, and #842 leaves the computation aligner-defined. So this is an internal-consistency argument, plus divergence from other implementations.

bwa-mem (bwamem.c:887-892) and minibwa (format.c:243-262) both compute TLEN per emitted record from that record's own 5′ position against the mate primary, and emit 0 when the two are on different references — including for supplementary records. minibwa is the closest structural analogue to this method, taking the same (record, mate primary) pair as input:

int tlen = 0;
if (this_tid >= 0 && r_next) {
    if (this_tid == r_next->tid) {                              // cross-reference guard
        if (r) {
            int this_pos5 = r->rev? r->te - 1 : this_pos;       // the record's own coordinates
            int next_pos5 = r_next->rev? r_next->te - 1 : r_next->ts;
            tlen = next_pos5 - this_pos5;
        }

samtools fixmate leaves supplementaries untouched, and htslib's CRAM decoder uses chain extents. Only htsjdk copies the value from a different record.

Proposed fix

Route through computeInsertSize, which already handles both cases:

supplemental.setInferredInsertSize(computeInsertSize(supplemental, matePrimary));

This keeps htsjdk internally consistent and matches bwa-mem/minibwa. I would deliberately not adopt htslib's chain-wide leftmost-to-rightmost basis, since that would also change the primaries' TLEN whenever a supplementary falls outside the pair's span.

Happy to open a PR if the approach looks right.

---

Investigation and reproduction assisted by Claude Code (Anthropic). All code references and the htsjdk 5.0.0 reproduction above were verified by hand.