# Scrna Seq Analysis

> Process 10x/Smart-seq scRNA-seq count matrices in scanpy/AnnData end-to-end — QC filtering, normalization, HVG selection, PCA/UMAP, Leiden clustering, marker genes, CellTypist annotation, and scVelo trajectory. Use when analyzing single-cell RNA-seq data, working with .h5ad/AnnData objects or 10x matrix.mtx output, clustering cells, annotating cell types, or computing RNA velocity/pseudotime.

- Skill: `pavel-kravchenko/scrna-seq-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/scrna-seq-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/scrna-seq-analysis/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/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.

```python
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.

```python
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.

```python
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:

```r
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