Chapter 13 — WGCNA: Weighted Gene Co-expression Network Analysis
Where this chapter sits. Companion to Lecture 13 and Tutorial 13. Prerequisites: Chapter 6 (count modeling, normalization).
13.1 The big idea
Differential expression and functional analysis (Chapters 6–7) ask “which genes change between conditions?”. WGCNA — Weighted Gene Co-expression Network Analysis — asks a different question: “which genes change together across samples?”.
Genes that consistently rise and fall together across samples are likely functionally related — they participate in the same pathway, are co-regulated by the same TF, or are responding to the same upstream signal. This “guilt by association” premise is the basis of all co-expression analysis1. WGCNA2,3 finds groups of such co-expressed genes (called modules) and connects each module to phenotypes via correlation with sample traits.
This is a systems-level analytical lens, complementary to per-gene DE. A module of 200 genes co-expressed and correlated with disease severity is often more informative than a list of 200 individually-significant DE genes. A DE list of 800 genes is hard to interpret; a handful of co-regulated modules — each summarized by an eigengene and tied to a trait — is actionable. WGCNA reframes “which genes change?” as “which transcriptional programs change?”
The key distinction between WGCNA and differential expression is that DE operates gene-by-gene — each gene is tested independently for a mean shift between conditions — whereas WGCNA operates across the covariance structure of all genes simultaneously. Two genes can each fail a DE test yet be informative together if their co-expression pattern changes with disease. Conversely, a DE gene that belongs to no co-expression module is an orphan whose mechanism is harder to infer. The two approaches are therefore genuinely complementary, not redundant: start with DE to find what changes; use WGCNA to find the programs within which those changes are embedded.
A co-expression module is a group of genes whose expression levels rise and fall together across samples, measured by Pearson correlation. Genes in the same module are assumed to be functionally linked — co-regulated by the same transcription factor, acting in the same pathway, or responding to the same upstream signal — even without prior knowledge of their function. This “guilt by association” logic1 is what makes WGCNA biologically productive: once a module is linked to a phenotype, every gene in it becomes a candidate for that biology, including uncharacterized genes co-expressed with well-known ones.
A module eigengene (ME) is the first principal component of the expression matrix formed by the module’s member genes across all samples. It is a single number per sample per module — the dominant axis of co-variation — that captures the average “activity” of the module in each sample. Because it is a PC score, the ME can be positive or negative; the sign is arbitrary, but the relative values across samples are meaningful. Once computed, an ME can be correlated with any sample-level trait (disease severity, age, treatment, etc.), converting hundreds of gene-level measurements into one tractable correlation test. Modules with similar MEs are candidates for merging (controlled by mergeCutHeight); modules with divergent MEs represent distinct biological programs.
The tutorial works through a bulk RNA-seq dataset (GSE152418, COVID-19 PBMC bulk RNA-seq from Arunachalam et al. 2020 — GEO accession GSE152418, comprising ~34 samples across healthy controls, convalescent patients, and ICU-severity COVID-19 patients) as the foundational worked example, because classic WGCNA was designed for bulk data and the interpretive logic is most transparent there. This dataset is an ideal teaching case: the severity gradient provides a clear trait axis, the sample size (~31 post-outlier-exclusion) is near the practical minimum for reliable module–trait correlation tests, and the biology (inflammatory vs. interferon modules tracking COVID severity) is well-characterized in the primary paper. For single-cell adaptations see §13.6.
13.1b Pre-processing steps before network construction
The workflow in Table 1 covers eight steps, but the first three — data pre-processing — are where subtle mistakes propagate most silently into the final network. Understanding why each step is done in this specific order is as important as knowing the function calls.
Three pre-processing steps precede network building and each has a non-obvious choice:
Outlier removal with
goodSamplesGenes(), then visual confirmation via hierarchical clustering and sample-level PCA. Remove samples that sit far from the main cluster before normalization, because outlier samples bias the variance-stabilizing transform.Variance-stabilizing transformation (VST) via DESeq2. WGCNA’s correlation is Pearson-based and assumes roughly homoscedastic data; raw counts violate this assumption because variance scales with mean. VST removes the mean–variance trend. Use
design = ~ 1(intercept only) for the DESeq2 object: WGCNA is an unsupervised method, and encoding the trait of interest in the design formula would mean the normalization step sees that variable, making the later module–trait correlation partially circular (“double dipping”).Gene filtering before running
pickSoftThreshold. The tutorial retains genes with at least 15 counts in three-quarters of the samples (rowSums(counts(dds) >= 15) >= 24). Lowly expressed genes contribute noise to pairwise correlations; filtering improves scale-free network fit and reduces the size of the TOM matrix.
Note that WGCNA expects the expression matrix as samples × genes (rows = samples, columns = genes) — the transpose of the standard RNA-seq convention of genes × samples. This is the most common orientation error and the tutorial transposes explicitly before passing to pickSoftThreshold.
Note also that WGCNA operates on a gene-by-gene correlation matrix, so statistical power depends on the number of samples, not cells. As a rule of thumb, the downstream module–trait correlation tests are meaningless below about N = 10 samples; the soft-threshold scan often fails to reach R² > 0.8 with fewer than ~15 samples; and reliable modules typically require N ≥ 20–303. If you have fewer samples than this, WGCNA will still run but its results should be treated as exploratory. This is one of the most common misapplications: a lab with 6 samples runs WGCNA, obtains colorful module-trait heatmaps with small p-values, and reports findings that do not replicate because N = 6 gives essentially no power for a Pearson correlation test.
WGCNA builds a gene × gene correlation matrix where each entry is the Pearson correlation of two genes across samples. Statistical confidence in a Pearson r with N = 10 samples is very low — the 95% CI for r = 0.8 with N = 10 is approximately ±0.25. Adding more cells per sample does not help: those cells are aggregated into one expression vector per sample, so N still equals the number of biological samples. The module–trait correlation tests that follow have the same N = samples constraint. Plan for at least 15–20 samples before committing to a WGCNA analysis, and include an independent replication cohort wherever possible.
Real biological networks — protein–protein interactions, metabolic networks, and gene co-regulatory networks — tend to be scale-free: a small number of “hub” genes have many connections, and most genes have few. A scale-free network’s degree distribution follows a power law, so the log of degree vs. the log of frequency is approximately linear with a negative slope. WGCNA’s soft-thresholding power transform is designed to push the observed co-expression network toward this scale-free property. The pickSoftThreshold() scan measures scale-free fit (R²) as a function of the chosen power β and returns the smallest β at which the network is approximately scale-free (R² ≥ 0.8). If no β achieves R² ≥ 0.8, the data does not fit the scale-free assumption — this is a signal that the dataset has unresolved batch structure, too few samples, or a biology that is not organized into hub-and-spoke modules.
The three pre-processing steps and the four network-construction/downstream steps are summarized in Table 1.
| Step | Function / action | Key choice | Common pitfall |
|---|---|---|---|
| 1. Outlier removal | goodSamplesGenes() + hierarchical clustering |
Remove outlier samples visually | Outliers bias VST downstream |
| 2. VST normalization | DESeq2::vst(design = ~ 1) |
Intercept-only design | Encoding the trait in design = circular correlation |
| 3. Gene filtering | rowSums(counts >= 15) >= 0.75 × N |
Retain genes with reliable counts | Too many low-count genes inflate TOM memory |
| 4. Soft-threshold power | pickSoftThreshold() |
Smallest \(\beta\) with \(R^2 > 0.8\) | Forcing \(\beta\) when \(R^2\) never reaches 0.8 |
| 5. Network construction | blockwiseModules() |
maxBlockSize; minModuleSize |
Memory overflow without block-wise mode |
| 6. Module eigengenes | moduleEigengenes() |
None; deterministic | Wrong orientation (samples × genes required) |
| 7. Module–trait correlation | cor() + corPvalueStudent() |
Pearson; N ≥ 10 | WGCNA::cor() shadows base::cor() globally |
| 8. Hub genes | kME or intramodularConnectivity() |
kME preferred for interpretability | TF analysis not shown; cross-reference to Chapter 7 |
Reading Table 1: work down the rows in order — each step’s output feeds the next. The most consequential user choices are in the “Key choice” column: the intercept-only design formula (row 2) prevents double-dipping; the gene-count filter (row 3) directly controls memory for the TOM step; and the soft-threshold power (row 4) determines the network’s sparsity. If you change any one of these, re-run all subsequent steps. The “Common pitfall” column maps to the error table in §13.9.
Why the intercept-only design prevents double-dipping. If you encode the severity trait in the DESeq2 design formula before running VST normalization, the normalization step will learn and partially remove variation associated with severity. The resulting normalized counts will then show reduced correlation with severity in the module–trait step — not because the biology is weak, but because you subtracted it. Using design = ~ 1 ensures the VST is agnostic to all traits and the module–trait correlations are genuine discoveries rather than artifacts of the normalization model.
The relationship between gene filtering and memory. The TOM matrix is genes × genes: for 20,000 genes it is a 400 million entry matrix, which at 8 bytes per double is ~3.2 GB. By filtering to genes with ≥ 15 counts in ≥ 75% of samples you typically reduce to 7,000–12,000 genes, shrinking the TOM to 49–144 million entries (0.4–1.2 GB), which fits comfortably in 8–16 GB RAM. If memory is still a bottleneck, blockwiseModules’s maxBlockSize argument splits the gene set into sub-blocks and computes the TOM block-by-block, accepting a small loss in cross-block module detection accuracy.
13.2 The network construction
Three steps build the gene co-expression network:
13.2.1 Pairwise correlations
Compute the Pearson correlation of every pair of genes across samples. The result is a gene × gene similarity matrix.
13.2.2 Power transform — the “weighted” in WGCNA
Raise the absolute correlation to a power \(\beta\):
\[a_{ij} = |\mathrm{cor}(x_i, x_j)|^\beta\]
This is the adjacency. Why power-transform? Real biological networks tend to be scale-free — a small number of genes have high connectivity (hubs), and most genes have low connectivity. Without a power transform, every pair gets a similarity score and the network is uniformly dense, hiding hub structure. With a power \(\beta = 6\) to \(12\), weak correlations are suppressed and strong ones are amplified, producing a network whose degree distribution looks scale-free.2 showed empirically that this soft-thresholding approach produces more biologically meaningful modules than hard-thresholded (binary) networks because it preserves the correlation gradient rather than discarding all information below a cutoff.
The right \(\beta\) is data-dependent. WGCNA’s pickSoftThreshold scans a range of \(\beta\) values and reports the scale-free fit (\(R^2\) of the linear fit between log degree and log frequency). Pick the smallest \(\beta\) at which \(R^2 > 0.8\). For typical bulk RNA-seq, that’s \(\beta = 6\)–\(12\) for unsigned networks or \(12\)–\(20\) for signed networks (which use the sign of the correlation to distinguish activation from repression).
blockwiseModules(power = ...) sets the soft-threshold exponent β. No universal default — you must choose it using pickSoftThreshold(): scan powers 1–20 (or up to 30 for signed networks) and pick the smallest β at which the scale-free topology fit R² exceeds 0.8. Typical values: β = 6–12 for unsigned networks, β = 12–20 for signed networks. Raising β further than necessary suppresses connectivity, making modules harder to detect; choosing a β at which R² < 0.8 means the network is not scale-free and downstream module biology may be unreliable.
blockwiseModules(mergeCutHeight = 0.25) controls how similar two module eigengenes must be before the modules are fused into one. Default: 0.25 (merge modules whose eigengenes have a dissimilarity ≤ 0.25, i.e., a correlation ≥ 0.75). Lowering it (e.g. 0.15) creates more, smaller, more distinct modules; raising it (e.g. 0.35–0.40) merges more aggressively, producing fewer but broader modules. If your heatmap shows many modules with near-identical trait correlations, raise mergeCutHeight. If you lose a biologically interesting sub-module, lower it.
Loading library(WGCNA) replaces base R’s cor() function with WGCNA’s own multi-threaded version in the global namespace. Default after loading: WGCNA’s cor() is active. This causes silent failures in downstream code that expects base::cor() — for example, Seurat’s FindTransferAnchors() will use WGCNA’s version and may return unexpected results. The tutorial pattern is to save and restore: temp_cor <- cor; cor <- WGCNA::cor; <WGCNA calls>; cor <- temp_cor. Always restore base::cor before any non-WGCNA code.
13.2.3 Topological overlap
Convert the adjacency to topological overlap matrix (TOM):
\[\mathrm{TOM}_{ij} = \frac{\sum_k a_{ik} a_{kj} + a_{ij}}{\min(\sum_k a_{ik}, \sum_k a_{kj}) + 1 - a_{ij}}\]
Two genes are TOM-similar if they have many shared connections, not just a high pairwise correlation. This makes the network more robust to spurious correlations. The intuition is that two genes i and j may have a modest pairwise correlation \(a_{ij}\), but if they both connect strongly to the same set of genes k, their TOM-similarity will be high — and they likely participate in the same biological process. This shared-neighborhood logic is what makes TOM better than raw adjacency for module detection.
The dissimilarity is \(1 - \mathrm{TOM}\), which is what’s clustered in the next step.
13.3 Module detection — hierarchical clustering with dynamic cutting
Run hierarchical clustering on the TOM dissimilarity matrix. The result is a dendrogram. To call modules, cut the dendrogram dynamically — cutreeDynamic finds branches whose internal connectivity is high relative to their connection to the rest of the tree. The key advantage over a fixed-height cut is that it handles branches of different sizes and connectivity levels naturally: a tight, highly-connected small module and a loose large module will both be detected even though they would require very different fixed-height thresholds.
Typical parameters:
minModuleSize = 30— modules smaller than 30 genes are merged into a “grey” module of un-grouped genesdeepSplit = 2— moderate sensitivity; higher values give more, smaller modules
The output is a colored module assignment per gene: each module gets a color label (turquoise, blue, brown, …) and the un-grouped genes get grey. The color labels themselves are arbitrary — they are simply WGCNA’s internal naming convention, and the same network run twice with a different random seed may assign a different color label to the same biological module. Never rely on a color label across datasets; always identify modules by their member genes or their eigengene correlation profile.
blockwiseModules(minModuleSize = 30) sets the minimum number of genes required for a branch to be called a module rather than assigned to grey. Default: 30. Lowering it (e.g. to 15 or 20) recovers smaller but potentially noisier modules; raising it (e.g. to 50–100) produces larger modules with better-supported eigengenes. If many genes accumulate in grey, lower minModuleSize. If you get dozens of tiny modules with no clear biological theme, raise it. For single-cell pseudobulk data with fewer genes, values of 15–20 are common.
13.3b Reading the gene dendrogram
blockwiseModules (or cutreeDynamic) produces a dendrogram with two colour bars below it: the top bar shows unmerged module assignments (many fine-grained modules as found by cutreeDynamic) and the bottom bar shows the merged assignments after modules with similar eigengenes are fused by mergeCutHeight. Reading the figure:
- Genes on nearby leaves are highly co-expressed (close in TOM dissimilarity).
- A clean result shows 5–20 distinct coloured branches in the merged bar; many more suggests
mergeCutHeightcould be raised to reduce redundant modules. - A single giant turquoise module with no further structure usually means the soft-threshold power is too low — the network is not sparse enough to resolve sub-modules.
- Grey genes at the base of the tree were not assigned to any module (
minModuleSizetoo small to form a module or below the connectivity threshold). If many genes are grey, consider loweringminModuleSize.
13.4 Module eigengenes — the per-module summary
Each module’s expression profile across samples is summarized by its module eigengene (ME) — the first principal component of the module’s expression. The ME is one number per sample per module, capturing the dominant co-variation.
In R:
MEs <- moduleEigengenes(expr, colors = moduleColors)$eigengenesThe result is a samples × modules matrix. Now you can correlate each module against any sample-level trait (treatment, disease severity, time point, etc.) — the module-trait correlation:
trait_cor <- cor(MEs, traits)
trait_pval <- corPvalueStudent(trait_cor, nSamples)Plot as a heatmap with significance asterisks. This is the canonical WGCNA figure.
13.5 Picking a soft-threshold power \(\beta\)
Practical workflow:
pickSoftThreshold(expr, powerVector = c(1:10, seq(12, 20, 2)))— runs the scan- Inspect the output: at each \(\beta\), what’s the scale-free fit \(R^2\)?
- Pick the smallest \(\beta\) at which \(R^2 > 0.8\)
- If \(R^2\) never reaches 0.8 across the scanned range, your data is too noisy or has technical structure that isn’t being modeled. Don’t force it; investigate
Typical bulk RNA-seq: \(\beta = 6\)–\(12\). Bigger studies (more samples) can use higher \(\beta\) without losing connectivity.
13.6 WGCNA on single-cell data — why classic WGCNA struggles
WGCNA was designed for bulk RNA-seq, where you have, say, 50 samples each with one expression vector. Modules emerge from inter-sample correlation.
Applied directly to scRNA-seq, classic WGCNA struggles because:
- Sparsity: at the cell level, most pairs of (gene, cell) entries are zero. Correlation of two sparse vectors is unstable.
- Scale: 50,000 cells × 20,000 genes = 1 × 10^9 entries. The TOM matrix is genes × genes = 4 × 10^8 entries. Memory becomes a problem.
- Heterogeneity: scRNA-seq has many cell types, each with its own co-expression structure. Mixing them gives a mash that doesn’t represent any single biology.
Two practical strategies:
- Pseudobulk first. Aggregate scRNA-seq counts to (sample × cell type) pseudobulk (Chapter 6) and run classic WGCNA on the per-cell-type matrices. Often works well for 8+ samples.
hdWGCNA4 — a single-cell-aware reimplementation that uses metacells (groups of ~10-20 cells aggregated by cluster) as the unit and adapts the network construction for sparsity. Modern, recommended for scRNA-seq specifically.
Table 2 compares the three approaches to co-expression network analysis across key practical dimensions.
| Approach | Input unit | Memory demand | Min recommended N | scRNA-seq sparsity handled? | R package |
|---|---|---|---|---|---|
| Classic WGCNA2,3 | Bulk sample | Moderate (TOM ∝ genes²) | ~15 samples | No — use pseudobulk first | WGCNA |
| Pseudobulk WGCNA | Sample × cell-type pseudobulk | Same as bulk | ~8 samples per cell type | Partially (aggregation reduces zeros) | WGCNA |
| hdWGCNA4 | Metacells (~10–20 cells each) | Lower per-cell-type | ~5 donors | Yes — designed for single-cell data | hdWGCNA |
Lecture 13 and the associated tutorial use the bulk workflow on the GSE152418 COVID-19 PBMC bulk RNA-seq dataset, since classic WGCNA on bulk is the foundational technique. For sc-specific WGCNA, see hdWGCNA.
13.7 Reading module-trait correlation heatmaps
A typical figure: rows = modules (turquoise, blue, brown, …), columns = traits (severity score, age, sex, time point, …), each cell colored by Pearson correlation, with \(-\log_{10}p\) as text annotation. The color scale runs from blue (strong negative correlation) through white (no correlation) to red (strong positive correlation).
Reading rules:
- A module with strong correlation to one trait + weak to others → that module is “the X-response module.” In the GSE152418 COVID dataset, modules strongly correlated with the ICU trait but not with sex or age are candidate inflammation or immune-response modules.
- A module correlated with multiple traits → either a confounded module (severity and age are correlated in many disease datasets, so a single module may reflect both) or a true general-stress signature. Partial correlation or regression adjustment can disentangle these.
- No modules correlated with anything → either no signal in the dataset, or your module detection settings need tuning. Try lowering
minModuleSize, increasing sample N, or inspecting whether a batch effect dominates the network. - Multiple modules with very similar trait patterns →
mergeCutHeightmay be too low; these modules are biologically redundant and can likely be merged. - Modules with opposite correlations to the same trait → a positive and a negative module for the same axis. These are typically reciprocal programs — one being induced as the other is repressed during the biological process. Both are informative.
Numerical confidence: a Pearson r ≈ 0.5 with N = 30 has a p-value near 0.005 (FDR-significant); r ≈ 0.7 would have p ≈ 0.00003. With N = 10 the same r = 0.5 gives p ≈ 0.14 — not significant. Always print the sample size alongside the heatmap so readers can judge effect size vs. power.
13.8 Hub genes and module export
Each module has hub genes — those most-connected within the module. Two complementary measures identify them:
- Intramodular connectivity (kWithin): total adjacency weight of a gene to all others in the same module. High kWithin = topologically central in the TOM network.
- Module membership / kME: the Pearson correlation of each gene’s expression profile with the module eigengene. kME close to 1.0 means the gene tracks the module’s average activity almost perfectly. This is the preferred hub metric in the tutorial because it is simple to compute and directly interpretable.
# kME approach (as in the tutorial)
mm <- as.data.frame(cor(norm.counts, module_eigengenes, use = "p"))
# Rank by kME for the module of interest (e.g. turquoise):
mm |> rownames_to_column("gene") |> arrange(desc(MEturquoise)) |> head(15)
# Intramodular connectivity (requires the adjacency matrix explicitly):
intramodConn <- intramodularConnectivity(adjacency, moduleColors)
hubs <- intramodConn |> arrange(desc(kWithin)) |> head(10)Hub genes are biologically interesting candidates: a module’s hub is often the TF or upstream regulator coordinating the rest. Cross-reference hub genes against TF databases for hypothesis generation. The choice between kME and kWithin is a matter of question: kME is simpler to compute and directly interpretable (“how much does this gene track the module average?”), whereas kWithin reflects topological centrality in the TOM network. For most applications kME is the preferred starting point; kWithin is more informative when you want to understand network topology rather than eigengene tracking. For functional follow-up, export the ranked hub gene list to Chapter 7’s clusterProfiler::enrichGO() pipeline — a module’s top pathway enrichment gives it a biological label that is more interpretable than its color name.
Module color names are arbitrary identifiers. WGCNA assigns colors (turquoise, blue, brown, …) as arbitrary labels to distinguish modules — the same analysis run with a different random seed may assign a different color to what is biologically the same module. Never cross-reference module colors across datasets or publications without first verifying that the member genes are the same. Always identify modules by their member gene lists or eigengene correlation profiles, not by their color labels.
cor shadowing. The WGCNA package replaces base R’s cor() function in the global environment when it is loaded. Code run after library(WGCNA) will silently use WGCNA’s cor() rather than base::cor(), which can produce surprising results in unrelated downstream code. The tutorial pattern saves base cor before calling blockwiseModules and restores it afterward: temp_cor <- cor; cor <- WGCNA::cor; ...; cor <- temp_cor.
For downstream functional analysis, export the module gene list and run Chapter 7’s ORA / GSEA pipeline on it. WGCNA modules are gene sets — they enrich for biology like any other gene set.
13.9 Common errors
Table 3 collects the errors most frequently encountered in the WGCNA workflow, their causes, and the fix.
| You see | What’s wrong | What to do |
|---|---|---|
WGCNA’s expr matrix has wrong orientation |
WGCNA expects samples in rows, genes in columns (opposite of most RNA-seq workflows) | Transpose before calling WGCNA functions |
| Scale-free fit \(R^2\) never reaches 0.8 | Data is noisy; or you have batch structure | Don’t force; inspect; consider stricter pre-filtering |
blockwiseModules runs out of memory |
TOM matrix doesn’t fit | Set maxBlockSize = 5000; accept the loss in cross-block module assignments |
| All genes end up in the grey module | minModuleSize too high or no real co-expression structure |
Lower minModuleSize to 20; investigate whether the data has clear biology |
| Module-trait correlations all p > 0.05 | Too few samples (N < 10) for the test to have power | Need bigger N; WGCNA’s downstream tests are not magic |
13.10 Where this material is also discussed
- Lecture: Lecture 13
- Tutorial: Tutorial 13
- Functional follow-up on module gene lists: Chapter 7
- Single-cell adaptation: hdWGCNA paper
13.11 Going further
The foundational papers are the weighted co-expression framework2 and the WGCNA R package3; the conceptual review of co-expression / guilt-by-association is van Dam et al.1; and the single-cell adaptation is hdWGCNA4. The full curated list is on Key Papers, Reviews & Benchmarks.