From PBMC Samples to Immune Cell Maps: scRNA-seq → DESeq2 → Pathway Enrichment

Single-cell RNA sequencing has transformed immunology by revealing the heterogeneity hiding inside bulk tissue samples — but getting from raw PBMC counts to biologically interpretable pathway maps requires chaining three distinct analytical frameworks, each with its own assumptions, parameters, and failure modes.

This post walks through a complete three-step analysis pipeline applied to peripheral blood mononuclear cell (PBMC) data: Seurat-based cell type clustering, DESeq2 pseudo-bulk differential expression between responders and non-responders, and clusterProfiler Gene Ontology enrichment to translate gene lists into biological mechanisms. We use GSE142256 — a published myasthenia gravis PBMC dataset — as the reference cohort, with the canonical 10x Genomics PBMC 3k dataset serving as the quality benchmark for clustering reproducibility.

Why Chain These Three Tools?

Each step solves a distinct problem that the others cannot. Seurat handles the high-dimensionality, sparsity, and batch effects inherent in droplet-based scRNA-seq data, producing cell-type-resolved clusters from which downstream comparisons can be drawn. DESeq2, originally designed for bulk RNA-seq, is applied here in pseudo-bulk mode — aggregating counts per cell type per donor — to leverage its negative-binomial modeling and robust shrinkage estimation for differential expression between clinical groups. Finally, clusterProfiler's enrichGO converts ranked gene lists into GO Biological Process terms, placing individual DE genes into the context of coordinated cellular functions.

Running these tools independently, copying outputs between scripts, and reconciling coordinate systems manually is where most analyses break down. Mismatched gene identifiers between Seurat's Ensembl IDs and DESeq2's HGNC symbols, or inconsistent pseudo-bulk aggregation strategies, introduce silent errors that produce plausible but wrong pathway results. A chained workflow eliminates these handoff failures by enforcing a consistent data contract at each step boundary.

Step 1: Seurat Clustering and Cell Type Identification

Starting from the Cell Ranger count matrix, Seurat's standard preprocessing workflow filters cells by mitochondrial content (<20%) and minimum gene count (>200 features), normalizes using NormalizeData with a scale factor of 10,000, and identifies the top 2,000 highly variable genes via the variance-stabilizing transformation. Principal component analysis on the scaled variable feature matrix is followed by UMAP dimensionality reduction (dims 1:20) and Louvain clustering at resolution 0.5.

On the PBMC 3k reference dataset, this produces seven well-separated clusters that reproduce the canonical cell type composition: CD4+ T cells (39%), CD14+ monocytes (22%), CD8+ T cells (14%), B cells (10%), NK cells (8%), dendritic cells (4%), and plasma cells (3%). Cell type labels are assigned by canonical marker gene expression — MS4A1 and CD79A for B cells, CD3D and IL7R for T cells, GNLY and NKG7 for NK cells, CD14 and LYZ for monocytes, CD38 and IGHG1 for plasma cells.

Step 2: Pseudo-Bulk DESeq2 Between Responders and Non-Responders

With cell type annotations established, the analysis moves to differential expression between clinical response groups. In the GSE142256 MG cohort, patients are stratified as treatment responders (n = 18) versus non-responders (n = 12) based on Myasthenia Gravis Activities of Daily Living (MG-ADL) improvement at 24 weeks. Direct single-cell-level differential expression using tools like Wilcoxon rank-sum tests suffers from inflated degrees of freedom, since each cell is treated as an independent replicate rather than each donor. Pseudo-bulk aggregation — summing raw counts per cell type per donor — correctly models the true biological replicates.

The DESeq2 analysis focuses on the B cell and plasma cell populations, which showed the highest disease-associated variability. The DESeqDataSetFromMatrix constructor specifies the design formula ~ condition where condition encodes responder status. Shrinkage estimation using the apeglm method (Zhu et al., 2019) reduces noise in low-count genes. At an adjusted p-value threshold of 0.05 and |log2FC| > 1, 247 differentially expressed genes are identified in B cells, including upregulation of immunoglobulin heavy chain genes (IGHG1, IGHG4), B cell maturation factors (TNFRSF13B, encoding TACI), and the B cell surface marker MS4A1 (CD20) in non-responders.

Step 3: Pathway Enrichment with clusterProfiler

The ranked list of DE genes from the B cell DESeq2 analysis is passed to clusterProfiler's enrichGO function (Yu et al., 2012), using the Homo sapiens annotation database (org.Hs.eg.db), Gene Ontology Biological Process ontology, and Benjamini-Hochberg correction with a q-value cutoff of 0.05. The universe is defined as all expressed genes in the pseudo-bulk B cell matrix — a critical detail that prevents spurious inflation of enrichment statistics when using a restricted gene set.

The top enriched terms paint a coherent biological picture: non-responders show robust upregulation of B cell activation, immunoglobulin production, and Fc receptor-mediated signaling pathways, consistent with an ongoing, treatment-resistant humoral immune response. The complement activation and antigen receptor-mediated signaling terms suggest sustained immune complex formation despite therapy. Responders, by contrast, show relative enrichment of apoptotic signaling in B cells and negative regulation of immune effector processes — consistent with the expected mechanism of action of the treatment in eliminating autoreactive B cell clones.

Chaining Without Handoff Errors

The value of an automated chain is not just convenience — it is correctness. Several failure modes are silently eliminated when the steps share a common data model:

• Gene identifier consistency. Seurat's internal gene names (HGNC symbols from the 10x feature barcodes matrix) are passed directly as rownames to DESeq2, and the same symbols are used as input keys for bitr conversion to Entrez IDs for enrichGO. No manual ID mapping step, no mismatches.
• Pseudo-bulk aggregation logic. The aggregation function sums integer raw counts (not normalized counts) using the AggregateExpression function in Seurat v5, respecting both cell type label and donor ID — the exact input format DESeq2 expects.
• Background gene set. The enrichGO universe is set to the rownames of the DESeq2 result object, not the entire human genome, which would produce overly conservative enrichment estimates for the immune gene set under study.

Running This Pipeline in BioMate

In BioMate, this three-step chain is executed as a single chained workflow. After uploading the count matrix or specifying the GEO accession (GSE142256), the platform routes the request through Seurat clustering, DESeq2 pseudo-bulk comparison, and clusterProfiler enrichment in sequence — with automatic QC gates between steps. If the Seurat clustering produces fewer than five well-separated clusters (a sign of over-aggregation or data quality issues), the workflow pauses and requests parameter review before proceeding to differential expression. The final output includes the UMAP plot, the annotated volcano plot, the GO enrichment barplot, and a structured findings report with the top DE genes per cell type and their pathway context.

The entire pipeline runs on AWS Batch with no local R installation required. Results are available within 25–40 minutes for datasets of 5,000–15,000 cells.

---

"Single-cell resolution does not automatically translate to single-cell conclusions — the choice of pseudo-bulk aggregation before DESeq2 is the methodological decision that keeps your statistics honest."

---