how to index a read in a circular sequence?

[yfarjoun]

I've been (re-)reviewing a PR to HTSJDK that enables TP:circular in @SQ and that took me back to the unfinished business of how to index and random-access circular contigs.

If we extended the index into the greater-than-LEN space, we could easily index the reads that "hang off" the end of the contig. Then, in order to query them, one should also add LEN and see if there are any reads at query_pos+LEN (when querying for query_pos). Since the assumption is that we only have linear alignments that have POS≤LEN, we only need to look at query_pos+LEN and not, e.g., query_pos+2*LEN.

Thoughts?

[jkbonfield]

Follow up from our file formats monthly catch up:

We determined two potential strategies, with the second being preferable.

1. In BAI/CSI the index maps genomic coordinates to file offsets (ie read start points), so we could potentially add a read at that extends beyond the end of the chromosome in two index bins. However it's not going to work in CRAM.

2. For circular genomes, we may just want to query two locations. If a user asks for Z:1-10 and chr Z is circular with LN 1,000,000, then we can do an index lookup for Z:1 and Z:1,000,000. Where this gets tricky is if the user queries Z:100,000. Do we still look up pos Z:1,000,000? If we have no prior knowledge of the maximum read length then yes! However potentially a smart implementation could look at read-group platforms and know what's possible and what's not, or implementations could also offer command-line tunables such as the maximum read length to permit for overlaps.

Neither system is perfect, but 2 has fewer problems.

There is also potentially an option for indexing to store additional bins that are beyond the end of the specified length (BAI/CSI). CRAM may infact already do this when writing out the container headers which specify the extents that things map to however I suspect it's specified to clip to the stated length (as buggy aligners sometimes do produce such data in error and it's a problem when attempting to do reference based compression if your reads overhang the reference), but I'm not 100% sure there is anything explicit to forbid the true alignment end point being listed in the index. It needs research.

The original PR #405 has a lot of discussion on this topic, but contained the comical get-out clause of "The impact of this representation on indexing and random access is yet to be explored by implementations". I guess the chickens are coming home to roost.