Chapter 7: Marker Gene Discovery & Single-Cell Differential Expression

Chapter 7: Marker Gene Discovery & Single-Cell Differential Expression

Identifying genes that define distinct clusters (marker discovery) and detecting transcriptional changes between conditions (e.g. diseased vs. healthy controls) require robust statistical modeling that handles the high sparsity, dropout, and biological variability of single-cell sequencing.

1. Cluster Marker Discovery (FindAllMarkers)

# Fast Wilcoxon Rank-Sum test for all cluster markers
markers <- FindAllMarkers(
    pbmc, 
    only.pos = TRUE, 
    min.pct = 0.25,        # Detected in at least 25% of cells in the cluster
    logfc.threshold = 0.25 # Minimum 1.28-fold difference
)

# Filter top 3 markers per cluster
library(dplyr)
top_markers <- markers %>% 
    group_by(cluster) %>% 
    slice_max(n = 3, order_by = avg_log2FC)

# Visualize with DotPlot
DotPlot(pbmc, features = unique(top_markers$gene)) + RotatedAxis()

2. Single-Cell Hurdle Modeling with MAST

The MAST framework (Model-based Analysis of Single-cell Transcriptomics) models scRNA-seq expression using a two-part generalized linear model: a logistic regression component modeling the probability that a gene is expressed (zero vs non-zero), and a Gaussian component modeling continuous positive expression levels.

# Running MAST differential expression in Seurat
cd4_vs_cd8 <- FindMarkers(
    pbmc, 
    ident.1 = "CD4 T cells", 
    ident.2 = "CD8 T cells", 
    test.use = "MAST",
    latent.vars = c("percent.mt") # Regress out cellular quality artifacts
)
head(cd4_vs_cd8)

3. Pseudobulk Differential Expression Workflow

In multi-donor experimental designs, sum raw UMI counts across all cells belonging to the same cell type within each individual biological donor, creating a pseudobulk count matrix suitable for standard negative binomial testing in DESeq2 or edgeR.

# Pseudobulk aggregation using Seurat v5 AggregateExpression
pseudo_counts <- AggregateExpression(
    pbmc, 
    assays = "RNA", 
    group.by = c("cell_type", "donor_id", "treatment_condition"),
    return.seurat = FALSE
)

# Extract counts for specific cell type (e.g. Monocytes)
mono_counts <- pseudo_counts$RNA[, grepl("^Monocytes_", colnames(pseudo_counts$RNA))]

# Build sample metadata table
colData <- data.frame(
    sample_id = colnames(mono_counts),
    condition = ifelse(grepl("treated", colnames(mono_counts)), "Treated", "Control")
)

# Run DESeq2 at true donor-level replication
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = mono_counts, colData = colData, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Treated", "Control"))
summary(res)

4. Visualization: Volcano Plots & Complex Heatmaps

# Plotting single-cell expression heatmap
DoHeatmap(pbmc, features = top_markers$gene, size = 3, angle = 45) + 
    scale_fill_gradient2(low = "blue", mid = "white", high = "red")