API reference 

residual_effects(data, perturbation_column, control_label=None, *, covariates, null='model', min_cells=10, min_control_mean=0.1, perturbations=None, gene_name_column=None, chunk_size=None, n_jobs=None, memory_limit_gb=None, format_mismatch_policy='auto', data_name=None, output_path=None, verbose=True, force=False)[source]

Perturbation effects adjusted for covariates by a model fitted on controls.

See crispyx.residual.residual_effects() for the full parameter documentation.

Parameters:

data (str | Path | AnnData)
perturbation_column (str | Series)
control_label (str | None)
covariates (Sequence[str])
null (Literal['model', 'robust'])
min_cells (int)
min_control_mean (float)
perturbations (Iterable[str] | None)
gene_name_column (str | None)
chunk_size (int | None)
n_jobs (int | None)
memory_limit_gb (float | None)
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off'])
data_name (str | None)
output_path (str | Path | None)
verbose (int | bool)
force (bool)

Return type:

AnnData

compute_overlap(sets_dict, *, metric='both')[source]: Compute pairwise overlap statistics. See crispyx.compute_overlap().

estimate_disk_usage(func, data, **kwargs)[source]: Estimate the disk space a crispyx function will need. See crispyx.estimate_disk_usage().

Plotting (`cx.pl`)

class crispyx._namespaces._PlottingNamespace[source]

Scanpy-style plotting entry points (cx.pl).

rank_genes_groups(data, **kwargs)[source]

rank_genes_groups_df(data, group, **kwargs)[source]

volcano(**kwargs)[source]

ma(**kwargs)[source]

top_genes_bar(**kwargs)[source]

qc_perturbation_counts(**kwargs)[source]

qc_summary(qc_result, **kwargs)[source]

materialize_rank_genes_groups(data, **kwargs)[source]

pca(data, **kwargs)[source]: Plot PCA scatter.

pca_variance_ratio(data, **kwargs)[source]: Plot PCA variance ratio.

pca_loadings(data, **kwargs)[source]: Plot PCA loadings.

umap(data, **kwargs)[source]: Plot UMAP embedding.

overlap_heatmap(result, **kwargs)[source]: Plot pairwise overlap heatmap. See crispyx.plot_overlap_heatmap().

Module Reference

The sections below document the underlying modules used by the namespace API. These are useful for advanced usage or for understanding parameter details.

Data loading and utilities

Helpers for working with AnnData .h5ad files in a streaming friendly way.

crispyx.data.drop_file_cache(path)[source]

Advise the kernel to drop page cache for path (Linux only).

On cgroup-limited systems (SLURM), page-cache pages count toward the memory limit. Calling this after a streaming read prevents the cached file data from consuming the cgroup budget.

The call is a no-op on non-Linux platforms or when the file cannot be opened.

Parameters:: path (str | Path)
Return type:: None

crispyx.data.is_dense_storage(path)[source]

Check if h5ad file stores X matrix as dense array.

Parameters:: path (str | Path) – Path to h5ad file.
Returns:: True if X is stored as dense array, False if sparse (CSR/CSC).
Return type:: bool

crispyx.data.get_matrix_storage_format(path)[source]

Detect matrix storage format in h5ad file.

Parameters:: path (str | Path) – Path to h5ad file.
Returns:: One of: ‘csr’, ‘csc’, ‘dense’
Return type:: str

class crispyx.data.AnnData(path, *, mode='r')[source]

Bases: object

Thin wrapper around a backed anndata.AnnData handle.

Parameters:

path (str | Path)
mode (str)

property filename: str: Return the underlying filename for compatibility with Scanpy.

property backed: AnnData: Return the lazily opened backed AnnData handle.

close()[source]

Close the underlying file handle if it is open.

Return type:: None

to_memory()[source]

Materialise the backed AnnData into memory.

Return type:: AnnData

property obs: _LazyFrameAccessor: Lazy accessor for observation (cell) metadata.

property var: _LazyFrameAccessor: Lazy accessor for variable (gene) metadata.

property uns: _LazyUnsMapping: Lazy accessor for unstructured annotations.

crispyx.data.read_backed(path, mode='r')[source]

Open an .h5ad file in backed mode, its matrices left on disk.

X and every layer stay on disk and are read only as they are sliced; obs, var, uns, obsm, varm, obsp, varp and .raw’s var are read into memory. This is anndata.read_h5ad(path, backed=mode) except for the layers, which anndata reads into memory in full – each as large as X.

Close the file with adata.file.close(); the layers, like X, cannot be read after that. adata.to_memory() loads everything.

Parameters:

path (str | Path)
mode (Literal['r', 'r+'])

Return type:

AnnData

crispyx.data.write_obsm_to_h5ad(path, key, data)[source]

Write a dense array to obsm/{key} in an h5ad file.

Parameters:

path (str | Path | File) – Path to the h5ad file, or an already-open h5py.File handle.
key (str) – Key under obsm (e.g., ‘X_pca’).
data (ndarray) – Dense numpy array of shape (n_obs, n_dims).

Return type:

None

crispyx.data.write_varm_to_h5ad(path, key, data)[source]

Write a dense array to varm/{key} in an h5ad file.

Parameters:

path (str | Path | File) – Path to the h5ad file, or an already-open h5py.File handle.
key (str) – Key under varm (e.g., ‘PCs’).
data (ndarray) – Dense numpy array of shape (n_vars, n_dims).

Return type:

None

crispyx.data.write_uns_dict_to_h5ad(path, key, data)[source]

Write a dict to uns/{key} in an h5ad file.

Handles scalar values, numpy arrays, and nested dicts. Uses AnnData-compatible encoding for proper round-trip compatibility.

Parameters:

path (str | Path) – Path to the h5ad file.
key (str) – Key under uns (e.g., ‘pca’).
data (dict) – Dictionary with string keys and scalar/array values.

Return type:

None

crispyx.data.write_obsp_to_h5ad(path, key, data)[source]

Write a sparse matrix to obsp/{key} in an h5ad file.

Stores in CSR format following AnnData conventions.

Parameters:

path (str | Path | File) – Path to the h5ad file, or an already-open h5py.File handle.
key (str) – Key under obsp (e.g., ‘distances’, ‘connectivities’).
data (spmatrix) – Sparse matrix of shape (n_obs, n_obs).

Return type:

None

crispyx.data.write_varp_to_h5ad(path, key, data)[source]

Write a sparse matrix to varp/{key} in an h5ad file.

Stores in CSR format following AnnData conventions.

Parameters:

path (str | Path | File) – Path to the h5ad file, or an already-open h5py.File handle.
key (str) – Key under varp (e.g., a gene-gene network/correlation matrix).
data (spmatrix) – Sparse matrix of shape (n_vars, n_vars).

Return type:

None

crispyx.data.resolve_data_path(data, *, require_exists=True)[source]

Resolve the on-disk path for a backed AnnData object or path-like input.

This utility supports flexible input types for crispyx functions, allowing users to pass either a file path or an AnnData object.

Parameters:

data (str | Path | 'AnnData' | ad.AnnData) – One of: - A string or Path to an h5ad file - A crispyx.AnnData wrapper (has .path attribute) - A backed anndata.AnnData object (has .filename attribute)
require_exists (bool) – If True (default), verify the resolved path exists.

Returns:

The resolved file path.

Return type:

Path

Raises:

TypeError – If data is an in-memory (non-backed) AnnData or unsupported type.
FileNotFoundError – If require_exists is True and the path does not exist.

Examples

>>> from crispyx.data import resolve_data_path
>>> path = resolve_data_path("data/counts.h5ad")
>>> path = resolve_data_path(adata_wrapper)
>>> path = resolve_data_path(backed_adata)

crispyx.data.resolve_output_path(input_path, *, suffix, data_name=None, output_path=None, output_dir=None)[source]

Construct an informative output path for an intermediate .h5ad file.

Output names follow the pattern {stem}_cx_{suffix}.h5ad where stem is the input file’s stem (or data_name when provided). The _cx_ token sits between the input identity and the processing step so that:

filenames sort naturally by input first, then by operation;
the producing tool is unambiguously marked without cluttering the start of the filename.

When output_path is provided it is returned directly, bypassing all name computation. This lets callers enforce an exact output filename. Either way, a path that is the input file itself raises ValueError.

Deprecated since version ``output_dir``: is kept for backward compatibility and will be removed in the next major version. Use output_path instead.

Parameters:

input_path (str | Path)
suffix (str)
data_name (str | None)
output_path (str | Path | None)
output_dir (str | Path | None)

Return type:

Path

crispyx.data.ensure_gene_symbol_column(adata, gene_name_column)[source]

Return a vector of gene symbols and verify they look like symbols, not Ensembl IDs.

Parameters:

adata (ad.AnnData | ad._core.anndata.AnnDataMixin)
gene_name_column (str | None)

Return type:

pd.Index

crispyx.data.resolve_control_label(labels, control_label, *, verbose=True)[source]

Return an explicit control label, inferring one when necessary.

Parameters:

labels (Sequence[str])
control_label (str | None)
verbose (int | bool)

Return type:

str

crispyx.data.read_h5ad_ondisk(path, *, n_obs=5, n_vars=5)[source]

Open an .h5ad file on disk, print a preview, and return a read-only view.

Parameters:

path (str | Path)
n_obs (int)
n_vars (int)

Return type:

AnnData

crispyx.data.validate_format_mismatch_policy(policy)[source]

Raise ValueError unless policy is one of the four accepted values.

Every public entry point that takes format_mismatch_policy calls this first, so a typo fails immediately – before a cached result is returned or a disk estimate silently drops its "scratch" entry.

Parameters:: policy (str)
Return type:: None

crispyx.data.resolve_auto_format_mismatch_policy(path, *, axis, chunk_size)[source]

Resolve format_mismatch_policy="auto" for a mismatched source.

Converts only when the projected cost of not converting – one full re-read of X per chunk beyond the first – is both structurally amortizable (_SLOW_AXIS_MIN_CHUNKS chunks) and large in absolute terms (_SLOW_AXIS_MIN_SAVING_SECONDS). Declining resolves to "warn" when the projected cost is material anyway (so the user hears about it and can convert the file once, up front) and to "off" when it is negligible, which keeps small local runs silent.

Callers must have established that path’s storage format really is mismatched with axis; estimate_disk_usage() calls this too, so its "scratch" entry reflects what the real call will do. Within one process that holds exactly, because the throughput measurement is cached per file; a query and a run in separate processes can still resolve differently if the filesystem’s speed changed between them.

Parameters:

path (Path)
axis (int)
chunk_size (int)

Return type:

_AutoPolicy

crispyx.data.scratch_copy_bytes(path, *, axis, policy, chunk_size)[source]

Disk the temporary fast-axis copy will need, or None for no copy.

The one place the disk-usage resolvers in de.py, batch.py and this module ask “will this call convert?”, so an estimate agrees with what the run then does – including under "auto", where the answer depends on measured throughput rather than on the policy alone, and so requires the same bounded probe read of the source that the run makes (cached, so the two see one measurement and the estimate is not charged twice).

Parameters:

path (str | Path)
axis (int)
policy (str)
chunk_size (int)

Return type:

float | None

crispyx.data.stream_on_fast_axis(path, *, axis, policy, fn_name, scratch_dir, chunk_size, memory_limit_gb=None, verbose=False, resumable=False)[source]

Yield the path to stream axis chunks of X from.

A CSR matrix streamed by columns (axis=1), or a CSC one streamed by rows (axis=0), makes anndata read the whole data/indices arrays for every chunk – one full read of the file per chunk instead of one in total, plus the matrix’s full size in transient memory. policy (the public format_mismatch_policy) decides what to do about that:

"auto": measure what the re-reads would actually cost on this filesystem and convert only when that is worth an extra read and write (see resolve_auto_format_mismatch_policy()), then behave as the resolved policy below. This is the default everywhere.
"convert": write a temporary fast-axis copy under scratch_dir (convert_to_csc() / convert_to_csr(), output buffers bounded by memory_limit_gb), yield its path, and delete it on exit. When the free space under scratch_dir cannot hold that copy, the call falls back to "warn" (one UserWarning naming the shortfall) rather than failing with ENOSPC after minutes of I/O.
"warn": yield path after one UserWarning quantifying the cost (and the one-time crispyx.data logger line).
"off": yield path silently.

A source already stored on its fast axis, or dense, is yielded unchanged. The yielded path holds X only (plus obs/var): read uns, layers, obsm etc. from the original path. Callers have made the format decision here, so they stream the result with iter_matrix_chunks(..., warn_slow_axis=False).

resumable says the caller was invoked with resume=True. The copy made here lives only for the duration of one call, so a computation that takes more restarts than one – exactly what resume exists for – rebuilds it from scratch every time. That is worth one warning: converting the source once, up front, is strictly cheaper for such a run.

Parameters:

path (Path)
axis (int)
policy (str)
fn_name (str)
scratch_dir (Path)
chunk_size (int)
memory_limit_gb (float | None)
verbose (int | bool)
resumable (bool)

Return type:

Iterator[Path]

crispyx.data.iter_matrix_chunks(adata, *, axis=0, chunk_size=1024, convert_to_dense=True, matrix=None, start_chunk=0, warn_slow_axis=True, selected=None, read_ahead=False)[source]

Yield chunks of matrix (defaults to adata.X).

Pass e.g. matrix=adata.layers[key] to stream a layer with the same chunking and slow-axis-format dispatch used for X. start_chunk skips directly to that 0-indexed chunk (e.g. for resuming a previously interrupted run) without reading the skipped chunks from disk. warn_slow_axis=False suppresses the one-time logger warning about streaming a compressed matrix off its fast axis; callers that already resolved that decision through stream_on_fast_axis() pass it. selected, a boolean mask along axis, skips (without reading) the chunks that hold no selected entry; the chunks yielded are whole. read_ahead=True reads the next chunk on a background thread while the caller works on the current one – worth it when that work releases the GIL (a numba kernel), so reading and computing overlap.

Parameters:

adata (ad.AnnData | ad._core.anndata.AnnDataMixin)
axis (int)
chunk_size (int)
convert_to_dense (bool)
matrix (np.ndarray | sp.spmatrix | None)
start_chunk (int)
warn_slow_axis (bool)
selected (np.ndarray | None)
read_ahead (bool)

Return type:

Iterator[tuple[slice, np.ndarray | sp.spmatrix]]

crispyx.data.chunk_starts(length, chunk_size, *, start_chunk=0, selected=None)[source]

First index of each chunk iter_matrix_chunks() reads: every chunk_size-th from chunk start_chunk on, less those holding no selected entry.

Parameters:

length (int)
chunk_size (int)
start_chunk (int)
selected (ndarray | None)

Return type:

list[int]

crispyx.data.normalize_total_block(block, *, library_size=None, target_sum=10000.0, dtype=<class 'numpy.float64'>)[source]

Return a library-size normalised dense view of block.

Parameters:

block (ndarray | spmatrix) – A slice of the expression matrix with shape (n_cells, n_genes).
library_size (ndarray | None) – Optional precomputed library sizes for the cells in block. When None the library size is computed from block directly.
target_sum (float) – Target total counts per cell after normalisation, matching the default used by scanpy.pp.normalize_total().
dtype (dtype | type) – Data type for the output dense array. Defaults to float64.

Returns:

A tuple (normalised, library_size) where normalised is a dense array containing the normalised counts and library_size contains the per-cell library sizes that were used.

Return type:

tuple

crispyx.data.write_filtered_subset(source_path, *, cell_mask, gene_mask, output_path, chunk_size=4096, var_assignments=None, row_nnz=None, total_nnz=None, data_dtype=None, chunk_cache=None, format_mismatch_policy='auto')[source]

Stream a filtered AnnData view to disk without materialising X.

Every slot is carried through the filter, not just X/obs/var: layers are streamed the same way as X (chunked, never fully materialised); obsm/varm/obsp/varp/uns are copied in memory (they are small relative to X — dense embeddings, sparse cell-cell/gene-gene graphs, or global metadata) and subset on whichever axis applies (obsm/obsp by cell_mask, varm/varp by gene_mask). A source .raw is not copied — a warning is emitted instead, since regenerating it from the source is the caller’s responsibility.

Parameters:

source_path (str | Path) – Path to source h5ad file.
cell_mask (ndarray) – Boolean mask for cells to include.
gene_mask (ndarray) – Boolean mask for genes to include.
output_path (str | Path) – Path for output h5ad file.
chunk_size (int) – Number of cells to process per chunk.
var_assignments (dict[str, Sequence] | None) – Optional dict of column assignments for var DataFrame.
row_nnz (ndarray | None) – Optional pre-computed non-zero counts per row for X. When provided along with total_nnz and data_dtype, skips X’s counting pass (layers always run their own counting pass).
total_nnz (int | None) – Optional pre-computed total non-zero count for X.
data_dtype (dtype | None) – Optional pre-computed data type for X’s sparse matrix.
chunk_cache (Any | None) – Optional _ChunkCache object from qc module. When provided, reads X’s CSR data from cache instead of re-reading the source matrix.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) – What to do when X is CSC, which streaming by rows re-reads in full for every chunk; see stream_on_fast_axis(). A temporary CSR copy goes beside output_path. Layers are streamed as stored.

Return type:

None

crispyx.data.concat_h5ad(paths, *, output_path, shared_cells='drop', label=None, chunk_size=4096, format_mismatch_policy='auto', verbose=True, force=False)[source]

Concatenate h5ad files along cells, on disk, dropping cells they share.

Screens are often delivered as several files that each carry the same control cells; concatenating them as they are counts every shared cell once per file. Cells are matched by obs_names: a barcode that occurs in more than one file must be the same cell there – the same number of non-zeros and the same total over the shared genes – or the call raises. obs values that differ between the copies (e.g. QC metrics computed per file) are warned about, and the first file’s are kept.

Each file’s matrix is read twice, row chunk by row chunk: once to size the output and check the shared cells, once to write it into one CSR matrix. It is never held in memory.

Parameters:

paths (Sequence[str | Path]) – The h5ad files, in output order.
output_path (str | Path) – Where to write the result (not one of paths).
shared_cells (Literal['drop', 'error']) – "drop" (default) keeps each shared cell where it first occurs and drops its later copies, after checking they are the same cell. "error" raises if any cell occurs in more than one file.
label (str | None) – If given, the name of an obs column recording each cell’s source file, in the order of paths: the file’s stem, or the path as given if two files share a stem.
chunk_size (int) – Number of cells read per streaming block.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) – What to do with a file whose X is CSC, which streaming by rows re-reads in full for every chunk; see stream_on_fast_axis(). Temporary CSR copies go beside output_path. Layers are streamed as stored.
verbose (int | bool) – Level 1 reports the cells and genes dropped and what was not carried over.
force (bool) – Recompute even when output_path holds a result of this call: same files (path, size, modification time), in the same order, with the same arguments, as recorded in its uns["crispyx"].

Returns:

The concatenated file. var holds the genes every file has, in the first file’s order, with the first file’s var columns; genes missing from any file are dropped. X (always CSR) and every layer present in all files are concatenated; obs is the union of the files’ columns (missing values are NA; a shared cell takes the values its first file lacks from a later copy); obsm keys present in all files as the same kind (array, sparse or DataFrame) and width are concatenated; uns keeps the keys every file has with equal values. obsp (cell-cell graphs), varm/varp and .raw describe a single file and are not carried over.

Return type:

AnnData

Examples

>>> cx.pp.concat_h5ad(
...     ["pilot_a.h5ad", "pilot_b.h5ad"], output_path="pilot.h5ad", label="source",
... )

crispyx.data.normalize_total_log1p(data, output_path=None, *, normalize=True, log1p=True, target_sum=10000.0, chunk_size=4096, output_dir=None, data_name=None, format_mismatch_policy='auto', verbose=True)[source]

Stream normalize and/or log-transform an h5ad file without loading it fully into memory.

This function processes the source file in chunks, optionally applying: 1. Total-count normalization (scanpy.pp.normalize_total equivalent) 2. Log1p transformation (scanpy.pp.log1p equivalent)

The output is written as a sparse CSR matrix. This is the streaming equivalent of calling scanpy.pp.normalize_total followed by scanpy.pp.log1p.

Only X is transformed. layers/obsm/varm/obsp/varp/ uns are carried through unchanged (cell/gene identity, order, and count never change here), same as write_filtered_subset() and downsample_counts(); a source .raw is dropped with a warning.

Parameters:

data (str | Path | 'AnnData' | ad.AnnData) – Path to source h5ad file, or a backed AnnData object.
output_path (str | Path | None) – Path for output h5ad file. If None, uses output_dir/data_name pattern.
normalize (bool) – Whether to apply total-count normalization. Default True.
log1p (bool) – Whether to apply log1p transformation. Default True.
target_sum (float) – Target total counts per cell after normalization. Default 1e4. Only used if normalize=True.
chunk_size (int) – Number of cells to process per chunk. Default 4096.
output_dir (str | Path | None) – Directory for output file. Defaults to input file’s directory.
data_name (str | None) – Custom name for output file. If None, uses “normalized” suffix.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) –
How to handle a source stored as CSC, whose row-(cell-)streaming here is O(total_nnz) per chunk and can be ~100x slower than CSR:
- "auto" (default): measure the source’s read throughput and convert only when the re-reads would actually cost more than one conversion does – roughly, a large file on a slow (networked) filesystem streamed in many chunks. See resolve_auto_format_mismatch_policy().
- "warn": proceed but log a single actionable warning.
- "convert": always convert the source to CSR in a temporary file (bounded-memory two-pass streaming) and stream from that; the temporary file is removed before returning. This calls convert_to_csr() internally, which temporarily needs ~2x the source file’s size in free disk space and warns automatically if that looks tight.
- "off": proceed silently with no warning.
verbose (int | bool) – Print progress information.

Returns:

Read-only AnnData wrapper pointing to the output file.

Return type:

AnnData

Examples

>>> # Full normalization + log1p (default)
>>> adata_norm = cx.pp.normalize_total_log1p(adata_ro, output_dir=OUTPUT_DIR, data_name="normalized")

>>> # Only log1p (no normalization)
>>> adata_log = cx.pp.normalize_total_log1p(adata_ro, normalize=False, output_dir=OUTPUT_DIR)

>>> # Only normalization (no log1p)
>>> adata_norm = cx.pp.normalize_total_log1p(adata_ro, log1p=False, output_dir=OUTPUT_DIR)

>>> # Use explicit output path
>>> adata_norm = cx.pp.normalize_total_log1p(adata_ro, "results/normalized.h5ad")

>>> # Auto-convert a CSC source to CSR for fast cell-streaming
>>> adata_norm = cx.pp.normalize_total_log1p(csc_path, format_mismatch_policy="convert")

crispyx.data.downsample_counts(data, output_path=None, *, counts_per_cell, chunk_size=4096, data_name=None, random_state=0, verbose=True)[source]

Stream-thin every cell’s counts down to at most counts_per_cell.

This is the streaming, dependency-free equivalent of scanpy.pp.downsample_counts(adata, counts_per_cell=..., replace=False): a cell already at or below the target is left unchanged; a cell above it has its counts thinned via exact sampling without replacement (see _thin_csr_chunk()), never loading the full matrix into memory.

Unlike crispyx.pp.subsample(), this never changes which cells or genes are present – it thins the values within every surviving cell. Use subsample first if you also want to reduce cell counts; the two are independent, composable streaming passes.

Only X is thinned, matching scanpy.pp.downsample_counts. layers/obsm/varm/obsp/varp/uns are carried through unchanged (cell/gene identity, order, and count never change here, only the values in X); a source .raw is dropped with a warning, same as write_filtered_subset(). X must hold non-negative integer counts.

Parameters:

data (str | Path | 'AnnData' | ad.AnnData) – Path to source h5ad file, or a backed AnnData object.
output_path (str | Path | None) – Path for output h5ad file. If None, derived from data_name.
counts_per_cell (int) – Target total count per cell. Must be a positive integer.
chunk_size (int) – Number of cells to process per chunk. Default 4096.
data_name (str | None) – Custom name for output file.
random_state (int) – Seed. Thinning is deterministic for a fixed seed and independent of chunk_size.
verbose (int | bool) – Print progress information.

Returns:

Read-only AnnData wrapper pointing to the output file.

Return type:

AnnData

crispyx.data.convert_to_csc(data, *, output_path=None, chunk_size=None, memory_limit_gb=None, output_dir=None, data_name=None, verbose=True)[source]

Convert a backed h5ad file’s matrix from CSR (or dense) to CSC format.

CSC format allows O(nnz_in_chunk) column-slicing instead of O(total_nnz) that CSR requires. This is required for efficient Wilcoxon rank-sum testing, which iterates over gene chunks with axis=1 access patterns.

The conversion streams the source by rows: one counting pass, then the non-zeros are scattered into CSC output buffers one column band at a time and each finished band is written into the pre-sized output file. A band is as many contiguous columns as fit half of memory_limit_gb, so peak memory is bounded by that budget (plus one row-chunk working buffer) regardless of the matrix size. When the whole output fits – the common case – there is a single band and hence a single second pass; a matrix K times larger than the budget costs K passes over the source instead of running out of memory.

Peak disk usage is not similarly bounded: the source and destination files coexist until the caller deletes the source, so this temporarily requires roughly 2x the source file’s size in free disk space (e.g. ~54 GB for a 27 GB screen). A warning is emitted automatically if the output filesystem looks too tight for that.

If the input file is already CSC, no file is written; the function returns a backed AnnData pointing to the original file.

Parameters:

data (str | Path | 'AnnData' | ad.AnnData) – Path to source h5ad file (CSR or dense), or a backed AnnData.
output_path (str | Path | None) – Explicit path for the output file. If None, a path is derived from output_dir/data_name with "_csc" appended to the stem.
chunk_size (int | None) – Number of rows (cells) to read at a time during both passes. Default 4096.
memory_limit_gb (float | None) – Memory budget in GB; the output buffers use at most half of it, and the matrix is converted in several column bands when it does not fit. When None (default), detected via psutil. On HPC nodes, pass the SLURM --mem value.
output_dir (str | Path | None) – Directory for the output file. Defaults to the source file’s directory.
data_name (str | None) – Custom name used when building the output filename.
verbose (int | bool) – Print progress messages.

Returns:

Backed (read-only) AnnData pointing to the written CSC h5ad file, or to the source file if it was already CSC.

Return type:

AnnData

Examples

>>> adata_csc = cx.pp.convert_to_csc(preprocessed_path, output_dir=OUTPUT_DIR)

crispyx.data.convert_to_csr(data, *, output_path=None, chunk_size=None, memory_limit_gb=None, output_dir=None, data_name=None, verbose=True)[source]

Convert a backed h5ad file’s matrix from CSC (or dense) to CSR format.

CSR format allows O(nnz_in_chunk) row-slicing instead of O(total_nnz) that CSC requires. This is needed for efficient NB-GLM, size factor computation, and any operation that iterates over cell (row) chunks.

The conversion mirrors convert_to_csc(): a counting pass, then the non-zeros are scattered into CSR output buffers one row band at a time, each band sized to fit half of memory_limit_gb, so peak memory is bounded by that budget plus one chunk working buffer. A dense source is already row-ordered and is copied chunk by chunk with no whole-matrix buffer at all.

Peak disk usage is not similarly bounded: the source and destination files coexist until the caller deletes the source, so this temporarily requires roughly 2x the source file’s size in free disk space. A warning is emitted automatically if the output filesystem looks too tight for that. This also applies to normalize_total_log1p()’s format_mismatch_policy="convert", which calls this function internally.

If the input file is already CSR, no file is written; the function returns a backed AnnData pointing to the original file.

Parameters:

data (str | Path | 'AnnData' | ad.AnnData) – Path to source h5ad file (CSC or dense), or a backed AnnData.
output_path (str | Path | None) – Explicit path for the output file. If None, a path is derived from output_dir/data_name with "_csr" appended to the stem.
chunk_size (int | None) – Number of rows (cells) to read at a time during both passes. Default is calculated automatically.
memory_limit_gb (float | None) – Memory budget in GB; the output buffers use at most half of it, and a CSC source is converted in several row bands when it does not fit. When None (default), detected via psutil.
output_dir (str | Path | None) – Directory for the output file. Defaults to the source file’s directory.
data_name (str | None) – Custom name used when building the output filename.
verbose (int | bool) – Print progress messages.

Returns:

Backed (read-only) AnnData pointing to the written CSR h5ad file, or to the source file if it was already CSR.

Return type:

AnnData

crispyx.data.compress_h5ad(src, dst, *, level=4, verify=True, n_jobs=-1, overwrite=False, verbose=False)[source]

Losslessly re-encode an .h5ad (or any HDF5) file with gzip + shuffle.

Every dataset of at least 1024 fixed-width elements – X in any layout, layers, obsm/varm/obsp/varp, dataframe columns, uns arrays – is rewritten with HDF5’s built-in shuffle and deflate filters in ~1 MiB chunks. Nothing else changes: the same groups, datasets, dtypes, shapes, sparse layout, links and attributes (values, types and order). The walk is independent of AnnData, so it applies to any .h5ad layout, including entries AnnData itself would not round-trip. Small datasets and variable-length strings are copied as they are.

Any HDF5 reader decodes the result without plugins: anndata/scanpy (in memory and backed="r"), crispyx, R (rhdf5, zellkonverter, anndata-R), h5ls. Chunks are compressed on n_jobs threads, so a single file uses every core; memory stays at a few MiB per thread regardless of the file size. Rewriting the file also reclaims space freed by in-place edits (e.g. rerunning pp.pca on the same file), which HDF5 never returns on its own.

Typical ratios: ~0.15-0.3 for a single-cell matrix and 0.55-0.75 for dense effect or DE matrices. Reading a compressed file costs decompression time. crispyx’s streaming readers inflate its chunks on a thread pool, several times faster than h5py’s single-threaded decoding, but still slower than an uncompressed read from the page cache. On fast local disks, compress files for archiving or sharing rather than inputs you are about to stream repeatedly; on HDD or network storage the smaller file often reads faster.

Parameters:

src (str | Path) – Source .h5ad file. It is only read.
dst (str | Path) – Destination path. Passing dst=src (with overwrite=True) replaces the source in place, but only after the new file is complete and verified; the source is never left half-written.
level (int) – Deflate level, 1-9.
verify (bool) – Re-read both files through h5py and compare every group, attribute and dataset bit for bit before the destination is put in place. Costs roughly one extra read of both files.
n_jobs (int) – Compression threads; -1 uses every core.
overwrite (bool) – Allow replacing an existing dst.
verbose (int | bool) – Show progress.

Returns:

src_bytes, dst_bytes, ratio (dst_bytes / src_bytes) and verified.

Return type:

Examples

>>> cx.compress_h5ad(OUTPUT_DIR / "screen.h5ad", OUTPUT_DIR / "screen.gz.h5ad")
>>> cx.compress_h5ad(path, path, overwrite=True)  # in place, after verification

crispyx.data.calculate_optimal_chunk_size(n_obs, n_vars, available_memory_gb=None, safety_factor=8.0, min_chunk=512, max_chunk=4096)[source]

Calculate optimal chunk size based on dataset dimensions and available memory.

Parameters:

n_obs (int) – Number of observations (cells) in the dataset.
n_vars (int) – Number of variables (genes) in the dataset.
available_memory_gb (float | None) – Available memory in gigabytes. If None, auto-detects using psutil.
safety_factor (float) – Safety multiplier to account for overhead (default 8.0 for backed operations).
min_chunk (int) – Minimum chunk size to return (default 512).
max_chunk (int) – Maximum chunk size to return (default 4096).

Returns:

Recommended chunk size, clamped to [min_chunk, max_chunk].

Return type:

Examples

>>> calculate_optimal_chunk_size(100000, 20000, available_memory_gb=32)
2000

crispyx.data.calculate_optimal_gene_chunk_size(n_obs, n_vars, n_groups=None, available_memory_gb=None, safety_factor=8.0, memory_fraction=0.5, min_chunk=32, max_chunk=512)[source]

Calculate optimal gene chunk size for column-wise operations.

For operations that iterate over genes (columns), such as Wilcoxon tests, each chunk loads all cells for a subset of genes. Memory usage is dominated by n_obs × chunk_size rather than chunk_size × n_vars.

Enhanced to account for the number of perturbation groups, which significantly impacts memory usage due to output array allocation.

Parameters:

n_obs (int) – Number of observations (cells) in the dataset.
n_vars (int) – Number of variables (genes) in the dataset.
n_groups (int | None) – Number of perturbation groups. If provided, used to estimate memory for output arrays. Large group counts require smaller chunks.
available_memory_gb (float | None) – Available memory in gigabytes. If None, auto-detects using psutil.
safety_factor (float) – Safety multiplier to account for overhead (default 8.0).
memory_fraction (float) – Fraction of available memory to use (default 0.5). Leave headroom for memory-mapped arrays and system overhead.
min_chunk (int) – Minimum chunk size to return (default 32).
max_chunk (int) – Maximum chunk size to return (default 512).

Returns:

Recommended gene chunk size, clamped to [min_chunk, max_chunk].

Return type:

Examples

>>> calculate_optimal_gene_chunk_size(100000, 20000, available_memory_gb=32)
512
>>> calculate_optimal_gene_chunk_size(4000000, 38000, n_groups=18000, available_memory_gb=128)
64

crispyx.data.calculate_wilcoxon_chunk_size(n_obs, n_vars, *, available_memory_gb=None, min_chunk=32, max_chunk=4096)[source]

Calculate optimal gene chunk size for Wilcoxon rank-sum tests.

Unlike calculate_optimal_gene_chunk_size(), this function has no n_groups cap. Wilcoxon writes all output arrays (effect, pvalue, z-score, etc.) to on-disk memory-mapped files immediately, so peak RAM per chunk is independent of the number of perturbation groups. The only effective cap is a per-chunk transient-memory budget:

transient ≈ chunk_size × (n_obs × 4 + n_ctrl × 8 + n_pert × 4)
≈ chunk_size × n_obs × 12 bytes

The budget is set to 15 % of available_memory_gb so that a single chunk never exceeds ~1/7th of the node RAM.

Parameters:

n_obs (int) – Number of cells in the dataset.
n_vars (int) – Number of genes in the dataset (used only for logging).
available_memory_gb (float | None) – Available memory in GB. When None, detected via psutil. On HPC nodes, pass the SLURM --mem value so the cap reflects the actual job allocation rather than system-wide free memory.
min_chunk (int) – Floor for the returned chunk size (default 32).
max_chunk (int) – Ceiling for the returned chunk size (default 4096). The cell-budget cap is applied before this ceiling, so max_chunk is only active for small/sparse datasets where the cell cap would be very large.

Returns:

Recommended gene chunk size, clamped to [min_chunk, max_chunk].

Return type:

Examples

>>> # Feng-gwsnf: 396K cells, 128 GB → 4067
>>> calculate_wilcoxon_chunk_size(396458, 32373, available_memory_gb=128)
4067
>>> # Feng-ts: 1.16M cells, 128 GB → 1378
>>> calculate_wilcoxon_chunk_size(1161864, 33165, available_memory_gb=128)
1378

crispyx.data.calculate_nb_glm_chunk_size(n_obs, n_vars, n_groups=None, available_memory_gb=None, memory_limit_gb=None, safety_factor=8.0, memory_fraction=0.5, min_chunk=32, max_chunk=256)[source]

Calculate optimal gene chunk size for NB-GLM operations.

NB-GLM iterates over genes (columns) and for each chunk: - Loads dense count data: n_obs × chunk_size - Computes dispersion: requires cell-level statistics - Fits GLM: design matrix operations

This function is specifically tuned for NB-GLM memory patterns, which differ from Wilcoxon (ranking) and t-test (simple statistics).

Parameters:

n_obs (int) – Number of observations (cells) in the dataset.
n_vars (int) – Number of variables (genes) in the dataset.
n_groups (int | None) – Number of perturbation groups. If provided, used to estimate memory for output arrays and design matrix overhead.
available_memory_gb (float | None) – Available memory in gigabytes. If None, auto-detects using psutil.
memory_limit_gb (float | None) – Optional hard memory limit in GB. If provided, uses the minimum of available memory and this limit.
safety_factor (float) – Safety multiplier to account for overhead (default 8.0).
memory_fraction (float) – Fraction of available memory to use (default 0.5).
min_chunk (int) – Minimum chunk size to return (default 32).
max_chunk (int) – Maximum chunk size to return (default 256).

Returns:

Recommended gene chunk size, clamped to [min_chunk, max_chunk]. For datasets where memory is sufficient, returns max_chunk (256). Only reduces chunk size when memory would be exceeded.

Return type:

Examples

>>> calculate_nb_glm_chunk_size(100000, 20000, n_groups=100, available_memory_gb=128)
256
>>> calculate_nb_glm_chunk_size(1200000, 36000, n_groups=500, available_memory_gb=128)
143

crispyx.data.calculate_pca_chunk_size(n_obs, n_vars, n_comps=50, available_memory_gb=None, method='auto', memory_fraction=0.5, min_chunk=256, max_chunk=4096)[source]

Calculate optimal chunk size for streaming PCA.

PCA memory usage depends on the method: - sparse_cov: O(genes²) for covariance matrix, fast for small gene counts - incremental: O(chunk × genes) for data chunks, better for large gene counts

Parameters:

n_obs (int) – Number of observations (cells) in the dataset.
n_vars (int) – Number of variables (genes) in the dataset.
n_comps (int) – Number of principal components to compute. Default 50.
available_memory_gb (float | None) – Available memory in gigabytes. If None, auto-detects using psutil.
method (str) – PCA method: ‘auto’, ‘sparse_cov’, or ‘incremental’. ‘auto’ selects based on gene count and available memory.
memory_fraction (float) – Fraction of available memory to use (default 0.5).
min_chunk (int) – Minimum chunk size to return (default 256).
max_chunk (int) – Maximum chunk size to return (default 4096).

Returns:

(chunk_size, selected_method) where selected_method is ‘sparse_cov’ or ‘incremental’.

Return type:

tuple[int, str]

Examples

>>> calculate_pca_chunk_size(100000, 8000, available_memory_gb=32)
(2048, 'sparse_cov')
>>> calculate_pca_chunk_size(100000, 50000, available_memory_gb=32)
(1024, 'incremental')

crispyx.data.calculate_adaptive_qc_thresholds(adata, perturbation_column, mode='conservative', sample_size=10000, chunk_size=None)[source]

Calculate adaptive QC thresholds based on data distribution.

Uses percentile-based approach to determine appropriate QC parameters that retain most of the data while filtering outliers.

Parameters:

adata (AnnData) – AnnData object (can be backed).
perturbation_column (str) – Column in adata.obs containing perturbation labels.
mode (str) – ‘conservative’ (10th percentile, retains ~90%) or ‘aggressive’ (5th percentile, retains ~95%).
sample_size (int) – Maximum number of cells to sample for gene expression analysis.
chunk_size (int | None) – Optional fixed chunk size to use. If None, calculated adaptively.

Returns:

Dictionary with keys: min_genes, min_cells_per_perturbation, min_cells_per_gene, chunk_size.

Return type:

Examples

>>> adata = ad.read_h5ad("data.h5ad", backed='r')
>>> thresholds = calculate_adaptive_qc_thresholds(adata, "perturbation")
>>> adata.file.close()

crispyx.data.standardize_dataset(dataset_path, perturbation_column, control_label, gene_name_column, output_dir, force=False)[source]

Standardize dataset column names and control labels with caching.

Creates a standardized copy of the dataset with: - perturbation_column renamed to ‘perturbation’ - control labels standardized to ‘control’ - gene_name_column set as var.index if specified

This function uses a streaming approach to avoid loading the X matrix into memory, making it suitable for very large datasets (>1M cells).

Standardized files are cached in {output_dir}/.cache/ and reused unless force=True.

Parameters:

dataset_path (Path | str) – Path to original dataset (.h5ad file).
perturbation_column (str) – Name of perturbation column in original dataset.
control_label (str | None) – Control label to standardize. If None, auto-detects.
gene_name_column (str | None) – Gene name column to use as var.index. If None, uses existing var.index.
output_dir (Path | str) – Directory for cached standardized files.
force (bool) – If True, regenerate standardized file even if cache exists.

Returns:

Path to standardized dataset (either cached or newly created).

Return type:

Path

Examples

>>> standardized_path = standardize_dataset(
...     "data/original.h5ad",
...     perturbation_column="gene",
...     control_label=None,
...     gene_name_column="gene_symbols",
...     output_dir="results",
...     force=False
... )

crispyx.data.needs_sorting_for_nbglm(path, perturbation_column='perturbation', *, min_cells=360000, min_perturbations=100, contiguity_threshold=0.5)[source]

Check if a dataset would benefit from sorting by perturbation for NB-GLM.

Large datasets with scattered cells benefit from having cells sorted by perturbation label, as this enables contiguous I/O reads instead of random access. This is especially important when the data is stored on HDD (rotational disk).

The default thresholds are based on I/O overhead analysis: - At ~100 IOPS (typical HDD), 360K cells = 1 hour of random I/O overhead - min_perturbations=100 ensures sufficient parallel workload to benefit - contiguity_threshold=0.5 catches scattered datasets

Parameters:

path (str | Path) – Path to h5ad file.
perturbation_column (str) – Column in obs containing perturbation labels.
min_cells (int) – Minimum number of cells for sorting to be recommended. Default: 360,000 (~1 hour of random I/O on HDD at 100 IOPS).
min_perturbations (int) – Minimum number of perturbations for sorting to be recommended.
contiguity_threshold (float) – If average contiguity is below this threshold, sorting is recommended. Contiguity is the fraction of a perturbation’s cells that are in a contiguous block (1.0 = perfectly contiguous, 0.0 = completely scattered).

Returns:

True if sorting is recommended, False otherwise.

Return type:

bool

crispyx.data.sort_by_perturbation(path, perturbation_column='perturbation', control_label=None, *, output_path=None, chunk_size=4096, force=False)[source]

Sort cells by perturbation label for contiguous I/O access.

Creates a new h5ad file with cells reordered so that all cells from each perturbation are contiguous. Control cells are placed first, followed by each perturbation group in alphabetical order. This enables efficient sequential reads when processing perturbations in parallel.

The function works in streaming mode to handle datasets larger than memory. Sorting information is stored in uns[‘sorting_metadata’].

Parameters:

path (str | Path) – Path to input h5ad file.
perturbation_column (str) – Column in obs containing perturbation labels.
control_label (str | None) – Label for control cells. If None, auto-detected from common patterns.
output_path (str | Path | None) – Path for output sorted file. If None, appends ‘_sorted’ to input name.
chunk_size (int) – Number of cells to process per chunk during streaming write.
force (bool) – If True, recreate sorted file even if it already exists.

Returns:

Path to the sorted h5ad file.

Return type:

Path

Examples

>>> sorted_path = sort_by_perturbation(
...     "data/large_dataset.h5ad",
...     perturbation_column="perturbation",
...     control_label="control",
... )
>>> # sorted_path is now "data/large_dataset_sorted.h5ad"

Notes

The sorted file contains additional metadata in uns[‘sorting_metadata’]: - original_path: Path to the original unsorted file - sort_order: Array mapping new indices to original indices - perturbation_boundaries: Dict mapping perturbation labels to (start, end) indices - sorted_by: The column used for sorting - timestamp: When the sorting was performed

For sparse inputs the output is always written in CSR format, since sorting benefits row-wise (per-perturbation) access patterns used by NB-GLM. CSC files used by Wilcoxon do not need perturbation sorting.

crispyx.data.get_perturbation_slice(adata_or_path, perturbation_label, perturbation_column='perturbation')[source]

Get slice for a perturbation’s cells, checking if file is sorted.

If the file is sorted by perturbation, returns a contiguous slice. Otherwise, returns None for the slice (caller should use boolean mask).

Parameters:

adata_or_path (str | Path | AnnData) – Path to h5ad file, or an already-opened AnnData object.
perturbation_label (str) – Label of the perturbation to get slice for.
perturbation_column (str) – Column containing perturbation labels.

Returns:

(slice object or None, is_sorted flag). If is_sorted is True, slice is valid for contiguous access. If is_sorted is False, slice is None and caller should use mask.

Return type:

tuple[slice | None, bool]

crispyx.data.load_obs(path)[source]

Load the obs metadata table from an h5ad file without reading X.

Parameters:: path (str | Path) – Path to the h5ad file.
Returns:: Full obs DataFrame in memory.
Return type:: pd.DataFrame

crispyx.data.load_var(path)[source]

Load the var metadata table from an h5ad file without reading X.

Parameters:: path (str | Path) – Path to the h5ad file.
Returns:: Full var DataFrame in memory.
Return type:: pd.DataFrame

crispyx.data.write_obs(path, df)[source]

Overwrite the obs metadata table in an h5ad file without touching X.

Parameters:

path (str | Path) – Path to the h5ad file (modified in-place).
df (DataFrame) – New obs DataFrame. Must have the same number of rows as the existing obs table. Index values are written as cell barcodes.

Raises:

ValueError – If the DataFrame length does not match the existing n_obs.

Return type:

None

crispyx.data.write_var(path, df)[source]

Overwrite the var metadata table in an h5ad file without touching X.

Parameters:

path (str | Path) – Path to the h5ad file (modified in-place).
df (DataFrame) – New var DataFrame. Must have the same number of rows as the existing var table.

Raises:

ValueError – If the DataFrame length does not match the existing n_vars.

Return type:

None

crispyx.data.standardise_gene_names(path, *, column=None, strip_version=True, normalise_mt_prefix=True, lookup_symbols=False, species='human', unmapped_action='warn', inplace=True)[source]

Standardise gene identifiers in the var metadata table.

Applies a deterministic normalisation pipeline:

Strip Ensembl version suffixes (ENSG00000123.4 → ENSG00000123).
Normalise mt- prefix to MT- (human mitochondrial convention).
Optionally resolve Ensembl IDs to HGNC symbols via mygene (requires pip install mygene). A tqdm progress bar is shown during batched lookups.

Parameters:

path (str | Path) – Path to the h5ad file.
column (str | None) – var column to normalise. None normalises the index (var_names).
strip_version (bool) – Strip ".N" Ensembl version suffixes.
normalise_mt_prefix (bool) – Convert lower-case mt- prefix to MT-.
lookup_symbols (bool) – If True, query mygene to map Ensembl IDs → gene symbols.
species (str) – Species string passed to mygene (default "human").
unmapped_action (Literal['keep', 'error', 'warn']) – What to do for IDs not found by mygene: "keep" leaves them unchanged, "warn" emits a warning, "error" raises.
inplace (bool) – If True, write the result back to the file and return None. If False, return a Series without modifying the file.

Returns:

Normalised gene names when inplace=False, else None.

Return type:

pd.Series or None

crispyx.data.normalise_perturbation_labels(path, column, *, strip_prefixes=None, strip_suffixes=None, strip_suffix_regex=None, control_aliases=None, canonical_control='NTC', inplace=True)[source]

Normalise perturbation labels stored in an obs column.

Applies transformations in order:

Strip specified prefixes via vectorised pd.Series.str.replace.
Strip specified suffixes via vectorised pd.Series.str.replace.
Apply a custom regex substitution (strip_suffix_regex).
Map known control aliases to canonical_control.

Parameters:

path (str | Path) – Path to the h5ad file.
column (str) – obs column containing perturbation labels.
strip_prefixes (list[str] | None) – List of prefix strings to remove (e.g. ["sg-", "sg"]).
strip_suffixes (list[str] | None) – List of suffix strings to remove (e.g. ["_KO", "_KD", "_P1P2"]).
strip_suffix_regex (str | None) – A Python regex applied via pd.Series.str.replace after prefix/suffix stripping.
control_aliases (list[str] | None) – Additional strings (case-insensitive) treated as control labels. The built-in aliases (ntc, ctrl, scramble, …) are always included.
canonical_control (str) – Canonical control label substituted for all matched aliases.
inplace (bool) – If True, write result back and return None. If False, return a Series without modifying the file.

Returns:

Normalised labels when inplace=False, else None.

Return type:

pd.Series or None

crispyx.data.detect_perturbation_column(adata, *, control_label=None, min_unique=2, verbose=True)[source]

Heuristically identify the obs column containing perturbation labels.

Scoring:

+3 if column name matches known aliases (perturbation, gene_target, …).
+2 if dtype is categorical or object.
+1 if unique-value count is in [min_unique, 5000].
+2 if at least one value matches a known control synonym or control_label when provided.

Parameters:

adata (str | Path | AnnData | AnnData) – Backed AnnData, AnnData, or path to h5ad file.
control_label (str | None) – Known control label; boosts the score for columns containing it.
min_unique (int) – Minimum number of unique values required for the +1 bonus.
verbose (int | bool) – Log the detected column name.

Returns:

Column name with the highest score, or None if no column scores above zero.

Return type:

str or None

crispyx.data.detect_gene_symbol_column(adata, *, verbose=True)[source]

Heuristically identify the var column containing gene symbols.

Scoring:

+3 if column name matches known aliases (gene_symbols, symbol, …).
+2 if values pass _validate_gene_symbols() without error.
+1 if values do not start with Ensembl prefixes.

Returns None when no column qualifies, which signals that var_names should be used as a fallback.

Parameters:

adata (str | Path | AnnData | AnnData) – Backed AnnData, AnnData, or path to h5ad file.
verbose (int | bool) – Log the detected column name.

Return type:

str or None

crispyx.data.infer_columns(adata, *, control_label=None, verbose=True)[source]

Detect perturbation and gene-symbol columns in a single call.

Parameters:

adata (str | Path | AnnData | AnnData) – Backed AnnData, AnnData, or path to h5ad file.
control_label (str | None) – Known control label forwarded to detect_perturbation_column().
verbose (int | bool) – Log detected column names.

Returns:

{"perturbation_column": ..., "gene_name_column": ...} where each value is the detected column name or None.

Return type:

class crispyx.data.OverlapResult(count_matrix, jaccard_matrix, set_sizes)[source]

Bases: object

Pairwise overlap statistics between named sets.

Parameters:

count_matrix (DataFrame)
jaccard_matrix (DataFrame)
set_sizes (Series)

count_matrix

(n_sets × n_sets) DataFrame of pairwise intersection sizes.

Type:: pandas.core.frame.DataFrame

jaccard_matrix

(n_sets × n_sets) DataFrame of Jaccard similarity coefficients.

Type:: pandas.core.frame.DataFrame

set_sizes

Series of sizes for each input set.

Type:: pandas.core.series.Series

count_matrix: DataFrame

jaccard_matrix: DataFrame

set_sizes: Series

crispyx.data.compute_overlap(sets_dict, *, metric='both')[source]

Compute pairwise overlap statistics between named sets.

Parameters:

sets_dict (dict[str, set | list]) – Mapping of name → set (or list, converted to set).
metric (Literal['count', 'jaccard', 'both']) – Which matrices to populate: "count", "jaccard", or "both" (default).

Returns:

Object with count_matrix, jaccard_matrix, and set_sizes.

Return type:

OverlapResult

Examples

>>> result = cx.tl.compute_overlap({
...     "dataset_A": {"BRCA1", "TP53", "EGFR"},
...     "dataset_B": {"TP53", "KRAS"},
... })
>>> result.jaccard_matrix

Disk-usage estimation

crispyx automatically warns – without blocking the call – when a streaming write’s disk footprint looks large relative to free space. To check usage before committing to a run, call estimate_disk_usage with the function you intend to run, its input file, and the same keyword arguments you plan to pass. Also available as cx.tl.estimate_disk_usage (see the Tools namespace above). See Disk space for the full explanation and a worked example.

crispyx.estimate_disk_usage(func, data, **kwargs)[source]

Estimate the disk space a crispyx function will need, before running it.

This is a standalone, on-demand query – separate from the automatic warning crispyx already emits from inside a real call when its own estimate looks tight. Resolvers read only cheap obs/uns metadata in backed mode to reproduce the group/batch counts the real function would compute in its own preamble; they never read X in full. The one exception is format_mismatch_policy="auto" (the default) on a source stored off the axis the function streams: deciding whether that call would write a temporary fast-axis copy means measuring the filesystem, so the query makes the same bounded 64 MB probe read of X the real call makes. The measurement is cached per file, so a call made later in the same process does not pay for it again and cannot resolve the decision differently than this query did.

Also available as cx.tl.estimate_disk_usage for Scanpy-style namespace discovery (the same pattern as compute_overlap, which is both top-level and under cx.tl); the two are the same function and always agree.

Parameters:

func (str | Callable) – Either the function’s name (e.g. "compute_normalized_effects", "t_test", "convert_to_csc") or the function object itself (e.g. crispyx.t_test). Both are accepted: the string form is the primary interface, and accepting the callable too costs one getattr while removing a class of typo’d-string bugs.
data (str | Path | AnnData) – Path to an h5ad file, or a backed/wrapped AnnData – same meaning as every other crispyx function’s data parameter.
**kwargs – The subset of the target function’s keyword arguments that affect group/shape counts (e.g. perturbation_column, control_label, batch_column, groupby, perturbations, min_cells). format_mismatch_policy matters too, as do chunk_size and memory_limit_gb through it: those two set how many chunks the run streams, which is what decides whether it converts the source and so whether "scratch" appears at all. Pass the ones you will pass to the real call; omitted, they resolve the way the real call resolves them. Keywords that genuinely cannot change the estimate (verbose, n_jobs, …) are accepted and ignored.

Returns:

Keyed by filesystem location – "tempdir" for disk-backed intermediate accumulators in $TMPDIR, "output" for the final result file, and "scratch" for the temporary fast-axis copy a format_mismatch_policy="convert" call writes beside the output. Not every function uses all of them; a conversion function like convert_to_csc only has "output". "output" and "scratch" are assessed where the target function would write: pass the same output_path (or output_dir) you will pass to it; with neither, the source file’s directory.

Return type:

dict[str, DiskEstimate]

Examples

>>> import crispyx as cx
>>> cx.estimate_disk_usage(
...     "compute_normalized_effects", "screen.h5ad",
...     perturbation_column="guide_target", batch_column="gem_group",
... )
{'tempdir': ..., 'output': ...}

Quality control

Quality control utilities for large .h5ad datasets.

class crispyx.qc.QualityControlResult(cell_mask, gene_mask, perturbation_keep, filtered, cell_gene_counts, gene_cell_counts)[source]

Bases: object

Result of quality control filtering.

Parameters:

cell_mask (ndarray)
gene_mask (ndarray)
perturbation_keep (Dict[str, bool])
filtered (AnnData | None)
cell_gene_counts (ndarray)
gene_cell_counts (ndarray)

cell_mask: ndarray

gene_mask: ndarray

perturbation_keep: Dict[str, bool]

filtered: AnnData | None

cell_gene_counts: ndarray

gene_cell_counts: ndarray

property filtered_path: Path | None

Compatibility accessor exposing the on-disk filename.

Returns None if no output file was written (output_dir was None).

crispyx.qc.filter_cells_by_gene_count(data, *, min_genes=100, gene_name_column=None, chunk_size=2048, return_counts=False, return_full_result=False, verbose=True)[source]

Return a boolean mask selecting cells with at least min_genes expressed genes.

This function can optionally compute both genes-per-cell (row nnz) AND cells-per-gene (column nnz) in a single matrix pass, avoiding a separate gene counting pass later.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad file, or a crispyx AnnData / backed anndata.AnnData.
min_genes (int) – Minimum number of expressed genes per cell.
gene_name_column (str | None) – Column in var containing gene names.
chunk_size (int) – Number of cells to process per chunk.
return_counts (bool) – If True, return (mask, counts) tuple instead of just mask. Ignored if return_full_result is True.
return_full_result (bool) – If True, return a _CellFilterResult containing cell_mask, gene_counts_per_cell, and gene_cell_counts_all (cells per gene for all cells, before any perturbation filtering).
verbose (int | bool) – Print how many cells passed filtering, and warn if fewer than half did.

Returns:

Boolean mask, optionally with counts, or full result dataclass.

Return type:

mask or (mask, counts) or _CellFilterResult

crispyx.qc.filter_perturbations_by_cell_count(data, *, perturbation_column, control_label=None, min_cells=50, base_mask=None, return_counts=False, verbose=True)[source]

Return a mask keeping cells whose perturbation has sufficient representation.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad file, or a crispyx AnnData / backed anndata.AnnData.
perturbation_column (str) – Column in obs containing perturbation labels.
control_label (str | None) – Label identifying control cells. If None, auto-detected.
min_cells (int) – Minimum number of cells required per perturbation.
base_mask (ndarray | None) – Optional mask for cells to consider (e.g., from prior filtering).
return_counts (bool) – If True, return (mask, cell_counts_per_perturbation) tuple.
verbose (int | bool) – Print how many perturbations passed filtering, and warn if fewer than half did.

Returns:

Boolean mask, optionally with cell counts per perturbation label.

Return type:

mask or (mask, counts)

crispyx.qc.filter_genes_by_cell_count(data, *, min_cells=100, cell_mask=None, gene_name_column=None, chunk_size=2048, return_counts=False, verbose=True)[source]

Return a boolean mask selecting genes expressed in at least min_cells cells.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad file, or a crispyx AnnData / backed anndata.AnnData.
min_cells (int) – Minimum number of cells expressing each gene.
cell_mask (ndarray | None) – Optional mask for cells to consider.
gene_name_column (str | None) – Column in var containing gene names.
chunk_size (int) – Number of cells to process per chunk.
return_counts (bool) – If True, return (mask, counts) tuple instead of just mask.
verbose (int | bool) – Print how many genes passed filtering, and warn if fewer than half did.

Returns:

Boolean mask, optionally with the raw cell counts per gene.

Return type:

mask or (mask, counts)

crispyx.qc.quality_control_summary(data, *, min_genes=100, min_cells_per_perturbation=50, min_cells_per_gene=100, perturbation_column, control_label=None, gene_name_column=None, chunk_size=None, memory_limit_gb=None, data_name=None, output_path=None, output_dir=None, cache_mode='memmap', delta_threshold=0.3, force_streaming=False, verbose=True)[source]

Run QC with automatic strategy selection for optimal performance.

This function automatically selects the best QC strategy based on: 1. Small data (any format): In-memory processing (fastest) 2. Large CSC data: Column-oriented streaming (memory efficient for CSC) 3. Large CSR/dense data: Row-oriented streaming (current behavior)

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad file, or a crispyx AnnData / backed anndata.AnnData.
min_genes (int) – Minimum number of expressed genes per cell.
min_cells_per_perturbation (int) – Minimum number of cells required per perturbation.
min_cells_per_gene (int) – Minimum number of cells expressing each gene.
perturbation_column (str) – Column in obs containing perturbation labels.
control_label (str | None) – Label identifying control cells. If None, auto-detected.
gene_name_column (str | None) – Column in var containing gene names.
chunk_size (int | None) – Number of cells to process per chunk. If None, automatically calculated based on available memory.
memory_limit_gb (float | None) – Optional memory limit in GB for strategy selection and chunk size. If None, auto-detected from system memory.
data_name (str | None) – Base name for output files.
output_path (str | Path | None) – Exact path for the output h5ad file. When provided, output_dir and data_name are ignored.
output_dir (str | Path | None) – Directory for output files. If None, returns QC masks without writing a filtered h5ad file (QualityControlResult.filtered will be None). Deprecated; use output_path instead. Will be removed in the next major version.
cache_mode (Literal['memory', 'memmap', 'none']) – Cache strategy for row-oriented streaming: ‘memory’ (fast, high RAM), ‘memmap’ (low RAM, disk-based), or ‘none’ (no caching, requires re-reading source during write). Default is ‘memmap’ for better memory efficiency.
delta_threshold (float) – Threshold for delta adjustment in row-oriented streaming. Default 0.3 (30%).
force_streaming (bool) – If True, always use streaming path regardless of data size. Useful for testing or memory-constrained environments.
verbose (int | bool)

Returns:

Dataclass containing masks, filtered AnnData (or None if output_dir is None), and QC statistics.

Return type:

QualityControlResult

Subsampling

Streaming, stratified/cluster subsampling for large .h5ad datasets.

crispyx.sample.subsample(data, *, n=None, frac=None, groupby=None, unit='cell', drop_insufficient=True, random_state=0, chunk_size=4096, output_path=None, data_name=None, verbose=True)[source]

Stream a stratified or cluster-sampled subset of an h5ad file to disk.

Cells are never loaded to decide the subset — the sampling mask is built entirely from .obs (fast, metadata-only), then handed to crispyx.data.write_filtered_subset() to stream X, layers, obsm/obsp/varm, and uns into the output the same way every other crispyx filtering function does.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, crispyx.AnnData, or a backed AnnData.
n (int | None) – Exactly one must be given. n draws exactly this many items from every stratum; frac draws round(frac * pool_size) items, matching pandas.DataFrameGroupBy.sample(n=, frac=) semantics — both are per-stratum quantities, not a global total.
frac (float | None) – Exactly one must be given. n draws exactly this many items from every stratum; frac draws round(frac * pool_size) items, matching pandas.DataFrameGroupBy.sample(n=, frac=) semantics — both are per-stratum quantities, not a global total.
groupby (str | Sequence[str] | None) – One obs column, a sequence of obs columns, or None (default). Defines the strata that are sampled independently of each other. None means a single global stratum (the whole dataset).
unit (str) – "cell" (default) draws individual cells. Any other value must name an obs column (e.g. "batch"): each unique value of that column within a stratum becomes one indivisible sampling item — a chosen unit keeps ALL of its cells in that stratum, an unchosen unit keeps NONE. n/frac then count units, not cells.
drop_insufficient (bool) – When a stratum’s item pool is smaller than the requested count: if True (default), the stratum is excluded entirely from the output; if False, all of its cells are kept instead. Either way, every affected stratum is reported via a warning, independent of verbose.
random_state (int) – Seed for the per-stratum RNG. Sampling is deterministic for a fixed seed and independent of stratum iteration order or chunk_size.
chunk_size (int) – Cells streamed per chunk when writing the output.
output_path (str | Path | None) – Explicit output path. If None, derived from data_name.
data_name (str | None) – Custom output name suffix.
verbose (int | bool) – Print how many cells were kept and across how many strata.

Returns:

Backed AnnData wrapper pointing to the subsampled output file.

Return type:

AnnData

Pseudo-bulk aggregation

Pseudo-bulk effect size estimators operating directly on .h5ad files.

crispyx.pseudobulk.aggregate_pseudobulk(data, *, groupby, method, layer=None, gene_name_column=None, perturbations=None, min_cells=5, bootstrap_size=None, random_state=0, chunk_size=None, memory_limit_gb=None, data_name=None, output_path=None, output_dir=None, verbose=True, force=False)[source]

Aggregate cells into absolute pseudo-bulk profiles.

Each output observation represents one observed combination of the columns in groupby. For example, groupby=["perturbation", "batch"] retains repeated batch-level measurements for every perturbation instead of averaging those batches together.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, crispyx.AnnData, or a backed AnnData.
groupby (str | Sequence[str]) – One observation column, or an ordered sequence of columns defining a pseudo-bulk sample.
method (Literal['mean_log1p', 'sum']) – Required; there is no default because each choice is wrong for the other use. "sum" sums raw counts – what count models need (DESeq2, edgeR, nb_glm_test(), causarray) – and therefore requires finite, non-negative, integer-valued input. "mean_log1p" computes the mean of per-cell log1p expression, for contrasts on the log scale. Integer count input is transformed with log1p first; continuous input is treated as already log1p-transformed and averaged unchanged.
layer (str | None) – Expression layer to aggregate. Uses X when omitted.
gene_name_column (str | None) – Column in adata.var containing gene symbols. Uses var_names when omitted.
perturbations (Iterable[str] | None) – Optional subset of grouping values to retain. A profile is kept when any of its groupby values matches, so the argument selects on whichever column holds the labels regardless of its position, and every observed combination of the remaining columns is preserved. With groupby=["perturbation", "batch"], perturbations=["KO_A"] therefore returns one profile per batch for KO_A.
min_cells (int) – Minimum number of source cells required for an observed group. Defaults to 5.
bootstrap_size (int | None) – If provided, draw exactly this many cells once with replacement from every retained group before aggregation. Sampling is deterministic for a fixed random_state and independent of streaming chunk sizes.
random_state (int) – Seed for bootstrap sampling.
chunk_size (int | None) – Number of cells read per streaming block. Automatically selected when omitted.
memory_limit_gb (float | None) – Soft memory budget used for automatic chunk sizing.
data_name (str | None) – Output naming controls. output_path takes precedence.
output_path (str | Path | None) – Output naming controls. output_path takes precedence.
output_dir (str | Path | None) – Output naming controls. output_path takes precedence.
verbose (int | bool) – Level 1 reports input interpretation and save summaries.
force (bool) – Recompute an existing output instead of reloading a matching result. A result matches when it was computed from the same input file (path, size, modification time) with the same result-shaping arguments, as recorded in its uns["crispyx"].

Returns:

On-disk result with absolute profiles in X, one row per profile named by its group key (the value for one groupby column, the values joined with "|" for several, e.g. "donor1|L2/3 IT"). obs holds the grouping columns, n_cells and n_cells_aggregated, and every other source obs column that is constant within each profile (donor sex, age, disease, …); columns that vary within a profile are left out. n_cells pivoted over two grouping columns is their composition table, e.g. cells per donor and cell type. uns['crispyx_pseudobulk'] describes the aggregation and uns['crispyx'] is the provenance stamp.

Return type:

AnnData

crispyx.pseudobulk.compute_normalized_effects(data, *, perturbation_column=None, groupby=None, control_label=None, reference=None, method='mean_log1p', baseline_count=1.0, gene_name_column=None, perturbations=None, batch_column=None, chunk_size=None, memory_limit_gb=None, data_name=None, output_path=None, output_dir=None, verbose=True)[source]

Compute library-size-normalised target-minus-reference effects from cells.

This is the one-command path from a cell-level .h5ad to an effect size: it normalises each cell’s counts to 10,000, aggregates per perturbation, and subtracts the reference, in a single streaming pass over the matrix.

Unlike compute_pseudobulk_effects(), which works on whatever scale its input already carries, this function normalises library size itself. Do not pre-normalise the input as well.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, or a backed AnnData object.
perturbation_column (str | Series | None) – Column in adata.obs identifying perturbation groups, or a pd.Series of labels indexed by obs_names (nothing is written to the file). Cells whose label is missing are not used.
groupby (str | None) – Scanpy-compatible alias for perturbation_column.
control_label (str | None) – Label of the reference group. Inferred from common patterns when omitted.
reference (str | None) – Scanpy-compatible alias for control_label.
method (Literal['mean_log1p', 'log_mean']) –
Which scale the averaging happens on. Both normalise library size first:
- "mean_log1p" averages per-cell log1p values (mean of logs).
- "log_mean" averages normalised counts, then takes log1p(baseline_count * mean) (log of mean).
These are different estimators, not different spellings: averaging before or after the log gives different answers whenever expression varies within a group. "sum" is not accepted here; it belongs to aggregate_pseudobulk().
baseline_count (float) – Scaling applied inside the log for method="log_mean". Ignored by "mean_log1p". Must be positive.
gene_name_column (str | None) – Column in adata.var with gene symbols. Uses var_names when omitted.
perturbations (Iterable[str] | None) – Subset of perturbation labels to score. All non-reference groups when omitted.
batch_column (str | Series | None) – Column in adata.obs identifying each cell’s batch, or a pd.Series indexed by obs_names (cells with a missing batch are not used). Combine several columns into one beforehand (e.g. lane and cell line). When given, effects are computed within each batch and combined with harmonic target/reference cell-count weights, w_b = n_tb * n_rb / (n_tb + n_rb): the estimate of a regression with one intercept per batch and one perturbation coefficient, which removes batch-driven confounding. Batches without reference cells support no contrast and are skipped. When omitted, a single pooled effect is computed. There is no flag to disable the correction: supplying the column is the request for it.
chunk_size (int | None) – Cells streamed per chunk. Auto-selected when omitted.
memory_limit_gb (float | None) – Soft memory budget. It sizes the chunk, and decides whether the per-(perturbation, batch) accumulators (pairs x genes float64 each; two with standard errors) stay in memory, which they do when they fit in half of it, or go to temporary files on disk. Results are identical regardless.
data_name (str | None) – Custom stem for the output filename.
output_path (str | Path | None) – Exact output h5ad path. Takes precedence over output_dir / data_name.
output_dir (str | Path | None) – Deprecated output directory override; use output_path.
verbose (int | bool) – 1 / True prints a summary line.

Returns:

On-disk result with perturbations in observations and the effect in X. layers['perturbation_profile'] holds the per-perturbation profile, uns['control_profile'] the pooled reference profile, and uns['method'] the estimator used. With batch_column, X is the batch-corrected effect, layers['perturbation_profile'] its batch-corrected profile, and layers['control_profile_matched'] the weight-matched reference, so that X == perturbation_profile - control_profile_matched exactly. uns['batch_column'] and uns['batch_ids'] record the stratification. method="log_mean" additionally records uns['baseline_count'].

method="mean_log1p" also writes two standard errors of the effect (Liu et al., 2026, doi:10.64898/2026.10.05.756839), treating data without batch_column as one batch:

layers['standard_error']: cells independent within a batch, with one within-batch variance per arm pooled over the batches with at least two cells on that side (w_b-weighted population variances); SE^2 = (s2_t sum_b w_b^2 / n_tb + s2_r sum_b w_b^2 / n_rb) / W^2, W = sum_b w_b. Without batches this is sqrt(v_t / n_t + v_r / n_r). NaN when no batch has two cells on a side.
layers['standard_error_robust']: from the spread of the per-batch effects d_b around the effect E, with no shared variance; SE^2 = B / (B - 1) sum_b w_b^2 (d_b - E)^2 / W^2 over the B contributing batches. NaN with fewer than two.

Return type:

AnnData

Generic batch statistics

Generic streaming statistics stratified by experimental batch.

class crispyx.batch.BatchReducer(initialize, update, finalize, compare=None, channels=None)[source]

Bases: object

Callbacks defining a mergeable streaming statistic.

A reducer must supply three callbacks, plus a fourth required only for mode="comparison":

initialize(width) → a fresh state object.
update(state, block) → None, or a replacement state.
finalize(state) → BatchStatistic, or a plain array.
compare(group_state, reference_state) → BatchStatistic, or a plain array. Optional; omit it for mode="group".

crispyx streams the matrix one gene chunk at a time and, inside a gene chunk, one cell chunk at a time. The callbacks therefore only ever see a slice of the data and must carry whatever running state the statistic needs.

Shapes and types

widthint

Number of genes in the current gene chunk. Equal to the run’s chunk_size except for the last chunk, which may be shorter, so size state from width rather than assuming it is constant.

blocknp.ndarray, shape (n_cells, width)

Dense array holding just the cells of one (group, batch) combination that fall inside the current cell chunk. Sparse input is densified first. n_cells is at least 1 and differs between calls.

How the cells of one combination are split across calls is an implementation detail – it changes with cell_chunk_size, with how the cells happen to be ordered in the file, and between crispyx releases – so a reducer must not depend on it. The dtype likewise follows the stored matrix, which for a typical counts or normalized file is float32: accumulate in float64, or the statistic’s accuracy will vary with the block sizes crispyx happens to hand you (summing a few thousand float32 rows at once loses roughly an order of magnitude of precision against summing a few hundred). One np.asarray(block, dtype=np.float64) at the top of update is enough; the state arrays below are float64 already.

stateany object

Created fresh for each gene chunk and never shared across gene chunks. Within one gene chunk, one state is held per (group, batch) and is updated by every cell chunk containing that combination. update may mutate it in place and return None, or return a replacement.

finalize / compare return value

One value per gene, so shape (width,). Returning a bare array is equivalent to BatchStatistic(values, weight=1.0).

Examples

A within-batch mean, weighted across batches by cell count:

import numpy as np
import crispyx as cx

def initialize(width):
    # One accumulator per gene in this chunk, plus a scalar cell count.
    return {"total": np.zeros(width, dtype=np.float64), "n": 0}

def update(state, block):
    # Accumulate in float64: block follows the file's dtype, and how
    # many cells arrive per call is not part of the contract.
    block = np.asarray(block, dtype=np.float64)
    # block is (n_cells, width); summing over cells leaves (width,).
    state["total"] += block.sum(axis=0)
    state["n"] += block.shape[0]
    # Returning None keeps the state that was mutated above.

def finalize(state):
    # values must be (width,); the weight makes the cross-batch
    # average cell-count weighted.
    return cx.BatchStatistic(state["total"] / state["n"], weight=state["n"])

reducer = cx.BatchReducer(initialize, update, finalize)

Adding compare enables mode="comparison". It receives the two states accumulated independently for the group and for the reference within the same batch:

def compare(group_state, reference_state):
    difference = (
        group_state["total"] / group_state["n"]
        - reference_state["total"] / reference_state["n"]
    )
    n_group, n_reference = group_state["n"], reference_state["n"]
    weight = n_group * n_reference / (n_group + n_reference)
    return cx.BatchStatistic(difference, weight=weight)

reducer = cx.BatchReducer(initialize, update, finalize, compare)

Multiple named channels from one pass

A statistic and something derived from the same state — e.g. a mean and its standard error, so a caller can form t = mean / se — can be packed into a single streaming pass instead of two. Set channels to the tuple of names finalize/compare will return, and have them return {name: BatchStatistic | ArrayLike, ...} with exactly those keys instead of one bare value:

def compare_mean_se(group_state, reference_state):
    diff = group_state["mean"] - reference_state["mean"]
    se = np.sqrt(group_state["var"] / group_state["n"]
                 + reference_state["var"] / reference_state["n"])
    n_group, n_reference = group_state["n"], reference_state["n"]
    weight = n_group * n_reference / (n_group + n_reference)
    return {
        "mean_diff": cx.BatchStatistic(diff, weight=weight),
        "se": cx.BatchStatistic(se, weight=weight),
    }

reducer = cx.BatchReducer(
    initialize, update, finalize, compare_mean_se,
    channels=("mean_diff", "se"),
)

Each channel is combined across batches independently (with its own weight). The first channel ("mean_diff" above) is written to result.X and every later one to result.layers[name]; each channel’s weights go to result.layers[f"{name}_weight_sum"]. Leave channels as None (the default) for the single-value form above.

Notes

The statistic must be expressible through state that merges across cell chunks. One that needs all of a group’s cells at once, such as an exact median, cannot be written this way.

initialize: Callable[[int], Any]

update: Callable[[Any, ndarray], Any | None]

finalize: Callable[[Any], BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | dict[str, BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]]]

compare: Callable[[Any, Any], BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | dict[str, BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]]] | None = None

channels: tuple[str, ...] | None = None

Parameters:

initialize (Callable[[int], Any])
update (Callable[[Any, ndarray], Any | None])
finalize (Callable[[Any], BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | dict[str, BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]]])
compare (Callable[[Any, Any], BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | dict[str, BatchStatistic | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]]] | None)
channels (tuple[str, ...] | None)

class crispyx.batch.BatchStatistic(values, weight=1.0)[source]

Bases: object

A finalized within-batch statistic and its aggregation weight.

Parameters:

values (numpy._typing._array_like._Buffer | numpy._typing._array_like._SupportsArray[numpy.dtype[Any]] | numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]] | bool | int | float | complex | str | bytes | numpy._typing._nested_sequence._NestedSequence[bool | int | float | complex | str | bytes]) – One value per gene in the current gene chunk, so shape (width,). Converted to float64.
weight (float | numpy._typing._array_like._Buffer | numpy._typing._array_like._SupportsArray[numpy.dtype[Any]] | numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]] | bool | int | complex | str | bytes | numpy._typing._nested_sequence._NestedSequence[bool | int | float | complex | str | bytes]) – How much this batch counts when the per-batch statistics are combined. Either a scalar applied to every gene, or one value per gene with shape (width,). Must be finite and non-negative; a weight of zero drops the batch for the genes concerned.

Notes

Batches are combined per gene as sum(weight * values) / sum(weight). The weight is therefore what makes a cross-batch average cell-count weighted rather than a plain mean of batch means. Genes whose weights sum to zero are reported as NaN.

values: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]

weight: float | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] = 1.0

crispyx.batch.batch_process(data, reducer, *, perturbation_column=None, groupby=None, control_label=None, reference=None, gene_name_column=None, perturbations=None, batch_column, mode='group', statistic_name, chunk_size=None, cell_chunk_size=None, data_name=None, output_path=None, output_dir=None, verbose=True, memory_limit_gb=None, force=False, resume=False, checkpoint_interval=None, format_mismatch_policy='auto')[source]

Compute a generic gene-wise statistic within experimental batches.

Parameters follow the differential-expression API: groupby aliases perturbation_column and reference aliases control_label. chunk_size controls genes per chunk, while cell_chunk_size controls rows supplied to the reducer at once.

In "group" mode, finalize is called for every observed (group, batch) state. In "comparison" mode, compare is called for group/reference states from shared batches. Finalized batch statistics are combined as sum(weight * values) / sum(weight).

The matrix is streamed gene-chunk by gene-chunk in a single pass (crispyx.data.iter_matrix_chunks(), axis=1) – the same access pattern wilcoxon_test uses. This is native and cheap for a CSC-stored source; see format_mismatch_policy for a CSR-stored one. Each gene chunk’s result is also the unit of resumable progress; see resume.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an existing h5ad file, crispyx.AnnData, or a backed anndata.AnnData. A normal in-memory anndata.AnnData is not accepted because the operation is designed to stream from disk. X must be a numeric cells-by-genes matrix, and obs must contain the grouping and batch columns.
reducer (BatchReducer) – A BatchReducer whose state can be merged across cell chunks. initialize(width) receives the current gene-chunk width; update(state, block) receives a dense cells-by-width array; and finalize(state) must return one value per gene, or – with reducer.channels set – a dict of named channels (see BatchReducer). Comparison mode also requires compare(group_state, reference_state).
perturbation_column (str | None) – Column in adata.obs containing group labels.
groupby (str | None) – Scanpy-compatible alias for perturbation_column. The two arguments are mutually exclusive.
control_label (str | None) – Reference group used in comparison mode. If omitted, crispyx infers a common control label.
reference (str | None) – Scanpy-compatible alias for control_label.
gene_name_column (str | None) – Column in adata.var containing gene symbols. Uses var_names when omitted.
perturbations (Iterable[str] | None) – Ordered subset of groups to return. All applicable groups are used when omitted.
batch_column (str) – Column in adata.obs containing experimental batch/stratum labels.
mode (Literal['group', 'comparison']) – "group" finalizes each group independently. "comparison" compares each group with the reference within shared batches.
statistic_name (str) – Stable identifier used in output naming and cache metadata.
chunk_size (int | None) – Genes processed per chunk. Automatically selected when omitted, so that both the densified cell slabs and the per-chunk (n_groups, chunk_size) accumulators (two float64 arrays per output layer) fit the memory budget. Also the unit of resumable progress – see resume.
cell_chunk_size (int | None) – Cells densified at a time within a gene chunk, i.e. the largest block a reducer update receives. Each slab is a dense cell_chunk_size × chunk_size float64 array, so this is the per-chunk working set. Automatically selected when omitted.
data_name (str | None) – Optional input stem override used to construct the output filename.
output_path (str | Path | None) – Exact output h5ad path. Takes precedence over naming arguments.
output_dir (str | Path | None) – Deprecated output directory override; use output_path instead.
verbose (int | bool) – Verbosity level. Level 1 prints run/reload and save summaries.
memory_limit_gb (float | None) – Soft memory budget used for automatic chunk sizing.
force (bool) – Recompute even when a matching output already exists. Edits to the input file are detected automatically through its path and modification time, so a regenerated source invalidates a cached result on its own. Use force when a reducer’s implementation changes without changing statistic_name, since the reducer itself cannot be fingerprinted.
resume (bool) – If True, resume an interrupted run from its checkpoint (<output>.progress.json), or – if that checkpoint is missing or corrupted – by scanning the partially-written output file itself for the last completed gene chunk. If False (default), always start over.
checkpoint_interval (int | None) – Gene chunks between checkpoint writes. Auto-selected from the number of gene chunks when omitted.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) –
How to handle a source stored as CSR, whose gene-(column-)streaming here re-reads the whole matrix once per gene chunk and holds it in memory – typically ~100x more I/O than CSC:
- "auto" (default): measure the source’s read throughput and convert only when the re-reads would actually cost more than one conversion does – roughly, a large file on a slow (networked) filesystem streamed in many chunks. A file that re-reads in milliseconds is streamed as-is, silently. See crispyx.data.resolve_auto_format_mismatch_policy().
- "convert": always convert the source to CSC in a temporary file beside the output (bounded-memory streaming via crispyx.data.convert_to_csc(), honouring memory_limit_gb) and stream from that; the temporary file is removed before returning. Needs ~2x the source file’s size in free disk space there; when that is not available the call falls back to "warn" behaviour with a warning naming the shortfall. Run cx.pp.convert_to_csc once instead if several steps will reuse the file.
- "warn": proceed on the CSR source after one warning that quantifies the cost.
- "off": proceed silently with no warning.

Returns:

On-disk result with groups in observations, genes in variables, the corrected statistic in X, and accumulated weights in layers['weight_sum']. With reducer.channels set, the first channel’s combined values are X and every later channel’s are layers[name]; each channel’s weights are layers[f"{name}_weight_sum"]. uns['channels'] lists the channels in order.

Return type:

AnnData

Notes

The reducer sees the values stored in X exactly as provided. Normalize or log-transform the input beforehand when the statistic requires it. A statistic must be expressible through mergeable state; callbacks that need all cells simultaneously are not suitable for this streaming interface.

Examples

A streaming mean keeps only a sum and count. Returning the count as the weight produces a cell-count-weighted mean across biological batches:

import numpy as np
import crispyx as cx

def initialize_mean(width):
    return {"sum": np.zeros(width, dtype=np.float64), "count": 0}

def update_mean(state, block):
    # Accumulate in float64: block follows the file's dtype, so
    # block.sum() on a float32 file would sum in float32.
    state["sum"] += np.asarray(block, dtype=np.float64).sum(axis=0)
    state["count"] += block.shape[0]

def finalize_mean(state):
    values = state["sum"] / state["count"]
    return cx.BatchStatistic(values, weight=state["count"])

mean_reducer = cx.BatchReducer(
    initialize=initialize_mean,
    update=update_mean,
    finalize=finalize_mean,
)
result = cx.tl.batch_process(
    "screen.h5ad",
    mean_reducer,
    groupby="perturbation",
    batch_column="batch",
    statistic_name="mean",
)

Comparison mode adds a callback receiving the independently accumulated group and reference states. This example computes within-batch mean differences with the harmonic-count weighting used by crispyx pseudobulk:

def compare_means(group_state, reference_state):
    n_group = group_state["count"]
    n_reference = reference_state["count"]
    difference = (
        group_state["sum"] / n_group
        - reference_state["sum"] / n_reference
    )
    weight = n_group * n_reference / (n_group + n_reference)
    return cx.BatchStatistic(difference, weight=weight)

comparison_reducer = cx.BatchReducer(
    initialize=initialize_mean,
    update=update_mean,
    finalize=finalize_mean,
    compare=compare_means,
)
result = cx.tl.batch_process(
    "screen.h5ad",
    comparison_reducer,
    groupby="perturbation",
    reference="control",
    batch_column="batch",
    mode="comparison",
    statistic_name="mean_difference",
)

Differential expression

All three DE functions reload an existing output .h5ad instead of rerunning when it was computed from the same input and arguments; pass force=True to recompute. See Auto-reload and re-run control.

Differential expression testing utilities.

class crispyx.de.DifferentialExpressionResult(genes, logfoldchanges, statistic, pvalue, method, perturbation, pvalue_adj=None, auc=None, mean_difference=None, result=None)[source]

Bases: object

One perturbation’s per-gene results; see RankGenesGroupsResult for what each field holds.

Parameters:

genes (Index)
logfoldchanges (ndarray)
statistic (ndarray)
pvalue (ndarray)
method (str)
perturbation (str)
pvalue_adj (ndarray | None)
auc (ndarray | None)
mean_difference (ndarray | None)
result (AnnData | None)

genes: Index

logfoldchanges: ndarray

statistic: ndarray

pvalue: ndarray

method: str

perturbation: str

pvalue_adj: ndarray | None = None

auc: ndarray | None = None

mean_difference: ndarray | None = None

result: AnnData | None = None

property result_path: Path

Path to the on-disk result file.

Returns:: Absolute path to the .h5ad file backing the result.
Return type:: Path
Raises:: AttributeError – If the result AnnData has not been initialised.

class crispyx.de.RankGenesGroupsResult(genes, groups, statistics, pvalues, pvalues_adj, logfoldchanges, u_statistics, pts, pts_rest, order, groupby, method, control_label, tie_correct, pvalue_correction, auc=None, mean_difference=None, result=None)[source]

Bases: Mapping[str, DifferentialExpressionResult]

Per-perturbation x per-gene DE results; rows follow groups.

logfoldchanges is the log2 fold change of mean expression for every method: for wilcoxon_test and t_test, log2((expm1(mean_pert) + 1e-9) / (expm1(mean_ctrl) + 1e-9)) on the log1p input (Scanpy’s convention); for nb_glm_test, the fitted perturbation coefficient on the log2 scale (lfc_base="ln": natural log). The method-specific effect sizes are auc (wilcoxon_test: U / (n_pert * n_ctrl), the probability that a perturbed cell exceeds a control cell, 0.5 for no shift) and mean_difference (t_test: the difference of means on the input’s own scale); each is None for the other methods.

Parameters:

genes (Index)
groups (list[str])
statistics (ndarray)
pvalues (ndarray)
pvalues_adj (ndarray)
logfoldchanges (ndarray)
u_statistics (ndarray)
pts (ndarray)
pts_rest (ndarray)
order (ndarray)
groupby (str)
method (str)
control_label (str)
tie_correct (bool)
pvalue_correction (Literal['benjamini-hochberg', 'bonferroni'])
auc (ndarray | None)
mean_difference (ndarray | None)
result (AnnData | None)

genes: Index

groups: list[str]

statistics: ndarray

pvalues: ndarray

pvalues_adj: ndarray

logfoldchanges: ndarray

u_statistics: ndarray

pts: ndarray

pts_rest: ndarray

order: ndarray

groupby: str

method: str

control_label: str

tie_correct: bool

pvalue_correction: Literal['benjamini-hochberg', 'bonferroni']

auc: ndarray | None = None

mean_difference: ndarray | None = None

result: AnnData | None = None

property result_path: Path

items()[source]: Return (group_name, DifferentialExpressionResult) pairs.

to_rank_genes_groups_dict()[source]

Convert results to Scanpy-compatible rank_genes_groups dict.

Returns:: Dictionary with keys 'params', 'names', 'scores', 'logfoldchanges', 'pvals', 'pvals_adj', 'pts', 'pts_rest', 'u_stat', 'full', and 'auc' when the method has one (wilcoxon_test). Each value (except 'params' and 'full') is a NumPy record array sorted by descending absolute score.
Return type:: dict

to_full_order_dict()[source]

Return unsorted matrices keyed by statistic name.

Unlike to_rank_genes_groups_dict(), values are in the original gene order (not sorted by score).

Returns:: Keys: 'scores', 'pvals', 'pvals_adj', 'logfoldchanges', 'u_stat', 'pts', 'pts_rest', and 'auc' when the method has one. Each value is a 2-D (n_groups, n_genes) NumPy array (copy).
Return type:: dict

crispyx.de.t_test(data, *, perturbation_column=None, groupby=None, control_label=None, reference=None, gene_name_column=None, perturbations=None, min_cells_expressed=0, min_pct_ctrl=0.01, min_pct_pert=0.002, min_mean_ctrl=0.05, min_mean_pert=0.005, cell_chunk_size=None, data_name=None, output_path=None, output_dir=None, n_jobs=None, verbose=True, resume=False, checkpoint_interval=None, corr_method='benjamini-hochberg', scanpy_format=False, memory_limit_gb=None, force=False)[source]

Perform a t-test comparing log-expression means for each perturbation.

Returns a RankGenesGroupsResult containing differential expression statistics. Results are stored in an h5ad file with layers containing the statistics.

The RankGenesGroupsResult implements the Mapping interface, so it can be used like a dict: result[perturbation_label] returns a DifferentialExpressionResult for that perturbation.

Input data should already be normalized and log-transformed (for example using scanpy.pp.normalize_total followed by scanpy.pp.log1p). To maintain backward compatibility with Scanpy-style workflows, count-like inputs are automatically normalized and log-transformed in streaming fashion, with a warning to encourage explicit preprocessing upstream.

A perturbation (or control) with a single cell has no variance estimate, so its comparison is left untested (NaN), as in scanpy.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, or a crispyx AnnData / backed anndata.AnnData containing log-transformed expression data.
perturbation_column (str | None) – Column in adata.obs indicating perturbation labels.
groupby (str | None) – Alias for perturbation_column (Scanpy-compatible). Mutually exclusive with perturbation_column.
control_label (str | None) – Label for the control/reference group. If None, infers from common patterns.
reference (str | None) – Alias for control_label (Scanpy-compatible). Mutually exclusive with control_label.
gene_name_column (str | None) – Column in adata.var with gene symbols. If None, uses adata.var_names.
perturbations (Iterable[str] | None) – Specific perturbations to test. If None, tests all non-control groups.
min_cells_expressed (int) – Minimum total cells (control + perturbation) expressing a gene for testing. Genes below this threshold are untested: NaN in score / pvalue / logfoldchanges and the effect size, with pts / pts_rest still populated.
min_pct_ctrl (float) – Minimum fraction of expressing cells for the control side. A gene is excluded only when both the control side and the perturbed side are jointly low. Default 0.01.
min_pct_pert (float) – Minimum fraction of expressing cells for the perturbed side. Default 0.002 (lower than ctrl; induction from near-zero baseline is biologically valid). Set to 0.0 to disable the pct check on pert.
min_mean_ctrl (float) – Minimum mean expression (log1p units) for the control side. Default 0.05. Excluded genes are written as NaN in score / pvalue / logfoldchanges / mean_difference; pts and pts_rest remain populated.
min_mean_pert (float) – Minimum mean expression for the perturbed side. Default 0.005. Together with min_pct_pert this forms a dual condition that is more robust to doublet / ambient-RNA artefacts than pct alone.
cell_chunk_size (int | None) – Number of cells to process per chunk (memory vs. speed tradeoff). This controls streaming along the cell axis and is distinct from any future perturbation_chunk_size option that would batch perturbations. Data must already be normalized/log-transformed before chunking.
data_name (str | None) – Custom name for output file. If None, uses “t_test” suffix.
output_path (str | Path | None) – Exact path for the output h5ad file. When provided, output_dir and data_name are ignored.
output_dir (str | Path | None) – Directory for output h5ad file. Defaults to input file’s directory. Deprecated; use output_path instead. Will be removed in the next major version.
n_jobs (int | None) – Number of parallel workers for computing statistics across perturbations. If None, uses all available cores. If 1, runs sequentially. If -1, uses all available cores.
verbose (int | bool) – If True, show a progress bar for perturbation processing. Requires tqdm.
resume (bool) – If True, continue an interrupted run from its last checkpoint, skipping the perturbations it had finished. Only a checkpoint from a call with the same input and result-affecting parameters is used; partial results live in a hidden .<output name>.resume directory until the output is written.
checkpoint_interval (int | None) – Number of perturbations between checkpoint saves. Auto-determined if None.
corr_method (Literal['benjamini-hochberg', 'bonferroni']) – Method for p-value correction: "benjamini-hochberg" or "bonferroni".
scanpy_format (bool) – If True, write Scanpy-compatible uns['rank_genes_groups'] structure in addition to the layer-based storage. Adds ~2-6 seconds of I/O overhead for large datasets. Default False for performance.
memory_limit_gb (float | None) – Maximum memory budget in GB. Used for automatic chunk size calculation. When None (default), detects available system memory via psutil. For HPC environments, set this to your SLURM --mem value (e.g., memory_limit_gb=128).
force (bool) – If True, rerun the analysis even when the output h5ad file already exists. If False (default), an existing result is loaded instead of rerun when it was computed from the same input file (path, size, modification time) with the same result-shaping arguments, as recorded in its uns["crispyx"]; otherwise it is recomputed.

Returns:

Differential expression results. Access results via dict-like interface: result[label].logfoldchanges, result[label].mean_difference, result[label].pvalue, etc. The h5ad file path is available at result.result_path.

The h5ad file holds one row per perturbation: the difference of means (mean_difference) in X and the layers z_score, pvalue, pvalue_adj, logfoldchanges and pts. pts_rest describes the control arm, so it is the per-gene column var["pts_rest"]. obs["n_genes_tested"] counts the genes each perturbation was tested on, after the filters.

Return type:

Notes

logfoldchanges follows scanpy’s formula, log2((expm1(mean_group) + 1e-9) / (expm1(mean_rest) + 1e-9)). When a gene has no expressing cells in one arm, that arm’s term is the constant, so the magnitude reported is set by the constant – its ceiling is log2(1 / 1e-9) = 29.9 – and not by the data. The p-values are unaffected, but ranking on logfoldchanges puts such genes above every real effect. pts / pts_rest tell them apart: a complete knockdown has pts = 0 beside a high pts_rest, an undetectable gene has both near zero.

crispyx.de.nb_glm_test(data, *, perturbation_column=None, groupby=None, control_label=None, reference=None, covariates=None, gene_name_column=None, perturbations=None, size_factors=None, size_factor_method='sparse', size_factor_scope='global', scale_size_factors=True, dispersion=None, dispersion_method='cox-reid', dispersion_scope='global', use_map_dispersion=True, shrink_dispersion=True, optimization_method='lbfgsb', max_iter=25, tol=1e-06, min_mu=0.5, poisson_init_iter=5, chunk_size=None, irls_batch_size=128, min_cells_expressed=0, min_pct_ctrl=0.01, min_pct_pert=0.002, min_cells_ctrl=1, min_cells_pert=1, min_mean_ctrl=0.05, min_mean_pert=0.005, min_total_count=1.0, cook_filter=False, lfc_shrinkage_type='none', lfc_base='log2', corr_method='benjamini-hochberg', se_method='sandwich', data_name=None, output_path=None, output_dir=None, scanpy_format=False, verbose=True, profiling=False, resume=False, checkpoint_interval=None, memory_limit_gb=None, max_dense_fraction=0.3, n_jobs=None, use_control_cache=True, freeze_control=None, force=False)[source]

Perform negative binomial GLM differential expression test.

Returns a RankGenesGroupsResult containing differential expression statistics. Uses a negative binomial GLM framework that can incorporate covariates. Results are stored in an h5ad file with layers containing the statistics.

The RankGenesGroupsResult implements the Mapping interface, so it can be used like a dict: result[perturbation_label] returns a DifferentialExpressionResult for that perturbation.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, or a crispyx AnnData / backed anndata.AnnData containing raw count data.
perturbation_column (str | None) – Column in adata.obs indicating perturbation labels.
groupby (str | None) – Alias for perturbation_column (Scanpy-compatible). Mutually exclusive with perturbation_column.
control_label (str | None) – Label for the control/reference group. If None, infers from common patterns.
reference (str | None) – Alias for control_label (Scanpy-compatible). Mutually exclusive with control_label.
covariates (Iterable[str] | None) – Additional columns in adata.obs to include as covariates in the GLM.
gene_name_column (str | None) – Column in adata.var with gene symbols. If None, uses adata.var_names.
perturbations (Iterable[str] | None) – Specific perturbations to test. If None, tests all non-control groups.
size_factors (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – Optional array of per-cell size factors. If None, computes size factors using the method specified by size_factor_method.
size_factor_method (Literal['sparse', 'deseq2']) –
Method for computing size factors when size_factors is None.
- "sparse": Sparse-aware median-of-ratios (default). Computes geometric means using only non-zero values, suitable for sparse single-cell data.
- "deseq2": Classic DESeq2/PyDESeq2 style. Uses only genes expressed in ALL cells (typically ~50-100 genes). Provides better numerical alignment with PyDESeq2 results but may be less robust for very sparse data.
size_factor_scope (Literal['global', 'per_comparison']) –
Scope for size factor computation.
- "global" (default): Compute size factors once on the full dataset. Recommended for CRISPR screens where all cells come from the same experiment and share a common sequencing depth distribution. Faster when combined with use_control_cache=True. Note: produces different results from PyDESeq2 (rho ~ 0.7-0.8) which uses per-comparison normalization.
- "per_comparison": Compute size factors separately for each control + perturbation comparison. This matches PyDESeq2’s behavior exactly, leading to near-perfect LFC, statistic, and p-value concordance (rho > 0.97 on Tian-crispra, rho > 0.99 on Adamson_subset). Use this when PyDESeq2 compatibility is required or for bulk RNA-seq style analysis.
scale_size_factors (bool) – If True (default), scale size factors so their geometric mean equals 1. This is the standard DESeq2/crispyx behavior. If False, use raw median-of-ratios size factors without rescaling, which matches PyDESeq2’s default behavior and can improve numerical alignment.
dispersion (float | None) – Fixed NB dispersion used for every gene and comparison. If None, it is estimated (see dispersion_scope).
dispersion_method (Literal['moments', 'cox-reid']) –
Method for estimating dispersion when dispersion is None.
- "moments": Method-of-moments (fast but less accurate).
- "cox-reid": Cox-Reid adjusted profile likelihood (slower but more accurate, similar to DESeq2). This is the default.
dispersion_scope (Literal['global', 'per_comparison']) –
Scope for dispersion estimation.
- "global" (default): Estimate each gene’s dispersion once and use it in every comparison. Without covariates it is estimated from all cells; with covariates, from the control cells fitted with the covariates (DESeq2’s MAP step), and each comparison is then one fit at that fixed dispersion. This is many times faster than per-comparison estimation and suits CRISPR screens, where perturbation effects are small relative to baseline expression. Per-comparison size factors keep per-comparison dispersion.
- "per_comparison": Estimate dispersion separately for each control + perturbation comparison. More accurate when perturbations cause large changes in gene expression variance, but significantly slower.
use_map_dispersion (bool) – If True (default), use MAP dispersion estimation with mean-dispersion trend. If False, use MLE dispersion without trend-based shrinkage.
shrink_dispersion (bool) – If True, fit a mean-dispersion trend and shrink gene-wise dispersions toward the trend using an empirical Bayes prior.
optimization_method (Literal['irls', 'lbfgsb']) –
Method for coefficient optimization.
- "lbfgsb": L-BFGS-B optimization (PyDESeq2 style, default). Directly optimizes the negative binomial log-likelihood.
- "irls": Iteratively Reweighted Least Squares (Fisher scoring). The classic GLM fitting approach.
max_iter (int) – Maximum iterations for GLM fitting.
tol (float) – Convergence tolerance for GLM fitting.
min_mu (float) – Minimum mean threshold for IRLS numerical stability. Predicted means are clamped to max(min_mu, predicted) during fitting to prevent numerical instability from very small means. Default: 0.5, matching PyDESeq2’s default. Set to 0.0 to disable clamping.
poisson_init_iter (int) – Initial Poisson iterations before switching to negative binomial.
chunk_size (int | None) – Number of genes to process per chunk (memory vs. speed tradeoff). Smaller values stream more, reducing peak memory at the cost of additional I/O.
irls_batch_size (int | None) – Maximum number of genes to densify per IRLS step. Keep this small to limit per-iteration memory when working with large sparse matrices. Set to None to process each chunk without additional batching.
min_cells_expressed (int) – Minimum total cells (control + perturbation) expressing a gene for testing. Genes below this threshold are untested: NaN in score / pvalue / logfoldchanges and the effect size, with pts / pts_rest still populated.
min_total_count (float) – Minimum total count across all cells for a gene to be tested.
min_pct_ctrl (float) – Minimum fraction of expressing cells for the control side. A gene is excluded only when both sides are jointly low. Default 0.01.
min_pct_pert (float) – Minimum fraction of expressing cells for the perturbed side. Default 0.002. Combined with min_mean_pert this forms a dual condition that is more robust to doublet / ambient-RNA artefacts.
min_mean_ctrl (float) – Minimum mean (size-factor-normalised) expression for the control side. Default 0.05. Excluded genes appear as NaN in pvalue / effect / logfc / se; pts and mean remain populated.
min_mean_pert (float) – Minimum mean expression for the perturbed side. Default 0.005.
min_cells_ctrl (int) –
Minimum number of expressing cells required in the control and perturbed arms for a (gene, perturbation) pair to be tested. Both default to 1: a pair with no counts at all in one arm has no finite effect to estimate, so it is reported as untested (NaN effect, statistic and p-value) and excluded from the multiple-testing correction rather than given a bound-determined estimate. See crispyx._statistics._nonestimable_glm_mask() for why, and for why the two sides are set separately (CRISPRi and CRISPRa want opposite asymmetries). Set either to 0 to disable that side.

The observation itself is not lost: pts and pts_rest are populated for every gene, so a pair expressed in 0% of perturbed and 95% of control cells is plainly visible there. lfc_shrinkage_type does not change this: an excluded pair is dropped before any fit, so it stays NaN under "apeglm" as well.
min_cells_pert (int) –
Minimum number of expressing cells required in the control and perturbed arms for a (gene, perturbation) pair to be tested. Both default to 1: a pair with no counts at all in one arm has no finite effect to estimate, so it is reported as untested (NaN effect, statistic and p-value) and excluded from the multiple-testing correction rather than given a bound-determined estimate. See crispyx._statistics._nonestimable_glm_mask() for why, and for why the two sides are set separately (CRISPRi and CRISPRa want opposite asymmetries). Set either to 0 to disable that side.

The observation itself is not lost: pts and pts_rest are populated for every gene, so a pair expressed in 0% of perturbed and 95% of control cells is plainly visible there. lfc_shrinkage_type does not change this: an excluded pair is dropped before any fit, so it stays NaN under "apeglm" as well.
cook_filter (bool) – Whether to apply Cook’s distance outlier filtering when available.
lfc_shrinkage_type (Literal['apeglm', 'none']) –
Type of log-fold change shrinkage to apply.
- "none": No shrinkage (default).
- "apeglm": Adaptive shrinkage using Cauchy prior (PyDESeq2-compatible). Preserves large effects while shrinking small/uncertain effects toward zero. Also updates standard errors to reflect posterior uncertainty.
lfc_base (Literal['log2', 'ln']) –
Log base for fold change output.
- "log2" (default): Output log2 fold change, matching PyDESeq2/edgeR.
- "ln": Output natural log fold change (raw GLM coefficients).
Standard error is also converted to match the selected log base. Wald statistics remain unchanged since both LFC and SE are scaled equally.
corr_method (Literal['benjamini-hochberg', 'bonferroni']) – Method for p-value correction: "benjamini-hochberg" or "bonferroni".
se_method (Literal['sandwich', 'fisher']) –
Method for computing standard errors.
- "sandwich" (default): Sandwich estimator SE = sqrt(c’ @ H @ M @ H @ c). More robust to model misspecification.
- "fisher": Standard Fisher information SE = sqrt(diag(inv(X’WX + ridge*I))). Matches PyDESeq2’s approach for better p-value parity.
data_name (str | None) – Custom name for output file. If None, uses “nb_glm” suffix.
output_path (str | Path | None) – Exact path for the output h5ad file. When provided, output_dir and data_name are ignored.
output_dir (str | Path | None) – Directory for output h5ad file. Defaults to input file’s directory. Deprecated; use output_path instead. Will be removed in the next major version.
scanpy_format (bool) – If True, write Scanpy-compatible uns['rank_genes_groups'] structure in addition to the layer-based storage. Adds ~2-6 seconds of I/O overhead for large datasets. Default False for performance.
verbose (int | bool) – If True, show a progress bar for perturbation fitting and log per-perturbation completion at DEBUG level. Requires tqdm to be installed for progress bar.
profiling (bool) – If True, enable timing and memory profiling. When enabled, stores profiling data in adata.uns["profiling"] with fields fit_seconds, fit_peak_memory_mb, and profiling_enabled. When False (default), adata.uns["profiling"] is set to "NA" to avoid profiling overhead in production.
resume (bool) – If True, continue an interrupted run from its last checkpoint, skipping the perturbations it had finished. Only a checkpoint from a call with the same input and result-affecting parameters is used; partial results live in a hidden .<output name>.resume directory until the output is written.
checkpoint_interval (int | None) – Number of perturbations to process between checkpoint saves. If None, auto-determined based on dataset size (1 for <100 perturbations, 10 for <1000, 50 for larger).
memory_limit_gb (float | None) – Optional memory limit in GB. If provided, this is used together with available system memory to determine when to switch to streaming mode for global dispersion estimation. The effective limit is min(available_memory, memory_limit_gb). Default is None (use system memory only).
max_dense_fraction (float) – Maximum fraction of available memory to use for dense matrix operations. If the estimated memory for densifying the full cell×gene matrix exceeds max_dense_fraction × min(available_memory, memory_limit_gb), the function switches to streaming mode. Default is 0.3 (30% of available memory).
n_jobs (int | None) – Number of parallel workers for fitting GLMs across perturbations. If None, uses all available cores. If 1, runs sequentially. If -1, uses all available cores.
use_control_cache (bool) – If True (default), precompute control cell statistics (intercept, weights, XᵀWX contributions) once and reuse them across all perturbation comparisons. This can significantly reduce memory and computation time when there are many perturbations and the control group is large. Only applies when no covariates are specified and size_factor_scope=”global”.
freeze_control (bool | None) –
Whether to use frozen control sufficient statistics instead of raw control matrix for parallel fitting. This dramatically reduces per-worker memory from ~5GB to ~1MB for large datasets, enabling more parallel workers.
- None (default): Auto-detect based on dataset size. Frozen control is enabled when control_matrix serialization would limit workers to <4, AND the required settings (dispersion_scope=’global’, shrink_dispersion=True) are met. For most large datasets (>500K cells), this auto-enables.
- True: Force frozen control mode. Raises ValueError if requirements not met.
- False: Disable frozen control (use raw control matrix).
The mode used is recorded in uns["frozen_control"] of the result. With None, an existing result is reused whichever mode the run that wrote it chose; pass True or False to require one.

Memory efficiency: Per-worker pickle size is reduced from (control_n × n_genes × 8) bytes to just ~1MB of sufficient statistics (W_sum, Wz_sum arrays).

Example: For Replogle-GW-k562 (75K control cells × 8K genes): - Without freeze_control: ~4.7 GB per worker → 2 workers max @ 128GB - With freeze_control: ~1 MB per worker → 32 workers @ 128GB - Time reduction: ~300 hours → ~10 hours

Requirements (enforced when True, auto-checked when None): - dispersion_scope=”global” (per-comparison dispersion needs raw matrix) - shrink_dispersion=True (ensures global dispersion is computed) - use_control_cache=True (required for caching)

Technical note: The intercept (β₀) is frozen to the control-only estimate. This is valid because control cells have perturbation indicator = 0, so μ_control = exp(β₀ + offset) is independent of the perturbation effect.
force (bool) – If True, rerun the analysis even when the output h5ad file already exists. If False (default), an existing result is loaded instead of rerun when it was computed from the same input file (path, size, modification time) with the same result-shaping arguments, as recorded in its uns["crispyx"]; otherwise it is recomputed.

Returns:

Differential expression results. Access results via dict-like interface: result[label].logfoldchanges, result[label].pvalue, etc. The h5ad file path is available at result.result_path.

The h5ad file holds one row per perturbation. X is the fold change in lfc_base (shrunk under lfc_shrinkage_type="apeglm"). Layers: z_score, pvalue, pvalue_adj, standard_error, intercept (ln scale), pts, converged (bool) and iterations (integer); logfoldchange_raw holds the MLE fold change only when X is shrunk, and with a log2 base logfoldchange_raw_ln / standard_error_ln keep the ln-scale values shrink_lfc() needs. Per-gene quantities are var columns: pts_rest, and, when a global dispersion is fitted, dispersion and dispersion_trend; a per-comparison dispersion is stored as the layers dispersion, dispersion_raw and dispersion_trend instead. obs["n_genes_tested"] counts the genes each perturbation was tested on, after the filters.

Return type:

crispyx.de.wilcoxon_test(data, *, perturbation_column=None, groupby=None, control_label=None, reference=None, gene_name_column=None, perturbations=None, batch_column=None, min_cells_expressed=0, min_pct_ctrl=0.01, min_pct_pert=0.002, min_mean_ctrl=0.05, min_mean_pert=0.005, chunk_size=None, tie_correct=True, corr_method='benjamini-hochberg', data_name=None, output_path=None, output_dir=None, verbose=True, resume=False, checkpoint_interval=None, scanpy_format=False, memory_limit_gb=None, force=False, format_mismatch_policy='auto')[source]

Perform a Wilcoxon rank-sum (Mann-Whitney U) test for each gene.

Input data must already be library-size normalised and log-transformed. The function operates directly on the provided matrix without additional preprocessing. As a safeguard, the first sparse chunk is inspected and a warning is emitted if the data appear to be raw counts (integer or count-like floats), encouraging explicit preprocessing upstream.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, or a crispyx AnnData / backed anndata.AnnData containing normalised, log-transformed data.
perturbation_column (str | Series | None) – Column in adata.obs indicating perturbation labels, or the labels themselves as a pd.Series indexed by cell (obs_names), e.g. shuffled labels for a permutation null, so no relabelled copy of the data has to be written. The Series is aligned to obs_names by index, so its order does not matter, but every cell needs a label; its name (default "perturbation") names the column in the result. A different Series is a different call: a result computed with one is not reused for another. On a CSR input each call makes a temporary CSC copy (see format_mismatch_policy); for many shuffles, convert once with crispyx.pp.convert_to_csc() and pass the converted file.
groupby (str | Series | None) – Alias for perturbation_column (Scanpy-compatible). Mutually exclusive with perturbation_column.
control_label (str | None) – Label for the control/reference group. If None, infers from common patterns.
reference (str | None) – Alias for control_label (Scanpy-compatible). Mutually exclusive with control_label.
gene_name_column (str | None) – Column in adata.var with gene symbols. If None, uses adata.var_names.
perturbations (Iterable[str] | None) – Specific perturbations to test. If None, tests all non-control groups.
batch_column (str | None) – Column in adata.obs identifying the batch / stratum of each cell (e.g. 10x lane / gem-group). When provided, a batch-stratified (van Elteren) Wilcoxon test is performed: cells are ranked within each batch separately and the per-stratum rank statistics are combined with unit weights (equivalent to summing the Mann-Whitney U statistics across strata). This removes rank inflation caused by batch effects. The low-expression filter, log-fold changes and pts remain pooled across all cells; only the rank test is stratified. When None (default) the standard pooled Wilcoxon test is used.
min_cells_expressed (int) – Minimum total cells (control + perturbation) expressing a gene for testing. Genes below this threshold are untested: NaN in score / pvalue / logfoldchanges and the effect size, with pts / pts_rest still populated. A gene whose values are all tied across both groups (e.g. zero in every cell, with the filters off) has no rank test and is untested in the same way; Scanpy reports score=0, pvalue=1 for it instead.
min_pct_ctrl (float) – Minimum fraction of expressing cells for the control side. A gene is excluded only when both sides are jointly low. Default 0.01.
min_pct_pert (float) – Minimum fraction of expressing cells for the perturbed side. Default 0.002 (lower than ctrl; induction from near-zero baseline is biologically valid). Combined with min_mean_pert this forms a dual condition more robust than pct alone.
min_mean_ctrl (float) – Minimum mean log1p expression for the control side. Default 0.05. Excluded genes are written as NaN in score / pvalue / logfoldchanges / auc; pts / pts_rest remain populated. Set to 0.0 together with pct thresholds to disable.
min_mean_pert (float) – Minimum mean expression for the perturbed side. Default 0.005.
chunk_size (int | None) – Number of genes to process per chunk (memory vs. speed tradeoff). Smaller values stream more, reducing peak memory at the cost of additional I/O.
tie_correct (bool) – Whether to apply tie correction to the U statistic. Default True for more accurate p-values when ties are present in the data.
corr_method (Literal['benjamini-hochberg', 'bonferroni']) – Method for p-value correction: “benjamini-hochberg” or “bonferroni”.
data_name (str | None) – Custom name for output file. If None, uses “wilcoxon” suffix.
output_path (str | Path | None) – Exact path for the output h5ad file. When provided, output_dir and data_name are ignored.
output_dir (str | Path | None) – Directory for output h5ad file. Defaults to input file’s directory. Deprecated; use output_path instead. Will be removed in the next major version.
verbose (int | bool) – If True, show a progress bar for gene chunk processing. Requires tqdm.
resume (bool) – If True, continue an interrupted run from its last checkpoint, at the first unfinished gene chunk (or perturbation batch, on the streaming path). The resolved chunk_size is part of the call’s identity, so pass the same value (or leave both on auto with the same memory). Only a checkpoint from a call with the same input and result-affecting parameters is used; partial results live in a hidden .<output name>.resume directory until the output is written.
checkpoint_interval (int | None) – Number of gene chunks between checkpoint saves. If None, auto-determined based on dataset size.
scanpy_format (bool) – If True, write Scanpy-compatible uns['rank_genes_groups'] structure in addition to the layer-based storage. Adds ~2-6 seconds of I/O overhead for large datasets. Default False for performance.
memory_limit_gb (float | None) – Maximum memory budget in GB. Controls whether the streaming path is used for very large datasets. When None (default), detects available system memory via psutil. For HPC environments with fixed allocations, set this to your SLURM --mem value (e.g., memory_limit_gb=128).
force (bool) – If True, rerun the analysis even when the output h5ad file already exists. If False (default), an existing result is loaded instead of rerun when it was computed from the same input file (path, size, modification time) with the same result-shaping arguments, as recorded in its uns["crispyx"]; otherwise it is recomputed.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) –
How to handle a source stored as CSR. The test streams the matrix by gene (column) chunks, which on CSR re-reads the whole matrix once per chunk and holds it in memory – typically ~100x more I/O than CSC:
- "auto" (default): measure the source’s read throughput and convert only when the re-reads would actually cost more than one conversion does – roughly, a large file on a slow (networked) filesystem streamed in many chunks. A file that re-reads in milliseconds is streamed as-is, silently. See crispyx.data.resolve_auto_format_mismatch_policy().
- "convert": always convert the source to CSC in a temporary file beside output_path (bounded-memory streaming via crispyx.data.convert_to_csc(), honouring memory_limit_gb) and stream from that; the temporary file is removed before returning. Needs ~2x the source file’s size in free disk space there; when that is not available the call falls back to "warn" behaviour with a warning naming the shortfall. Run cx.pp.convert_to_csc once instead if several steps will reuse the file.
- "warn": proceed on the CSR source after one warning that quantifies the cost.
- "off": proceed silently.

Returns:

Differential expression results. Access results via dict-like interface: result[label].logfoldchanges, result[label].auc, result[label].pvalue, etc. The h5ad file path is available at result.result_path.

The h5ad file holds one row per perturbation: the AUC (U / (n_pert * n_ctrl), 0.5 for no shift) in X and the layers z_score, pvalue, pvalue_adj, logfoldchanges, u_statistic and pts. pts_rest describes the control arm, so it is the per-gene column var["pts_rest"]. obs["n_genes_tested"] counts the genes each perturbation was tested on, after the filters.

Return type:

Notes

Differences from scanpy.tl.rank_genes_groups(method="wilcoxon"). Scanpy tests every gene; crispyx applies four low-expression filters by default (min_pct_ctrl=0.01, min_pct_pert=0.002, min_mean_ctrl=0.05, min_mean_pert=0.005) and reports the genes they exclude as untested (NaN), which also keeps them out of the multiple-testing correction. obs["n_genes_tested"] says how many genes each comparison kept. Setting all four to 0 tests Scanpy’s gene set; a gene whose values are all tied is still untested here, where Scanpy reports score=0, pvalue=1.

The rank test itself runs in a numba prange kernel over perturbations and uses every CPU the process is allowed to run on; there is no separate worker-count parameter.

logfoldchanges follows scanpy’s formula, log2((expm1(mean_group) + 1e-9) / (expm1(mean_rest) + 1e-9)). When a gene has no expressing cells in one arm, that arm’s term is the constant, so the magnitude reported is set by the constant – its ceiling is log2(1 / 1e-9) = 29.9 – and not by the data. The p-values are unaffected, but ranking on logfoldchanges puts such genes above every real effect. pts / pts_rest tell them apart: a complete knockdown has pts = 0 beside a high pts_rest, an undetectable gene has both near zero.

crispyx.de.shrink_lfc(data, *, output_dir=None, data_name=None, method='stats', prior_scale_mode='global', min_mu=0.0, n_jobs=-1, batch_size=128, profiling=False, memory_limit_gb=None, verbose=True)[source]

Apply apeGLM log-fold change shrinkage to existing NB-GLM results.

This function applies apeGLM shrinkage using a Cauchy prior on the LFC coefficient. Two methods are available:

stats (default, recommended): Uses pre-computed MLE LFC and SE from the h5ad file with vectorized Newton-Raphson optimization. ~35× faster than full and maintains consistency with stored MLE coefficients.
full: Re-loads original count data and runs per-gene L-BFGS-B optimization. May produce different results from stored MLE for lowly expressed genes due to min_mu clamping differences in the likelihood.

Note

The “stats” method is recommended because CRISPYx NB-GLM fitting uses min_mu=0.5 clamping for numerical stability, which affects the stored MLE coefficients. The “full” method re-evaluates the likelihood without this constraint, potentially finding different optima for lowly expressed genes. The “stats” method preserves shrinkage direction (always toward zero) by working directly with the stored statistics.

This enables separating the base NB-GLM fitting from shrinkage for: - Benchmarking: measure base fitting and shrinkage times separately - Flexibility: apply shrinkage to existing results - Speed: use stats method for production

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file, or a crispyx AnnData / backed anndata.AnnData containing NB-GLM results from nb_glm_test. Must have required layers and metadata in uns.
output_dir (str | Path | None) – Directory for output h5ad file. Defaults to input file’s directory.
data_name (str | None) – Custom name for output file. If None, appends “_shrunk” to input name.
method (Literal['stats', 'full']) –
Shrinkage method to use:
- "stats" (default): Fast vectorized shrinkage using pre-computed MLE statistics. Uses Newton-Raphson optimization across all genes simultaneously. ~35× faster than “full” and maintains consistency with stored MLE coefficients.
- "full": Full model re-fitting with L-BFGS-B per gene. Re-loads original count data. Note: May produce different results for lowly expressed genes due to min_mu clamping differences. Use “stats” for consistent shrinkage behavior.
prior_scale_mode (Literal['global', 'per_comparison']) –
How to estimate the Cauchy prior scale parameter:
- "global" (default): Estimate prior scale once from all perturbations’ MLE LFCs. Faster and often more stable.
- "per_comparison": Estimate prior scale separately for each perturbation vs control comparison. Matches PyDESeq2’s behavior exactly. Use for benchmarking to demonstrate parity.
min_mu (float) – Minimum mean threshold for shrinkage likelihood evaluation. Default: 0.0 (no clamping), matching PyDESeq2’s lfc_shrink which omits min_mu entirely. PyDESeq2 uses min_mu=0.5 for NB-GLM fitting but does NOT pass min_mu to the shrinkage optimizer. This is intentional: the MLE coefficients represent the best fit, and shrinkage should evaluate the same likelihood surface.
n_jobs (int) – Number of parallel jobs for per-gene optimization (only used when method=”full”). Default -1 uses all available cores.
batch_size (int) – Number of genes per joblib batch (only used when method=”full”). Default: 128, matching PyDESeq2.
profiling (bool) – If True, enable timing and memory profiling. When enabled, stores profiling data in adata.uns[“profiling”] with fields: - shrinkage_seconds: Time for lfcShrink operation - shrinkage_peak_memory_mb: Peak memory during shrinkage - profiling_enabled: True When False (default), adata.uns[“profiling”] is set to “NA”.
memory_limit_gb (float | None) – Optional memory budget in gigabytes. With method="full", the number of parallel workers is capped so that they fit in the budget (None, the default, uses the memory currently available). method="stats" holds no per-worker data and ignores it.
verbose (int | bool) – Print basic progress/completion messages. Defaults to True.

Returns:

Updated differential expression results with shrunken LFCs. The result h5ad has: - logfoldchange: shrunken LFC values - logfoldchange_raw: original MLE LFC values (preserved) - standard_error: posterior SE reflecting shrinkage uncertainty - X: updated to the shrunken LFC

Return type:

Examples

>>> # Fast default (recommended for production)
>>> shrunk = crispyx.de.shrink_lfc("nb_glm_result.h5ad")

>>> # Benchmark-accurate (matches PyDESeq2 exactly)
>>> shrunk = crispyx.de.shrink_lfc(
...     "nb_glm_result.h5ad",
...     method="full",
...     prior_scale_mode="per_comparison",
... )

>>> # First run NB-GLM without shrinkage
>>> result = crispyx.de.nb_glm_test(
...     "data.h5ad",
...     perturbation_column="perturbation",
...     lfc_shrinkage_type="none",  # No shrinkage during fitting
... )
>>> # Then apply shrinkage as a separate step
>>> shrunk_result = crispyx.de.shrink_lfc(result.result_path)

Covariate-adjusted effects

residual_effects() fits each gene’s log1p(CP10k) expression on covariates (lane, donor, depth, mitochondrial share, …) over the control cells and reports each perturbation’s mean residual, in one streaming pass over the counts.

Covariate-adjusted perturbation effects: control-fitted linear model residuals.

crispyx.residual.residual_effects(data, perturbation_column, control_label=None, *, covariates, null='model', min_cells=10, min_control_mean=0.1, perturbations=None, gene_name_column=None, chunk_size=None, n_jobs=None, memory_limit_gb=None, format_mismatch_policy='auto', data_name=None, output_path=None, verbose=True, force=False)[source]

Perturbation effects adjusted for covariates by a model fitted on controls.

The method of Feng et al. (2026, Cell Genomics, doi:10.1016/j.xgen.2025.101076). Each gene’s log1p(CP10k) expression is fitted to the covariates by least squares over the control cells only; a perturbation’s effect is the mean, over its cells, of observed minus fitted expression. Differences in donor, lane or batch between a perturbation’s cells and the controls are thereby removed rather than read as an effect.

Everything comes from one pass over the rows of the raw counts. The mean residual of perturbation t is mean_t(y) - mean_t(x) @ B, so only per-perturbation sums and the controls’ cross-products are kept, never a per-cell residual; memory is one perturbations x genes array.

Parameters:

data (str | Path | AnnData | AnnData) – Path to an h5ad file of raw counts, crispyx.AnnData, or an anndata.AnnData opened with backed='r'. The function normalises (log1p of counts per 10,000) itself and refuses input with non-integer or negative values in any chunk it reads.
perturbation_column (str | Series) – Column in adata.obs naming each cell’s perturbation, or a pd.Series of labels indexed by obs_names. Cells whose label is missing (NaN) are not used, which restricts a file to a subset of cells without writing a copy; row chunks without a used cell are not read.
control_label (str | None) – Label of the control cells the model is fitted on. Inferred from common patterns when omitted.
covariates (Sequence[str]) –
Terms of the model besides the intercept; required, with no default. Each is an obs column (one-hot coded, first level dropped, when categorical, string or boolean; used as is when numeric) or one of two values computed from the counts:
- "log_total_counts": natural log of the cell’s total counts;
- "mito_fraction": share of the cell’s counts in genes whose symbol starts with MT-.
Adjust for how cells were processed (donor, lane, batch), which no perturbation changes. Avoid what the perturbation itself can change (depth, mitochondrial share, cell cycle): adjusting for it removes part of the effect.

[] gives the plain difference from the controls. Interactions (lane x donor) and other scores are obs columns made beforehand. Perturbed cells in a categorical level no control shares cannot be adjusted and are left out (counted at verbose=1); numeric covariates must be finite in used cells, and an obs column named like a computed value raises. Collinear covariates are fine unless a perturbation’s effect would then depend on how the fit splits them, which raises.
null (Literal['model', 'robust']) –
Spread the z-score divides by:
- "model" (the paper’s test): sigma_g * sqrt(1/n_t + h_t), with sigma_g the residual standard deviation of the control fit and h_t = mean_t(x) (X^T X)^+ mean_t(x)^T the control fit’s own error at the perturbation’s mean design, which every perturbation shares. h_t is negligible with many controls (the paper’s sigma_g / sqrt(n_t)) and keeps z calibrated with few. Treats cells as independent; cells sharing a guide or donor are not, so it calls too many hits.
- "robust": z = (effect * sqrt(n_t) - m_g) / s_g, with m_g the median over perturbations of effect * sqrt(n_t) and s_g 1.4826 x their median absolute deviation. It assumes most perturbations leave gene g unchanged – true for genome-wide screens, false for a focused panel where every target is active.
h_t is written to obs (leverage) and sigma_g, m_g and s_g to var, so either z can be recomputed.
min_cells (int) – Perturbations with fewer usable cells are left out (listed at verbose=1).
min_control_mean (float) – Genes are tested when the controls’ mean log1p(CP10k) exceeds this (Feng et al. used 0.1). Untested genes keep their effect but have NaN z-scores and p-values.
perturbations (Iterable[str] | None) – Subset of perturbation labels to report. All non-control labels when omitted. The fit uses the controls either way.
gene_name_column (str | None) – Column in adata.var with gene symbols. Uses var_names when omitted.
chunk_size (int | None) – Cells read per chunk. Auto-selected when omitted.
n_jobs (int | None) – Threads for the per-chunk computation (a chunk’s perturbations are split between them); the next chunk is read meanwhile. All cores when omitted. Results do not depend on it.
memory_limit_gb (float | None) – Memory budget. Three perturbations x genes float64 arrays (accumulator, effect, robust spread) and one genes x design-columns array per thread must fit in it, else the call refuses before reading the matrix. Detected from the system when omitted.
format_mismatch_policy (Literal['auto', 'warn', 'convert', 'off']) –
How to handle a source stored as CSC, whose row-(cell-)streaming here is O(total_nnz) per chunk and can be ~100x slower than CSR:
- "auto" (default): measure the source’s read throughput and convert only when the re-reads would cost more than one conversion. See crispyx.data.resolve_auto_format_mismatch_policy().
- "warn": proceed but log a single actionable warning.
- "convert": stream from a temporary CSR copy written beside the output (about the source’s size in free disk), removed before returning.
- "off": proceed silently with no warning.
data_name (str | None) – Custom stem for the output filename.
output_path (str | Path | None) – Exact output h5ad path. Defaults to <stem>_cx_residual_effects.h5ad next to the input.
verbose (int | bool) – 1 / True prints progress and a summary.
force (bool) – Recompute even when an existing output was made by the same call.

Returns:

On-disk result with perturbations in observations and genes in variables. X is the effect: the mean residual, a difference of mean log1p(CP10k) values (natural log) – not a log2 fold change. Layers z_score, pvalue (two-sided normal) and pvalue_adj (Benjamini-Hochberg within each perturbation). obs holds n_cells, n_genes_tested and leverage; var holds control_mean, residual_sd, robust_center, robust_sd and tested; uns records the null, the covariates and the design’s column names.

Return type:

AnnData

GLM fitting

Batched IRLS for negative binomial and Poisson GLMs. NBGLMBatchFitter fits an arbitrary dense design; StructuredGLMBatchFitter exploits the block structure of a [covariates | one-hot groups] design, which is the shape of every perturbation screen and of any design carrying a many-level categorical covariate. fit_nb_glm_batch_auto() chooses between them.

Generalized linear models utilities for differential expression.

class crispyx.glm.NBGLMResult(coef, se, dispersion, converged, n_iter, deviance, max_cooks=None)[source]

Bases: object

Result of fitting a negative binomial GLM for a single gene.

Parameters:

coef (ndarray)
se (ndarray)
dispersion (float)
converged (bool)
n_iter (int)
deviance (float)
max_cooks (float | None)

coef: ndarray

se: ndarray

dispersion: float

converged: bool

n_iter: int

deviance: float

max_cooks: float | None = None

class crispyx.glm.NBGLMBatchResult(coef, se, dispersion, converged, n_iter, deviance)[source]

Bases: object

Result of fitting NB GLM for multiple genes in a batch.

Parameters:

coef (ndarray)
se (ndarray)
dispersion (ndarray)
converged (ndarray)
n_iter (ndarray)
deviance (ndarray)

coef: ndarray

se: ndarray

dispersion: ndarray

converged: ndarray

n_iter: ndarray

deviance: ndarray

class crispyx.glm.ControlStatisticsCache(control_matrix, control_n, control_offset, beta_intercept, control_dispersion, control_xtwx_intercept, control_xtwz_intercept, control_mean_expr, control_expr_counts, pts_rest, global_size_factors=None, global_dispersion=None, global_dispersion_trend=None, global_disp_prior_var=None, dispersion_scope=None, frozen_control_W_sum=None, frozen_control_W_wald_sum=None, frozen_control_Wz_sum=None, frozen_control_mu_sum=None, frozen_control_resid_sq_sum=None, frozen_control_Y_sum=None, use_frozen_control=False)[source]

Bases: object

Cached statistics for control cells to avoid redundant computation.

When fitting independent NB-GLM models for multiple perturbations, each comparison includes the same control cells. This cache precomputes and stores control cell contributions to the IRLS normal equations, allowing them to be reused across all perturbation comparisons instead of being redundantly computed.

The cache stores: - Control cell intercept (β₀): baseline log-expression per gene - Control dispersion: estimated from control cells only - XᵀWX contribution from control cells (for intercept column) - XᵀWz contribution from control cells - Global size factors (optional): precomputed on all cells for consistency - Global dispersion (optional): precomputed MAP dispersion using all cells - Global dispersion prior variance: for MAP shrinkage

Memory optimization: We store control_matrix as dense (not mu/weights) since mu and weights change during IRLS and must be recomputed anyway. Storing the dense matrix avoids repeated .toarray() calls in workers.

This reduces IRLS complexity from O(n_perturbations × n_control × n_genes × n_iter) to O(n_control × n_genes × n_iter) for control-related computations.

Parameters:

control_matrix (ndarray)
control_n (int)
control_offset (ndarray)
beta_intercept (ndarray)
control_dispersion (ndarray)
control_xtwx_intercept (ndarray)
control_xtwz_intercept (ndarray)
control_mean_expr (ndarray)
control_expr_counts (ndarray)
pts_rest (ndarray)
global_size_factors (ndarray | None)
global_dispersion (ndarray | None)
global_dispersion_trend (ndarray | None)
global_disp_prior_var (float | None)
dispersion_scope (str | None)
frozen_control_W_sum (ndarray | None)
frozen_control_W_wald_sum (ndarray | None)
frozen_control_Wz_sum (ndarray | None)
frozen_control_mu_sum (ndarray | None)
frozen_control_resid_sq_sum (ndarray | None)
frozen_control_Y_sum (ndarray | None)
use_frozen_control (bool)

control_matrix: ndarray

control_n: int

control_offset: ndarray

beta_intercept: ndarray

control_dispersion: ndarray

control_xtwx_intercept: ndarray

control_xtwz_intercept: ndarray

control_mean_expr: ndarray

control_expr_counts: ndarray

pts_rest: ndarray

global_size_factors: ndarray | None = None

global_dispersion: ndarray | None = None

global_dispersion_trend: ndarray | None = None

global_disp_prior_var: float | None = None

dispersion_scope: str | None = None

frozen_control_W_sum: ndarray | None = None

frozen_control_W_wald_sum: ndarray | None = None

frozen_control_Wz_sum: ndarray | None = None

frozen_control_mu_sum: ndarray | None = None

frozen_control_resid_sq_sum: ndarray | None = None

frozen_control_Y_sum: ndarray | None = None

use_frozen_control: bool = False

crispyx.glm.precompute_control_statistics(control_matrix, control_offset, *, max_iter=10, tol=1e-06, min_mu=0.5, dispersion_method='moments', global_size_factors=None, freeze_control=False)[source]

Precompute control cell statistics for reuse across perturbation comparisons.

This function fits a simple intercept-only model to control cells to estimate the baseline expression level (β₀) per gene. The resulting intercept, mean expression, weights, and XᵀWX/XᵀWz contributions are cached for reuse.

Parameters:

control_matrix (ndarray | csr_matrix) – Expression matrix for control cells, shape (n_control, n_genes).
control_offset (ndarray) – Log size factors for control cells, shape (n_control,).
max_iter (int) – Maximum IRLS iterations for intercept estimation.
tol (float) – Convergence tolerance.
min_mu (float) – Minimum fitted mean value.
dispersion_method (Literal['moments', 'cox-reid']) – Method for dispersion estimation.
global_size_factors (ndarray | None) – Optional precomputed size factors for all cells (n_cells_total,). When provided, stored in cache for use across all comparisons.
freeze_control (bool) – If True, compute frozen sufficient statistics (W_sum, Wz_sum, etc.) and set control_matrix to None to save memory. This reduces per-worker pickle size from ~5GB to ~1MB for large datasets. Workers can use the frozen stats with fit_batch_with_frozen_control() instead of the raw matrix. Default False.

Returns:

Cached statistics for control cells.

Return type:

crispyx.glm.precompute_control_statistics_streaming(path, control_mask, control_offset, *, max_iter=10, tol=1e-06, min_mu=0.5, global_size_factors=None, freeze_control=True, chunk_size=4096)[source]

Streaming version of precompute_control_statistics for very large control groups.

Instead of densifying the full control matrix (which can exceed 100+ GiB for large datasets), this function reads control cells from disk in chunks and accumulates sufficient statistics for the intercept-only IRLS model.

Peak memory is O(chunk_size × n_genes) instead of O(n_control × n_genes).

Parameters:

path (str | Path) – Path to the h5ad file containing the count matrix.
control_mask (np.ndarray) – Boolean mask over all cells indicating control cells, shape (n_cells,).
control_offset (np.ndarray) – Log size factors for control cells, shape (n_control,).
max_iter (int) – Maximum IRLS iterations for intercept estimation.
tol (float) – Convergence tolerance.
min_mu (float) – Minimum fitted mean value.
global_size_factors (np.ndarray | None) – Optional precomputed size factors for all cells (n_cells_total,).
freeze_control (bool) – Must be True for streaming mode (raw control_matrix is never materialised).
chunk_size (int) – Number of control cells to process per chunk. Default 4096.

Returns:

Cached statistics with frozen control sufficient statistics.

Return type:

crispyx.glm.precompute_global_dispersion(control_cache, all_cell_matrix, all_cell_offset, *, n_grid=25, min_disp=1e-08, max_disp=10.0, fit_type='parametric', fast_mode=True, max_dense_fraction=0.3, memory_limit_gb=None)[source]

Precompute global dispersion trend using all cells and update cache.

This function computes a global dispersion trend using all cells in the dataset (control + all perturbations), similar to how DESeq2/PyDESeq2 estimates dispersion from all samples. The trend is then used for MAP shrinkage in all per-perturbation comparisons.

Using global dispersion has several advantages: 1. More stable estimates (larger sample size for trend fitting) 2. ~10× faster per-perturbation fitting (skips MAP estimation) 3. More consistent results across perturbations

Memory-adaptive behavior: If the estimated memory for densifying the matrix exceeds max_dense_fraction × min(available_memory, memory_limit_gb), the function switches to chunk-wise streaming processing.

Parameters:

control_cache (ControlStatisticsCache) – Precomputed control cell statistics cache.
all_cell_matrix (ndarray | csr_matrix) – Count matrix for all cells, shape (n_cells, n_genes).
all_cell_offset (ndarray) – Log size factors for all cells, shape (n_cells,).
n_grid (int) – Number of grid points for dispersion MAP estimation.
min_disp (float) – Minimum dispersion value.
max_disp (float) – Maximum dispersion value.
fit_type (Literal['parametric', 'local', 'mean']) – Type of trend fitting (“parametric”, “local”, or “mean”).
fast_mode (bool) – If True, use simple trend shrinkage (MoM → trend) instead of expensive MAP estimation. This is ~50× faster and suitable for most use cases. If False, use full MAP estimation with grid search + refinement.
max_dense_fraction (float) – Maximum fraction of available memory to use for dense matrix. If estimated memory exceeds this, switch to streaming mode. Default is 0.3 (30% of available memory).
memory_limit_gb (float | None) – Optional explicit memory limit in GB. If provided, the effective memory budget is min(available_memory, memory_limit_gb).

Returns:

Updated cache with global_dispersion, global_dispersion_trend, and global_disp_prior_var fields populated.

Return type:

crispyx.glm.precompute_global_dispersion_from_path(path, control_cache, all_cell_offset, *, min_disp=1e-08, max_disp=10.0, fit_type='parametric', chunk_size=4096)[source]

Path-based streaming global dispersion for very large datasets.

This function estimates dispersion by streaming directly from an h5ad file, avoiding the need to load the entire matrix into memory. Uses method-of-moments estimation accumulated across chunks read from disk.

This is the preferred method for datasets that exceed available memory (e.g., Replogle-GW-k562 with ~2M cells × 8K genes = ~131 GB).

Parameters:

path (str | Path) – Path to the h5ad file containing the count matrix.
control_cache (ControlStatisticsCache) – Precomputed control cell statistics cache.
all_cell_offset (np.ndarray) – Log size factors for all cells, shape (n_cells,).
min_disp (float) – Minimum dispersion value.
max_disp (float) – Maximum dispersion value.
fit_type (Literal['parametric', 'local', 'mean']) – Type of trend fitting (“parametric”, “local”, or “mean”).
chunk_size (int) – Number of cells to process per chunk. Larger chunks are faster but use more memory. Default 4096 balances speed and memory.

Returns:

Updated cache with global_dispersion, global_dispersion_trend, and global_disp_prior_var fields populated.

Return type:

class crispyx.glm.NBGLMFitter(design, *, offset=None, dispersion=None, max_iter=50, tol=1e-06, poisson_init_iter=20, ridge_penalty=1e-06, min_mu=0.0, min_total_count=1.0, compute_cooks=False, dispersion_method='cox-reid')[source]

Bases: object

L-BFGS-B solver for negative binomial GLMs.

Parameters:

design (ArrayLike) – The design matrix with shape (n_samples, n_features).
offset (ArrayLike | None) – Optional log-scale offset (e.g. library size) per sample.
dispersion (float | None) – Optional dispersion (alpha) for the negative binomial. If None the dispersion is re-estimated at each iteration using a method-of-moments update similar to the one used in statsmodels.
max_iter (int) – Maximum number of outer iterations for alternating optimization.
tol (float) – Absolute tolerance on the coefficient updates for convergence.
poisson_init_iter (int) – Maximum number of iterations for the Poisson initialisation stage. If set to 0 the Poisson warm start is skipped and coefficients are initialised at zero.
ridge_penalty (float) – Small diagonal ridge penalty added to the weighted normal equations to improve numerical stability. This does not change the estimator when the system is well conditioned but prevents failures when the Hessian is nearly singular.
min_mu (float) – Lower bound on the fitted mean to avoid issues with extremely small predicted counts for lowly expressed genes.
min_total_count (float) – Genes with a total count below this threshold are not fitted (the resulting NBGLMResult will report converged=False).
dispersion_method (Literal['moments', 'cox-reid']) –
Method for estimating dispersion when dispersion is None.
- "moments": Method-of-moments (fast but less accurate).
- "cox-reid": Cox-Reid adjusted profile likelihood (slower but more accurate, similar to DESeq2).
compute_cooks (bool)

fit_gene(counts)[source]

Fit NB GLM for a single gene using L-BFGS-B optimization.

Parameters:: counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes])
Return type:: NBGLMResult

fit_matrix(matrix, *, batch_size=None)[source]

Fit NB GLM for every gene (column) in a count matrix.

Parameters:

matrix (array-like of shape (n_samples, n_genes)) – Raw count matrix. Sparse (CSC) and dense formats are accepted.
batch_size (int or None, optional) – Number of genes to densify at once when matrix is sparse. None processes all genes in one batch.

Returns:

One result per gene, in column order.

Return type:

list of NBGLMResult

estimate_dispersion_cox_reid(y, mu, weights, *, initial_alpha=0.1, bounds=(1e-08, 1000.0))[source]

Estimate dispersion using Cox-Reid adjusted profile likelihood.

This method maximizes the adjusted profile log-likelihood for the dispersion parameter, which includes a bias correction term based on the Cox-Reid adjustment. This approach is similar to DESeq2’s dispersion estimation.

Parameters:

y (ndarray) – Count vector of shape (n_samples,).
mu (ndarray) – Current fitted mean values of shape (n_samples,).
weights (ndarray) – IRLS weights from the current fit.
initial_alpha (float) – Starting value for the optimization.
bounds (tuple[float, float]) – Lower and upper bounds for the dispersion parameter.

Returns:

Estimated dispersion parameter.

Return type:

float

crispyx.glm.build_design_matrix(obs_frame, *, covariate_columns, perturbation_indicator, intercept=True)[source]

Construct a design matrix from covariates and a perturbation indicator.

Parameters:

obs_frame – Pandas DataFrame (preferred) or structured numpy array containing covariate columns.
covariate_columns (Sequence[str]) – Columns that should be included in the design matrix. Categorical columns are expanded using one-hot encoding (dropping the first level).
perturbation_indicator (ndarray) – Binary array of shape (n_samples,) marking perturbed cells.
intercept (bool) – Whether to prepend an intercept column to the design matrix.

Returns:

design – The numeric design matrix as a numpy.ndarray of float64.
column_names – The column names corresponding to the design matrix.

Return type:

tuple[ndarray, list[str]]

crispyx.glm.fit_dispersion_trend(means, dispersions, *, min_mean=0.5, fit_type='parametric', n_iter=10)[source]

Fit a smooth mean-dispersion trend using DESeq2/PyDESeq2-style Gamma GLM.

The parametric trend models dispersion as dispersion = asymptDisp + extraPois / mean.

This matches PyDESeq2’s fitDispersionTrend which uses iteratively reweighted least squares (IRLS) with a Gamma family and log link, fitting E[disp] = a0 + a1 / mean.

Outliers are iteratively removed based on prediction ratio bounds.

Parameters:

means (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Mean expression values per gene (normalized counts).
dispersions (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Raw genewise dispersion estimates.
min_mean (float) – Minimum mean value for fitting (genes below this are excluded).
fit_type (Literal['parametric', 'local', 'mean']) – Type of trend fitting: - “parametric”: DESeq2/PyDESeq2-style Gamma GLM (recommended) - “local”: Weighted local regression (LOWESS-like) - “mean”: Simple median (fallback)
n_iter (int) – Number of iterations for IRLS fitting.

Returns:

Fitted trend values for each gene.

Return type:

np.ndarray

crispyx.glm.shrink_dispersions(raw, trend, *, prior_df=None, min_prior_df=1.0, max_prior_df=100.0, outlier_sigma=2.0, n_iter=5)[source]

Shrink dispersions toward fitted trend using empirical Bayes.

This implements a DESeq2/PyDESeq2-style empirical Bayes shrinkage where the prior variance is estimated from the distribution of log-dispersion residuals around the trend using an iterative trimmed variance estimator.

Genes with dispersions more than outlier_sigma standard deviations above the trend keep their MLE (not shrunken), matching PyDESeq2’s outlier handling.

The shrinkage formula is:: log(shrunk) = (prior_df * log(trend) + log(raw)) / (prior_df + 1)

This is equivalent to a posterior mean estimate under a log-normal prior.

Parameters:

raw (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Raw MLE dispersion estimates.
trend (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Fitted mean-dispersion trend values.
prior_df (float | None) – Prior degrees of freedom controlling shrinkage strength. If None, estimated empirically from the data using iterative trimming.
min_prior_df (float) – Minimum allowed prior degrees of freedom.
max_prior_df (float) – Maximum allowed prior degrees of freedom.
outlier_sigma (float) – Number of standard deviations above trend beyond which genes keep their MLE dispersion (not shrunken). Set to inf to disable.
n_iter (int)

Returns:

Shrunken dispersion estimates.

Return type:

np.ndarray

crispyx.glm.estimate_dispersion_map(Y, mu, trend, *, prior_var=None, min_disp=1e-08, max_disp=10.0, n_grid=25, refine=True, n_jobs=-1)[source]

Estimate MAP dispersion using vectorized grid search + optional refinement.

This implements PyDESeq2-style MAP estimation where the dispersion is estimated by maximizing log L(Y | mu, alpha) + log prior(alpha | trend, prior_var).

The prior is log-normal: log(alpha) ~ N(log(trend), prior_var).

Optimization (v5): Uses fused Numba kernel that combines grid search with Brent’s method refinement in a single parallel pass over genes. This eliminates joblib process spawning overhead, achieving ~2-3× speedup over the previous joblib-based refinement while maintaining identical accuracy to scipy’s minimize_scalar.

Default is n_grid=25, refine=True which provides optimal balance of speed and accuracy. The Brent refinement makes grid size largely irrelevant for accuracy (all grid sizes 15-50 achieve perfect correlation).

Parameters:

Y (ndarray) – Count matrix of shape (n_cells, n_genes).
mu (ndarray) – Fitted mean matrix of shape (n_cells, n_genes).
trend (ndarray) – Dispersion trend values of shape (n_genes,).
prior_var (float | None) – Variance of the log-normal prior. If None, estimated from data using the variance of log-dispersion residuals around trend.
min_disp (float) – Minimum allowed dispersion value.
max_disp (float) – Maximum allowed dispersion value.
n_grid (int) – Number of grid points for initial search. More points = better initial estimate but slower grid search. Default is 50 for good accuracy (96% Top-100 overlap with PyDESeq2 without refinement).
refine (bool) – If True, refine the grid search result using Brent’s method. Default is True for best accuracy (99% Top-100 overlap with PyDESeq2). Set to False for ~2× speedup if slight accuracy loss is acceptable.
n_jobs (int) – Number of parallel jobs for refinement. -1 uses all cores.

Returns:

MAP dispersion estimates of shape (n_genes,).

Return type:

np.ndarray

crispyx.glm.estimate_glm_map_dispersion(design, counts, *, offset, size_factors, max_iter=25, tol=1e-06, poisson_init_iter=5, dispersion_method='cox-reid', min_mu=0.0, min_total_count=1.0)[source]

Per-gene MAP dispersion of an NB GLM, and the trend it is shrunk toward.

The dispersion step of DESeq2: fit the GLM to every gene (with dispersion_method estimates), fit a mean-dispersion trend to those estimates, then take each gene’s maximum a posteriori dispersion at the fitted means under a log-normal prior centred on the trend. A gene whose fit did not converge is evaluated at its normalised counts, and one with no finite estimate takes the trend value.

Parameters:

design (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Design matrix of shape (n_samples, n_features).
counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
offset (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Log size factors, shape (n_samples,).
size_factors (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Size factors, shape (n_samples,), for the normalised means.
max_iter (int) – As for fit_nb_glm_batch_auto().
tol (float) – As for fit_nb_glm_batch_auto().
poisson_init_iter (int) – As for fit_nb_glm_batch_auto().
dispersion_method (Literal['moments', 'cox-reid']) – As for fit_nb_glm_batch_auto().
min_mu (float) – As for fit_nb_glm_batch_auto().
min_total_count (float) – As for fit_nb_glm_batch_auto().

Returns:

Arrays of shape (n_genes,).

Return type:

dispersion, trend

crispyx.glm.shrink_lfc_apeglm(counts, design_matrix, size_factors, dispersion, mle_coef, mle_se, *, shrink_index=1, prior_scale=None, prior_no_shrink_scale=15.0, max_iter=100, tol=1e-06, n_jobs=-1, batch_size=128, min_mu=0.0)[source]

Apply apeGLM LFC shrinkage using Cauchy prior (PyDESeq2-compatible).

This implements the apeGLM (approximate posterior estimation for GLM) shrinkage method used by DESeq2/PyDESeq2. The method re-fits the NB-GLM model with a Cauchy prior penalty on the LFC coefficient.

Key features matching PyDESeq2: - L-BFGS-B optimization with analytical gradients - Grid search fallback over [-5, 5] with 50 points when optimization fails - Parallel per-gene optimization using joblib with batch_size=128 - Proper NB likelihood formulation with numerical stability - min_mu clamping for consistency with NB-GLM fitting (if min_mu > 0)

Parameters:

counts (ndarray) – Raw count matrix (n_cells, n_genes).
design_matrix (ndarray) – Design matrix (n_cells, n_params).
size_factors (ndarray) – Size factors for each cell (n_cells,).
dispersion (ndarray) – Gene-wise dispersion estimates (n_genes,).
mle_coef (ndarray) – MLE coefficient matrix (n_params, n_genes) for warm-starting.
mle_se (ndarray) – Standard errors of MLE LFC estimates (n_genes,).
shrink_index (int) – Index of the coefficient to shrink (default: 1, the LFC coefficient).
prior_scale (float | None) – Scale parameter for Cauchy prior. If None, estimated globally from the MLE LFC distribution (matching PyDESeq2’s approach).
prior_no_shrink_scale (float) – Scale for normal prior on non-shrunk coefficients (default: 15.0).
max_iter (int) – Maximum iterations for L-BFGS-B optimization.
tol (float) – Convergence tolerance.
n_jobs (int) – Number of parallel jobs. Default -1 uses all available cores.
batch_size (int) – Number of genes per joblib batch (default: 128, matching PyDESeq2).
min_mu (float) – Minimum mean value for mu clamping in NB log-likelihood. If > 0, mu is clamped to be at least min_mu. This should match the min_mu used during NB-GLM fitting to ensure the stored coefficients are consistent with the likelihood surface (default: 0.0 = no clamping).

Returns:

shrunk_coef – Shrunken coefficient matrix (n_params, n_genes).
shrunk_se – Approximate standard errors from inverse Hessian (n_genes,).
converged – Boolean array indicating convergence for each gene (n_genes,).

Return type:

tuple[ndarray, ndarray, ndarray]

crispyx.glm.shrink_lfc_apeglm_from_stats(mle_lfc, mle_se, xtwx_diag=None, *, base_mean=None, prior_scale=None, max_iter=100, tol=1e-08, use_gene_specific_prior=True, hybrid_fallback=True, hybrid_mle_se_threshold=3.0, hybrid_base_mean_threshold=10.0)[source]

Apply apeGLM-style shrinkage using pre-computed MLE statistics.

This is a memory-efficient, fully vectorized version of apeGLM shrinkage that uses pre-computed MLE coefficients and standard errors without requiring access to the full count matrix. It applies a Cauchy prior and finds the MAP estimate using a damped Newton-Raphson optimization across all genes simultaneously.

Accuracy improvements (v2):

Gene-specific prior scales based on expression level (sqrt(base_mean)).
Moment-corrected SE using observed Fisher information approximation.
Hybrid fallback: marks genes with abs(MLE)/SE > threshold or low expression for full NB-GLM re-fitting to achieve Eff rho >= 0.98.

The posterior mode is found by solving beta_MAP = argmin { (beta - mle_lfc)^2 / (2 * se^2) + log(1 + (beta/s)^2) } which is the negative log posterior with a Cauchy prior.

This implementation uses a robust damped Newton method with Hessian clamping to ensure convergence even when the standard Hessian becomes negative (which can happen when abs(beta) > prior_scale). This matches PyDESeq2’s behavior which uses L-BFGS-B for robustness.

Parameters:

mle_lfc (ndarray) – MLE log-fold change estimates (n_genes,).
mle_se (ndarray) – Standard errors of MLE estimates (n_genes,).
xtwx_diag (ndarray | None) – Diagonal elements of X’WX for Fisher information (n_genes,). If provided, used for posterior SE computation. If None, uses approximation based on MLE SE.
base_mean (ndarray | None) – Mean normalized expression per gene (n_genes,). Used for gene-specific prior scaling. If None, uniform prior scale is used.
prior_scale (float | None) – Base scale parameter for Cauchy prior. If None, estimated adaptively. When use_gene_specific_prior=True, this is scaled per-gene.
max_iter (int) – Maximum iterations for optimization.
tol (float) – Convergence tolerance.
use_gene_specific_prior (bool) – If True, scale prior_scale by 1/sqrt(base_mean) per gene. This accounts for the fact that lowly-expressed genes have higher variance and should be shrunk more aggressively. Default True.
hybrid_fallback (bool) – If True, return a mask indicating genes that need full NB-GLM re-fitting due to problematic statistics (large abs(MLE)/SE or low expression).
hybrid_mle_se_threshold (float) – Genes with abs(MLE)/SE > this threshold are marked for full re-fitting.
hybrid_base_mean_threshold (float) – Genes with base_mean < this threshold are marked for full re-fitting.

Returns:

shrunk_lfc – Shrunken log-fold change estimates (n_genes,).
shrunk_se – Posterior standard errors (n_genes,).
converged – Boolean array indicating convergence for each gene (n_genes,).
needs_full_refit – Boolean array indicating genes that need full NB-GLM re-fitting (n_genes,). Only populated when hybrid_fallback=True, otherwise all False.

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

crispyx.glm.compute_cooks_distance_batch(Y, mu, dispersion, n_params=2)[source]

Compute Cook’s distance for each observation in a batch of genes.

Cook’s distance measures the influence of each observation on the fitted model. Large values indicate potential outliers that disproportionately affect the estimates.

For GLMs, Cook’s distance is computed as:: D_i = (r_i^2 / p) * (h_ii / (1 - h_ii)^2)

where r_i is the Pearson residual, p is the number of parameters, and h_ii is the leverage (diagonal of the hat matrix).

Parameters:

Y (ndarray) – Count matrix of shape (n_cells, n_genes).
mu (ndarray) – Fitted mean matrix of shape (n_cells, n_genes).
dispersion (ndarray) – Dispersion estimates of shape (n_genes,).
n_params (int) – Number of model parameters (default 2 for intercept + treatment).

Returns:

Cook’s distance matrix of shape (n_cells, n_genes).

Return type:

np.ndarray

crispyx.glm.filter_outliers_cooks(Y, mu, dispersion, *, n_params=2, threshold_quantile=0.99)[source]

Identify and replace outlier counts based on Cook’s distance.

Following DESeq2’s approach: 1. Compute Cook’s distance for each observation 2. Identify outliers where Cook’s D > F(threshold_quantile, p, n-p) 3. Replace outlier counts with trimmed mean from non-outlier samples

Parameters:

Y (ndarray) – Count matrix of shape (n_cells, n_genes).
mu (ndarray) – Fitted mean matrix of shape (n_cells, n_genes).
dispersion (ndarray) – Dispersion estimates of shape (n_genes,).
n_params (int) – Number of model parameters.
threshold_quantile (float) – Quantile of F distribution for outlier threshold.

Returns:

Y_filtered: Count matrix with outliers replaced
outlier_mask: Boolean matrix indicating outliers (n_cells, n_genes)

Return type:

Tuple[np.ndarray, np.ndarray]

crispyx.glm.estimate_covariate_effects_streaming(backed_adata, *, obs_df, perturbation_labels, control_label, covariate_columns, size_factors, chunk_size=2048, poisson_iter=10, tol=1e-06, return_intercept=False)[source]

Estimate global intercept and covariate effects using control cells only.

This function fits a Poisson regression using only control cells to estimate: - An intercept (baseline expression for control cells) - Covariate effects (if any covariates are specified)

By using only control cells, the intercept represents the true control baseline that can then be used as an offset in per-perturbation fitting. This ensures that perturbation effects are properly estimated as deviations from control.

For the intercept-only case (no covariates), the closed-form MLE is used:: intercept = log(sum(Y) / sum(size_factors))

When covariates are present, IRLS is used with proper per-gene weighting.

Parameters:

backed_adata – Backed AnnData object opened in read mode.
obs_df (pd.DataFrame) – Full obs DataFrame with all cells (already loaded).
perturbation_labels (np.ndarray) – Array of perturbation labels for all cells.
control_label (str) – The label identifying control cells.
covariate_columns (Sequence[str]) – List of covariate column names to include.
size_factors (np.ndarray) – Per-cell size factors (length n_cells).
chunk_size (int) – Number of cells to process per chunk.
poisson_iter (int) – Number of Poisson IRLS iterations.
tol (float) – Convergence tolerance.
return_intercept (bool) – If True, also return the global intercept coefficients. This is useful for joint fitting where the intercept should be shared across all perturbation comparisons.

Returns:

If return_intercept is False:: Covariate effects of shape (n_covariates, n_genes). These are the log-scale regression coefficients for each covariate.
If return_intercept is True:: Tuple of (covariate_effects, intercept) where: - covariate_effects has shape (n_covariates, n_genes) - intercept has shape (n_genes,) representing the control baseline

Return type:

np.ndarray or tuple

crispyx.glm.estimate_global_dispersion_streaming(backed_adata, *, obs_df, perturbation_labels, control_label, covariate_columns, size_factors, beta_intercept, beta_cov=None, beta_perturbation=None, chunk_size=2048, dispersion_method='cox-reid', poisson_iter=10, tol=1e-06)[source]

Estimate global per-gene dispersion using all cells via streaming.

This function streams through all cells to estimate dispersion for each gene using a full design matrix. The dispersion is estimated using all conditions together, which provides more stable estimates than per-perturbation estimation, similar to how PyDESeq2 estimates dispersion from all samples.

The function uses the pre-estimated intercept, perturbation effects, and covariate effects to compute fitted values (mu), then estimates dispersion from the residuals using method-of-moments.

Parameters:

backed_adata – Backed AnnData object opened in read mode.
obs_df (pd.DataFrame) – Full obs DataFrame with all cells (already loaded).
perturbation_labels (np.ndarray) – Array of perturbation labels for all cells.
control_label (str) – The label identifying control cells (used as reference level).
covariate_columns (Sequence[str]) – List of covariate column names to include.
size_factors (np.ndarray) – Per-cell size factors (length n_cells).
beta_intercept (np.ndarray) – Pre-estimated global intercept coefficients, shape (n_genes,).
beta_cov (np.ndarray | None) – Pre-estimated covariate effects, shape (n_covariates, n_genes). If None or empty, no covariate adjustment is applied.
beta_perturbation (np.ndarray | None) – Pre-estimated perturbation effects, shape (n_perturbations, n_genes). If None, perturbation effects are estimated via Poisson IRLS.
chunk_size (int) – Number of cells to process per chunk.
dispersion_method (Literal['moments', 'cox-reid']) – Method for dispersion estimation: - “moments”: Method-of-moments (fast but less accurate) - “cox-reid”: Cox-Reid adjusted profile likelihood (more accurate)
poisson_iter (int) – Number of Poisson IRLS iterations for refining mu estimates. Only used if beta_perturbation is None.
tol (float) – Convergence tolerance.

Returns:

Dispersion estimates of shape (n_genes,). These are the alpha values for the negative binomial distribution: Var(Y) = mu + alpha * mu^2.

Return type:

np.ndarray

class crispyx.glm.NBGLMBatchFitter(design, *, offset=None, max_iter=25, tol=1e-06, poisson_init_iter=5, dispersion_method='cox-reid', min_mu=0.0, min_total_count=1.0, ridge_penalty=1e-06, family='nb')[source]

Bases: object

Vectorized batch fitter for NB GLM across multiple genes.

This fitter processes all genes simultaneously using vectorized operations, providing significant speedup compared to per-gene fitting. It uses IRLS (Iteratively Reweighted Least Squares) with batched matrix operations.

Parameters:

design (ArrayLike) – Design matrix with shape (n_samples, n_features).
offset (ArrayLike | None) – Log-scale offset (e.g., log size factors) per sample.
max_iter (int) – Maximum IRLS iterations.
tol (float) – Convergence tolerance on coefficient updates.
poisson_init_iter (int) – Initial Poisson iterations for warm start.
dispersion_method (Literal['moments', 'cox-reid']) – Method for dispersion estimation: “moments” or “cox-reid”.
min_mu (float) – Minimum fitted mean to avoid numerical issues.
min_total_count (float) – Minimum total count for a gene to be fitted.
family (Literal['nb', 'poisson']) – "nb" (default) fits a negative binomial and estimates the dispersion; "poisson" fits a genuine Poisson GLM, whose variance is the mean. A Poisson fit is not an NB fit with a small estimated dispersion: the weights are mu rather than mu / (1 + alpha mu) and no dispersion is estimated. dispersion is reported as 0.
ridge_penalty (float)

fit_batch(counts, gene_batch_size='auto', use_numba=True, fixed_dispersion=None)[source]

Fit NB GLM for all genes in the count matrix.

Memory-optimized implementation with optional gene batching and Numba acceleration for the 2-parameter model (intercept + perturbation).

Parameters:

counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
gene_batch_size (int | Literal['auto'] | None) – Number of genes to process per batch. If “auto”, calculated based on memory constraints (~100 MB per batch). If None, process all genes at once (legacy behavior).
use_numba (bool) – Whether to use Numba-accelerated IRLS for 2-feature models. Default True for better memory efficiency.
fixed_dispersion (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – Per-gene NB dispersion of shape (n_genes,) to hold fixed instead of estimating it. Ignored when family="poisson".

Returns:

Results for all genes with vectorized arrays.

Return type:

fit_batch_with_covariate_offset(counts, covariate_offset)[source]

Fit NB GLM with pre-computed covariate offset.

This method is used in the joint fitting approach where covariate effects are estimated globally and then held fixed during per-perturbation fitting. The covariate offset is subtracted from the working response during IRLS.

Parameters:

counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
covariate_offset (ndarray) – Pre-computed covariate offset of shape (n_samples, n_genes), representing X_cov @ beta_cov for the covariate portion of the design.

Returns:

Results for all genes with vectorized arrays.

Return type:

fit_batch_with_joint_offsets(counts, *, intercept_offset=None, covariate_offset=None, fixed_dispersion=None)[source]

Fit NB GLM with pre-computed intercept and covariate offsets.

This method is used in the joint fitting approach where the global intercept (and optionally covariates and dispersion) are estimated using all cells and then held fixed during per-perturbation fitting.

The design matrix should NOT include an intercept column when using intercept_offset, as the global intercept is added as an offset.

Parameters:

counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
intercept_offset (ndarray | None) – Pre-computed global intercept of shape (n_genes,). If provided, the design matrix should not include an intercept column.
covariate_offset (ndarray | None) – Pre-computed covariate offset of shape (n_samples, n_genes), representing X_cov @ beta_cov for the covariate portion.
fixed_dispersion (ndarray | None) – If provided, use these dispersion values instead of estimating. Shape (n_genes,).

Returns:

Results for all genes with vectorized arrays.

Return type:

fit_batch_with_control_cache(perturbation_matrix, perturbation_offset, control_cache, *, perturbation_indicator, valid_mask=None)[source]

Fit NB GLM using precomputed control cell statistics.

This method provides significant speedup by reusing control cell contributions (XᵀWX, XᵀWz) from a precomputed cache instead of redundantly computing them for each perturbation comparison.

The design matrix is [1, perturbation_indicator] where: - Control cells have indicator = 0 - Perturbation cells have indicator = 1

The control contribution to XᵀWX and XᵀWz is taken from the cache, and only perturbation cell contributions are computed fresh.

Parameters:

perturbation_matrix (ndarray | csr_matrix) – Expression matrix for perturbation cells only, shape (n_pert, n_genes).
perturbation_offset (ndarray) – Log size factors for perturbation cells, shape (n_pert,).
control_cache (ControlStatisticsCache) – Precomputed control cell statistics from precompute_control_statistics.
perturbation_indicator (ndarray) – Binary indicator for perturbation cells in the combined design. Should be shape (n_control + n_pert,) with 0 for control, 1 for perturbation.
valid_mask (ndarray | None) – Optional boolean mask for genes to fit, shape (n_genes,).

Returns:

Fitted coefficients and statistics.

Return type:

fit_batch_with_frozen_control(perturbation_matrix, perturbation_offset, control_cache, *, valid_mask=None)[source]

Fit NB GLM using frozen control sufficient statistics (memory-efficient).

This method uses precomputed sufficient statistics from control cells instead of the raw control_matrix, reducing per-worker memory from ~5GB to ~1MB for large datasets.

Key differences from fit_batch_with_control_cache: - β₀ (intercept) is FROZEN to the value estimated from control cells - Only β₁ (perturbation effect) is estimated - Control contributions (W_sum, Wz_sum) are pre-computed constants - No access to raw control_matrix (control_cache.control_matrix is None)

The design matrix is [1, perturbation_indicator] where: - Control cells have indicator = 0 (contributions are frozen) - Perturbation cells have indicator = 1

Parameters:

perturbation_matrix (ndarray | csr_matrix) – Expression matrix for perturbation cells only, shape (n_pert, n_genes).
perturbation_offset (ndarray) – Log size factors for perturbation cells, shape (n_pert,).
control_cache (ControlStatisticsCache) – Precomputed control cell statistics with use_frozen_control=True. Must have frozen_control_W_sum and frozen_control_Wz_sum set.
valid_mask (ndarray | None) – Optional boolean mask for genes to fit, shape (n_genes,).

Returns:

Fitted coefficients and statistics.

Return type:

Notes

This method is designed for parallel processing where each worker handles a subset of perturbations. By using frozen control stats:

Per-worker pickle size: ~5GB → ~1MB (control_matrix not needed)
Memory enables: 2 workers → 32 workers (for 128GB memory limit)
Time reduction: ~300h → ~10h (for genome-wide screens)

Mathematical justification: With global dispersion and fixed β₀, the control cells’ contribution to XᵀWX and XᵀWz is constant across all perturbation comparisons:

μ_control = exp(β₀ + offset) is fixed (no perturbation indicator)
W_control = μ²/(μ + α*μ²) depends only on μ and global α
z_control = η + (Y - μ)/μ - offset = β₀ + (Y - μ)/μ

Therefore, sum(W_control) and sum(W_control * z_centered) are constants that can be pre-computed once and reused across all comparisons.

crispyx.glm.detect_onehot_block(design, min_block=2)[source]

Find a maximal block of disjoint one-hot columns in a design matrix.

A column qualifies if it is binary and its support does not overlap the support of any column already selected. Columns are scanned from the smallest support upwards, because a wide binary covariate (sex, say, covering about half the cells) would otherwise be selected first and block every narrow group column it overlaps.

This exists for callers that hold a design matrix in which the structure has already been flattened. A caller that still knows which columns are the group indicators should pass them to StructuredGLMBatchFitter directly rather than rediscover them.

Parameters:

design (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Design matrix of shape (n_samples, n_features).
min_block (int) – Minimum number of qualifying columns; below this an empty array is returned.

Returns:

Indices of the block, sorted in the original column order, or an empty array when fewer than min_block columns qualify.

Return type:

numpy.ndarray

class crispyx.glm.StructuredGLMResult(coef, se, dispersion, deviance, converged, n_iter, mu=None, dev_resid=None)[source]

Bases: object

Result of fitting a batch of GLMs on a [covariates | groups] design.

coef and se are laid out with the covariate columns first, in the order of design, followed by the group columns in the order of groups.

Parameters:

coef (ndarray)
se (ndarray)
dispersion (ndarray)
deviance (ndarray)
converged (ndarray)
n_iter (ndarray)
mu (ndarray | None)
dev_resid (ndarray | None)

coef: ndarray

se: ndarray

dispersion: ndarray

deviance: ndarray

converged: ndarray

n_iter: ndarray

mu: ndarray | None = None

dev_resid: ndarray | None = None

class crispyx.glm.StructuredGLMBatchFitter(design, groups, *, offset=None, family='nb', max_iter=50, tol=1e-08, ridge=1e-06, ridge_group=0.0001, clip_group=10.0, min_mu=0.0, dispersion_method='cox-reid')[source]

Bases: object

Batched IRLS for designs [X | G] with G one-hot.

Every perturbation screen has this shape: a few dense covariates (an intercept, library size, continuous covariates) beside many indicator columns, one per perturbation or per batch/donor level, at most one of which is set for any cell. Because the indicator supports are disjoint, G' W G is diagonal and each Newton step can be taken through the Schur complement of that diagonal block, at O(n p d_X^2 + p a d_X^2) instead of the O(n p (d_X + a)^2) a dense solver pays. The result is the exact IRLS solution of the same GLM, not an approximation.

When to prefer this over NBGLMBatchFitter: the advantage grows with a / d_X and is not automatic. Measured on one Newton step at n=3,000, p=8,563: 20.9x at d_X=31, a=200; 9.4x at d_X=12, a=100; 4.4x at d_X=2, a=50; but only 1.5x at d_X=12, a=20. The crossover is near a ~ d_X. Memory is the other argument: the dense per-gene Hessian at d=231, p=8,563 is 3.66 GB against 0.42 GB for the structured blocks.

Parameters:

design (ArrayLike) – Dense covariates X, shape (n_samples, n_features). Include the intercept here.
groups – One-hot indicators G, shape (n_samples, n_groups), dense or sparse. Rows that are all zero are cells belonging to no group.
offset (ArrayLike | None) – Log-scale offset (e.g. log size factors) per cell.
family (Literal['nb', 'poisson']) – "nb" or "poisson".
max_iter (int) – IRLS iteration cap and convergence tolerance; a gene is converged when both its relative deviance change and its largest coefficient change fall below tol.
tol (float) – IRLS iteration cap and convergence tolerance; a gene is converged when both its relative deviance change and its largest coefficient change fall below tol.
ridge (float) – L2 penalty on the covariate block.
ridge_group (float) – L2 penalty on the group block. A group whose cells carry no counts for a gene has an unbounded coefficient (the MLE is -inf); this penalty plus clip_group is what keeps it finite and defined rather than leaving it wherever the iteration happened to stop.
clip_group (float) – Bound on the magnitude of a group coefficient, on the log scale.
min_mu (float) – Floor on the fitted mean. 0 (the default) applies no floor.
dispersion_method (Literal['moments', 'cox-reid']) – How the NB dispersion is estimated when it is not supplied: "moments" or "cox-reid". Matches NBGLMBatchFitter, so the two solvers are interchangeable on a design both can express.

fit_batch(counts, *, dispersion=None, return_mu=True, return_dev_resid=True, gene_batch_size='auto')[source]

Fit one GLM per gene.

Parameters:

counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
dispersion (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – Per-gene NB dispersion to hold fixed. Ignored for Poisson. When None and family="nb", it is estimated by method of moments from a Poisson fit and the model is then refitted with it, so the returned coefficients and dispersion describe the same model.
return_mu (bool) – Whether to return the (n_samples, n_genes) fitted means and deviance residuals. Both are large; skip what is not needed.
return_dev_resid (bool) – Whether to return the (n_samples, n_genes) fitted means and deviance residuals. Both are large; skip what is not needed.
gene_batch_size (int | Literal['auto'] | None) – Genes fitted at a time. A Newton step holds a dozen (n_samples, batch) work arrays, so this is what bounds peak memory; "auto" sizes it for ~100 MB of them, as NBGLMBatchFitter does. None fits every gene in one batch. (counts itself is held whole either way.)

Return type:

StructuredGLMResult

counterfactual_means(result, *, design=None, offset=None)[source]

Fitted means with every group off, and with each group on in turn.

Answers “what would these cells have expressed under each perturbation, and under none”, which is the quantity a counterfactual or log-fold-change estimate is built from.

Returns:

baseline – Shape (n_cells, n_genes) – all group indicators zero.
per_group – Shape (n_cells, n_genes, n_groups) – group k set for all cells. This is the large one; n_cells * n_genes * n_groups float64 values.

Parameters:

result (StructuredGLMResult)
design (ndarray | None)
offset (ndarray | None)

Return type:

tuple[ndarray, ndarray]

crispyx.glm.fit_glm_onehot(counts, design, groups, *, family='nb', dispersion=None, offset=None, **kwargs)[source]

Fit [design | groups] GLMs for every gene in counts.

One-shot wrapper around StructuredGLMBatchFitter; see it for the parameters and for when this solver is the right choice.

Parameters:

counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes])
design (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes])
family (Literal['nb', 'poisson'])
dispersion (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None)
offset (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None)

Return type:

StructuredGLMResult

crispyx.glm.fit_nb_glm_batch_auto(design, counts, *, offset=None, protected_columns=(), max_iter=25, tol=1e-06, poisson_init_iter=5, dispersion_method='cox-reid', min_mu=0.0, min_total_count=1.0, ridge_penalty=1e-06, fixed_dispersion=None)[source]

Fit a batch NB GLM, using the structured solver when it pays.

A design built from a categorical covariate with many levels – batch, donor, lane – carries a block of disjoint indicator columns, which makes the per-gene Hessian arrowhead-structured and lets StructuredGLMBatchFitter take each Newton step far more cheaply than the dense path. The advantage grows with the number of groups relative to the number of remaining covariates and is not automatic (it was measured at only 1.5x for 20 groups beside 12 covariates), so the structured solver is used only when the groups outnumber the covariates.

Results are returned in the column order of design either way, so this is a drop-in replacement for NBGLMBatchFitter(design, ...).fit_batch().

Parameters:

design (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Design matrix of shape (n_samples, n_features).
counts (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]) – Count matrix of shape (n_samples, n_genes).
protected_columns (Sequence[int]) – Columns that must stay in the covariate block even if they would qualify as group indicators – the intercept and any column whose coefficient is the quantity being estimated. clip_group bounds group coefficients, which is appropriate for nuisance level effects and not for an effect under test.
offset (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – As for NBGLMBatchFitter.
max_iter (int) – As for NBGLMBatchFitter.
tol (float) – As for NBGLMBatchFitter.
poisson_init_iter (int) – As for NBGLMBatchFitter.
dispersion_method (Literal['moments', 'cox-reid']) – As for NBGLMBatchFitter.
min_mu (float) – As for NBGLMBatchFitter.
min_total_count (float) – As for NBGLMBatchFitter.
ridge_penalty (float) – As for NBGLMBatchFitter.
fixed_dispersion (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None) – Per-gene dispersion of shape (n_genes,) to hold fixed instead of estimating it; dispersion_method is then unused.

Return type:

NBGLMBatchResult

Highly variable genes

Highly variable gene (HVG) selection for streaming, disk-backed AnnData.

Producer half of the var["highly_variable"] contract that dimred.py’s PCA entry points already consume. Reimplemented as a streaming, format-aware pass rather than a thin scanpy.pp.highly_variable_genes wrapper because every scanpy flavor computes gene moments via _get_mean_var(X) against a fully realized X – there is no bounded-memory, chunked-from-disk code path, the same reason normalize_total_log1p and the QC/DE streaming paths in this package are crispyx’s own reimplementations.

crispyx.hvg.highly_variable_genes(data, *, flavor='seurat_v3', n_top_genes=2000, min_mean=0.0125, max_mean=3.0, min_disp=0.5, n_bins=20, span=0.3, layer=None, perturbation_column=None, control_label=None, cell_mask='control', gene_name_column=None, chunk_size=None, inplace=True, verbose=True)[source]

Select highly variable genes from a backed AnnData, streamed off disk.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad file, or a crispyx AnnData / backed anndata.AnnData.
flavor (Literal['seurat_v3', 'mean_dispersion']) –
HVG selection method:
- "seurat_v3" (default; Stuart et al. 2019, scanpy’s flavor="seurat_v3"): ranks genes by standardized variance fit via LOESS. Expects raw counts; requires n_top_genes.
- "mean_dispersion" (Satija et al. 2015, scanpy’s flavor="seurat"): bins genes by mean expression and z-normalizes dispersion within each bin. Expects log1p-normalized data; needs no extra dependency.
n_top_genes (int | None) – Number of top genes to select by normalized variance/dispersion. Required for flavor="seurat_v3". For flavor="mean_dispersion", if None, min_mean/max_mean/min_disp thresholds are used instead.
min_mean (float) – mean_dispersion-only threshold cutoffs used when n_top_genes is None.
max_mean (float) – mean_dispersion-only threshold cutoffs used when n_top_genes is None.
min_disp (float) – mean_dispersion-only threshold cutoffs used when n_top_genes is None.
n_bins (int) – mean_dispersion-only: number of equal-width bins used to normalize dispersion by mean expression level.
span (float) – seurat_v3-only: LOESS smoothing span (fraction of genes used to estimate the variance at each point).
layer (str | None) – If provided, use adata.layers[layer] instead of adata.X.
perturbation_column (str | None) – Column in obs containing perturbation labels. Required when cell_mask="control" (its default); ignored when cell_mask is None or an explicit array.
control_label (str | None) – Label identifying control cells. If None, auto-detected from perturbation_column. Only used when cell_mask="control".
cell_mask (Literal['control'] | ~numpy.ndarray | None) –
Which cells to compute gene moments over:
- "control" (default): control cells only, resolved from perturbation_column/control_label. This is a CRISPR/Perturb-seq-specific default – selecting HVGs from all cells lets on-target perturbation effects dominate the variable gene set, structuring downstream PCA around which perturbation a cell received rather than baseline cell-state heterogeneity.
- None: use all cells (the scanpy/Seurat default).
- An explicit boolean array: a custom cell subset.
gene_name_column (str | None) – Column in var containing gene names, used only to sanity-check that gene identifiers look like symbols, not Ensembl IDs.
chunk_size (int | None) – Number of cells (or, for CSC-stored data, genes) to process per chunk. If None, automatically calculated based on available memory.
inplace (bool) – If True on a backed AnnData, write var["highly_variable"] (bool), var["means"], var["variances"], var["variances_norm"] back to the file and return None. If False, return a DataFrame with those same columns, indexed by gene, without modifying the file.
verbose (int | bool) – Print progress and a completion summary.

Returns:

DataFrame indexed by gene when inplace=False, else None.

Return type:

pd.DataFrame or None

Raises:

ValueError – If cell_mask="control" is requested without perturbation_column, if the resolved/explicit cell mask selects zero cells, or if flavor="seurat_v3" is used without n_top_genes.

Dimension reduction

Dimension reduction module for crispyx.

Provides streaming/on-disk PCA and neighbor computation to avoid memory issues with large datasets. Follows Scanpy-style API patterns.

crispyx.dimred.pca(adata, n_comps=50, method='auto', use_highly_variable=True, chunk_size=None, random_state=0, copy=False, show_progress=True)[source]

Compute Principal Component Analysis (PCA) on backed AnnData.

Streaming implementation that works with on-disk data to avoid memory issues. Automatically selects the optimal method based on dataset characteristics.

For backed data (crispyx.AnnData wrapper), results are written directly to the h5ad file using a close-write-reopen pattern, keeping .X on disk.

Parameters:

adata (AnnData) – The annotated data matrix.
n_comps (int) – Number of principal components to compute. Default 50.
method (str) – PCA method to use: - ‘auto’: Automatically select based on gene count and memory - ‘sparse_cov’: Use sparse covariance method (fast for ≤15K genes) - ‘incremental’: Use IncrementalPCA (memory-efficient for >15K genes)
use_highly_variable (bool) – If True and ‘highly_variable’ exists in var, restrict to HVGs. Default True.
chunk_size (int | None) – Number of cells to process per chunk. If None, automatically calculated based on available memory.
random_state (int) – Random seed (not currently used, for API compatibility).
copy (bool) – If True, return a copy of adata with PCA results. If False, modify adata in place and return None.
show_progress (bool) – Show progress bars during computation. Default True.

Returns:

adata (AnnData | None) – If copy=True, returns modified AnnData. Otherwise modifies in place.
Modifies adata
————–
obsm[‘X_pca’] – PCA-transformed data (n_obs × n_comps).
varm[‘PCs’] – Principal components (n_vars × n_comps). Only includes selected genes.
uns[‘pca’] – Dict with ‘variance’, ‘variance_ratio’, ‘use_highly_variable’.

Return type:

‘AnnData’ | None

Examples

>>> import crispyx as cx
>>> adata = cx.read_backed("data.h5ad")
>>> cx.pp.pca(adata, n_comps=50)
>>> adata.obsm['X_pca'].shape
(n_obs, 50)

crispyx.dimred.neighbors(adata, n_neighbors=15, n_pcs=None, use_rep='X_pca', metric='euclidean', method='umap', random_state=0, copy=False, show_progress=True)[source]

Compute k-nearest neighbors graph from embeddings.

Uses pre-computed embeddings (typically PCA) to build a KNN graph. The embeddings are loaded into memory for efficient distance computation.

For backed data (crispyx.AnnData wrapper), results are written directly to the h5ad file using a close-write-reopen pattern, keeping .X on disk.

Parameters:

adata (AnnData) – The annotated data matrix with embeddings in .obsm.
n_neighbors (int) – Number of neighbors in the KNN graph. Default 15.
n_pcs (int | None) – Number of PCs to use from the embedding. If None, uses all.
use_rep (str) – Key in .obsm to use for distance computation. Default ‘X_pca’.
metric (str) – Distance metric. Default ‘euclidean’. Supports ‘euclidean’, ‘cosine’, ‘manhattan’, etc.
method (str) – KNN algorithm: ‘umap’ (uses pynndescent, fast approximate) or ‘sklearn’ (exact but slower). Default ‘umap’.
random_state (int) – Random seed for reproducibility.
copy (bool) – If True, return a copy with results. Otherwise modify in place.
show_progress (bool) – Show progress information. Default True.

Returns:

adata (AnnData | None) – If copy=True, returns modified AnnData. Otherwise modifies in place.
Modifies adata
————–
obsp[‘distances’] – Sparse distance matrix (n_obs × n_obs).
obsp[‘connectivities’] – Sparse connectivity matrix (n_obs × n_obs).
uns[‘neighbors’] – Dict with parameters: n_neighbors, method, metric, use_rep.

Return type:

‘AnnData’ | None

Examples

>>> import crispyx as cx
>>> adata = cx.read_backed("data.h5ad")
>>> cx.pp.pca(adata, n_comps=50)
>>> cx.pp.neighbors(adata, n_neighbors=15)
>>> adata.obsp['connectivities']
<sparse matrix (n_obs, n_obs)>

crispyx.dimred.umap(adata, min_dist=0.5, spread=1.0, n_components=2, neighbors_key='neighbors', random_state=0, copy=False)[source]

Compute UMAP embedding from pre-computed neighbor graph.

Uses the neighbor graph stored in adata.obsp (computed by cx.pp.neighbors) to create a 2D UMAP embedding. This is memory-efficient because only the neighbor graph and embedding need to be in memory, not the full expression matrix.

For backed data (crispyx.AnnData wrapper), results are written directly to the h5ad file using a close-write-reopen pattern, keeping .X on disk.

Parameters:

adata (AnnData) – The annotated data matrix with a neighbor graph computed.
min_dist (float) – The effective minimum distance between embedded points. Smaller values result in a more clustered embedding. Default 0.5.
spread (float) – The effective scale of embedded points. In combination with min_dist this determines how clustered/clumped the embedding is. Default 1.0.
n_components (int) – Number of UMAP dimensions. Default 2.
neighbors_key (str) – Key in .uns where neighbor information is stored. Default ‘neighbors’.
random_state (int) – Random seed for reproducibility.
copy (bool) – If True, return a copy with results. Otherwise modify in place.

Returns:

adata (AnnData | None) – If copy=True, returns modified AnnData. Otherwise modifies in place.
Modifies adata
————–
obsm[‘X_umap’] – UMAP embedding (n_obs × n_components).
uns[‘umap’] – Dict with UMAP parameters.

Return type:

‘AnnData’ | None

Examples

>>> import crispyx as cx
>>> adata = cx.read_backed("data.h5ad")
>>> cx.pp.pca(adata, n_comps=50)
>>> cx.pp.neighbors(adata, n_neighbors=15)
>>> cx.tl.umap(adata)
>>> adata.obsm['X_umap'].shape
(n_obs, 2)

Notes

This function wraps scanpy.tl.umap. The neighbor graph must be computed first using cx.pp.neighbors(). Memory requirements scale linearly with the number of cells: approximately 0.75MB per 1000 cells for 15 neighbors.

Plotting

Plotting utilities for crispyx with Scanpy-style helpers.

These functions are designed to work with on-disk AnnData objects and avoid loading full count matrices into memory.

crispyx.plotting.materialize_rank_genes_groups(data, *, key='rank_genes_groups', groups=None, n_genes=None, gene_symbols=None)[source]

Create a minimal in-memory AnnData with Scanpy-style rank_genes_groups.

This helper is intended for plotting only. It never loads the expression matrix from disk; instead it constructs an empty sparse matrix and injects Scanpy-compatible uns['rank_genes_groups'].

Parameters:

data (str | Path | AnnData | AnnData)
key (str)
groups (Sequence[str] | None)
n_genes (int | None)
gene_symbols (str | None)

Return type:

AnnData

crispyx.plotting.rank_genes_groups(data, groups=None, *, n_genes=20, gene_symbols=None, key='rank_genes_groups', **kwargs)[source]

Scanpy-style rank_genes_groups plot from on-disk crispyx results.

Parameters:

data (str | Path | AnnData | AnnData)
groups (Sequence[str] | None)
n_genes (int)
gene_symbols (str | None)
key (str)

crispyx.plotting.rank_genes_groups_df(data, group, *, key='rank_genes_groups', n_genes=None)[source]

Return a tidy DataFrame for rank_genes_groups results from disk.

Parameters:

data (str | Path | AnnData | AnnData)
group (str | Sequence[str])
key (str)
n_genes (int | None)

Return type:

DataFrame

crispyx.plotting.plot_volcano(*, data=None, group=None, de_df=None, key='rank_genes_groups', p_cut=0.05, lfc_cut=1.0, ax=None, show=None, savepath=None)[source]

Volcano plot for a single group.

Provide de_df directly or supply data and group to read from disk.

Parameters:

data (str | Path | AnnData | AnnData | None)
group (str | None)
de_df (DataFrame | None)
key (str)
p_cut (float)
lfc_cut (float)
show (bool | None)
savepath (str | Path | None)

crispyx.plotting.plot_top_genes_bar(*, data=None, group=None, de_df=None, key='rank_genes_groups', topn=15, ax=None, show=None, savepath=None)[source]

Horizontal bar plot of top-ranked genes for a group.

Parameters:

data (str | Path | AnnData | AnnData | None)
group (str | None)
de_df (DataFrame | None)
key (str)
topn (int)
show (bool | None)
savepath (str | Path | None)

crispyx.plotting.plot_ma(*, data, de_result=None, group, reference=None, perturbation_column=None, key='rank_genes_groups', de_df=None, mean_mode='raw', target_sum=10000.0, n_genes=None, p_cut=0.05, lfc_cut=1.0, chunk_size=1024, ax=None, show=None, savepath=None)[source]

MA plot using raw counts or normalized log1p means.

Parameters:

data (str | Path | AnnData | AnnData) – Path or backed AnnData containing raw counts.
de_result (str | Path | AnnData | AnnData | None) – Path or AnnData with rank_genes_groups results. Defaults to data.
mean_mode (str) – “raw” or “log1p” (normalized log1p means).
group (str)
reference (str | None)
perturbation_column (str | None)
key (str)
de_df (DataFrame | None)
target_sum (float)
n_genes (int | None)
p_cut (float)
lfc_cut (float)
chunk_size (int)
show (bool | None)
savepath (str | Path | None)

crispyx.plotting.plot_qc_perturbation_counts(*, data, perturbation_column, cell_mask=None, top_n=None, ax=None, show=None, savepath=None)[source]

Plot per-perturbation cell counts (optionally after QC filtering).

Parameters:

data (str | Path | AnnData | AnnData)
perturbation_column (str)
cell_mask (ndarray | None)
top_n (int | None)
show (bool | None)
savepath (str | Path | None)

crispyx.plotting.plot_qc_summary(qc_result, *, bins=50, min_genes=None, min_cells_per_gene=None, ax=None, show=None, savepath=None)[source]

Plot QC summary distributions from a QualityControlResult.

Parameters:

qc_result (QualityControlResult)
bins (int)
min_genes (int | None)
min_cells_per_gene (int | None)
show (bool | None)
savepath (str | Path | None)

crispyx.plotting.plot_pca(data, *, color=None, use_raw=None, layer=None, sort_order=True, groups=None, projection='2d', components=None, palette=None, na_color='lightgray', na_in_legend=True, size=None, frameon=None, legend_fontsize=None, legend_fontweight=None, legend_loc='right margin', legend_fontoutline=None, colorbar_loc='right', ncols=4, wspace=None, hspace=0.25, title=None, show=None, save=None, ax=None, return_fig=None, **kwargs)[source]

Plot PCA scatter from on-disk crispyx/backed AnnData or in-memory AnnData.

Wrapper around scanpy.pl.pca that works with backed and in-memory AnnData. Loads only the PCA embeddings and specified color columns.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad, crispyx.AnnData, or backed anndata.AnnData with X_pca computed.
color (str | Sequence[str] | None) – Keys for annotations of observations in .obs or variables in .var.
components (str | Sequence[str] | None) – e.g. ‘1,2’ or [‘1,2’, ‘3,4’]. Default first 2 components.
projection (str) – ‘2d’ or ‘3d’.
palette – Color palette for categorical annotations.
size (float | None) – Point size.
show (bool | None) – Show the figure.
save (str | bool | None) – Save the figure.
**kwargs – Passed to scanpy.pl.pca.
use_raw (bool | None)
layer (str | None)
sort_order (bool)
groups (str | Sequence[str] | None)
na_color (str)
na_in_legend (bool)
frameon (bool | None)
legend_fontsize (int | float | str | None)
legend_fontweight (int | str | None)
legend_loc (str)
legend_fontoutline (int | None)
colorbar_loc (str | None)
ncols (int)
wspace (float | None)
hspace (float)
title (str | Sequence[str] | None)
return_fig (bool | None)

Return type:

matplotlib.axes.Axes or list of Axes, or None if show=True.

crispyx.plotting.plot_pca_variance_ratio(data, *, n_pcs=None, log=False, show=None, save=None)[source]

Plot variance ratio explained by each PC.

Wrapper around scanpy.pl.pca_variance_ratio that works with backed AnnData.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad, crispyx.AnnData, or backed anndata.AnnData with PCA computed.
n_pcs (int | None) – Number of PCs to show. Default shows all computed.
log (bool) – Plot on log scale.
show (bool | None) – Show the figure.
save (str | bool | None) – Save the figure.

Return type:

matplotlib.axes.Axes or None if show=True.

crispyx.plotting.plot_pca_loadings(data, *, components=None, include_lowest=True, show=None, save=None)[source]

Plot gene loadings for principal components.

Wrapper around scanpy.pl.pca_loadings that works with backed and in-memory AnnData.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad, crispyx.AnnData, or backed anndata.AnnData with PCA computed.
components (int | str | Sequence[int] | None) – Which PCs to plot loadings for. e.g. [1, 2, 3] or ‘1,2,3’. Default shows first few components.
include_lowest (bool) – Show genes with lowest loadings (most negative) as well.
show (bool | None) – Show the figure.
save (str | bool | None) – Save the figure.

Return type:

matplotlib.axes.Axes or None if show=True.

crispyx.plotting.plot_umap(data, *, color=None, use_raw=None, layer=None, sort_order=True, groups=None, components=None, palette=None, na_color='lightgray', na_in_legend=True, size=None, frameon=None, legend_fontsize=None, legend_fontweight=None, legend_loc='right margin', legend_fontoutline=None, colorbar_loc='right', ncols=4, wspace=None, hspace=0.25, title=None, show=None, save=None, ax=None, return_fig=None, **kwargs)[source]

Plot UMAP embedding from on-disk crispyx/backed AnnData or in-memory AnnData.

Wrapper around scanpy.pl.umap that works with backed and in-memory AnnData. Loads only the UMAP embeddings and specified color columns.

Parameters:

data (str | Path | AnnData | AnnData) – Path to h5ad, crispyx.AnnData, or backed anndata.AnnData with X_umap computed.
color (str | Sequence[str] | None) – Keys for annotations of observations in .obs or variables in .var.
components (str | Sequence[int] | None) – Which dimensions to use (e.g. [0, 1] for first two). Default first 2.
palette – Color palette for categorical annotations.
size (float | None) – Point size.
show (bool | None) – Show the figure.
save (str | bool | None) – Save the figure.
**kwargs – Passed to scanpy.pl.umap.
use_raw (bool | None)
layer (str | None)
sort_order (bool)
groups (str | Sequence[str] | None)
na_color (str)
na_in_legend (bool)
frameon (bool | None)
legend_fontsize (int | float | str | None)
legend_fontweight (int | str | None)
legend_loc (str)
legend_fontoutline (int | None)
colorbar_loc (str | None)
ncols (int)
wspace (float | None)
hspace (float)
title (str | Sequence[str] | None)
return_fig (bool | None)

Return type:

matplotlib.axes.Axes or list of Axes, or None if show=True.

Examples

>>> import crispyx as cx
>>> adata = cx.read_backed("data.h5ad")
>>> cx.pl.umap(adata, color="perturbation")

Profiling

Profiling utilities for timing and memory measurement.

This module provides a unified Profiler class for measuring execution time and memory usage of code sections. It supports:

Timing measurement with start/stop labels
Memory tracking via tracemalloc (Python objects) or RSS (process-level)
Continuous background memory sampling
Visualization utilities for timing and memory plots

Example usage:

from crispyx.profiling import Profiler

# Basic timing
profiler = Profiler(timing=True)
profiler.start("data_loading")
# ... load data ...
profiler.stop("data_loading")
print(profiler.get_report())

# Combined timing and memory with context manager
with Profiler(timing=True, memory=True) as p:
    p.start("processing")
    # ... process data ...
    p.snapshot("after_processing")
    p.stop("processing")

# Access results programmatically
stats = p.get_stats()
print(f"Total time: {stats['timing']['total_seconds']}s")
print(f"Peak memory: {stats['memory']['peak_mb']}MB")

class crispyx.profiling.Profiler(timing=False, memory=False, memory_method='tracemalloc', sampling=False, sample_interval=0.1, top_n=10)[source]

Bases: object

Unified profiler for timing and memory measurement.

This class combines timing profiling, memory snapshots, and continuous memory sampling into a single interface. It can be used as a context manager or with explicit start/stop calls.

Parameters:

timing (bool, default=False) – Enable timing measurement of code sections.
memory (bool, default=False) – Enable memory tracking (snapshots at labeled points).
memory_method ({"tracemalloc", "rss"}, default="tracemalloc") – Method for memory measurement: - “tracemalloc”: Python object memory via tracemalloc (more detailed) - “rss”: Process resident set size via psutil (total process memory)
sampling (bool, default=False) – Enable continuous background memory sampling. Runs a daemon thread that records memory usage at sample_interval intervals.
sample_interval (float, default=0.1) – Seconds between samples when sampling=True.
top_n (int, default=10) – Number of top allocations to report (only for tracemalloc).

Examples

>>> with Profiler(timing=True, memory=True) as p:
...     p.start("section1")
...     # ... code ...
...     p.stop("section1")
>>> print(p.get_report())

start(label)[source]

Start timing a labeled section.

Parameters:: label (str)
Return type:: None

stop(label)[source]

Stop timing a labeled section and return elapsed time.

Parameters:: label (str)
Return type:: float

get_total_time()[source]

Get total elapsed time since first start() call.

Return type:: float

snapshot(label)[source]

Take a memory snapshot at the current point.

Parameters:: label (str)
Return type:: None

reset_peak()[source]

Reset peak memory tracking to current memory level.

This is useful when profiling a specific operation after loading data, so that peak memory measures only the operation’s memory usage, not the memory used by previously loaded data.

Return type:: None

start_sampling()[source]

Start background memory sampling thread.

Return type:: None

stop_sampling()[source]

Stop background memory sampling thread.

Return type:: None

get_stats()[source]

Get profiling statistics as a dict.

Returns:

Dictionary with keys "timing" and "memory".

The "timing" entry contains "total_seconds" (float) and "sections" (a mapping of label to {"seconds": float, "percent": float}).

The "memory" entry contains "peak_mb" (float), "snapshots" (a mapping of label to snapshot data), "samples" (list of timestamp/memory pairs if sampling is enabled), and "top_allocations" (if tracemalloc is used).

Return type: