From 01a1a1f3261bcaa14c6c7c8d9689324b1a16eb3a Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 29 Jul 2026 10:17:14 +0200 Subject: [PATCH] feat(solo): --runMode soloCellFiltering Cell-calls an existing raw count matrix without aligning anything, taking the raw directory and an output prefix as STAR does: `--runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix`. Cell calling is a decision about a matrix, not about reads. Re-calling with different `--soloCellFilter` parameters should not mean re-aligning 400 million reads, and a matrix produced by another tool should be callable too. The matrix is streamed into the same temp-body form the align path builds, so `called_cells` and `emptydrops_called` are the identical code here and there rather than a second implementation free to drift. Counts are rounded on the way in: a multimapper matrix carries real values, and the filters work on UMI totals. `--runMode` becomes a token list, because that is what STAR's is: the mode followed by its arguments. The mode itself is now validated rather than falling back to `alignReads`, so a typo is refused instead of quietly running something else. The standalone `emptydrops` binary still exists and still carries its own copy of the algorithm, which no longer matches this one. Removing it means moving `test/solo_genefull_compare.py` and `test/solo_genefull_h5_compare.py` to the new mode first, so it is left alone here rather than broken. --- CHANGELOG.md | 8 ++ src/lib.rs | 5 +- src/params/mod.rs | 63 +++++++++--- src/solo/count.rs | 191 ++++++++++++++++++++++++++++++++++++ tests/alignment_features.rs | 90 +++++++++++++++++ 5 files changed, 342 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b6e4893..31b6b73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -97,6 +97,14 @@ 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. +- **`--runMode soloCellFiltering `** cell-calls + an existing raw count matrix without aligning anything. Cell calling + is a decision about a matrix, not about reads: re-calling with + different `--soloCellFilter` parameters should not mean re-aligning, + and a matrix produced elsewhere should be callable too. It streams + the matrix into the same form the align path produces, so the filters + are the identical code rather than a second implementation. + ### Bug fixes - **STARsolo `Gene` assignment now requires exon concordance**, matching diff --git a/src/lib.rs b/src/lib.rs index 6fce173..d38274d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -55,7 +55,7 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { if let Some(hint) = cpu::upgrade_hint() { info!("{hint}"); } - info!("runMode: {}", params.run_mode); + info!("runMode: {}", params.run_mode_in.join(" ")); info!("runThreadN: {}", params.run_thread_n); // Configure the rayon global pool from `--runThreadN` **before** @@ -74,11 +74,12 @@ pub fn run(params: &Parameters) -> anyhow::Result<()> { .build_global(); } - match params.run_mode { + match params.run_mode() { RunMode::GenomeGenerate => genome_generate(params), RunMode::AlignReads => align_reads(params), RunMode::InputAlignmentsFromBAM => bam_dedup::run(params), RunMode::LiftOver => liftover::run(params), + RunMode::SoloCellFiltering => crate::solo::count::run_cell_filtering(params), } } diff --git a/src/params/mod.rs b/src/params/mod.rs index 0536a85..fc5c0e5 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -37,6 +37,8 @@ pub enum RunMode { GenomeGenerate, InputAlignmentsFromBAM, LiftOver, + /// Cell-call an existing raw count matrix, without aligning anything. + SoloCellFiltering, } impl std::str::FromStr for RunMode { @@ -47,9 +49,10 @@ impl std::str::FromStr for RunMode { "genomeGenerate" => Ok(Self::GenomeGenerate), "inputAlignmentsFromBAM" => Ok(Self::InputAlignmentsFromBAM), "liftOver" => Ok(Self::LiftOver), + "soloCellFiltering" => Ok(Self::SoloCellFiltering), _ => Err(format!( "unknown runMode '{s}'; expected 'alignReads', 'genomeGenerate', \ - 'inputAlignmentsFromBAM', or 'liftOver'" + 'inputAlignmentsFromBAM', 'liftOver', or 'soloCellFiltering'" )), } } @@ -62,6 +65,7 @@ impl std::fmt::Display for RunMode { Self::GenomeGenerate => write!(f, "genomeGenerate"), Self::InputAlignmentsFromBAM => write!(f, "inputAlignmentsFromBAM"), Self::LiftOver => write!(f, "liftOver"), + Self::SoloCellFiltering => write!(f, "soloCellFiltering"), } } } @@ -465,9 +469,11 @@ impl std::fmt::Display for SoloType { )] pub struct Parameters { // ── Run ───────────────────────────────────────────────────────────── - /// Run mode: alignReads or genomeGenerate - #[arg(long = "runMode", default_value = "alignReads")] - pub run_mode: RunMode, + /// Run mode, plus its arguments. `--runMode soloCellFiltering` takes two + /// more tokens: the raw count-matrix directory and the output prefix + /// (STAR `SoloFeature_loadRawMatrix.cpp`). + #[arg(long = "runMode", num_args = 1.., default_values_t = vec!["alignReads".to_string()])] + pub run_mode_in: Vec, /// Number of threads #[arg(long = "runThreadN", default_value_t = NonZeroUsize::new(1).unwrap())] @@ -1162,6 +1168,16 @@ pub struct Parameters { } impl Parameters { + /// The run mode. Validation guarantees it parses, so this cannot fail + /// after `validate()`; before it, an unknown mode reads as `alignReads` + /// and validation is what rejects it. + pub fn run_mode(&self) -> RunMode { + self.run_mode_in + .first() + .and_then(|m| m.parse().ok()) + .unwrap_or(RunMode::AlignReads) + } + /// Build an output path by concatenating `suffix` onto `out_file_name_prefix`. pub fn output_path(&self, suffix: &str) -> PathBuf { PathBuf::from(format!("{}{suffix}", self.out_file_name_prefix)) @@ -1353,8 +1369,29 @@ impl Parameters { shlex::try_join(args.iter().map(AsRef::as_ref)).ok() }; + // The run mode itself must be one this build knows: an unrecognised + // one would otherwise fall through to alignReads and silently do + // something the user did not ask for. + if let Some(mode) = params.run_mode_in.first() + && mode.parse::().is_err() + { + return Err(command.error( + ErrorKind::InvalidValue, + mode.parse::().unwrap_err(), + )); + } + + // `--runMode soloCellFiltering `. + if params.run_mode() == RunMode::SoloCellFiltering && params.run_mode_in.len() < 3 { + return Err(command.error( + ErrorKind::WrongNumberOfValues, + "--runMode soloCellFiltering needs the raw count-matrix directory and the \ + output prefix: --runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix", + )); + } + // genomeGenerate requires FASTA files - if params.run_mode == RunMode::GenomeGenerate && params.genome_fasta_files.is_empty() { + if params.run_mode() == RunMode::GenomeGenerate && params.genome_fasta_files.is_empty() { return Err(command.error( ErrorKind::MissingRequiredArgument, "--genomeFastaFiles is required when --runMode genomeGenerate", @@ -1371,7 +1408,7 @@ impl Parameters { // alignReads requires read files — except SmartSeq, which gets its reads // from --readFilesManifest instead. - if params.run_mode == RunMode::AlignReads + if params.run_mode() == RunMode::AlignReads && params.read_files_in.is_empty() && params.solo_type != SoloType::SmartSeq { @@ -1408,7 +1445,7 @@ impl Parameters { } // inputAlignmentsFromBAM: only --bamRemoveDuplicatesType is implemented so far - if params.run_mode == RunMode::InputAlignmentsFromBAM { + if params.run_mode() == RunMode::InputAlignmentsFromBAM { let dedup = params.bam_remove_duplicates_type.as_str(); if dedup == "-" { return Err(command.error( @@ -1431,7 +1468,7 @@ impl Parameters { } // liftOver requires a chain file and a GTF to lift - if params.run_mode == RunMode::LiftOver { + if params.run_mode() == RunMode::LiftOver { if params.genome_chain_files.is_empty() { return Err(command.error( ErrorKind::MissingRequiredArgument, @@ -1541,7 +1578,7 @@ impl Parameters { // validation time we can only enforce the genomeGenerate rule; // for alignReads, GenomeIndex::load checks for the on-disk files // and surfaces a clear error if neither source is available. - if params.run_mode == RunMode::GenomeGenerate + if params.run_mode() == RunMode::GenomeGenerate && params.quant_transcriptome_sam() && params.sjdb_gtf_file.is_none() { @@ -1552,7 +1589,7 @@ impl Parameters { } // ── STARsolo validation ───────────────────────────────────────── - if params.run_mode == RunMode::AlignReads && params.solo_enabled() { + if params.run_mode() == RunMode::AlignReads && params.solo_enabled() { // CB_UMI_Complex needs one CB position + whitelist per segment. if params.solo_type == SoloType::CbUmiComplex { if params.solo_cb_position.is_empty() { @@ -1884,7 +1921,7 @@ mod tests { #[test] fn defaults() { let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap(); - assert_eq!(p.run_mode, RunMode::AlignReads); + assert_eq!(p.run_mode(), RunMode::AlignReads); assert_eq!(p.run_thread_n, NonZeroUsize::new(1).unwrap()); assert_eq!(p.run_rng_seed, 777); assert_eq!(p.genome_dir, PathBuf::from("./GenomeDir")); @@ -1978,7 +2015,7 @@ mod tests { "11", ]) .unwrap(); - assert_eq!(p.run_mode, RunMode::GenomeGenerate); + assert_eq!(p.run_mode(), RunMode::GenomeGenerate); assert_eq!(p.genome_dir, PathBuf::from("/data/genome")); assert_eq!( p.genome_fasta_files, @@ -2018,7 +2055,7 @@ mod tests { "Basic", ]) .unwrap(); - assert_eq!(p.run_mode, RunMode::AlignReads); + assert_eq!(p.run_mode(), RunMode::AlignReads); assert_eq!(p.genome_dir, PathBuf::from("/idx/hg38")); assert_eq!( p.read_files_in, diff --git a/src/solo/count.rs b/src/solo/count.rs index 9af3edf..576462c 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -1881,6 +1881,197 @@ fn write_barcodes_subset( Ok(()) } +/// `--runMode soloCellFiltering `: cell-call an +/// existing raw count matrix without aligning anything. +/// +/// STAR `SoloFeature_loadRawMatrix.cpp` plus the same `--soloCellFilter` used +/// at the end of a solo run. The point of the mode is that cell calling is a +/// decision about a matrix, not about reads: re-calling with different +/// parameters should not mean re-aligning, and a matrix produced elsewhere +/// should be callable too. +pub fn run_cell_filtering(params: &crate::params::Parameters) -> anyhow::Result<()> { + let raw_dir = Path::new(¶ms.run_mode_in[1]); + let out_prefix = ¶ms.run_mode_in[2]; + + let find = |base: &str| -> Result { + for name in [base.to_string(), format!("{base}.gz")] { + let p = raw_dir.join(&name); + if p.exists() { + return Ok(p); + } + } + Err(Error::Parameter(format!( + "{}: no {base} (or {base}.gz) in the raw matrix directory", + raw_dir.display() + ))) + }; + + let barcodes = read_first_column(&find("barcodes.tsv")?)?; + let features_path = find("features.tsv")?; + let matrix_path = find("matrix.mtx")?; + + // Stream the matrix into the same temp-body form the align path produces, + // so the filters below are the identical code rather than a second + // implementation that can drift from it. + let mut body = tempfile::NamedTempFile::new().map_err(|e| Error::io(e, raw_dir))?; + let mut totals: HashMap = HashMap::default(); + let mut n_features = 0usize; + { + let mut out = std::io::BufWriter::new(body.as_file_mut()); + let reader = open_maybe_gz(&matrix_path)?; + let mut header_seen = false; + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, &matrix_path))?; + if line.starts_with('%') { + continue; + } + let mut fields = line.split_whitespace(); + let (Some(first), Some(second), Some(third)) = + (fields.next(), fields.next(), fields.next()) + else { + continue; + }; + if !header_seen { + // ` ` + n_features = first.parse().unwrap_or(0); + header_seen = true; + continue; + } + let (gene, cb, count) = ( + first.parse::().unwrap_or(0), + second.parse::().unwrap_or(0), + third.parse::().unwrap_or(0.0), + ); + if gene == 0 || cb == 0 { + continue; + } + // Counts can be real-valued in a multimapper matrix; the cell + // filters work on UMI totals, so round rather than refuse. + let count = count.round().max(0.0) as u64; + let e = totals.entry(cb - 1).or_insert((0, 0)); + e.0 += count; + e.1 += 1; + writeln!(out, "{gene} {cb} {count}").map_err(|e| Error::io(e, raw_dir))?; + } + out.flush().map_err(|e| Error::io(e, raw_dir))?; + } + if n_features == 0 { + return Err(Error::Parameter(format!( + "{}: no MatrixMarket header, so the matrix shape is unknown", + matrix_path.display() + )) + .into()); + } + + let mut cells: Vec = totals + .iter() + .map(|(&cb, &(n_umis, n_genes))| CellStat { + cb, + n_reads: n_umis, + n_umis, + n_genes, + }) + .collect(); + cells.sort_unstable_by_key(|c| c.cb); + log::info!( + "soloCellFiltering: {} barcodes with counts, {n_features} features", + cells.len() + ); + + let called = if params + .solo_cell_filter + .first() + .is_some_and(|m| m == "EmptyDrops_CR") + { + Some(emptydrops_called( + &cells, + &body, + n_features, + ¶ms.solo_cell_filter, + )?) + } else { + called_cells(&cells, ¶ms.solo_cell_filter) + }; + let Some(cbs) = called.filter(|c| !c.is_empty()) else { + log::warn!("soloCellFiltering: no cells called; writing nothing"); + return Ok(()); + }; + + let out_dir = Path::new(out_prefix); + let (dir, prefix): (&Path, &str) = if out_prefix.ends_with('/') { + (out_dir, "") + } else { + ( + out_dir.parent().unwrap_or_else(|| Path::new(".")), + out_dir.file_name().and_then(|s| s.to_str()).unwrap_or(""), + ) + }; + std::fs::create_dir_all(dir).map_err(|e| Error::io(e, dir))?; + + let remap: HashMap = cbs + .iter() + .enumerate() + .map(|(i, &cb)| (cb, i as u32 + 1)) + .collect(); + + let bc_path = dir.join(format!("{prefix}barcodes.tsv")); + let mut bc = String::new(); + for &cb in &cbs { + let Some(name) = barcodes.get(cb as usize) else { + continue; + }; + bc.push_str(name); + bc.push('\n'); + } + 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))?; + + let mtx_path = dir.join(format!("{prefix}matrix.mtx")); + let nnz = finalize_matrix( + &body, + &mtx_path, + false, + n_features, + cbs.len(), + 0, + Some(&remap), + )?; + log::info!( + "soloCellFiltering: {} cells, {nnz} entries -> {}", + cbs.len(), + dir.display() + ); + Ok(()) +} + +/// First whitespace-separated column of a (possibly gzipped) file. +fn read_first_column(path: &Path) -> Result, Error> { + let reader = open_maybe_gz(path)?; + let mut out = Vec::new(); + for line in reader.lines() { + let line = line.map_err(|e| Error::io(e, path))?; + if let Some(first) = line.split_whitespace().next() { + out.push(first.to_string()); + } + } + Ok(out) +} + +/// Open a file, transparently decompressing `.gz`. +fn open_maybe_gz(path: &Path) -> Result, Error> { + let file = std::fs::File::open(path).map_err(|e| Error::io(e, path))?; + if path + .extension() + .is_some_and(|e| e.eq_ignore_ascii_case("gz")) + { + Ok(Box::new(BufReader::new(flate2::read::GzDecoder::new(file)))) + } else { + Ok(Box::new(BufReader::new(file))) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 2db7d39..62ee32a 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2024,3 +2024,93 @@ fn test_wasp_samtag() { "all 10 unique reads overlapping the het SNV should pass WASP (vW:i:1)" ); } + +// --------------------------------------------------------------------------- +// --runMode soloCellFiltering +// --------------------------------------------------------------------------- + +/// Cell calling is a decision about a matrix, not about reads: the mode takes +/// an existing raw matrix and writes the called subset, with no genome and no +/// FASTQ involved. +#[test] +fn test_run_mode_solo_cell_filtering_calls_cells_from_a_raw_matrix() { + let tmpdir = TempDir::new().unwrap(); + let raw = tmpdir.path().join("raw"); + fs::create_dir_all(&raw).unwrap(); + + // 20 barcodes: five real cells with 1000 UMIs each, fifteen with 2. + let n_features = 50usize; + let n_cells = 20usize; + let mut entries: Vec<(usize, usize, u64)> = Vec::new(); + for cb in 1..=n_cells { + let per_gene = if cb <= 5 { 100 } else { 1 }; + let n_genes = if cb <= 5 { 10 } else { 2 }; + for gene in 1..=n_genes { + entries.push((gene, cb, per_gene)); + } + } + { + let mut f = fs::File::create(raw.join("matrix.mtx")).unwrap(); + writeln!(f, "%%MatrixMarket matrix coordinate integer general").unwrap(); + writeln!(f, "%").unwrap(); + writeln!(f, "{} {} {}", n_features, n_cells, entries.len()).unwrap(); + for (g, c, v) in &entries { + writeln!(f, "{g} {c} {v}").unwrap(); + } + } + { + let mut f = fs::File::create(raw.join("barcodes.tsv")).unwrap(); + for cb in 0..n_cells { + writeln!(f, "{:016b}", cb).unwrap(); + } + } + { + let mut f = fs::File::create(raw.join("features.tsv")).unwrap(); + for g in 0..n_features { + writeln!(f, "gene{g}\tGENE{g}\tGene Expression").unwrap(); + } + } + + let out = tmpdir.path().join("filtered/"); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "soloCellFiltering", + raw.to_str().unwrap(), + out.to_str().unwrap(), + "--soloCellFilter", + "TopCells", + "5", + ]) + .assert() + .success(); + + let barcodes = fs::read_to_string(out.join("barcodes.tsv")).unwrap(); + assert_eq!( + barcodes.lines().count(), + 5, + "the five deep barcodes are the called cells" + ); + + let matrix = fs::read_to_string(out.join("matrix.mtx")).unwrap(); + let header = matrix.lines().nth(2).unwrap(); + let fields: Vec<&str> = header.split_whitespace().collect(); + assert_eq!(fields[0], "50", "features are carried through"); + assert_eq!(fields[1], "5", "columns are the called cells"); + assert_eq!(fields[2], "50", "10 genes × 5 cells"); + + assert!( + out.join("features.tsv").exists(), + "the feature list travels with the matrix" + ); +} + +/// The mode needs both paths; asking for it without them is refused rather +/// than run against a guess. +#[test] +fn test_run_mode_solo_cell_filtering_requires_its_paths() { + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "soloCellFiltering"]) + .assert() + .failure(); +}