haystackfm

A GPU-accelerated FM-index for DNA sequence alignment. It runs on native targets (Vulkan / Metal / DX12) through wgpu, and compiles to WebAssembly for in-browser WebGPU use — the same index code, on the desktop or in a tab.

haystackfm supports the full 16-symbol IUPAC ambiguity alphabet (A C G T N R Y S W K M B D H V) in both its CPU and GPU query paths, with a pluggable Alphabet trait for swapping in custom matching semantics.

🔎 Try it in your browser: the live demo builds an index and runs count / locate queries entirely client-side via WebGPU — no server, no upload.

What is an FM-index?

An FM-index is a compressed full-text substring index built on the Burrows–Wheeler Transform (BWT). It answers two core questions over a text of length n without keeping the text uncompressed:

• count(P) — how many times pattern P occurs, in O(|P|) time.
• locate(P) — where each occurrence is, in O(|P| + occ·k) time.

That makes it the workhorse behind short-read aligners (BWA, Bowtie) and other bioinformatics tooling. haystackfm pushes the construction and the query hot loops onto the GPU while keeping a CPU implementation as the correctness ground truth.

Feature overview

Where to go next

• New here? Start with Installation and the Quick Start.
• Want the mental model? Read Concepts and Architecture.
• Looking for the type-level API? It lives on docs.rs.

Installation

haystackfm is a Rust crate published on crates.io. Its minimum supported Rust version (MSRV) is 1.87.

CPU only (default)

cargo add haystackfm

The default cpu feature builds and queries indices entirely on the CPU — no GPU driver or wgpu dependency required.

With GPU acceleration

cargo add haystackfm --features gpu

The gpu feature pulls in wgpu and enables GPU-accelerated construction plus the batched GPU query paths. On native targets it uses Vulkan, Metal, or DX12 depending on your platform.

For the browser (WebAssembly)

cargo add haystackfm --features wasm

Build the WASM package with wasm-pack:

wasm-pack build --target web --features wasm

This emits a pkg/ directory with haystackfm.js + haystackfm_bg.wasm and TypeScript type definitions. See WebAssembly for the JS/TS API.

Feature flags at a glance

See Feature Flags for details.

Quick Start

The fastest path: build an index on the CPU and run a couple of queries.

Build and query (CPU)

use haystackfm::{DnaSequence, FmIndex, FmIndexConfig};

let seq = DnaSequence::from_str("ACGTACGT")?;
let config = FmIndexConfig::default();
let index = FmIndex::build_cpu(&[seq], &config)?;

// A DnaSequence holds IUPAC-encoded bytes; encode the query the same way.
let query = DnaSequence::from_str("ACGT")?;

assert_eq!(index.count(query.as_slice()), 2);
let positions = index.locate(query.as_slice()); // Vec<(SeqId, u32)> — (seq id, offset)

count returns how many times the pattern occurs; locate returns the (SeqId, offset) of each occurrence. A SeqId is the reference’s 0-based position in build order, not its FASTA name — resolve one with index.seq_header(id), and go the other way with index.seq_id("chr1"). See Sequence ids vs. headers for why, and Count & Locate for the full surface.

Build on the GPU

GPU construction is async and returns the same FmIndex type:

use haystackfm::{DnaSequence, FmIndex, FmIndexConfig};
use haystackfm::gpu::locate::locate_batch_gpu;

let seqs = vec![DnaSequence::from_str("ACGTACGT")?];
let config = FmIndexConfig::default();

let index = FmIndex::build(&seqs, &config).await?;   // async GPU build

// Batched, IUPAC-aware GPU locate
let ctx = haystackfm::gpu::GpuContext::new().await?;
let query = DnaSequence::from_str("ACGT")?;
let queries: Vec<&[u8]> = vec![query.as_slice()];
let hits = locate_batch_gpu(&ctx, &index, &queries).await?;
// hits[i] = Vec<(seq_id, offset_within_seq)>

Requires the gpu feature. See GPU Acceleration.

Prefer the browser?

The live demo does all of this — build (CPU or GPU), search, serialize — without leaving the page. It’s the quickest way to see haystackfm run.

Configuration

FmIndexConfig::default() is a sensible starting point. The knobs that matter most:

let config = FmIndexConfig {
    sa_sample_rate: 4,
    ..Default::default()
};

Concepts

A quick tour of the moving parts, so the rest of the guide has names to hang on.

The building blocks

An FM-index is assembled from a few classic components:

• Suffix Array (SA) — the sorted order of all suffixes of the text. It’s the backbone that turns “where does P occur” into a contiguous range.
• Burrows–Wheeler Transform (BWT) — a reversible permutation of the text derived from the SA. It clusters similar contexts together, which is what makes the index compressible and searchable.
• C array — for each symbol, the number of text characters that sort strictly before it. A tiny histogram (16 values for the IUPAC alphabet).
• Occ table — rank support: Occ(c, i) = how many times symbol c appears in the BWT up to position i. This is the hot data structure every query step touches.

Together, the LF-mapping (C[c] + Occ(c, i)) lets a backward search walk the text one symbol at a time, shrinking an SA interval until it exactly covers every occurrence of the pattern.

Backward search

count and locate both start with a backward search: process the pattern right-to-left, and at each step apply LF-mapping to narrow the [lo, hi) SA interval. When the interval is empty the pattern is absent; otherwise its size is the occurrence count.

Because haystackfm is IUPAC-aware, a single query symbol can match several BWT symbols, so the search tracks a union of intervals rather than one — see IUPAC Ambiguity.

Locate

count only needs the final interval size. locate must turn each SA-interval position into a text coordinate. Storing the whole SA would defeat compression, so haystackfm samples it every k positions (sa_sample_rate) and LF-walks from an unsampled position until it lands on a sample — trading a little query time for a much smaller index.

Sequence ids vs. headers

An index is built over several reference sequences concatenated into one text, so every result has to say which reference it came from. haystackfm keeps two distinct things apart here:

• A header is the FASTA name of a reference — "chr1", "NC_045512.2". It is what a human reads, and it is a String.
• A SeqId is that reference’s position in build order — the first sequence you pass to the builder is SeqId(0). It is a u32.

Queries only ever report SeqId. Locating an occurrence yields (SeqId, offset), never a header. This matters because locate cost is per occurrence: a seed conserved across hundreds of references produces hundreds of hits, and cloning a header string for each one — only for the caller to look it up again — is pure overhead. Returning the integer lets a caller index a Vec directly.

Translating between the two is a separate, O(1) step:

index.seq_header(id)   // Option<&str>   — id  -> header
index.seq_id(header)   // Option<SeqId>  — header -> id
index.seq_headers()    // &[String]      — every header, indexed by SeqId::index()

The usual pattern is to build whatever label table you need once at load time from seq_headers(), then index it by SeqId::index() per hit.

Ids are stable: they are assigned in build order and preserved across to_bytes / from_bytes, so an id recorded against one load of an index still means the same reference after a reload. They are not portable across differently-built indexes — reordering or adding sequences renumbers everything.

Because seq_id must be an unambiguous inverse, headers must be unique: a collision fails the build with FmIndexError::DuplicateHeader. Sequences supplied without a header are auto-named seq_{i} in build order.

Bidirectional index

Maximal Exact Match (MEM) and Super-Maximal Exact Match (SMEM) finding need to extend a match in both directions. That requires a bidirectional FM-index — BidirFmIndex — which maintains synchronized intervals over the forward and reverse texts. See MEM / SMEM Finding.

Where the GPU comes in

Construction (SA, BWT, Occ) and the query hot loops are all data-parallel, which maps well to GPU compute. haystackfm implements each stage as WGSL compute shaders and keeps the CPU implementation as the correctness oracle — every GPU path is parity-tested against it. The full pipeline is laid out in Architecture.

Count & Locate

The two core exact-match queries. Both operate on IUPAC-encoded byte slices (DnaSequence::as_slice()).

count

let n = index.count(query.as_slice()); // u32 — number of occurrences

Runs a backward search and returns the final SA-interval size — O(|P|), independent of how many times the pattern occurs.

locate

let hits = index.locate(query.as_slice()); // Vec<(SeqId, u32)>

Returns every occurrence as (seq_id, offset_within_seq), where seq_id is a SeqId — the reference’s 0-based position in build order, not its FASTA name — and offset is a u32 within that reference. Cost is O(|P| + occ·k), where occ is the number of hits and k is sa_sample_rate — each hit LF-walks at most k steps to the nearest SA sample.

No header string is allocated per hit. Since locate cost scales with occurrence count, that matters most exactly where it hurts: a seed conserved across many references.

Sequence ids and headers

SeqId and the FASTA header are separate things; convert between them explicitly. Both directions are O(1).

index.seq_header(id)?;         // Option<&str>   — id -> header
index.seq_id("chr1");          // Option<SeqId>  — header -> id
index.seq_headers();           // &[String]      — all headers, indexed by SeqId::index()

To label results, resolve per hit:

for (id, offset) in index.locate(query.as_slice()) {
    println!("{}:{offset}", index.seq_header(id).unwrap());
}

Or, in a hot loop, build your own table once and index it by SeqId::index():

// Done once at load time, not per occurrence.
let labels: Vec<MyLabel> = index.seq_headers().iter().map(MyLabel::parse).collect();

for (id, offset) in index.locate(query.as_slice()) {
    let label = &labels[id.index()];   // plain Vec indexing — no hashing
}

Headers must be unique: a duplicate fails the build with FmIndexError::DuplicateHeader, which is what makes seq_id an exact inverse of seq_header. Sequences built without an explicit header are named seq_{i} in build order.

Ids survive serialization — see Serialization.

Index metadata

index.text_len();        // u32 — total length of the concatenated text
index.num_sequences();   // u32 — number of indexed sequences

GPU batch locate

When you have many patterns, resolve them together on the GPU. It is IUPAC-aware and returns results per query:

use haystackfm::gpu::locate::locate_batch_gpu;

let ctx = haystackfm::gpu::GpuContext::new().await?;
let queries: Vec<&[u8]> = vec![q1.as_slice(), q2.as_slice()];
let hits = locate_batch_gpu(&ctx, &index, &queries).await?;
// hits[i] = Vec<(seq_id, offset_within_seq)> for queries[i]

Requires the gpu feature. GPU results carry the same occurrences as the CPU locate output (order is not guaranteed), reported with the same (SeqId, offset) shape — so the same label-resolution code works for both. FmIndex::locate_gpu wraps this for a whole batch. For when the GPU path is worthwhile, see GPU Acceleration.

MEM / SMEM Finding

Exact-match queries answer “is this whole pattern present?”. Seeding a read aligner instead needs the longest exact stretches shared between a query and the reference — Maximal Exact Matches (MEMs) and Super-Maximal Exact Matches (SMEMs). These require the bidirectional index, BidirFmIndex.

Build a bidirectional index

use haystackfm::{DnaSequence, BidirFmIndex, FmIndexConfig};

let refs = vec![DnaSequence::from_str("ACGTACGTACGT")?];
let config = FmIndexConfig::default();
let bidir = BidirFmIndex::build_cpu(&refs, &config)?;

CPU MEM / SMEM

let query = DnaSequence::from_str("ACGT")?;

// Vec<Mem>. `min_len` filters short matches; `locate` resolves positions.
let smems = bidir.find_smems(query.as_slice(), /*min_len=*/18, /*locate=*/true);
let mems  = bidir.find_mems(query.as_slice(),  /*min_len=*/18, /*locate=*/true);

Both are IUPAC-aware (N matches any of A/C/G/T). Passing locate = false skips position resolution and leaves Mem::positions empty — cheaper when you only need match extents and counts.

Mem

pub struct Mem {
    pub query_start: usize,              // 0-based inclusive
    pub query_end:   usize,              // 0-based exclusive
    pub match_count: u32,                // number of occurrences
    pub positions:   Vec<(SeqId, u32)>,  // (seq id, offset) — empty when locate=false
}

positions identifies each reference by SeqId — its 0-based build order, not its FASTA name. Resolve one with bidir.seq_header(id), or build a label table once from bidir.seq_headers() and index it by id.index(). This is the hot path the id representation exists for: a conserved seed can occur in hundreds of references, and positions are resolved per occurrence.

GPU MEM / SMEM

For batches of queries, run the GPU pipeline. Queries are passed as &[DnaSequence], and reference boundaries come from the index:

#[cfg(feature = "gpu")]
async fn run(bidir: &haystackfm::BidirFmIndex, query: haystackfm::DnaSequence) -> Result<(), Box<dyn std::error::Error>> {
let boundaries = bidir.seq_boundaries();     // reference-sequence boundaries
let queries = [query.clone()];               // &[DnaSequence]

// Vec<Vec<MemHit>> — the GPU context is drawn from a process-wide cache.
let smem_hits = bidir.find_smems_gpu(&queries, /*min_len=*/18, boundaries, /*max_hits_per_mem=*/1024).await?;
let mem_hits  = bidir.find_mems_gpu(&queries,  18, boundaries, 1024).await?;
// smem_hits[query_i] = Vec<MemHit> with resolved reference positions
Ok(()) }

MemHit

pub struct MemHit {                     // GPU result type
    pub query_start: u32,
    pub query_end:   u32,
    pub match_count: u32,
    pub positions:   Vec<(SeqId, u32)>, // same shape as Mem::positions
    pub truncated:   bool,              // true if positions capped at max_hits_per_mem
}

positions has the same (SeqId, offset) shape as the CPU Mem::positions, so the same label-resolution code works against either path.

max_hits_per_mem caps how many positions each MEM resolves; when a MEM has more occurrences than the cap, truncated is set. GPU results are parity-tested against the CPU find_mems / find_smems output.

IUPAC Ambiguity

Queries and references may contain any of the 16 IUPAC nucleotide symbols. Two symbols match when their base sets share at least one nucleotide.

For example, a query N matches any base; a query R matches A or G (and any ambiguity code whose base set includes one of them).

CPU/GPU parity

Both paths use the same compatibility lookup — compatible_symbols on the CPU and the COMPAT table in WGSL on the GPU. Parity tests enforce that the two stay in sync, so a query returns the same matches whether it runs on the CPU or the GPU.

Opting out

If you want strict A/C/G/T matching where ambiguity codes never match (useful for peer-comparable benchmarks), build with the ExactDna alphabet instead of the default. See Custom Alphabets.

Custom Alphabets

Matching semantics are pluggable through the Alphabet trait (src/alphabet.rs). Rather than carrying a generic type parameter, FmIndex and BidirFmIndex store a runtime AlphabetFns bundle (function pointers plus a serialization tag), so the index type itself stays alphabet-agnostic while the match rules are chosen at build time.

Built-in alphabets

Choosing one

The default build_cpu / build use IupacDna. To pick a different alphabet, use the _with constructors:

use haystackfm::alphabet::ExactDna;
use haystackfm::{FmIndex, BidirFmIndex};

let index = FmIndex::build_cpu_with::<ExactDna>(&seqs, &config)?;
let bidir = BidirFmIndex::build_cpu_with::<ExactDna>(&seqs, &config)?;

Implementing your own

Implement Alphabet for a custom type to define your own symbol set and match rules. The trait carries a safety contract — stable function pointers and a unique serialization tag (≥ 128) so serialized indices can be matched back to their alphabet. See the trait docs in src/alphabet.rs on docs.rs for the exact requirements.

GPU Acceleration

With the gpu feature, construction and the query hot loops offload to the GPU via wgpu. On native targets that means Vulkan, Metal, or DX12; in the browser it means WebGPU (see WebAssembly).

Building on the GPU

let index = FmIndex::build(&seqs, &config).await?; // async; returns a normal FmIndex

The suffix array, BWT, and Occ table are all built with WGSL compute shaders. The result is the same FmIndex you’d get from build_cpu — you can query it on either the CPU or the GPU.

GpuContext caching

Adapter and device initialization is not free, so haystackfm caches a process-wide GpuContext in an OnceLock (src/gpu/context_cache.rs). GPU functions accept an optional pre-initialized context: pass None on the first call and reuse the cached context afterward. This is what keeps benchmarks and tests from paying init overhead on every call.

GPU query paths

• Batch locate — locate_batch_gpu runs backward search and SA resolution on the GPU for a batch of patterns. See Count & Locate.
• MEM / SMEM — find_mems_gpu / find_smems_gpu run the bidirectional-extension pipeline on the GPU. See MEM / SMEM Finding.

Every GPU path is validated against the CPU implementation, which remains the correctness oracle. The multi-pass shader pipelines are described in Architecture.

When the GPU helps

The GPU wins on batches of queries over large indices, where there’s enough data-parallel work to amortize dispatch and transfer costs. For a single short pattern on a small index, the CPU path is typically faster — the query itself is cheaper than moving data to the device. Measure with your own workload; the benchmarks in the repo (cargo bench --features gpu) are a good starting template.

WebAssembly

The same index runs in the browser. With the wasm feature, haystackfm exposes a wasm-bindgen JS/TS API and uses the browser’s WebGPU for GPU construction and queries.

Building the package

cargo add haystackfm --features wasm
wasm-pack build --target web --features wasm

This produces a pkg/ directory containing haystackfm.js, haystackfm_bg.wasm, and .d.ts type definitions.

JS/TS API

import init, { FmIndexBuilder, FmIndexHandle } from "./pkg/haystackfm.js";

await init();

const builder = new FmIndexBuilder(/*sa_sample_rate=*/32);
builder.add_fasta(`>seq1\nACGTACGT\n>seq2\nTGCATGCA`);

const handle = await builder.build_gpu();    // or builder.build_cpu()

console.log(handle.count("ACGT"));           // number of occurrences
console.log(handle.locate("ACGT"));          // Array of [seqId, offset] pairs

// seqId is the sequence's 0-based build order, not its FASTA header.
console.log(handle.seq_header(0));           // "seq1"  — id -> header, O(1)
console.log(handle.seq_id("seq1"));          // 0       — header -> id, O(1)
console.log(handle.seq_headers());           // ["seq1", "seq2"] — indexed by seqId

console.log(handle.text_len());              // total text length
console.log(handle.num_sequences());         // number of indexed sequences

const bytes = handle.to_bytes();             // Uint8Array — serialize
const restored = FmIndexHandle.from_bytes(bytes);

Requirements

GPU builds and queries require a browser with WebGPU. CPU builds (build_cpu) work anywhere WASM runs. See Browser Support, and try the live demo to see it end-to-end.

Serialization

A built index can be serialized to bytes and reloaded, so you can construct once and cache the result rather than rebuilding on every run.

Rust

let bytes = index.to_bytes()?;              // Vec<u8>
let restored = FmIndex::from_bytes(&bytes)?; // FmIndex

BidirFmIndex exposes the same pair. The serialized form records the alphabet’s serialization tag, so a reloaded index keeps the matching semantics it was built with (see Custom Alphabets).

Sequence ids are stable

The headers are stored in build order, so a SeqId means the same reference before and after a round trip — you can record ids against a serialized index and reuse them on reload. (The header -> id map is derived from the header list rather than stored, so from_bytes rebuilds it; this is invisible to callers.)

Ids are only meaningful within one index. A differently-built index — sequences added, removed, or reordered — renumbers them; match on headers when crossing that boundary.

JavaScript / TypeScript

const bytes = handle.to_bytes();             // Uint8Array
const restored = FmIndexHandle.from_bytes(bytes);

Trusting input

Deserialization reconstructs an index from a byte buffer. Treat from_bytes input as you would any deserialized data — only load payloads you produced or otherwise trust. See the project security policy for the threat model.

Architecture

How the pieces fit together, from construction through the GPU query pipelines.

Construction pipeline

Input: [DnaSequence]
  → concatenate + add sentinels
  → GPU: build suffix array          radix_sort.wgsl + sa_*.wgsl
  → GPU: derive BWT                  bwt_gather.wgsl
  → CPU: compute C array             (trivial histogram, 16 values)
  → GPU: build Occ table             occ_scan.wgsl + prefix_sum.wgsl
  → sample SA at rate k

CPU query pipeline

O(m) count, O(m + occ·k) locate:

backward_search (IUPAC multi-interval) → SA interval set [lo, hi)
locate: LF-walk from each position to nearest SA sample, fused symbol+rank
        lookup per step (OccTable::lf_step — one block-plane read instead of two)

GPU locate pipeline (2 passes)

Pass 1  locate_search.wgsl   — backward search, IUPAC multi-interval → (match_count, intervals)
Pass 2  locate_resolve.wgsl  — LF-walk to sampled SA position → (seq_id, offset)

GPU MEM / SMEM pipeline (3 passes)

Pass 1  mem_find.wgsl        — bidirectional extension, IUPAC multi-interval → raw SA intervals
Pass 2  mem_resolve.wgsl     — resolve each SA position via LF-walk
Pass 3  ref_map.wgsl         — binary-search reference boundaries → (ref_id, offset)

Key modules

Compact Occ table

No resident Bwt is kept after construction. The OccTable compacts to the effective symbol alphabet (bitplane-encoded lanes, not the full 16-symbol IUPAC space) and reconstructs BWT rows/bytes on demand for GPU upload. This keeps the index small while still feeding the fixed-lane layout the shaders expect.

Feature Flags

cpu

The default. Everything builds and queries on the CPU with no GPU dependency. This is the correctness ground truth for every GPU path.

gpu

Adds wgpu and the GPU construction + query paths (FmIndex::build, locate_batch_gpu, find_mems_gpu / find_smems_gpu). Required for all GPU benchmarks and tests:

cargo test --features gpu
cargo bench --features gpu

wasm

Adds the wasm-bindgen bindings for browser use and implies gpu (browser WebGPU). Build the package with wasm-pack build --target web --features wasm. See WebAssembly.

Browser Support

GPU construction and queries in the browser require WebGPU.

CPU builds (build_cpu) work in any browser that runs WebAssembly — WebGPU is only needed for the GPU build path (build_gpu).

Checking at runtime

if ("gpu" in navigator) {
  // WebGPU available — build_gpu will work
} else {
  // Fall back to build_cpu
}

The live demo does exactly this: it detects WebGPU, shows a badge, and enables or disables the GPU build button accordingly.

API Reference

The complete, always-current type-level API reference is generated from the source by rustdoc and hosted on docs.rs:

👉 docs.rs/haystackfm

There you’ll find every public type, trait, and function with its signatures and doc comments, including:

• FmIndex — build_cpu, build_cpu_with, build, count, locate, to_bytes, from_bytes, the id/header accessors seq_headers, seq_header, seq_id, and the base accessors sequence, sequence_by_header.
• BidirFmIndex — build_cpu, build_cpu_with, find_mems, find_smems, the GPU variants, and the same id/header and base accessors.
• SeqId — a reference’s stable 0-based id, reported by every query in place of its FASTA header. See Sequence ids vs. headers.
• Mem / MemHit — result types for MEM/SMEM finding.
• FmIndexConfig — construction knobs.
• alphabet — the Alphabet trait, IupacDna, ExactDna, DnaSequence, and compatible_symbols.
• gpu — GpuContext, locate_batch_gpu, and the GPU MEM/SMEM functions (behind the gpu feature).

This guide covers the how and why; docs.rs is the exhaustive what.

Live Demo

haystackfm runs entirely in your browser. The demo builds an FM-index (on the CPU or the GPU via WebGPU), runs count / locate queries, and serializes the index — all client-side, with nothing uploaded to a server.

What you can do

1. Input sequences — paste plain ACGT lines or FASTA, or load the bundled example.
2. Build — construct the index on the CPU or, if your browser supports WebGPU, on the GPU.
3. Search — enter a pattern and see the occurrence count and located positions.
4. Serialize — export the built index to bytes and reload it, exercising to_bytes / from_bytes.

Requirements

The GPU build path needs a browser with WebGPU (Chrome/Edge 113+, Safari 18+, Firefox behind a flag — see Browser Support). The CPU build works in any browser that runs WebAssembly; the demo detects WebGPU and enables the GPU button only when it’s available.

Contributing

Contributions are welcome. The authoritative guides live in the repository:

• CONTRIBUTING.md — build / test / lint workflow and PR conventions.
• CODE_OF_CONDUCT.md — community expectations.
• SECURITY.md — how to report a security issue (please don’t open a public issue for vulnerabilities).

Local checks

Before opening a PR, make sure these pass:

cargo fmt --check
cargo clippy --all-features -- -D warnings
cargo test --features gpu

GPU tests and benchmarks require the gpu feature and a working GPU/driver.

License

haystackfm is licensed under the Apache License, Version 2.0. Unless you state otherwise, any contribution you submit for inclusion is licensed as above, without additional terms or conditions.