Chapter 13: Correlation & Regression Analysis in RNA-Seq: WGCNA Co-Expression Modules

1. Weighted Gene Co-Expression Network Analysis (WGCNA)

While differential expression analysis evaluates each gene individually, Weighted Gene Co-expression Network Analysis (WGCNA) reconstructs the global architecture of transcriptomic co-regulation. Instead of applying hard correlation cutoffs (which lose subtle connectivity), WGCNA raises the continuous correlation \(s_{ij}\) to a soft-thresholding power \(\beta\):

$$a_{ij} = |s_{ij}|^\beta$$

The parameter \(\beta\) is selected so the resulting network satisfies scale-free network topology (\(P(k) \sim k^{-\gamma}\)), where \(R^2 \ge 0.85\).

2. WGCNA Analysis Pipeline

library(WGCNA)
enableWGCNAThreads()

# Input: Transposed Variance-Stabilized Expression Matrix (Samples as rows, Genes as cols)
datExpr <- t(assay(vsd)[1:1000, ]) # Top 1000 variable genes for demonstration

# 1. Pick Soft-Thresholding Power Beta
powers <- c(seq(1, 10, by = 1), seq(12, 20, by = 2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 0)

# Optimal beta where scale-free topology fit index R^2 > 0.85
best_beta <- sft$powerEstimate
if (is.na(best_beta)) best_beta <- 6
cat("Selected soft-thresholding power beta:", best_beta, "\n")

# 2. Compute Adjacency and Topological Overlap Matrix (TOM)
adjacency <- adjacency(datExpr, power = best_beta)
TOM <- TOMsimilarity(adjacency)
dissTOM <- 1 - TOM

# 3. Hierarchical Clustering and Dynamic Tree Cut for Module Detection
geneTree <- hclust(as.dist(dissTOM), method = "average")
dynamicMods <- cutreeDynamic(
    dendro = geneTree,
    distM = dissTOM,
    deepSplit = 2,
    pamRespectsDendro = FALSE,
    minClusterSize = 30
)

# Convert numeric labels to colors
moduleColors <- labels2colors(dynamicMods)
table(moduleColors)

# 4. Correlate Module Eigengenes (MEs) with Clinical Traits
MEs <- moduleEigengenes(datExpr, colors = moduleColors)$eigengenes
trait_vector <- as.numeric(clinical_clean$Condition == "Treated")

moduleTraitCor <- cor(MEs, trait_vector, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(datExpr))

cat("Module-Trait Correlations with Treatment Status:\n")
print(data.frame(Correlation = round(moduleTraitCor, 3), P_Value = moduleTraitPvalue))