WGCNA tutorial: soft threshold, module detection, module-trait association and hub genes

Input data

Sample size, gene filtering and data transformation

WGCNA results are driven mainly by the input matrix. The official FAQ advises against WGCNA on fewer than 15 samples and recommends at least 20. The input is a matrix with samples in rows and genes in columns that has been normalized and log-transformed or variance-stabilized.

# R >= 4.3; WGCNA 1.74 is on CRAN; its dependencies impute, preprocessCore and GO.db are on Bioconductor
install.packages("BiocManager")
BiocManager::install(c("impute", "preprocessCore", "GO.db", "AnnotationDbi",
                       "DESeq2", "org.Hs.eg.db"))
install.packages("WGCNA")

# Load DESeq2 and other Bioconductor packages first and WGCNA last, so S4Vectors does not mask cor
suppressMessages({library(DESeq2); library(org.Hs.eg.db); library(WGCNA)})
options(stringsAsFactors = FALSE, timeout = 600)
acc <- "GSE130970"   # 78 NAFLD liver biopsies, RNA-seq, with fibrosis stage and other histology scores
dir.create("geo", showWarnings = FALSE)
f_cnt <- file.path("geo", paste0(acc, "_raw_counts.tsv.gz"))
f_mat <- file.path("geo", paste0(acc, "_series_matrix.txt.gz"))
if (!file.exists(f_cnt)) download.file(paste0(
  "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=", acc,
  "&format=file&file=", acc, "_raw_counts_GRCh38.p13_NCBI.tsv.gz"), f_cnt, mode = "wb")
if (!file.exists(f_mat)) download.file(paste0(
  "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE130nnn/", acc, "/matrix/", acc,
  "_series_matrix.txt.gz"), f_mat, mode = "wb")

cnt <- as.matrix(read.delim(f_cnt, row.names = 1, check.names = FALSE))  # 39376 x 78, row names are Entrez IDs
# 1. Drop low-count genes: count < 10 in more than 90% of samples (the example in the official FAQ)
keep <- rowSums(cnt >= 10) >= 0.1 * ncol(cnt)
sym  <- mapIds(org.Hs.eg.db, rownames(cnt), "SYMBOL", "ENTREZID")
cnt  <- cnt[keep & !is.na(sym), ]
rownames(cnt) <- make.unique(sym[rownames(cnt)])                      # 23337 genes
# 2. Variance-stabilizing transformation; never feed raw counts or un-logged FPKM into WGCNA
dds  <- DESeqDataSetFromMatrix(cnt, data.frame(row.names = colnames(cnt), x = rep(1, ncol(cnt))), ~1)
expr <- assay(vst(dds, blind = TRUE))
# 3. Keep the top 5000 genes by MAD and transpose to rows = samples, columns = genes
datExpr <- t(expr[order(apply(expr, 1, mad), decreasing = TRUE)[1:5000], ])
dim(datExpr)   # 78 5000

# Clinical traits: parse the characteristics rows of the Series Matrix
sm  <- readLines(gzfile(f_mat))
row <- function(p) lapply(sm[startsWith(sm, p)], function(x) gsub('"', "", strsplit(x, "\t")[[1]][-1]))
ch  <- row("!Sample_characteristics_ch1")
tr  <- as.data.frame(lapply(ch, function(v) sub("^[^:]+: ", "", v)))
names(tr) <- sapply(ch, function(v) sub(":.*", "", v[1]))
rownames(tr) <- row("!Sample_geo_accession")[[1]]
tr <- tr[rownames(datExpr), ]
traits <- data.frame(fibrosis = as.numeric(tr$`fibrosis stage`),
                     NAS = as.numeric(tr$`nafld activity score`),
                     steatosis = as.numeric(tr$`steatosis grade`),
                     ballooning = as.numeric(tr$`cytological ballooning grade`),
                     inflammation = as.numeric(tr$`lobular inflammation grade`),
                     age = as.numeric(tr$`age at biopsy`),
                     male = as.numeric(tr$Sex == "M"), row.names = rownames(tr))

Outlier samples

Outlier detection: sample tree plus standardized connectivity

The official tutorial sets the cut height on the sample tree by hand (15 for the female liver data, removing 1 sample). With more samples and no clearly isolated branch, standardized connectivity Z.k is more objective; the common cutoff is Z.k < -2.5.

The highest merge height in the GSE130970 sample tree is 84.4 and there is no single isolated sample; Z.k flags GSM3758009 (-6.14) and GSM3758047. Removing them leaves 76 samples. In this dataset keeping the two samples gives similar results (β is 7 in both cases, 14 modules, 1334 grey genes); outliers matter more in small datasets.

When the sample tree shows two large branches, first plot clinical and technical variables (batch, sex, sequencing date) under the tree to see whether the split comes from biological groups or technical factors, then decide whether to remove samples, correct the data or build separate networks.

Forgetting to transpose is an error that raises no error: a matrix with genes in rows and samples in columns still returns TRUE from goodSamplesGenes, and the network is then built with samples as genes. Check dim(datExpr) before analysis; the number of rows should equal the number of samples.

gsg <- goodSamplesGenes(datExpr, verbose = 0); gsg$allOK      # FALSE if genes have missing values or zero variance
sampleTree <- hclust(dist(datExpr), method = "average")
plot(sampleTree, cex = 0.6, main = "Sample clustering")       # look for samples hanging off on their own

# Standardized connectivity Z.k: samples with very low connectivity in the sample network are outliers (common cutoff -2.5)
A  <- adjacency(t(datExpr), type = "distance")
Zk <- as.numeric(scale(colSums(A) - 1))
rownames(datExpr)[Zk < -2.5]          # GSE130970: GSM3758009 GSM3758047
datExpr <- datExpr[Zk >= -2.5, ]; traits <- traits[rownames(datExpr), ]
nSamples <- nrow(datExpr)             # 76

Soft threshold

Soft-threshold β: R² cutoff, mean connectivity and empirical values

pickSoftThreshold returns the first β whose scale-free fit R² (signed R²) exceeds RsquaredCut; the default RsquaredCut is 0.85 and the default networkType is "unsigned". The networkType and correlation used to choose β must be identical to those used later in blockwiseModules.

When no power reaches the cutoff, the official FAQ (changed to more conservative values in December 2017) gives empirical β values: for unsigned and signed hybrid networks, 9 for fewer than 20 samples, 8 for 20–30, 7 for 30–40 and 6 for more than 40; for signed networks 18, 16, 14 and 12. The precondition is that batch effects, outlier samples and DEG-based gene filtering have been ruled out. On this page the FAQ value β = 12 for a signed network on 76 samples gives 13 modules, 1661 grey genes and a strongest fibrosis correlation of 0.53 (0.58 at β = 7).

The FAQ considers β below 15 reasonable for unsigned and signed hybrid networks and below 30 for signed networks. If R² never reaches 0.8 in that range and mean connectivity stays in the hundreds, the data usually contain a factor that pulls a subset of samples apart globally. This page simulated that case (adding 1.5 to 40% of genes in half of the samples): R² reaches only 0.68 at β = 30, mean connectivity is still 26 at β = 30, and β = 12 produces one large module of 2088 genes.

Signed or unsigned: the FAQ recommends signed or signed hybrid networks. An unsigned network puts positively and negatively correlated genes into the same module, so the direction of the module eigengene no longer matches every gene in the module. For correlation the FAQ recommends bicor (biweight midcorrelation) with maxPOutliers = 0.05 or 0.10; without it, when expression is strongly driven by a binary variable (disease status, genotype), bicor treats one group of samples as outliers.

allowWGCNAThreads(4)    # if parallel runs fail inside RStudio, use disableWGCNAThreads()
powers <- c(1:10, seq(12, 20, 2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed",
                         corFnc = bicor, corOptions = list(maxPOutliers = 0.05),
                         RsquaredCut = 0.85, verbose = 0)
fi <- sft$fitIndices
data.frame(power = fi$Power, R2 = round(-sign(fi$slope) * fi$SFT.R.sq, 3), mean.k = round(fi$mean.k., 1))
beta <- sft$powerEstimate      # GSE130970 (76 samples): 7

# If no power reaches the cutoff, use the empirical values from the official FAQ (December 2017 version)
if (is.na(beta)) {
  n <- nSamples; signed <- TRUE
  beta <- if (n < 20) 9 else if (n < 30) 8 else if (n < 40) 7 else 6
  if (signed) beta <- beta * 2
}

par(mfrow = c(1, 2))
plot(fi$Power, -sign(fi$slope) * fi$SFT.R.sq, type = "n", xlab = "power", ylab = "signed R^2")
text(fi$Power, -sign(fi$slope) * fi$SFT.R.sq, fi$Power, col = "red"); abline(h = 0.85, col = "red")
plot(fi$Power, fi$mean.k., type = "n", xlab = "power", ylab = "mean connectivity")
text(fi$Power, fi$mean.k., fi$Power, col = "red")

Module detection

Key blockwiseModules parameters and their measured effects

The defaults of the one-step blockwiseModules function differ from what the official tutorial uses. The table lists WGCNA 1.74 defaults, the tutorial values and the effects measured on this page, with a baseline of 76 samples, 5000 genes, signed + bicor and β = 7 (14 modules, 1369 grey genes).

This page reproduced official Tutorial I: 3600 female liver probes, 134 samples, β = 6, TOMType = "unsigned", minModuleSize = 30 and mergeCutHeight = 0.25 give 18 modules and 99 grey genes, with every module size identical to the tutorial table, in 7.9 seconds on a single thread.

β also changes the module count. On the same data, β = 2, 3, 4 and 7 give 8, 8, 12 and 14 modules. There is no correct number of modules; reporting the parameters used together with a sensitivity analysis holds up better than repeatedly tuning until one module correlates strongly with a trait.

net <- blockwiseModules(datExpr, power = beta,
         networkType = "signed", TOMType = "signed",     # networkType must match pickSoftThreshold
         corType = "bicor", maxPOutliers = 0.05,
         maxBlockSize = 6000,                             # larger than the gene count, so there is a single block
         minModuleSize = 30, deepSplit = 2, mergeCutHeight = 0.25,
         pamRespectsDendro = FALSE, numericLabels = TRUE, verbose = 3)
moduleColors <- labels2colors(net$colors)
table(moduleColors)       # GSE130970: 14 modules, 1369 genes in grey
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], "Module",
                    dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)

Memory

maxBlockSize, the cost of block splitting, and memory needs

TOM is a genes × genes matrix, so memory grows with the square of the gene count. When the gene count exceeds maxBlockSize (default 5000), blockwiseModules first pre-clusters genes into blocks with a k-means-like method, builds a network in each block and finally merges similar modules. Block splitting misassigns genes whose natural module spans blocks, and the official tutorial recommends using as few blocks as possible.

Effect of block splitting on results: comparing 15,000 genes in 4 blocks with 1 block, the adjusted Rand index (ARI) of the module assignments is 0.42; of the 27 single-block modules, 15 keep fewer than 70% of their genes together in the blockwise result. The genes of the fibrosis-related yellow module from the main analysis fall mainly into 2 modules (190 and 128 genes) in the 15,000-gene single-block result, and are spread across 3 modules (123, 92 and 70 genes) in the blockwise result. The official tutorial compared blockwise and single-block analysis on 3600 probes and found them very similar; that conclusion does not hold at tens of thousands of genes.

Memory guidance in the official tutorial: about 8000–10000 probes with 4 GB, about 20000 with 16 GB, about 30000 with 32 GB. On the 16 GB machine used here, 20,000 genes in one block completed with a 6.5 GB peak, and about 22,500 genes in one block crashed (other programs were using memory at the same time). The hard upper limit of maxBlockSize is sqrt(2^31) ≈ 46340.

When memory is short, first reduce the gene count (the top 5000–10000 genes by MAD usually capture the main co-expression structure) so that all genes fit in one block; if you really need every gene, use a machine with more memory. Describe the block boundaries in the paper if you split into blocks.

Module-trait

Computing and reading the module-trait heatmap

Each module is represented by its module eigengene (ME, the first principal component of the module expression matrix), which is then correlated with clinical traits. Each heatmap cell shows the correlation and its P value; red is positive and blue is negative.

MEs <- orderMEs(moduleEigengenes(datExpr, moduleColors)$eigengenes)
moduleTraitCor <- cor(MEs, traits, use = "p")
moduleTraitP   <- corPvalueStudent(moduleTraitCor, nSamples)
textMatrix <- paste0(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitP, 1), ")")
par(mar = c(6, 9, 3, 3))
labeledHeatmap(moduleTraitCor, xLabels = names(traits), yLabels = names(MEs), ySymbols = names(MEs),
               colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix,
               setStdMargins = FALSE, cex.text = 0.6, zlim = c(-1, 1), main = "Module-trait relationships")
# bicor also works for binary or ordinal traits, but turn off robust handling of the trait:
# bicor(MEs, traits, maxPOutliers = 0.05, robustY = FALSE)

Hub genes

GS, MM, intramodular connectivity and hub gene criteria

GS (gene significance) is the correlation of a gene's expression with the trait, MM (module membership, i.e. kME) is the correlation of a gene's expression with the module eigengene, and kWithin is intramodular connectivity. In a module significantly correlated with a trait, a significant positive correlation between GS and MM means the module's core genes are more likely to be related to the trait.

Whether |GS| > 0.2 is significant depends on sample size. The |GS| needed for two-sided P < 0.05 is 0.444 with 20 samples, 0.361 with 30, 0.279 with 50, 0.226 with 76 and 0.197 with 100; 0.2 becomes significant only at 97 samples. With fewer than 97 samples, genes passing |GS| > 0.2 may have no significant correlation with the trait; convert the cutoff according to sample size or require a significant GS P value directly.

WGCNA author Peter Langfelder wrote on the Bioconductor support site that kME and intramodular connectivity usually give very similar rankings; his group usually uses kME because it is simpler to calculate and has an associated P value; the ranking is only a guide for prioritizing follow-up, several genes with similarly high kME can all be considered hubs, and rank 3 versus rank 7 makes no real difference. The top 10 genes by kWithin in the yellow module on this page are TNFRSF12A, YWHAH, KPNA2, ARPC5, THBS2, ANXA2, TIGAR, ANXA2P2, LGALS3 and CCL20, which largely coincide with the 12 genes passing |MM| > 0.8 and |GS| > 0.2.

Within the yellow module, |MM| and |GS| correlate at 0.715 (P = 1.5e-64). A high correlation means the module's core genes are also the genes most correlated with fibrosis. When this correlation is close to 0, the module's hub genes are not necessarily related to the trait even if the ME-trait correlation is significant.

Langfelder et al. (PLoS One, 2013) compared hub gene selection with standard meta-analysis: hub genes in consensus modules were more useful for biological interpretation, while standard meta-analysis did as well as or better than the network method in validation success. When the goal is a biomarker, put GS significance ahead of MM.

module <- "yellow"; trait <- "fibrosis"
MM <- cor(datExpr, MEs[, paste0("ME", module)], use = "p")[, 1]          # module membership (kME)
GS <- cor(datExpr, traits[[trait]], use = "p")[, 1]                       # gene significance
inMod <- moduleColors == module
verboseScatterplot(abs(MM[inMod]), abs(GS[inMod]), xlab = paste("MM in", module),
                   ylab = paste("GS for", trait), col = module, abline = TRUE)
abline(v = 0.8, h = 0.2, lty = 2)

# |GS| needed for two-sided P < 0.05 at sample size n
tq <- qt(0.975, nSamples - 2); sqrt(tq^2 / (tq^2 + nSamples - 2))   # 0.226 at n = 76

# Intramodular connectivity kWithin
adj <- adjacency(datExpr, power = beta, type = "signed",
                 corFnc = "bicor", corOptions = "maxPOutliers = 0.05")
kIM <- intramodularConnectivity(adj, moduleColors)
hub <- data.frame(gene = colnames(datExpr), MM = MM, GS = GS, kWithin = kIM$kWithin)[inMod, ]
hub <- hub[order(-hub$kWithin), ]
head(hub, 10)
subset(hub, abs(MM) > 0.8 & abs(GS) > 0.2)      # GSE130970: 12 genes
cor(hub$kWithin, hub$MM, method = "spearman")    # 0.95: the two rankings are almost identical

Errors and pitfalls

Common errors, oversized grey modules and too few modules

Most error messages below were reproduced on WGCNA 1.74; the rest come from the official FAQ.

# Error: unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE)
find("cor")                 # fails when IRanges / S4Vectors come before WGCNA
cor <- WGCNA::cor           # set it temporarily, restore after blockwiseModules
net <- blockwiseModules(datExpr, power = beta, networkType = "signed", numericLabels = TRUE)
cor <- stats::cor
# Or restart R, load DESeq2 and the other packages first, then library(WGCNA) last

# Error: REAL() can only be applied to a 'numeric', not a 'integer'
storage.mode(datExpr) <- "double"   # triggered by an integer matrix (e.g. raw counts), which should not be used anyway

Peer review

Common reviewer objections and how to answer them

The most frequent objections to WGCNA papers are sample size and unvalidated hub genes. Validation can be done entirely with public data: test with modulePreservation whether modules are preserved in an independent dataset, then check the correlation of hub genes with the same trait.

# External data processed the same way into a vst matrix (rows = samples, columns = genes): GSE135251, 216 samples
common <- intersect(colnames(datExpr), colnames(datExpr2))
mp <- modulePreservation(
  list(ref = list(data = datExpr[, common]), test = list(data = datExpr2[, common])),
  list(ref = setNames(moduleColors, colnames(datExpr))[common]),
  referenceNetworks = 1, nPermutations = 50, networkType = "signed",
  corFnc = "bicor", randomSeed = 1, maxGoldModuleSize = 300, maxModuleSize = 1000, verbose = 0)
z <- mp$preservation$Z$ref.ref$inColumnsAlsoPresentIn.test
z[, c("moduleSize", "Zsummary.pres")]     # Zsummary > 10 strong, 2-10 moderate, < 2 not preserved

Chinese community

Verified practices and counterexamples from Chinese community posts

The following come from WGCNA tutorials on CSDN, Jianshu, Tencent Cloud and Zhihu columns; each has an error message or source code, and each was reproduced or checked on this page.

Hand off to an agent

Hand this workflow off to an agent

Downloading, transforming, choosing β, building the network, plotting and external validation follow fixed steps and suit a science agent; choosing the dataset, defining traits and interpreting modules remain your job.

Example one-line instruction: "Run WGCNA on GSE130970: filter low counts in the NCBI raw counts and apply vst, keep the top 5000 genes by MAD, remove outlier samples with Z.k < -2.5; choose β for a signed network with bicor (maxPOutliers = 0.05), run single-block blockwiseModules (minModuleSize 30, mergeCutHeight 0.25); make a module-trait heatmap with fibrosis, NAS, sex and age; for the most correlated module output the GS-MM scatter plot and a hub table ranked by kWithin; then run modulePreservation and hub gene-fibrosis correlations in GSE135251; add a sensitivity analysis over deepSplit 0–4."

1. 01

   Data and preprocessing

   Download counts and the Series Matrix, parse traits, apply low-count filtering, vst and MAD selection, and output the sample tree and list of outliers.

2. 02

   Choose β and build the network

   Run pickSoftThreshold and plot it, confirm that networkType matches, use a single block if memory allows, and record time and memory.

3. 03

   Module-trait and hubs

   Output the heatmap (with corrected P values), GS-MM scatter plot, hub gene table and module gene lists.

4. 04

   Validation and sensitivity analysis

   Run module preservation and hub gene correlations in an independent dataset; check whether key modules are stable across deepSplit and β.

5. 05

   Adversarial review

   Check for DEG-based filtering, whether β and network type match, whether blocks were split, whether the |GS| cutoff is significant, and whether the heatmap contains confounding modules such as sex.

• The workspace keeps the raw data, R scripts, sessionInfo, β selection plots, module assignment table, heatmap and validation results.
• You still need to check that traits are coded correctly, that outlier removal is justified and how the selected modules are interpreted biologically.
• Wet-lab validation such as qPCR and immunohistochemistry is up to you.