Refactor/output direct to h5ad

CHANGELOG.md

@@ -0,0 +1,51 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [Unreleased]
+
+### Changed
+
+- `cyto workflow gex`/`crispr` (and `cyto ibu count --h5ad`) now write `.h5ad` output
+  natively from Rust instead of staging through MTX and shelling out to `pycyto convert`.
+  `cyto-ibu-count` builds the sparse count matrix directly from its in-memory dedup
+  table and writes it via the `anndata`/`anndata-hdf5` crates, so h5ad output no longer
+  round-trips through on-disk MTX text files or spawns a Python subprocess.
+- `cell-filter` (GEX droplet filtering) and `geomux` (CRISPR guide assignment) are
+  unchanged and still run as external Python tools via `uv`, invoked on the natively
+  written `.h5ad`.
+- `pycyto` is no longer installed or invoked anywhere in `cyto workflow`. It remains
+  useful standalone for aggregating multiple `cyto` outputs into one `h5ad`
+  (`pycyto aggregate`), documented in the README.
+- The h5ad `X` matrix is now stored as unsigned 32-bit integers rather than the
+  float32 that `pycyto convert` produced by default. Both `cell-filter` and `geomux`
+  already cast counts back to `int` internally, so this has no behavioral effect
+  downstream.
+
+### Added
+
+- `cyto ibu count --h5ad`: writes a native AnnData `.h5ad` file directly (mutually
+  exclusive with `--mtx`).
+- New `cyto-ibu-count` dependencies: `anndata`, `anndata-hdf5` (statically builds
+  `libhdf5` from source via `hdf5-metno-sys`, requires `cmake` at build time — no
+  system `libhdf5` needed at runtime), and `nalgebra-sparse`.
+
+### Fixed
+
+- A file-locking race: `cyto-workflow` processes probes concurrently via Rayon, and
+  HDF5's advisory file locking could spuriously report a just-written, already-closed
+  `.h5ad` as locked while a sibling probe's `.h5ad` was still open elsewhere in the
+  same process. `cell-filter`/`geomux` subprocess invocations now set
+  `HDF5_USE_FILE_LOCKING=FALSE` to avoid this.
+
+### Removed
+
+- `convert_to_h5ad()` and the `pycyto` external-tool install path
+  (`VERSION_PYCYTO`, the `uv tool install pycyto` step in
+  `ArgsWorkflow::validate_requirements()`).
+- `Module::ConversionH5ad` timing entry (h5ad write time is now folded into the
+  `Counting` step) and `ArgsWorkflow::mtx()` (no longer needed now that h5ad output
+  doesn't stage through MTX).

README.md

@@ -434,10 +434,10 @@ ACGTACGT   ENSG00000000005   12
 
 #### Matrix Market Format
 
-For downstream analysis with scanpy/Seurat:
+For downstream analysis from MTX with scanpy/Seurat:
 
 ```bash
-cyto ibu count -i sample.ibu -f features.tsv -o counts_mtx --format mtx
+cyto ibu count -i sample.ibu -f features.tsv -o counts_mtx --mtx
 ```
 
 Generates:
@@ -446,15 +446,17 @@ Generates:
 - `barcodes.tsv` - Cell barcodes
 - `features.tsv` - Feature names
 
-#### Convert to h5ad
+#### h5ad Format
 
-Use [pycyto](https://github.com/arcinstitute/pycyto) utilities for format conversion and aggregation:
+Writes a native AnnData file directly without MTX staging. This is the default format for `cyto workflow`:
 
 ```bash
-# Convert MTX to h5ad
-pycyto mtx-to-h5ad counts_mtx/ output.h5ad
+cyto ibu count -i sample.ibu -f features.tsv -o output.h5ad --h5ad
+```
+
+To aggregate multiple cyto outputs into a single h5ad per sample, use [pycyto](https://github.com/arcinstitute/pycyto):
 
-# Aggregate cyto output into a single h5ad per sample
+```bash
 pycyto aggregate <config>.json <cyto_output_dir> <aggr_dir>
 ```
 

crates/cyto-cli/CLAUDE.md

@@ -31,8 +31,9 @@ Defines all CLI argument structures using Clap. This crate is a pure definition
 
 - Geometry presets: V2 presets force `remap_window=5`, V1 uses default of 1
 - `MultiPairedInput.is_binseq()` auto-detects format by file extension
-- `ArgsWorkflow.validate_requirements()` transparently installs Python tools (`pycyto`, `cell-filter`, `geomux`) via `uv tool install` at pinned versions
-- External tool versions are pinned as constants: `VERSION_GEOMUX`, `VERSION_CELL_FILTER`, `VERSION_PYCYTO`
+- `ArgsWorkflow.validate_requirements()` transparently installs Python tools (`cell-filter`, `geomux`) via `uv tool install` at pinned versions, only when format is h5ad and the corresponding step isn't skipped. h5ad itself is written natively by `cyto-ibu-count` — no Python tool required.
+- External tool versions are pinned as constants: `VERSION_GEOMUX`, `VERSION_CELL_FILTER`
+- `ArgsCount` (`src/ibu/count.rs`) has mutually-exclusive `mtx`/`h5ad` flags (`conflicts_with` each other); `from_wf_path()` takes a `CountFormat` and sets the corresponding bool
 
 ## Dependencies (within workspace)
 

crates/cyto-cli/src/ibu/count.rs

@@ -1,5 +1,7 @@
 use std::path::Path;
 
+use crate::workflow::CountFormat;
+
 use super::IbuInput;
 
 #[derive(clap::Parser, Debug)]
@@ -17,9 +19,25 @@ pub struct ArgsCount {
     /// (1) barcodes.txt.gz
     /// (2) features.txt.gz
     /// (3) matrix.mtx.gz
-    #[clap(long, requires = "output", requires = "features")]
+    #[clap(
+        long,
+        requires = "output",
+        requires = "features",
+        conflicts_with = "h5ad"
+    )]
     pub mtx: bool,
 
+    /// Output h5ad format directly.
+    ///
+    /// Writes a native .h5ad `AnnData` file
+    #[clap(
+        long,
+        requires = "output",
+        requires = "features",
+        conflicts_with = "mtx"
+    )]
+    pub h5ad: bool,
+
     /// Number of threads to use in counting
     #[clap(short = 't', long, default_value = "1")]
     pub num_threads: usize,
@@ -50,14 +68,15 @@ impl ArgsCount {
         out_path: P,
         features_path: P,
         num_threads: usize,
-        mtx: bool,
+        format: CountFormat,
         suffix: Option<String>,
     ) -> Self {
         Self {
             input: IbuInput::from_path(sort_path),
             output: Some(out_path.as_ref().to_str().unwrap().to_string()),
             features: Some(features_path.as_ref().to_str().unwrap().to_string()),
-            mtx,
+            mtx: format == CountFormat::Mtx,
+            h5ad: format == CountFormat::H5ad,
             compressed: false,
             feature_col: 1,
             num_threads,

crates/cyto-cli/src/workflow/mod.rs

@@ -8,7 +8,6 @@ use crate::{ArgsCrispr, ArgsGex};
 
 pub const VERSION_GEOMUX: &str = "0.5.5";
 pub const VERSION_CELL_FILTER: &str = "0.1.2";
-pub const VERSION_PYCYTO: &str = "0.1.14";
 
 #[derive(Subcommand, Debug)]
 pub enum WorkflowCommand {
@@ -126,8 +125,19 @@ pub struct ArgsWorkflow {
     pub format: CountFormat,
 }
 impl ArgsWorkflow {
+    /// Checks for the external Python tools needed to post-process an h5ad file.
+    ///
+    /// h5ad itself is written natively (no external tool required); `cell-filter`
+    /// (GEX droplet filtering) and `geomux` (CRISPR guide assignment) are only
+    /// needed when producing h5ad output and their respective step isn't skipped.
     pub fn validate_requirements(&self, mode: WorkflowMode) -> Result<()> {
-        if self.format == CountFormat::H5ad || !self.no_filter {
+        let needs_cell_filter =
+            mode == WorkflowMode::Gex && self.format == CountFormat::H5ad && !self.no_filter;
+        let needs_geomux = mode == WorkflowMode::Crispr
+            && self.format == CountFormat::H5ad
+            && !self.skip_assignment;
+
+        if needs_cell_filter || needs_geomux {
             debug!("Checking if `uv` exists in $PATH");
             match Command::new("uv").args(["--version"]).output() {
                 Ok(_) => debug!("Found `uv` in $PATH"),
@@ -136,27 +146,16 @@ impl ArgsWorkflow {
                     bail!("Encountered an unexpected error checking for `uv`: {e}");
                 }
             }
-            transparent_uv_install("pycyto", VERSION_PYCYTO)?;
         }
-        if mode == WorkflowMode::Gex && !self.no_filter {
+        if needs_cell_filter {
             transparent_uv_install("cell-filter", VERSION_CELL_FILTER)?;
         }
-        if mode == WorkflowMode::Crispr {
+        if needs_geomux {
             transparent_uv_install("geomux", VERSION_GEOMUX)?;
         }
         Ok(())
     }
 
-    /// Check whether the workflow should output mtx files
-    ///
-    /// This is true if the format is mtx or h5ad but mtx is consumed by h5ad
-    pub fn mtx(&self) -> bool {
-        match self.format {
-            CountFormat::H5ad | CountFormat::Mtx => true,
-            CountFormat::Tsv => false,
-        }
-    }
-
     /// Check whether the workflow should output h5ad files
     pub fn to_h5ad(&self) -> bool {
         match self.format {

crates/cyto-ibu-count/CLAUDE.md

@@ -2,7 +2,7 @@
 
 ## Purpose
 
-Creates barcode-by-feature count matrices from sorted IBU files. Deduplicates UMIs by selecting the most abundant index per UMI (discarding ties), optionally aggregates counts at a higher level (e.g., gene level from probe level), and writes output in TSV or MTX format.
+Creates barcode-by-feature count matrices from sorted IBU files. Deduplicates UMIs by selecting the most abundant index per UMI (discarding ties), optionally aggregates counts at a higher level (e.g., gene level from probe level), and writes output in TSV, MTX, or h5ad format.
 
 ## Key Source Files
 
@@ -14,11 +14,13 @@ Creates barcode-by-feature count matrices from sorted IBU files. Deduplicates UM
   - `DeduplicateError` — Custom errors: `UnsortedIbu`, `MaxIndexExceeded`, `EmptyStream`
   - Extensive unit tests (~400 lines) covering deduplication edge cases
 - `src/lib.rs` — Output formatting and aggregation:
-  - `run()` — Entry point: loads features, deduplicates, optionally aggregates, writes output
+  - `run()` — Entry point: loads features, deduplicates, optionally aggregates, dispatches to h5ad/mtx/tsv writer based on `ArgsCount.h5ad`/`.mtx`
   - `aggregate_unit()` — Aggregates counts by feature name (e.g., probe -> gene level)
   - `write_counts_tsv()` — TSV output (encoded 2-bit or decoded nucleotide barcodes)
   - `write_counts_mtx()` — Matrix Market format: `matrix.mtx.gz`, `barcodes.tsv.gz`, `features.tsv.gz` (parallel gzip via `gzp`)
   - `load_features()` — Reads feature names from TSV at specified column index
+- `src/h5ad.rs` — Native AnnData output:
+  - `write_counts_h5ad()` — Builds a CSR sparse matrix directly from `BarcodeIndexCounts` (no MTX intermediate) and writes it as a `.h5ad` file via the `anndata`/`anndata-hdf5` crates. Forces gzip compression (instead of the crate's zstd default) so the file is readable by plain `h5py`/`anndata` without the optional `hdf5plugin` package.
 
 ## Key Types
 
@@ -30,6 +32,7 @@ Creates barcode-by-feature count matrices from sorted IBU files. Deduplicates UM
 
 - Barcodes can be output as 2-bit encoded integers (`--compressed`) or decoded nucleotide strings (default)
 - MTX output uses parallel gzip compression via `gzp::ParCompress`
+- h5ad output writes directly from the in-memory count table — no MTX staging and no external Python conversion tool (previously `pycyto convert`)
 - Barcode suffix support (e.g., `-ProbeA`) for multiplexed experiments
 - Feature aggregation happens post-deduplication: probe-level counts are summed to gene-level
 
@@ -38,6 +41,11 @@ Creates barcode-by-feature count matrices from sorted IBU files. Deduplicates UM
 - `cyto-cli` — `ArgsCount` argument struct
 - `cyto-io` — I/O handle creation
 
+## External Dependencies
+
+- `anndata` / `anndata-hdf5` — Native h5ad read/write; statically builds `libhdf5` from source via `hdf5-metno-sys` (requires `cmake` at build time, no system libhdf5 needed at runtime)
+- `nalgebra-sparse` — CSR/COO sparse matrix construction for h5ad's `X`
+
 ## Testing
 
 ```bash

crates/cyto-ibu-count/Cargo.toml

@@ -20,7 +20,10 @@ log = { workspace = true }
 hashbrown = { workspace = true }
 serde = { workspace = true }
 
+anndata = "0.7"
+anndata-hdf5 = "0.5"
 gzp = "2.0.2"
+nalgebra-sparse = "0.11"
 thiserror = "2.0.18"
 
 [lints]

crates/cyto-ibu-count/src/h5ad.rs

@@ -0,0 +1,86 @@
+use std::path::Path;
+
+use anndata::backend::{Compression, WriteConfig, set_default_write_config};
+use anndata::data::DataFrameIndex;
+use anndata::{AnnData, AnnDataOp};
+use anndata_hdf5::H5;
+use anyhow::Result;
+use hashbrown::HashMap;
+use ibu::Header;
+use log::info;
+use nalgebra_sparse::coo::CooMatrix;
+use nalgebra_sparse::csr::CsrMatrix;
+
+use crate::dedup::BarcodeIndexCounts;
+use crate::extend_suffix;
+
+/// Writes a barcode-by-feature count matrix directly to a `.h5ad` file.
+///
+/// Uses gzip rather than this crate's zstd default so the file is readable by
+/// plain `h5py`/`anndata` without requiring the optional `hdf5plugin` package.
+pub fn write_counts_h5ad<P: AsRef<Path>>(
+    path: P,
+    counts: &BarcodeIndexCounts,
+    features: &[String],
+    header: Header,
+    suffix: Option<&str>,
+) -> Result<()> {
+    set_default_write_config(WriteConfig {
+        compression: Some(Compression::Gzip(6)),
+        block_size: None,
+    });
+
+    let n_var = features.len();
+    let n_obs = counts.get_num_barcodes();
+    let nnz = counts.get_nnz();
+
+    let mut obs_names = Vec::with_capacity(n_obs);
+    let mut bc_idx_map = HashMap::with_capacity(n_obs);
+    let mut row_idx = Vec::with_capacity(nnz);
+    let mut col_idx = Vec::with_capacity(nnz);
+    let mut vals: Vec<u32> = Vec::with_capacity(nnz);
+    let mut dbuf = Vec::default();
+
+    for record in counts.iter_counts() {
+        let obs_idx = if let Some(idx) = bc_idx_map.get(&record.barcode()) {
+            *idx
+        } else {
+            dbuf.clear();
+            bitnuc::from_2bit(record.barcode(), header.bc_len as usize, &mut dbuf)?;
+            extend_suffix(&mut dbuf, suffix);
+
+            let obs_idx = obs_names.len();
+            obs_names.push(std::str::from_utf8(&dbuf)?.to_string());
+            bc_idx_map.insert(record.barcode(), obs_idx);
+            obs_idx
+        };
+
+        row_idx.push(obs_idx);
+        col_idx.push(record.index() as usize);
+        vals.push(record.count() as u32);
+    }
+
+    let coo = CooMatrix::try_from_triplets(n_obs, n_var, row_idx, col_idx, vals)
+        .map_err(|e| anyhow::anyhow!("Unable to build sparse count matrix: {e}"))?;
+    let csr = CsrMatrix::from(&coo);
+
+    if let Some(parent) = path.as_ref().parent() {
+        std::fs::create_dir_all(parent)?;
+    }
+    if path.as_ref().exists() {
+        std::fs::remove_file(path.as_ref())?;
+    }
+
+    let adata = AnnData::<H5>::new(path.as_ref())?;
+    adata.set_x(csr)?;
+    adata.set_obs_names(DataFrameIndex::from(obs_names))?;
+    adata.set_var_names(DataFrameIndex::from(features.to_vec()))?;
+    adata.close()?;
+
+    info!(
+        "Finished writing h5ad counts to {}",
+        path.as_ref().display()
+    );
+
+    Ok(())
+}

crates/cyto-ibu-count/src/lib.rs

@@ -13,7 +13,9 @@ use ibu::{Header, Reader};
 use log::{debug, error, info};
 
 mod dedup;
+mod h5ad;
 pub use dedup::{BarcodeIndexCount, BarcodeIndexCounts, deduplicate_umis};
+use h5ad::write_counts_h5ad;
 
 /// Extends a barcode buffer with an optional suffix
 fn extend_suffix(buffer: &mut Vec<u8>, suffix: Option<&str>) {
@@ -326,7 +328,19 @@ pub fn run(args: &ArgsCount) -> Result<()> {
         }
     }
 
-    if args.mtx {
+    if args.h5ad {
+        write_counts_h5ad(
+            args.output
+                .as_ref()
+                .expect("Must provide an output path to write h5ad"),
+            &counts,
+            features
+                .expect("Must provide a feature file to write h5ad")
+                .as_slice(),
+            header,
+            args.suffix.as_deref(),
+        )
+    } else if args.mtx {
         write_counts_mtx(
             args.output
                 .as_ref()