mod column_subset;
mod handlers;
mod hdf5_io;
mod interactive;
mod qc;
mod sparse_backend;
mod sparse_data_visitors;
mod sparse_io;
mod sparse_io_vector;
mod sparse_util;
mod utilities;
#[allow(dead_code)]
mod zarr_io;
use crate::handlers::analysis::{run_histogram, run_stat, RunHistogramArgs, RunStatArgs};
#[cfg(feature = "hdf5")]
use crate::handlers::builders::{
run_build_from_10x_matrix, run_build_from_10x_molecule, run_build_from_h5ad, From10xMatrixArgs,
From10xMoleculeArgs, FromH5adArgs,
};
use crate::handlers::builders::{
run_build_from_fragments, run_build_from_mtx, run_build_from_zarr_triplets, FromFragmentsArgs,
FromMtxArgs, FromZarrArgs,
};
#[cfg(feature = "hdf5")]
use crate::handlers::exporters::{run_export_to_h5ad, ToH5adArgs};
use crate::handlers::exporters::{run_export_to_mtx, ToMtxArgs};
use crate::handlers::inspection::{
show_info, take_column_names, take_columns, take_row_names, take_rows, InfoArgs,
TakeColumnNamesArgs, TakeColumnsArgs, TakeRowNamesArgs, TakeRowsArgs,
};
use crate::handlers::listing::{list_h5, list_zarr, ListH5Args, ListZarrArgs};
use crate::handlers::merging::{
align_backends, run_merge_backend, run_merge_mtx, AlignDataArgs, MergeBackendArgs, MergeMtxArgs,
};
use crate::handlers::transformation::{
reorder_rows, run_convert, run_split, run_squeeze, run_subsample, subset_columns, subset_rows,
ConvertArgs, ReorderRowsArgs, RunSqueezeArgs, SplitArgs, SubsampleArgs, SubsetColumnsArgs,
SubsetRowsArgs,
};
use clap::{Parser, Subcommand};
fn main() -> anyhow::Result<()> {
let cli = Cli::parse();
if cli.verbose {
std::env::set_var("RUST_LOG", "info");
}
env_logger::init();
if let Some(n) = cli.n_threads {
if n == 0 {
anyhow::bail!("--n-threads must be >= 1");
}
rayon::ThreadPoolBuilder::new()
.num_threads(n)
.build_global()
.ok();
}
match run(&cli) {
Err(e) if is_broken_pipe(&e) => {
std::process::exit(0);
}
other => other,
}
}
fn is_broken_pipe(err: &anyhow::Error) -> bool {
err.chain().any(|cause| {
cause
.downcast_ref::<std::io::Error>()
.is_some_and(|io| io.kind() == std::io::ErrorKind::BrokenPipe)
})
}
fn run(cli: &Cli) -> anyhow::Result<()> {
match &cli.commands {
Commands::FromMtx(args) => {
run_build_from_mtx(args)?;
}
Commands::ListH5(args) => {
list_h5(args)?;
}
Commands::ListZarr(args) => {
list_zarr(args)?;
}
Commands::FromZarr(args) => {
run_build_from_zarr_triplets(args)?;
}
#[cfg(feature = "hdf5")]
Commands::From10xMatrix(args) => {
run_build_from_10x_matrix(args)?;
}
#[cfg(feature = "hdf5")]
Commands::FromH5ad(args) => {
run_build_from_h5ad(args)?;
}
#[cfg(feature = "hdf5")]
Commands::From10xMolecule(args) => {
run_build_from_10x_molecule(args)?;
}
#[cfg(feature = "hdf5")]
Commands::ToH5ad(args) => {
run_export_to_h5ad(args)?;
}
Commands::ToMtx(args) => {
run_export_to_mtx(args)?;
}
Commands::Subsample(args) => {
run_subsample(args)?;
}
Commands::Split(args) => {
run_split(args)?;
}
Commands::Convert(args) => {
run_convert(args)?;
}
Commands::FromFragments(args) => {
run_build_from_fragments(args)?;
}
Commands::Info(args) => {
show_info(args)?;
}
Commands::Statistics(args) => {
run_stat(args)?;
}
Commands::Histogram(args) => {
run_histogram(args)?;
}
Commands::Squeeze(args) => {
run_squeeze(args)?;
}
Commands::Columns(args) => {
take_columns(args)?;
}
Commands::Rows(args) => {
take_rows(args)?;
}
Commands::ColumnNames(args) => {
take_column_names(args)?;
}
Commands::RowNames(args) => {
take_row_names(args)?;
}
Commands::SubsetColumns(args) => {
subset_columns(args)?;
}
Commands::SubsetRows(args) => {
subset_rows(args)?;
}
Commands::AlignData(args) => {
align_backends(args)?;
}
Commands::ReorderRows(args) => {
reorder_rows(args)?;
}
Commands::MergeMtx(args) => {
run_merge_mtx(args)?;
}
Commands::MergeBackend(args) => {
run_merge_backend(args)?;
}
}
Ok(())
}
#[derive(Parser, Debug)]
#[command(
version,
about = "Data Backend for Expedited Acquisition and Neighbourhood Search.",
long_about = "Data Backend for Expedited Acquisition and Neighbourhood Search (data-beans).
We assume non-negative sparse matrices were generated by single-cell omics (feature x cell).
This tool creates a data structure for faster access, organized as follows:
(root)
├── nrow
├── ncell
├── by_column
│ ├── data
│ ├── indices (row indices)
│ └── indptr (column pointers)
└── by_row
├── data
├── indices (column indices)
└── indptr (row pointers)
For more details on each command, use '--help' after the command name."
)]
struct Cli {
#[arg(short = 'v', long, global = true)]
verbose: bool,
#[arg(
long = "n-threads",
visible_aliases = ["threads", "num-threads"],
global = true,
value_name = "N",
help = "Limit the number of CPU threads",
long_help = "Limit the global rayon thread pool and HDF5/blosc compression threads.\n\
Useful for bounding peak memory on `from-*` and `merge-*` subcommands,\n\
since each rayon worker may hold its own intermediate buffers.\n\
Defaults to all logical CPUs when unset."
)]
n_threads: Option<usize>,
#[command(subcommand)]
commands: Commands,
}
#[derive(Subcommand, Debug)]
enum Commands {
#[command(about = "Build backend from `mtx` file and associated `tsv` files")]
FromMtx(FromMtxArgs),
#[cfg(feature = "hdf5")]
#[command(
name = "from-10x-matrix",
about = "Build backend from 10X Genomics feature-barcode matrix `h5`",
long_about = "Build a backend from a 10X Genomics Cell Ranger feature-barcode matrix HDF5 file.\n\
Expected layout: matrix/{data,indices,indptr,barcodes,features/...}\n\
\n\
This is Cell Ranger's `filtered_feature_bc_matrix.h5`,\n\
or its `raw_feature_bc_matrix.h5`, from count or multi.\n\
\n\
Use --root-group-name, --data-field, etc. to customize field paths.\n\
Use --select-row-type to filter features (default: 'gene').",
visible_alias = "from-h5"
)]
From10xMatrix(From10xMatrixArgs),
#[cfg(feature = "hdf5")]
#[command(
about = "Build backend from AnnData `h5ad` file (CELLxGENE schema)",
long_about = "Build a backend from an AnnData h5ad file (CELLxGENE schema v7, AnnData spec v0.1.0).\n\
\n\
Auto-detects sparse format (CSR/CSC) and transposes to (features x cells).\n\
Prefers raw/X (raw counts) over X (processed) when available.\n\
\n\
Additionally outputs:\n\
- {output}.cell_metadata.tsv.gz (all obs columns per cell)\n\
- {output}.barcode_to_donor.tsv.gz (barcode-to-donor mapping, if donor_id exists)\n\
- {output}.sample_metadata.tsv.gz (one row per donor, if donor_id exists)\n\
\n\
Column names become barcode@donor_id for multi-donor data.\n\
Use --select-row-type to filter by biotype (e.g., 'protein_coding')."
)]
FromH5ad(FromH5adArgs),
#[cfg(feature = "hdf5")]
#[command(
name = "from-10x-molecule",
about = "Build backend from 10X molecule_info.h5",
long_about = "Build a backend by aggregating per-molecule counts.\n\
The input is a 10X Genomics molecule_info.h5 file.\n\
The output is a feature x cell sparse matrix.\n\
\n\
Each molecule has (barcode_idx, feature_idx, count, gem_group, library_idx).\n\
Molecules are filtered by --library-type. Cell Ranger's pass_filter,\n\
its valid cell calls, can filter too.\n\
Survivors are aggregated into count triplets.\n\
\n\
Barcode names are formatted as SEQUENCE-GEMGROUP (e.g., AAACCTGA-1).\n\
Handles multi-sample (cellranger aggr) and multi-library (CITE-seq) data."
)]
From10xMolecule(From10xMoleculeArgs),
#[cfg(feature = "hdf5")]
#[command(
name = "to-h5ad",
about = "Export backend to an AnnData `h5ad` file (scanpy-readable)",
long_about = "Export a backend to an AnnData h5ad file readable by scanpy / anndata.\n\
This is the inverse of `from-h5ad`.\n\
\n\
The (features x cells) backend is transposed back to AnnData's (cells x features) layout,\n\
and written as a gzip-compressed CSR `X`. Gzip, not Blosc,\n\
so h5py needs no filter plugin.\n\
\n\
Row names map to var/_index; column names to obs/_index.\n\
Attach cell and feature annotations with --obs and --var.\n\
Feeding back the `*.cell_metadata.tsv.gz` that `from-h5ad` emitted is the usual round trip.",
visible_alias = "to-anndata"
)]
ToH5ad(ToH5adArgs),
#[command(
name = "to-mtx",
about = "Export backend to a 10x-style MatrixMarket triplet (MEX)",
long_about = "Export a backend to a 10x Genomics Cell Ranger MEX triplet directory,\n\
loadable with `scanpy.read_10x_mtx(dir)` or Seurat `Read10X(dir)`.\n\
This is the inverse of `from-mtx` and needs no HDF5.\n\
\n\
Writes into the output directory (Cell Ranger v3 gzipped layout):\n\
- matrix.mtx.gz (integer MatrixMarket, features x barcodes, 1-based)\n\
- features.tsv.gz (id, name, feature_type; composite names split on '_')\n\
- barcodes.tsv.gz (one barcode per line)\n\
\n\
Pass --no-gzip for the uncompressed layout.",
visible_alias = "to-10x"
)]
ToMtx(ToMtxArgs),
#[command(
about = "Build backend from triplets in 10X Xenium `zarr`",
long_about = "Build a backend from triplets in `zarr` format.\n\
Supports conversion and indexing for fast access."
)]
FromZarr(FromZarrArgs),
#[command(
name = "from-fragments",
about = "Build ATAC/histone backend from a fragments TSV file",
long_about = "Build a (feature x cell) sparse backend by streaming a TSV.\n\
The input is an scATAC or histone fragments file.\n\
\n\
Expected format (tab-separated, one fragment per line):\n\
chr<TAB>start<TAB>end<TAB>barcode[<TAB>count]\n\
\n\
Both plain gzip and bgzipped files are accepted;\n\
'#' header lines (e.g. cellranger-arc metadata) are skipped.\n\
\n\
Features are user-supplied peaks, via --peaks <bed>.\n\
Otherwise they are fixed-width genome tiles,\n\
discovered on the fly at --bin-size, which defaults to 5000.\n\
\n\
Each fragment contributes 1 to every feature it overlaps.\n\
With --use-count it contributes its column-5 count instead."
)]
FromFragments(FromFragmentsArgs),
#[command(
about = "List contents of `h5` file",
long_about = "List what are included in the `h5` file. Shows datasets, groups,\n\
and metadata."
)]
ListH5(ListH5Args),
#[command(
about = "List contents of `zarr` file",
long_about = "List what are included in the `zarr` file.\n\
Displays structure and available arrays."
)]
ListZarr(ListZarrArgs),
#[command(
about = "Sort rows by name order",
long_about = "Sort rows according to the order of row names specified in a row name file.\n\
Useful for aligning datasets and ensuring consistent row order."
)]
ReorderRows(ReorderRowsArgs),
#[command(
about = "Take columns and output dense matrix",
visible_aliases = ["take-columns"],
long_about = "Take columns from the sparse matrix,\n\
and save them to an `output` file as a dense matrix,\n\
for quick examination.\n\
Useful for extracting subsets for visualization or analysis."
)]
Columns(TakeColumnsArgs),
#[command(
about = "Take rows and output dense matrix (transposed)",
visible_aliases = ["take-rows"],
long_about = "Take rows from the sparse matrix,\n\
and save them to an `output` file as a dense matrix,\n\
for quick examination. For convenience,\n\
it will output a transposed (`column x selected_row`) matrix."
)]
Rows(TakeRowsArgs),
#[command(
about = "List column names",
long_about = "List all column names in the backend.\n\
Useful for inspecting available features."
)]
ColumnNames(TakeColumnNamesArgs),
#[command(
about = "List row names",
long_about = "List all row names in the backend.\n\
Useful for inspecting available samples or observations."
)]
RowNames(TakeRowNamesArgs),
#[command(
about = "Subset columns and create new backend",
long_about = "Take columns from the sparse matrix and create a new sparse matrix backend.\n\
Allows for focused analysis on selected features."
)]
SubsetColumns(SubsetColumnsArgs),
#[command(
about = "Subset rows and create new backend",
long_about = "Take rows from the sparse matrix and create a new sparse matrix backend.\n\
Allows for focused analysis on selected samples/observations."
)]
SubsetRows(SubsetRowsArgs),
#[command(
about = "Randomly subsample cells and/or genes into a smaller backend",
long_about = "Draw a random subset of cells and/or genes into a new backend.\n\
Cells are columns; genes are rows.\n\
This is handy for quick test and demo datasets.\n\
\n\
Specify counts with --cells and --genes.\n\
Or specify fractions with --cell-frac and --gene-frac.\n\
An unset dimension keeps all of its entries.\n\
\n\
Sampling is reproducible via --seed. It reads only the selected columns.\n\
Cost therefore scales with the output size, not the input.",
visible_alias = "downsample"
)]
Subsample(SubsampleArgs),
#[command(
about = "Split cells into train/test halves, or into K cross-validation folds",
long_about = "Split a backend by CELL into disjoint halves, or into K folds.\n\
\n\
Every cell lands on exactly one side, and the halves together\n\
cover the input. Only the selected columns are read,\n\
so the cost scales with the output, not the input.\n\
\n\
On spatial data, split by REGION, not by cell.\n\
A random cell split leaves each test cell ringed by training cells,\n\
and adjacent cells are near-duplicates,\n\
so the model has effectively seen the test cell already.\n\
\x20 --coord positions.csv --coord-columns 4,5 --grid 8\n\
tiles the coordinate bounding box 8x8 and assigns whole tiles,\n\
so a test region gets a training boundary, not a training interior.\n\
\n\
--groups is the escape hatch for a non-geometric grouping:\n\
donor, sample, slide. Two columns: cell name, group label,\n\
as .parquet or delimited text (.tsv/.csv, optionally .gz).\n\
Cells sharing a label stay on the same side. Names must\n\
match EXACTLY -- a fuzzy match would put a cell in the\n\
wrong half with nothing downstream able to detect it.\n\
\n\
Outputs:\n\
\x20 {out}.train.zarr.zip and {out}.test.zarr.zip\n\
\x20 or {out}.fold{k}.train / .test under --folds\n\
\n\
THE WHOLE SCHEME, end to end:\n\
\x20 1. data-beans split data.zarr -o cv --test-frac 0.2\n\
\x20 (add --coord/--grid on spatial data; --folds K for CV)\n\
\x20 2. train each method on cv.train.zarr.zip\n\
\x20 3. senna predict cv.test.zarr.zip --model M -o pred\n\
\x20 pinto predict cv.test.zarr.zip --model M -o pred\n\
\x20 4. compare pred.predictive.parquet across methods, on\n\
\x20 eval_llik_per_count minus eval_null_llik_per_count\n\
\x20 (pinto: llik_per_count minus null_llik_per_count).\n\
\x20 NOT on the plain llik column -- that one is the\n\
\x20 backend's own and differs by decoder. Filter the\n\
\x20 count column > 0 first; empty rows carry NaN.\n\
\n\
To ablate features, do it at PREDICT, not here:\n\
\x20 senna predict ... --ablate-features hide.txt\n\
hides those genes from the encoder and scores on them, so\n\
the score is a prediction rather than a reconstruction.\n\
Ablating the TRAINING half instead (subsample --gene-frac)\n\
answers a different question -- how much the dictionary\n\
needed those genes -- and does NOT remove the advantage a\n\
larger latent has when it is fitted on the cell it scores.\n\
\n\
Selection: with 10+ arms on one test half, the winner's\n\
score is biased upward by the max over arms. Use --folds,\n\
pick the arm on one fold, and report on another."
)]
Split(SplitArgs),
#[command(
about = "Convert a backend between on-disk formats (zarr <-> h5)",
long_about = "Convert a backend to a different on-disk format,\n\
preserving the matrix and row/column names:\n\
- zarr <-> h5 (re-encodes the data)\n\
- .zarr <-> .zarr.zip (unzip / re-zip)\n\
\n\
The output format follows --backend and the output path (pass\n\
--no-zip to keep a .zarr directory instead of a .zarr.zip archive).",
visible_alias = "to-backend"
)]
Convert(ConvertArgs),
#[command(
about = "Align data backends",
long_about = "To ensure that column names are aligned for multimodal analysis.\n\
We will only keep columns and rows matched across files. 1st row:\n\
`D(1,1)-D(1,2)`, 2nd row: `D(2,1)-D(2,2)`, etc.",
visible_alias = "align"
)]
AlignData(AlignDataArgs),
#[command(
about = "Merge multiple `.mtx` files",
long_about = "Merge multiple 10x `.mtx` files into one fileset.\n\
Useful for combining datasets from different sources."
)]
MergeMtx(MergeMtxArgs),
#[command(
about = "Merge multiple backend files",
long_about = "Merge multiple backend file(sets) into one.\n\
Supports various formats and options for merging.",
visible_alias = "merge"
)]
MergeBackend(MergeBackendArgs),
#[command(
about = "Squeeze out sparse rows/columns",
long_about = "Squeeze out rows and columns with too few non-zeros.\n\
It will overwrite the original (be careful) and save the indices kept."
)]
Squeeze(RunSqueezeArgs),
#[command(
about = "Show basic matrix info",
long_about = "Show basic information of a sparse matrix. If output header is provided,\n\
row and column names will be saved."
)]
Info(InfoArgs),
#[command(
about = "Take matrix statistics",
long_about = "Take basic statistics from a sparse matrix.\n\
The output file will contain columns of (1) `nnz` - number of non-zero elements,\n\
(2) `tot` - total sum, (3) `mu` - average `μ`,\n\
(4) `sig` - standard deviation `σ`.",
visible_alias = "stat"
)]
Statistics(RunStatArgs),
#[command(
about = "ASCII log-scale histogram of a row/column statistic",
long_about = "Print an ASCII log10(x+1) histogram of a statistic.\n\
It covers the per-feature (row) and per-cell (column) axes.\n\
This is the summary `squeeze --show-histogram` shows,\n\
exposed as a standalone command.\n\
\n\
Choose the statistic with --stat, default `nnz`, also `sum`, `mean`, `sd`;\n\
and the margin with --dim, default `both`.\n\
Pass -o/--output to also dump the raw per-unit values.",
visible_alias = "hist"
)]
Histogram(RunHistogramArgs),
}