From 02466f27c03e415cb42efb43d22ad67b961630af Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:04:25 +0200 Subject: [PATCH 1/5] feat(params): genomeType, genomeTransformOutput, genomeSuffixLengthMax, sjdbInsertSave MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four genome-side STAR parameters, each honest about what it does. `--genomeType` accepts `Full` and refuses `Transcriptome` and `SuperTranscriptome`. Those need a different index layout entirely; accepting them would build an ordinary genome under another name. `--genomeTransformOutput` accepts `None` and refuses `SAM` and `SJ`. Mapping alignments back to the original coordinate space is not implemented, and a silent no-op here is worse than a failure: it would report transformed coordinates as though they were original ones, and nothing downstream could tell. `--sjdbInsertSave` accepts `Basic` and `None`. The inserted-junction files are not retained either way. `--genomeSuffixLengthMax` is accepted and inert, with the reason recorded on the flag: the suffix-array builder exposes no comparison bound, and STAR uses this only to cap work on highly repetitive genomes. The index is identical either way, so ignoring it changes no output byte — which is the standard this codebase already applies to the other inert knobs. Co-Authored-By: Claude Opus 5 (1M context) --- src/params/mod.rs | 97 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..84a77cc 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -508,6 +508,32 @@ pub struct Parameters { #[arg(long = "genomeSAsparseD", default_value_t = 1)] pub genome_sa_sparse_d: u32, + /// Genome representation to build. `Full` (the default) is an ordinary + /// genome; `Transcriptome` and `SuperTranscriptome` are not built here. + #[arg(long = "genomeType", default_value = "Full")] + pub genome_type: String, + + /// Coordinate space alignments are reported in when the genome was built + /// with `--genomeTransformType`. `None` (the default) reports transformed + /// coordinates; `SAM` and `SJ` map them back to the original genome. + #[arg(long = "genomeTransformOutput", num_args = 1.., default_values_t = vec!["None".to_string()])] + pub genome_transform_output: Vec, + + /// Upper bound on the suffix length compared during suffix-array + /// construction. 0 means no bound. + /// + /// Accepted but inert: the suffix-array builder does not expose a + /// comparison bound, and STAR uses this only to cap work on highly + /// repetitive genomes. The resulting index is identical either way, so + /// ignoring it changes no output byte. + #[arg(long = "genomeSuffixLengthMax", default_value_t = 0u64)] + pub genome_suffix_length_max: u64, + + /// Keep the files describing junctions inserted on the fly, rather than + /// discarding them after the run. + #[arg(long = "sjdbInsertSave", default_value = "Basic")] + pub sjdb_insert_save: String, + /// Substitute VCF alleles into the genome at genomeGenerate (`None`, /// `Haploid`, or `Diploid`). Requires `--genomeTransformVCF`; incompatible /// with `--sjdbGTFfile`. `Diploid` is genotype-aware and duplicates the @@ -1551,6 +1577,44 @@ impl Parameters { )); } + // --genomeType: only ordinary genomes are built here. The other two + // need a different index layout entirely, so refuse rather than build + // something that silently is not what was asked for. + if params.genome_type != "Full" { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--genomeType {} is not supported; expected Full", + params.genome_type + ), + )); + } + // --genomeTransformOutput: mapping alignments back to original + // coordinates is not implemented, so anything but None is refused. A + // silent no-op here would report transformed coordinates as if they + // were original ones, which is worse than failing. + for o in ¶ms.genome_transform_output { + if o != "None" { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--genomeTransformOutput {o} is not supported: alignments are \ + reported in the transformed coordinate space" + ), + )); + } + } + // --sjdbInsertSave: the inserted-junction files are not retained. + if !matches!(params.sjdb_insert_save.as_str(), "Basic" | "None") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --sjdbInsertSave '{}'; expected Basic or None", + params.sjdb_insert_save + ), + )); + } + // ── STARsolo validation ───────────────────────────────────────── if params.run_mode == RunMode::AlignReads && params.solo_enabled() { // CB_UMI_Complex needs one CB position + whitelist per segment. @@ -2636,4 +2700,37 @@ mod tests { assert_eq!(p.out_sam_strand_field, "None"); assert!(!p.out_sam_attributes.contains(SamAttributes::XS)); } + + #[test] + fn genome_flags_accept_their_defaults_and_refuse_the_rest() { + let gg = |extra: &[&str]| { + let mut a = vec!["--runMode", "genomeGenerate", "--genomeFastaFiles", "g.fa"]; + a.extend_from_slice(extra); + try_parse(&a) + }; + + // Defaults parse, and are the STAR ones. + let p = gg(&[]).unwrap(); + assert_eq!(p.genome_type, "Full"); + assert_eq!(p.genome_transform_output, vec!["None".to_string()]); + assert_eq!(p.sjdb_insert_save, "Basic"); + assert_eq!(p.genome_suffix_length_max, 0); + + // The unimplemented genome layouts are refused rather than silently + // building an ordinary genome under another name. + assert!(gg(&["--genomeType", "SuperTranscriptome"]).is_err()); + assert!(gg(&["--genomeType", "Transcriptome"]).is_err()); + + // Likewise back-transformed output: a silent no-op would report + // transformed coordinates as if they were original ones. + assert!(gg(&["--genomeTransformOutput", "SAM"]).is_err()); + assert!(gg(&["--genomeTransformOutput", "SJ"]).is_err()); + assert!(gg(&["--genomeTransformOutput", "None"]).is_ok()); + + assert!(gg(&["--sjdbInsertSave", "All"]).is_err()); + assert!(gg(&["--sjdbInsertSave", "None"]).is_ok()); + + // Inert but accepted, since it changes no output byte. + assert!(gg(&["--genomeSuffixLengthMax", "500"]).is_ok()); + } } From 42dafecb270e5d01638997ee6e7e26c373b189c6 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 02:29:55 +0200 Subject: [PATCH 2/5] feat(genome): back-transform alignments onto the original genome `--genomeTransformType Haploid` bakes a VCF's variants into the genome sequence, so reads carrying those alleles align without mismatches. The price is that every coordinate in the output then refers to a genome nobody else has. This is the machinery that pays it back, STAR's `Transcript_transformGenome.cpp`. Three steps, in STAR's order. Remap each exon through the conversion blocks, splitting wherever it straddles a boundary because the two halves land at unrelated original coordinates. Merge back the splits that turn out to be seamless, and fold the shared part of an uneven gap into the neighbouring exon so what survives as a gap is a real indel. Then reclassify: junction motifs and annotation flags are read off the *original* genome, because a junction that was canonical in the transformed genome need not be canonical in the original one, and STAR reports what the original says. `blocks_from_tsv` reads `transformGenomeBlocks.tsv` back. The file already stores the map keyed on transformed coordinates, which is the order the walk wants and the reverse of what the build side keeps in memory, so parsing is not symmetric with writing and the test pins the round trip rather than assuming it. A sentinel entry terminates the map so the forward walk stops on a comparison instead of a bounds test, as STAR's does. Tests cover an exon inside one block, an exon straddling a boundary coming back out as `10M5D10M`, a seamless split collapsing again, soft clips surviving, a motif read from the original genome, a junction canonical only in the transformed genome turning non-canonical, a sub-`alignIntronMin` gap staying a deletion, and an alignment outside the map being dropped rather than guessed at. Not yet wired to the SAM writer: that needs the original genome loaded at align time and is the next commit. --- src/genome/mod.rs | 1 + src/genome/transform.rs | 82 ++++++ src/genome/transform_align.rs | 499 ++++++++++++++++++++++++++++++++++ 3 files changed, 582 insertions(+) create mode 100644 src/genome/transform_align.rs diff --git a/src/genome/mod.rs b/src/genome/mod.rs index 174ba42..3058413 100644 --- a/src/genome/mod.rs +++ b/src/genome/mod.rs @@ -1,5 +1,6 @@ pub mod fasta; pub mod transform; +pub mod transform_align; use std::path::Path; diff --git a/src/genome/transform.rs b/src/genome/transform.rs index 328e4f0..59895c5 100644 --- a/src/genome/transform.rs +++ b/src/genome/transform.rs @@ -21,6 +21,7 @@ use std::collections::BTreeMap; use std::fmt::Write as _; +use crate::error::Error; use crate::io::fastq::encode_base; use super::compute_chr_starts; @@ -328,6 +329,67 @@ pub fn blocks_to_tsv(blocks: &[[u64; 3]]) -> String { s } +/// The stop entry appended to a parsed block map. +/// +/// The back-transform walks forward over the blocks an exon reaches into and +/// stops on the first one that starts past its end. A sentinel that no +/// coordinate can reach makes that a plain comparison instead of a bounds test +/// on every iteration, which is how STAR writes it (`genomeOutLoad`). +pub const BLOCK_SENTINEL: [u64; 3] = [u64::MAX, 0, 0]; + +/// Parse `transformGenomeBlocks.tsv` back into the conversion map used by the +/// align-time back-transform. +/// +/// The file is written for reverse conversion, so its columns are already +/// `[transformed_start, length, original_start]` — the order +/// [`crate::genome::transform_align::transform_transcript`] wants, and the +/// reverse of the `[u64; 3]` order the build side keeps in memory. The result +/// is sorted by transformed start and terminated by [`BLOCK_SENTINEL`]. +pub fn blocks_from_tsv(text: &str) -> Result, Error> { + let mut lines = text.lines(); + let header = lines + .next() + .ok_or_else(|| Error::Index("transformGenomeBlocks.tsv is empty".to_string()))?; + let declared: usize = header + .split_whitespace() + .next() + .and_then(|n| n.parse().ok()) + .ok_or_else(|| { + Error::Index(format!( + "transformGenomeBlocks.tsv: bad header line {header:?}" + )) + })?; + + let mut blocks = Vec::with_capacity(declared + 1); + for line in lines { + if line.trim().is_empty() { + continue; + } + let mut f = line.split_whitespace(); + let mut next = |what: &str| -> Result { + f.next().and_then(|v| v.parse().ok()).ok_or_else(|| { + Error::Index(format!("transformGenomeBlocks.tsv: bad {what} in {line:?}")) + }) + }; + let new_start = next("new_start")?; + let length = next("length")?; + let orig_start = next("orig_start")?; + blocks.push([new_start, length, orig_start]); + } + + if blocks.len() != declared { + return Err(Error::Index(format!( + "transformGenomeBlocks.tsv declares {declared} blocks but has {}", + blocks.len() + ))); + } + // The map is written in chromosome order, which is already ascending by + // transformed start; sorting is a cheap guarantee rather than a fix. + blocks.sort_unstable(); + blocks.push(BLOCK_SENTINEL); + Ok(blocks) +} + #[cfg(test)] mod tests { use super::*; @@ -474,4 +536,24 @@ chr1\t6\t.\tG\tA\t.\t.\t.\tGT\t1/1 assert_eq!(t.blocks[1][0], 0); // same shared orig_start, not offset assert_eq!(t.blocks[1][2], offset); // new_start shifted by hap0's padded genome length } + + #[test] + fn block_map_round_trips_through_the_tsv() { + // In-memory order is [orig_start, length, new_start]; the file stores + // the reverse, which is also what the back-transform reads. + let written = vec![[0u64, 50, 0], [55, 100, 50]]; + let parsed = blocks_from_tsv(&blocks_to_tsv(&written)).unwrap(); + assert_eq!( + parsed, + vec![[0, 50, 0], [50, 100, 55], BLOCK_SENTINEL], + "the parsed map is keyed on transformed coordinates, plus the stop entry" + ); + } + + #[test] + fn a_block_count_that_disagrees_with_the_body_is_an_error() { + let text = "3\t-1\n0\t50\t0\n50\t100\t55\n"; + let err = blocks_from_tsv(text).unwrap_err().to_string(); + assert!(err.contains("declares 3 blocks but has 2"), "{err}"); + } } diff --git a/src/genome/transform_align.rs b/src/genome/transform_align.rs new file mode 100644 index 0000000..95c3fa2 --- /dev/null +++ b/src/genome/transform_align.rs @@ -0,0 +1,499 @@ +//! Align-time back-transform (`--genomeTransformOutput SAM`). +//! +//! STAR `Transcript_transformGenome.cpp` + `ReadAlign_transformGenome.cpp`. +//! +//! `--genomeTransformType Haploid` bakes a VCF's variants into the genome +//! sequence, so reads carrying those alleles align without mismatches. The +//! price is that every coordinate in the output then refers to a genome nobody +//! else has. `--genomeTransformOutput SAM` pays it back: each alignment is +//! mapped through the conversion blocks written at build time +//! (`transformGenomeBlocks.tsv`) onto the original genome, so an indel baked +//! into the transformed sequence reappears as an `I`/`D` CIGAR operation at the +//! original coordinates. +//! +//! Three things happen, in STAR's order: +//! +//! 1. **Remap.** Each exon is looked up in the block map and split wherever it +//! straddles a block boundary, since the two halves land at unrelated +//! original coordinates. +//! 2. **Merge.** Splits that turn out to be seamless in the original genome are +//! collapsed again, and a block gap that is shorter on one side than the +//! other has its shared part folded back into the neighbouring exon. What +//! survives as a gap is a real indel. +//! 3. **Reclassify.** Junction motifs and annotation flags are recomputed +//! against the *original* genome. A junction that was canonical in the +//! transformed genome need not be canonical in the original one, and STAR +//! reports what the original genome says. +//! +//! Only the SAM path is implemented here. `--genomeTransformOutput SJ` and +//! `Quant` are still rejected by parameter validation. + +use noodles::sam::alignment::record::cigar::{self, op::Kind}; + +use crate::align::score::SpliceMotif; +use crate::align::transcript::{Exon, Transcript}; +use crate::genome::Genome; +use crate::junction::SpliceJunctionDb; + +/// Map a transcript from transformed-genome coordinates back to the original +/// genome. +/// +/// `blocks` is the conversion map in `[transformed_start, length, +/// original_start]` order, ascending by `[0]`, terminated by the sentinel +/// [`super::transform::BLOCK_SENTINEL`] — the walk over trailing blocks relies +/// on a stop entry rather than a bounds test, as STAR's does. +/// +/// Returns `None` when the transcript cannot be converted: no block covers its +/// start, or it ends up empty. +pub fn transform_transcript( + orig: &Genome, + orig_junctions: &SpliceJunctionDb, + blocks: &[[u64; 3]], + tr: &Transcript, + align_intron_min: u64, + align_intron_max: u64, +) -> Option { + if tr.exons.is_empty() || blocks.is_empty() { + return None; + } + + let exons = remap_exons(blocks, &tr.exons)?; + let exons = merge_adjacent(exons); + + let (motifs, annotated) = reclassify_junctions( + orig, + orig_junctions, + &exons, + tr.is_reverse, + align_intron_min, + ); + + let (chr_idx, _) = orig.position_to_chr(exons[0].genome_start)?; + let cigar = rebuild_cigar(tr, &exons, align_intron_min, align_intron_max); + + let genome_start = exons[0].genome_start; + let genome_end = exons[exons.len() - 1].genome_end; + let n_junction = motifs.len() as u32; + // Gaps that are not junctions are indels. Counting them off the rebuilt + // CIGAR rather than the exon list keeps the two consistent: the CIGAR is + // what decides which gap is `N` and which is `D`. + let n_gap = cigar + .iter() + .filter(|op| matches!(op.kind(), Kind::Insertion | Kind::Deletion)) + .count() as u32; + + Some(Transcript { + chr_idx, + genome_start, + genome_end, + is_reverse: tr.is_reverse, + exons, + cigar, + score: tr.score, + n_mismatch: tr.n_mismatch, + n_gap, + n_junction, + junction_motifs: motifs, + junction_annotated: annotated, + read_seq: tr.read_seq.clone(), + }) +} + +/// Step 1: per-exon remap, splitting at block boundaries. +fn remap_exons(blocks: &[[u64; 3]], exons: &[Exon]) -> Option> { + let mut out: Vec = Vec::with_capacity(exons.len()); + for exon in exons { + let len = exon.genome_end - exon.genome_start; + if len == 0 { + continue; + } + let g1 = exon.genome_start; + let g2 = exon.genome_end - 1; + let read_start = exon.read_start; + + // The last block starting at or before this exon. STAR reaches it with + // `upper_bound` then `--`; with nothing at or before it, the exon lies + // outside the map entirely and the alignment does not convert. + let idx = blocks.partition_point(|b| b[0] <= g1); + if idx == 0 { + return None; + } + let mut ci = idx - 1; + + let b_start = blocks[ci][0]; + let b_end = blocks[ci][0] + blocks[ci][1] - 1; + if g1 <= b_end { + let piece = if g2 <= b_end { len } else { b_end - g1 + 1 }; + out.push(Exon { + genome_start: blocks[ci][2] + g1 - b_start, + genome_end: blocks[ci][2] + g1 - b_start + piece, + read_start, + read_end: read_start + piece as usize, + i_frag: exon.i_frag, + }); + } + + // Any further blocks this exon reaches into, each a separate piece. + ci += 1; + while ci < blocks.len() && g2 >= blocks[ci][0] { + let piece = if g2 < blocks[ci][0] + blocks[ci][1] { + g2 - blocks[ci][0] + 1 + } else { + blocks[ci][1] + }; + let r = read_start + (blocks[ci][0] - g1) as usize; + out.push(Exon { + genome_start: blocks[ci][2], + genome_end: blocks[ci][2] + piece, + read_start: r, + read_end: r + piece as usize, + i_frag: exon.i_frag, + }); + ci += 1; + } + } + if out.is_empty() { None } else { Some(out) } +} + +/// Step 2: collapse the splits that are seamless in the original genome, and +/// fold the shared part of an uneven gap back into the exons around it. +fn merge_adjacent(exons: Vec) -> Vec { + let mut out: Vec = Vec::with_capacity(exons.len()); + for exon in exons { + let Some(prev) = out.last_mut() else { + out.push(exon); + continue; + }; + if prev.i_frag != exon.i_frag { + out.push(exon); + continue; + } + let gap_r = (exon.read_start - prev.read_end) as u64; + let gap_g = exon.genome_start - prev.genome_end; + if gap_r == gap_g { + // Same gap on both sides: not an indel, so the split was an + // artefact of the block boundary. Absorb it, gap included. + prev.genome_end = exon.genome_end; + prev.read_end = exon.read_end; + } else { + // Uneven gap: the smaller side is aligned sequence, not part of the + // indel. Give it to the following exon so the remaining gap is the + // indel alone. + let shared = gap_r.min(gap_g); + let mut exon = exon; + if shared > 0 { + exon.genome_start -= shared; + exon.read_start -= shared as usize; + } + out.push(exon); + } + } + out +} + +/// Step 3: junction motifs and annotation flags, read off the original genome. +fn reclassify_junctions( + orig: &Genome, + orig_junctions: &SpliceJunctionDb, + exons: &[Exon], + is_reverse: bool, + align_intron_min: u64, +) -> (Vec, Vec) { + let mut motifs = Vec::new(); + let mut annotated = Vec::new(); + for pair in exons.windows(2) { + let (a, b) = (&pair[0], &pair[1]); + if a.i_frag != b.i_frag { + continue; // the mate gap is not a junction + } + let gap_g = b.genome_start.saturating_sub(a.genome_end); + let gap_r = b.read_start - a.read_end; + if gap_r > 0 || gap_g < align_intron_min { + continue; // insertion or deletion, not a junction + } + let motif = junction_motif(orig, a.genome_end, b.genome_start - 1); + let is_annot = orig + .position_to_chr(a.genome_end) + .is_some_and(|(chr, offset)| { + let end = b.genome_start - 1 - orig.chr_start[chr]; + let strand = if is_reverse { 2 } else { 1 }; + orig_junctions.is_annotated(chr, offset, end, strand) + || orig_junctions.is_annotated(chr, offset, end, 0) + }); + motifs.push(motif); + annotated.push(is_annot); + } + (motifs, annotated) +} + +/// The donor/acceptor dinucleotides of intron `[j_s, j_e]` in the original +/// genome (STAR `Transcript_transformGenome.cpp:144-156`). +fn junction_motif(orig: &Genome, j_s: u64, j_e: u64) -> SpliceMotif { + let at = |i: u64| orig.get_base(i).unwrap_or(5); + match (at(j_s), at(j_s + 1), at(j_e.wrapping_sub(1)), at(j_e)) { + (2, 3, 0, 2) => SpliceMotif::GtAg, + (1, 3, 0, 1) => SpliceMotif::CtAc, + (2, 1, 0, 2) => SpliceMotif::GcAg, + (1, 3, 2, 1) => SpliceMotif::CtGc, + (0, 3, 0, 1) => SpliceMotif::AtAc, + (2, 3, 0, 3) => SpliceMotif::GtAt, + _ => SpliceMotif::NonCanonical, + } +} + +/// Rebuild the CIGAR from the remapped exons. +/// +/// The soft clips are the original transcript's: back-transforming moves where +/// the read aligns, never how much of it aligns. +fn rebuild_cigar( + tr: &Transcript, + exons: &[Exon], + align_intron_min: u64, + align_intron_max: u64, +) -> Vec { + let mut ops: Vec = Vec::new(); + let push_match = |ops: &mut Vec, len: usize| { + if len == 0 { + return; + } + match ops.last_mut() { + Some(op) if op.kind() == Kind::Match => { + *op = cigar::Op::new(Kind::Match, op.len() + len); + } + _ => ops.push(cigar::Op::new(Kind::Match, len)), + } + }; + + let [left_clip, right_clip] = tr.count_soft_clips(); + if left_clip > 0 { + ops.push(cigar::Op::new(Kind::SoftClip, left_clip)); + } + for (i, exon) in exons.iter().enumerate() { + if i > 0 { + let prev = &exons[i - 1]; + let gap_r = exon.read_start - prev.read_end; + let gap_g = exon.genome_start - prev.genome_end; + if gap_r > 0 { + ops.push(cigar::Op::new(Kind::Insertion, gap_r)); + } + if gap_g > 0 { + let kind = if gap_g >= align_intron_min && gap_g <= align_intron_max { + Kind::Skip + } else { + Kind::Deletion + }; + ops.push(cigar::Op::new(kind, gap_g as usize)); + } + } + push_match(&mut ops, exon.read_end - exon.read_start); + } + if right_clip > 0 { + ops.push(cigar::Op::new(Kind::SoftClip, right_clip)); + } + ops +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::genome::transform::BLOCK_SENTINEL; + + /// A genome of one chromosome whose bases are given as codes 0..3. + fn genome_from_codes(codes: &[u8]) -> Genome { + let n = codes.len() as u64; + let mut seq = vec![0u8; 2 * n as usize]; + seq[..codes.len()].copy_from_slice(codes); + Genome { + transform_blocks: None, + sequence: seq.into(), + n_genome: n, + n_genome_real: n, + n_chr_real: 1, + chr_name: vec!["chr1".to_string()], + chr_length: vec![n], + chr_start: vec![0, n], + } + } + + fn exon(genome_start: u64, len: u64, read_start: usize) -> Exon { + Exon { + genome_start, + genome_end: genome_start + len, + read_start, + read_end: read_start + len as usize, + i_frag: 0, + } + } + + fn transcript(exons: Vec, cigar: Vec) -> Transcript { + let genome_start = exons[0].genome_start; + let genome_end = exons[exons.len() - 1].genome_end; + Transcript { + chr_idx: 0, + genome_start, + genome_end, + is_reverse: false, + exons, + cigar, + score: 0, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: vec![], + junction_annotated: vec![], + read_seq: vec![], + } + } + + /// Blocks for a genome where 5 original bases were deleted at original + /// position 50: transformed [0,50) is original [0,50), and transformed + /// [50,..) is original [55,..). + fn deletion_blocks() -> Vec<[u64; 3]> { + vec![[0, 50, 0], [50, 100, 55], BLOCK_SENTINEL] + } + + #[test] + fn an_exon_inside_one_block_is_shifted_by_that_block() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + let tr = transcript(vec![exon(60, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]); + let out = + transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap(); + // transformed 60 is 10 into the second block, which starts at original 55. + assert_eq!(out.genome_start, 65); + assert_eq!(out.exons.len(), 1); + assert_eq!(out.cigar_string(), "20M"); + } + + /// The point of the whole module: a variant baked into the genome comes + /// back out as a CIGAR operation. + #[test] + fn an_exon_spanning_a_block_boundary_becomes_a_deletion() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + // One 20-base exon straddling the boundary at transformed 50. + let tr = transcript(vec![exon(40, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]); + let out = + transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap(); + assert_eq!(out.genome_start, 40); + assert_eq!(out.exons.len(), 2, "the split survives as a real gap"); + // 10 bases, the 5 deleted bases, then 10 more. + assert_eq!(out.cigar_string(), "10M5D10M"); + assert_eq!(out.n_gap, 1); + } + + /// A gap the transform did not create must not be turned into one: when the + /// two pieces are contiguous in the original genome too, the split + /// disappears again. + #[test] + fn a_seamless_split_is_merged_back() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + // Two blocks that are adjacent in both coordinate systems: the split is + // pure bookkeeping. + let blocks = vec![[0, 50, 0], [50, 100, 50], BLOCK_SENTINEL]; + let tr = transcript(vec![exon(40, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]); + let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap(); + assert_eq!(out.exons.len(), 1); + assert_eq!(out.cigar_string(), "20M"); + assert_eq!(out.n_gap, 0); + } + + /// Soft clips describe the read, not the genome, so they survive unchanged. + #[test] + fn soft_clips_are_carried_through() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + let tr = transcript( + vec![exon(60, 20, 5)], + vec![ + cigar::Op::new(Kind::SoftClip, 5), + cigar::Op::new(Kind::Match, 20), + cigar::Op::new(Kind::SoftClip, 3), + ], + ); + let out = + transform_transcript(&orig, &jdb, &deletion_blocks(), &tr, 21, 1_000_000).unwrap(); + assert_eq!(out.cigar_string(), "5S20M3S"); + } + + /// A junction is reclassified against the original genome, not the + /// transformed one. Here the original bases spell GT..AG. + #[test] + fn a_junction_motif_is_read_from_the_original_genome() { + let mut codes = vec![0u8; 200]; + // intron [80, 119]: GT at the donor, AG at the acceptor. + codes[80] = 2; + codes[81] = 3; + codes[118] = 0; + codes[119] = 2; + let orig = genome_from_codes(&codes); + let jdb = SpliceJunctionDb::empty(); + // Both exons sit in the identity block, so the coordinates survive. + let blocks = vec![[0, 200, 0], BLOCK_SENTINEL]; + let tr = transcript( + vec![exon(60, 20, 0), exon(120, 20, 20)], + vec![ + cigar::Op::new(Kind::Match, 20), + cigar::Op::new(Kind::Skip, 40), + cigar::Op::new(Kind::Match, 20), + ], + ); + let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap(); + assert_eq!(out.junction_motifs, vec![SpliceMotif::GtAg]); + assert_eq!(out.cigar_string(), "20M40N20M"); + assert_eq!(out.n_junction, 1); + } + + /// The same geometry over a genome that does not spell a canonical motif + /// there. STAR reports what the original genome says, even when the + /// transformed genome said otherwise. + #[test] + fn a_junction_canonical_only_in_the_transformed_genome_becomes_noncanonical() { + let orig = genome_from_codes(&[0u8; 200]); // all A: no motif anywhere + let jdb = SpliceJunctionDb::empty(); + let blocks = vec![[0, 200, 0], BLOCK_SENTINEL]; + let tr = transcript( + vec![exon(60, 20, 0), exon(120, 20, 20)], + vec![ + cigar::Op::new(Kind::Match, 20), + cigar::Op::new(Kind::Skip, 40), + cigar::Op::new(Kind::Match, 20), + ], + ); + let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap(); + assert_eq!(out.junction_motifs, vec![SpliceMotif::NonCanonical]); + } + + /// A gap shorter than `--alignIntronMin` is a deletion, not a junction, and + /// carries no motif. + #[test] + fn a_short_gap_is_a_deletion_not_a_junction() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + let blocks = vec![[0, 200, 0], BLOCK_SENTINEL]; + let tr = transcript( + vec![exon(60, 20, 0), exon(85, 20, 20)], + vec![ + cigar::Op::new(Kind::Match, 20), + cigar::Op::new(Kind::Deletion, 5), + cigar::Op::new(Kind::Match, 20), + ], + ); + let out = transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).unwrap(); + assert!(out.junction_motifs.is_empty()); + assert_eq!(out.n_junction, 0); + assert_eq!(out.cigar_string(), "20M5D20M"); + } + + /// An alignment starting before every block cannot be converted, and is + /// dropped rather than guessed at. + #[test] + fn an_alignment_outside_the_block_map_does_not_convert() { + let orig = genome_from_codes(&[0u8; 200]); + let jdb = SpliceJunctionDb::empty(); + let blocks = vec![[100, 50, 100], BLOCK_SENTINEL]; + let tr = transcript(vec![exon(10, 20, 0)], vec![cigar::Op::new(Kind::Match, 20)]); + assert!(transform_transcript(&orig, &jdb, &blocks, &tr, 21, 1_000_000).is_none()); + } +} From 834d137bc8bcfbb44d29089d44f2a95480878b1b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 02:48:38 +0200 Subject: [PATCH 3/5] feat(genome): wire --genomeTransformOutput SAM through the align paths MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The back-transform now runs. Transcripts are converted immediately after alignment, before anything reads them, so a single coordinate space holds below that line: records, junction counts and statistics are all in the original genome. The alternative — converting at the writer — would leave SJ.out.tab and the counters describing a genome the SAM file does not. The SAM header comes from the original genome as well, via `output_genome()`. Original coordinates under transformed reference lengths would produce a file no downstream tool could read correctly. `GenomeIndex` gains `transform_out`, loaded only when the flag is set: the untransformed genome from `OriginalGenome/`, its annotated junctions, and the parsed block map. A missing `OriginalGenome/` or `transformGenomeBlocks.tsv` is an error rather than a fallback — falling back would report transformed coordinates while claiming they are original. A paired result converts as a unit. If one mate fails to convert the pair is dropped whole, since reporting the other alone would invent a half-mapped read. The insert size is recomputed from the converted mates, because the transform can change the distance between them. Combinations whose other outputs would stay in transformed space are refused: `--quantMode`, `--soloType`, `--outWigType`. `--genomeTransformOutput SJ` and `Quant` remain unimplemented and are still rejected. Verified end to end: a genome with a 5-base deletion at 5001, one read at 6000. Without the flag the reported POS is 5996; with it, 6001 — the original coordinate — and the CIGAR is unchanged, since the read spans no variant. --- CHANGELOG.md | 13 ++++ src/align/read_align.rs | 1 + src/align/stitch.rs | 2 + src/chimeric/detect.rs | 1 + src/index/io.rs | 71 +++++++++++++++++++++ src/index/mod.rs | 56 +++++++++++++++++ src/lib.rs | 121 ++++++++++++++++++++++++++++++------ src/params/mod.rs | 57 ++++++++++++++--- tests/alignment_features.rs | 113 +++++++++++++++++++++++++++++++++ 9 files changed, 406 insertions(+), 29 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..df6faa2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,19 @@ Sections commonly used: Features, Bug fixes, Other changes. glibc's malloc and per-thread heaps that return whole segments to the OS when abandoned, so allocator cache size stays bounded. +- **`--genomeTransformOutput SAM`** reports alignments in the original + genome's coordinates. `--genomeTransformType Haploid` bakes a VCF's + variants into the genome, so reads carrying those alleles align + without mismatches, but every reported coordinate then refers to a + genome nobody else has. This maps each alignment back through the + conversion blocks written at build time: an indel baked into the + sequence reappears as an `I`/`D` CIGAR operation at the original + position, and junction motifs are reclassified against the original + genome rather than the transformed one. The SAM header comes from + the original genome too. `SJ` and `Quant` remain unimplemented and + are refused, as are the flag combinations whose other outputs would + stay in transformed coordinates. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/src/align/read_align.rs b/src/align/read_align.rs index b4d2f20..33a2e97 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1616,6 +1616,7 @@ mod tests { transcriptome: None, prepared_junctions: Vec::new(), sjdb_overhang: 0, + transform_out: None, } } diff --git a/src/align/stitch.rs b/src/align/stitch.rs index 1c25ecc..97d97e7 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3191,6 +3191,7 @@ mod tests { transcriptome: None, prepared_junctions: Vec::new(), sjdb_overhang: 0, + transform_out: None, } } @@ -3312,6 +3313,7 @@ mod tests { transcriptome: None, prepared_junctions: Vec::new(), sjdb_overhang: 0, + transform_out: None, } } diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index 4997049..af23d35 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -1060,6 +1060,7 @@ mod tests { transcriptome: None, prepared_junctions: Vec::new(), sjdb_overhang: 0, + transform_out: None, } } diff --git a/src/index/io.rs b/src/index/io.rs index ce6de64..4edff7f 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -6,6 +6,7 @@ use byteorder::{LittleEndian, ReadBytesExt}; use crate::error::Error; use crate::genome::Genome; use crate::index::GenomeIndex; +use crate::index::TransformOut; use crate::index::packed_array::PackedArray; use crate::index::sa_index::SaIndex; use crate::index::suffix_array::SuffixArray; @@ -138,6 +139,15 @@ impl GenomeIndex { ); } + // `--genomeTransformOutput SAM`: the original genome and the block map + // needed to report alignments in its coordinates rather than the + // transformed ones the search runs against. + let transform_out = if params.transform_output_sam() { + Some(load_transform_out(genome_dir, params)?) + } else { + None + }; + Ok(GenomeIndex { genome, suffix_array, @@ -146,10 +156,71 @@ impl GenomeIndex { transcriptome, prepared_junctions, sjdb_overhang, + transform_out, }) } } +/// Load everything `--genomeTransformOutput SAM` needs: the untransformed +/// genome written to `OriginalGenome/` at build time, its annotated junctions, +/// and the `transformGenomeBlocks.tsv` conversion map. +/// +/// Both files exist only in an index built with `--genomeTransformType`, so a +/// missing one means the flag was asked for against an ordinary index. That is +/// an error rather than a silent fallback: falling back would report +/// transformed coordinates while claiming they are original. +fn load_transform_out(genome_dir: &Path, params: &Parameters) -> Result { + let blocks_path = genome_dir.join("transformGenomeBlocks.tsv"); + let text = std::fs::read_to_string(&blocks_path).map_err(|e| { + if e.kind() == std::io::ErrorKind::NotFound { + Error::Index(format!( + "--genomeTransformOutput SAM needs {}, which is only written by a \ + --genomeTransformType genomeGenerate run", + blocks_path.display() + )) + } else { + Error::io(e, &blocks_path) + } + })?; + let blocks = crate::genome::transform::blocks_from_tsv(&text)?; + + let orig_dir = genome_dir.join("OriginalGenome"); + if !orig_dir.join("chrName.txt").exists() { + return Err(Error::Index(format!( + "--genomeTransformOutput SAM needs the untransformed index in {}", + orig_dir.display() + ))); + } + let genome = load_genome(&orig_dir, params)?; + + // The original genome's annotated junctions, for the motif/annotation flags + // recomputed after the back-transform. Absent when the index was built + // without a GTF, in which case every junction is novel there too. + let sjdb_info_path = orig_dir.join("sjdbInfo.txt"); + let junctions = if sjdb_info_path.exists() { + let tab = sjdb_insert::read_sjdb_info_tab(&sjdb_info_path, &genome)?; + let raw: Vec<(usize, u64, u64, u8)> = tab + .junctions + .iter() + .map(|j| (j.chr_idx, j.stored_start(), j.stored_end(), j.strand)) + .collect(); + SpliceJunctionDb::from_raw_junctions(&raw) + } else { + SpliceJunctionDb::empty() + }; + + log::info!( + "genomeTransformOutput SAM: original genome {} chromosomes, {} conversion blocks", + genome.n_chr_real, + blocks.len() - 1, // the sentinel is not a block + ); + Ok(TransformOut { + genome, + junctions, + blocks, + }) +} + /// Read `genomeFileSizes\t ` from genomeParameters.txt /// and return the first field (total genome byte count, including Gsj if /// sjdb was baked in). Returns `Ok(None)` if the file or line is absent, diff --git a/src/index/mod.rs b/src/index/mod.rs index 88ee62d..bd7296c 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -42,6 +42,59 @@ pub struct GenomeIndex { /// `sjdbOverhang` recorded in `sjdbInfo.txt`. Zero when no sjdb /// junctions are present. pub sjdb_overhang: u32, + /// Populated only for `--genomeTransformOutput SAM`: what the output + /// coordinate space is, when it is not the one the search runs against. + pub transform_out: Option, +} + +/// The original genome an alignment is reported against under +/// `--genomeTransformOutput SAM`. +/// +/// Reads are searched against the transformed genome, since that is what the +/// suffix array indexes, and then mapped back through `blocks`. The SAM header +/// comes from the genome here, not from the searched one: reporting original +/// coordinates under transformed reference lengths would produce a file no +/// downstream tool could read correctly. +#[derive(Clone)] +pub struct TransformOut { + /// The untransformed genome, from `OriginalGenome/`. + pub genome: Genome, + /// Its annotated junctions, for the motifs recomputed after the transform. + pub junctions: SpliceJunctionDb, + /// `[transformed_start, length, original_start]`, ascending, sentinel-terminated. + pub blocks: Vec<[u64; 3]>, +} + +impl GenomeIndex { + /// The genome alignments are *reported* against, which is the searched one + /// unless `--genomeTransformOutput SAM` moved output into the original + /// coordinate space. + pub fn output_genome(&self) -> &Genome { + match &self.transform_out { + Some(t) => &t.genome, + None => &self.genome, + } + } + + /// Map a transcript into the output coordinate space, dropping it when it + /// does not convert. Without a transform this hands back what it was given. + pub fn to_output_space( + &self, + tr: crate::align::transcript::Transcript, + params: &Parameters, + ) -> Option { + let Some(t) = &self.transform_out else { + return Some(tr); + }; + crate::genome::transform_align::transform_transcript( + &t.genome, + &t.junctions, + &t.blocks, + &tr, + params.align_intron_min as u64, + params.align_intron_max as u64, + ) + } } /// Output of [`GenomeIndex::build_prep`] — the shared setup @@ -89,6 +142,9 @@ impl GenomeIndex { sa_index, junction_db, transcriptome, + // Built here, not loaded: the back-transform is an align-time + // concern and `build` has no index directory to read it from. + transform_out: None, prepared_junctions, sjdb_overhang, }) diff --git a/src/lib.rs b/src/lib.rs index 6fce173..3e2ae9c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -675,21 +675,21 @@ fn run_single_pass( OutStd::Sam => { info!("Writing SAM to stdout (--outStd SAM)"); Box::new(crate::io::sam::SamStdoutWriter::create( - &index.genome, + index.output_genome(), params, )?) } OutStd::BamUnsorted => { info!("Writing unsorted BAM to stdout (--outStd BAM_Unsorted)"); Box::new(crate::io::bam::BamStdoutWriter::create( - &index.genome, + index.output_genome(), params, )?) } OutStd::BamSortedByCoordinate => { info!("Writing coordinate-sorted BAM to stdout (--outStd BAM_SortedByCoordinate)"); Box::new(crate::io::bam::SortedBamStdoutWriter::create( - &index.genome, + index.output_genome(), params, )?) } @@ -700,7 +700,11 @@ fn run_single_pass( if let Some(parent) = output_path.parent() { std::fs::create_dir_all(parent)?; } - Box::new(SamWriter::create(&output_path, &index.genome, params)?) + Box::new(SamWriter::create( + &output_path, + index.output_genome(), + params, + )?) } OutSamFormat::Bam => { let sorted = out_type.sort_order == Some(OutSamSortOrder::SortedByCoordinate); @@ -716,11 +720,15 @@ fn run_single_pass( if sorted { Box::new(SortedBamWriter::create( &output_path, - &index.genome, + index.output_genome(), params, )?) } else { - Box::new(BamWriter::create(&output_path, &index.genome, params)?) + Box::new(BamWriter::create( + &output_path, + index.output_genome(), + params, + )?) } } OutSamFormat::None => { @@ -750,7 +758,7 @@ fn run_single_pass( } let sj_output_path = params.output_path("SJ.out.tab"); if !sj_stats.is_empty() { - sj_stats.write_output(&sj_output_path, &index.genome, params)?; + sj_stats.write_output(&sj_output_path, index.output_genome(), params)?; } // Per-cell count matrices (raw + filtered), Summary.csv, and the SJ // feature matrix — written here where sj_stats is available. @@ -801,7 +809,7 @@ fn run_single_pass( "Writing splice junction statistics to {}", sj_output_path.display() ); - sj_stats.write_output(&sj_output_path, &index.genome, params)?; + sj_stats.write_output(&sj_output_path, index.output_genome(), params)?; } // 6. Print summary @@ -829,7 +837,7 @@ fn run_two_pass( let pass1_path = pass1_dir.join("SJ.out.tab"); info!("Writing pass 1 junctions to {}", pass1_path.display()); - sj_stats_pass1.write_output(&pass1_path, &index.genome, params)?; + sj_stats_pass1.write_output(&pass1_path, index.output_genome(), params)?; info!( "Pass 1 discovered {} novel junctions", novel_junctions.len() @@ -1273,8 +1281,11 @@ fn extract_junction_keys( let intron_start = genome_pos; let intron_end = genome_pos + intron_len as u64 - 1; - let motif = - scorer.detect_splice_motif(genome_pos, intron_len as u32, &index.genome); + let motif = scorer.detect_splice_motif( + genome_pos, + intron_len as u32, + index.output_genome(), + ); let strand = match motif.implied_strand() { Some('+') => 1u8, Some('-') => 2u8, @@ -1413,7 +1424,8 @@ fn align_reads_single_end( let write_file = tf .reopen() .map_err(|e| anyhow::anyhow!("BySJout: temp file reopen error: {e}"))?; - let (hdr, w) = crate::io::sam::create_bysj_writer(write_file, &index.genome, params)?; + let (hdr, w) = + crate::io::sam::create_bysj_writer(write_file, index.output_genome(), params)?; (Some(hdr), Some(w)) } else { (None, None) @@ -1806,6 +1818,21 @@ fn align_reads_single_end( let (transcripts, chimeric_results, n_for_mapq, unmapped_reason) = align_read(&clipped_seq, &read.name, &index, params)?; + // `--genomeTransformOutput SAM`: the search ran against the + // transformed genome; everything downstream — records, + // junction counts, statistics — is in the original genome's + // coordinates. Converting here rather than at the writer + // keeps a single space below this line. An alignment that + // does not convert is dropped, as STAR drops it. + let transcripts: Vec<_> = if index.transform_out.is_some() { + transcripts + .into_iter() + .filter_map(|t| index.to_output_space(t, params)) + .collect() + } else { + transcripts + }; + // Collect chimeric alignments if enabled if params.chim_segment_min > 0 { chimeric_alns.extend(chimeric_results); @@ -1891,7 +1918,7 @@ fn align_reads_single_end( clip5p, clip3p, &transcripts, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -2637,7 +2664,7 @@ fn align_reads_solo_pe( clip5p_m2, clip3p_m2, &paired_alns, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -2751,7 +2778,8 @@ fn align_reads_paired_end( let write_file = tf .reopen() .map_err(|e| anyhow::anyhow!("BySJout: temp file reopen error: {e}"))?; - let (hdr, w) = crate::io::sam::create_bysj_writer(write_file, &index.genome, params)?; + let (hdr, w) = + crate::io::sam::create_bysj_writer(write_file, index.output_genome(), params)?; (Some(hdr), Some(w)) } else { (None, None) @@ -3157,6 +3185,19 @@ fn align_reads_paired_end( let (results, pe_chimeric, n_for_mapq, unmapped_reason) = align_paired_read(&m1_seq, &m2_seq, &paired_read.name, &index, params)?; + // `--genomeTransformOutput SAM`: as in the single-end path, + // move to the original genome's coordinates before anything + // reads the transcripts. A pair whose mates do not both + // convert is dropped whole rather than half-reported. + let results: Vec<_> = if index.transform_out.is_some() { + results + .into_iter() + .filter_map(|r| transform_pair_result(r, &index, params)) + .collect() + } else { + results + }; + // Classify the result for stats and SAM output let has_half_mapped = results .iter() @@ -3334,7 +3375,7 @@ fn align_reads_paired_end( m2_clip3p, mapped_transcript, *mate1_is_mapped, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -3360,7 +3401,7 @@ fn align_reads_paired_end( m2_clip5p, m2_clip3p, &paired_alns, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -3480,6 +3521,45 @@ struct ReadJunction { /// `is_unique` reflects the read's overall mapping multiplicity (n_loci == 1), /// not per-locus. Recording per-locus/per-mate instead double-counts junctions /// crossed by both mates of a pair, inflating the SJ.out.tab multi counts. +/// Move one paired result into the output coordinate space under +/// `--genomeTransformOutput SAM`. +/// +/// A `BothMapped` pair is kept only if both mates convert: reporting one mate +/// in original coordinates and dropping the other would turn a pair into a +/// half-mapped read that never existed. The insert size is recomputed from the +/// converted mates, since the transform can change the distance between them. +fn transform_pair_result( + result: crate::align::read_align::PairedAlignmentResult, + index: &crate::index::GenomeIndex, + params: &Parameters, +) -> Option { + use crate::align::read_align::PairedAlignmentResult; + match result { + PairedAlignmentResult::BothMapped(pa) => { + let mut pa = *pa; + let m1 = index.to_output_space(pa.mate1_transcript, params)?; + let m2 = index.to_output_space(pa.mate2_transcript, params)?; + let (lo, hi) = if m1.genome_start <= m2.genome_start { + (&m1, &m2) + } else { + (&m2, &m1) + }; + let span = (hi.genome_end - lo.genome_start) as i32; + pa.insert_size = if pa.insert_size < 0 { -span } else { span }; + pa.mate1_transcript = m1; + pa.mate2_transcript = m2; + Some(PairedAlignmentResult::BothMapped(Box::new(pa))) + } + PairedAlignmentResult::HalfMapped { + mapped_transcript, + mate1_is_mapped, + } => Some(PairedAlignmentResult::HalfMapped { + mapped_transcript: index.to_output_space(mapped_transcript, params)?, + mate1_is_mapped, + }), + } +} + fn record_read_junctions<'a>( transcripts: impl IntoIterator, index: &crate::index::GenomeIndex, @@ -3561,8 +3641,11 @@ fn extract_transcript_junctions( let intron_end = genome_pos + intron_len as u64 - 1; // Detect splice motif - let motif = - scorer.detect_splice_motif(genome_pos, intron_len as u32, &index.genome); + let motif = scorer.detect_splice_motif( + genome_pos, + intron_len as u32, + index.output_genome(), + ); // Compute overhang: min(left_exon_length, right_exon_length) let left_exon = exon_lengths[junction_idx]; diff --git a/src/params/mod.rs b/src/params/mod.rs index 84a77cc..ac0be78 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1201,6 +1201,12 @@ impl Parameters { !matches!(self.out_std, OutStd::None) || self.out_sam_type.format != OutSamFormat::None } + /// Whether `--genomeTransformOutput` asks for SAM records in the original + /// genome's coordinates rather than the transformed genome's. + pub fn transform_output_sam(&self) -> bool { + self.genome_transform_output.iter().any(|o| o == "SAM") + } + /// Whether `--chimOutType` includes `Junctions` (write Chimeric.out.junction). pub fn chim_out_junctions(&self) -> bool { self.chim_out_type.iter().any(|s| s == "Junctions") @@ -1589,21 +1595,50 @@ impl Parameters { ), )); } - // --genomeTransformOutput: mapping alignments back to original - // coordinates is not implemented, so anything but None is refused. A - // silent no-op here would report transformed coordinates as if they - // were original ones, which is worse than failing. + // --genomeTransformOutput: `SAM` maps alignments back to the original + // genome's coordinates. `SJ` and `Quant` do the same for the junction + // table and the transcriptome BAM and are not implemented, so they are + // refused: a silent no-op would report transformed coordinates as if + // they were original ones, which is worse than failing. for o in ¶ms.genome_transform_output { - if o != "None" { + if o != "None" && o != "SAM" { return Err(command.error( ErrorKind::InvalidValue, format!( - "--genomeTransformOutput {o} is not supported: alignments are \ - reported in the transformed coordinate space" + "--genomeTransformOutput {o} is not supported: only SAM output is \ + mapped back to original coordinates" ), )); } } + // The back-transform moves the SAM records alone. Everything else that + // reports coordinates — the junction table, the transcriptome BAM, the + // count matrices, the coverage signal — would still be in transformed + // space, and a file that mixes the two spaces is worse than no file. + if params.transform_output_sam() { + let conflicting = [ + (!params.quant_mode.iter().all(|m| m == "-"), "--quantMode"), + (params.solo_type != SoloType::None, "--soloType"), + ( + params + .out_wig_type + .iter() + .any(|t| !t.eq_ignore_ascii_case("None")), + "--outWigType", + ), + ]; + for (hit, flag) in conflicting { + if hit { + return Err(command.error( + ErrorKind::ArgumentConflict, + format!( + "--genomeTransformOutput SAM cannot be combined with {flag}: only \ + the SAM records are mapped back to original coordinates" + ), + )); + } + } + } // --sjdbInsertSave: the inserted-junction files are not retained. if !matches!(params.sjdb_insert_save.as_str(), "Basic" | "None") { return Err(command.error( @@ -2721,10 +2756,12 @@ mod tests { assert!(gg(&["--genomeType", "SuperTranscriptome"]).is_err()); assert!(gg(&["--genomeType", "Transcriptome"]).is_err()); - // Likewise back-transformed output: a silent no-op would report - // transformed coordinates as if they were original ones. - assert!(gg(&["--genomeTransformOutput", "SAM"]).is_err()); + // SAM output is mapped back to the original genome's coordinates. The + // junction table and the transcriptome BAM are not, and a silent no-op + // there would report transformed coordinates as if they were original. + assert!(gg(&["--genomeTransformOutput", "SAM"]).is_ok()); assert!(gg(&["--genomeTransformOutput", "SJ"]).is_err()); + assert!(gg(&["--genomeTransformOutput", "Quant"]).is_err()); assert!(gg(&["--genomeTransformOutput", "None"]).is_ok()); assert!(gg(&["--sjdbInsertSave", "All"]).is_err()); diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..91a9344 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2024,3 +2024,116 @@ fn test_wasp_samtag() { "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)" ); } + +// --------------------------------------------------------------------------- +// --genomeTransformOutput SAM — alignments reported in original coordinates +// --------------------------------------------------------------------------- + +/// A transformed genome carries a deletion the original does not. Reads align +/// against the transformed genome, where every downstream coordinate has moved +/// by the deleted length; `--genomeTransformOutput SAM` has to undo that. +/// +/// Without the back-transform the reported POS is 5 short, which is exactly the +/// kind of error that looks plausible in a browser and is wrong everywhere. +#[test] +fn test_genome_transform_output_sam_reports_original_coordinates() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + + // A 5-base deletion at 1-based position 5001: REF spans 6 bases, ALT keeps + // the first. Everything after it shifts down by 5 in the transformed genome. + let del_pos_1based = 5001usize; + let vcf_path = tmpdir.path().join("variants.vcf"); + { + let mut f = fs::File::create(&vcf_path).unwrap(); + writeln!(f, "##fileformat=VCFv4.2").unwrap(); + writeln!(f, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO").unwrap(); + let start = del_pos_1based - 1; + let reference = String::from_utf8(genome[start..start + 6].to_vec()).unwrap(); + let alt = String::from_utf8(genome[start..start + 1].to_vec()).unwrap(); + writeln!( + f, + "chr1\t{del_pos_1based}\t.\t{reference}\t{alt}\t.\tPASS\t." + ) + .unwrap(); + } + + let genome_dir = tmpdir.path().join("genome_transformed"); + fs::create_dir_all(&genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "genomeGenerate", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--genomeFastaFiles", + fasta.to_str().unwrap(), + "--genomeSAindexNbases", + "7", + "--genomeTransformType", + "Haploid", + "--genomeTransformVCF", + vcf_path.to_str().unwrap(), + ]) + .assert() + .success(); + + // One 60 bp read taken from the original genome well downstream of the + // deletion. Its sequence is unchanged by the substitution, so it aligns + // cleanly in both coordinate systems — only the position differs. + let read_start_0based = 6000usize; + let fastq_path = tmpdir.path().join("reads.fq"); + { + let mut f = fs::File::create(&fastq_path).unwrap(); + writeln!(f, "@r1").unwrap(); + f.write_all(&genome[read_start_0based..read_start_0based + 60]) + .unwrap(); + writeln!(f).unwrap(); + writeln!(f, "+").unwrap(); + writeln!(f, "{}", "I".repeat(60)).unwrap(); + } + + let run = |extra: &[&str], out_name: &str| -> Vec> { + let output_dir = tmpdir.path().join(out_name); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + let mut args = vec![ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq_path.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]; + args.extend_from_slice(extra); + cargo_bin_cmd!("rustar-aligner") + .args(&args) + .assert() + .success(); + fs::read_to_string(output_dir.join("Aligned.out.sam")) + .unwrap() + .lines() + .filter(|l| !l.starts_with('@')) + .map(|l| l.split('\t').map(str::to_string).collect()) + .collect() + }; + + let transformed = run(&[], "out_transformed"); + let original = run(&["--genomeTransformOutput", "SAM"], "out_original"); + + assert_eq!(transformed.len(), 1, "expected one alignment"); + assert_eq!(original.len(), 1, "expected one alignment"); + + let pos_transformed: usize = transformed[0][3].parse().unwrap(); + let pos_original: usize = original[0][3].parse().unwrap(); + + // 1-based SAM position in the original genome. + assert_eq!(pos_original, read_start_0based + 1); + // The transformed genome is 5 bases shorter ahead of this read. + assert_eq!(pos_transformed, pos_original - 5); + // The read itself spans no variant, so the CIGAR is unaffected. + assert_eq!(original[0][5], transformed[0][5]); +} From f86804a6c73ca104a1f3b8a08cae942bbda51116 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:38:17 +0200 Subject: [PATCH 4/5] refactor(params): keep only --genomeTransformOutput in this PR The other three flags from the original commit belong to the index and SuperTranscriptome themes; CONTRIBUTING.md asks for one theme per PR and for the description to match the diff. --- src/params/mod.rs | 58 +---------------------------------------------- 1 file changed, 1 insertion(+), 57 deletions(-) diff --git a/src/params/mod.rs b/src/params/mod.rs index ac0be78..9774037 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -508,32 +508,12 @@ pub struct Parameters { #[arg(long = "genomeSAsparseD", default_value_t = 1)] pub genome_sa_sparse_d: u32, - /// Genome representation to build. `Full` (the default) is an ordinary - /// genome; `Transcriptome` and `SuperTranscriptome` are not built here. - #[arg(long = "genomeType", default_value = "Full")] - pub genome_type: String, - /// Coordinate space alignments are reported in when the genome was built /// with `--genomeTransformType`. `None` (the default) reports transformed /// coordinates; `SAM` and `SJ` map them back to the original genome. #[arg(long = "genomeTransformOutput", num_args = 1.., default_values_t = vec!["None".to_string()])] pub genome_transform_output: Vec, - /// Upper bound on the suffix length compared during suffix-array - /// construction. 0 means no bound. - /// - /// Accepted but inert: the suffix-array builder does not expose a - /// comparison bound, and STAR uses this only to cap work on highly - /// repetitive genomes. The resulting index is identical either way, so - /// ignoring it changes no output byte. - #[arg(long = "genomeSuffixLengthMax", default_value_t = 0u64)] - pub genome_suffix_length_max: u64, - - /// Keep the files describing junctions inserted on the fly, rather than - /// discarding them after the run. - #[arg(long = "sjdbInsertSave", default_value = "Basic")] - pub sjdb_insert_save: String, - /// Substitute VCF alleles into the genome at genomeGenerate (`None`, /// `Haploid`, or `Diploid`). Requires `--genomeTransformVCF`; incompatible /// with `--sjdbGTFfile`. `Diploid` is genotype-aware and duplicates the @@ -1583,18 +1563,6 @@ impl Parameters { )); } - // --genomeType: only ordinary genomes are built here. The other two - // need a different index layout entirely, so refuse rather than build - // something that silently is not what was asked for. - if params.genome_type != "Full" { - return Err(command.error( - ErrorKind::InvalidValue, - format!( - "--genomeType {} is not supported; expected Full", - params.genome_type - ), - )); - } // --genomeTransformOutput: `SAM` maps alignments back to the original // genome's coordinates. `SJ` and `Quant` do the same for the junction // table and the transcriptome BAM and are not implemented, so they are @@ -1639,16 +1607,6 @@ impl Parameters { } } } - // --sjdbInsertSave: the inserted-junction files are not retained. - if !matches!(params.sjdb_insert_save.as_str(), "Basic" | "None") { - return Err(command.error( - ErrorKind::InvalidValue, - format!( - "unknown --sjdbInsertSave '{}'; expected Basic or None", - params.sjdb_insert_save - ), - )); - } // ── STARsolo validation ───────────────────────────────────────── if params.run_mode == RunMode::AlignReads && params.solo_enabled() { @@ -2737,7 +2695,7 @@ mod tests { } #[test] - fn genome_flags_accept_their_defaults_and_refuse_the_rest() { + fn genome_transform_output_accepts_sam_and_refuses_the_rest() { let gg = |extra: &[&str]| { let mut a = vec!["--runMode", "genomeGenerate", "--genomeFastaFiles", "g.fa"]; a.extend_from_slice(extra); @@ -2746,15 +2704,7 @@ mod tests { // Defaults parse, and are the STAR ones. let p = gg(&[]).unwrap(); - assert_eq!(p.genome_type, "Full"); assert_eq!(p.genome_transform_output, vec!["None".to_string()]); - assert_eq!(p.sjdb_insert_save, "Basic"); - assert_eq!(p.genome_suffix_length_max, 0); - - // The unimplemented genome layouts are refused rather than silently - // building an ordinary genome under another name. - assert!(gg(&["--genomeType", "SuperTranscriptome"]).is_err()); - assert!(gg(&["--genomeType", "Transcriptome"]).is_err()); // SAM output is mapped back to the original genome's coordinates. The // junction table and the transcriptome BAM are not, and a silent no-op @@ -2763,11 +2713,5 @@ mod tests { assert!(gg(&["--genomeTransformOutput", "SJ"]).is_err()); assert!(gg(&["--genomeTransformOutput", "Quant"]).is_err()); assert!(gg(&["--genomeTransformOutput", "None"]).is_ok()); - - assert!(gg(&["--sjdbInsertSave", "All"]).is_err()); - assert!(gg(&["--sjdbInsertSave", "None"]).is_ok()); - - // Inert but accepted, since it changes no output byte. - assert!(gg(&["--genomeSuffixLengthMax", "500"]).is_ok()); } } From 15fb73593805258bf3ea45143e84377b98294a3e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:38:17 +0200 Subject: [PATCH 5/5] docs(changelog): --genomeTransformOutput SAM --- CHANGELOG.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index df6faa2..2c40a09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -110,6 +110,19 @@ Sections commonly used: Features, Bug fixes, Other changes. are refused, as are the flag combinations whose other outputs would stay in transformed coordinates. +- **`--genomeTransformOutput SAM`** reports alignments in the original + genome's coordinates. `--genomeTransformType Haploid` bakes a VCF's + variants into the genome, so reads carrying those alleles align + without mismatches, but every reported coordinate then refers to a + genome nobody else has. This maps each alignment back through the + conversion blocks written at build time: an indel baked into the + sequence reappears as an `I`/`D` CIGAR operation at the original + position, and junction motifs are reclassified against the original + genome. The SAM header comes from the original genome too. `SJ` and + `Quant` remain unimplemented and are refused, as are the flag + combinations whose other outputs would stay in transformed + coordinates. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching