diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..571f2f2 --- /dev/null +++ b/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). diff --git a/README.md b/README.md index 9413578..ed97096 100644 --- a/README.md +++ b/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 .json ``` diff --git a/crates/cyto-cli/CLAUDE.md b/crates/cyto-cli/CLAUDE.md index b791839..7095572 100644 --- a/crates/cyto-cli/CLAUDE.md +++ b/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) diff --git a/crates/cyto-cli/src/ibu/count.rs b/crates/cyto-cli/src/ibu/count.rs index 9d0fbdf..366173c 100644 --- a/crates/cyto-cli/src/ibu/count.rs +++ b/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, ) -> 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, diff --git a/crates/cyto-cli/src/workflow/mod.rs b/crates/cyto-cli/src/workflow/mod.rs index 58531b6..1776327 100644 --- a/crates/cyto-cli/src/workflow/mod.rs +++ b/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 { diff --git a/crates/cyto-ibu-count/CLAUDE.md b/crates/cyto-ibu-count/CLAUDE.md index 4a4b2b5..3f14ee7 100644 --- a/crates/cyto-ibu-count/CLAUDE.md +++ b/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 diff --git a/crates/cyto-ibu-count/Cargo.toml b/crates/cyto-ibu-count/Cargo.toml index ade126d..7c71310 100644 --- a/crates/cyto-ibu-count/Cargo.toml +++ b/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] diff --git a/crates/cyto-ibu-count/src/h5ad.rs b/crates/cyto-ibu-count/src/h5ad.rs new file mode 100644 index 0000000..b03576a --- /dev/null +++ b/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>( + 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 = 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::
::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(()) +} diff --git a/crates/cyto-ibu-count/src/lib.rs b/crates/cyto-ibu-count/src/lib.rs index 9427581..f25efca 100644 --- a/crates/cyto-ibu-count/src/lib.rs +++ b/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, 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() diff --git a/crates/cyto-map/src/detect.rs b/crates/cyto-map/src/detect.rs index a6577dc..ac5d027 100644 --- a/crates/cyto-map/src/detect.rs +++ b/crates/cyto-map/src/detect.rs @@ -128,9 +128,7 @@ impl PositionAccumulator { let hi = best_pos.saturating_add(window); self.counts .iter() - .filter(|((c, m, p), _)| { - *c == component && *m == mate && *p >= lo && *p <= hi - }) + .filter(|((c, m, p), _)| *c == component && *m == mate && *p >= lo && *p <= hi) .map(|(_, &count)| count) .sum() } @@ -2367,10 +2365,7 @@ mod tests { // Per-file W=1: [65,67] on file A's acc = 0 + 5000 + 0 = 5000. windowed_match_count: 5000, windowed_match_proportion: 0.5, - top_positions: vec![ - (ReadMate::R2, 66, 5000), - (ReadMate::R2, 68, 500), - ], + top_positions: vec![(ReadMate::R2, 66, 5000), (ReadMate::R2, 68, 500)], }, ], total_reads_sampled: 10_000, @@ -2433,7 +2428,10 @@ mod tests { ]) .unwrap(); - assert_eq!(aggregated.remap_window, 2, "max_remap_window = max(1,2) = 2"); + assert_eq!( + aggregated.remap_window, 2, + "max_remap_window = max(1,2) = 2" + ); assert_eq!(aggregated.total_reads_sampled, 20_000); // Merged accumulator at R2: 66->8000, 67->2000, 68->2500. diff --git a/crates/cyto-workflow/CLAUDE.md b/crates/cyto-workflow/CLAUDE.md index 4183ec1..74459a4 100644 --- a/crates/cyto-workflow/CLAUDE.md +++ b/crates/cyto-workflow/CLAUDE.md @@ -2,20 +2,19 @@ ## Purpose -Orchestrates end-to-end analysis pipelines. Runs the full sequence: map -> sort -> umi-correct -> reads -> count -> convert -> filter/assign. Parallelizes post-mapping steps across probes using Rayon. Invokes external Python tools (`pycyto`, `cell-filter`, `geomux`) via `std::process::Command`. +Orchestrates end-to-end analysis pipelines. Runs the full sequence: map -> sort -> umi-correct -> reads -> count -> filter/assign. Parallelizes post-mapping steps across probes using Rayon. h5ad is written natively by `cyto-ibu-count` (no external tool); `cell-filter` and `geomux` are still invoked via `std::process::Command` for GEX droplet filtering and CRISPR guide assignment. ## Key Source Files - `src/gex.rs` — `run()`: GEX workflow entry point. Calls `cyto_map::run_gex()`, then parallelizes `ibu_steps()` across all per-probe IBU files. Distributes threads proportionally across files. - `src/crispr.rs` — `run()`: CRISPR workflow entry point. Same structure as GEX but passes `ArgsGeomux` for guide assignment step. - `src/utils.rs` — Core workflow utilities: - - `ibu_steps()` — Orchestrates per-IBU pipeline: sort -> umi-correct (optional) -> reads stats (optional) -> count -> h5ad conversion (optional) -> filter/assign. Cleans up intermediate files. + - `ibu_steps()` — Orchestrates per-IBU pipeline: sort -> umi-correct (optional) -> reads stats (optional) -> count (writes h5ad/mtx/tsv directly based on format) -> filter/assign. Cleans up intermediate files. - `identify_ibu_files()` — Globs `outdir/ibu/*.ibu`, excludes `.sort.ibu` - - `convert_to_h5ad()` — Calls `pycyto convert`, removes MTX directory on success - `filter_h5ad()` — Calls `cell-filter` (EmptyDrops), handles missing filtered output gracefully - `assign_guides()` — Calls `geomux` with full parameter passthrough, handles known warning conditions - `write_done_file()` / `write_timings_file()` — Writes workflow completion marker and timing TSV -- `src/timing.rs` — `ModuleTiming` (ibu_name, module, elapsed_secs), `Module` enum (Mapping, InitialSort, UmiCorrection, ReadsDump, Counting, ConversionH5ad, DropletFiltering, GuideAssignment) +- `src/timing.rs` — `ModuleTiming` (ibu_name, module, elapsed_secs), `Module` enum (Mapping, InitialSort, UmiCorrection, ReadsDump, Counting, DropletFiltering, GuideAssignment) ## Key Types @@ -28,6 +27,7 @@ Orchestrates end-to-end analysis pipelines. Runs the full sequence: map -> sort - Thread distribution: total threads divided evenly across IBU files, minimum 1 per file - Intermediate IBU files are removed as they're consumed (unsorted -> sorted -> umi-corrected) - External tool errors surface as `bail!()` with stdout/stderr logged, except for known warnings (e.g., "No guides passed the cell threshold") which are logged as warnings +- `filter_h5ad()`/`assign_guides()` set `HDF5_USE_FILE_LOCKING=FALSE` on the `cell-filter`/`geomux` subprocess env: since probes are counted concurrently via Rayon within this same process, HDF5's advisory file locking can spuriously report a just-written, already-closed h5ad file as locked when a sibling probe's h5ad is still open elsewhere in-process - Output directory structure: `ibu/`, `counts/`, `stats/`, `assignments/`, `metadata/` ## Dependencies (within workspace) diff --git a/crates/cyto-workflow/src/timing.rs b/crates/cyto-workflow/src/timing.rs index 9dac64e..68a01fb 100644 --- a/crates/cyto-workflow/src/timing.rs +++ b/crates/cyto-workflow/src/timing.rs @@ -26,7 +26,6 @@ pub enum Module { UmiCorrection, ReadsDump, Counting, - ConversionH5ad, DropletFiltering, GuideAssignment, } @@ -38,7 +37,6 @@ impl Display for Module { Module::UmiCorrection => write!(f, "UmiCorrection"), Module::ReadsDump => write!(f, "ReadsDump"), Module::Counting => write!(f, "Counting"), - Module::ConversionH5ad => write!(f, "ConversionH5ad"), Module::DropletFiltering => write!(f, "DropletFiltering"), Module::GuideAssignment => write!(f, "GuideAssignment"), } diff --git a/crates/cyto-workflow/src/utils.rs b/crates/cyto-workflow/src/utils.rs index c5d9f60..8954f94 100644 --- a/crates/cyto-workflow/src/utils.rs +++ b/crates/cyto-workflow/src/utils.rs @@ -6,7 +6,7 @@ use std::time::Instant; use anyhow::bail; use anyhow::{Context, Result}; use cyto_cli::ibu::ArgsReads; -use cyto_cli::workflow::{ArgsGeomux, CrisprMappingCommand, GexMappingCommand}; +use cyto_cli::workflow::{ArgsGeomux, CountFormat, CrisprMappingCommand, GexMappingCommand}; use cyto_cli::{ ibu::{ArgsCount, ArgsSort, ArgsUmi}, workflow::{ArgsWorkflow, WorkflowMode}, @@ -42,44 +42,6 @@ fn strip_ibu_basename(ibu_path: &str) -> Result<&str> { Ok(base_ibu_path) } -fn convert_to_h5ad>(count_path: P) -> Result<()> { - info!( - "Converting MTX {} -> {}.h5ad", - count_path.as_ref().display(), - count_path.as_ref().display() - ); - - let output = Command::new("pycyto") - .arg("convert") - .arg(count_path.as_ref().display().to_string()) - .arg(format!("{}.h5ad", count_path.as_ref().display())) - .output()?; - if output.status.success() { - debug!( - "Successfully converted {} to h5ad", - count_path.as_ref().display() - ); - debug!("Removing MTX directory"); - std::fs::remove_dir_all(&count_path).context(format!( - "Unable to remove directory {}", - count_path.as_ref().display() - ))?; - } else { - error!( - "Unable to run h5ad conversion for {}", - count_path.as_ref().display() - ); - error!("stdout: {}", std::str::from_utf8(&output.stdout)?); - error!("stderr: {}", std::str::from_utf8(&output.stderr)?); - bail!( - "Unable to convert {} to h5ad", - count_path.as_ref().display() - ); - } - - Ok(()) -} - fn filter_h5ad>( count_path: P, stats_outdir: P, @@ -96,6 +58,10 @@ fn filter_h5ad>( .arg(&out_h5ad) .arg("--logfile") .arg(logfile) + // Our own process writes h5ad files concurrently across probes via Rayon; + // HDF5's advisory file locking can spuriously report a still-open file as + // locked under that concurrency, so disable it for the reading process. + .env("HDF5_USE_FILE_LOCKING", "FALSE") .output() .context("Unable to run cell-filter")?; @@ -184,6 +150,9 @@ pub fn assign_guides>( } let output = Command::new("geomux") .args(&geomux_args_vec) + // See the comment in `filter_h5ad`: avoids spurious lock errors from + // concurrent h5ad writes across probes within our own process. + .env("HDF5_USE_FILE_LOCKING", "FALSE") .output() .context("Unable to run geomux")?; @@ -298,13 +267,16 @@ pub fn ibu_steps>( // Locate the expected feature path let feature_path = outdir.as_ref().join("metadata").join("features.tsv"); // Build the expected count path - let count_path = if wf_args.mtx() { - outdir.as_ref().join("counts").join(base_ibu_path) - } else { - outdir + let count_path = match wf_args.format { + CountFormat::Mtx => outdir.as_ref().join("counts").join(base_ibu_path), + CountFormat::H5ad => outdir .as_ref() .join("counts") - .join(format!("{base_ibu_path}.counts.tsv")) + .join(format!("{base_ibu_path}.h5ad")), + CountFormat::Tsv => outdir + .as_ref() + .join("counts") + .join(format!("{base_ibu_path}.counts.tsv")), }; // Create the argument struct let count_args = ArgsCount::from_wf_path( @@ -312,7 +284,7 @@ pub fn ibu_steps>( &count_path, &feature_path, 1, - wf_args.mtx(), + wf_args.format, if base_ibu_path == DEFAULT_OUTPUT_BASENAME { None } else { @@ -320,7 +292,7 @@ pub fn ibu_steps>( }, ); - // Run the counting step + // Run the counting step (writes h5ad directly when format is h5ad) info!("Counting {sort_path} -> {}", count_path.display()); let start = Instant::now(); cyto_ibu_count::run(&count_args)?; @@ -332,17 +304,8 @@ pub fn ibu_steps>( std::fs::remove_file(&sort_path).context("Unable to remove IBU file")?; } - // Convert to h5ad if required + // Post-process the h5ad if required if wf_args.to_h5ad() { - let start = Instant::now(); - convert_to_h5ad(&count_path)?; - let elapsed = start.elapsed(); - timings.push(ModuleTiming::new( - base_ibu_path, - Module::ConversionH5ad, - elapsed, - )); - match wf_mode { WorkflowMode::Gex => { if !wf_args.no_filter {