feat: WASP allele-specific-mapping filter (SE + PE)

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
+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.

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.

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<noodles ...::cigar::Op>` (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`.

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<W: AlignmentWriter + ?Sized>(
     };
     let mut bysj_meta: Vec<BySJReadMeta> = 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<crate::wasp::WaspContext> =
+        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<W: AlignmentWriter + ?Sized>(
                     }
                 } 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<W: AlignmentWriter + ?Sized>(
                         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);
                     }
@@ -1331,6 +1372,34 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
     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<crate::wasp::WaspContext> =
+        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(),
@@ -1654,7 +1723,7 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                         .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,
@@ -1665,6 +1734,19 @@ fn align_reads_paired_end<W: AlignmentWriter + ?Sized>(
                         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);
                     }

src/params/mod.rs

@@ -262,6 +262,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<Self, Self::Err> {
+        match s {
+            "None" => Ok(Self::None),
+            "SAMtag" => Ok(Self::SAMtag),
+            _ => Err(format!("unknown waspOutputMode value: '{s}'")),
+        }
+    }
+}
+
 // ---------------------------------------------------------------------------
 // Parameters struct
 // ---------------------------------------------------------------------------
@@ -661,6 +681,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<PathBuf>,
+
     // ── Chimeric ────────────────────────────────────────────────────────
     // ── Debug ───────────────────────────────────────────────────────
     /// Filter for debug logging: only log detailed alignment info for reads matching this name
@@ -962,6 +991,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)
     }
 

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.",
                 ),
         )
     }