bio-expression-matrix-sparse-handling

What this skill does

## Version Compatibility

Reference examples tested with: numpy 1.26+, scipy 1.12+, pandas 2.2+, anndata 0.10+, scanpy 1.10+, Matrix R package 1.6+, HDF5Array 1.30+ (Bioconductor), DelayedArray 0.28+, zellkonverter 1.12+, zarr-python 2.18+, dask 2024.1+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Sparse Matrix Handling

**"Store / compute on a single-cell matrix without blowing up memory"** -> Pick the sparse format that matches the access pattern (CSC for column ops, CSR for row ops), respect the R/Python convention difference (Bioconductor stores cells in columns; AnnData stores cells in rows), and use HDF5/Zarr backed mode for matrices too large for RAM.

## The Single Most Important Modern Insight -- The R/Python interop transpose is silent and catastrophic

R Bioconductor (Seurat, SingleCellExperiment) stores cells in COLUMNS and uses dgCMatrix (Column-Compressed Sparse). Python scverse (AnnData, scanpy) stores cells in ROWS and defaults to CSR (Compressed Sparse Row). Round-tripping with `anndata2ri`, `zellkonverter`, or `rpy2` triggers a TRANSPOSE under the hood -- once per direction. Two flips silently cancel. A debugging session that converts back and forth multiple times can end up with mysteriously transposed data and no error.

| Conversion | Implicit transpose |
|------------|-------------------|
| Python CSR -> R dgCMatrix | YES (rows <-> cols) |
| AnnData `.X` (cells x genes) -> SingleCellExperiment counts (genes x cells) | YES |
| Seurat `@assays$RNA@counts` (genes x cells) -> AnnData `.X` (cells x genes) | YES |
| `scipy.sparse.csr_matrix(dense)` | None (just format conversion) |
| `csr.tocsc()` / `csc.tocsr()` | None (semantically same matrix, different layout) |

The safe pattern for one-shot conversions: file-based intermediate. `adata.write('file.h5ad')` then `zellkonverter::readH5AD('file.h5ad')` -- avoids the in-memory rpy2/reticulate gymnastics, and the file roundtrip makes orientation explicit.

For a 1M-cell single-cell matrix, the implicit transpose is non-trivial -- minutes of wall time and a temporary memory peak roughly equal to nnz x 12 bytes (CSC) or x 16 bytes (CSR with int64). Avoid unnecessary transposes by aligning the format to the consumer.

Two adjacent insights:

1. **Sparse becomes inefficient above ~10-15% density.** Sparse iteration has cache-unfriendly indirection; dense iteration is sequential. Above ~10-15% nonzero, dense is often faster for most operations even though it uses more memory.
2. **AnnData backed mode quietly differs from in-memory in important ways.** `sc.read_h5ad('file.h5ad', backed='r')` returns an AnnData where `.X` is a wrapped HDF5 dataset, read-only. Many scanpy functions silently load to memory; some functions error or hang on backed mode.

## Algorithmic Taxonomy

| Format | Layout | Fast for | Slow for |
|--------|--------|----------|----------|
| dgCMatrix (CSC, R) | i (row indices), p (col pointers), x (values) | Column slicing; per-cell ops in single-cell (cells in cols); matrix-vector with col vector | Row slicing |
| dgRMatrix (CSR, R) | j (col indices), p (row pointers), x (values) | Row slicing; per-gene ops when genes in rows | Column slicing |
| dgTMatrix (COO, R triplet) | i, j, x | Random insertion when building; reading MTX | Most operations -- convert to dgCMatrix after build |
| scipy `csc_matrix` | Same as dgCMatrix | Column ops in Python | Row ops |
| scipy `csr_matrix` | Same as dgRMatrix | Row ops (NumPy convention); matmul; sklearn defaults | Column ops |
| scipy `coo_matrix` | Triplet (i, j, data) | Construction; MTX I/O | Most ops -- convert after build |
| HDF5 (h5ad, h5) | Single-file binary chunked | Random access via chunks; compression; widely supported | Cloud / parallel writes |
| Zarr | Chunked array, per-chunk file (or S3 object) | Cloud-native; parallel writes; Dask integration | Single-file simplicity |
| DelayedArray + HDF5Array | Bioc lazy evaluation over HDF5 | Out-of-memory bulk ops in R | Speed of in-memory |

## Decision Tree by Scenario

| Scenario | Recommended approach |
|----------|---------------------|
| Single-cell (>10k cells), per-cell ops | dgCMatrix in R (Bioconductor); CSR `adata.X` in Python (scanpy) |
| Bulk RNA-seq (60-80% density) | Dense -- sparse overhead exceeds benefit |
| Single-cell pseudobulk (after donor aggregation) | Dense -- now 60-80% density typically |
| 1M+ cells, can't fit in RAM | scanpy backed mode OR HDF5SummarizedExperiment + DelayedArray |
| TCGA + GTEx + recount3 scale (100k+ samples) | HDF5Array / Zarr + Dask |
| 10X CellRanger 3.0+ output | `sc.read_10x_h5()` or `Read10X_h5()` -- the .h5 is faster than the .mtx triplet |
| Cloud-native (anndata on S3, dask compute) | Zarr |
| Local workstation, single-machine | HDF5 -- faster, more widely supported |
| Building sparse matrix incrementally | COO (dgTMatrix / coo_matrix); convert to CSC/CSR after |
| R <-> Python conversion | File-based intermediate (`adata.write` then `zellkonverter::readH5AD`); aware of the transpose |

## Check Sparsity

**Goal:** Decide whether sparse is the right format given the actual data density.

**Approach:** Compute nonzero fraction; rule of thumb is sparse > ~85% (single-cell). Below that, dense often wins.

```python
import numpy as np
import scipy.sparse as sp

def sparsity(m):
    if sp.issparse(m):
        return 1 - m.nnz / (m.shape[0] * m.shape[1])
    return (m == 0).mean()

s = sparsity(adata.X)
print(f'{s:.1%} sparse')
```

Memory math:

| Data | Format | Bytes |
|------|--------|-------|
| 30k x 100k single-cell matrix, 5% density | dgCMatrix | (5% * 3e9) * 12 bytes ~= 1.8 GB |
| Same | Dense double | 24 GB |
| 60k x 100k bulk, 70% density | dgCMatrix | 50 GB (worse than dense!) |
| Same | Dense double | 48 GB |

For single-cell (typically 90-95% sparse), sparse is essential. For bulk RNA-seq (typically 60-80% density), dense is faster and not appreciably larger.

## dgCMatrix / scipy CSC / CSR

**Goal:** Construct, query, and convert sparse matrices in the format matching the consumer's expected layout.

**Approach:** `Matrix::sparseMatrix(i, j, x, dims=...)` (R) or `scipy.sparse.csr_matrix((data, (i, j)))` (Python); preserve row/column names; convert layout (CSC <-> CSR) without changing semantics.

```python
import scipy.sparse as sp
import pandas as pd

dense_df = pd.read_csv('counts.csv', index_col=0)
sparse_csr = sp.csr_matrix(dense_df.values)
sparse_csc = sp.csc_matrix(dense_df.values)

gene_names = dense_df.index.tolist()
sample_names = dense_df.columns.tolist()

sparse_csr.tocsc()
sparse_csc.tocsr()
```

```r
library(Matrix)

dense_mat <- as.matrix(read.csv('counts.csv', row.names = 1))
sparse_dgc <- as(dense_mat, 'CsparseMatrix')

class(sparse_dgc)
rownames(sparse_dgc) <- rownames(dense_mat)
colnames(sparse_dgc) <- colnames(dense_mat)
```

dgTMatrix is best for building matrices incrementally (reading MTX, parsing per-row); convert to dgCMatrix for downstream ops:

```r
mat_t <- as(triplet_data, 'TsparseMatrix')
mat_c <- as(mat_t, 'CsparseMatrix')
```

## HDF5 vs Zarr -- The Cloud-Native Shift

HDF5 (Hierarchical Data Format 5): hierarchical, single-file binary. Random access via chunks; supports compression (gzip, blosc, lz4). On-disk format for AnnData `.h5ad`, MuData `.h5mu`, 10x Genomics `.h5`, HDF5SummarizedExperiment.

Zarr: cloud-native, chunked array storage. Each chunk is a separate file (or S3 object). Parallel-write friendly; splittable by Dask. Format used by recent AnnData (`anndata.write_zarr`), SpatialData, and large-cohort consortia.

| Criterion | HDF5 | Zarr |
|-----------|------|------|
| File structure | Single binary file | Directory of