TLEN meaning when leftmost and rightmost are the same for a pair, or the same for a single record

[jkbonfield]

This came about from samtools/htslib#1948 in what I consider to be an implementation bug in my CRAM decoder.
See also #522 and #48, but I'm ignoring for now the issue of primary vs secondary and simply looking at primary alignments

Although this is arguably an implementation bug, while looking at it I also see it's a specification weakness too.

The CRAM specification defines how to compute TLEN when we have a read-pair stored as a pair (ie not "detached"). Here the reads are in the same slice and we can derive TLEN on the fly without storing it verbatim, but to do so we have to tightly define the meaning of TLEN. CRAM 3.1 section 10.4 defines this as:

27: Find leftmost and rightmost mapped coordinate in records this_record and next_record.
28: For leftmost of this_record and next_record: template_size ← rightmost − lef tmost + 1
29: For rightmost of this_record and next_record: template_size ← −(rightmost − lef tmost + 1)

All well and good, but what if leftmost and rightmost are the same? Ie a containment:

------------>
    <-----

That's obviously under-specified. A better definition may be to define how we assign the positive (rightmost-leftmost+1) value and use an else clause for the other alignments. We lack the statement that SAM has about the two ends having differing sign:

If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing sign.

This obviously makes the explicit assumption of "two ends". I guess we could define "ends" to be primary READ1 and READ2 and any middle reads (READ1+READ2) or supplementaries are excluded from that statement, but it's open to interpretation and I can clearly see cases where this is challenging when, for example, the supplementary or "middle" read happens to be outside the extents of the primary READ1/READ2 ends and so modifies the definition of leftmost and rightmost.

What about where we share the same start but differing ends, or the same end but differing starts? That's clearly not covered by "if segments cover the same coordinates", but it does confuse the definition of leftmost and rightmost. Similarly containments where the leftmost and rightmost are the same alignment and the "other end" is shorter and somewhere between those two points. Such cases don't happen for properly mapped data, but it can always be triggered by base calling errors and/or alignment failures. The specification needs to be clear.

If we have 3 reads, with a READ1+READ2, then that read could potentially be the leftmost or rightmost outright (not tied) due to misalignments. So maybe leftmost and rightmost there refer only to the READ1 and READ2 coordinates? Frankly triplets (or more) are just a disaster, but they could occur with supplementary alignments, so we have READ1, READ1(supp), READ2 for expample. (Ignored for now as the language is confusing enough.)

eg

----->  1+2
   -------> 1
        <-------- 2

In the above case, read 1 and 2 form a nice pair, but 1+2 is the leftmost and so we can (and do) get read1 and read2 having the same TLEN sign. Technically this is spec compliant, but I doubt it was intended to be. Sigh.

Focusing purely on pairs I can see a number of disparities between samtools, picard and noodles CRAM implementations and also samtools fixmates. I have a PR brewing for this, but I've been producing test data too. However it's also a SAM issue.

I have a work in progress in testing this, and some of my correct vs incorrect assignments are themselves wrong due to the nuance of "arbitrary" in the SAM spec above. I'll continue to fix this. For now however my on-going test files and spreadsheet of results are attached.

https://docs.google.com/spreadsheets/d/11pAe0I9794ItNKrPEw3N0F27LVPafr9Cq_F-RNBNkzs/edit?usp=sharing

TLEN.test.tar.gz

My thoughts on how to fix these issues are still in flux, but:

• For SAM: it's pretty good, but we need to be explicit about the containment case where leftmost and rightmost are the same read and we have an extra read which is neither leftmost nor rightmost. Even if it's something such as "If segments cover the same coordinates or the leftmost and rightmost coordinates are from the same segment, then..."

• For CRAM: we need to explain leftmost and rightmost in a far better fashion, again compensating for the situation of segments with the same start and/or end coordinates, and containments. As it's written in pseudocode, defining the positive one and using an else clause for others feels the most natural solution.

• We may also wish to add a recommendation about ties rather than "arbitrary". Eg if we sort a file we shouldn't expect the signs to flip, which is why htslib looks at BAM flags to resolve them. However it's implementation defined and you can be spec compliant without doing that (eg noodles) so maybe that's not necessary.

PS. Once I've explored this more I'll look into creating more test files for hts-specs, but it may need multiple "valid" decodes to cover the arbitrary definitions.

I'll probably also work on a samtools fixmates --tlen1 / --tlen2 option at some point so it can be spec compliant with the current definitions, and also fix it so it's not so broken on more than 2 reads, maybe. However that's an implementation issue.

[jkbonfield]

On further consideration:

If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs

The intent is clear - we must have different signs for the two "ends" even if we're choosing which is which arbitrarily due to equal coverage. However, what's not clear is what an "end" is. It talks about the mapped start of the template and the mapped end of the template, which could be interpretted purely as observations on where any segment maps, or it could be interpreted as having knowledge of the sequencing chemistry and computing where READ1 / READ2 extents go to.

It does talk about "primary alignments", but that is defined earlier to be a single "read alignment" out of a potential set of multiple alignments, and a "read alignment" is defined to be either a "linear alignment" or "chimeric alignment". So by definition, when TLEN talks about the start and end of a primary alignment it is including any chimeric alignments (ie supplementaries) that are part of the primary alignment. This totally makes sense, and is also clearly why it was listed in a chain or PNEXT rather than just "other end". So far so good (if complex).

However the inclusion of supplementaries now confuses the definition of what the "start of the template" and "end of the template" really means. It could strictly be defined as leftmost / rightmost mapping, or it could be restricted purely to the first and last base of the segment reported for READ1 and READ2. Note this doesn't necessarily mean non-supplementary either.

. Typically, one
of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the
others are called “supplementary” and are distinguished by the supplementary alignment flag. All the
SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80
flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

Hence we could have (1/2 is for READ1 and READ2, "=" is for non-supplementary and "~" is for supplementary, and "." is soft-clips):

1====>...                        <=====2
     ......1~>

or

1~>......                        <=====2
      ...1=====>

Here the aligner has used the decision to define the "representative" alignment as the longest match, but it is defined to be arbitrary and it could easily have also chosen to use the first / last matches as the representative one instead.

In both cases, the desired TLEN is clear and the same thing - all the way from left ro right.

---

Given the soft-clips, we know which bases of the original sequence reported from the sequencing instrument are being used in the alignment. So base 1(ish) may be interpreted to be the start of the alignment. Consider this case though:

......1~>                        <=====2
           1=====>...

Here the representative alignment is the longest one and it represents base 1 to most of the way through the called sequence. We have a short segment of softclip at the 3' end, and that short segment is recorded as a supplementary alignment. However that alignment is before the non-supplementary.

The definition of the template extents feels arbitrary and I think it's reasonable to say it's up to the aligner to decide whether it views the template extents as all the way from left to right, or only from the left end of non-supplementary read to the right. This is probably also mired in the leftmost/rightmost vs 5' ends argument which we've already decided is implementation specific.

---

What does this mean for TLEN?

• We shouldn't blithly ignore supplementaries that are part of the primary alignment (ie part of the RNEXT/PNEXT chain on a non-secondary alignment)
• However we also shouldn't state TLEN must be computed in this specific way and the calculation of it is up to the aligner.
• With that in mind, the exact interpretation of what records constitute an "end" is also subjective.

We can write our own implementations obviously (eg samtools fixmate), but that's outside the bounds of the specification.

Frankly what it's all coming down to is the specification can probably just throw its hands up and say "you're on your own!" lol.

When it comes to CRAM however, we painted ourselves into a corner by defining sets of alignments within a slice as forming template alignments from which we can compute TLEN rather than storing it. We attempted to codify this, but it is very clear that different implementations codify it in different ways. Yes, we're back to TLEN#1 and TLEN#2 interpretations between htsjdk and htslib for CRAM decoding.

So I think I need to be brutal and modify the CRAM spec to be equally vague and arbitrary, but perhaps with a footnote to give recommendations to encoders to side-step the issue by using "detached" reads whenever there are scenarios that may lead to differing interpretations. (And be equally vague as to what those scenarios are as codifying that is also a nightmare)

[colinhercus]

I think you should just make TLEN obsolete

[jkbonfield]

It's hard to do that, but I'd be in favour of essentially saying it's aligner specific (with recommendations on meaning) as that's essentially where we're at and have always been at.

For CRAM it's trickier as it's formally defined in the specification. Furthermore, if we wish to get the compression benefits of not duplicating data between read-pairs, then defining it has a purpose. It just needs some nuance for the contentious cases and maybe recommending anything that's not "obvious" is stored verbatim so we get most of the compression gains while avoiding the myriad of weird corner cases. Hopefully...

Edit: a demonstration of using mate-flags to encode NovaSeq read-pair fields via a defined algorithm vs using "detached" throughout and storing verbatim (standard CRAM 3.1 output from htslib):

-rw-r--r-- 1 jkb team117 175754187 Sep  8 09:29 /tmp/_paired.cram
-rw-r--r-- 1 jkb team117 204027635 Sep  8 09:38 /tmp/_detached.cram

Ie not exploiting this data redundancy during encode would cost us 16% storage for this NovaSeq sample. The details of the specific CRAM data series are CF (cram flags), MF (flate flags), NS (RNEXT equiv), TS (template size), NF (next fragment, how many records later). Fully detached we have a very big TS field as expected and MF increases too. With pairing data we vastly reduce TS/MF but add some back in fields like NF, NP, and CF.

The flags,RNEXT,PNEXT,TLEN fields in the two files consume:

paired:    7.4Mb
detached: 34.4Mb

So that's a 78% reduction in storage for those fields.

Attempting to always use read-pairing data even on the contentious cases (as htslib right now stores verbatim rnext/pnext/tlen when the meaning of tlen is ambiguous) only saves a marginal amount more, as most (70%) of the existing detached reads come from read-pairs spanning CRAM slices.

Obviously it depends on chemistry, library, instrument, etc too.