While we're here, I notice that samtools-sort.1 (line 78) says that -M produces output that cannot be indexed. I haven't tested it, but I don't see why this would be the case…?

While we're here, I notice that samtools-sort.1 (line 78) says that -M produces output that cannot be indexed. I haven't tested it, but I don't see why this would be the case…?

I think it's wrong and probably written forgetting that it can do partially mapped files still.

Rather curiously, --write-index does fail with an inappropriate "Ignoring --write-index as it only works for position sorted files" message. A subsequent samtools index works just fine. That's also a buglet of samtools sort.

It's not something we can determine up-front whether it is entirely unmapped, but actually that's not a problem anyway as we can still build an index and still query it via things like idxstats so it has some benefit. Technically we could even still use the index for random access and sub-setting of the data, albeit non-trivially. So there's no reason to forbid indexing of unmapped data.

(As an aside, bizarrely samtools import accepts --write-index but doesn't make an index.)

Also sadly, the decision of making minhash either coordinate or unsorted is making the index statement invalid. Fortunately our index doesn't check as it validates on the alignments rather than trusting the header (but it could and maybe some do). However the man page for sort explicitly states we must use coordinate sort if we want an index, which isn't true as we can build an index on the unsorted data (because technically it's still coordinate sorted, as noted above).

I'll try picard. Once I find a version of java that works.

Edit: as I expected. Samtools sort -M on develop produces a bam file which picard can index, while this PR produces a file that it cannot index as it's not SO:coordinate.

This is the problem really with adjusting things to look "nice". We now have a decision to make. Do we want this data to be obvious to users, or do we want to retain the (admittedly minor) ability of being able to query index stats on an unaligned file? I think probably it's OK, as samtools import claims unsorted too (correct) but also falls under the "technically sorted" definition of there being no chromosomes that are out of order. So it's a really minor thing, but something that needs a concious decision rather than an accidental change of compatibility.

Edit: ah there's another reason for indexing unaligned data. In CRAM it directly translates to the ability to sub-set using samtools cat, as CRAM's indexing lends itself neatly to fractional extraction of unaligned data simply by taking the Nth-Mth containers. However the fact that htsjdk can no longer index the minhash data is irrelevant for CRAM, as it cannot index normal aligned CRAMs either! (Or at least if it can it's not exposed via picard)

Updated as a separate commit for ease of review. I'll squash later (it won't make it into this upcoming release, sorry).

Also this talk of indexing reminds me of the original inspiration which came from Illumina's Tony Cox.

His code was computing a minimiser and actually having that minimiser as a coordinate on a fake chromosome (which sadly limits the size due to BAM indexing capabilities, unless we have many fake chromosomes). An index then is obviously more sensible. The benefit of this is then you can provide a query sequence, produce minimisers from that query, and then do random access to pull out sequences that share that minimiser. It's a crude way of finding data that may potentially align to a query sequence without first assembling or aligning it all. The hit rate wasn't that great though as you can only place a sequence in one location, so you can't have proper sequence sketches consisting of multiple minimisers. That also limits it to short read data.

A cool idea though. :) I just took the locality => better compression component of it and ran with that.