Differential Expression and Enrichment

Differential Expression and Enrichment

Once single-cell RNA-seq data have been clustered and cell types assigned, the natural next question is: what genes define each cluster, and which biological pathways are they involved in? Differential expression (DE) analysis answers the first question by comparing gene expression between groups of cells—typically one cluster versus all others, or one condition versus another within the same cell type. Because single-cell count data are sparse (many zeros) and overdispersed, standard bulk RNA-seq methods must be adapted. Tools such as Seurat's FindMarkers() and Scanpy's rank_genes_groups() use statistical models (negative binomial or logistic regression) suited to this data structure, and they report log fold change, a p-value, and an adjusted p-value that corrects for testing thousands of genes simultaneously. Interpreting a list of hundreds of significant genes is difficult, so we turn to enrichment analysis . Instead of examining genes one at a time, enrichment methods ask whether a predefined gene set—such as a KEGG pathway, a Gene Ontology (GO) biological process, or a Reactome pathway—is represented more often among your significant genes than expected by chance. Over-representation analysis (ORA) uses only the significant genes and a hypergeometric or Fisher's exact test, while Gene Set Enrichment Analysis (GSEA) ranks all genes by their DE statistic and detects coordinated shifts in pathway membership, making it more sensitive to subtle but consistent changes. For pharmacy and biotech applications, this step is where data becomes insight: an upregulated "cholesterol biosynthesis" signature in a drug-treated hepatocyte cluster, for example, directly suggests a mechanism of action or an adaptive stress response. The following example performs cluster-based differential expression in R using Seurat and then runs over-representation analysis with clusterProfiler: # Differential expression with Seurat (R) library(Seurat) library(clusterProfiler) library(org.Hs.eg.db) # seurat_obj: a processed Seurat object with clusters in metadata # Find markers for cluster 2 vs. all other cells markers <- FindMarkers(seurat_obj, ident.1 = "2", only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) # Keep significant genes and convert symbols to Entrez IDs sig_genes <- rownames(markers[markers$p_val_adj < 0.05, ]) entrez_ids <- bitr(sig_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID # Over-representation analysis against KEGG pathways ego_kegg <- enrichKEGG(gene = entrez_ids, organism = "hsa", p.adjust.method = "BH") head(ego_kegg@result[, c("Description", "p.adjust", "Count")]) The equivalent workflow in Python uses Scanpy for DE and GSEA-based enrichment via gseapy: # Differential expression with Scanpy (Python) import scanpy as sc import gseapy as gp # adata: a processed AnnData object with cluster labels in adata.obs['leiden'] # Rank genes for each cluster using a Wilcoxon rank-sum test sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon") # Extract significant genes for cluster 2 result = sc.get.rank_genes_groups_df(adata, group="2") sig_genes = result[(result["pvals_adj"] < 0.05) & (result["logfoldchanges"] > 0.25)]["names"].tolist() # Gene Set Enrichment Analysis against MSigDB Hallmark pathways enr = gp.enrichr(gene_list=sig_genes, gene_sets="MSigDB_Hallmark_2020", organism="Human", outdir=None) print(enr.results[["Term", "Adjusted P-value", "Overlap"]].head()) When interpreting results, always consider the biology behind the statistics. A low adjusted p-value tells you a pathway is overrepresented, but not whether the direction of change is biologically meaningful. Separating upregulated from downregulated gene sets before enrichment, checking that the leading-edge genes are expressed in the expected cell type, and validating top markers with known literature or a dot plot are essential habits. Beware of common pitfalls too: enrichment results depend heavily on the background gene set chosen, and highly expressed housekeeping genes or mitochondrial genes can dominate DE lists if not filtered appropriately. Differential expression identifies marker genes that define clusters or respond to experimental conditions, using statistics adapted for sparse single-cell data. Multiple-testing correction (e.g., Benjamini–Hochberg adjusted p-values) is mandatory when testing thousands of genes. Enrichment analysis translates gene lists into biological meaning by testing overrepresentation of curated pathways (GO, KEGG, Reactome, MSigDB). ORA tests only significant genes, while GSEA uses the full ranked list and detects subtle, coordinated pathway shifts. Always interpret enrichment results in context: check direction of change, background gene sets, and consistency with known cell-type biology. Practice exercise: Using the PBMC 3k dataset you processed in previous lessons, run sc.tl.rank_genes_groups() on your Leiden clusters, select the top 50 positive markers for the cluster you annotated as monocytes, and perform an Enrichr analysis against the MSigDB Hallmark gene sets. Then answer the following: which three pathways are most significantly enriched, do the enriched terms match the expected biology of monocytes (e.g., inflammatory and innate immune responses), and how would excluding ribosomal and mitochondrial genes from your marker list change the results? Write a short summary of your findings as if reporting them in a research notebook.