From bac20af1b96efbe6aab312917e92d0fa2487e381 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 00:13:03 +0100 Subject: [PATCH] Solo: refuse what the read loops cannot do, check barcode and UMI input Options the solo loops have no code for, and more than one --soloUMIdedup, are refused instead of writing empty or zero files with exit 0. A whitelist of the wrong length is fatal. CB_UMI_Complex positions must parse and set the UMI width. --soloUMIlen above 32 is refused. soloCellFiltering decompresses features.tsv. The gene name tag reaches a GTF given at mapping. CellReads.stats is counted per feature. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 16 ++ DIVERGENCE.md | 10 ++ src/lib.rs | 2 +- src/params/mod.rs | 63 +++++++ src/quant/mod.rs | 18 +- src/quant/transcriptome.rs | 2 +- src/solo/count.rs | 14 +- src/solo/mod.rs | 136 ++++++++++---- src/solo/whitelist.rs | 38 +++- tests/solo_input_checks.rs | 358 +++++++++++++++++++++++++++++++++++++ 10 files changed, 609 insertions(+), 48 deletions(-) create mode 100644 tests/solo_input_checks.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ab8f0e..03c09e9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,22 @@ Sections commonly used: Features, Bug fixes, Other changes. ## [Unreleased] +### Bug fixes + +- Solo runs refuse options their read loops do not implement (`--quantMode + GeneCounts`/`TranscriptomeSAM`, `--quantGeneCoverage`, `--quantGeneSplicing`, + `--outReadsUnmapped Fastx`, `--outFilterType BySJout`, `--chimSegmentMin`, + `--waspOutputMode`, and more than one `--soloUMIdedup`) instead of writing + empty or zero output with exit 0. +- A whitelist whose barcodes are not as long as the barcode read gives is now + fatal, as in STAR. It used to give an empty matrix, or false matches. +- `CB_UMI_Complex`: a missing or bad position is fatal, and the UMI width comes + from `--soloUMIposition` rather than `--soloUMIlen`. +- `--soloUMIlen` above 32 is refused rather than corrupting UMIs. +- `--runMode soloCellFiltering` no longer writes gzip bytes into `features.tsv`. +- `--sjdbGTFtagExonParentGeneName` applies when the GTF is given at mapping. +- `CellReads.stats` counts each feature's reads in that feature's file. + ### Other changes - Third-party notices completed. `LICENSES/` gains the licences of code diff --git a/DIVERGENCE.md b/DIVERGENCE.md index f0890fb..e9326cb 100644 --- a/DIVERGENCE.md +++ b/DIVERGENCE.md @@ -171,6 +171,16 @@ A rustar-aligner extension with no STAR equivalent, off by default (`No`). With **Source.** `src/solo/count.rs` (`write_cellranger_summary`), `src/solo/gene.rs` (`Region`). +### 1.10 Solo options that are refused rather than implemented + +**What STAR does.** STARsolo runs accept `--quantMode GeneCounts TranscriptomeSAM`, `--outReadsUnmapped Fastx`, `--outFilterType BySJout`, `--chimSegmentMin`, `--waspOutputMode`, and several `--soloUMIdedup` methods (one matrix per method). These descriptions come from STAR's documentation; solo runs of the Homebrew binary could not be compared on macOS, where it stops at `geneInfo.tab`. + +**What rustar-aligner does.** With `--soloType` it exits with an error naming the option. The same goes for this tool's own `--quantGeneCoverage` and `--quantGeneSplicing`. + +**Why.** The solo read loops have no code for these. They were accepted and wrote empty or all-zero files with exit 0, or dropped all methods but the first. An error is better than a result that looks real. Each one can be lifted when the loop implements it. + +**Source.** `src/params/mod.rs` (STARsolo validation). + --- ## 2. Cases where rustar-aligner outperforms STAR diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..d6b2135 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1051,7 +1051,7 @@ fn run_single_pass( } if let Some(info) = sctx.tag_read_info() { let info = info.lock().unwrap(); - let umi_len = params.solo_umi_len as usize; + let umi_len = sctx.umi_len; writer.finish_with(&mut |record| { crate::io::sam::apply_solo_read_info(record, &info, &sctx.whitelist, umi_len); })?; diff --git a/src/params/mod.rs b/src/params/mod.rs index 47e4ecf..40d592d 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -2311,8 +2311,57 @@ impl Parameters { // ── STARsolo validation ───────────────────────────────────────── if params.run_mode() == RunMode::AlignReads && params.solo_enabled() { + // The solo read loops have no code for these, so they would write + // empty or all-zero files and exit 0 (DIVERGENCE §1.10). + let unsupported = [ + (params.quant_gene_counts(), "--quantMode GeneCounts"), + ( + params.quant_transcriptome_sam(), + "--quantMode TranscriptomeSAM", + ), + (params.quant_gene_coverage(), "--quantGeneCoverage"), + (params.quant_gene_splicing(), "--quantGeneSplicing"), + ( + params.out_reads_unmapped == OutReadsUnmapped::Fastx, + "--outReadsUnmapped Fastx", + ), + ( + params.out_filter_type == OutFilterType::BySJout, + "--outFilterType BySJout", + ), + (params.chim_segment_min > 0, "--chimSegmentMin"), + ( + params.wasp_output_mode != WaspOutputMode::None, + "--waspOutputMode", + ), + ]; + if let Some((_, flag)) = unsupported.iter().find(|(on, _)| *on) { + return Err(command.error( + ErrorKind::ArgumentConflict, + format!("{flag} is not implemented with --soloType, so it cannot be used here"), + )); + } + // Only the first method reaches the matrix (DIVERGENCE §1.10). + if params.solo_umi_dedup.len() > 1 { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--soloUMIdedup takes one method here, got {}", + params.solo_umi_dedup.len() + ), + )); + } // CB_UMI_Complex needs one CB position + whitelist per segment. if params.solo_type == SoloType::CbUmiComplex { + if params.solo_umi_position.is_empty() { + return Err(command.error( + ErrorKind::MissingRequiredArgument, + "--soloType CB_UMI_Complex requires --soloUMIposition", + )); + } + if let Err(e) = crate::solo::check_positions(¶ms) { + return Err(command.error(ErrorKind::InvalidValue, e.to_string())); + } if params.solo_cb_position.is_empty() { return Err(command.error( ErrorKind::MissingRequiredArgument, @@ -2537,6 +2586,20 @@ impl Parameters { "--soloCBlen and --soloUMIlen must be > 0 for soloType CB_UMI_Simple", )); } + // UMIs are packed in a u64 too. CB_UMI_Complex takes its width from + // --soloUMIposition, checked above. + if params.solo_type != SoloType::CbUmiComplex + && params.solo_umi_len as usize > crate::solo::whitelist::CB_LEN_MAX + { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--soloUMIlen {} exceeds the maximum of {}", + params.solo_umi_len, + crate::solo::whitelist::CB_LEN_MAX + ), + )); + } // Cell barcode cannot exceed a u64 packing (32 bases). if params.solo_cb_len as usize > crate::solo::whitelist::CB_LEN_MAX { return Err(command.error( diff --git a/src/quant/mod.rs b/src/quant/mod.rs index f3c39f6..db955d0 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -69,7 +69,12 @@ impl GeneAnnotation { /// Build from GTF exon records using `gene_tag` as the gene grouping attribute. /// /// `gene_tag` is STAR's `sjdbGTFtagExonParentGene` (default `"gene_id"`). - pub fn from_gtf_exons_configured(exons: &[GtfRecord], genome: &Genome, gene_tag: &str) -> Self { + pub fn from_gtf_exons_configured( + exons: &[GtfRecord], + genome: &Genome, + gene_tag: &str, + gene_name_tags: &[String], + ) -> Self { let n_chrs = genome.n_chr_real; let mut gene_ids: Vec = Vec::new(); let mut gene_names: Vec = Vec::new(); @@ -96,11 +101,9 @@ impl GeneAnnotation { gene_is_reverse.push(is_rev); // gene_name (symbol) for features.tsv col 2; STAR falls back to // the gene_id when the GTF record omits gene_name. - let gene_name = exon - .attributes - .get("gene_name") - .cloned() - .unwrap_or_else(|| gene_id.clone()); + let gene_name = + crate::quant::transcriptome::find_attr_list(&exon.attributes, gene_name_tags) + .unwrap_or_else(|| gene_id.clone()); gene_names.push(gene_name); gene_ids.push(gene_id); gene_span.push(None); @@ -290,7 +293,7 @@ impl GeneAnnotation { /// Build from GTF exon records using default `"gene_id"` attribute (backward-compatible). pub fn from_gtf_exons(exons: &[GtfRecord], genome: &Genome) -> Self { - Self::from_gtf_exons_configured(exons, genome, "gene_id") + Self::from_gtf_exons_configured(exons, genome, "gene_id", &["gene_name".to_string()]) } pub fn n_genes(&self) -> usize { @@ -629,6 +632,7 @@ pub fn resolve_gene_annotation( &exons, genome, ¶ms.sjdb_gtf_tag_exon_parent_gene, + ¶ms.sjdb_gtf_tag_exon_parent_gene_name, )); } match transcriptome { diff --git a/src/quant/transcriptome.rs b/src/quant/transcriptome.rs index 0a02059..d9bc526 100644 --- a/src/quant/transcriptome.rs +++ b/src/quant/transcriptome.rs @@ -31,7 +31,7 @@ const MISSING_GENE_TYPE: &str = "MissingGeneType"; /// STAR's `find_attr_list` (`GTF.cpp`): the last of `tags`, in list order, that has an entry in /// `attrs`. `--sjdbGTFtagExonParentGeneName`/`GeneType` accept several candidate keys; when more /// than one matches, the last one in the list wins. -fn find_attr_list(attrs: &HashMap, tags: &[String]) -> Option { +pub(crate) fn find_attr_list(attrs: &HashMap, tags: &[String]) -> Option { tags.iter() .filter_map(|t| attrs.get(t)) .next_back() diff --git a/src/solo/count.rs b/src/solo/count.rs index 2acfb8a..a408667 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1534,7 +1534,7 @@ pub fn write_gene_matrix( } else { 0.0 }; - let umi_len = params.solo_umi_len as usize; + let umi_len = ctx.umi_len; let solo_dir = params .solo_out_file_names @@ -1654,7 +1654,7 @@ pub fn write_gene_matrix( // `--soloCellReadStats CB`: the per-cell read summary, alongside the // raw matrix because its UMI and gene columns are the raw totals. - if let Some(cell_stats) = &ctx.cell_read_stats { + if let Some(cell_stats) = ctx.cell_read_stats.as_ref().map(|s| &s[fi]) { let umi_gene: std::collections::BTreeMap = mstats .cells .iter() @@ -1875,7 +1875,7 @@ pub fn write_gene_matrix( sorted.len(), gzip, )?; - let umi_len = params.solo_umi_len as usize; + let umi_len = ctx.umi_len; // STAR's per-line filters (`SoloReadFeature_inputRecords.cpp:23-104`), // counted only for a junction with a row: STAR skips the others first. let has_row = |start: u64, end: u64| row.contains_key(&(start, end)); @@ -2799,7 +2799,13 @@ pub fn run_cell_filtering(params: &crate::params::Parameters) -> anyhow::Result< std::fs::write(&bc_path, bc).map_err(|e| Error::io(e, &bc_path))?; let feat_path = dir.join(format!("{prefix}features.tsv")); - std::fs::copy(&features_path, &feat_path).map_err(|e| Error::io(e, &feat_path))?; + // The filtered matrix is written plain, whatever the raw one was. + let mut feat_out = std::fs::File::create(&feat_path).map_err(|e| Error::io(e, &feat_path))?; + std::io::copy( + &mut crate::solo::whitelist::open_maybe_gzip(&features_path)?, + &mut feat_out, + ) + .map_err(|e| Error::io(e, &feat_path))?; let mtx_path = dir.join(format!("{prefix}matrix.mtx")); let nnz = finalize_matrix( diff --git a/src/solo/mod.rs b/src/solo/mod.rs index b91d1fb..9096b1a 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -203,7 +203,51 @@ fn parse_position_anchored(spec: &str) -> Result<(Anchored, Anchored), Error> { } } +/// Refuse a `CB_UMI_Complex` position that cannot be parsed, or a UMI too wide to +/// pack. [`SoloBarcodeLayout::from_params`] assumes this has passed. +pub fn check_positions(params: &Parameters) -> Result<(), Error> { + let anchored = params.solo_adapter_sequence != "-"; + for spec in params + .solo_cb_position + .iter() + .chain(std::iter::once(¶ms.solo_umi_position)) + { + if anchored { + parse_position_anchored(spec)?; + } else { + parse_position(spec)?; + } + } + match SoloBarcodeLayout::from_params(params).umi_len() { + Some(len) if len > whitelist::CB_LEN_MAX => Err(invalid_pos( + ¶ms.solo_umi_position, + "the UMI is longer than 32 bases", + )), + _ => Ok(()), + } +} + impl SoloBarcodeLayout { + /// UMI width, or `None` when two different anchors make it vary by read. + pub fn umi_len(&self) -> Option { + match self { + Self::Simple { umi_len, .. } => Some(*umi_len), + Self::Complex { umi, .. } => Some(umi.1), + Self::Anchored { umi, .. } => anchored_len(umi), + } + } + + /// Cell barcode width of each segment, `None` where it varies by read. + pub fn cb_segment_lens(&self) -> Vec> { + match self { + Self::Simple { cb_len, .. } => vec![Some(*cb_len)], + Self::Complex { cb_segments, .. } => { + cb_segments.iter().map(|&(_, l)| Some(l)).collect() + } + Self::Anchored { cb_segments, .. } => cb_segments.iter().map(anchored_len).collect(), + } + } + /// Build the layout from CLI parameters. `CB_UMI_Complex` parses /// `--soloCBposition`/`--soloUMIposition`; otherwise fixed Simple geometry. pub fn from_params(params: &Parameters) -> Self { @@ -393,6 +437,13 @@ pub fn barcode_q30( } } +fn anchored_len(field: &(Anchored, Anchored)) -> Option { + let (start, end) = field; + (start.anchor == end.anchor) + .then(|| usize::try_from(end.dist - start.dist + 1).ok()) + .flatten() +} + fn slice_or_empty(data: &[u8], start: usize, len: usize) -> Vec { if start + len <= data.len() { data[start..start + len].to_vec() @@ -1291,7 +1342,9 @@ pub struct SoloContext { /// `--soloCellReadStats CB`: the per-cell read summary, or `None` when the /// flag is off. Behind a mutex like the other per-read collections; the /// lock is only taken when the flag asks for the file. - pub cell_read_stats: Option>, + pub cell_read_stats: Option>>, + /// UMI width: `--soloUMIposition`'s for `CB_UMI_Complex`, else `--soloUMIlen`. + pub umi_len: usize, /// Chromosome indices named by `--genomeChrSetMitochondrial`, for the /// `mito` column. pub mito_chr: std::collections::HashSet, @@ -1528,6 +1581,8 @@ impl SoloContext { genome: &crate::genome::Genome, transcriptome: Option<&crate::quant::transcriptome::TranscriptomeIndex>, ) -> Result { + let layout = SoloBarcodeLayout::from_params(params); + let cb_lens = layout.cb_segment_lens(); let whitelist = if params.solo_type == SoloType::CbUmiComplex { // One whitelist per CB segment → combined cartesian-product whitelist. let paths: Vec = params @@ -1539,7 +1594,7 @@ impl SoloContext { "STARsolo CB_UMI_Complex: combining {} segment whitelists", paths.len() ); - let wl = CbWhitelist::load_complex(&paths)?; + let wl = CbWhitelist::load_complex(&paths, &cb_lens)?; log::info!("STARsolo: {} combined whitelist barcodes", wl.len()); wl } else { @@ -1550,6 +1605,13 @@ impl SoloContext { path.display() ); let wl = CbWhitelist::load(&path)?; + if cb_lens[0].is_some_and(|want| want != wl.barcode_len()) { + return Err(whitelist::wrong_length( + &path, + wl.barcode_len(), + cb_lens[0].unwrap_or(0), + )); + } log::info!("STARsolo: {} whitelist barcodes loaded", wl.len()); wl } @@ -1621,8 +1683,10 @@ impl SoloContext { (tags || velocyto).then(|| Mutex::new(Vec::new())) }) .collect(); + let n_features = features.len(); Ok(Self { - layout: SoloBarcodeLayout::from_params(params), + umi_len: layout.umi_len().unwrap_or(params.solo_umi_len as usize), + layout, whitelist, match_type: params.solo_cb_match_type(), strand, @@ -1641,8 +1705,11 @@ impl SoloContext { velocyto_enabled, velocyto_records: Mutex::new(VelocytoRecords::default()), want_multi, - cell_read_stats: (params.solo_cell_read_stats == "CB") - .then(|| Mutex::new(crate::solo::cell_reads::CellReadStats::new())), + cell_read_stats: (params.solo_cell_read_stats == "CB").then(|| { + (0..n_features) + .map(|_| Mutex::new(crate::solo::cell_reads::CellReadStats::new())) + .collect() + }), mito_chr: params .genome_chr_set_mitochondrial .iter() @@ -2183,34 +2250,37 @@ impl SoloContext { let Some(stats) = &self.cell_read_stats else { return; }; - let counted_u = out.per_feature.iter().any(|f| f.record.is_some()); - let counted_m = out.per_feature.iter().any(|f| f.multi_gene.is_some()); - let feature_u = counted_u || out.per_feature.iter().any(|f| f.multi.is_some()); - let flag = crate::solo::cell_reads::CellReadFlag { - cb_perfect: matches!(cb_match, CbMatch::Exact(_)), - cb_mm_unique: matches!(cb_match, CbMatch::Corrected(_)), - cb_mm_multiple: matches!(cb_match, CbMatch::Multi(_)), - genome_u: n_loci == 1, - genome_m: n_loci > 1, - feature_u, - feature_m: counted_m, - // The region columns split by strand: STAR reports an antisense - // read under `exonicAS`/`intronicAS`, not under `exonic`/`intronic`. - exonic: !class.antisense && class.region == Some(Region::Exonic), - intronic: !class.antisense && class.region == Some(Region::Intronic), - exonic_as: class.antisense && class.region == Some(Region::Exonic), - intronic_as: class.antisense && class.region == Some(Region::Intronic), - mito: !self.mito_chr.is_empty() - && transcripts - .iter() - .any(|t| self.mito_chr.contains(&t.chr_idx)), - counted_u, - counted_m, - }; - let mut stats = stats.lock().unwrap(); - match cb_resolved { - Some(cb) => stats.add_cell(cb, &flag), - None => stats.add_no_cb(&flag), + // STAR keeps one file per feature, so what a read counted to stays there. + for (stats, f) in stats.iter().zip(&out.per_feature) { + let counted_u = f.record.is_some(); + let counted_m = f.multi_gene.is_some(); + let feature_u = counted_u || f.multi.is_some(); + let flag = crate::solo::cell_reads::CellReadFlag { + cb_perfect: matches!(cb_match, CbMatch::Exact(_)), + cb_mm_unique: matches!(cb_match, CbMatch::Corrected(_)), + cb_mm_multiple: matches!(cb_match, CbMatch::Multi(_)), + genome_u: n_loci == 1, + genome_m: n_loci > 1, + feature_u, + feature_m: counted_m, + // The region columns split by strand: STAR reports an antisense + // read under `exonicAS`/`intronicAS`, not under `exonic`/`intronic`. + exonic: !class.antisense && class.region == Some(Region::Exonic), + intronic: !class.antisense && class.region == Some(Region::Intronic), + exonic_as: class.antisense && class.region == Some(Region::Exonic), + intronic_as: class.antisense && class.region == Some(Region::Intronic), + mito: !self.mito_chr.is_empty() + && transcripts + .iter() + .any(|t| self.mito_chr.contains(&t.chr_idx)), + counted_u, + counted_m, + }; + let mut stats = stats.lock().unwrap(); + match cb_resolved { + Some(cb) => stats.add_cell(cb, &flag), + None => stats.add_no_cb(&flag), + } } } diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index a5287ed..142020a 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -356,11 +356,19 @@ impl CbWhitelist { /// Matching the assembled CB against this is equivalent to STARsolo's /// per-segment matching for both Exact and 1MM (a 1MM in the concatenation is /// a 1MM in exactly one segment). Errors if the combined length exceeds 32. - pub fn load_complex(paths: &[std::path::PathBuf]) -> Result { + pub fn load_complex( + paths: &[std::path::PathBuf], + segment_lens: &[Option], + ) -> Result { let segs: Vec<(Vec, usize)> = paths .iter() .map(|p| Self::load_packed(p)) .collect::>()?; + for ((path, (_, len)), want) in paths.iter().zip(&segs).zip(segment_lens) { + if want.is_some_and(|w| w != *len) { + return Err(wrong_length(path, *len, want.unwrap_or(0))); + } + } let total_len: usize = segs.iter().map(|(_, l)| l).sum(); if total_len == 0 || total_len > CB_LEN_MAX { return Err(Error::from(std::io::Error::new( @@ -497,6 +505,11 @@ impl CbWhitelist { /// `cb_seq`). On an exact hit the whitelist's exact-count is incremented. pub fn match_cb(&self, cb_seq: &[u8], cb_qual: &[u8], match_type: CbMatchType) -> CbMatch { let len = cb_seq.len(); + // Packing drops the length, so a shorter barcode would equal a longer + // entry that starts with `A`s. + if matches!(self, Self::List { len: want, .. } if *want != len) { + return CbMatch::NoMatch; + } match self { Self::NoWhitelist { .. } => match pack_barcode(cb_seq) { // No whitelist: every N-free barcode is its own "cell". We @@ -580,7 +593,17 @@ fn qual_at(qual: &[u8], pos: usize) -> u8 { } /// Open a file, transparently decompressing `.gz`. -fn open_maybe_gzip(path: &Path) -> Result, Error> { +pub(crate) fn wrong_length(path: &Path, found: usize, want: usize) -> Error { + Error::from(std::io::Error::new( + std::io::ErrorKind::InvalidData, + format!( + "{}: whitelist barcodes have length {found}, but the barcode read gives {want}", + path.display() + ), + )) +} + +pub(crate) fn open_maybe_gzip(path: &Path) -> Result, Error> { let file = File::open(path).map_err(|e| Error::io(e, path))?; let is_gz = path .extension() @@ -711,6 +734,17 @@ mod tests { f } + #[test] + fn barcode_of_another_length_never_matches() { + // Packed, `AACCC` and the 3-base entry `CCC` are the same number. + let PackResult::NoN(packed) = pack_barcode(&[1, 1, 1]) else { + unreachable!() + }; + let wl = CbWhitelist::from_packed_list(vec![packed], 3); + let matched = wl.match_cb(&[0, 0, 1, 1, 1], &[b'I'; 5], "Exact".parse().unwrap()); + assert!(matches!(matched, CbMatch::NoMatch)); + } + #[test] fn pack_roundtrip() { let s = "ACGTACGT"; diff --git a/tests/solo_input_checks.rs b/tests/solo_input_checks.rs new file mode 100644 index 0000000..ebf0f7d --- /dev/null +++ b/tests/solo_input_checks.rs @@ -0,0 +1,358 @@ +//! Solo options that used to be accepted and then give an empty or wrong result +//! with exit 0, and the per-feature and gene-name outputs that depend on them. + +use assert_cmd::cargo::cargo_bin_cmd; +use std::fs; +use std::io::{Read, Write}; +use std::path::{Path, PathBuf}; +use tempfile::TempDir; + +const CB: &str = "AAAACCCCGGGGTTTT"; +const READ_LEN: usize = 60; +/// One gene with an intron: exons at 2001-2300 and 2901-3200. +const EXONIC: usize = 2_020; +const INTRONIC: usize = 2_500; + +fn genome() -> Vec { + let bases = *b"ACGT"; + let mut state = 12_345u32; + (0..12_000) + .map(|_| { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12_345); + bases[((state >> 16) & 3) as usize] + }) + .collect() +} + +struct Fixture { + dir: TempDir, + genome_dir: PathBuf, +} + +impl Fixture { + fn new() -> Self { + let dir = TempDir::new().unwrap(); + let fasta = dir.path().join("g.fa"); + let mut f = fs::File::create(&fasta).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(&genome()).unwrap(); + writeln!(f).unwrap(); + let genome_dir = dir.path().join("idx"); + fs::create_dir_all(&genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "genomeGenerate", "--genomeSAindexNbases", "5"]) + .arg("--genomeDir") + .arg(&genome_dir) + .arg("--genomeFastaFiles") + .arg(&fasta) + .arg("--outFileNamePrefix") + .arg(dir.path().join("gen_")) + .assert() + .success(); + let gtf = dir.path().join("a.gtf"); + let attrs = r#"gene_id "G1"; transcript_id "T1"; gene_name "NAME"; gene_symbol "SYMBOL";"#; + fs::write( + >f, + format!( + "chr1\tt\texon\t2001\t2300\t.\t+\t.\t{attrs}\nchr1\tt\texon\t2901\t3200\t.\t+\t.\t{attrs}\n" + ), + ) + .unwrap(); + Self { dir, genome_dir } + } + + /// Write `barcodes` paired with reads at `at` (0-based), and return the paths. + fn reads(&self, at: &[usize], barcodes: &[String]) -> (PathBuf, PathBuf) { + let g = genome(); + let (cdna, bc) = (self.dir.path().join("c.fq"), self.dir.path().join("b.fq")); + let (mut fc, mut fb) = ( + fs::File::create(&cdna).unwrap(), + fs::File::create(&bc).unwrap(), + ); + for (i, (&p, b)) in at.iter().zip(barcodes).enumerate() { + let seq = std::str::from_utf8(&g[p..p + READ_LEN]).unwrap(); + writeln!(fc, "@r{i}\n{seq}\n+\n{}", "I".repeat(READ_LEN)).unwrap(); + writeln!(fb, "@r{i}\n{b}\n+\n{}", "I".repeat(b.len())).unwrap(); + } + (cdna, bc) + } + + fn whitelist(&self, name: &str, barcode: &str) -> PathBuf { + let p = self.dir.path().join(name); + fs::write(&p, format!("{barcode}\n")).unwrap(); + p + } + + fn run(&self, cdna: &Path, bc: &Path, args: &[&str]) -> assert_cmd::assert::Assert { + cargo_bin_cmd!("rustar-aligner") + .arg("--genomeDir") + .arg(&self.genome_dir) + .arg("--readFilesIn") + .arg(cdna) + .arg(bc) + .arg("--sjdbGTFfile") + .arg(self.dir.path().join("a.gtf")) + .arg("--outFileNamePrefix") + .arg(self.out()) + .args(args) + .assert() + } + + fn out(&self) -> PathBuf { + self.dir.path().join("out/") + } +} + +fn stderr(a: &assert_cmd::assert::Assert) -> String { + String::from_utf8_lossy(&a.get_output().stderr).into_owned() +} + +/// Refused while parsing, so nothing needs to exist on disk. +fn refused(args: &[&str]) -> String { + let a = cargo_bin_cmd!("rustar-aligner") + .args(["--genomeDir", "none", "--readFilesIn", "a", "b"]) + .args(args) + .assert() + .failure(); + stderr(&a) +} + +#[test] +fn solo_refuses_side_outputs_its_read_loops_do_not_write() { + for (args, named) in [ + (&["--quantMode", "GeneCounts"][..], "GeneCounts"), + (&["--quantMode", "TranscriptomeSAM"][..], "TranscriptomeSAM"), + (&["--outReadsUnmapped", "Fastx"][..], "outReadsUnmapped"), + (&["--outFilterType", "BySJout"][..], "BySJout"), + (&["--chimSegmentMin", "12"][..], "chimSegmentMin"), + ] { + let mut all = vec!["--soloType", "CB_UMI_Simple"]; + all.extend_from_slice(args); + let err = refused(&all); + assert!(err.contains(named), "{args:?}: {err}"); + } +} + +#[test] +fn umi_longer_than_a_u64_is_refused() { + let err = refused(&["--soloType", "CB_UMI_Simple", "--soloUMIlen", "33"]); + assert!(err.contains("--soloUMIlen 33"), "{err}"); +} + +#[test] +fn several_umi_dedup_methods_are_refused() { + let err = refused(&[ + "--soloType", + "CB_UMI_Simple", + "--soloUMIdedup", + "1MM_All", + "Exact", + ]); + assert!(err.contains("--soloUMIdedup"), "{err}"); +} + +#[test] +fn complex_needs_a_valid_umi_position() { + let base = [ + "--soloType", + "CB_UMI_Complex", + "--soloCBwhitelist", + "wl", + "--soloCBposition", + "0_0_0_15", + ]; + assert!(refused(&base).contains("--soloUMIposition")); + let mut bad = base.to_vec(); + bad.extend_from_slice(&["--soloUMIposition", "0_16_nonsense"]); + assert!(refused(&bad).contains("0_16_nonsense")); + let mut wide = base.to_vec(); + wide.extend_from_slice(&["--soloUMIposition", "0_16_0_48"]); + assert!(refused(&wide).contains("longer than 32")); + let mut bad_cb = base[..4].to_vec(); + bad_cb.extend_from_slice(&["--soloCBposition", "oops", "--soloUMIposition", "0_16_0_25"]); + assert!(refused(&bad_cb).contains("oops")); +} + +#[test] +fn whitelist_of_another_length_is_fatal() { + let fx = Fixture::new(); + let (c, b) = fx.reads(&[EXONIC], &[format!("{CB}AAAAACCCCC")]); + let wl = fx.whitelist("wl14.txt", &CB[..14]); + let a = fx + .run( + &c, + &b, + &[ + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl.to_str().unwrap(), + ], + ) + .failure(); + assert!(stderr(&a).contains("length 14"), "{}", stderr(&a)); +} + +#[test] +fn complex_segment_whitelist_of_another_length_is_fatal() { + let fx = Fixture::new(); + let (c, b) = fx.reads(&[EXONIC], &[format!("{CB}AAAAACCCCC")]); + let (w1, w2) = ( + fx.whitelist("a.txt", "AAAAAA"), + fx.whitelist("b.txt", "CCCCCCCCCC"), + ); + // 6 + 10 is the right total; each file is the wrong way round. + let a = fx + .run( + &c, + &b, + &[ + "--soloType", + "CB_UMI_Complex", + "--soloCBwhitelist", + w1.to_str().unwrap(), + w2.to_str().unwrap(), + "--soloCBposition", + "0_0_0_9", + "0_10_0_15", + "--soloUMIposition", + "0_16_0_25", + ], + ) + .failure(); + assert!(stderr(&a).contains("length 6"), "{}", stderr(&a)); +} + +#[test] +fn complex_umi_width_comes_from_its_position() { + let fx = Fixture::new(); + let (c, b) = fx.reads(&[EXONIC], &[format!("{CB}ACGTACGT")]); + let wl = fx.whitelist("wl.txt", CB); + fx.run( + &c, + &b, + &[ + "--soloType", + "CB_UMI_Complex", + "--soloCBwhitelist", + wl.to_str().unwrap(), + "--soloCBposition", + "0_0_0_15", + "--soloUMIposition", + "0_16_0_23", + "--soloBarcodeReadLength", + "0", + "--outSAMtype", + "BAM", + "SortedByCoordinate", + "--outSAMattributes", + "NH", + "CB", + "UB", + ], + ) + .success(); + let bam = fs::File::open(fx.out().join("Aligned.sortedByCoord.out.bam")).unwrap(); + let mut decoded = Vec::new(); + flate2::read::MultiGzDecoder::new(bam) + .read_to_end(&mut decoded) + .unwrap(); + let tag = b"UBZACGTACGT\0"; + assert!( + decoded.windows(tag.len()).any(|w| w == tag), + "the UB tag must be the eight bases of --soloUMIposition" + ); +} + +/// An intronic read counts to GeneFull and not to Gene, so the two features' +/// files must disagree about it. +#[test] +fn cell_reads_stats_are_per_feature() { + let fx = Fixture::new(); + let (c, b) = fx.reads( + &[EXONIC, INTRONIC], + &[format!("{CB}AAAAACCCCC"), format!("{CB}GGGGGTTTTT")], + ); + let wl = fx.whitelist("wl.txt", CB); + fx.run( + &c, + &b, + &[ + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl.to_str().unwrap(), + "--soloFeatures", + "Gene", + "GeneFull", + "--soloCellReadStats", + "CB", + ], + ) + .success(); + let counted_u = |feature: &str| -> String { + let text = fs::read_to_string(fx.out().join(format!("Solo.out/{feature}/CellReads.stats"))) + .unwrap(); + let mut lines = text.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + let col = header.iter().position(|h| *h == "countedU").unwrap(); + let row = lines.find(|l| l.starts_with(CB)).unwrap(); + row.split('\t').nth(col).unwrap().to_string() + }; + assert_eq!(counted_u("Gene"), "1"); + assert_eq!(counted_u("GeneFull"), "2"); +} + +#[test] +fn gene_name_tag_applies_to_a_gtf_given_at_mapping() { + let fx = Fixture::new(); + let (c, b) = fx.reads(&[EXONIC], &[format!("{CB}AAAAACCCCC")]); + let wl = fx.whitelist("wl.txt", CB); + fx.run( + &c, + &b, + &[ + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl.to_str().unwrap(), + "--sjdbGTFtagExonParentGeneName", + "gene_symbol", + ], + ) + .success(); + let features = fs::read_to_string(fx.out().join("Solo.out/Gene/raw/features.tsv")).unwrap(); + assert!(features.contains("G1\tSYMBOL"), "{features}"); +} + +/// `soloCellFiltering` writes a plain matrix, so a gzipped feature list has to be +/// decompressed on the way through, not copied into a `.tsv`. +#[test] +fn cell_filtering_writes_plain_features_from_gzipped_ones() { + let dir = TempDir::new().unwrap(); + let raw = dir.path().join("raw"); + fs::create_dir_all(&raw).unwrap(); + fs::write( + raw.join("matrix.mtx"), + "%%MatrixMarket matrix coordinate integer general\n%\n2 2 2\n1 1 500\n2 2 1\n", + ) + .unwrap(); + fs::write(raw.join("barcodes.tsv"), "AAAA\nCCCC\n").unwrap(); + let mut gz = flate2::write::GzEncoder::new( + fs::File::create(raw.join("features.tsv.gz")).unwrap(), + flate2::Compression::default(), + ); + gz.write_all(b"g1\tG1\tGene Expression\ng2\tG2\tGene Expression\n") + .unwrap(); + gz.finish().unwrap(); + let out = dir.path().join("filtered/"); + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "soloCellFiltering"]) + .arg(&raw) + .arg(&out) + .args(["--soloCellFilter", "TopCells", "1"]) + .assert() + .success(); + let written = fs::read(out.join("features.tsv")).unwrap(); + assert!(written.starts_with(b"g1\tG1"), "{written:?}"); +}