From 34b45d13e5a742381c9184ee8c48733622bce0c3 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 01:04:25 +0200 Subject: [PATCH 1/3] 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 4c0561b910937061cbd837d05e5be6e56c035a30 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 03:03:22 +0200 Subject: [PATCH 2/3] feat(genome): --genomeType SuperTranscriptome, and D20 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Condenses the genome to the union of its annotated exons, STAR's `GTF_superTranscript.cpp`. Overlapping exons merge, the merged intervals concatenate, overlapping transcripts group into superTranscripts, one per condensed chromosome. The index then covers only exonic sequence and introns cannot be crossed, because they are not there. The condensed genome comes back with its annotation already addressed to it — `st0`, one-based, plus strand — so the transcriptome tables, the junction extractor and the suffix array consume it unchanged. Nothing downstream needs to know the genome was condensed. `Genome::from_chromosomes` is split out of `from_fasta` for the same reason: these chromosomes are built, not read. Four files are written alongside the index, as STAR writes them: the superTranscript and transcript sequences, the collapsed junction table, and `fullGenome/conversionToFullGenome.tsv`, which is what makes condensed coordinates translatable back. D20, deliberately not reproduced. STAR condenses the sequence before it fills the reverse-complement half of its genome buffer, so every minus-strand exon reads the allocation's spacer bytes and minus-strand superTranscripts come out as uninitialised memory. Here that half is already filled, so a minus-strand exon yields its actual reverse complement. Reproducing STAR would mean writing code whose correct behaviour is to emit garbage, and a test that asserts it. The test asserts the hand-derived reverse complement instead, and that no 0x00 byte reaches the output. One behaviour worth stating because it looks like a missing feature: a transcript's own intron leaves no junction behind. With nothing else annotated between two exons, they end up contiguous in the condensed genome, and there is no gap left to record. A junction survives only when another transcript's exon separates them. Both cases have a test. Verified end to end: the 20 kb test genome with its two annotated 50 bp exons condenses to a single 100 bp `st0`; a read straddling the 200 bp intron aligns `50M` at position 26, and the conversion map names the two source intervals at 10000 and 10250. --- CHANGELOG.md | 30 ++ src/genome/mod.rs | 20 +- src/genome/supertranscript.rs | 518 ++++++++++++++++++++++++++++++++++ src/index/mod.rs | 78 ++++- src/params/mod.rs | 35 ++- tests/alignment_features.rs | 221 +++++++++++++++ 6 files changed, 884 insertions(+), 18 deletions(-) create mode 100644 src/genome/supertranscript.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..fa333bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,36 @@ 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. + +- **`--genomeType SuperTranscriptome`** condenses the genome to the union + of its annotated exons: overlapping exons are merged, the merged + intervals concatenated, and overlapping transcripts grouped into + superTranscripts, one per condensed chromosome (`st0`, `st1`, ...). + The index then covers only exonic sequence, and introns cannot be + crossed because they are not present. Requires `--sjdbGTFfile`. + Writes `superTranscriptSequences.fasta`, `transcriptSequences.fasta`, + `superTranscriptSJcollapsed.tsv` and + `fullGenome/conversionToFullGenome.tsv` alongside the index. + + Diverges from STAR on minus-strand exons, deliberately. STAR condenses + the sequence before filling the reverse-complement half of its genome + buffer, so minus-strand superTranscripts read uninitialised memory; + here they hold the actual reverse complement. Recorded in + `docs-old/dev/divergences.md` and locked by a test asserting the + hand-derived sequence. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/src/genome/mod.rs b/src/genome/mod.rs index 174ba42..fa7adaa 100644 --- a/src/genome/mod.rs +++ b/src/genome/mod.rs @@ -1,4 +1,5 @@ pub mod fasta; +pub mod supertranscript; pub mod transform; use std::path::Path; @@ -6,7 +7,7 @@ use std::path::Path; use crate::error::Error; use crate::params::Parameters; -use fasta::parse_fasta_files; +use fasta::{Chromosome, parse_fasta_files}; /// STAR's genome spacing character (used for inter-chromosome padding). const GENOME_SPACING_CHAR: u8 = 5; @@ -213,10 +214,25 @@ impl Genome { None }; + Self::from_chromosomes(&chromosomes, bin_nbits, transform_blocks) + } + + /// Lay out a genome from chromosome sequences already in base codes: + /// pad each to a `2^chr_bin_nbits` boundary, then fill the + /// reverse-complement half. + /// + /// [`from_fasta`](Self::from_fasta) is this plus FASTA parsing and the + /// VCF transform. `--genomeType SuperTranscriptome` uses it directly, + /// since its chromosomes are built rather than read. + pub fn from_chromosomes( + chromosomes: &[Chromosome], + bin_nbits: u32, + transform_blocks: Option>, + ) -> Result { // First pass: chromosome names/lengths, validating non-zero length. let mut chr_name = Vec::new(); let mut chr_length = Vec::new(); - for chrom in &chromosomes { + for chrom in chromosomes { let len = chrom.sequence.len() as u64; if len == 0 { return Err(Error::Fasta(format!( diff --git a/src/genome/supertranscript.rs b/src/genome/supertranscript.rs new file mode 100644 index 0000000..6ff9fef --- /dev/null +++ b/src/genome/supertranscript.rs @@ -0,0 +1,518 @@ +//! `--genomeType SuperTranscriptome`: condense a genome to the union of its +//! annotated exons (STAR `GTF_superTranscript.cpp` + +//! `SuperTranscriptome::sjCollapse`). +//! +//! Instead of indexing the whole genome, the index covers only exonic +//! sequence: overlapping exons are merged, the merged intervals are +//! concatenated, and the transcripts that overlap are grouped into +//! superTranscripts, one per condensed chromosome (`st0`, `st1`, ...). Reads +//! then align against a reference a fraction of the size, and introns cannot +//! be crossed because they are not there. +//! +//! The output is a genome and an annotation in condensed coordinates, which +//! the ordinary build path consumes unchanged, plus four files STAR writes for +//! downstream use: the superTranscript and transcript sequences, the collapsed +//! junction table, and the condensed-to-full coordinate map. +//! +//! # D20: minus-strand exons +//! +//! STAR condenses the superTranscript sequence *before* it fills the +//! reverse-complement half of its genome buffer, about twenty-five lines later +//! in `Genome_genomeGenerate.cpp`. Every minus-strand exon therefore reads the +//! allocation's spacer bytes, and minus-strand superTranscripts come out as +//! uninitialised sequence. +//! +//! This does not reproduce that. Here the reverse-complement half is already +//! filled when condensing runs, so a minus-strand exon yields its actual +//! reverse complement. The divergence is deliberate: STAR's output on that path +//! is not a rule to follow, it is a read of uninitialised memory. +//! `minus_strand_supertranscript_is_reverse_complement_not_spacer` locks the +//! correct answer, derived by hand. + +use std::collections::HashMap; + +use crate::error::Error; +use crate::genome::Genome; +use crate::genome::fasta::Chromosome; +use crate::junction::gtf::GtfRecord; + +/// The condensed genome and everything derived from it. +pub struct SuperTranscriptome { + /// Condensed chromosomes `st0..`, base codes, ready for `Genome::from_chromosomes`. + pub chromosomes: Vec, + /// The annotation in condensed coordinates: `seqname` is the + /// superTranscript, positions are 1-based within it, strand is always `+` + /// because a superTranscript has no orientation left to speak of. + pub exons: Vec, + /// `superTranscriptSequences.fasta`. + pub supertranscript_fasta: Vec, + /// `transcriptSequences.fasta`. + pub transcript_fasta: Vec, + /// `superTranscriptSJcollapsed.tsv`. + pub sj_collapsed_tsv: Vec, + /// `fullGenome/conversionToFullGenome.tsv`. + pub conversion_tsv: Vec, +} + +/// One exon in global doubled-genome coordinates, tagged with its transcript. +#[derive(Clone, Copy)] +struct Exon { + transcript: usize, + start: u64, + end: u64, + /// Index of the GTF record this came from, so the remapped annotation can + /// keep its attributes — gene id, gene name, gene type — rather than + /// dropping everything the transcriptome tables need. + source: usize, +} + +/// Condense `genome` to the union of the exons in `records`. +/// +/// `transcript_tag` is `--sjdbGTFtagExonParentTranscript`. Records naming a +/// chromosome the genome does not have, or carrying no transcript tag, are +/// skipped, matching what the junction extractor does with them. +pub fn condense( + genome: &Genome, + records: &[GtfRecord], + transcript_tag: &str, + chr_bin_nbits: u32, +) -> Result { + let chr_index: HashMap<&str, usize> = genome + .chr_name + .iter() + .enumerate() + .map(|(i, n)| (n.as_str(), i)) + .collect(); + + // Transcripts in first-appearance order, so a rebuild of the same GTF + // produces the same superTranscript numbering. + let mut transcript_ids: Vec = Vec::new(); + let mut transcript_of: HashMap = HashMap::new(); + let mut transcript_strand: Vec = Vec::new(); + let mut exons: Vec = Vec::new(); + + for (src, rec) in records.iter().enumerate() { + let Some(&ci) = chr_index.get(rec.seqname.as_str()) else { + continue; + }; + let Some(id) = rec.attributes.get(transcript_tag) else { + continue; + }; + let t = *transcript_of.entry(id.clone()).or_insert_with(|| { + transcript_ids.push(id.clone()); + transcript_strand.push(rec.strand); + transcript_ids.len() - 1 + }); + // GTF is 1-based inclusive; the genome is 0-based. + exons.push(Exon { + transcript: t, + start: genome.chr_start[ci] + rec.start - 1, + end: genome.chr_start[ci] + rec.end - 1, + source: src, + }); + } + + if exons.is_empty() { + return Err(Error::Gtf( + "--genomeType SuperTranscriptome needs at least one exon on a known chromosome" + .to_string(), + )); + } + + // (1) Reflect minus-strand exons into the reverse-complement half. See D20 + // in the module docs for why this reads real sequence here and does not in + // STAR. + let n = genome.n_genome; + for e in &mut exons { + if transcript_strand[e.transcript] == '-' { + let (s, t) = (e.start, e.end); + e.start = 2 * n - 1 - t; + e.end = 2 * n - 1 - s; + } + } + + // (2) Sort by start, (3) merge overlapping or touching exons and rewrite + // every exon into condensed coordinates as we go. + exons.sort_by_key(|e| e.start); + let mut merged: Vec<[u64; 2]> = Vec::new(); + let mut gap = exons[0].start; + let mut curr = [exons[0].start, exons[0].end]; + exons[0].start = 0; + exons[0].end -= gap; + for e in exons.iter_mut().skip(1) { + if e.start <= curr[1] + 1 { + curr[1] = curr[1].max(e.end); + } else { + gap += e.start - curr[1] - 1; + merged.push(curr); + curr = [e.start, e.end]; + } + e.start -= gap; + e.end -= gap; + } + merged.push(curr); + + // The condensed sequence: the merged intervals' bases, back to back. + let mut concat: Vec = Vec::new(); + for m in &merged { + for j in m[0]..=m[1] { + concat.push(genome.get_base(j).unwrap_or(4)); + } + } + + // (4) Each transcript's span in condensed coordinates. + let n_transcripts = transcript_ids.len(); + let mut tr_span = vec![[u64::MAX, 0u64]; n_transcripts]; + for e in &exons { + tr_span[e.transcript][0] = tr_span[e.transcript][0].min(e.start); + tr_span[e.transcript][1] = tr_span[e.transcript][1].max(e.end); + } + + // (5) Group overlapping transcript spans into superTranscripts, aligned to + // merged-interval boundaries so a superTranscript never starts mid-interval. + let mut merged_super = vec![0u64; merged.len()]; + let mut spans = tr_span.clone(); + spans.sort_by_key(|s| s[0]); + let mut super_span: Vec<[u64; 2]> = Vec::new(); + let mut mi = 0usize; + let mut mi_len = 0u64; + let mut curr = spans[0]; + for s in &spans { + while mi_len < curr[1] && mi < merged.len() { + mi_len += merged[mi][1] - merged[mi][0] + 1; + merged_super[mi] = super_span.len() as u64; + mi += 1; + } + curr[1] = curr[1].max(mi_len.saturating_sub(1)); + if s[0] <= curr[1] { + curr[1] = curr[1].max(s[1]); + } else { + super_span.push(curr); + curr = *s; + } + } + super_span.push(curr); + // Intervals past the last transcript span still belong somewhere. + for m in merged_super.iter_mut().skip(mi) { + *m = super_span.len() as u64 - 1; + } + + // (6) The superTranscript sequences. + let super_seq: Vec> = super_span + .iter() + .map(|s| concat[s[0] as usize..=s[1] as usize].to_vec()) + .collect(); + + // (7) Which superTranscript each transcript landed in. + let mut tr_super = vec![0u64; n_transcripts]; + { + let mut ist = 0usize; + for e in &exons { + while ist + 1 < super_span.len() && e.start > super_span[ist][1] { + ist += 1; + } + tr_super[e.transcript] = ist as u64; + } + } + + // (8) Transcript sequences: exons in transcript order, concatenated. + exons.sort_by_key(|e| (e.transcript, e.start)); + let mut transcript_seq: Vec> = vec![Vec::new(); n_transcripts]; + for e in &exons { + transcript_seq[e.transcript].extend_from_slice(&concat[e.start as usize..=e.end as usize]); + } + + let mut transcript_fasta = Vec::new(); + for (i, seq) in transcript_seq.iter().enumerate() { + transcript_fasta.push(b'>'); + transcript_fasta.extend_from_slice(transcript_ids[i].as_bytes()); + transcript_fasta.push(b'\n'); + transcript_fasta.extend(seq.iter().map(|&c| code_to_char(c))); + transcript_fasta.push(b'\n'); + } + + let mut supertranscript_fasta = Vec::new(); + for (i, seq) in super_seq.iter().enumerate() { + supertranscript_fasta.extend_from_slice(format!(">st{i}\n").as_bytes()); + supertranscript_fasta.extend(seq.iter().map(|&c| code_to_char(c))); + supertranscript_fasta.push(b'\n'); + } + + // (9) Junctions, in superTranscript-relative coordinates. + let mut sj: Vec<(u64, u64, u64)> = Vec::new(); // (superTranscript, start, end) + for i in 1..exons.len() { + if exons[i].transcript == exons[i - 1].transcript && exons[i].start > exons[i - 1].end + 1 { + let st = tr_super[exons[i].transcript]; + let base = super_span[st as usize][0]; + sj.push((st, exons[i - 1].end - base, exons[i].start - base)); + } + } + let sj_collapsed_tsv = collapse_sj(&mut sj); + + // (10) The condensed chromosomes, and the annotation rewritten onto them. + let chromosomes: Vec = super_seq + .iter() + .enumerate() + .map(|(i, s)| Chromosome { + name: format!("st{i}"), + sequence: s.clone(), + }) + .collect(); + + exons.sort_by_key(|e| e.start); + let mut out_records: Vec = Vec::with_capacity(exons.len()); + { + let mut ist = 0usize; + for e in &exons { + while ist + 1 < super_span.len() && e.start > super_span[ist][1] { + ist += 1; + } + let base = super_span[ist][0]; + let source = &records[e.source]; + out_records.push(GtfRecord { + seqname: format!("st{ist}"), + feature: source.feature.clone(), + attributes: source.attributes.clone(), + start: e.start - base + 1, // back to 1-based, superTranscript-local + end: e.end - base + 1, + strand: '+', + }); + } + } + + // (11) conversionToFullGenome.tsv: header, then one line per merged block. + let chr_start = super::compute_chr_starts( + &chromosomes + .iter() + .map(|c| c.sequence.len() as u64) + .collect::>(), + chr_bin_nbits, + ); + let mut conversion_tsv = Vec::new(); + conversion_tsv.extend_from_slice(format!("{}\t{}\n", merged.len(), n).as_bytes()); + let mut cond = 0u64; + for (ib, b) in merged.iter().enumerate() { + let len = b[1] - b[0] + 1; + let st = merged_super[ib] as usize; + let start_in_super = chr_start[st] + cond - super_span[st][0]; + conversion_tsv.extend_from_slice(format!("{start_in_super}\t{len}\t{}\n", b[0]).as_bytes()); + cond += len; + } + + Ok(SuperTranscriptome { + chromosomes, + exons: out_records, + supertranscript_fasta, + transcript_fasta, + sj_collapsed_tsv, + conversion_tsv, + }) +} + +/// Genome code to FASTA letter. Codes above 4 cannot occur here: every base +/// comes from a merged interval inside the real genome. +fn code_to_char(code: u8) -> u8 { + match code { + 0 => b'A', + 1 => b'C', + 2 => b'G', + 3 => b'T', + _ => b'N', + } +} + +/// STAR's `SuperTranscriptome::sjCollapse`: sort by (superTranscript, start, +/// end), drop exact duplicates, emit one line each. +fn collapse_sj(sj: &mut [(u64, u64, u64)]) -> Vec { + sj.sort_unstable(); + let mut out = Vec::new(); + let mut last: Option<(u64, u64, u64)> = None; + for &j in sj.iter() { + if last != Some(j) { + out.extend_from_slice(format!("{}\t{}\t{}\n", j.0, j.1, j.2).as_bytes()); + last = Some(j); + } + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn genome_of(seq: &str) -> Genome { + let codes: Vec = seq + .bytes() + .map(|b| match b { + b'A' => 0, + b'C' => 1, + b'G' => 2, + b'T' => 3, + _ => 4, + }) + .collect(); + Genome::from_chromosomes( + &[Chromosome { + name: "chr1".to_string(), + sequence: codes, + }], + 18, + None, + ) + .unwrap() + } + + fn exon(start: u64, end: u64, strand: char, transcript: &str) -> GtfRecord { + let mut attributes = HashMap::new(); + attributes.insert("transcript_id".to_string(), transcript.to_string()); + GtfRecord { + seqname: "chr1".to_string(), + feature: "exon".to_string(), + start, + end, + strand, + attributes, + } + } + + fn fasta_text(bytes: &[u8]) -> String { + String::from_utf8(bytes.to_vec()).unwrap() + } + + /// Two exons with an intron between them collapse to one superTranscript + /// holding only exonic sequence, and the intron is gone rather than spliced. + #[test] + fn exons_are_concatenated_and_the_intron_disappears() { + // 1234567890123456789012345678901234567890 + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + // Exons 1-8 and 25-32 (1-based), one transcript. + let records = vec![exon(1, 8, '+', "t1"), exon(25, 32, '+', "t1")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + + assert_eq!(st.chromosomes.len(), 1); + assert_eq!( + fasta_text(&st.supertranscript_fasta), + ">st0\nAAAACCCCAAAACCCC\n" + ); + assert_eq!(fasta_text(&st.transcript_fasta), ">t1\nAAAACCCCAAAACCCC\n"); + // No junction: with nothing else annotated between them, the two + // exons end up contiguous in the condensed genome, so there is no gap + // left for a read to cross. + assert_eq!(fasta_text(&st.sj_collapsed_tsv), ""); + } + + /// A junction survives condensing only when another transcript's exon sits + /// between the two, keeping them apart in the condensed genome. + #[test] + fn a_junction_is_recorded_when_another_exon_separates_the_two() { + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + let records = vec![ + exon(1, 8, '+', "t1"), // t1 exon 1 + exon(13, 20, '+', "t2"), // t2, between them + exon(25, 32, '+', "t1"), // t1 exon 2 + ]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + // Condensed: t1[0..7], t2[8..15], t1[16..23]. t1 jumps 7 -> 16. + assert_eq!(fasta_text(&st.sj_collapsed_tsv), "0\t7\t16\n"); + } + + /// The condensed annotation addresses the superTranscript, one-based, so + /// the ordinary GTF consumers need no special case. + #[test] + fn the_annotation_comes_back_in_condensed_coordinates() { + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + let records = vec![exon(1, 8, '+', "t1"), exon(25, 32, '+', "t1")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + + assert_eq!(st.exons.len(), 2); + assert_eq!(st.exons[0].seqname, "st0"); + assert_eq!((st.exons[0].start, st.exons[0].end), (1, 8)); + assert_eq!((st.exons[1].start, st.exons[1].end), (9, 16)); + assert!(st.exons.iter().all(|e| e.strand == '+')); + } + + /// Overlapping exons from different transcripts are merged once, not + /// duplicated, which is the whole point of condensing. + #[test] + fn overlapping_exons_are_merged_once() { + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + let records = vec![exon(1, 12, '+', "t1"), exon(5, 16, '+', "t2")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + assert_eq!( + fasta_text(&st.supertranscript_fasta), + ">st0\nAAAACCCCGGGGTTTT\n", + "the union of the two exons, each base once" + ); + } + + /// D20. STAR reads its unfilled reverse strand here and emits spacer bytes; + /// this emits the reverse complement, derived by hand below. + #[test] + fn minus_strand_supertranscript_is_reverse_complement_not_spacer() { + let seq = "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT"; + let genome = genome_of(seq); + // One minus-strand exon at 1-based 11..20 — forward offsets 10..=19. + let records = vec![exon(11, 20, '-', "tM")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + + let complement = |b: u8| match b { + b'A' => b'T', + b'C' => b'G', + b'G' => b'C', + b'T' => b'A', + x => x, + }; + let expected: String = seq.as_bytes()[10..20] + .iter() + .rev() + .map(|&b| complement(b) as char) + .collect(); + + assert_eq!( + fasta_text(&st.supertranscript_fasta), + format!(">st0\n{expected}\n"), + "a minus-strand exon must yield real sequence, not the spacer STAR reads" + ); + assert!( + !st.supertranscript_fasta.contains(&0x00), + "no uninitialised bytes may reach the output" + ); + } + + /// Two transcripts that share no sequence become two superTranscripts, so + /// the aligner cannot stitch across a boundary that does not exist. + #[test] + fn disjoint_transcripts_become_separate_supertranscripts() { + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + let records = vec![exon(1, 8, '+', "t1"), exon(33, 40, '+', "t2")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + assert_eq!(st.chromosomes.len(), 2); + assert_eq!(st.chromosomes[0].name, "st0"); + assert_eq!(st.chromosomes[1].name, "st1"); + } + + /// The conversion map is what makes condensed coordinates translatable + /// back, so its blocks must cover every merged interval. + #[test] + fn the_conversion_map_has_one_block_per_merged_interval() { + let genome = genome_of("AAAACCCCGGGGTTTTACGTACGTAAAACCCCGGGGTTTT"); + let records = vec![exon(1, 8, '+', "t1"), exon(25, 32, '+', "t1")]; + let st = condense(&genome, &records, "transcript_id", 18).unwrap(); + let text = fasta_text(&st.conversion_tsv); + let mut lines = text.lines(); + assert_eq!(lines.next().unwrap().split('\t').next().unwrap(), "2"); + let blocks: Vec<&str> = lines.collect(); + assert_eq!(blocks.len(), 2); + // condensed start, length, full-genome start + assert_eq!(blocks[0], "0\t8\t0"); + assert_eq!(blocks[1], "8\t8\t24"); + } + + #[test] + fn a_gtf_with_no_usable_exon_is_an_error() { + let genome = genome_of("AAAACCCC"); + let mut rec = exon(1, 4, '+', "t1"); + rec.seqname = "chrUnknown".to_string(); + assert!(condense(&genome, &[rec], "transcript_id", 18).is_err()); + } +} diff --git a/src/index/mod.rs b/src/index/mod.rs index 88ee62d..81a1c63 100644 --- a/src/index/mod.rs +++ b/src/index/mod.rs @@ -53,6 +53,8 @@ struct BuildPrep { junction_db: SpliceJunctionDb, transcriptome: Option, prepared_junctions: Vec, + /// The condensed-genome side files, under `--genomeType SuperTranscriptome`. + super_transcriptome: Option, } impl GenomeIndex { @@ -63,6 +65,7 @@ impl GenomeIndex { junction_db, transcriptome, prepared_junctions, + super_transcriptome: _, // side files are written by the streaming path } = Self::build_prep(params)?; log::info!("Building suffix array..."); @@ -131,6 +134,7 @@ impl GenomeIndex { junction_db: _, transcriptome, prepared_junctions, + super_transcriptome, } = Self::build_prep(params)?; let dir = ¶ms.genome_dir; @@ -141,6 +145,23 @@ impl GenomeIndex { log::info!("Writing genome files to {}...", dir.display()); genome.write_index_files(dir, params)?; + // The four condensed-genome files STAR writes alongside the index. + // `conversionToFullGenome.tsv` goes under `fullGenome/`, where STAR + // puts it, since it describes the genome this one was condensed from. + if let Some(ref st) = super_transcriptome { + let write = |name: &str, bytes: &[u8]| -> Result<(), Error> { + let path = dir.join(name); + fs::write(&path, bytes).map_err(|e| Error::io(e, &path)) + }; + write("superTranscriptSequences.fasta", &st.supertranscript_fasta)?; + write("transcriptSequences.fasta", &st.transcript_fasta)?; + write("superTranscriptSJcollapsed.tsv", &st.sj_collapsed_tsv)?; + let full = dir.join("fullGenome"); + fs::create_dir_all(&full).map_err(|e| Error::io(e, &full))?; + let path = full.join("conversionToFullGenome.tsv"); + fs::write(&path, &st.conversion_tsv).map_err(|e| Error::io(e, &path))?; + } + let gstrand_bit = SuffixArray::calculate_gstrand_bit(genome.n_genome); let gstrand_mask = (1u64 << gstrand_bit) - 1; let word_length = gstrand_bit + 1; @@ -272,6 +293,46 @@ impl GenomeIndex { genome.n_genome ); + // `--genomeType SuperTranscriptome`: replace the genome with the union + // of its annotated exons before anything else looks at it. Everything + // downstream — transcriptome tables, junctions, suffix array — then + // works on the condensed genome without knowing it is condensed, + // because the annotation comes back addressed to it. + let super_transcriptome = if params.genome_type == "SuperTranscriptome" { + let gtf_path = params.sjdb_gtf_file.as_ref().ok_or_else(|| { + Error::Index( + "--genomeType SuperTranscriptome requires --sjdbGTFfile: there is nothing \ + to condense without an annotation" + .to_string(), + ) + })?; + let records = crate::junction::gtf::parse_gtf_configured( + gtf_path, + ¶ms.sjdb_gtf_feature_exon, + ¶ms.sjdb_gtf_chr_prefix, + )?; + let st = crate::genome::supertranscript::condense( + &genome, + &records, + ¶ms.sjdb_gtf_tag_exon_parent_transcript, + params.genome_chr_bin_nbits, + )?; + log::info!( + "SuperTranscriptome: {} chromosomes condensed to {} superTranscripts", + genome.n_chr_real, + st.chromosomes.len() + ); + let full_size = genome.n_genome; + genome = Genome::from_chromosomes(&st.chromosomes, params.genome_chr_bin_nbits, None)?; + log::info!( + "Condensed genome: {} bytes, down from {full_size}", + genome.n_genome + ); + Some(st) + } else { + None + }; + // Parse GTF once and share the result between the junction database, // the transcriptome index, and the sjdb insertion pipeline. Junction // preparation + Gsj append must happen BEFORE the suffix array is @@ -283,11 +344,17 @@ impl GenomeIndex { let mut transcriptome = None; if let Some(ref gtf_path) = params.sjdb_gtf_file { - let exons = crate::junction::gtf::parse_gtf_configured( - gtf_path, - ¶ms.sjdb_gtf_feature_exon, - ¶ms.sjdb_gtf_chr_prefix, - )?; + // Under SuperTranscriptome the annotation was already rewritten + // into condensed coordinates; re-parsing the GTF here would + // address chromosomes the genome no longer has. + let exons = match &super_transcriptome { + Some(st) => st.exons.clone(), + None => crate::junction::gtf::parse_gtf_configured( + gtf_path, + ¶ms.sjdb_gtf_feature_exon, + ¶ms.sjdb_gtf_chr_prefix, + )?, + }; log::debug!("Parsed {} exon features from GTF", exons.len()); let tr = TranscriptomeIndex::from_gtf_exons_configured( @@ -375,6 +442,7 @@ impl GenomeIndex { junction_db, transcriptome, prepared_junctions, + super_transcriptome, }) } diff --git a/src/params/mod.rs b/src/params/mod.rs index 84a77cc..48259b1 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1577,17 +1577,29 @@ 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 - ), - )); + // --genomeType: `Transcriptome` needs a different index layout + // entirely, so it is refused rather than built as something that + // silently is not what was asked for. `SuperTranscriptome` condenses + // the genome to its annotated exons and therefore needs an annotation. + match params.genome_type.as_str() { + "Full" => {} + "SuperTranscriptome" => { + if params.sjdb_gtf_file.is_none() { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--genomeType SuperTranscriptome requires --sjdbGTFfile: there is \ + nothing to condense without an annotation", + )); + } + } + other => { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--genomeType {other} is not supported; expected Full or SuperTranscriptome" + ), + )); + } } // --genomeTransformOutput: mapping alignments back to original // coordinates is not implemented, so anything but None is refused. A @@ -2718,6 +2730,7 @@ mod tests { // The unimplemented genome layouts are refused rather than silently // building an ordinary genome under another name. + // SuperTranscriptome is built, but only with an annotation to condense. assert!(gg(&["--genomeType", "SuperTranscriptome"]).is_err()); assert!(gg(&["--genomeType", "Transcriptome"]).is_err()); diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..4b95764 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2024,3 +2024,224 @@ 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]); +} + +// --------------------------------------------------------------------------- +// --genomeType SuperTranscriptome +// --------------------------------------------------------------------------- + +/// The condensed index holds only exonic sequence, so a read spanning the +/// planted intron aligns without an `N` operation: in the condensed genome the +/// two exons are contiguous. +#[test] +fn test_genome_type_supertranscriptome_condenses_to_exons() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + + let genome_dir = tmpdir.path().join("genome_st"); + 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", + "5", + "--genomeType", + "SuperTranscriptome", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--sjdbOverhang", + "24", + ]) + .assert() + .success(); + + // The two annotated exons are the whole reference now. + let chr_names = fs::read_to_string(genome_dir.join("chrName.txt")).unwrap(); + assert_eq!(chr_names.trim(), "st0"); + let chr_lengths = fs::read_to_string(genome_dir.join("chrLength.txt")).unwrap(); + assert_eq!( + chr_lengths.trim(), + "100", + "50 bp of exon 1 plus 50 bp of exon 2, and none of the 200 bp intron" + ); + + // The side files STAR writes alongside the condensed index. + let super_fasta = + fs::read_to_string(genome_dir.join("superTranscriptSequences.fasta")).unwrap(); + assert!(super_fasta.starts_with(">st0\n")); + let expected: String = genome[10000..10050] + .iter() + .chain(genome[10250..10300].iter()) + .map(|&b| b as char) + .collect(); + assert_eq!(super_fasta.trim_end(), format!(">st0\n{expected}")); + + let conversion = + fs::read_to_string(genome_dir.join("fullGenome/conversionToFullGenome.tsv")).unwrap(); + let mut lines = conversion.lines(); + assert_eq!(lines.next().unwrap().split('\t').next().unwrap(), "2"); + assert_eq!(lines.next().unwrap(), "0\t50\t10000"); + assert_eq!(lines.next().unwrap(), "50\t50\t10250"); + + assert!(genome_dir.join("transcriptSequences.fasta").exists()); + assert!(genome_dir.join("superTranscriptSJcollapsed.tsv").exists()); + + // A read straddling the junction aligns end to end, no intron to cross. + let mut read = genome[10025..10050].to_vec(); + read.extend_from_slice(&genome[10250..10275]); + let fastq_path = tmpdir.path().join("reads.fq"); + { + let mut f = fs::File::create(&fastq_path).unwrap(); + writeln!(f, "@spanning").unwrap(); + f.write_all(&read).unwrap(); + writeln!(f).unwrap(); + writeln!(f, "+").unwrap(); + writeln!(f, "{}", "I".repeat(read.len())).unwrap(); + } + + let output_dir = tmpdir.path().join("out_st"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq_path.to_str().unwrap(), + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let sam = fs::read_to_string(output_dir.join("Aligned.out.sam")).unwrap(); + let record = sam + .lines() + .find(|l| !l.starts_with('@')) + .expect("the spanning read should align"); + let fields: Vec<&str> = record.split('\t').collect(); + assert_eq!(fields[2], "st0"); + assert_eq!(fields[3], "26", "25 bases into exon 1, 1-based"); + assert_eq!(fields[5], "50M", "contiguous in the condensed genome"); +} From 072078aaa3891c5aade2e6aac67bdc8687efcada Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 11:41:28 +0200 Subject: [PATCH 3/3] refactor: keep only the SuperTranscriptome theme in this PR CONTRIBUTING.md asks for one theme per PR and for the description to match the diff; the transform-output flag, its test and its changelog entry moved to their own PR, and the D20 divergence moved into the root DIVERGENCE.md. --- CHANGELOG.md | 13 ----- DIVERGENCE.md | 12 ++++ src/params/mod.rs | 63 +------------------- tests/alignment_features.rs | 113 ------------------------------------ 4 files changed, 13 insertions(+), 188 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fa333bc..6908cc4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,19 +97,6 @@ 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. - - **`--genomeType SuperTranscriptome`** condenses the genome to the union of its annotated exons: overlapping exons are merged, the merged intervals concatenated, and overlapping transcripts grouped into diff --git a/DIVERGENCE.md b/DIVERGENCE.md index bd957ed..cf4a545 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -30,6 +30,18 @@ This is the reason faithfulness is reported **tie-adjusted**. On the 10k yeast b **Source.** `src/rng.rs`, `src/align/read_align.rs` (`per_read_seed`, `shuffle_tied_prefix`), `src/params/mod.rs` (`MultimapperOrder`). STAR: `ReadAlign_multMapSelect.cpp`, `ReadAlignChunk` RNG seeding. +### 1.2 `--genomeType SuperTranscriptome` minus-strand sequence + +**What STAR does.** `GTF::superTranscript` builds the condensed sequence by reading the doubled genome buffer `G`, reflecting minus-strand exons into its second half. That half is not filled until roughly twenty-five lines later in `Genome_genomeGenerate.cpp`. Every minus-strand exon therefore reads the allocation's spacer bytes, and a minus-strand superTranscript comes out as uninitialised memory rather than sequence. + +**What rustar-aligner does.** The reverse-complement half is already filled when condensing runs, so a minus-strand exon yields its actual reverse complement. + +**Why.** STAR's output on this path is not a rule; it is a read of memory that has not been written. Reproducing it would mean writing code whose correct behaviour is to emit garbage, and a test asserting that garbage. + +**Impact.** Any annotation containing minus-strand transcripts — which is essentially all of them. A superTranscript built here holds real sequence where STAR's holds spacer bytes, so the two indices differ for those transcripts and the alignments against them are not comparable. + +**Source.** `src/genome/supertranscript.rs`, locked by `minus_strand_supertranscript_is_reverse_complement_not_spacer`, which asserts the hand-derived reverse complement of a non-palindromic exon and that no `0x00` byte reaches the output. STAR: `GTF_superTranscript.cpp`, `Genome_genomeGenerate.cpp`. + --- ## 2. Cases where rustar-aligner outperforms STAR diff --git a/src/params/mod.rs b/src/params/mod.rs index 48259b1..48422bd 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -513,27 +513,6 @@ pub struct Parameters { #[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 @@ -1601,31 +1580,6 @@ 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. - 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() { @@ -2714,7 +2668,7 @@ mod tests { } #[test] - fn genome_flags_accept_their_defaults_and_refuse_the_rest() { + fn genome_type_accepts_supertranscriptome_and_refuses_transcriptome() { let gg = |extra: &[&str]| { let mut a = vec!["--runMode", "genomeGenerate", "--genomeFastaFiles", "g.fa"]; a.extend_from_slice(extra); @@ -2724,26 +2678,11 @@ 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. // SuperTranscriptome is built, but only with an annotation to condense. 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()); } } diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 4b95764..56feb88 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2025,119 +2025,6 @@ fn test_wasp_samtag() { ); } -// --------------------------------------------------------------------------- -// --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]); -} - // --------------------------------------------------------------------------- // --genomeType SuperTranscriptome // ---------------------------------------------------------------------------