From 7fe6ffa33ca6fb1380cc711f1a277a70f9fb34c6 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 23 Jul 2026 19:49:47 +0200 Subject: [PATCH 1/5] fix(ci): resolve clippy 1.97 lints and RUSTSEC advisories CI's clippy advanced to 1.97 (new lints fail the `-D warnings` gate) and rustsec flagged several advisories in the dependency tree. Both block all PRs. Clippy: - src/index/suffix_array.rs: assert!(x == 0) -> assert_eq! (manual_assert_eq) - src/io/sam.rs, tests/alignment_features.rs: byte-slice literals -> byte strings, e.g. [b'S', b'M'] -> *b"SM" (byte_char_slices) - src/quant/transcriptome.rs: if-let/else-return-None -> `?` (question_mark) Security audit (Cargo.lock only, no manifest change): - memmap2 0.9.10 -> 0.9.11 (RUSTSEC-2026-0186, out-of-bounds pointer offset) - crossbeam-epoch 0.9.18 -> 0.9.20 (RUSTSEC-2026-0204, invalid pointer deref) - anyhow 1.0.102 -> 1.0.104 (RUSTSEC-2026-0190, unsound downcast_mut) No behavioural change. Unblocks CI for open and future PRs. Co-Authored-By: Claude Opus 4.8 (1M context) --- Cargo.lock | 12 ++++++------ src/index/suffix_array.rs | 2 +- src/io/sam.rs | 4 ++-- src/quant/transcriptome.rs | 11 ++++------- tests/alignment_features.rs | 2 +- 5 files changed, 14 insertions(+), 17 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 7936618..6a6c361 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -78,9 +78,9 @@ dependencies = [ [[package]] name = "anyhow" -version = "1.0.102" +version = "1.0.104" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7f202df86484c868dbad7eaa557ef785d5c66295e41b460ef922eca0723b842c" +checksum = "330a5ed07fa54e4702c9d6c4174f74427fc0ef6e214bbd677ae50a5099946470" [[package]] name = "assert_cmd" @@ -288,9 +288,9 @@ dependencies = [ [[package]] name = "crossbeam-epoch" -version = "0.9.18" +version = "0.9.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" dependencies = [ "crossbeam-utils", ] @@ -639,9 +639,9 @@ checksum = "f52b00d39961fc5b2736ea853c9cc86238e165017a493d1d5c8eac6bdc4cc273" [[package]] name = "memmap2" -version = "0.9.10" +version = "0.9.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "714098028fe011992e1c3962653c96b2d578c4b4bce9036e15ff220319b1e0e3" +checksum = "d1219ed1b7f229ee7104d281dd01d6802fe28bb6e95d292942c4daacdeb798c0" dependencies = [ "libc", ] diff --git a/src/index/suffix_array.rs b/src/index/suffix_array.rs index 570a7dd..e1a6017 100644 --- a/src/index/suffix_array.rs +++ b/src/index/suffix_array.rs @@ -275,7 +275,7 @@ mod tests { let first_base = genome.sequence[first_pos as usize]; // In "AAB", the first suffix lexicographically is "A" (from pos 0 or 1) - assert!(first_base == 0); // A + assert_eq!(first_base, 0); // A } /// Differential: the caps-sa-backed builder must produce a `PackedArray` diff --git a/src/io/sam.rs b/src/io/sam.rs index 29716ea..3c94b40 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -1516,8 +1516,8 @@ mod tests { assert!(rgs.contains_key(&b"rg0"[..])); let map = rgs.get(&b"rg0"[..]).unwrap(); // SM and LB should be present as other_fields - let sm_tag = HeaderOtherTag::<_>::try_from([b'S', b'M']).unwrap(); - let lb_tag = HeaderOtherTag::<_>::try_from([b'L', b'B']).unwrap(); + let sm_tag = HeaderOtherTag::<_>::try_from(*b"SM").unwrap(); + let lb_tag = HeaderOtherTag::<_>::try_from(*b"LB").unwrap(); let sm: &[u8] = map.other_fields().get(&sm_tag).unwrap().as_ref(); let lb: &[u8] = map.other_fields().get(&lb_tag).unwrap().as_ref(); assert_eq!(sm, b"sample0"); diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index c4c86ae..d332418 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -1036,13 +1036,10 @@ fn align_to_one_transcript( } else { // Coalesce: extend the last projected exon by this block's length // (STAR adds block EX_L directly — does NOT update start position). - if let Some(last) = proj_exons.last_mut() { - let len = block_end - block_start; - last.genome_end += len; - last.read_end = block.read_end; - } else { - return None; - } + let last = proj_exons.last_mut()?; + let len = block_end - block_start; + last.genome_end += len; + last.read_end = block.read_end; } // Advance `ex` across any splice boundary BEFORE the next block, diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 27e3a8d..4dfe9c7 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -16,7 +16,7 @@ use tempfile::TempDir; /// LCG pseudo-random sequence generator (identical LCG to existing tests). fn lcg_seq(seed: u32, length: usize) -> Vec { - let bases: [u8; 4] = [b'A', b'C', b'G', b'T']; + let bases: [u8; 4] = *b"ACGT"; let mut state = seed; let mut seq = Vec::with_capacity(length); for _ in 0..length { From 242b770ec40c33f6beddd2bd97ca2b0eb044e5ce Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 23 Jul 2026 19:14:54 +0200 Subject: [PATCH 2/5] docs: add STAR-rs porting groundwork (attribution + conventions) Add LICENSES/STAR_RS_LICENSE (STAR-rs MIT notice, Copyright (c) 2026 Benjamin Demaille) and note it in LICENSES/README.md, for attribution of modules ported from STAR-rs. STAR-rs is dual MIT-OR-Apache; ported code is taken under its MIT option, compatible with this repo's MIT license. Add docs-old/dev/porting-from-star-rs.md documenting the STAR-rs to rustar-aligner type/module mapping, the licensing and cellranger reference-only caveat, the determinism decision (adopt STAR-rs's thread-invariant, no-rand model), and how to validate ports against native STAR via the existing test/ harness with DIVERGENCES.md as the acceptance spec. This is the base of a series of themed, dependency-stacked PRs bringing the STAR-rs feature surface (WASP, clipping, signal tracks, genome transforms, STARlong, STARsolo, ...) into rustar-aligner. Co-Authored-By: Claude Opus 4.8 (1M context) --- LICENSES/README.md | 8 ++- LICENSES/STAR_RS_LICENSE | 21 ++++++++ docs-old/dev/porting-from-star-rs.md | 75 ++++++++++++++++++++++++++++ 3 files changed, 103 insertions(+), 1 deletion(-) create mode 100644 LICENSES/STAR_RS_LICENSE create mode 100644 docs-old/dev/porting-from-star-rs.md diff --git a/LICENSES/README.md b/LICENSES/README.md index cb6189e..8afa7ce 100644 --- a/LICENSES/README.md +++ b/LICENSES/README.md @@ -1,3 +1,9 @@ The original MIT STAR license is included here for attribution to the original author Alexander Dobin -The LICENSE of this repo has been chosen to also be MIT to match that of the STAR repo https://github.com/alexdobin/STAR \ No newline at end of file +The LICENSE of this repo has been chosen to also be MIT to match that of the STAR repo https://github.com/alexdobin/STAR + +`STAR_RS_LICENSE` is the MIT license of [STAR-rs](https://github.com/BenjaminDEMAILLE/STAR-rs) +(Copyright (c) 2026 Benjamin Demaille), included for attribution of modules ported +from STAR-rs into rustar-aligner. STAR-rs is dual-licensed `MIT OR Apache-2.0`; the +ported code is taken under its MIT option, which is compatible with this repo's MIT +license. See `docs-old/dev/porting-from-star-rs.md` for the porting conventions. \ No newline at end of file diff --git a/LICENSES/STAR_RS_LICENSE b/LICENSES/STAR_RS_LICENSE new file mode 100644 index 0000000..fc97d60 --- /dev/null +++ b/LICENSES/STAR_RS_LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Benjamin Demaille + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/docs-old/dev/porting-from-star-rs.md b/docs-old/dev/porting-from-star-rs.md new file mode 100644 index 0000000..8cc57a8 --- /dev/null +++ b/docs-old/dev/porting-from-star-rs.md @@ -0,0 +1,75 @@ +# Porting features from STAR-rs into rustar-aligner + +This note records the conventions for bringing features over from +[STAR-rs](https://github.com/BenjaminDEMAILLE/STAR-rs) (a clean-room Rust re-port of +STAR 2.7.11b, validated byte-identical against native STAR) into rustar-aligner. +It is the shared reference for the themed, dependency-stacked porting PRs. + +## Why + +rustar-aligner and STAR-rs are two independent Rust reimplementations of STAR. +rustar-aligner already covers single/paired-end alignment, splice junctions, +two-pass, 4-tier chimeric, GeneCounts, TranscriptomeSAM and sorted BAM. STAR-rs +additionally implements STARsolo, STARlong, WASP, genome transforms, signal tracks, +adapter clipping, PE mate-overlap, SuperTranscriptome, liftOver and more. These +notes cover moving that surface over, one theme per PR. + +## Licensing + +STAR-rs is dual-licensed `MIT OR Apache-2.0` (Copyright (c) 2026 Benjamin Demaille). +Ported code is taken under the **MIT** option, compatible with this repo's MIT +license. `LICENSES/STAR_RS_LICENSE` preserves the STAR-rs MIT notice for attribution. + +**Never copy `vendor/cellranger-upstream` code from STAR-rs.** It is 10x Genomics +source-available (non-OSI) and kept in STAR-rs as algorithm *reference only*. When +porting STARsolo, reimplement the algorithms from the STAR C++ source, as STAR-rs did. + +## Type / module mapping + +STAR-rs is a 10-crate workspace (`star_core`, `star_index`, `star_seed`, +`star_align`, `star_cli`, ...); rustar-aligner is a single crate (`rustar_aligner`) +with a flat module tree. A ported file's `use star_*::...` header is repointed to +`crate::...`, and the core types are adapted field-by-field: + +| STAR-rs | rustar-aligner (`src/align/transcript.rs`) | +| --- | --- | +| `star_align::star_model::Transcript` | `crate::align::transcript::Transcript` | +| `Exon { r, g, len, i_frag }` | `Exon { read_start, genome_start, genome_end, i_frag }` (start+end, not start+len) | +| `tr.str_: u8` (0 fwd / 1 rev) | `tr.is_reverse: bool` | +| `tr.chr` | `tr.chr_idx` | +| `tr.n_exons()` | `tr.exons.len()` | +| `star_core::Cigar` (derived from exons late) | `Vec` (materialized in `Transcript`) | +| `star_index::load::StarGenome` (unified) | split `genome::Genome` + `index::GenomeIndex` + `index::suffix_array::SuffixArray` | + +Base encoding is identical in both: `A=0, C=1, G=2, T=3, N=4`, reverse complement in +the second half of the genome array. + +## Determinism + +rustar-aligner currently uses `rand` (`StdRng`, seeded by `--runRNGseed`) for +multimapper-order / tie-break randomness. STAR-rs instead uses a hand-rolled, +per-read-seeded `splitmix64` shuffle so output is **byte-identical regardless of +thread count** (a stronger correctness property than STAR's per-thread `mt19937`). + +**Decision for the porting effort:** adopt STAR-rs's thread-invariant model. Ported +code should not introduce new `rand` usage; the migration of the existing tie-break +path lands in its own PR before STARsolo (whose UMI dedup depends on deterministic +order). This is a behavioural change from STAR's per-thread RNG and is flagged here +for team discussion. + +## Validation + +STAR-rs's `DIVERGENCES.md` (entries D1-D24) enumerates every intentional difference +from the native-STAR oracle, each locked by a named test. Treat it as the acceptance +spec: a ported feature should match native STAR except where DIVERGENCES.md documents +an intentional difference. + +rustar-aligner already has a differential harness under `test/` (`compare_sam.py`, +`compare_pe.py`, `assess_faithfulness.py`, `compare_junctions.py`, +`compare_chimeric.py`, driven by `test/ci.sh`) that runs native STAR and diffs +outputs. Use it to validate each port against `STAR` on the yeast test set; it is a +local/manual harness and is not part of the GitHub `alls-green` CI gate. + +Every PR must still pass the standard gate: `cargo build`, `cargo test`, +`cargo clippy --all-targets -- -D warnings`, `cargo fmt --check`, and MSRV +`cargo +1.89 check`. From 3e84d9716503039dcbbec2247f20979754539b53 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 23 Jul 2026 19:47:43 +0200 Subject: [PATCH 3/5] feat: WASP allele-specific-mapping filter (single-end) Port STAR's WASP filter (--waspOutputMode SAMtag, --varVCFfile) from STAR-rs `crates/star-align/src/star_wasp.rs`. A uniquely-mapped read overlapping heterozygous SNVs is re-mapped with every alternative-allele combination; if all re-maps land on the identical locus it passes (vW:i:1), otherwise a failure code (2..7) is reported. Emits vW, and vA/vG when requested in --outSAMattributes. - src/wasp/mod.rs: load_vcf, wasp_type (re-map decision), wasp_variants (vG/vA), and record annotation. Works in base-code space; variant/read overlap is computed by walking the transcript CIGAR (same convention as build_md_tag), so it does not depend on the exon read-coordinate frame. - src/params: add --waspOutputMode (None|SAMtag) and --varVCFfile, validation (SAMtag requires a VCF), and vW/vA/vG in --outSAMattributes. - src/lib.rs: load the VCF once per run, re-map with a WASP-relaxed parameter set (match_nmin=0, multimap_score_range=1, multimap_nmax=10, matching STAR), and stamp vW/vA/vG onto the SAM records in the single-end path. Paired-end WASP is a follow-up. See docs-old/dev/porting-from-star-rs.md. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/lib.rs | 43 ++++- src/params/mod.rs | 41 ++++ src/params/sam.rs | 9 +- src/wasp/mod.rs | 483 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 574 insertions(+), 2 deletions(-) create mode 100644 src/wasp/mod.rs diff --git a/src/lib.rs b/src/lib.rs index c9a3aa9..73e1931 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -34,6 +34,7 @@ pub mod junction; pub mod mapq; pub mod quant; pub mod stats; +pub mod wasp; use log::info; use noodles::sam::alignment::record::cigar; @@ -962,6 +963,34 @@ fn align_reads_single_end( }; let mut bysj_meta: Vec = Vec::new(); + // WASP allele-specific filtering context: load the VCF once, and build the + // relaxed re-map parameter set. `None` when --waspOutputMode is not SAMtag. + let wasp_ctx: Option = + if params.wasp_output_mode == params::WaspOutputMode::SAMtag { + let vcf = params + .var_vcf_file + .as_ref() + .expect("validated: SAMtag requires --varVCFfile"); + let ctx = crate::wasp::WaspContext::load( + vcf, + &index.genome.chr_name, + &index.genome.chr_start, + params, + ) + .map_err(|source| error::Error::Io { + source, + path: vcf.clone(), + })?; + info!( + "WASP: loaded {} heterozygous SNVs from {}", + ctx.snps.len(), + vcf.display() + ); + Some(ctx) + } else { + None + }; + info!("Aligning reads..."); loop { // Sequential FASTQ reading (unavoidable bottleneck) @@ -1099,7 +1128,7 @@ fn align_reads_single_end( } } else if transcripts.len() <= max_multimaps { // Mapped (within multimap limit) - let records = SamWriter::build_alignment_records( + let mut records = SamWriter::build_alignment_records( &read.name, &clipped_seq, &clipped_qual, @@ -1108,6 +1137,18 @@ fn align_reads_single_end( params, n_for_mapq, )?; + // WASP allele-specific filtering: stamp vW/vA/vG on the records. + if let Some(ctx) = wasp_ctx.as_ref() { + crate::wasp::annotate_records_se( + &mut records, + &transcripts, + &clipped_seq, + &read.name, + &index, + ctx, + params.out_sam_attributes, + )?; + } for record in records { buffer.push(record); } diff --git a/src/params/mod.rs b/src/params/mod.rs index a63b5e8..76cbfdf 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -221,6 +221,26 @@ impl std::str::FromStr for TwopassMode { } } +/// STAR's `--waspOutputMode`: WASP allele-specific-mapping filter output. +#[derive(Debug, Clone, PartialEq, Eq, Default)] +pub enum WaspOutputMode { + #[default] + None, + /// Emit the `vW` SAM tag (and `vA`/`vG` when requested in `--outSAMattributes`). + SAMtag, +} + +impl std::str::FromStr for WaspOutputMode { + type Err = String; + fn from_str(s: &str) -> Result { + match s { + "None" => Ok(Self::None), + "SAMtag" => Ok(Self::SAMtag), + _ => Err(format!("unknown waspOutputMode value: '{s}'")), + } + } +} + // --------------------------------------------------------------------------- // Parameters struct // --------------------------------------------------------------------------- @@ -616,6 +636,15 @@ pub struct Parameters { #[arg(long = "twopass1readsN", default_value_t = -1, allow_hyphen_values = true)] pub twopass1_reads_n: i64, + // ── WASP allele-specific filtering ────────────────────────────────── + /// WASP output mode: None or SAMtag (emit the vW tag) + #[arg(long = "waspOutputMode", default_value = "None")] + pub wasp_output_mode: WaspOutputMode, + + /// VCF of heterozygous SNVs for WASP filtering (required with --waspOutputMode SAMtag) + #[arg(long = "varVCFfile")] + pub var_vcf_file: Option, + // ── Chimeric ──────────────────────────────────────────────────────── // ── Debug ─────────────────────────────────────────────────────── /// Filter for debug logging: only log detailed alignment info for reads matching this name @@ -917,6 +946,18 @@ impl Parameters { )); } + // WASP SAMtag mode requires a VCF of heterozygous SNVs; fold the vW bit + // into out_sam_attributes so the writer emits it (vA/vG stay opt-in). + if params.wasp_output_mode == WaspOutputMode::SAMtag { + if params.var_vcf_file.is_none() { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--waspOutputMode SAMtag requires --varVCFfile", + )); + } + params.out_sam_attributes |= SamAttributes::VW; + } + Ok(params) } diff --git a/src/params/sam.rs b/src/params/sam.rs index 43006f7..a1ed7b3 100644 --- a/src/params/sam.rs +++ b/src/params/sam.rs @@ -20,6 +20,10 @@ bitflags::bitflags! { const JI = 1 << 6; const XS = 1 << 7; const RG = 1 << 8; + // WASP allele-specific-mapping tags (require --waspOutputMode SAMtag). + const VW = 1 << 9; + const VA = 1 << 10; + const VG = 1 << 11; const STANDARD = Self::NH.bits() | Self::HI.bits() | Self::AS.bits() @@ -49,6 +53,9 @@ impl FromStr for SamAttributes { "jI" => Self::JI, "XS" => Self::XS, "RG" => Self::RG, + "vW" => Self::VW, + "vA" => Self::VA, + "vG" => Self::VG, other => return Err(format!("unknown --outSAMattributes token '{other}'")), }) } @@ -94,7 +101,7 @@ impl clap::Args for SamAttributes { .default_values(["Standard"]) .help( "SAM optional tags: Standard, All, None, or any combination of \ - NH HI AS NM nM MD jM jI XS RG.", + NH HI AS NM nM MD jM jI XS RG vW vA vG.", ), ) } diff --git a/src/wasp/mod.rs b/src/wasp/mod.rs new file mode 100644 index 0000000..792477a --- /dev/null +++ b/src/wasp/mod.rs @@ -0,0 +1,483 @@ +//! WASP allele-specific-mapping filter. +//! +//! Ported from STAR-rs `crates/star-align/src/star_wasp.rs`, which itself ports +//! STAR's `ReadAlign_waspMap.cpp`, `Transcript::variationAdjust` and +//! `Variation::loadVCF`. See `docs-old/dev/porting-from-star-rs.md`. +//! +//! A read overlapping heterozygous SNP(s) is re-mapped with every alternative +//! allele combination; if all re-maps land on the identical locus the read passes +//! (`vW:i:1`), otherwise a failure code is reported. Reads overlapping no variant +//! get no `vW` tag (code `-1`). +//! +//! This revision covers **single-end**; paired-end WASP is a follow-up. +//! +//! Everything works in rustar-aligner's base-code space (`A=0,C=1,G=2,T=3,N=4`), +//! and variant/read overlap is computed by walking the transcript CIGAR in +//! SAM/forward-genomic order (the same convention as `io::sam::build_md_tag`), +//! so it does not depend on the exon read-coordinate frame. + +use std::path::Path; + +use noodles::sam::alignment::record::cigar; +use noodles::sam::alignment::record::data::field::Tag; +use noodles::sam::alignment::record_buf::RecordBuf; +use noodles::sam::alignment::record_buf::data::field::Value; +use noodles::sam::alignment::record_buf::data::field::value::Array; + +use crate::align::read_align::align_read; +use crate::align::transcript::Transcript; +use crate::error::Error; +use crate::index::GenomeIndex; +use crate::io::fastq::complement_base; +use crate::params::{Parameters, SamAttributes}; + +/// One heterozygous SNV: absolute 0-based genomic `loci`, and `nt = [ref, allele0, +/// allele1]` as base codes. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct Snp { + pub loci: u64, + pub nt: [u8; 3], +} + +fn nt_code(b: u8) -> u8 { + match b { + b'A' | b'a' => 0, + b'C' | b'c' => 1, + b'G' | b'g' => 2, + b'T' | b't' => 3, + _ => 4, + } +} + +/// Reverse-complement a base-code sequence (`N`=4 preserved). +fn rc_codes(v: &[u8]) -> Vec { + v.iter().rev().map(|&b| complement_base(b)).collect() +} + +/// STAR `scanVCF`: single-char REF and ALT, genotype not homozygous-reference and +/// heterozygous (`gt0 != gt1`); `nt = [REF, altV[gt0], altV[gt1]]` with +/// `altV = [REF, ALT...]`. Returns SNPs sorted by absolute genomic position (only +/// ACGT variants on known chromosomes). `chr_names`/`chr_starts` come from the +/// loaded genome. +pub fn load_vcf( + path: &Path, + chr_names: &[String], + chr_starts: &[u64], +) -> std::io::Result> { + let text = std::fs::read_to_string(path)?; + let mut snps = Vec::new(); + for line in text.lines() { + if line.starts_with('#') || line.is_empty() { + continue; + } + let f: Vec<&str> = line.split('\t').collect(); + if f.len() < 10 { + continue; + } + let (chr, pos, r#ref, alt, sample) = (f[0], f[1], f[3], f[4], f[9]); + if r#ref.len() != 1 { + continue; + } + let alt_v: Vec<&str> = std::iter::once(r#ref).chain(alt.split(',')).collect(); + if alt_v.iter().skip(1).any(|a| a.len() != 1) { + continue; + } + let sb = sample.as_bytes(); + if sb.len() < 3 { + continue; + } + let (Some(gt0), Some(gt1)) = ((sb[0] as char).to_digit(10), (sb[2] as char).to_digit(10)) + else { + continue; + }; + // Skip homozygous-reference / homozygous genotypes (heteroOnly default). + if (gt0 == 0 && gt1 == 0) || gt0 == gt1 { + continue; + } + let Some(chr_idx) = chr_names.iter().position(|c| c == chr) else { + continue; + }; + let (Some(a0), Some(a1)) = ( + alt_v.get(gt0 as usize).map(|s| s.as_bytes()[0]), + alt_v.get(gt1 as usize).map(|s| s.as_bytes()[0]), + ) else { + continue; + }; + let nt = [nt_code(r#ref.as_bytes()[0]), nt_code(a0), nt_code(a1)]; + if nt.iter().any(|&c| c >= 4) { + continue; + } + let Ok(pos1) = pos.parse::() else { + continue; + }; + snps.push(Snp { + loci: pos1 - 1 + chr_starts[chr_idx], + nt, + }); + } + snps.sort_by_key(|s| s.loci); + Ok(snps) +} + +/// The matched-allele code for a read base at a SNP: `1`/`2` for the SNP's two +/// alleles, `3` if neither, `4` if the read base is `N`. +fn classify_allele(snp: &Snp, read_code: u8) -> u8 { + if read_code > 3 { + 4 + } else if snp.nt[1] == read_code { + 1 + } else if snp.nt[2] == read_code { + 2 + } else { + 3 + } +} + +/// Variants overlapping an alignment, by walking the CIGAR in SAM/forward-genomic +/// order. `sam_codes` is the read in that frame (reverse-complemented for a reverse +/// alignment). Returns `(snp_index, read_pos_in_sam_frame, allele)` per overlapping +/// SNP, in genomic order. `snps` must be sorted by `loci`. +fn variation_overlap(snps: &[Snp], sam_codes: &[u8], tr: &Transcript) -> Vec<(usize, usize, u8)> { + use cigar::op::Kind; + let mut out = Vec::new(); + let mut genome_pos = tr.genome_start; + let mut read_pos: usize = 0; + for op in &tr.cigar { + match op.kind() { + Kind::Match | Kind::SequenceMatch | Kind::SequenceMismatch => { + for _ in 0..op.len() { + if let Ok(isnp) = snps.binary_search_by(|s| s.loci.cmp(&genome_pos)) { + let read_code = sam_codes.get(read_pos).copied().unwrap_or(4); + out.push((isnp, read_pos, classify_allele(&snps[isnp], read_code))); + } + genome_pos += 1; + read_pos += 1; + } + } + Kind::Deletion | Kind::Skip => { + genome_pos += op.len() as u64; + } + Kind::Insertion | Kind::SoftClip => { + read_pos += op.len(); + } + Kind::HardClip | Kind::Pad => {} + } + } + out +} + +/// STAR `Transcript::variationAdjust`: the `vG`/`vA` tag values for an alignment +/// overlapping heterozygous SNP(s). `vg` = each SNP's chromosome-relative 0-based +/// genomic coordinate; `va` = the matched-allele code, in genomic order. Empty when +/// the alignment overlaps no variant. `read_codes` is the forward read (base codes); +/// `chr_start` is the alignment chromosome's genome offset. +pub fn wasp_variants( + chr_start: u64, + snps: &[Snp], + read_codes: &[u8], + tr: &Transcript, +) -> (Vec, Vec) { + let sam_codes = if tr.is_reverse { + rc_codes(read_codes) + } else { + read_codes.to_vec() + }; + let mut vg = Vec::new(); + let mut va = Vec::new(); + for (isnp, _rp, allele) in variation_overlap(snps, &sam_codes, tr) { + vg.push((snps[isnp].loci - chr_start) as i32); + va.push(allele); + } + (vg, va) +} + +fn same_exons(a: &Transcript, b: &Transcript) -> bool { + a.exons.len() == b.exons.len() + && a.exons.iter().zip(b.exons.iter()).all(|(x, y)| { + x.read_start == y.read_start + && x.genome_start == y.genome_start + && x.genome_end == y.genome_end + }) +} + +/// STAR `waspMap` (single-end): the `vW` code for one alignment, or `-1` for no `vW` +/// tag (read overlaps no variant). `1` = passed; `2` multimaps; `3` variant base is +/// `N`; `4` a remap is unmapped; `5` a remap multimaps; `6` a remap lands elsewhere; +/// `7` too many variants. `n_tr` is the read's mapped-locus count; `remap_params` +/// is the WASP-relaxed [`wasp_remap_params`]. +pub fn wasp_type( + index: &GenomeIndex, + snps: &[Snp], + read_codes: &[u8], + read_name: &str, + tr: &Transcript, + n_tr: usize, + remap_params: &Parameters, +) -> Result { + let sam_codes = if tr.is_reverse { + rc_codes(read_codes) + } else { + read_codes.to_vec() + }; + let vars = variation_overlap(snps, &sam_codes, tr); + if vars.is_empty() { + return Ok(-1); + } + if n_tr > 1 { + return Ok(2); + } + if vars.len() > 10 { + return Ok(7); + } + if vars.iter().any(|&(_, _, allele)| allele > 3) { + return Ok(3); + } + + let actual: Vec = vars.iter().map(|&(_, _, a)| a).collect(); + let n = vars.len(); + // All 2^n allele combinations of {1, 2}. + for mask in 0..(1u32 << n) { + let combo: Vec = (0..n) + .map(|i| if mask & (1 << i) != 0 { 2 } else { 1 }) + .collect(); + if combo == actual { + continue; + } + // Flip each variant's base to the combination's allele in the SAM frame, + // then map back to the forward read and re-align. + let mut modified = sam_codes.clone(); + for (iv, &(isnp, read_pos, _)) in vars.iter().enumerate() { + modified[read_pos] = snps[isnp].nt[combo[iv] as usize]; + } + let remap_read = if tr.is_reverse { + rc_codes(&modified) + } else { + modified + }; + let (trs, _chim, _n_for_mapq, _reason) = + align_read(&remap_read, read_name, index, remap_params)?; + if trs.is_empty() { + return Ok(4); // remap unmapped or too-many-loci + } + if trs.len() > 1 { + return Ok(5); + } + if !same_exons(&trs[0], tr) { + return Ok(6); + } + } + Ok(1) +} + +/// The parameter set used for WASP re-mapping: relaxed match filter and a fixed +/// multimap window, matching STAR's `waspMap` (STAR-rs passes `match_nmin=0`, +/// `multimap_score_range=1`, `multimap_nmax=10`). Built once and reused per read. +pub fn wasp_remap_params(params: &Parameters) -> Parameters { + let mut p = params.clone(); + p.out_filter_match_nmin = 0; + p.out_filter_match_nmin_over_lread = 0.0; + p.out_filter_multimap_score_range = 1; + p.out_filter_multimap_nmax = 10; + p +} + +/// Insert `vW`/`vA`/`vG` tags on a record, gated on the requested attributes. +fn insert_wasp_tags(record: &mut RecordBuf, vw: i32, vg: &[i32], va: &[u8], attrs: SamAttributes) { + let data = record.data_mut(); + if attrs.contains(SamAttributes::VW) { + data.insert(Tag::new(b'v', b'W'), Value::from(vw)); + } + if attrs.contains(SamAttributes::VA) && !va.is_empty() { + data.insert( + Tag::new(b'v', b'A'), + Value::Array(Array::Int8(va.iter().map(|&x| x as i8).collect())), + ); + } + if attrs.contains(SamAttributes::VG) && !vg.is_empty() { + data.insert( + Tag::new(b'v', b'G'), + Value::Array(Array::Int32(vg.to_vec())), + ); + } +} + +/// Loaded WASP state for a run: the heterozygous SNVs and the relaxed re-map +/// parameter set. Built once and shared read-only across the per-read parallel loop. +pub struct WaspContext { + pub snps: Vec, + pub remap_params: Parameters, +} + +impl WaspContext { + /// Load the VCF and build the WASP re-map parameters for a run. + pub fn load( + vcf_path: &Path, + chr_names: &[String], + chr_starts: &[u64], + params: &Parameters, + ) -> std::io::Result { + Ok(Self { + snps: load_vcf(vcf_path, chr_names, chr_starts)?, + remap_params: wasp_remap_params(params), + }) + } +} + +/// Compute WASP tags for a single-end read and stamp them onto its already-built +/// SAM records. `records[i]` corresponds to `transcripts[i]`. A uniquely-mapped read +/// is re-mapped (computing `vW`); a multi-mapped read overlapping variants gets +/// `vW:i:2`. Reads overlapping no variant are left untagged. +pub fn annotate_records_se( + records: &mut [RecordBuf], + transcripts: &[Transcript], + read_codes: &[u8], + read_name: &str, + index: &GenomeIndex, + ctx: &WaspContext, + attrs: SamAttributes, +) -> Result<(), Error> { + let snps = &ctx.snps; + if transcripts.len() == 1 { + let tr = &transcripts[0]; + let vw = wasp_type(index, snps, read_codes, read_name, tr, 1, &ctx.remap_params)?; + if vw != -1 { + let chr_start = index.genome.chr_start[tr.chr_idx]; + let (vg, va) = wasp_variants(chr_start, snps, read_codes, tr); + if let Some(rec) = records.get_mut(0) { + insert_wasp_tags(rec, vw, &vg, &va, attrs); + } + } + } else { + for (rec, tr) in records.iter_mut().zip(transcripts.iter()) { + let chr_start = index.genome.chr_start[tr.chr_idx]; + let (vg, va) = wasp_variants(chr_start, snps, read_codes, tr); + if !vg.is_empty() { + insert_wasp_tags(rec, 2, &vg, &va, attrs); + } + } + } + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::align::transcript::{Exon, Transcript}; + use noodles::sam::alignment::record::cigar::{Op, op::Kind}; + + fn tr_fwd_50m(genome_start: u64) -> Transcript { + Transcript { + chr_idx: 0, + genome_start, + genome_end: genome_start + 50, + is_reverse: false, + exons: vec![Exon { + genome_start, + genome_end: genome_start + 50, + read_start: 0, + read_end: 50, + i_frag: 0, + }], + cigar: vec![Op::new(Kind::Match, 50)], + score: 49, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: vec![], + junction_annotated: vec![], + read_seq: vec![], + } + } + + #[test] + fn nt_code_maps_bases() { + assert_eq!(nt_code(b'A'), 0); + assert_eq!(nt_code(b'c'), 1); + assert_eq!(nt_code(b'G'), 2); + assert_eq!(nt_code(b'T'), 3); + assert_eq!(nt_code(b'N'), 4); + } + + #[test] + fn rc_codes_reverses_and_complements() { + // A C G T N (0 1 2 3 4) -> reverse-complement -> N A C G T (4 0 1 2 3) + assert_eq!(rc_codes(&[0, 1, 2, 3, 4]), vec![4, 0, 1, 2, 3]); + } + + #[test] + fn classify_allele_codes() { + let snp = Snp { + loci: 10, + nt: [0, 0, 2], + }; // ref A, alleles A/G + assert_eq!(classify_allele(&snp, 0), 1); // A -> allele 1 + assert_eq!(classify_allele(&snp, 2), 2); // G -> allele 2 + assert_eq!(classify_allele(&snp, 1), 3); // C -> neither + assert_eq!(classify_allele(&snp, 4), 4); // N + } + + #[test] + fn load_vcf_parses_het_snp() { + use std::io::Write as _; + let mut f = tempfile::NamedTempFile::new().unwrap(); + writeln!(f, "##fileformat=VCFv4.2").unwrap(); + writeln!( + f, + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1" + ) + .unwrap(); + // het 0/1 at chrI:100 (A->G); hom-ref skipped; multi-char REF skipped + writeln!(f, "chrI\t100\t.\tA\tG\t.\t.\t.\tGT\t0|1").unwrap(); + writeln!(f, "chrI\t200\t.\tA\tG\t.\t.\t.\tGT\t0|0").unwrap(); + writeln!(f, "chrI\t300\t.\tAC\tG\t.\t.\t.\tGT\t0|1").unwrap(); + let names = vec!["chrI".to_string()]; + let starts = vec![0u64]; + let snps = load_vcf(f.path(), &names, &starts).unwrap(); + assert_eq!(snps.len(), 1); + assert_eq!(snps[0].loci, 99); // 0-based, chr_start 0 + assert_eq!(snps[0].nt, [0, 0, 2]); // ref A, allele0 A, allele1 G + } + + #[test] + fn variation_overlap_forward() { + // SNP at genomic 110, transcript maps read[0..50] to genome[100..150] fwd. + let snps = vec![Snp { + loci: 110, + nt: [0, 0, 2], + }]; + let tr = tr_fwd_50m(100); + let mut read = vec![0u8; 50]; + read[10] = 2; // read base at the SNP is G (allele 2) + let out = variation_overlap(&snps, &read, &tr); + assert_eq!(out, vec![(0, 10, 2)]); + } + + #[test] + fn wasp_variants_reverse_strand_frame() { + // Reverse alignment: sam frame is rc(read). SNP at genome 110. + let snps = vec![Snp { + loci: 110, + nt: [0, 0, 2], + }]; + let mut tr = tr_fwd_50m(100); + tr.is_reverse = true; + // sam_codes[10] must be the read base at the SNP; build read so rc(read)[10]=G(2). + let mut sam = vec![0u8; 50]; + sam[10] = 2; + let read = rc_codes(&sam); // forward read whose rc is `sam` + let (vg, va) = wasp_variants(0, &snps, &read, &tr); + assert_eq!(vg, vec![110]); // chr-relative coord + assert_eq!(va, vec![2]); + } + + #[test] + fn wasp_variants_empty_when_no_overlap() { + let snps = vec![Snp { + loci: 500, + nt: [0, 0, 2], + }]; + let tr = tr_fwd_50m(100); + let (vg, va) = wasp_variants(0, &snps, &[0u8; 50], &tr); + assert!(vg.is_empty() && va.is_empty()); + } +} From 24047b94f5c3b5509b85dfd5caa1b4c07d683b68 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 24 Jul 2026 01:52:58 +0200 Subject: [PATCH 4/5] docs: condense CLAUDE.md status + source-layout sections Replace the single giant phase-log paragraph and the stale hardcoded test count / known-issues list with a short faithfulness summary and a pointer to ROADMAP.md / CHANGELOG.md / docs-old for per-phase history. Refresh the source layout (params/ split, cpu.rs, sa_build.rs, sjdb_insert.rs, transcriptome.rs, log.rs) and the dependency list (caps-sa, mimalloc, rand 0.10, tempfile as a runtime dep) to match the current tree. Documentation only. Co-Authored-By: Claude Opus 4.8 (1M context) --- CLAUDE.md | 116 +++++++++++++++++++++++------------------------------- 1 file changed, 49 insertions(+), 67 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 3ff1099..b00c57f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -32,7 +32,11 @@ Always run `cargo clippy --all-targets`, `cargo fmt --check`, and `cargo test` b ## Current Status -**396 tests passing, 0 clippy warnings.** SE: 8613/8926 compare_sam.py (96.5%; note: lower due to seeded-RNG tie-break PR diverging from STAR's mt19937), **99.815% faithfulness (tie-adjusted)** (8611/8627 non-tie reads exact), 299 tie-breaking diffs excluded. 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 rustar-aligner-only SE reads**. PE: **8390 both-mapped** (STAR: 8390), **0 half-mapped**, 0 MAPQ inflations / 0 deflations, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded), **0 proper-pair diffs**, **0 NH diffs**. Phase 17.A: `scoreSeedBest` pre-extension. Phase 17.B: per-mate seeding. Phase 17.C: STAR-faithful SCORE-GATE + mappedFilter. Phase 17.D: combined-span penalty fix + dedup ordering. Phase 17.8: `--quantMode GeneCounts`. Phase E fix (2026-04-21): mate_id-aware diagonal dedup. Phase E2 (2026-04-22): STAR-faithful combined-read seeding. Phase E3 (2026-04-22): combined-threshold half-mapped fallback. Phase E4 (2026-04-22): PE-CHECK2 unconditional. Phase E5 (2026-04-23): split_combined_wt n_mismatch propagation. Phase E6 (2026-04-24): tie-adjusted faithfulness metric in assess_faithfulness.py. Phase F1: --runRNGseed + seeded primary tie-break (PR #5). Phase F2: --outSAMattrRGline (PR #6). Phase F3: --quantMode TranscriptomeSAM (PR #7). Phase F4: SJDB insertion into Genome+SA at genomeGenerate (PR #8). Phase G1 (2026-04-29): junction_shifts fix in split_combined_wt (rDNA cross-copy false-splice filter). Phase G2 (2026-04-29): MAX_RECURSION 10k→100k + sa_pos_to_forward overflow fix (ERR12389696.7118031 NH=3→9). Phase 17.2 (2026-04-29): coordinate-sorted BAM output (`--outSAMtype BAM SortedByCoordinate` → `Aligned.sortedByCoord.out.bam`). Phase 17.4 (2026-04-29): `--outReadsUnmapped Fastx` → `Unmapped.out.mate1` / `Unmapped.out.mate2`; writes unmapped + TooManyLoci reads; PE writes both mates for fully-unmapped and half-mapped pairs. Phase 17.6 (2026-05-01): `--outStd SAM/BAM_Unsorted/BAM_SortedByCoordinate` — routes primary alignment output to stdout via `Box` trait dispatch; `SamStdoutWriter`, `BamStdoutWriter`, `SortedBamStdoutWriter` in sam.rs/bam.rs; verified with samtools pipe (967 records). Phase G3 (2026-05-01): SA tie-breaking fix — `compare_suffixes` tie-breaker changed from `pos_b.cmp(&pos_a)` to `packed_a.cmp(&packed_b)` (ascending by packed SA value with strand bit); rustar-aligner SA is now **byte-for-byte identical** to STAR's SA for the yeast genome (10,862 → 0 entry diffs). diff AS: 6→4 cases (4 remaining are rustar-aligner improvements: .844151 VIII 0mm vs STAR VII 6mm, .4972950 spliced vs unspliced mate2). Phase 17.3 (2026-05-01): PE chimeric detection — `detect_inter_mate_chimeric` in `chimeric/detect.rs`; intra-mate multi-cluster chimeric via cluster splitting + mate2 read_pos adjustment; inter-mate chimeric for discordant pairs (diff chr, same strand, or >1Mb); `align_paired_read` returns 4-tuple including `Vec`; no benchmark regression (8390 both-mapped, 0 half-mapped). Phase 17.11 (2026-05-01): `--chimOutType WithinBAM` — chimeric alignments written as supplementary records (FLAG 0x800) in primary BAM; donor record has full SEQ + SA tag; acceptor has FLAG 0x800 + SA tag + empty SEQ; `build_within_bam_records` in `chimeric/output.rs`; `chim_out_junctions()` / `chim_out_within_bam()` helpers in params.rs; supports mixed `--chimOutType Junctions WithinBAM`. Phase 17.7 (2026-05-01): GTF tag parameters — `--sjdbGTFchrPrefix`, `--sjdbGTFfeatureExon`, `--sjdbGTFtagExonParentTranscript`, `--sjdbGTFtagExonParentGene`; `_configured` variants in `junction/gtf.rs`, `quant/mod.rs`, `quant/transcriptome.rs`, `junction/mod.rs`; all 4 production paths thread params; backward-compat wrappers preserve zero test disruption. Phase 17.9 (2026-05-01): `--outBAMcompression` (BGZF level -1–9, default 1; -1/0=NONE, 1-8=flate2 levels, ≥9=BEST) + `--limitBAMsortRAM` (bytes, 0=unlimited; aborts sort if ~400 bytes/record estimate exceeds limit); `bgzf_compression()` + `make_bgzf_writer()` helpers in `io/bam.rs`; threaded through all 4 BAM writers (unsorted file, sorted file, unsorted stdout, sorted stdout). PE chimericDetectionOld (2026-05-01): per-mate `detect_chimeric_old` called on `all_m1_transcripts` / `all_m2_transcripts` pools after `filter_paired_transcripts` in `read_align.rs`. Phase 17.12 (2026-05-01): BySJout disk buffering — `BySJReadMeta` struct + `NamedTempFile` SAM temp file replaces `Vec`; `create_bysj_writer` / `bysj_write_records` / `bysj_read_n_records` helpers in `io/sam.rs`; `tempfile` moved to `[dependencies]`. Phase 17.13 (2026-05-01): 8 integration tests in `tests/alignment_features.rs` — synthetic 20kb genome with planted GT-AG intron; tests cover BAM output, PE alignment, spliced reads, BySJout, GeneCounts, unmapped output, two-pass mode. Phase 12.2 (2026-05-04): SE chimeric Tier 1b soft-clip re-mapping — `detect_from_soft_clips` in `chimeric/detect.rs` re-seeds the primary alignment's soft-clipped bases when `detect_chimeric_old` finds no partner; `adjust_read_positions` helper shifts sub-seq coords into full-read space for right clips; called as Step 3c in `read_align.rs`. Phase 17.10 (2026-05-04): Chimeric Tier 3 — `detect_from_chimeric_residuals` in `chimeric/detect.rs` re-seeds outer uncovered read regions (before donor / after acceptor) of each found chimeric pair; enables 3-way gene-fusion detection; called as Step 3d in `read_align.rs`. See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs-old/](docs-old/) for per-phase development notes. The published Astro Starlight docs site is in [docs/](docs/). +End-to-end SE + PE RNA-seq alignment with splice-junction detection, two-pass mode, 4-tier chimeric detection, GeneCounts / TranscriptomeSAM quantification, SortedByCoordinate BAM, and multi-threaded processing. 0 clippy warnings. + +**Faithfulness vs STAR** (10k yeast reads, ERR12389696, 150 bp): SE **99.815% tie-adjusted** (8611/8627 non-tie reads exact; 0 STAR-only, 0 rustar-aligner-only, 1 CIGAR-only tie). PE **8390 both-mapped** (matches STAR), 0 half-mapped, **99.883% tie-adjusted** exact; 0 MAPQ inflations/deflations, 0 NH diffs, 0 proper-pair diffs. The SA is byte-for-byte identical to STAR's for the yeast genome. Remaining raw disagreements are verified genuine ties (SA-order / seeded-RNG tie-break — rustar-aligner uses `StdRng`, not STAR's mt19937). + +Phase-by-phase history (what changed, when, and why) lives in **[ROADMAP.md](ROADMAP.md)** and **[CHANGELOG.md](CHANGELOG.md)**; per-phase dev notes in **[docs-old/](docs-old/)**; the published Astro Starlight docs site is in **[docs/](docs/)**. Update those, not this section, when adding phases. ## Source Layout @@ -40,7 +44,10 @@ Always run `cargo clippy --all-targets`, `cargo fmt --check`, and `cargo test` b src/ main.rs -- Thin entry: parse CLI (clap), init logging, call lib::run() lib.rs -- run() dispatches on RunMode (AlignReads | GenomeGenerate) - params.rs -- ~52 STAR CLI params via clap derive, --camelCase long names + cpu.rs -- Runtime CPU/SIMD feature-detection guard; VERSION_BODY from build.rs + params/ + mod.rs -- STAR CLI params via clap derive, --camelCase long names, Parameters::validate() + sam.rs -- SAM-attribute / outSAMtype parsing helpers error.rs -- Error enum with thiserror (Parameter, Io, Fasta, Index, Alignment, Gtf) mapq.rs -- MAPQ calculation (STAR lookup table: n=1→255, n=2→3, n≥5→0) stats.rs -- Alignment statistics, Log.final.out writer, UnmappedReason enum @@ -50,7 +57,9 @@ src/ index/ mod.rs -- GenomeIndex (build + load + write) packed_array.rs-- Variable-width bit packing (1-64 bits per element) - suffix_array.rs-- SA construction, custom comparator, strand encoding + packed_stream.rs- Streaming writer for the PackedArray bit format (SA build) + suffix_array.rs-- SA data structure + strand encoding + naive oracle for tests + sa_build.rs -- STAR-faithful SA construction via the caps-sa crate (sentinel transform, in-mem vs ext-mem) sa_index.rs -- K-mer lookup table (35-bit entries with flags) io.rs -- Load index from disk (Genome, SA, SAindex) align/ @@ -63,14 +72,17 @@ src/ io/ mod.rs -- Module exports fastq.rs -- FASTQ reader (plain + gzip, noodles wrapper) - sam.rs -- SAM writer (header + records, noodles wrapper) + sam.rs -- SAM writer (header + records, noodles wrapper), BySJout disk buffering bam.rs -- BAM writer (BGZF compression, streaming unsorted + in-memory sorted output) + log.rs -- STAR-compatible Log.out / Log.progress.out writers junction/ mod.rs -- GTF parsing, junction database, motif detection, two-pass filtering sj_output.rs -- SJ.out.tab writer + sjdb_insert.rs -- SJDB insertion into Genome+SA at genomeGenerate (ports sjdbPrepare/sjdbBuildIndex.cpp) gtf.rs -- GTF parser (internal) quant/ mod.rs -- Gene-level read counting (--quantMode GeneCounts, ReadsPerGene.out.tab) + transcriptome.rs- Transcriptome-coordinate SAM output (--quantMode TranscriptomeSAM) chimeric/ mod.rs -- Module exports detect.rs -- Chimeric detection (Tier 1: transcript-pair, Tier 1b: soft-clip re-seed, Tier 2: multi-cluster, Tier 3: residual re-seed) @@ -101,79 +113,49 @@ src/ ## Dependencies +See `Cargo.toml` for the authoritative list; the load-bearing choices: + ```toml [dependencies] -clap = { version = "4", features = ["derive"] } -anyhow = "1" -thiserror = "2" -log = "0.4" -env_logger = "0.11" -memmap2 = "0.9" -byteorder = "1" -noodles = { version = "0.80", features = ["fastq", "sam", "bam", "bgzf"] } -bstr = "1" -flate2 = "1" -rayon = "1" -dashmap = "6" -chrono = "0.4" +clap = { version = "4", features = ["derive"] } # CLI parsing (derive) +anyhow = "1" # top-level error propagation +thiserror = "2" # Error enum +noodles = { version = "0.111", features = ["fastq", "sam", "bam", "bgzf"] } # SAM/BAM/FASTQ IO +memmap2 = "0.9" # mmap the genome/SA index +caps-sa = "0.6" # suffix-array construction (see index/sa_build.rs) +rayon = "1" # per-read parallel alignment +dashmap = "6" # concurrent junction/gene maps +rand = "0.10" # --runRNGseed tie-breaking (StdRng, NOT mt19937) +tempfile = "3" # BySJout disk buffering (now a runtime dep) +bitflags = "2" # SAM flag bitsets +shlex = "2" # split parametersFiles / command strings +flate2 = "1" # gzip FASTQ + BGZF levels +mimalloc = { version = "0.1", default-features = false } # global allocator (bounds peak RSS, faster small allocs) +# also: log, env_logger, byteorder, bstr, chrono + +[build-dependencies] # build.rs stamps OS/arch/SIMD/git-hash into VERSION_BODY +chrono = { version = "0.4", default-features = false, features = ["clock"] } [dev-dependencies] -tempfile = "3" assert_cmd = "2" predicates = "3" ``` -## Testing Pattern - -- Unit tests: `#[cfg(test)]` in each module -- Integration tests: `tests/` directory (future phases) -- Every phase uses differential testing against STAR where applicable -- Test data tiers: synthetic micro-genome → chr22 → full human genome - -**Current test status**: 364/364 tests passing (359 unit + 5 integration), 0 clippy warnings - -## Known Issues — Disagreement Root Causes (10k SE yeast) - -**299 total position disagreements — ALL verified as genuine ties** (SA-order ties + seeded-RNG tie-break divergence from STAR's mt19937): - -Both tools find identical alignment sets for all 299 disagreements. Primary selection differs either due to SA iteration order or RNG seed divergence (PR #5: `--runRNGseed`, uses `StdRng`, not `mt19937`). - -- **100+ diff-chr ties** — same set of alignments, different repeat copy chosen as primary. -- **27+ same-chr ties** — same alignment set, different primary due to tie-breaking. - -**1 CIGAR-only disagreement (same position, different CIGAR):** -- `ERR12389696.13573895`: both tools align to XV:218357 MAPQ=255, but rustar-aligner gives `100M1I45M4S` (insertion at read pos 100) while STAR gives `108M1I37M4S` (insertion at 108). Root cause: both alignments score AS=133. The 71-base seed is found at RC pos 29 (rustar-aligner) vs RC pos 37 (STAR) due to different Lmapped chain paths through a long homopolymer region. Same diagonal, different starting position → different insertion placement. Seed-level tie. +**Release profile is deliberately slow to build**: `lto = "fat"`, `codegen-units = 1`, `strip = true` in `[profile.release]` (tuned for the seed-search / DP inner loops). For iterating, prefer `cargo build` (debug) or `cargo test`; CI overrides via `CARGO_PROFILE_RELEASE_LTO=thin` / `CODEGEN_UNITS=16`. -**1 truly actionable SE issue:** - -1. **CIGAR-only insertion placement (1 read)** — `ERR12389696.13573895` (see above). Requires matching STAR's exact Lmapped chain to fix. - -Previously listed issues now resolved: -- `ERR12389696.18296181` (rustar-aligner-only false splice): filtered out by score threshold (score=95/92 < 98 threshold for 150bp read). Both the 1-exon soft-clip AND the 2-exon false splice fail the minimum score gate — read correctly unmapped. -- `ERR12389696.13766843` (STAR-only NM=10 read): rustar-aligner already maps this correctly (VII:24449, 28S121M1S, MAPQ=255) — it was incorrectly listed as "STAR-only". - -**Phase 16.38 (STAR-faithful filter ordering):** -- Moved dedup + score-range filter (multMapSelect) to run BEFORE quality filters (mappedFilter). STAR's ordering: `multMapSelect → mappedFilter`. Old code ran quality filters first, which could remove the high-scoring primary leaving a lower-scoring secondary as apparent "best". No benchmark change (the specific false splice read was already filtered by score threshold), but structurally correct for edge cases. - -**Phase 16.37 (alignIntronMax=0 fix) resolved:** -- `ERR12389696.5825571`: now aligns as `XV:80779 121M607028N13M16S` (exact match with STAR). Root cause: rustar-aligner computed a finite intron limit of 589824 when `alignIntronMax=0`, blocking the 607kb intron. STAR uses `>0` guard in `stitchAlignToTranscript.cpp` line 100 — alignIntronMax=0 means no limit. Fix: use `u32::MAX` sentinel in score.rs. -- `ERR12389696.16030539`: MAPQ inflation resolved. Both tools now find two alignments (XV:121224 128M925400N10M12S and XV:598336 128M448288N10M12S), both MAPQ=3. Primary selection differs (tie). MAPQ inflate: 1→0. - -See [ROADMAP.md](ROADMAP.md) and [docs-old/](docs-old/) for full issue tracking. - -## PE Status (Updated 2026-04-29 — Phase G2) - -**Current**: **PE both-mapped = 8390** (STAR: 8390, exact match!), **half-mapped = 0**, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded). MAPQ inflations: **0**, deflations: 0. NH diffs: **0**. 1 STAR-only mate (`.18919121`, SA-level diff), 1 rustar-aligner-only mate (`.6302610`, pre-existing FP). - -**Phase G2** (2026-04-29): `MAX_RECURSION` 10,000→100,000 + `sa_pos_to_forward` overflow fix. `ERR12389696.7118031` was the sole source of both NH diffs and MAPQ inflations (NH=3 vs STAR's NH=9, rustar-aligner MAPQ=1 vs STAR MAPQ=0). Root cause: the 47-WA rDNA cluster (4 copies × multiple seeds per mate) exhausted 10k recursions before exploring the 4th within-copy pair. Fix: raise the per-cluster recursion budget from 10k to 100k. Also fixed: `sa_pos_to_forward` underflow panic for reverse-strand seeds near genome boundary (now `saturating_sub`). Also added guard in `finalize_transcript` to reject WTs where `adjusted_genome_start + ref_len > n_genome`. - -**Phase G1** (2026-04-29): `split_combined_wt` junction_idx fix. Reduced `.16980960`'s pairs from 11 to 9, matching STAR's NH=9. +## Testing Pattern -**Note on faithfulness change**: Phase F1 (--runRNGseed PR) changed PE tie-breaking from SA-order to seeded StdRng, increasing tie-breaking diffs. Phase G1 improved faithfulness from 99.755% → 99.865%. Phase G2: 99.865% → 99.883%. +- Unit tests: `#[cfg(test)]` in each module (the bulk of the suite) +- Integration tests: `tests/` directory (`alignment_features.rs`, `phase9_threading.rs`, `transcriptome_sam.rs`) — synthetic genomes, run the built binary via `assert_cmd` +- Differential testing against STAR: Python/shell harness in `test/` (NOT `tests/`) — `compare_sam.py`, `compare_pe.py`, `assess_faithfulness.py`, `compare_junctions.py`, `compare_chimeric.py`, etc. `test/ci.sh` drives them. +- Test data tiers: synthetic micro-genome → 10k yeast reads (ERR12389696) → full human genome +- **Don't hardcode a test count in prose** — it goes stale fast. Get the current number from `cargo test`. (`git grep -c '#\[test\]'` currently reports ~455 test fns.) -## Remaining Limitations (Top 2) +## Known Issues & Limitations -- No STARsolo single-cell features — Phase 14 (deferred) -- 4 PE AS diffs (rustar-aligner improvements, not bugs): `.844151` finds VIII:451791 0mm vs STAR's VII:1001391 6mm; `.4972950` finds correct spliced mate2 vs STAR's unspliced. Both cases: STAR's combined-window approach fails to stitch a PE pair at the better location. +Full issue tracking is in **[ROADMAP.md](ROADMAP.md)**; feature status in **[docs-old/phase17_features.md](docs-old/phase17_features.md)**. The durable takeaways a future session needs: -See [docs-old/phase17_features.md](docs-old/phase17_features.md) for full feature status. +- **Remaining SE/PE disagreements vs STAR are almost all genuine ties.** Both tools find the identical alignment set; the primary differs only because of SA-iteration order or RNG divergence (rustar-aligner tie-breaks with seeded `StdRng` via `--runRNGseed`, not STAR's mt19937). When a diff appears, **first check whether it is a tie** before treating it as a bug — the alignment sets matching is the signal. +- **One actionable SE CIGAR-only tie**: `ERR12389696.13573895` — same position (XV:218357), same AS=133, but insertion placed at read pos 100 vs STAR's 108, because the 71-base seed anchors at a different offset through a homopolymer run. Fixing it requires reproducing STAR's exact Lmapped chain. +- **4 PE AS diffs are rustar-aligner improvements, not bugs** (`.844151`, `.4972950`, …): rustar-aligner stitches a better-scoring pair than STAR's combined-window approach finds. Do not "fix" these toward STAR. +- **Not implemented**: STARsolo single-cell features (deferred). From b5b09dee53c5ba3f0a962306da8fb34fb6a79a1f Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Fri, 24 Jul 2026 02:03:40 +0200 Subject: [PATCH 5/5] feat: WASP allele-specific-mapping filter (paired-end) Complete the WASP filter for paired-end reads, the documented follow-up to the single-end port (3e84d97). Ports STAR-rs `wasp_type_pe` / `wasp_variants_pe` into rustar-aligner's two-transcript-per-pair model. - src/wasp/mod.rs: add wasp_variants_pe (combined vG/vA over mate1 then mate2, STAR's `mate1 ++ rc(mate2)` combined-transcript order), wasp_type_pe (single vW code: walk each mate's variants in its SAM frame, apply every allele combination, map back to forward mates, and re-map the pair via align_paired_read with the WASP-relaxed params; a combination passes only if exactly one both-mapped pair lands on the same exons for both mates), and annotate_records_pe (both mate records share vW/vA/vG; multimapped pairs overlapping variants get vW:i:2). - src/lib.rs: build the WaspContext once in align_reads_paired_end and stamp vW/vA/vG onto the built pair records, mirroring the single-end path. Gate green: cargo build, clippy -D warnings, fmt, 464 tests (2 new PE variant-frame unit tests). Byte-level validation against native STAR via the test/ harness is a separate step. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/lib.rs | 43 ++++++++- src/wasp/mod.rs | 251 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 291 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 73e1931..ee39cc0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1372,6 +1372,34 @@ fn align_reads_paired_end( let quant = quant_ctx.map(Arc::clone); let tr = tr_idx.map(Arc::clone); + // WASP allele-specific filtering context: load the VCF once, and build the + // relaxed re-map parameter set. `None` when --waspOutputMode is not SAMtag. + let wasp_ctx: Option = + if params.wasp_output_mode == params::WaspOutputMode::SAMtag { + let vcf = params + .var_vcf_file + .as_ref() + .expect("validated: SAMtag requires --varVCFfile"); + let ctx = crate::wasp::WaspContext::load( + vcf, + &index.genome.chr_name, + &index.genome.chr_start, + params, + ) + .map_err(|source| error::Error::Io { + source, + path: vcf.clone(), + })?; + info!( + "WASP: loaded {} heterozygous SNVs from {}", + ctx.snps.len(), + vcf.display() + ); + Some(ctx) + } else { + None + }; + info!( "Reading paired-end from {} and {}", params.read_files_in[0].display(), @@ -1695,7 +1723,7 @@ fn align_reads_paired_end( .iter() .map(|pa| PairedAlignment::clone(pa)) .collect(); - let records = SamWriter::build_paired_records( + let mut records = SamWriter::build_paired_records( &paired_read.name, &m1_seq, &m1_qual, @@ -1706,6 +1734,19 @@ fn align_reads_paired_end( params, n_for_mapq, )?; + // WASP allele-specific filtering: stamp vW/vA/vG on the pair. + if let Some(ctx) = wasp_ctx.as_ref() { + crate::wasp::annotate_records_pe( + &mut records, + &paired_alns, + &m1_seq, + &m2_seq, + &paired_read.name, + &index, + ctx, + params.out_sam_attributes, + )?; + } for record in records { buffer.push(record); } diff --git a/src/wasp/mod.rs b/src/wasp/mod.rs index 792477a..8b6ebe8 100644 --- a/src/wasp/mod.rs +++ b/src/wasp/mod.rs @@ -9,7 +9,11 @@ //! (`vW:i:1`), otherwise a failure code is reported. Reads overlapping no variant //! get no `vW` tag (code `-1`). //! -//! This revision covers **single-end**; paired-end WASP is a follow-up. +//! Covers **single-end** ([`annotate_records_se`]) and **paired-end** +//! ([`annotate_records_pe`]). For a pair, variants overlapping mate1 and mate2 +//! are walked in mate1-then-mate2 order (STAR's combined `mate1 ++ rc(mate2)` +//! transcript frame), a single `vW` code is computed by re-mapping the whole +//! pair with each allele combination, and both mate records share `vW`/`vA`/`vG`. //! //! Everything works in rustar-aligner's base-code space (`A=0,C=1,G=2,T=3,N=4`), //! and variant/read overlap is computed by walking the transcript CIGAR in @@ -24,7 +28,9 @@ use noodles::sam::alignment::record_buf::RecordBuf; use noodles::sam::alignment::record_buf::data::field::Value; use noodles::sam::alignment::record_buf::data::field::value::Array; -use crate::align::read_align::align_read; +use crate::align::read_align::{ + PairedAlignment, PairedAlignmentResult, align_paired_read, align_read, +}; use crate::align::transcript::Transcript; use crate::error::Error; use crate::index::GenomeIndex; @@ -359,6 +365,202 @@ pub fn annotate_records_se( Ok(()) } +/// STAR `Transcript::variationAdjust` for a paired-end alignment: the combined +/// `vG`/`vA` tag values from variants overlapping mate1 then mate2 (STAR's +/// `mate1 ++ rc(mate2)` combined-transcript order). Both mate records share this +/// array. `chr_start1`/`chr_start2` are the mates' chromosome genome offsets; +/// `m1_codes`/`m2_codes` are the forward mate reads (base codes). +#[allow(clippy::too_many_arguments)] +pub fn wasp_variants_pe( + chr_start1: u64, + chr_start2: u64, + snps: &[Snp], + m1_codes: &[u8], + m2_codes: &[u8], + tr1: &Transcript, + tr2: &Transcript, +) -> (Vec, Vec) { + let mut vg = Vec::new(); + let mut va = Vec::new(); + for (chr_start, codes, tr) in [(chr_start1, m1_codes, tr1), (chr_start2, m2_codes, tr2)] { + let sam_codes = if tr.is_reverse { + rc_codes(codes) + } else { + codes.to_vec() + }; + for (isnp, _rp, allele) in variation_overlap(snps, &sam_codes, tr) { + vg.push((snps[isnp].loci - chr_start) as i32); + va.push(allele); + } + } + (vg, va) +} + +/// STAR `waspMap` (paired-end): the `vW` code for a pair, or `-1` for no `vW` tag +/// (neither mate overlaps a variant). Codes match [`wasp_type`]: `1` passed, `2` +/// pair multimaps, `3` a variant base is `N`, `4` a re-map is not both-mapped, +/// `5` a re-map multimaps, `6` a re-map lands elsewhere, `7` too many variants. +/// Each allele combination is applied per mate in its SAM frame, mapped back to +/// forward mates, and re-mapped as a pair with the WASP-relaxed [`remap_params`]. +#[allow(clippy::too_many_arguments)] +pub fn wasp_type_pe( + index: &GenomeIndex, + snps: &[Snp], + m1_codes: &[u8], + m2_codes: &[u8], + read_name: &str, + tr1: &Transcript, + tr2: &Transcript, + n_tr: usize, + remap_params: &Parameters, +) -> Result { + let sam1 = if tr1.is_reverse { + rc_codes(m1_codes) + } else { + m1_codes.to_vec() + }; + let sam2 = if tr2.is_reverse { + rc_codes(m2_codes) + } else { + m2_codes.to_vec() + }; + let vars1 = variation_overlap(snps, &sam1, tr1); + let vars2 = variation_overlap(snps, &sam2, tr2); + if vars1.is_empty() && vars2.is_empty() { + return Ok(-1); + } + if n_tr > 1 { + return Ok(2); + } + let n1 = vars1.len(); + let n = n1 + vars2.len(); + if n > 10 { + return Ok(7); + } + if vars1 + .iter() + .chain(vars2.iter()) + .any(|&(_, _, allele)| allele > 3) + { + return Ok(3); + } + + // Combined actual alleles, mate1 then mate2 (matches STAR's combined read). + let actual: Vec = vars1 + .iter() + .chain(vars2.iter()) + .map(|&(_, _, a)| a) + .collect(); + for mask in 0..(1u32 << n) { + let combo: Vec = (0..n) + .map(|i| if mask & (1 << i) != 0 { 2 } else { 1 }) + .collect(); + if combo == actual { + continue; + } + // Flip each variant's base to the combination's allele in each mate's SAM + // frame, then map back to the forward mate reads. + let mut mod1 = sam1.clone(); + for (iv, &(isnp, read_pos, _)) in vars1.iter().enumerate() { + mod1[read_pos] = snps[isnp].nt[combo[iv] as usize]; + } + let mut mod2 = sam2.clone(); + for (iv, &(isnp, read_pos, _)) in vars2.iter().enumerate() { + mod2[read_pos] = snps[isnp].nt[combo[n1 + iv] as usize]; + } + let remap1 = if tr1.is_reverse { + rc_codes(&mod1) + } else { + mod1 + }; + let remap2 = if tr2.is_reverse { + rc_codes(&mod2) + } else { + mod2 + }; + + let (results, _chim, _n_for_mapq, _reason) = + align_paired_read(&remap1, &remap2, read_name, index, remap_params)?; + let both: Vec<&PairedAlignment> = results + .iter() + .filter_map(|r| match r { + PairedAlignmentResult::BothMapped(pa) => Some(pa.as_ref()), + PairedAlignmentResult::HalfMapped { .. } => None, + }) + .collect(); + if both.is_empty() { + return Ok(4); // re-map unmapped / half-mapped / too-many-loci + } + if both.len() > 1 { + return Ok(5); + } + let pa = both[0]; + if !same_exons(&pa.mate1_transcript, tr1) || !same_exons(&pa.mate2_transcript, tr2) { + return Ok(6); + } + } + Ok(1) +} + +/// Compute WASP tags for a paired-end read and stamp them onto its already-built +/// SAM records. `records` is `[mate1, mate2, mate1, mate2, ...]` (two per pair, as +/// [`crate::io::sam::SamWriter::build_paired_records`] emits); `paired[i]` +/// corresponds to `records[2i]`/`records[2i+1]`. A uniquely-mapped pair is +/// re-mapped (computing a single `vW` shared by both mate lines); a multi-mapped +/// pair overlapping variants gets `vW:i:2`. Pairs overlapping no variant are left +/// untagged. +#[allow(clippy::too_many_arguments)] +pub fn annotate_records_pe( + records: &mut [RecordBuf], + paired: &[PairedAlignment], + m1_codes: &[u8], + m2_codes: &[u8], + read_name: &str, + index: &GenomeIndex, + ctx: &WaspContext, + attrs: SamAttributes, +) -> Result<(), Error> { + let snps = &ctx.snps; + if paired.len() == 1 { + let pa = &paired[0]; + let (tr1, tr2) = (&pa.mate1_transcript, &pa.mate2_transcript); + let vw = wasp_type_pe( + index, + snps, + m1_codes, + m2_codes, + read_name, + tr1, + tr2, + 1, + &ctx.remap_params, + )?; + if vw != -1 { + let cs1 = index.genome.chr_start[tr1.chr_idx]; + let cs2 = index.genome.chr_start[tr2.chr_idx]; + let (vg, va) = wasp_variants_pe(cs1, cs2, snps, m1_codes, m2_codes, tr1, tr2); + for rec in records.iter_mut().take(2) { + insert_wasp_tags(rec, vw, &vg, &va, attrs); + } + } + } else { + for (i, pa) in paired.iter().enumerate() { + let (tr1, tr2) = (&pa.mate1_transcript, &pa.mate2_transcript); + let cs1 = index.genome.chr_start[tr1.chr_idx]; + let cs2 = index.genome.chr_start[tr2.chr_idx]; + let (vg, va) = wasp_variants_pe(cs1, cs2, snps, m1_codes, m2_codes, tr1, tr2); + if !vg.is_empty() { + for k in 0..2 { + if let Some(rec) = records.get_mut(2 * i + k) { + insert_wasp_tags(rec, 2, &vg, &va, attrs); + } + } + } + } + } + Ok(()) +} + #[cfg(test)] mod tests { use super::*; @@ -480,4 +682,49 @@ mod tests { let (vg, va) = wasp_variants(0, &snps, &[0u8; 50], &tr); assert!(vg.is_empty() && va.is_empty()); } + + #[test] + fn wasp_variants_pe_combines_mate1_then_mate2() { + // SNP over mate1 (genome 110) and mate2 (genome 310); both forward. + let snps = vec![ + Snp { + loci: 110, + nt: [0, 0, 2], + }, + Snp { + loci: 310, + nt: [0, 0, 2], + }, + ]; + let tr1 = tr_fwd_50m(100); + let tr2 = tr_fwd_50m(300); + let mut m1 = vec![0u8; 50]; + m1[10] = 2; // allele 2 at mate1 SNP + let mut m2 = vec![0u8; 50]; + m2[10] = 2; // allele 2 at mate2 SNP + let (vg, va) = wasp_variants_pe(0, 0, &snps, &m1, &m2, &tr1, &tr2); + // mate1 variant first, then mate2 (STAR's mate1 ++ rc(mate2) order). + assert_eq!(vg, vec![110, 310]); + assert_eq!(va, vec![2, 2]); + } + + #[test] + fn wasp_variants_pe_reverse_mate2_frame() { + // mate2 aligns reverse: its SAM frame is rc(mate2 read). + let snps = vec![Snp { + loci: 310, + nt: [0, 0, 2], + }]; + let tr1 = tr_fwd_50m(100); + let mut tr2 = tr_fwd_50m(300); + tr2.is_reverse = true; + let m1 = vec![0u8; 50]; // no mate1 overlap + // Build mate2 forward read whose rc has allele 2 at the SNP offset (10). + let mut sam2 = vec![0u8; 50]; + sam2[10] = 2; + let m2 = rc_codes(&sam2); + let (vg, va) = wasp_variants_pe(0, 0, &snps, &m1, &m2, &tr1, &tr2); + assert_eq!(vg, vec![310]); + assert_eq!(va, vec![2]); + } }