Bio Applied Single Cell Scanpy

Single-Cell RNA-seq Analysis with Scanpy

When to Use

• Given a raw scRNA-seq count matrix (10x matrix.mtx/h5ad/h5ad) and need cells QC'd, normalized, clustered, and visualized end to end.
• Need to detect and remove low-quality cells (empty droplets, doublets, dying cells) before downstream analysis.
• Need to identify highly variable genes, reduce dimensionality (PCA/UMAP), and cluster cells with Leiden/Louvain.
• Need marker genes per cluster to annotate cell types (e.g. PBMC populations: T/B/NK/monocytes/DCs).
• Comparing clustering resolutions or marker-detection methods (wilcoxon vs t-test) on the same dataset.

Version Compatibility

scanpy >=1.10, anndata >=0.10, Python >=3.10. leidenalg (or igraph with flavor='igraph') required for sc.tl.leiden. Examples assume the 10x PBMC3k dataset via sc.datasets.pbmc3k().

Prerequisites

• pip install scanpy leidenalg python-igraph
• Familiarity with AnnData objects (adata.X, .obs, .var, .obsm) — see the anndata skill.
• Input as an AnnData object (cells x genes), raw (unnormalized) integer counts in adata.X.

Core Pipeline

Goal: go from a raw counts matrix to annotated clusters. Approach: QC-filter cells/genes, normalize + log-transform, select HVGs, scale, PCA, build a neighbor graph, embed with UMAP, cluster with Leiden, then rank marker genes per cluster.

import scanpy as sc
import numpy as np

sc.settings.verbosity = 1


def run_qc(adata, min_genes=200, min_cells=3, max_pct_mt=20.0):
    """Flag mitochondrial genes, compute QC metrics, and filter low-quality cells/genes.

    Returns the filtered AnnData (view materialized as a copy).
    """
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    sc.pp.calculate_qc_metrics(
        adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
    )
    sc.pp.filter_cells(adata, min_genes=min_genes)
    sc.pp.filter_genes(adata, min_cells=min_cells)
    adata = adata[adata.obs.pct_counts_mt < max_pct_mt, :].copy()
    return adata


def normalize_and_select_hvg(adata, target_sum=1e4, min_mean=0.0125, max_mean=3, min_disp=0.5):
    """Normalize counts, log-transform, then flag highly variable genes.

    IMPORTANT: HVG selection must run on normalized+logged (not scaled) data,
    and adata.raw should be saved before subsetting so DE can use full,
    unscaled expression later.
    """
    sc.pp.normalize_total(adata, target_sum=target_sum)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(
        adata, min_mean=min_mean, max_mean=max_mean, min_disp=min_disp
    )
    adata.raw = adata  # keep full normalized (unscaled) matrix for marker DE
    return adata


def cluster(adata, n_pcs=40, n_neighbors=15, resolution=0.5):
    """Subset to HVGs, scale, PCA, build neighbor graph, embed (UMAP), and cluster (Leiden)."""
    adata = adata[:, adata.var.highly_variable].copy()
    sc.pp.scale(adata, max_value=10)
    sc.tl.pca(adata, svd_solver="arpack", n_comps=50)
    sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)
    sc.tl.umap(adata)
    sc.tl.leiden(adata, resolution=resolution, flavor="igraph", n_iterations=2)
    return adata


adata = sc.datasets.pbmc3k()
adata = run_qc(adata)
adata = normalize_and_select_hvg(adata)
adata = cluster(adata)
print(f"{adata.obs['leiden'].nunique()} clusters over {adata.n_obs} cells")

Goal: find and interpret marker genes per cluster. Approach: run rank_genes_groups, then score canonical PBMC marker sets with a dotplot to assign cell-type labels.

def find_markers(adata, groupby="leiden", method="wilcoxon"):
    """Rank differentially expressed genes per cluster vs. all other clusters.

    method='wilcoxon' is the recommended default (robust, no distribution
    assumptions); use 't-test' only for a quick sanity check.
    Uses adata.raw (unscaled, normalized) automatically if present.
    """
    sc.tl.rank_genes_groups(adata, groupby, method=method)
    return adata


CANONICAL_PBMC_MARKERS = {
    "T cells (CD4+)": ["CD3D", "CD4", "IL7R"],
    "T cells (CD8+)": ["CD3D", "CD8A", "CD8B"],
    "B cells": ["CD79A", "MS4A1", "CD19"],
    "NK cells": ["GNLY", "NKG7", "NCAM1"],
    "Monocytes (CD14+)": ["CD14", "LYZ", "S100A8"],
    "Monocytes (FCGR3A+)": ["FCGR3A", "MS4A7"],
    "Dendritic cells": ["FCER1A", "CST3"],
    "Platelets": ["PPBP", "PF4"],
}

adata = find_markers(adata)
all_markers = [g for genes in CANONICAL_PBMC_MARKERS.values() for g in genes]
available = [g for g in all_markers if g in adata.raw.var_names]
sc.pl.dotplot(adata, available, groupby="leiden", standard_scale="var", figsize=(12, 4))

Leiden resolution guide: 0.2-0.4 -> few broad populations; 0.5-1.0 -> typical PBMC subtypes; 1.0-2.0 -> fine-grained subclusters (may overfit noise).

Pitfalls

• Scaling before HVG selection: inflates dispersion estimates — always normalize -> log1p -> select HVGs -> then scale.
• Losing raw counts: set adata.raw = adata right after normalize/log1p, before subsetting to HVGs and scaling; rank_genes_groups and dotplots need unscaled expression to be interpretable.
• Batch effects mistaken for biology: color the UMAP by n_genes_by_counts, total_counts, and any batch/sample covariate before trusting cluster boundaries; integrate (Harmony/BBKNN/scVI) if clusters track batch.
• Doublets not removed: an unusually high n_genes_by_counts or total_counts often indicates a doublet, not a real cell type — run explicit doublet detection (e.g. Scrublet) before clustering.
• Over/under-clustering: don't pick resolution blindly — compare 2-3 resolutions and check marker consistency before finalizing cluster count.

See Also

• anndata — underlying data structure (adata.X/.obs/.var/.obsm) manipulated throughout this pipeline.
• bio-single-cell-doublet-detection — run Scrublet/DoubletFinder before QC filtering.
• bio-single-cell-batch-integration — Harmony/BBKNN/scVI correction when samples span batches.
• bio-single-cell-cell-annotation — automated cell-type annotation beyond manual marker dotplots.