From 0a4cd83ee954fc1f47888fd67d1dae321b7371b8 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:04:25 +0200 Subject: [PATCH 1/7] 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 3b50731..3d4ad11 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -514,6 +514,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 @@ -1802,6 +1828,33 @@ 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" + ), + )); + } + } // --readFilesType: only Fastx is supported. Reading pre-aligned SAM // input is separate work; refuse rather than silently treating a SAM // file as FASTQ. @@ -1818,6 +1871,17 @@ 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() { // CB_UMI_Complex needs one CB position + whitelist per segment. @@ -2998,4 +3062,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 7303d62f89f2e018428b727594837bbd501145cb Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 02:29:55 +0200 Subject: [PATCH 2/7] 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 9a10f629acbce92b2e3738a5493917f8d3c92799 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 02:48:38 +0200 Subject: [PATCH 3/7] 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 | 12 ++++ src/align/read_align.rs | 1 + src/align/stitch.rs | 2 + src/chimeric/detect.rs | 1 + src/genome/transform_align.rs | 1 - src/index/io.rs | 71 +++++++++++++++++++++ src/index/mod.rs | 56 +++++++++++++++++ src/lib.rs | 113 +++++++++++++++++++++++++++------ src/params/mod.rs | 58 +++++++++++++---- tests/alignment_features.rs | 114 ++++++++++++++++++++++++++++++++++ 10 files changed, 398 insertions(+), 31 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..109e025 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -169,6 +169,18 @@ Sections commonly used: Features, Bug fixes, Other changes. union: they came from one molecule. Writes `matrix.mtx`, `features.tsv` and `transcriptEndDistanceDistribution.txt` under `Solo.out/Transcript3p/raw/`. +- **`--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 diff --git a/src/align/read_align.rs b/src/align/read_align.rs index ec66a16..84354eb 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1611,6 +1611,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 8dcca5c..9c93af4 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -3212,6 +3212,7 @@ mod tests { transcriptome: None, prepared_junctions: Vec::new(), sjdb_overhang: 0, + transform_out: None, } } @@ -3333,6 +3334,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 a32af77..5bf403a 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/genome/transform_align.rs b/src/genome/transform_align.rs index 95c3fa2..9d1e1f6 100644 --- a/src/genome/transform_align.rs +++ b/src/genome/transform_align.rs @@ -95,7 +95,6 @@ pub fn transform_transcript( n_junction, junction_motifs: motifs, junction_annotated: annotated, - read_seq: tr.read_seq.clone(), }) } 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 5086fcb..ca442b5 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -712,21 +712,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, )?) } @@ -737,7 +737,7 @@ 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); @@ -753,11 +753,11 @@ 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 => { @@ -794,7 +794,7 @@ fn run_single_pass( let sj_output_path = params.output_path("SJ.out.tab"); // `--outSJtype None` suppresses SJ.out.tab entirely. if params.out_sj_type != "None" && !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. @@ -846,7 +846,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 @@ -874,7 +874,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() @@ -1318,8 +1318,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, @@ -1459,7 +1462,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) @@ -1852,6 +1856,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); @@ -1937,7 +1956,7 @@ fn align_reads_single_end( clip5p, clip3p, &transcripts, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -2683,7 +2702,7 @@ fn align_reads_solo_pe( clip5p_m2, clip3p_m2, &paired_alns, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -2798,7 +2817,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) @@ -3204,6 +3224,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() @@ -3381,7 +3414,7 @@ fn align_reads_paired_end( m2_clip3p, mapped_transcript, *mate1_is_mapped, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -3407,7 +3440,7 @@ fn align_reads_paired_end( m2_clip5p, m2_clip3p, &paired_alns, - &index.genome, + index.output_genome(), params, n_for_mapq, )?; @@ -3527,6 +3560,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, @@ -3608,8 +3680,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 3d4ad11..4de4dce 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1368,6 +1368,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") @@ -1840,17 +1846,18 @@ 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" ), )); } @@ -1870,7 +1877,34 @@ impl Parameters { )); } } - + // 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( @@ -3083,10 +3117,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 86691d3..aad8351 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2204,3 +2204,117 @@ fn test_read_name_separator_cuts_the_qname_and_is_configurable() { ); } } + + +// --------------------------------------------------------------------------- +// --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 f9b879bdffab6dd42f6de3087ffd1c461e340f6d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:38:17 +0200 Subject: [PATCH 4/7] 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 4de4dce..a48b443 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -514,32 +514,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 @@ -1834,18 +1814,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 @@ -1905,16 +1873,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() { @@ -3098,7 +3056,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); @@ -3107,15 +3065,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 @@ -3124,11 +3074,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 5a4332ab1083e08fc66af096fbdbe1060733a3c1 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:38:17 +0200 Subject: [PATCH 5/7] docs(changelog): --genomeTransformOutput SAM --- CHANGELOG.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 109e025..675f3b0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -182,6 +182,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 - `--runThreadN 1` ran on every logical core instead of on one. The From c54c691e2adf88c144f1b9a0f7e76110a5961699 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 26 Aug 2026 23:46:01 +0200 Subject: [PATCH 6/7] test(surface): --genomeTransformOutput is accepted now --- src/genome/transform_align.rs | 1 - tests/parameter_surface.rs | 1 - 2 files changed, 2 deletions(-) diff --git a/src/genome/transform_align.rs b/src/genome/transform_align.rs index 9d1e1f6..ec5084a 100644 --- a/src/genome/transform_align.rs +++ b/src/genome/transform_align.rs @@ -340,7 +340,6 @@ mod tests { n_junction: 0, junction_motifs: vec![], junction_annotated: vec![], - read_seq: vec![], } } diff --git a/tests/parameter_surface.rs b/tests/parameter_surface.rs index df3dd94..7e17dda 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -112,7 +112,6 @@ const NOT_YET_ACCEPTED: &[&str] = &[ "clip5pAdapterSeq", // Genome index types and transforms. "genomeSuffixLengthMax", - "genomeTransformOutput", "genomeType", "sjdbInsertSave", // STARsolo barcode chemistry. From b6449f34552105779e2bc45784d8d08a2a719514 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 00:23:40 +0200 Subject: [PATCH 7/7] style: rustfmt after the rebase --- src/lib.rs | 12 ++++++++++-- tests/alignment_features.rs | 1 - 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index ca442b5..33beea8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -737,7 +737,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.output_genome(), params)?) + Box::new(SamWriter::create( + &output_path, + index.output_genome(), + params, + )?) } OutSamFormat::Bam => { let sorted = out_type.sort_order == Some(OutSamSortOrder::SortedByCoordinate); @@ -757,7 +761,11 @@ fn run_single_pass( params, )?) } else { - Box::new(BamWriter::create(&output_path, index.output_genome(), params)?) + Box::new(BamWriter::create( + &output_path, + index.output_genome(), + params, + )?) } } OutSamFormat::None => { diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index aad8351..2f39071 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2205,7 +2205,6 @@ fn test_read_name_separator_cuts_the_qname_and_is_configurable() { } } - // --------------------------------------------------------------------------- // --genomeTransformOutput SAM — alignments reported in original coordinates // ---------------------------------------------------------------------------