WGCNA Explained: Building Gene Co-Expression Networks and Finding Modules

A differential expression analysis hands you a list of genes that moved, but it says nothing about which of those genes move together. Two genes that rise and fall in lockstep across 40 samples are almost certainly wired into the same regulatory program, regardless of whether either crossed a significance cutoff. WGCNA — Weighted Gene Co-expression Network Analysis — is the R package built to find those coordinated groups. Instead of testing one gene at a time, it treats the expression matrix as a network, clusters genes into co-expression modules, collapses each module to a single representative profile, and then asks which modules track your sample traits. This post walks through the machinery that makes that work: soft-thresholding, the topological overlap matrix, and module eigengenes.

From a correlation matrix to a network

WGCNA starts where most people already are: a matrix of genes by samples, ideally variance-stabilized counts from DESeq2 or a log-normalized microarray table. The first move is to compute the pairwise correlation between all genes across samples. With N genes that is an N × N similarity matrix, and each entry is just a Pearson correlation.

The obvious next step would be to threshold that matrix — say, call two genes "connected" if their correlation exceeds 0.7 — turning it into an unweighted graph. WGCNA deliberately refuses to do this. A hard cutoff is arbitrary and throws away information: a correlation of 0.69 and one of 0.71 are treated as opposite when they are essentially identical. Instead WGCNA builds a weighted network where every pair keeps a continuous connection strength between 0 and 1. The question is how to transform a raw correlation into that weight — the job of the soft-thresholding power.

Soft-thresholding and scale-free topology

The core transformation is to raise the absolute correlation to a power: the adjacency between gene i and gene j is a_ij = |cor(i, j)|^β, where β is the soft-thresholding power. Because correlations live between 0 and 1, raising them to a power like 6 pushes weak correlations toward zero while barely touching strong ones. A correlation of 0.3 becomes 0.3^6 ≈ 0.0007 (effectively no edge), while 0.9 becomes 0.9^6 ≈ 0.53 (still a substantial edge). The power acts as a soft, continuous filter rather than a cliff.

How do you pick β? WGCNA anchors the choice in a biological expectation: real gene networks are approximately scale-free, meaning most genes have few connections and a handful of hub genes have many. The pickSoftThreshold function sweeps a range of candidate powers and, for each, computes how well the resulting network fits a scale-free topology, reported as a model fit R^2. You choose the lowest power at which the fit R^2 crosses roughly 0.8 to 0.9 and the mean connectivity has leveled off.

Signed versus unsigned networks

Before you compute adjacency you have to decide whether negative correlations count as connections. In an unsigned network the weight uses the absolute correlation, so a gene pair that is strongly anti-correlated (cor = -0.9) becomes just as connected as a co-activated pair. In a signed network, negative correlations are pushed toward zero, so modules contain only genes moving in the same direction.

Most current practice favors signed networks. Anti-correlated genes usually belong to functionally distinct programs (one repressed while another activates), and lumping them into one module blurs interpretation. The tradeoff is that signed networks need a higher power — roughly double — to achieve comparable scale-free fit, hence the β ≈ 12 figure in the table above.

The topological overlap matrix (TOM)

Adjacency alone is a noisy basis for clustering, because it captures only the direct correlation between two genes, and direct correlations are easily corrupted by noise in small sample sets. WGCNA's key refinement is the topological overlap matrix. The topological overlap between two genes asks not only "how connected are these two directly?" but "how many neighbors do they share?" Two genes both strongly connected to the same 30 other genes have high topological overlap even if their direct correlation is only moderate.

Formally the TOM combines the direct adjacency a_ij with the shared connection strength summed over all common neighbors, normalized by the smaller of the two genes' total connectivity. The effect is a denoised similarity: shared network context reinforces genuine relationships and washes out spurious correlations. WGCNA then converts TOM to a dissimilarity, 1 - TOM, and uses that as the distance for hierarchical clustering.

Cutting the dendrogram into modules

Clustering on 1 - TOM with average linkage produces a dendrogram where tight clusters of high-overlap genes hang together as distinct branches. Modules are those branches. Rather than cutting the tree at a single fixed height — which would badly split nested or unevenly sized clusters — WGCNA uses the Dynamic Tree Cut algorithm, which adaptively detects branch boundaries by shape.

Each detected module is assigned a color label (turquoise, blue, brown, and so on, ordered by size), and genes that fit no module land in grey. Two parameters govern the result: minModuleSize (often 30) sets the smallest allowed module, and deepSplit (0 to 4) tunes how aggressively branches are subdivided. The grey module is not a real module, just the leftover bin — never interpret it as a biological unit.

Module eigengenes and merging

A module might hold 500 genes, which is unwieldy to correlate against traits. WGCNA summarizes each module with a single profile, the module eigengene: the first principal component of that module's expression matrix across samples. Because module genes are by construction highly co-expressed, this first PC captures most of their shared variance, making the eigengene a faithful stand-in for the module's behavior.

Eigengenes also solve a practical problem. Dynamic Tree Cut sometimes splits one biological program into two near-identical modules. By correlating the eigengenes themselves and merging any pair above a threshold (via mergeCloseModules with a cutHeight such as 0.25, i.e. eigengene correlation above 0.75), you collapse redundant modules into one. The eigengene is the atom of everything downstream: what you correlate, plot, and test.

Linking modules to traits with a heatmap

The payoff step is the module-trait heatmap. You correlate each module eigengene against each sample trait — genotype, treatment, clinical score — producing a matrix of correlation coefficients and their p-values. Rendered with labeledHeatmap, each cell is colored by the correlation (red for positive, blue for negative) and annotated with the coefficient and p-value. A module that lights up strongly against "disease status" is your candidate co-regulated program driving that phenotype.

From there you drill into the significant module: rank member genes by module membership (kME, each gene's correlation to its own eigengene) and gene significance (each gene's correlation to the trait). Genes high in both are the intramodular hub genes most likely to be biologically central. Those hubs feed naturally into downstream analysis: run them through clusterProfiler for GO and KEGG enrichment to name the program, or push the module's genes into a STRING/Cytoscape PPI network to see the physical wiring.

Sample size and preprocessing caveats

WGCNA is correlation-driven, so it needs samples. Correlations across 6 samples are meaningless; the authors recommend at least 15 and prefer above 20. Filter aggressively first — drop low-count, low-variance genes to keep the informative 5000 to 15000, use goodSamplesGenes to flag genes and samples with too much missing data, and cluster samples to catch outliers. Feed WGCNA normalized input (variance-stabilized or log counts, never raw counts) and confirm no batch effect is masquerading as biology, since correlation structure faithfully reproduces whatever technical artifact you leave in.

Practical Takeaways

• WGCNA builds a weighted network — raise correlations to a soft-thresholding power β chosen for scale-free fit (R^2 ≈ 0.8+) instead of applying a hard correlation cutoff.
• Prefer a signed network so modules contain only co-directional genes; expect to roughly double the power (β ≈ 12) versus unsigned.
• The TOM denoises pairwise adjacency by rewarding shared neighbors; cluster on 1 - TOM and cut with Dynamic Tree Cut, not a fixed height.
• The module eigengene (first PC) is the atom of the analysis — you merge modules, correlate traits, and rank hub genes all through it.
• Read the module-trait heatmap for modules that track your phenotype, then use kME and gene significance to find intramodular hubs.
• WGCNA needs samples: aim for 20+, filter to variance-informative genes, and kill batch effects before you trust a single module.

WGCNA assumes bulk RNA-seq or microarray input, but the same soft-thresholding and TOM machinery now drives hdWGCNA, which builds modules from pseudobulked metacells in single-cell RNA-seq data. Whichever input you start from, treat a module as a hypothesis, not a conclusion — validate its hub genes independently before they anchor a paper's narrative.