Tutorial

Identifying Cell Types in Single-Cell RNA-Seq with Scanpy

Difficulty: Beginner Time: 15 min read

Introduction

Single-cell RNA sequencing (scRNA-seq) has transformed biology by allowing researchers to measure the expression levels of thousands of genes in individual cells. This makes it possible to discover new, rare cell types that would be masked in bulk tissue sequencing. In this tutorial, we will use the Python package Scanpy to process a raw count matrix, perform dimensionality reduction, cluster the cells, and identify distinct cell types using marker genes.

Architecture Overview


graph TD
    Raw[(Raw Count Matrix)] --> QC[Quality Control & Filtering]
    QC --> Norm[Normalization & Log Transform]
    Norm --> HVG[Highly Variable Gene Selection]
    HVG --> Scale[Scale & Center Data]
    Scale --> PCA[Principal Component Analysis]
    PCA --> Graph[K-Nearest Neighbors Graph]
    
    Graph --> UMAP[UMAP 2D Projection]
    Graph --> Leiden[Leiden Clustering]
    
    UMAP -.-> Visual[Visual Analysis]
    Leiden -.-> Visual
    Leiden --> Markers[Differential Expression: Rank Genes Groups]
    Markers --> Annotation(Cell Type Annotation)

Prerequisites

  • Python 3.8+
  • Basic familiarity with pandas DataFrames
  • A working Jupyter Notebook environment

Step 1: Install Dependencies

pip install scanpy anndata pandas matplotlib seaborn

Step 2: Loading the Data

Scanpy uses the AnnData (Annotated Data) object, which stores the main count matrix alongside annotations for observations (cells) and variables (genes). For this tutorial, we will use a built-in dataset of 3k Peripheral Blood Mononuclear Cells (PBMCs) from 10x Genomics.

import scanpy as sc
import pandas as pd

# Set verbosity and matplotlib settings
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=80, facecolor='white')

# Load the PBMC 3k dataset
adata = sc.datasets.pbmc3k()

# Make variable names (gene names) unique
adata.var_names_make_unique()
print(adata)
# Output: AnnData object with n_obs x n_vars = 2700 x 32738

Step 3: Quality Control and Preprocessing

Raw scRNA-seq data contains empty droplets and dead cells. Dead or dying cells typically leak their cytoplasmic RNA but retain their mitochondrial RNA, resulting in a high percentage of mitochondrial genes. We must filter these out.

# Calculate the percentage of mitochondrial genes
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)

# Filter cells with less than 200 genes or more than 2500 genes
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_cells(adata, max_genes=2500)

# Filter cells with > 5% mitochondrial counts (likely dying cells)
adata = adata[adata.obs.pct_counts_mt < 5, :]

# Filter out genes expressed in fewer than 3 cells
sc.pp.filter_genes(adata, min_cells=3)

# Save the raw data before normalization for differential expression testing later
adata.raw = adata

Step 4: Normalization and Feature Selection

Because different cells are sequenced to different depths, we must normalize the data. We then log-transform it and select the highly variable genes, which contain the most information about cell types.

# Normalize counts per cell to 10,000
sc.pp.normalize_total(adata, target_sum=1e4)

# Logarithmize the data
sc.pp.log1p(adata)

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

# Actually filter the data to only keep these genes
adata = adata[:, adata.var.highly_variable]

# Scale the data to unit variance and zero mean
sc.pp.scale(adata, max_value=10)

Step 5: Dimensionality Reduction and Clustering

We reduce the ~2000 highly variable genes down to 50 Principal Components (PCA). We then build a neighborhood graph and cluster the cells using the Leiden algorithm.

# Run PCA
sc.tl.pca(adata, svd_solver='arpack')

# Compute the neighborhood graph using the first 40 PCs
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

# Run UMAP for 2D visualization
sc.tl.umap(adata)

# Cluster the cells using the Leiden algorithm
sc.tl.leiden(adata, resolution=0.5)

# Plot the UMAP colored by the Leiden clusters
sc.pl.umap(adata, color=['leiden'])

Step 6: Marker Gene Identification and Cell Type Annotation

Now that we have distinct clusters, we need to figure out what cell type each cluster represents. We do this by finding the differentially expressed genes (marker genes) for each cluster.

# Rank genes to find markers for each Leiden cluster
# We use the raw, unscaled data for statistical testing
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', use_raw=True)

# Plot the top 5 marker genes for each cluster
sc.pl.rank_genes_groups(adata, n_genes=5, sharey=False)

# Based on literature, we can map markers to known cell types:
# Cluster 0: IL7R, CD3D (CD4 T Cells)
# Cluster 1: CD14, LYZ (CD14 Monocytes)
# Cluster 2: MS4A1 (B Cells)
# Cluster 3: CD8A (CD8 T Cells)

# Create a dictionary to map cluster numbers to cell types
new_cluster_names = {
    '0': 'CD4 T',
    '1': 'CD14 Monocytes',
    '2': 'B',
    '3': 'CD8 T',
    # ... add mappings for all clusters
}

# Apply the new annotations to the AnnData object
adata.rename_categories('leiden', new_cluster_names.values())

# Plot the final annotated UMAP
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='PBMC Cell Types')

Conclusion

With just a few lines of code, Scanpy allows us to take raw, noisy biological data of 2,700 cells across 32,000 genes, clean it, and project it into a beautiful 2D UMAP where distinct cell types emerge organically. By combining Leiden clustering with Wilcoxon rank-sum tests for differential expression, we can confidently annotate our clusters to discover the underlying biology of the tissue sample.