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

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

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.

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

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.

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.

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.