feat: output anndata

Cargo.toml

@@ -35,6 +35,7 @@ memmap2 = "0.9"
 byteorder = "1"
 noodles = { version = "0.113", features = ["fastq", "sam", "bam", "bgzf"] }
 bstr = "1"
+nalgebra-sparse = { version = "0.11", optional = true }
 # Use the pure-Rust `zlib-rs` backend instead of flate2's default `miniz_oxide`:
 # ~2-3x faster inflate/deflate on the FASTQ decode + BGZF (BAM) paths, with no C
 # toolchain dependency. Backend is chosen at build time; API is unchanged.
@@ -66,9 +67,28 @@ libmimalloc-sys = { version = "0.1.49", features = ["extended"] } # mi_option_se
 libdeflater = "1.25.2"
 noodles-bgzf = { version = "0.49", features = ["libdeflate"] }
 
+# --- AnnData/MuData output (--soloOutputFormat Zarr) -----------------------
+# scverse/anndata-rs. Optional because the tree is heavy (polars for the
+# obs/var DataFrames, zarrs for the Zarr v3 store). `--no-default-features`
+# gives back the lean aligner-only build; `--soloOutputFormat` then only
+# accepts `MTX`.
+#
+# `nalgebra-sparse` / `polars` are pinned to the versions anndata 0.7 depends
+# on — the `CsrMatrix` and `DataFrame` we hand to it must be the *same* types.
+anndata = { version = "0.7", optional = true }
+anndata-zarr = { version = "0.2", optional = true }
+polars = { version = "0.53", default-features = false, optional = true }
+
+[features]
+default = ["anndata-out"]
+# Shared plumbing for AnnData backend.
+anndata-out = ["dep:anndata", "dep:nalgebra-sparse", "dep:polars", "dep:anndata-zarr"]
+
 [dev-dependencies]
 assert_cmd = "2"
 predicates = "3"
+# Reading back the hand-written Zarr group metadata in the MuData output test.
+serde_json = "1"
 
 [build-dependencies]
 chrono = { version = "0.4", default-features = false, features = ["clock"] }
@@ -86,3 +106,10 @@ opt-level = 3
 lto = "fat"
 codegen-units = 1
 strip = true
+
+# `obsm`/`varm` DataFrames only get their index (the `obs`/`var` names) stamped
+# on by anndata-rs `ig/set_index`, which is what the solo Zarr writer relies on.
+# Drop this section once that lands in a crates.io release.
+[patch.crates-io]
+anndata = { path = "../anndata-rs/anndata" }
+anndata-zarr = { path = "../anndata-rs/anndata-zarr" }

README.md

@@ -216,7 +216,7 @@ resident; the 16 GB sparse index is stable at ~54 s.</sub>
 - Chimeric alignment detection — SE and PE, 4-tier pipeline: transcript-pair search, multi-cluster, soft-clip re-seeding, residual outer re-seeding for multi-junction fusions (`--chimSegmentMin`)
 - Gene-level read counting (`--quantMode GeneCounts` → `ReadsPerGene.out.tab`)
 - Transcriptome-coordinate SAM output (`--quantMode TranscriptomeSAM`)
-- **Single-cell quantification (STARsolo)** — `--soloType CB_UMI_Simple`, `CB_UMI_Complex` (multi-segment barcodes), and `SmartSeq` (plate-based, SE + PE); features `Gene`, `GeneFull` (pre-mRNA), `SJ`, and `Velocyto` (spliced/unspliced/ambiguous); barcode correction (`--soloCBmatchWLtype` Exact/1MM/1MM_multi/…), UMI dedup (`--soloUMIdedup` 1MM_All/1MM_CR/1MM_Directional/…), multi-gene UMI filtering, multi-mapper resolution (`--soloMultiMappers` Uniform/PropUnique/EM/Rescue), cell calling (`--soloCellFilter` CellRanger2.2/TopCells/EmptyDrops_CR), gzip output, and `Summary.csv` — writes STARsolo-compatible `Solo.out/<feature>/{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}`
+- **Single-cell quantification (STARsolo)** — `--soloType CB_UMI_Simple`, `CB_UMI_Complex` (multi-segment barcodes), and `SmartSeq` (plate-based, SE + PE); features `Gene`, `GeneFull` (pre-mRNA), `SJ`, and `Velocyto` (spliced/unspliced/ambiguous); barcode correction (`--soloCBmatchWLtype` Exact/1MM/1MM_multi/…), UMI dedup (`--soloUMIdedup` 1MM_All/1MM_CR/1MM_Directional/…), multi-gene UMI filtering, multi-mapper resolution (`--soloMultiMappers` Uniform/PropUnique/EM/Rescue), cell calling (`--soloCellFilter` CellRanger2.2/TopCells/EmptyDrops_CR), gzip output, and `Summary.csv` — writes STARsolo-compatible `Solo.out/<feature>/{raw,filtered}/{matrix.mtx, barcodes.tsv, features.tsv}`, or, with `--soloOutputFormat Zarr`, one MuData store (`Solo.out/matrix.zarr`) holding every feature as a layer of a `gex` modality plus a junction-indexed `sj` modality
 - WASP allele-specific-mapping filter (`--waspOutputMode SAMtag`, `--varVCFfile`) — vW/vA/vG tags, SE + PE
 - Paired-end mate-overlap merging (`--peOverlapNbasesMin`, `--peOverlapMMp`)
 - Coverage-track output (`--outWigType bedGraph` → `Signal.{Unique,UniqueMultiple}.str{1,2}.out.bg`)

ROADMAP.md

@@ -343,7 +343,7 @@ Single-cell quantification layered around the existing aligner: the cDNA read al
 
 **Live verification — PASS:** rustar-aligner's `Gene/raw` matrix is **byte-identical to real STARsolo's** for the CellRanger-style run, confirmed deterministically (3/3 runs). The reference STAR (2.7.10b) and a Linux build of rustar-aligner run in a consistent Linux container (`test/Dockerfile.solodiff` + `test/solo_diff_docker.sh`, via colima — no Docker Desktop). This was necessary because STAR 2.7.11b reads 0 input reads on Apple-Silicon macOS (a known STAR/macOS bug, `nextChar=-1`). 479 lib + 11 integration tests, 0 clippy warnings.
 
-**Phase 14.4 — MVP COMPLETE** (2026-06-11): UMI deduplication + raw count-matrix output. New `src/solo/count.rs`: `UmiDedup` (`--soloUMIdedup`: Exact / NoDedup / 1MM_All [default, connected-components within Hamming-1] / 1MM_Directional / 1MM_Directional_UMItools, `dirCountAdd` 0/−1); deferred 1MM_multi CB resolution via STAR's count+quality posterior (weight = `exactCount·10^(−q/10)`, prior from `whitelist.exact_count_snapshot()`); `build_matrix` groups reads by (cell,gene), collapses UMIs, and `write_gene_matrix` writes `Solo.out/Gene/raw/{matrix.mtx, barcodes.tsv, features.tsv}` (MatrixMarket `nFeatures nBarcodes nEntries`, entries `gene+1 cell+1 count`, 1-based; CellRanger-v3 3-column features.tsv; whitelist-sorted barcodes.tsv). Wired into `align_reads` post-alignment. `--soloUMIdedup` validation in params. End-to-end test (`test_starsolo_gene_matrix`): 8 reads, one cell, two Hamming-distant UMI clouds → 2 deduped molecules → matrix `1 1 2`. **A working 10x Chromium Gene count matrix.** 475 lib + 10 integration tests, 0 clippy warnings.
+**Phase 14.4 — MVP COMPLETE** (2026-06-11): UMI deduplication + raw count-matrix output. New `src/solo/count.rs`: `UmiDedup` (`--soloUMIdedup`: Exact / NoDedup / 1MM_All [default, connected-components within Hamming-1] / 1MM_Directional / 1MM_Directional_UMItools, `dirCountAdd` 0/−1); deferred 1MM_multi CB resolution via STAR's count+quality posterior (weight = `exactCount·10^(−q/10)`, prior from `whitelist.exact_count_snapshot()`); `build_matrix` groups reads by (cell,gene), collapses UMIs, and `write_matrix_market` writes `Solo.out/Gene/raw/{matrix.mtx, barcodes.tsv, features.tsv}` (MatrixMarket `nFeatures nBarcodes nEntries`, entries `gene+1 cell+1 count`, 1-based; CellRanger-v3 3-column features.tsv; whitelist-sorted barcodes.tsv). Wired into `align_reads` post-alignment. `--soloUMIdedup` validation in params. End-to-end test (`test_starsolo_gene_matrix`): 8 reads, one cell, two Hamming-distant UMI clouds → 2 deduped molecules → matrix `1 1 2`. **A working 10x Chromium Gene count matrix.** 475 lib + 10 integration tests, 0 clippy warnings.
 
 **Phase 14.5–14.11 + performance** (2026-07): completed the feature-parity set — `Summary.csv` (STARsolo-faithful, CellRanger funnel split to its own file), `--soloCellFilter` CellRanger2.2/TopCells/**EmptyDrops_CR** (Monte-Carlo ambient rescue in the `filtered/` writer), `--soloFeatures` **GeneFull/SJ/Velocyto** (spliced/unspliced/ambiguous per Sullivan 2025), `--soloMultiMappers` Uniform/PropUnique/EM/Rescue, chemistries **CB_UMI_Complex** (multi-segment) and **SmartSeq** (plate-based, SE + PE fragment counts), and a rustar-vs-STARsolo SJ + multi-mapper diff harness. Performance: pipelined solo FASTQ decode, parallelized matrix build + EmptyDrops MC, libdeflate/zlib-rs for matrix gzip + BGZF, and an **O(log n + k) segment-tree gene-overlap query** (replacing STAR's linear scan — the #1 solo hotspot, ~14% wall reduction). Sparse suffix array (`--genomeSAsparseD`, byte-identical to STAR's D=2) for a 31% smaller index. 516 tests, 0 clippy warnings.
 

docs-old/phase14_starsolo.md

@@ -238,7 +238,7 @@ usable single-cell output.
   when no candidate has positive weight.
 - **`build_matrix`** groups reads by `(cell, gene)` into UMI→multiplicity maps
   (resolved multi-CB records folded in), then dedups each.
-- **`write_gene_matrix`** writes `Solo.out/Gene/raw/`:
+- **`write_matrix_market`** writes `Solo.out/Gene/raw/`:
   - `matrix.mtx` — `%%MatrixMarket matrix coordinate integer general`; dims
     `nFeatures nBarcodes nEntries`; entries `gene+1 cell+1 count` (1-based),
     iterated in cell-column order.

src/junction/mod.rs

@@ -7,7 +7,7 @@
 /// - Junction statistics collection for SJ.out.tab output
 pub(crate) mod chr_start_end;
 pub(crate) mod gtf;
-mod sj_output;
+pub(crate) mod sj_output;
 pub mod sjdb_insert;
 
 pub use sj_output::SpliceJunctionStats;

src/junction/sj_output.rs

@@ -35,6 +35,15 @@ pub(crate) struct SjKey {
     pub motif: u8, // Encoded motif value
 }
 
+/// One junction's `SJ.out.tab` count columns, snapshotted out of [`SjCounts`].
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct SjRowCounts {
+    pub annotated: bool,
+    pub unique: u32,
+    pub multi: u32,
+    pub max_overhang: u32,
+}
+
 /// Counts for a single junction
 #[derive(Debug)]
 pub(crate) struct SjCounts {
@@ -250,75 +259,82 @@ impl SpliceJunctionStats {
     /// solo-feature matrix. Returns the (intron_start, intron_end) absolute-coord
     /// keys so the SJ recorder can be mapped to matrix rows.
     pub(crate) fn sj_feature_order(&self, params: &Parameters) -> Vec<(u64, u64)> {
-        let surviving = self.compute_surviving_junctions(params);
-        let mut keys: Vec<(usize, u64, u64)> = self
-            .junctions
-            .iter()
-            .filter(|e| surviving.contains(e.key()))
-            .map(|e| {
-                let k = e.key();
-                (k.chr_idx, k.intron_start, k.intron_end)
-            })
-            .collect();
-        keys.sort_unstable();
-        keys.into_iter().map(|(_, s, e)| (s, e)).collect()
+        self.sorted_junctions(params)
+            .into_iter()
+            .map(|(k, _)| (k.intron_start, k.intron_end))
+            .collect()
     }
 
-    /// Write the 9-column `SJ.out.tab` lines (sorted) to `writer`; returns the
-    /// number written. Shared by `write_output` and the SJ feature's
-    /// `features.tsv`, so both stay in the same order as the SJ matrix rows.
-    pub(crate) fn write_sj_lines(
-        &self,
-        writer: &mut dyn std::io::Write,
-        genome: &Genome,
-        params: &Parameters,
-    ) -> Result<u32, Error> {
+    /// The surviving junctions in `SJ.out.tab` order (chromosome, then intron
+    /// start/end), each with its `(annotated, unique, multi, max_overhang)`
+    /// counts. The single source of row order for `SJ.out.tab`, the SJ feature's
+    /// `features.tsv`/`var`, and the SJ matrix rows.
+    pub(crate) fn sorted_junctions(&self, params: &Parameters) -> Vec<(SjKey, SjRowCounts)> {
         let surviving = self.compute_surviving_junctions(params);
-        let mut output_junctions: Vec<_> = self
+        let mut rows: Vec<(SjKey, SjRowCounts)> = self
             .junctions
             .iter()
             .filter(|entry| surviving.contains(entry.key()))
             .map(|entry| {
-                let key = entry.key().clone();
                 let counts = entry.value();
                 (
-                    key,
-                    counts.annotated,
-                    counts.unique_count.load(Ordering::Relaxed),
-                    counts.multi_count.load(Ordering::Relaxed),
-                    counts.max_overhang.load(Ordering::Relaxed),
+                    entry.key().clone(),
+                    SjRowCounts {
+                        annotated: counts.annotated,
+                        unique: counts.unique_count.load(Ordering::Relaxed),
+                        multi: counts.multi_count.load(Ordering::Relaxed),
+                        max_overhang: counts.max_overhang.load(Ordering::Relaxed),
+                    },
                 )
             })
             .collect();
-
-        output_junctions.sort_by(|a, b| {
+        rows.sort_by(|a, b| {
             a.0.chr_idx
                 .cmp(&b.0.chr_idx)
                 .then(a.0.intron_start.cmp(&b.0.intron_start))
                 .then(a.0.intron_end.cmp(&b.0.intron_end))
         });
+        rows
+    }
+
+    /// `(chromosome name, 1-based intron start, 1-based intron end)` for a junction.
+    pub(crate) fn locus<'g>(key: &SjKey, genome: &'g Genome) -> Result<(&'g str, u64, u64), Error> {
+        let chr_name = genome
+            .chr_name
+            .get(key.chr_idx)
+            .ok_or_else(|| Error::Index("Invalid chromosome index in junction".to_string()))?;
+        let chr_start_pos = genome.chr_start[key.chr_idx];
+        Ok((
+            chr_name,
+            key.intron_start - chr_start_pos + 1,
+            key.intron_end - chr_start_pos + 1,
+        ))
+    }
 
+    /// Write the 9-column `SJ.out.tab` lines (sorted) to `writer`; returns the
+    /// number written. Shared by `write_output` and the SJ feature's
+    /// `features.tsv`, so both stay in the same order as the SJ matrix rows.
+    pub(crate) fn write_sj_lines(
+        &self,
+        writer: &mut dyn std::io::Write,
+        genome: &Genome,
+        params: &Parameters,
+    ) -> Result<u32, Error> {
         let mut written = 0u32;
-        for (key, annotated, unique, multi, max_overhang) in &output_junctions {
-            let chr_name = genome
-                .chr_name
-                .get(key.chr_idx)
-                .ok_or_else(|| Error::Index("Invalid chromosome index in junction".to_string()))?;
-            let chr_start_pos = genome.chr_start[key.chr_idx];
-            let chr_pos_start = key.intron_start - chr_start_pos + 1;
-            let chr_pos_end = key.intron_end - chr_start_pos + 1;
+        for (key, c) in self.sorted_junctions(params) {
+            let (chr_name, start, end) = Self::locus(&key, genome)?;
             writeln!(
                 writer,
                 "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
                 chr_name,
-                chr_pos_start,
-                chr_pos_end,
+                start,
+                end,
                 key.strand,
                 key.motif,
-                i32::from(*annotated),
-                unique,
-                multi,
-                max_overhang
+                i32::from(c.annotated),
+                c.unique,
+                c.multi,
+                c.max_overhang
             )
             .map_err(|e| Error::Index(format!("SJ write: {e}")))?;
             written += 1;

src/lib.rs

@@ -437,7 +437,29 @@ fn write_solo_output(
             recorder.n_multi_records(),
         );
     }
-    crate::solo::write_gene_matrix(sctx, params, stats, Some(&**sj_stats), &index.genome)?;
+    // Validated in `Parameters::validate`, so an unknown/unavailable format
+    // never reaches here.
+    match crate::solo::OutputFormat::parse(&params.solo_output_format) {
+        #[cfg(feature = "anndata-out")]
+        Some(crate::solo::OutputFormat::Zarr) => {
+            crate::solo::adata::write_mudata(
+                sctx,
+                params,
+                stats,
+                Some(&**sj_stats),
+                &index.genome,
+            )?;
+        }
+        _ => {
+            crate::solo::write_matrix_market(
+                sctx,
+                params,
+                stats,
+                Some(&**sj_stats),
+                &index.genome,
+            )?;
+        }
+    }
     Ok(())
 }
 

src/params/mod.rs

@@ -1133,6 +1133,16 @@ pub struct Parameters {
     #[arg(long = "soloOutGzip", default_value = "no")]
     pub solo_out_gzip: String,
 
+    /// Container format for the solo count matrices (rustar extension beyond
+    /// STARsolo). `MTX` (default) is STARsolo's `raw/`+`filtered/` MatrixMarket
+    /// triplet. `Zarr` writes one sharded Zarr v3 MuData store,
+    /// `<soloOutFileNames[0]>/<matrix stem>.zarr`, holding every feature as a
+    /// layer of a `gex` modality plus a junction-indexed `sj` modality — cells ×
+    /// features (AnnData's obs × var orientation, the transpose of `.mtx`).
+    /// One format per run: the two writers each consume the count records.
+    #[arg(long = "soloOutputFormat", default_value = "MTX")]
+    pub solo_output_format: String,
+
     /// Velocyto ambiguous-molecule handling (rustar extension beyond STARsolo).
     /// `yes` (default) writes the three `spliced`/`unspliced`/`ambiguous` matrices
     /// like STARsolo — exon-only molecules with no junction/intron evidence stay in
@@ -1642,6 +1652,26 @@ impl Parameters {
                     ));
                 }
             }
+            // soloOutputFormat value (and whether this build has the backend).
+            let f = &params.solo_output_format;
+            match crate::solo::OutputFormat::parse(f) {
+                Some(fmt) if fmt.is_available() => {}
+                Some(fmt) => {
+                    return Err(command.error(
+                        ErrorKind::InvalidValue,
+                        format!(
+                            "--soloOutputFormat {f} needs the `{}` cargo feature, which this binary was built without",
+                            fmt.cargo_feature(),
+                        ),
+                    ));
+                }
+                None => {
+                    return Err(command.error(
+                        ErrorKind::InvalidValue,
+                        format!("unsupported --soloOutputFormat '{f}'; expected MTX or Zarr"),
+                    ));
+                }
+            }
             // Gene-level features need a gene model (SJ does not — junctions come
             // from the alignments).
             let needs_gtf = params