From 92a406ec0c9d47e1f92be6e2e5648af8b84fc66a 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 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 d1fcb2e5299d480b4b7207e2b5e112fa68fa7af8 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 | 29 ++ src/genome/mod.rs | 20 +- src/genome/supertranscript.rs | 518 ++++++++++++++++++++++++++++++++++ src/index/mod.rs | 78 ++++- src/params/mod.rs | 35 ++- tests/alignment_features.rs | 108 +++++++ 6 files changed, 770 insertions(+), 18 deletions(-) create mode 100644 src/genome/supertranscript.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346..5066c80 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -169,6 +169,35 @@ 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. + +- **`--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 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 3d4ad11..8e23d2b 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1828,17 +1828,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 @@ -3080,6 +3092,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 86691d3..3b2552e 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2204,3 +2204,111 @@ fn test_read_name_separator_cuts_the_qname_and_is_configurable() { ); } } + +// --------------------------------------------------------------------------- +// --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 567e29b0e61ab49d1670b446f3cf758a0ab1044a 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. --- src/params/mod.rs | 62 -------------------------------------- tests/parameter_surface.rs | 1 - 2 files changed, 63 deletions(-) diff --git a/src/params/mod.rs b/src/params/mod.rs index 8e23d2b..f245cb7 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -519,27 +519,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 @@ -1852,21 +1831,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" - ), - )); - } - } // --readFilesType: only Fastx is supported. Reading pre-aligned SAM // input is separate work; refuse rather than silently treating a SAM // file as FASTQ. @@ -1883,17 +1847,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() { // CB_UMI_Complex needs one CB position + whitelist per segment. @@ -3086,26 +3039,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/parameter_surface.rs b/tests/parameter_surface.rs index df3dd94..4ae97c0 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -113,7 +113,6 @@ const NOT_YET_ACCEPTED: &[&str] = &[ // Genome index types and transforms. "genomeSuffixLengthMax", "genomeTransformOutput", - "genomeType", "sjdbInsertSave", // STARsolo barcode chemistry. "soloAdapterMismatchesNmax",