Where this lecture fits

• Bonus track — self-paced after the core workshop; order is: Lec 15 → Lec 13 → Lec 16 → Lec 14 → Lec 17
• You are here: Lec 13 — finding co-expressed gene modules and tying them to phenotype
• Next: Lec 16 — Spatial Transcriptomics
• Companion tutorial: Tutorial 13 — WGCNA on GSE152418
• Companion book chapter: Chapter 13 — WGCNA

Goals of this lecture

• Understand what a co-expression network is and what makes it “weighted”
• Pick a soft-thresholding power so the network is approximately scale-free
• Build modules with dynamic tree cut, summarize them with eigengenes
• Correlate modules with sample traits — and read the resulting heatmap
• Apply WGCNA to scRNA-seq via pseudobulk or metacells (hdWGCNA)

WGCNA — the concept

WGCNA concept — co-expression network → modules → traits

• Step back from individual genes — WGCNA groups co-expressed genes into modules that move together across samples
• Each module is summarized by an eigengene, then correlated with traits/conditions
• Outputs are programs, not gene lists — much easier to interpret and act on

From correlations to a network

• For each pair of genes (i, j) compute Pearson correlation \(r_{ij}\) across samples
• Raise to a soft-thresholding power \(\beta\): \(a_{ij} = |r_{ij}|^{\beta}\) — emphasizes strong correlations, dampens weak ones
• Result is a weighted adjacency matrix — every gene is connected to every other gene, but most edges are near zero
• Convert to a topological overlap matrix (TOM) so genes that share many neighbours are clustered together

Why “weighted”?

A hard threshold (binary network) loses information: any gene pair just barely above the cutoff looks identical to a perfectly correlated one. The weighted soft-thresholding preserves the gradient — and Zhang & Horvath (2005) showed it produces more biologically meaningful modules.

Picking the soft-threshold power \(\beta\)

• Choose the smallest \(\beta\) for which the network is approximately scale-free (\(R^2 \gtrsim 0.8\))
• Typical values: 6 for unsigned, 12 for signed networks (rule of thumb from Horvath)
• Use pickSoftThreshold() — it returns the \(R^2\) vs. \(\beta\) table

library(WGCNA)
sft <- pickSoftThreshold(t(expr), powerVector = c(1:10, seq(12, 20, by = 2)),
                         networkType = "signed", verbose = 5)
plot(sft$fitIndices$Power, -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,
     xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit, R^2")
abline(h = 0.80, col = "red")

From TOM to modules

1. Compute TOM: \(\textrm{TOM}_{ij} = \frac{\sum_u a_{iu} a_{uj} + a_{ij}}{\min(k_i, k_j) + 1 - a_{ij}}\)
2. Cluster genes by 1 - TOM with average linkage
3. Cut the dendrogram with dynamic tree cut (cutreeDynamic) → modules (one per “branch”)
4. Merge modules whose eigengenes are highly correlated (> 0.75) — avoids redundant near-duplicate modules

bwnet <- blockwiseModules(t(expr),
                          power            = 18,        # from pickSoftThreshold
                          networkType      = "signed",
                          TOMType          = "signed",
                          minModuleSize    = 30,
                          mergeCutHeight   = 0.25,
                          numericLabels    = FALSE,
                          maxBlockSize     = 14000)
table(bwnet$colors)         # genes per module — colours are the WGCNA convention

Module eigengenes — the per-module “summary gene”

For each module, the first principal component of the module’s expression matrix is the module eigengene (ME). It’s a single number per sample per module — and that’s what you correlate with traits.

MEs <- moduleEigengenes(t(expr), bwnet$colors)$eigengenes

Correlate modules with phenotypes

trait <- model.matrix(~ severity, data = phenoData)[, -1]   # one-hot encoded
moduleTraitCor    <- cor(MEs, trait, use = "pairwise.complete.obs")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, n = nrow(MEs))
labeledHeatmap(Matrix      = moduleTraitCor,
               xLabels     = colnames(trait),
               yLabels     = colnames(MEs),
               textMatrix  = signif(moduleTraitPvalue, 2),
               colors      = blueWhiteRed(50))

Reading the heatmap

• Each cell: a Pearson r between a module eigengene and a trait, with the p-value
• Modules with |r| > 0.5 and p < 0.05 are candidates
• Hub genes within a strong module = highest module membership (kME) — these are your priority candidates for follow-up
• Always cross-check that hub genes make biological sense (GO / KEGG enrichment of the module)

Why classic WGCNA struggles on raw scRNA

• Sparsity: ~95% zeros makes pairwise correlations dominated by joint zeros
• No replicates: each cell is its own column → noise drowns the structure
• Solution: aggregate first, network second

Two practical strategies

# hdWGCNA sketch (signed, per-celltype)
seu <- SetupForWGCNA(seu, gene_select = "fraction", fraction = 0.05,
                     wgcna_name = "Tcells")
seu <- MetacellsByGroups(seu, group.by = "celltype",
                         k = 25, max_shared = 10, ident.group = "celltype")
seu <- NormalizeMetacells(seu)
seu <- TestSoftPowers(seu, networkType = "signed")
# soft_power here (9) is what TestSoftPowers picked on the *metacell* data;
# it differs from the bulk example's 18 because each scan chooses the smallest
# beta that reaches scale-free R^2 > 0.8 for *its* data — not a fixed default.
seu <- ConstructNetwork(seu, soft_power = 9, setDatExpr = FALSE)

cor shadowing — a common R footgun

temp_cor <- cor
cor      <- WGCNA::cor
bwnet    <- blockwiseModules(norm.counts, ...)
cor      <- temp_cor   # restore — otherwise unrelated downstream code breaks

Step 0 — Data orientation gotcha

norm.counts <- assay(dds_norm) |> t()   # now: samples x genes

Sample-level outlier detection

Before networking, look at your samples:

# Hierarchical clustering of samples in raw-count space
htree <- hclust(dist(t(data)), method = "average")
plot(htree)

# PCA — more sensitive on small N
pca <- prcomp(t(data))
plot(pca$x[, 1:2])

• The dendrogram’s y-axis is merge height (Euclidean distance). A sample merging at 2–3× the height of all others is an outlier candidate.
• PCA is more sensitive on small N — a sample can look borderline in the dendrogram but clearly isolated on PC1/PC2. Trust the PCA.
• goodSamplesGenes() additionally flags genes or samples with too many missing values.
• Drop outliers before normalization — they bias VST dispersion estimates.

Recap — reading the module–trait heatmap

• Each cell shows Pearson r between a module eigengene and a severity trait, with a p-value
• Modules with |r| > 0.5 and p < 0.05 are candidates for follow-up
• Hub genes = highest kME (gene–module correlation) within a strong module — these are your priority follow-ups
• Always run GO/pathway enrichment on the full module gene list to give it a biological label
• Multiple testing applies: prefer FDR-adjusted p-values, and check that hits survive permutation of trait labels

Checklist

1. QC: goodSamplesGenes(), drop sample/gene outliers — before normalization (outliers bias VST)
2. Normalize: VST (DESeq2::vst, design ~ 1) for bulk; log-normalize per cell type for pseudobulk
3. Cluster samples (hclust + PCA) — drop visible outliers; PCA is more sensitive on small N
4. Pick power: pickSoftThreshold to scale-free \(R^2 \gtrsim 0.8\); smallest β that crosses the line
5. Build network: blockwiseModules (signed, TOMType = “signed”, randomSeed set)
6. Restore cor after blockwiseModules (see the cor-shadowing note)
7. Correlate module eigengenes with traits; flag |r| > 0.5
8. Hub genes by kME; functional enrichment per module
9. Validate in an independent dataset before claiming biology

What to remember from Lecture 13

• WGCNA finds modules of co-expressed genes; eigengenes summarize them per sample
• Choose β to make the network approximately scale-free, not the largest β you can fit
• The interesting output is module–trait correlations + hub genes
• Don’t run WGCNA directly on raw single-cell counts — pseudobulk or hdWGCNA metacells first
• Validate in an independent dataset

Companion tutorial

Tutorial 13 — WGCNA on GSE152418 walks through a worked example: COVID-19 patient PBMC bulk RNA-seq, building a signed network, and correlating modules with disease severity.

Coming up next

• Lec 16 — Spatial Transcriptomics: gene expression with tissue coordinates
• Hands-on: Tutorial 13 — WGCNA on GSE152418

Further reading

• Langfelder & Horvath (2008) BMC Bioinformatics — the original WGCNA paper
• Morabito et al. (2023) Cell Reports Methods — hdWGCNA
• WGCNA package tutorials — Horvath’s official walkthroughs