RNA golden tests, pipeline guide, changelog, contributing guide

.github/ISSUE_TEMPLATE/bug_report.yml

@@ -0,0 +1,43 @@
+name: Bug report
+description: Wrong output, a crash, or an error
+labels: [bug]
+body:
+  - type: textarea
+    id: what
+    attributes:
+      label: What happened
+      description: What you expected, and what you got instead (error message, wrong alignment, crash).
+    validations:
+      required: true
+  - type: input
+    id: version
+    attributes:
+      label: bonefish version
+      description: Output of `bonefish --version`, and how you installed it (release binary, source, Docker, conda).
+    validations:
+      required: true
+  - type: textarea
+    id: command
+    attributes:
+      label: Command line
+      render: shell
+    validations:
+      required: true
+  - type: input
+    id: system
+    attributes:
+      label: System
+      description: OS and CPU (e.g. Ubuntu 22.04, AMD EPYC 7763 / Apple M2), memory, threads used.
+    validations:
+      required: true
+  - type: textarea
+    id: data
+    attributes:
+      label: Data
+      description: Read type (DNA / direct RNA / cDNA, chemistry, basecaller), reference, and if possible a few reads that reproduce the problem (attach or link).
+  - type: textarea
+    id: log
+    attributes:
+      label: Log output
+      description: bonefish's stderr (add `-v` for stage timers).
+      render: text

.github/ISSUE_TEMPLATE/config.yml

@@ -0,0 +1,5 @@
+blank_issues_enabled: true
+contact_links:
+  - name: Usage documentation
+    url: https://github.com/DylanCosto/bonefish/blob/main/docs/usage.md
+    about: Presets, input and output options, memory

.github/ISSUE_TEMPLATE/feature_request.yml

@@ -0,0 +1,16 @@
+name: Feature request
+description: A new option, preset, output or workflow
+labels: [enhancement]
+body:
+  - type: textarea
+    id: need
+    attributes:
+      label: What would you like bonefish to do?
+      description: The use case, and how you do it today (e.g. with minimap2 options).
+    validations:
+      required: true
+  - type: textarea
+    id: data
+    attributes:
+      label: Data
+      description: Read type and a public dataset that shows the need, if any.

CHANGELOG.md

@@ -0,0 +1,90 @@
+# Changelog
+
+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.2.0 (2026-10-21)
+
+**New: direct RNA and cDNA.** Spliced alignment to a genome is available through two presets:
+
+- `-x rna`: direct RNA, as minimap2 `-ax splice -uf -k14`.
+- `-x cdna`: cDNA, as minimap2 `-ax splice`. `splice:hq` is also available.
+
+Spliced runs accept `--junc-bed` and minimap2's splice options, and add the `ts:A` transcript-strand
+tag.
+
+**Speed:** against minimap2 on HG002, direct RNA runs 2.6x faster and cDNA 5.4x.
+
+- Introns are aligned against only the reference near the chain's anchors, with a 10-mer search
+  for exons that have no anchor.
+- cDNA chains are aligned once instead of once per strand.
+- Chaining uses the SIMD chainers, and read-end extensions use reduced windows.
+
+**Accuracy:** 99.5-99.9% of primary alignments are identical to minimap2's, and GENCODE junction
+accuracy is unchanged. Report: `docs/validation/2026-10-21-rna-cdna/`.
+
+**DNA:** output unchanged (`9bc7ff3e…`).
+
+**Removed:** the unused GPU experiment and the WebAssembly browser demo.
+
+## v0.1.7 (2026-10-21)
+
+**AVX-512 kernels** for the wide and narrow DP, chaining, small sorts and minimizer sketching,
+chosen at run time.
+
+- On a Xeon Gold 5418Y, the genome-wide sample takes 27% less CPU time than with v0.1.6.
+- On AVX2 machines it takes 3% less.
+- Output unchanged (`9bc7ff3e…`).
+
+## v0.1.6 (2026-10-16)
+
+**Fix:** the BAI/CSI index written with `-o x.bam` had about 2.5x too many chunks, so region readers
+(Clair3, longphase, IGV) were slower on bonefish BAMs.
+
+- The index is now identical in content to `samtools index`.
+- Clair3 on HG002 chr20 runs in 181 s instead of 216 s.
+- Output unchanged.
+
+## v0.1.5 (2026-10-16)
+
+**Sub-anchored gap filling is back**, outside tandem repeats, recovering most of v0.1.4's cost.
+
+- Whole-genome re-validation on HG002 and HG005: SNP and indel F1 equal to minimap2's.
+- The tandem-repeat benchmark matches minimap2.
+- Output `9bc7ff3e9f498cc1530b2b1d78c47e10`.
+
+## v0.1.4 (2026-10-15)
+
+**Fix:** fewer false SV calls in tandem repeats.
+
+- Sub-anchored gap filling is off: it could align reads through the wrong repeat copy.
+- Costs about 6% CPU.
+- Output `8fbcd8ce9df67b78a978ef6497873caa`.
+
+## v0.1.3 (2026-10-15)
+
+**Fixes:**
+
+- Reads from a pipe now work.
+- A CSI index is written for references longer than 512 Mbp.
+- 7 crashes on corrupt BAM/BGZF input are fixed.
+
+**Changes:** the default batch is `-K 100M`, using 2.7 GB less memory. Output unchanged.
+
+## v0.1.2 (2026-10-14)
+
+**Input and CLI hardening:**
+
+- Malformed input is fatal, with a clear error.
+- `-o x.bam` implies sorting.
+- `-y` copies dorado BAM tags and `@RG` lines.
+
+## v0.1.1 (2026-10-14)
+
+**Platforms and validation:**
+
+- Faster arm64 chaining.
+- The DP kernel is ported to NEON and SSE4.1, so arm64 output equals x86.
+- First whole-genome validation (HG002).
+- Output `deed7949d35f15bd8eed935be47a1d81`.

CONTRIBUTING.md

@@ -0,0 +1,50 @@
+# Contributing
+
+Bug reports, test data and pull requests are welcome. For bugs, please use the issue template:
+the version, command line, CPU and a few reads that reproduce the problem make it fast to fix.
+
+## Build and test
+
+```bash
+cargo build --release --features jemalloc
+cargo test --release --features jemalloc       # unit tests, golden-output regression, doc tests
+cargo clippy --release --features jemalloc --all-targets
+```
+
+CI runs the same on Linux x86_64, Linux arm64 and macOS (arm64). On x86 runners with AVX-512,
+CI requires the AVX-512 tests to run rather than skip.
+
+## Output changes
+
+The golden-output test (`bonefish/tests/regression.rs`) pins the exact output of every preset on
+small synthetic datasets in `bonefish/tests/data/`. Every change therefore falls into one of two
+cases:
+
+- **Speed or refactoring:** the output must stay byte-identical, and the regression test must pass
+  unchanged. For larger changes, also compare a real dataset before and after (SAM with `@PG`
+  lines removed). The validation reports quote MD5s for the genome-wide HG001 sample.
+- **Intended output changes:** regenerate the expected hashes and explain the change in the pull
+  request:
+
+  ```bash
+  BONEFISH_BLESS=1 cargo test --release --features jemalloc --test regression -- --nocapture
+  ```
+
+  Paste the printed rows into `bonefish/tests/data/expected.tsv`. An output change also needs
+  evidence that accuracy does not get worse. For the `dna` preset, that means variant calling
+  against GIAB (Clair3 for small variants, Sniffles2 for SVs) on at least one chromosome. For the
+  spliced presets, it means junction accuracy against an annotation and against simulated reads
+  with known intron chains. `docs/validation/` shows how earlier changes were validated, with
+  scripts.
+
+The test datasets are generated deterministically by `scripts/make_regression_data.py` (DNA) and
+`scripts/make_regression_rna_data.py` (RNA/cDNA).
+
+## Code
+
+- **Style:** match the surrounding code.
+- **SIMD:** kernels need a scalar or portable fallback, selected at run time. Tests that compare
+  them with the reference implementation should skip through
+  `crate::align::dp::simd_test_available` when the CPU lacks the feature.
+- **DNA speed:** new presets or options must not slow the `dna` preset. Keep per-read and per-call
+  branches out of the DP and chaining inner loops.

README.md

@@ -55,7 +55,8 @@ in a few minutes.
 | `cdna` | cDNA (minimap2 `splice`) |
 
 Input formats, output options, memory use and CPU support are covered in
-[docs/usage.md](docs/usage.md).
+[docs/usage.md](docs/usage.md). To replace minimap2 in an existing pipeline, including ONT's
+wf-human-variation, see [docs/pipelines.md](docs/pipelines.md).
 
 ## Speed and accuracy
 
@@ -127,6 +128,11 @@ The rest comes from engineering:
 `bonefish-core` (imported as `bonefish`) exposes the aligner as a Rust library: `bonefish::Aligner`
 with `Preset::Dna`, `Preset::Rna` and `Preset::Cdna`. See `bonefish-core/examples/`.
 
+## Contributing
+
+Bug reports and pull requests are welcome; see [CONTRIBUTING.md](CONTRIBUTING.md). Changes in each
+release are listed in [CHANGELOG.md](CHANGELOG.md).
+
 ## Credits and license
 
 bonefish is a fork of [rammap](https://github.com/jwanglab/rammap), a pure-Rust reimplementation

bonefish/tests/data/expected.tsv

@@ -6,3 +6,11 @@ dna_paf_cs	bonefish	-c --cs	0dc5b802637a8b0f
 dna_md_eqx	bonefish	-a --MD --eqx	0d347ab709e3136b
 alias_sam	bonefish	-ax lr:hq-fast	81ec695369877c8c
 alias_map_ont	bonefish	-ax map-ont	81ec695369877c8c
+rna_sam	bonefish	-x rna -a	3ae4701b00bc65f8
+rna_paf_cs	bonefish	-x rna -c --cs	f779a1c786e11fa3
+rna_md_eqx	bonefish	-x rna -a --MD --eqx	f603af3054848bb4
+cdna_sam	bonefish	-x cdna -a	60d6d8c7c31cd299
+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

bonefish/tests/data/rna_junc.bed

@@ -0,0 +1,28 @@
+chrA	3000	27808	g000.0	0	+	3000	27808	0	6	274,96,100,201,210,82,	0,1165,11029,12943,24298,24726,
+chrA	3000	27808	g000.1	0	+	3000	27808	0	5	274,100,201,210,82,	0,11029,12943,24298,24726,
+chrA	35456	47240	g001.0	0	-	35456	47240	0	3	80,247,188,	0,8523,11596,
+chrA	49884	90342	g002.0	0	-	49884	90342	0	5	109,274,151,124,64,	0,15496,15883,27828,40394,
+chrA	93260	100977	g003.0	0	+	93260	100977	0	4	154,62,283,124,	0,535,938,7593,
+chrA	93260	100977	g003.1	0	+	93260	100977	0	3	154,283,124,	0,938,7593,
+chrA	107965	111167	g004.0	0	+	107965	111167	0	3	101,107,248,	0,1516,2954,
+chrA	115619	116229	g005.0	0	-	115619	116229	0	2	65,219,	0,391,
+chrA	118487	134147	g006.0	0	-	118487	134147	0	3	243,69,217,	0,562,15443,
+chrA	141963	158302	g007.0	0	+	141963	158302	0	7	167,171,137,179,134,134,102,	0,11618,11893,12276,12683,14844,16237,
+chrA	161228	164150	g008.0	0	+	161228	164150	0	3	235,167,196,	0,594,2726,
+chrA	166758	193519	g009.0	0	-	166758	193519	0	6	162,150,149,32,287,131,	0,300,9758,12447,24328,26630,
+chrA	166758	193519	g009.1	0	-	166758	193519	0	5	162,149,32,287,131,	0,9758,12447,24328,26630,
+chrA	198909	199392	g010.0	0	+	198909	199392	0	2	77,253,	0,230,
+chrA	203051	267166	g011.0	0	-	203051	267166	0	7	189,200,51,226,77,239,297,	0,18218,33887,53819,61928,63385,63818,
+chrB	3000	8456	g100.0	0	+	3000	8456	0	3	179,290,96,	0,3123,5360,
+chrB	16257	43864	g101.0	0	-	16257	43864	0	6	147,81,166,23,155,292,	0,504,8828,9206,9844,27315,
+chrB	47448	48127	g102.0	0	-	47448	48127	0	2	260,60,	0,619,
+chrB	55992	58423	g103.0	0	-	55992	58423	0	6	71,120,212,210,294,124,	0,242,1070,1435,1742,2307,
+chrB	55992	58423	g103.1	0	-	55992	58423	0	5	71,212,210,294,124,	0,1070,1435,1742,2307,
+chrB	60544	85044	g104.0	0	-	60544	85044	0	6	178,297,58,283,225,283,	0,368,868,12247,23111,24217,
+chrB	88220	108914	g105.0	0	-	88220	108914	0	5	50,177,22,109,194,	0,6535,19686,20052,20500,
+chrB	116030	116737	g106.0	0	+	116030	116737	0	2	245,107,	0,600,
+chrB	123988	125057	g107.0	0	-	123988	125057	0	3	241,253,156,	0,425,913,
+chrB	132503	136179	g108.0	0	+	132503	136179	0	3	275,169,135,	0,1699,3541,
+chrB	139218	141033	g109.0	0	-	139218	141033	0	2	93,122,	0,1693,
+chrB	143084	143360	g110.0	0	+	143084	143360	0	2	75,97,	0,179,
+chrB	145583	169642	g111.0	0	+	145583	169642	0	7	123,280,65,205,277,208,297,	0,3585,4476,22208,22766,23421,23762,

bonefish/tests/regression.rs

@@ -54,11 +54,23 @@ fn expected() -> Vec<(String, String, String, String)> {
         .collect()
 }
 
-fn map_args(case_args: &str) -> Vec<String> {
-    let mut a: Vec<String> = case_args.split_whitespace().map(String::from).collect();
+/// 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
+/// `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> {
+    let junc = data("rna_junc.bed").to_str().unwrap().to_string();
+    let mut a: Vec<String> = case_args.split_whitespace().map(|x| x.replace("{junc}", &junc)).collect();
     a.extend(["-t", "2"].map(String::from));
-    a.push(data("ref.fa.gz").to_str().unwrap().into());
-    a.push(data("reads.fq.gz").to_str().unwrap().into());
+    let (r, q) = if name.starts_with("rna_") {
+        ("rna_ref.fa.gz", "rna_reads.fq.gz")
+    } else if name.starts_with("cdna_") {
+        ("rna_ref.fa.gz", "cdna_reads.fq.gz")
+    } else {
+        ("ref.fa.gz", "reads.fq.gz")
+    };
+    a.push(data(r).to_str().unwrap().into());
+    a.push(data(q).to_str().unwrap().into());
     a
 }
 
@@ -67,7 +79,7 @@ fn outputs_match_expected() {
     let bless = std::env::var_os("BONEFISH_BLESS").is_some();
     let mut failures = Vec::new();
     for (name, source, args, hash) in expected() {
-        let a = map_args(&args);
+        let a = map_args(&name, &args);
         let out = run(&a.iter().map(String::as_str).collect::<Vec<_>>());
         let got = fnv1a(&out);
         if bless && source == "bonefish" {

docs/pipelines.md

@@ -0,0 +1,85 @@
+# Using bonefish in pipelines
+
+bonefish takes the same inputs as minimap2 and writes standard SAM or BAM, so it can replace
+minimap2 in a pipeline in two ways:
+
+1. **Swap the command.** Replace the aligner call itself.
+2. **Align first.** Align with bonefish, then hand the sorted BAM to a pipeline that accepts
+   aligned input.
+
+## Replacing a minimap2 command
+
+| minimap2 | bonefish |
+|---|---|
+| `minimap2 -ax lr:hq -y -t 16 ref.fa reads.fq \| samtools sort -o out.bam && samtools index out.bam` | `bonefish -t 16 -y -o out.bam ref.fa reads.fq` |
+| `samtools fastq -T '*' reads.bam \| minimap2 -ax lr:hq -y ...` (dorado BAM) | `bonefish -t 16 -y -o out.bam ref.fa reads.bam` |
+| `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 ...` |
+
+Notes:
+
+- **`-o out.bam`** writes the sorted BAM and its `.bai` index in one step: the same records and
+  order as `samtools sort --write-index`.
+- **`-y`** keeps each read's tags (methylation `MM`/`ML`, `RG`, ...) and the input's `@RG`
+  header lines. Those headers carry dorado's basecalling model, which tools such as Clair3 read.
+- **`--junc-bed`** takes BED12 (one transcript per line) or BED6 (one intron per line). minimap2's
+  `paftools.js gff2bed` converts a GTF or GFF3 file.
+- **Indexes:** build one once per preset for repeated runs (see [usage](usage.md#indexes)).
+
+## ONT's wf-human-variation
+
+[wf-human-variation](https://github.com/epi2me-labs/wf-human-variation) accepts aligned BAMs.
+Align with bonefish, then run the workflow on the result; it skips its own alignment step:
+
+```bash
+bonefish -t 16 -y -o sample.bam GRCh38.fa reads.bam
+nextflow run epi2me-labs/wf-human-variation --bam sample.bam --ref GRCh38.fa \
+    --snp --sv --mod --str --cnv --sex XY
+```
+
+- **Same reference:** use the same reference file for bonefish and the workflow.
+- **Basecalling model:** with `-y` and a dorado BAM as input, the `@RG` lines carry the model. From
+  FASTQ input, pass it explicitly with `--override_basecaller_cfg`, e.g.
+  `dna_r10.4.1_e8.2_400bps_sup@v5.0.0`.
+- **Validation:** on HG002 chr20, every analysis of the workflow (small variants, SVs,
+  methylation, STRs, CNVs) gave the same or slightly better results from bonefish's BAM than from
+  minimap2's ([report](validation/2026-10-14-wf-human-variation/README.md)).
+
+## Variant callers and methylation
+
+bonefish's sorted, indexed BAMs have been tested with these tools:
+
+- **Small variants:** Clair3.
+- **Structural variants:** Sniffles2.
+- **Phasing:** longphase, WhatsHap.
+- **Methylation:** modkit, which needs `-y` so the `MM`/`ML` tags are kept.
+
+The BAMs are standard (they pass Picard ValidateSamFile), so other tools that read sorted BAMs
+should work too.
+
+Accuracy against minimap2's BAMs is in [benchmarks.md](benchmarks.md).
+
+## RNA and cDNA
+
+Spliced BAMs from `-x rna` or `-x cdna` are the input that genome-based transcript tools take, for example:
+
+- **Isoform discovery and quantification:** IsoQuant, bambu, StringTie2 (`-L`), FLAIR.
+- **Direct-RNA modifications:** modkit with `-y`; dorado's m6A and other calls are carried as
+  `MM`/`ML` tags.
+
+**Transcript strand:** the `ts:A` tag gives the transcript strand of each spliced alignment, as
+minimap2 does.
+
+**Not yet validated:** the presets are checked against minimap2, GENCODE junctions and simulated
+reads ([report](validation/2026-10-21-rna-cdna/README.md)), but not yet against these downstream
+tools.
+
+## Containers
+
+```bash
+docker build -t bonefish https://github.com/DylanCosto/bonefish.git
+docker run --rm -v "$PWD":/data bonefish bonefish -t 16 -y -o /data/out.bam /data/ref.fa /data/reads.bam
+```
+
+The image includes `procps`, which Nextflow needs for task monitoring.