Bio Single Cell Data Io

Version Compatibility

Reference examples tested with: scanpy 1.10+, Seurat 5.0+, anndata 0.10+

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.

Single-Cell Data I/O

"Load my 10X data" -> Parse a Cell Ranger matrix into an annotated object (cells, genes, counts, metadata).

• Python: sc.read_10x_mtx() / sc.read_10x_h5() -> AnnData
• R: Read10X() / Read10X_h5() -> CreateSeuratObject()

Governing Principle

The dominant failure in single-cell I/O is not a crash; it is a silent semantic change to the matrix during read or conversion. Three traps drive almost every lost-data bug.

Orientation is opposite across ecosystems. AnnData is cells x genes (observations on rows, obs indexes rows, var indexes columns); Seurat and SingleCellExperiment are genes x cells (features on rows, cells on columns). So adata.X is the transpose of LayerData(seu) and assay(sce). A faithful conversion must transpose AND swap which axis the metadata annotates; getting the transpose right but the metadata axis wrong is the single most common silent conversion bug.

Sparse storage compounds the transpose. R Matrix::dgCMatrix is CSC; scanpy conventionally stores X as CSR. Transposing a CSR matrix yields CSC for free, so a correct AnnData->Seurat hop involves both a logical transpose and a CSR<->CSC change. Forcing dense (.toarray(), as_dense= on write) on a 500k-cell x 30k-gene float32 matrix materializes ~60 GB; keep X and layers sparse and check with scipy.sparse.issparse(adata.X).

Conversion is lossy by default, and the loss is silent. Cross-ecosystem hops drop layers, obsp/varp, nested uns, and coerce categoricals to character/NA. adata.raw has its own var (its purpose is to survive HVG subsetting) and tools disagree on whether to read X or .raw.X (use_raw=), so a mismatched expectation silently uses the wrong matrix. Always diff slot inventories before and after any cross-ecosystem conversion, and keep the original file.

One more governing fact: the Cell Ranger filtered matrix is cell-CALLED, not ambient-corrected. The widespread claim that the filtered matrix is "decontaminated" is false. Keep the RAW (unfiltered) matrix, because EmptyDrops, SoupX, CellBender, and DecontX all require it and filtered-only storage is irreversible.

Choosing a Storage Format

Methods and tool maturity move fast here. Before committing a conversion route, verify the chosen package is still maintained and matches installed versions (packageVersion, pip show).

Loading 10X Cell Ranger Output

Goal: Read a Cell Ranger matrix correctly, keeping the raw matrix and non-GEX features when present.

Approach: Read the raw (unfiltered) MEX/HDF5 matrix; select stable Ensembl IDs for reproducible joins; retain Antibody/CRISPR features by disabling gex_only.

import scanpy as sc

# raw_feature_bc_matrix has every barcode (needed by EmptyDrops/SoupX/CellBender); filtered_feature_bc_matrix has only called cells
adata = sc.read_10x_mtx('raw_feature_bc_matrix/', var_names='gene_ids', gex_only=False)
# gene_ids (Ensembl) is stable across annotation releases; gene_symbols (default) is ambiguous and non-unique
# gex_only=False keeps Antibody Capture / CRISPR Guide; split later by adata.var['feature_types']
adata.var_names_make_unique()

library(Seurat)
counts <- Read10X(data.dir = 'filtered_feature_bc_matrix/')          # list when multiple feature types present
seurat_obj <- CreateSeuratObject(counts = counts, project = 'PBMC', min.cells = 3, min.features = 200)

Read functions return symbols by default. sc.read_10x_h5 has no var_names argument (symbols by default; Ensembl IDs land in var['gene_ids']). make_unique appends -1/-2 to duplicate symbols, which can mask distinct paralog/PAR loci, so prefer IDs when joining datasets.

AnnData Object Structure

Goal: Place counts, normalized values, metadata, and embeddings in the conventional slots so downstream tools find them.

Approach: Keep integer counts in layers['counts'], log-normalized values in X, and a frozen full-gene snapshot in .raw before HVG subsetting.

import anndata as ad

# X is (n_obs, n_vars) = cells x genes; obs indexes rows, var indexes columns
adata.layers['counts'] = adata.X.copy()   # integer UMIs, kept to recompute or feed count models (scVI, DESeq2)
# ... normalize_total + log1p populate X ...
adata.raw = adata                         # frozen log-normalized full-gene snapshot; survives later var-subsetting
adata = adata[:, adata.var['highly_variable']].copy()

Slot roles: X/layers align to both axes (each exactly cells x genes); obs/obsm align to cells; var/varm align to genes; obsp/varp are square pairwise graphs; uns is unstructured. adata.raw.to_adata() reconstitutes the snapshot. Slicing the parent by obs also slices raw on obs, but var-slicing does NOT shrink raw.

Seurat v5 Object Structure

Goal: Read and write the right assay layer under the v5 layers API.

Approach: Use LayerData()/$-accessors; rejoin split layers after merge() before any function expecting one layer.

counts <- LayerData(seurat_obj, layer = 'counts')      # v5; GetAssayData(slot=) is the superseded v4 form
counts <- seurat_obj[['RNA']]$counts                   # shorthand
merged <- merge(obj1, y = c(obj2, obj3), add.cell.ids = c('S1', 'S2', 'S3'))
merged <- JoinLayers(merged)                           # merge() splits layers (counts.1, counts.2); rejoin first

Seurat v5 stores counts/data/scale.data as layers in an Assay5; v4 used fixed slots via GetAssayData(slot=). After merge(), layers split per object until JoinLayers().