Construction from slice of text

[RagnarGrootKoerkamp]

Hey,

Looking throught all the different fm-index crates on crates.io, AWRY looks like one of the more promising ones.

Unfortunately, you only seem to support construction from a FASTA file. It'd be super helpful if you could also support construction directly from a &[u8].

Could you also write a bit more about exactly which alphabets are then supported? Should input be ACTG, ACTGactg, or is alphabet size detected automatically and would any ASCII work?

Oh, and a brief comparison in the readme with all the other existing FM-index crates/implementations going around would also be helpful. At least I know that the one in rust-bio is missing some of the optimizations needed for low memory requirements, and I think some of the other implementations going around have similar inefficiencies, but it's hard to know as an outsider to the field of state-of-the-art FM index implementations :)

[Sawwave]

Hi Ragnar,

I've asked the maintainers of the suffix array generation library I'm using to implement this as a feature on their end. I could emulate this by simply writing the given text slice as a fasta and read it in on creation, but that seems a bit hackish and it'd be nice to not have to go to disk.

Assuming they get this done in a timely manner, I expect this to be a quick fix, but in the meantime I'm waiting to see when this will be implemented.

As for the alphabets, the current supported alphabets are DNA, RNA, and Amino Acid sequences, and are case insensitive. I'll add this to the documentation.

I'll consider adding comparisons against the other crates, but as those crates may be updated I may expect that the rankings could become out of date. Perhaps this comparison is a contribution you might be suited to add?

[RagnarGrootKoerkamp]

Hey, thanks for your replies :)

Ah so it looks like libsufr does not actually do linear time construction? Then I'll also ask the authors of that crate for a quick benchmark, since otherwise it's really hard to say anything about performance.

But anyway, the size of the suffix array is easy to get by just multiplying n by either 4 or 8 depending on whether it uses 32 of 64bit integers, right?

How hard would it be to support generic 8bit alphabets as well? For my current project I somewhat need that as well. (I could work around it by splitting 8-bit values over 4 2-bit bases, but that'd be somewhat annoying.)

[Sawwave]

We use sufr to generate the suffix array, but AWRY stores the data in a bespoke manner that allows for efficient computation of the occurrence function, so we aren't limited by the sufr SA storage methodology, just the runtime/mem for building the SA. I've already implemented the MemSize trait as you suggested, it should be included in the next version release (version 0.2.0).

I've given a lot of though to generic alphabets, and there's enough complication that would make for some distinct design challenges. Essentially, the issue with generic 8-bit alphabets is that at the beginning of every block there needs to be a "milestone" count of the number of occurrences of each character in the alphabet up to that point. For Nucleotide data, that's 6 values- A,G,C,T,N, and the sentinel $.
For a generic 8-bit alphabet that means that the milestone data for each block would need to be 256 64-bit values, or 2KiB of data for each 256 positions in the text.

My current thoughts on the matter are that generic text indexing would probably end up going something like this:
Break each ascii character into 2 nybbles, and encode each nybble as a 5-bit value that includes an extra indicator bit to denote if it's from the upper nybble or lower nybble. That way, each milestone array would only contain 33 64-bit values, (2^5 + the $ sentinel). The extra bit is necessary, otherwise the count() function will return extra values that map to invalid reading frames (i.e., a query matches something that starts halfway through a 8-bit element, and just happens to match the query).

[Sawwave]

The other potential option for generic 8-bit data is to deal with the larger milestone arrays by creating larger blocks in the BWT. AWRY has blocks of 256 values to map easily to 256-bit SIMD vectors, but a block could contain multiple 256-bit SIMD vectors. This would mean that the act of counting the character occurrences between the beginning of a block and the actual query position may require additional work to bitwise-smush the vectors together and use popcount operations, but these actions are currently a tiny runtime component compared to the comparatively enormous hit of the near-guaranteed cache miss of reading in the block.

[RagnarGrootKoerkamp]

Ok, it sounds like for now the easiest is then to simply encode my 8-bit data into DNA. Then the only limitation is that the FM-index may return matches that don't start at a multiple of 4 (aligned with the original data) but I can filter those out manually.

I'll consider adding comparisons against the other crates, but as those crates may be updated I may expect that the rankings could become out of date. Perhaps this comparison is a contribution you might be suited to add?

At the moment (for the past 3 months, really) I'm just looking for an off-the-shelf FM-index that's easily usable in Rust. While I could do an FM-index comparison, I'm really not the right person for that, never having written one myself.

Some implementations that I've come across (half encouraging you to compare them ;), half note to self -- maybe I turn this into a quick blogpost):

• sdsl lite
• FM-Index by mpetri, who recommends this over sdsl lite. But I don't have found/created rust bindings yet.
• https://github.com/TravisWheelerLab/AvxWindowFmIndex, the C++ this is based on.

Then there's a bunch on crates.io:

• https://crates.io/crates/fm-index
• https://crates.io/crates/awry (you :) )
• https://crates.io/crates/nucleic-acid
• https://crates.io/crates/qwt: quad-wavelet-tree only, not a full FM-index.
• https://crates.io/crates/lt-fm-index
• https://docs.rs/bio/2.0.3/bio/data_structures/fmindex/trait.FMIndexable.html

Really, the only way to choose between them is with some statistics. Say for human genome DNA input, and maybe also some generic text input (although most probably don't support that):

• construction time
• max construction space
• datastructure space
• negative query speed (random queries of some fixed length of say 10bp and 1000bp)
• positive query speed (substring queries of same fixed length)

Also, it would be good to have a little table writing for each the theoretical expected/worst-case construction time (whichever is more relevant) and theoretical space usage, and what the highlight of each method is.

[Sawwave]

I'll consider adding comparisons against the other crates, but as those crates may be updated I may expect that the rankings could become out of date. Perhaps this comparison is a contribution you might be suited to add?

I was just suggesting a benchmark comparison of the various FM-index crates, no implementation knowledge necessary :)

No worries if that's not something you're interested in, it's just not something I will have time for in the near future, but seems like something you might be able to do to contribute to the community.

Ok, it sounds like for now the easiest is then to simply encode my 8-bit data into DNA.

This should work, although it might end up slower than other fm-index implementations that are designed for full u8 compatibility. Just keep in mind that if some of the locations need to be filtered out, this will mean that the count() function will be unreliable - it'll return all the actual instances of your query, but will also contain the count of things you'd need to filter out as well.

[RagnarGrootKoerkamp]

I was just suggesting a benchmark comparison of the various FM-index crates, no implementation knowledge necessary :)

No worries if that's not something you're interested in, it's just not something I will have time for in the near future, but seems like something you might be able to do to contribute to the community.

(Sorry, I'm tired at the moment and have been bothered by this for a while now, and as you can see none of the other crates does this, so it's not specific to you, but)
I really think it's on library authors to sell and compare their work. Just adding some new crate without a 'why' is really just adding noise to crates.io and making it harder for people to choose and for the community to settle on one optimized implementation. I could compare things myself, but that'll only be relevant as long as no new crates pop up.

You write that your port a fastest-in-class FM-index, but until you show the numbers for this port specifically, it's just impossible as someone not in-the-know to know what that means. Especially for highly performant code, rewriting something using slightly different code (and in a different language also) could probably make up to 2x speedup or slowdown, so just referring to the numbers for the C++ implementation doesn't cut it for me.

I'm sure you have some benchmarks comparing AWRY to the original C++ implementation. Just showing those (and ideally that performance is close or better) already gives me much more confidence to go with AWRY :)

[Sawwave]

I'm sure you have some benchmarks comparing AWRY to the original C++ implementation. Just showing those (and ideally that performance is close or better) already gives me much more confidence to go with AWRY :)

This project is designed to be publicly usable, but it's a side-project in addition to my current job for our research group and affiliated groups. I will post benchmarking when I get it, but it will likely not be anywhere near the holistic comparison to the rest of the rust ecosystem that you would like to see.

The original version of this library is (if I can say so myself) incredibly well optimized, and the current benchmarks that I've run suggest that the AWRY is somewhere around 50% slower than the C implementation. This roughly agrees with most performance comparisons I've seen, where Rust code is fast, but can't compete with well written C or C++.

I'll post benchmarks when I get them, but I'm currently working with the writers of the sufr suffix array generation library to deal with some performance issues on their side when the input has large runs of the same character (as you get when analyzing human chromosomes, for instance). When they've gotten performance figured out on their end, we can get an accurate view of the overall performance of the system.

[RagnarGrootKoerkamp]

AWRY is somewhere around 50% slower than the C implementation.

Allright, that sounds reasonable. For final applications I suppose that might be enough of a difference to switch to the C++ implementation with FFI, but for me for a quick comparison that's good enough for sure.

Rust code is fast, but can't compete with well written C or C++.

I disagree ;) You just have to work slightly harder on it, especially in making sure to use unchecked indexing in hot loops (that's also what they do in that post and they get within 1% which may just be noise). In rust this feels slightly iffy, but keep in mind that array indexing operation in C++ is unsafe by default, so you really shouldn't feel bad about adding a bit of unsafe every once in a while to improve performance.
See also eg this recent blog I wrote:
https://curiouscoding.nl/posts/static-search-tree/

I'm currently working with the writers of the sufr suffix array generation library to deal with some performance issues on their side

Thanks!

[RagnarGrootKoerkamp]

(Just a headsup that as I'm watching the repo, I see all the version releases coming by. I'd say there isn't much need to tag new versions all the time. Once a week is plenty for sure. Either way I/others can also just use the git version here if the latest is needed. That should also keep the commit history a bit cleaner with just the actual commits/PRs, rather than a lot of version bumps and release tags.)