Transcriptome preset (-x tx), 2.3-2.8x faster than minimap2 for quantification

CHANGELOG.md

@@ -4,6 +4,22 @@ Each release's binaries and full notes are on the
 [releases page](https://github.com/DylanCosto/bonefish/releases). "Output" means the alignment
 records. Each entry quotes the MD5 of the genome-wide HG001 sample's SAM output, `@PG` lines removed.
 
+## v0.3.0 (unreleased)
+
+**New: `-x tx`.** The preset maps direct RNA or cDNA to a transcriptome for quantifiers such as
+oarfish and NanoCount (minimap2 `-ax map-ont -N 100`).
+
+- **Speed:** on HG002 it is 2.3x (cDNA) and 2.8x (direct RNA) faster than minimap2.
+- **Accuracy:** it maps the same reads as minimap2. oarfish gene counts agree at r = 0.9995 or
+  higher. Report: `docs/validation/2026-10-21-transcriptome/`.
+
+**Faster with many reference sequences:** the ALT-contig count is computed once per index, not
+by scanning every target for every read. This matters for transcriptomes and metagenome
+assemblies.
+
+**Faster chaining:** isolated anchors are dropped before chaining. DNA uses 1.2% less CPU, with
+output unchanged (`9bc7ff3e…`).
+
 ## v0.2.0 (2026-10-21)
 
 **New: direct RNA and cDNA.** Spliced alignment to a genome is available through two presets:

README.md

@@ -6,7 +6,8 @@ RNA and cDNA reads to a reference genome. It writes SAM, PAF or a sorted, indexe
 - **Accurate.** On whole human genomes (GIAB HG002 and HG005), SNP, indel and structural-variant
   calls match minimap2's. Methylation calls agree at r = 0.999.
 - **Fast.**
-  - Aligner alone: 2.6-8x faster than minimap2 on DNA, 2.6x on direct RNA and 5.4x on cDNA.
+  - Aligner alone: 2.6-8x faster than minimap2 on DNA, 2.6x on direct RNA, 5.4x on cDNA and
+    2.3-2.8x for transcriptome quantification.
   - Reads to a sorted BAM: about 3x faster than `minimap2 | samtools sort`.
 - **Drop-in.** It reads dorado's BAM output, keeping methylation and other tags, and accepts
   minimap2's options and preset names.
@@ -53,6 +54,7 @@ in a few minutes.
 | `dna` (default) | genomic DNA: R10.4.1 or R9.4.1, sup or hac, ultra-long (minimap2 `lr:hq`) |
 | `rna` | direct RNA (minimap2 `splice -uf -k14`) |
 | `cdna` | cDNA (minimap2 `splice`) |
+| `tx` | RNA or cDNA to a transcriptome, for quantification with oarfish or NanoCount (minimap2 `map-ont -N 100`) |
 
 Input formats, output options, memory use and CPU support are covered in
 [docs/usage.md](docs/usage.md). To replace minimap2 in an existing pipeline, including ONT's
@@ -81,9 +83,11 @@ Sniffles2 against GIAB v5.0q.
 | DNA, HG002 R9.4.1 chr20 (4.6 Gbp) | 412 s | 103 s | 4.0x |
 | cDNA, HG002 PCS114 (809k reads) | 1,223 s | 227 s | 5.4x |
 | direct RNA, HG002 RNA004 (902k reads) | 456 s | 177 s | 2.6x |
+| direct RNA to the transcriptome (`-x tx`, 200k reads) | 106 s | 37 s | 2.8x |
 
 On RNA, 99.5-99.9% of primary alignments are identical to minimap2's, and junction accuracy
-against GENCODE is the same.
+against GENCODE is the same. Transcript quantification from `-x tx` alignments agrees with
+minimap2's at gene level at r = 0.9995-0.99997.
 
 bonefish is also faster than dorado's aligner (4.3x), mm2plus (5.1x), mm2-fast (14x) and
 Winnowmap2 (38x).

bonefish-core/src/align/index.rs

@@ -135,6 +135,10 @@ pub struct Index {
     /// computed on first use for the ava self-skip fast path. See [`name_rank`].
     #[serde(skip)]
     name_rank: std::sync::OnceLock<NameRankData>,
+    /// Derived (not serialized): number of ALT contigs, counted on first use (after `--alt` has
+    /// marked them), so per-read code does not scan every target (500k+ for a transcriptome).
+    #[serde(skip)]
+    n_alt: std::sync::OnceLock<usize>,
 }
 
 /// Dense lexicographic ranking of target-sequence names, used to turn the
@@ -211,6 +215,7 @@ impl IndexBuilder {
                 backend: LookupBackend::BucketHash(super::index_bucket::BucketHashLookup::empty()),
                 packed_seqs: Default::default(),
                 name_rank: Default::default(),
+                n_alt: Default::default(),
             },
             w, k, max_occ,
             bucket_bits,
@@ -304,6 +309,11 @@ impl IndexBuilder {
 }
 
 impl Index {
+    /// Number of ALT contigs (computed once).
+    pub fn n_alt(&self) -> usize {
+        *self.n_alt.get_or_init(|| self.seqs.iter().filter(|s| s.is_alt).count())
+    }
+
     /// Dense lexicographic ranking of the target names, computed once on first use and cached
     pub fn name_rank(&self) -> &NameRankData {
         self.name_rank.get_or_init(|| {
@@ -538,6 +548,7 @@ impl Index {
                 backend: LookupBackend::BucketHash(lookup),
                 packed_seqs: Storage::mapped(map.clone(), qo, qn),
                 name_rank: Default::default(),
+                n_alt: Default::default(),
             });
             pos += hdr.part_len;
             pos += pad_to(pos, RAW_ALIGN);
@@ -707,6 +718,7 @@ impl Index {
             backend,
             packed_seqs: packed_seqs.into(),
             name_rank: Default::default(),
+            n_alt: Default::default(),
         };
 
         eprintln!("[*] Index loaded: {} seqs, {}M bases in {:.1}s",
@@ -727,6 +739,7 @@ impl Index {
             backend: LookupBackend::BucketHash(super::index_bucket::BucketHashLookup::empty()),
             packed_seqs: Default::default(),
             name_rank: Default::default(),
+            n_alt: Default::default(),
         }
     }
 
@@ -740,6 +753,7 @@ impl Index {
             backend: LookupBackend::BucketHash(super::index_bucket::BucketHashLookup::empty()),
             packed_seqs: Default::default(),
             name_rank: Default::default(),
+            n_alt: Default::default(),
         }
     }
 
@@ -769,6 +783,7 @@ impl Index {
             backend: LookupBackend::BucketHash(super::index_bucket::BucketHashLookup::empty()),
             packed_seqs: Default::default(),
             name_rank: Default::default(),
+            n_alt: Default::default(),
         };
 
         let mut offset = 0usize;

bonefish-core/src/align/map.rs

@@ -652,6 +652,36 @@ fn get_minimizer_index(qlen: usize, a: &Minimizer, mini_pos: &[u64]) -> Option<u
 }
 
 /// Estimate chain divergence from minimizer matching
+/// Remove anchors (sorted by `x`) with no other anchor of the same reference and strand within
+/// `max_dist` reference bases. Such an anchor is outside every other anchor's chaining window,
+/// so it can only form a one-anchor chain, which `min_cnt` (> 1) always drops: the chains are
+/// the same without it. A transcriptome yields many such scattered hits (a k-mer of the read
+/// on hundreds of short transcripts).
+fn drop_isolated_anchors(a: &mut Vec<Minimizer>, max_dist: i32, min_cnt: i32) {
+    let n = a.len();
+    if min_cnt <= 1 || n < 2 {
+        return;
+    }
+    let md = max_dist.max(0) as u64;
+    let near = |p: &Minimizer, q: &Minimizer| p.x >> 32 == q.x >> 32 && q.x.wrapping_sub(p.x) <= md;
+    let mut w = 0;
+    for i in 0..n {
+        if (i > 0 && near(&a[i - 1], &a[i])) || (i + 1 < n && near(&a[i], &a[i + 1])) {
+            a[w] = a[i];
+            w += 1;
+        }
+    }
+    a.truncate(w);
+}
+
+/// Copy each region's chain anchors (`anchor_count` from offset `original_as` of `chains`).
+fn fill_anchors(regs: &mut [Mapping], chains: &[Minimizer]) {
+    for r in regs.iter_mut() {
+        let (a, n) = (r.original_as, r.anchor_count as usize);
+        r.anchors = chains[a..a + n].to_vec();
+    }
+}
+
 fn estimate_divergence(regs: &mut [Mapping], _anchors: &[Minimizer], qlen: usize, mini_pos: &[u64], mi: &Index) {
     if mini_pos.is_empty() { return; }
 
@@ -839,6 +869,13 @@ pub fn map_query(
     // presets), otherwise DP
     let is_splice = opt.flags.contains(AlignFlags::SPLICE);
     let mut anchors = std::mem::take(&mut ctx.anchors);
+    // The rescue's density test below counts every anchor, isolated ones included.
+    let n_anchors = anchors.len();
+    if !is_splice && !opt.flags.intersects(AlignFlags::RMQ_CHAIN | AlignFlags::PAR_CHAIN) {
+        // Widest reference distance any chaining pass of this read uses (main, the -g re-chain).
+        let reach = opt.chaining.max_gap.max(max_chain_gap_ref).max(opt.chaining.bandwidth);
+        drop_isolated_anchors(&mut anchors, reach, opt.chaining.min_cnt);
+    }
     let (mut u, mut chains) = if is_splice {
         chain_anchors_cdna(&opt.chaining, max_chain_gap_ref, max_chain_gap_qry, &mut anchors, &mut ctx.chain_bufs)
     } else if opt.flags.contains(AlignFlags::RMQ_CHAIN) {
@@ -889,7 +926,7 @@ pub fn map_query(
             let needs_rescue = ((qlen.saturating_sub(coverage) > opt.chaining.rmq_rescue_size as usize)
                 || (coverage > (qlen as f32 * opt.chaining.rmq_rescue_ratio) as usize))
                 && crate::align::env_flags::EXP_RESCUE_COVER.is_none_or(|f| (coverage as f32) < f * qlen as f32)
-                && (opt.chaining.rescue_max_density == 0 || ctx.anchors.len() < opt.chaining.rescue_max_density as usize * qlen);
+                && (opt.chaining.rescue_max_density == 0 || n_anchors < opt.chaining.rescue_max_density as usize * qlen);
             let needs_rescue = needs_rescue && (!(opt.chaining.rescue_join || *crate::align::env_flags::EXP_RESCUE_JOIN) || {
                 let max_dist = opt.chaining.max_gap.max(opt.chaining.bandwidth_long);
                 let j = chains_joinable(&u, &chains, max_dist, opt.chaining.bandwidth_long);
@@ -952,12 +989,12 @@ pub fn map_query(
 
         if cnt == 0 { continue; }
 
-        // Collect chain anchors
-        let chain_anchors_vec = chains[k_offset..k_offset + cnt as usize].to_vec();
-
+        // The chain's anchors are copied into the region only if it survives the parent /
+        // secondary filtering below (`fill_anchors`): with many overlapping chains per read
+        // (repeats, transcript isoforms) most regions are dropped there.
         let is_qstrand = opt.flags.contains(AlignFlags::QSTRAND);
         let ChainCoords { rid, rev, rs, re, qs, qe, hash: chain_hash } =
-            chain_coords(&chain_anchors_vec, qseq.len(), is_qstrand, read_hash);
+            chain_coords(&chains[k_offset..k_offset + cnt as usize], qseq.len(), is_qstrand, read_hash);
 
         regs.push(Mapping {
             ref_id: rid,
@@ -968,7 +1005,7 @@ pub fn map_query(
             score: sc,
             initial_chain_score: sc,  // initially same as score
             anchor_count: cnt,
-            anchors: chain_anchors_vec,
+            anchors: Vec::new(),
             is_reverse: rev,
             s2_score: None,
             is_secondary: false,
@@ -1004,7 +1041,7 @@ pub fn map_query(
     }
 
     // Mark ALT-contig chains; then sort chains by descending score.
-    let n_alt = mi.seqs.iter().filter(|s| s.is_alt).count();
+    let n_alt = mi.n_alt();
     if n_alt > 0 {
         for r in regs.iter_mut() {
             if mi.seqs[r.ref_id].is_alt { r.is_alt = true; }
@@ -1027,6 +1064,7 @@ pub fn map_query(
 
     // Skip parent-based filtering when ALL_CHAINS is set (e.g. ava-ont, ava-pb)
     if opt.flags.contains(AlignFlags::ALL_CHAINS) {
+        fill_anchors(&mut regs, &chains);
         // Nearby-seed computation on post-filter chains
         // Build squeezed anchor array from surviving chains, sorted by position
         let mut pos_order: Vec<usize> = (0..regs.len()).collect();
@@ -1219,6 +1257,7 @@ pub fn map_query(
         compacted.push(r);
     }
     regs = compacted;
+    fill_anchors(&mut regs, &chains);
     let mut filtered_regs = regs;
     let pc3 = Instant::now();
     PC_PARENT_NS.fetch_add((pc3 - pc2).as_nanos() as u64, Relaxed);

bonefish-core/src/align/pipeline.rs

@@ -1769,7 +1769,7 @@ fn process_query_core(
     let squeezed = map_result.squeezed;
     let split_mode = out.split_mode;
     // Mark ALT contigs on chains before sorting
-    let n_alt = mi.seqs.iter().filter(|s| s.is_alt).count();
+    let n_alt = mi.n_alt();
     if n_alt > 0 {
         for r in regs.iter_mut() {
             if mi.seqs[r.ref_id].is_alt { r.is_alt = true; }

bonefish-core/src/api.rs

@@ -43,6 +43,8 @@ pub enum Preset {
     Cdna,
     /// Accurate cDNA, spliced (minimap2 `-x splice:hq`).
     SpliceHq,
+    /// RNA or cDNA reads to a transcriptome, for quantification (minimap2 `-x map-ont -N 100`).
+    Tx,
 }
 
 impl Preset {
@@ -53,6 +55,7 @@ impl Preset {
             Preset::Rna => "rna",
             Preset::Cdna => "cdna",
             Preset::SpliceHq => "splice:hq",
+            Preset::Tx => "tx",
         }
     }
 }
@@ -568,7 +571,7 @@ fn parse_cigar_string(s: &str) -> Vec<CigarOp> {
 // ---------------------------------------------------------------------------
 
 /// Apply a preset by string name ("dna" or one of its minimap2-name aliases, "rna", "cdna",
-/// "splice" (= cdna) or "splice:hq").
+/// "splice" (= cdna), "splice:hq" or "tx").
 ///
 /// This is the same function used by the CLI's `-x` flag. Modifies `opt`,
 /// `k` and `w` in place. Returns an error if the preset name
@@ -613,6 +616,18 @@ pub fn apply_preset_str(opt: &mut MapOptions, k: &mut usize, w: &mut usize, pres
             opt.alignment.gap_subanchor = true;
             opt.alignment.gap_subanchor_tr = 150;
         },
+        // Reads (direct RNA or cDNA) against a transcriptome, for quantifiers that take all
+        // alignments (oarfish, NanoCount, salmon): minimap2's map-ont seeding (k15 w10, no seed
+        // occurrence range, 5 kb query gap limit) for short and noisy RNA reads, up to 100
+        // secondary alignments per read (isoforms share most of their sequence), and otherwise
+        // the dna preset's cheaper chaining and DP (there are no introns to align).
+        "tx" => {
+            apply_preset_str(opt, k, w, "dna")?;
+            *k = 15; *w = 10;
+            opt.chaining.max_gap = 5000;
+            opt.seeding.min_mid_occ = 10; opt.seeding.max_mid_occ = 1_000_000;
+            opt.filtering.best_n = 100;
+        },
         // Spliced alignment of transcripts, as minimap2 2.30/2.31 `-x splice` (options.c): k15 w5,
         // both transcript strands tried (-ub) with splice-site flank scoring, introns up to 200 kb
         // (-G: chaining reference gap and bandwidth), cDNA-mode chaining, 2 kb query gap limit,
@@ -1088,7 +1103,7 @@ mod tests {
 
     #[test]
     fn test_preset_as_str_roundtrip() {
-        let presets = [Preset::Dna, Preset::Rna, Preset::Cdna, Preset::SpliceHq];
+        let presets = [Preset::Dna, Preset::Rna, Preset::Cdna, Preset::SpliceHq, Preset::Tx];
         for p in &presets {
             let s = p.as_str();
             assert!(!s.is_empty(), "preset {:?} should have a string form", p);

bonefish/src/main.rs

@@ -52,11 +52,13 @@ struct AlignArgs {
     ///   rna              — direct RNA vs a reference genome (spliced; minimap2 -x splice -uf -k14)
     ///   cdna             — cDNA vs a reference genome (spliced; minimap2 -x splice); alias: splice
     ///   splice:hq        — accurate cDNA (minimap2 -x splice:hq)
+    ///   tx               — RNA/cDNA vs a transcriptome, for quantification (minimap2 -x map-ont -N 100)
     /// dna = minimap2 lr:hq seeding and scoring plus --ext-bw 200 --ext-adaptive 64 --ext-rescue
     /// --ext-probe 2000,1000 --gap-bw-margin 500 --gap-bw-narrow 32 --ungapped-fill --max-gap-ref 5000
     /// --seg-rmq --clean-dp --rescue-join --rescue-max-density 2 --gap-subanchor-tr 150
     /// cdna = -k15 -w5 -g2k -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200,100 -ub --junc-bonus 9
     /// --junc-pen 5 --splice-flank yes; rna = cdna with -k14 -uf
+    /// tx = dna with -k15 -w10 -g5000 -U10,1000000 -N100
     #[arg(short = 'x', long, verbatim_doc_comment)]
     pub preset: Option<String>,
 

bonefish/tests/data/expected.tsv

@@ -14,3 +14,4 @@ cdna_paf_cs	bonefish	-x cdna -c --cs	21afca13cdff0b35
 cdna_junc_bed	bonefish	-x cdna -a --junc-bed {junc}	6292d532a2fcabc3
 cdna_splice_hq	bonefish	-x splice:hq -a	1988a08a9ab7353a
 cdna_alias_splice	bonefish	-x splice -a	60d6d8c7c31cd299
+tx_sam	bonefish	-x tx -a	d53bfcfb58fa2138

bonefish/tests/regression.rs

@@ -55,7 +55,8 @@ fn expected() -> Vec<(String, String, String, String)> {
 }
 
 /// Arguments for an expected.tsv case: its options (`{junc}` stands for the annotation BED), then
-/// the reference and reads. Cases named `rna_*` / `cdna_*` use the spliced dataset (made by
+/// the reference and reads (`tx_*` runs the transcriptome preset on the DNA dataset, for its
+/// seeding and secondary settings). Cases named `rna_*` / `cdna_*` use the spliced dataset (made by
 /// `scripts/make_regression_rna_data.py`): direct-RNA-like or cDNA-like reads on a genome of
 /// multi-exon genes.
 fn map_args(name: &str, case_args: &str) -> Vec<String> {

docs/benchmarks.md

@@ -96,6 +96,20 @@ the reference near the chain's anchors with a 10-mer search for unanchored exons
 aligned once for both transcript strands. Details:
 [validation/2026-10-21-rna-cdna](validation/2026-10-21-rna-cdna/README.md).
 
+## Transcriptome quantification
+
+`-x tx` against minimap2 `-ax map-ont -N 100` on GENCODE v49 transcripts (HG002, 200k reads each):
+
+| | minimap2 | bonefish `-x tx` |
+|---|---:|---:|
+| cDNA, wall | 185 s | 81 s |
+| direct RNA, wall | 106 s | 37 s |
+| reads mapped | | the same as minimap2 |
+| oarfish gene counts vs minimap2's: Pearson r | | 0.9995 (cDNA), 0.99997 (direct RNA) |
+
+On simulated reads, the primary lands on the true transcript at least as often as with minimap2.
+Details: [validation/2026-10-21-transcriptome](validation/2026-10-21-transcriptome/README.md).
+
 ## Metagenomes, duplex reads and rDNA
 
 Against minimap2 `lr:hq` on the same reads (v0.1.7): ONT's Zymo fecal reference (1.70 M R10.4.1 sup

docs/pipelines.md

@@ -16,6 +16,7 @@ minimap2 in a pipeline in two ways:
 | `minimap2 -ax splice -uf -k14 ref.fa drna.fq` | `bonefish -x rna -a ref.fa drna.fq` |
 | `minimap2 -ax splice ref.fa cdna.fq` | `bonefish -x cdna -a ref.fa cdna.fq` |
 | `minimap2 -ax splice --junc-bed anno.bed ...` | `bonefish -x cdna -a --junc-bed anno.bed ...` |
+| `minimap2 -ax map-ont -N 100 transcripts.fa reads.fq` (for quantification) | `bonefish -x tx -a transcripts.fa reads.fq` |
 
 Notes:
 
@@ -75,6 +76,22 @@ minimap2 does.
 reads ([report](validation/2026-10-21-rna-cdna/README.md)), but not yet against these downstream
 tools.
 
+## Transcript quantification
+
+oarfish, NanoCount and similar quantifiers take reads aligned to the transcriptome with many
+secondary alignments. `-x tx` writes up to 100 secondaries per read, as `minimap2 -N 100` does:
+
+```bash
+bonefish -t 16 -x tx -a --eqx transcripts.fa reads.fq | samtools view -b -o tx.bam -
+```
+
+- **Agreement with minimap2:** on HG002 direct RNA and cDNA, oarfish gene counts from these
+  alignments agree with counts from minimap2's at r = 0.9995-0.99997. Simulated reads land on
+  their true transcript at least as often as with minimap2
+  ([report](validation/2026-10-21-transcriptome/README.md)).
+- **oarfish needs a workaround for now:** it only accepts BAMs whose `@PG` header names an
+  aligner it has validated, and bonefish is not yet on its list.
+
 ## Containers
 
 ```bash

docs/usage.md

@@ -8,6 +8,7 @@
 | `rna` | direct RNA (RNA004), spliced to a genome | `-x splice -uf -k14` |
 | `cdna` | cDNA, spliced to a genome, transcript strand from splice signals | `-x splice` |
 | `splice:hq` | accurate cDNA | `-x splice:hq` |
+| `tx` | direct RNA or cDNA to a transcriptome, for quantification (oarfish, NanoCount) | `-x map-ont -N 100` |
 
 The minimap2 names `lr:hq`, `lr:hq-fast`, `map-ont` and `lr` are accepted for `dna`, and `splice`
 for `cdna`. PacBio, short-read and assembly presets are not supported.
@@ -21,6 +22,7 @@ build an index once:
 bonefish -d ref.raw --idx-raw ref.fa              # dna; loads in under a second (8.5 GB for GRCh38)
 bonefish -x cdna -d cdna.raw --idx-raw ref.fa     # spliced presets need their own index (k/w differ)
 bonefish -x rna  -d rna.raw  --idx-raw ref.fa
+bonefish -x tx   -d tx.raw   --idx-raw transcripts.fa
 ```
 
 Without `--idx-raw`, `-d` writes a smaller index that takes a few seconds to load. Indexes are