solo: --runMode soloCellFiltering

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 <raw dir> <output prefix>`** 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

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),
     }
 }
 

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<String>,
 
     /// 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::<RunMode>().is_err()
+        {
+            return Err(command.error(
+                ErrorKind::InvalidValue,
+                mode.parse::<RunMode>().unwrap_err(),
+            ));
+        }
+
+        // `--runMode soloCellFiltering <raw dir> <output prefix>`.
+        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,

src/solo/count.rs

@@ -1881,6 +1881,197 @@ fn write_barcodes_subset(
     Ok(())
 }
 
+/// `--runMode soloCellFiltering <raw dir> <output prefix>`: 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(&params.run_mode_in[1]);
+    let out_prefix = &params.run_mode_in[2];
+
+    let find = |base: &str| -> Result<PathBuf, Error> {
+        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<u32, (u64, u32)> = 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 {
+                // `<features> <barcodes> <entries>`
+                n_features = first.parse().unwrap_or(0);
+                header_seen = true;
+                continue;
+            }
+            let (gene, cb, count) = (
+                first.parse::<u32>().unwrap_or(0),
+                second.parse::<u32>().unwrap_or(0),
+                third.parse::<f64>().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<CellStat> = 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,
+            &params.solo_cell_filter,
+        )?)
+    } else {
+        called_cells(&cells, &params.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<u32, u32> = 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<Vec<String>, 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<Box<dyn BufRead>, 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::*;