# Bio Applied Scrna Preprocessing

> Build AnnData from 10x MEX, compute scanpy QC (pct_counts_mt, n_genes_by_counts), MAD-filter cells, normalize/log1p, select HVGs. Use for scRNA-seq preprocessing, doublet/dying-cell filtering, or AnnData QC before clustering.

- Skill: `pavel-kravchenko/bio-applied-scrna-preprocessing` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-scrna-preprocessing`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-scrna-preprocessing/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-scrna-preprocessing

---


# scRNA-seq: Quality Control and Preprocessing

## When to Use

- Loading raw 10x Genomics / Cell Ranger or STARsolo output (`matrix.mtx.gz` + `barcodes.tsv.gz` + `features.tsv.gz`) into an AnnData object.
- Computing per-cell QC metrics and deciding which cells/genes to discard before clustering.
- Choosing between fixed vs. MAD-based (median absolute deviation) QC thresholds.
- Normalizing counts (library-size CP10K, scran pooling, or Pearson residuals) and selecting highly variable genes (HVGs).
- Debugging AnnData `.obs`/`.var`/`.layers`/`.raw` structure or `var_names_make_unique()` errors.

## Version Compatibility

scanpy >= 1.10, anndata >= 0.10, Python >= 3.10. Examples below use the current scanpy API (`sc.pp.calculate_qc_metrics`, `sc.pp.highly_variable_genes`); pre-1.9 scanpy lacked `inplace=True` defaults and `sc.experimental.pp`.

## Prerequisites

```bash
pip install scanpy anndata scipy numpy pandas
```

Concepts: sparse matrices (`scipy.sparse.csr_matrix`), the AnnData data model. Related upstream step: `bio-expression-matrix-counts-ingest` (loading a raw count matrix); downstream steps: `bio-applied-dimensionality-reduction`, `bio-applied-cell-type-annotation`.

## AnnData Structure

| Slot | Content |
|------|---------|
| `adata.X` | count matrix (n_cells x n_genes), usually sparse CSR |
| `adata.obs` | per-cell metadata (QC metrics, cell type labels) |
| `adata.var` | per-gene metadata (gene symbols, `highly_variable` flag) |
| `adata.obsm` | embeddings (`X_pca`, `X_umap`) |
| `adata.uns` | unstructured metadata (color palettes, parameters) |
| `adata.layers` | additional matrices (e.g. `raw_counts` before normalization) |

Count matrix generation (FASTQ -> counts) via Cell Ranger/STARsolo: barcode correction against a whitelist (1 Hamming distance), splice-aware alignment (STAR), UMI deduplication (barcode+gene+UMI collapsed), then cell calling (EmptyDrops). Aim for 20-50k reads/cell for clustering, 500k+ for rare isoform detection.

**Goal:** Load a Cell Ranger MEX directory into AnnData and preserve raw counts for later DE testing.

**Approach:** Use `sc.read_10x_mtx`, deduplicate gene symbols immediately, then snapshot `.X` into a `raw_counts` layer before any filtering or normalization touches it.

```python
import scanpy as sc
import scipy.sparse as sp


def load_10x_counts(mex_dir: str) -> "sc.AnnData":
    """Load a Cell Ranger/STARsolo MEX directory into AnnData.

    Expects matrix.mtx.gz, barcodes.tsv.gz, features.tsv.gz in mex_dir.
    """
    adata = sc.read_10x_mtx(mex_dir, var_names="gene_symbols", cache=True)
    adata.var_names_make_unique()  # duplicate symbols -> cryptic downstream errors
    adata.layers["raw_counts"] = sp.csr_matrix(adata.X.copy())
    return adata
```

## QC Metrics and MAD-Based Filtering

| Metric | Low = problem | High = problem |
|--------|---------------|-----------------|
| `n_genes_by_counts` | empty droplet / dead cell | doublet |
| `total_counts` | dead/lysed cell | doublet |
| `pct_counts_mt` | -- | dying cell (cytoplasmic RNA leaked) |

MAD thresholds adapt to each dataset's own distribution instead of hardcoding e.g. `pct_counts_mt < 5`, which is too strict for metabolically active cell types (cardiomyocytes, hepatocytes) and too loose for others.

**Goal:** Flag low-quality cells using robust, per-dataset MAD thresholds instead of fixed cutoffs.

**Approach:** Compute standard QC metrics with `sc.pp.calculate_qc_metrics`, then flag outliers on total counts (low), genes detected (low), and mito percent (high) using median +/- n*MAD, and filter before normalizing.

```python
import numpy as np
import pandas as pd
import scanpy as sc


def mad_outlier(series: pd.Series, n_mad: int = 3, direction: str = "high") -> pd.Series:
    """Flag values more than n_mad median-absolute-deviations from the median.

    direction="high" flags values above median + n_mad*MAD (e.g. doublets, MT%).
    direction="low" flags values below median - n_mad*MAD (e.g. dead cells).
    """
    median = series.median()
    mad = (series - median).abs().median()
    if direction == "high":
        return series > median + n_mad * mad
    return series < median - n_mad * mad


def qc_filter(adata: "sc.AnnData", n_mad: int = 3) -> "sc.AnnData":
    """Compute QC metrics and drop outlier cells (filter BEFORE normalizing)."""
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True, percent_top=None)

    fail_total = mad_outlier(adata.obs["total_counts"], n_mad, direction="low")
    fail_genes = mad_outlier(adata.obs["n_genes_by_counts"], n_mad, direction="low")
    fail_mt = mad_outlier(adata.obs["pct_counts_mt"], n_mad, direction="high")

    keep = ~(fail_total | fail_genes | fail_mt)
    return adata[keep].copy()
```

Always inspect the joint scatter of `total_counts` vs. `pct_counts_mt` (and vs. `n_genes_by_counts`) before trusting the automatic cutoff — a bimodal MT% distribution usually means a real dying-cell population, not noise.

## Normalization and HVG Selection

**Goal:** Normalize for library-size differences, stabilize variance, and select the top variable genes for downstream PCA/clustering.

**Approach:** Library-size normalize to a fixed target sum (CP10K), log1p-transform (must run on normalized, not scaled, data), then pick HVGs with the Seurat-flavor mean/dispersion binning. Save the full gene set to `.raw` before subsetting to HVGs so DE testing later can use the complete matrix.

```python
import scanpy as sc

# 1. Library-size normalization (CP10K) + log transform
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)  # log(x+1): compresses dynamic range, symmetrizes fold-changes

# Alternatives for strong composition bias / variance stabilization:
# sc.external.pp.scran_normalization(adata)                     # scran pooling
# sc.experimental.pp.normalize_pearson_residuals(adata)          # SCTransform-equivalent

# 2. HVG selection (Seurat flavor: bins genes by mean expr, normalizes dispersion per bin)
sc.pp.highly_variable_genes(
    adata, flavor="seurat", n_top_genes=2000, batch_key=None
)
adata.raw = adata  # keep full matrix for DE testing before subsetting
adata_hvg = adata[:, adata.var["highly_variable"]].copy()
adata_hvg.write_h5ad("scrna_preprocessed.h5ad")
```

## Pitfalls

- **Mitochondrial threshold too strict**: a fixed `pct_counts_mt < 5%` discards real cells in metabolically active types (cardiomyocytes, hepatocytes). Always check the joint distribution with total counts, and prefer MAD-based thresholds.
- **Skipping `var_names_make_unique()`**: duplicate gene symbols cause cryptic shape-mismatch errors deep in scanpy/PCA/plotting code.
- **Normalizing before QC filtering**: low-quality cells (empty droplets, dying cells) bias the library-size estimate for every other cell. Filter first, normalize second.
- **HVG selection on raw or scaled data**: `sc.pp.highly_variable_genes(flavor="seurat")` expects log-normalized (not raw, not scaled) data — scaling first inflates dispersion estimates and shifts which genes are called variable.
- **Losing raw counts**: overwriting `.X` with normalized values without first saving to `.layers["raw_counts"]` or `.raw` makes downstream count-based DE (e.g. DESeq2/edgeR-style pseudobulk) impossible.

## See Also

- `bio-applied-single-cell-scanpy` — dimensionality reduction and clustering on the preprocessed object.
- `bio-applied-cell-type-annotation` — marker-based and reference-based cell type calling.
- `bio-applied-trajectory-analysis` — pseudotime/trajectory inference downstream of clustering.
- `bio-expression-matrix-counts-ingest` — loading raw count matrices from other formats.