Scrna seq analysis

scrna-seq-analysis

When to Use

• Loading raw 10x Genomics (matrix.mtx.gz/barcodes.tsv.gz/features.tsv.gz) or Smart-seq count data into an AnnData object
• QC-filtering cells (doublets, dying cells, empty droplets) before any downstream analysis
• Clustering cells (Leiden/Louvain), running UMAP, and finding marker genes per cluster
• Assigning cell type labels manually (canonical markers) or automatically (CellTypist, SingleR)
• Inferring differentiation trajectories or RNA velocity (scVelo) from spliced/unspliced counts

Version Compatibility

scanpy ≥1.10, anndata ≥0.10, Python ≥3.10, celltypist ≥1.6, scvelo ≥0.3, leidenalg ≥0.10. For R: Seurat ≥5.0 (R ≥4.3).

Prerequisites

• pip install scanpy anndata celltypist scvelo leidenalg igraph
• Familiarity with sparse matrices (scipy.sparse) and pandas DataFrames
• Related skills: bio-single-cell-preprocessing, bio-expression-matrix-sparse-handling

Goal: Load raw 10x output, remove low-quality cells/doublets, and normalize so downstream comparisons aren't confounded by library size.

Approach: Use MAD (median absolute deviation) thresholds rather than fixed cutoffs — they adapt to each dataset's distribution instead of assuming "200 genes" or "20% MT" is universal.

import scanpy as sc
import numpy as np

def load_and_qc_filter(path, mt_prefix="MT-", n_mad=3, min_cells_per_gene=3):
    """Load 10x mtx data, compute QC metrics, and filter cells by MAD outlier thresholds.

    Returns a filtered AnnData with raw counts preserved in .layers['counts'].
    """
    adata = sc.read_10x_mtx(path, var_names="gene_symbols", cache=True)
    adata.var_names_make_unique()
    adata.layers["counts"] = adata.X.copy()

    adata.var["mt"] = adata.var_names.str.startswith(mt_prefix)
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True, percent_top=None)

    def mad_outlier(series, direction="low"):
        median, mad = series.median(), (series - series.median()).abs().median()
        return series < median - n_mad * mad if direction == "low" else series > median + n_mad * mad

    fail = (
        mad_outlier(adata.obs["total_counts"], "low")
        | mad_outlier(adata.obs["n_genes_by_counts"], "low")
        | mad_outlier(adata.obs["pct_counts_mt"], "high")
    )
    adata = adata[~fail].copy()
    sc.pp.filter_genes(adata, min_cells=min_cells_per_gene)

    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    return adata

adata = load_and_qc_filter("data/filtered_feature_bc_matrix/")
print(f"{adata.n_obs} cells x {adata.n_vars} genes after QC")

Goal: Reduce dimensionality, cluster cells, and rank marker genes per cluster.

Approach: HVG-subset → scale → PCA → neighbor graph → UMAP → Leiden; always keep the full log-normalized matrix in .raw so marker/DE tests use all genes, not just HVGs.

import scanpy as sc

def cluster_and_rank_markers(adata, n_top_genes=2000, n_pcs=40, resolution=0.5):
    """Run the standard HVG->PCA->neighbors->UMAP->Leiden pipeline and rank cluster markers.

    Mutates and returns adata with .obs['leiden'], .obsm['X_umap'], and
    .uns['rank_genes_groups'] populated.
    """
    adata.raw = adata  # full gene set retained for DE / marker tests

    sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes)
    hvg = adata[:, adata.var.highly_variable].copy()
    sc.pp.scale(hvg, max_value=10)
    sc.tl.pca(hvg, svd_solver="arpack", n_comps=50)
    sc.pp.neighbors(hvg, n_neighbors=10, n_pcs=n_pcs)
    sc.tl.umap(hvg)
    sc.tl.leiden(hvg, resolution=resolution)

    adata.obsm["X_pca"] = hvg.obsm["X_pca"]
    adata.obsm["X_umap"] = hvg.obsm["X_umap"]
    adata.obs["leiden"] = hvg.obs["leiden"].values

    sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", use_raw=True)
    return adata

adata = cluster_and_rank_markers(adata)
sc.pl.umap(adata, color=["leiden", "CD3D", "CD79A", "LYZ"])
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)

Goal: Assign cell type labels automatically instead of manually curating marker gene lists per cluster.

Approach: CellTypist runs a pre-trained logistic regression classifier directly on log1p(CP10K)-normalized data; majority_voting=True pools per-cell predictions to the Leiden cluster level to reduce single-cell noise.

import celltypist
from celltypist import models

def annotate_with_celltypist(adata, model_name="Immune_All_Low.pkl"):
    """Predict cell types with a pre-trained CellTypist model and majority-vote per cluster.

    Requires adata.X to be log1p-normalized with target_sum=1e4 (CellTypist's training norm).
    """
    models.download_models(model=model_name, force_update=False)
    model = models.Model.load(model=model_name)
    predictions = celltypist.annotate(adata, model=model, majority_voting=True)
    return predictions.to_adata()

adata = annotate_with_celltypist(adata)
sc.pl.umap(adata, color=["majority_voting", "conf_score"])

Seurat (R) Equivalent

For teams standardizing on R/Bioconductor, the same pipeline in Seurat:

library(Seurat)

# Load 10x data and build Seurat object
counts <- Read10X(data.dir = "data/filtered_feature_bc_matrix/")
so <- CreateSeuratObject(counts = counts, min.cells = 3, min.features = 200)

# QC: mitochondrial percentage, then filter
so[["percent.mt"]] <- PercentageFeatureSet(so, pattern = "^MT-")
so <- subset(so, subset = nFeature_RNA > 200 & percent.mt < 20)

# Normalize, HVG, scale, PCA
so <- NormalizeData(so, normalization.method = "LogNormalize", scale.factor = 1e4)
so <- FindVariableFeatures(so, selection.method = "vst", nfeatures = 2000)
so <- ScaleData(so)
so <- RunPCA(so, npcs = 50)

# Cluster and embed
so <- FindNeighbors(so, dims = 1:40)
so <- FindClusters(so, resolution = 0.5)
so <- RunUMAP(so, dims = 1:40)

# Marker genes per cluster (Wilcoxon rank-sum, matches sc.tl.rank_genes_groups)
markers <- FindAllMarkers(so, only.pos = TRUE, test.use = "wilcox", logfc.threshold = 0.25)
DimPlot(so, reduction = "umap", label = TRUE)

Pitfalls

• Forgetting adata.raw = adata before subsetting to HVGs — you lose the full gene set needed for marker/DE tests
• Over-clustering: start at resolution=0.3–0.5; only increase if known subtypes aren't separating
• Batch effects: color UMAP by batch/sample before interpreting clusters — unintegrated batches masquerade as biology (see bio-single-cell-batch-integration)
• MT% threshold is tissue-dependent: 5–10% for neurons, 20–25% for heart/liver/muscle; MAD-based filtering adapts automatically
• Doublets: run scDblFinder/DoubletFinder (or sc.pp.scrublet per batch) before QC filtering, not after
• CellTypist requires the exact training normalization (log1p, target_sum=1e4) — mismatched normalization silently degrades predictions

Key Formats

• AnnData (.h5ad): adata.obs (cell metadata), adata.var (gene metadata), adata.X (cells x genes matrix), adata.obsm['X_umap']/X_pca' (embeddings), adata.layers (alternate matrices)
• 10x output: barcodes.tsv.gz, features.tsv.gz, matrix.mtx.gz (Market Exchange sparse format)

See Also

• bio-single-cell-preprocessing — deeper QC/doublet-detection detail
• bio-single-cell-clustering — Leiden/Louvain resolution tuning
• bio-single-cell-cell-annotation — manual marker-based and SingleR annotation
• bio-single-cell-trajectory-inference — pseudotime and scVelo RNA velocity