Chapter 6: Module–Trait Associations, Matrix-Wide BH FDR Correction & Linear Modeling
Chapter 6: Module–Trait Associations, Matrix-Wide BH FDR Correction & Linear Modeling
Calculating Module Eigengenes (MEs); computing Pearson correlations against continuous and binary physiological stress traits; applying Benjamini-Hochberg FDR correction across the entire test matrix; and fitting linear models (lm(ME ~ treatment * timepoint)).
1. Module Eigengenes (MEs) & Biological Trait Formulations
To test whether co-expression modules correlate with phenotypic stress responses, WGCNA summarizes the collective expression of all genes in module $q$ into a single characteristic synthetic vector: the Module Eigengene ($ME_q$).
Mathematically, $ME_q$ is defined as the first principal component ($PC1$) of the standardized expression matrix of module $q$. It captures the maximum possible variance (typically $40\text{--}65\%$) of intramodular expression trajectories across all experimental samples.
# Calculate module eigengenes across all samples
MEs <- moduleEigengenes(expr_mat, colors = module_colors[colnames(expr_mat)])$eigengenes
# Order eigengenes so closely related modules cluster adjacently
MEs <- orderMEs(MEs)
# Verify sample alignment between eigengenes and metadata
trait_numeric <- trait_numeric[rownames(MEs), , drop = FALSE]
stopifnot(identical(rownames(MEs), rownames(trait_numeric)))
cat("Calculated", ncol(MEs), "Module Eigengenes across", nrow(MEs), "samples.\n")2. Pairwise Trait Correlation with Sample-Specific $n$ & Matrix-Wide BH FDR Correction
In plant stress trials, physiological measurements (such as Relative Water Content or photosynthetic efficiency) occasionally contain missing observations due to technical constraints. Standard matrix correlation routines silently fail or distort sample size $n$.
The rigorous workflow computes Pearson correlation and Student's t-test p-values only across pairwise complete observations, and crucially applies Benjamini-Hochberg (BH) multiple-testing correction across the entire $n_{ME} \times n_{\text{trait}}$ matrix:
nME <- ncol(MEs)
nTrait <- ncol(trait_numeric)
module_trait_cor <- matrix(NA_real_, nME, nTrait, dimnames = list(colnames(MEs), colnames(trait_numeric)))
module_trait_p <- module_trait_cor
module_trait_n <- matrix(0L, nME, nTrait, dimnames = dimnames(module_trait_cor))
# Compute pairwise correlation with explicit observation size verification
for (i in seq_len(nME)) {
for (j in seq_len(nTrait)) {
x <- MEs[, i]
y <- trait_numeric[, j]
ok <- is.finite(x) & is.finite(y)
n <- sum(ok)
module_trait_n[i, j] <- n
# Require at least 4 valid pairwise samples and non-zero variance
if (n >= 4 && length(unique(x[ok])) > 1 && length(unique(y[ok])) > 1) {
test <- cor.test(x[ok], y[ok], method = "pearson")
module_trait_cor[i, j] <- unname(test$estimate)
module_trait_p[i, j] <- test$p.value
}
}
}
# Apply Benjamini-Hochberg (BH) False Discovery Rate across FULL matrix
module_trait_q <- matrix(
p.adjust(as.vector(module_trait_p), method = "BH"),
nrow = nrow(module_trait_p),
ncol = ncol(module_trait_p),
dimnames = dimnames(module_trait_p)
)
# Screening Rule: Identify modules significantly associated with Drought (|r| >= 0.60, q < 0.05)
trait_name <- "is_drought"
stopifnot(trait_name %in% colnames(module_trait_cor))
candidate_MEs <- rownames(module_trait_cor)[
abs(module_trait_cor[, trait_name]) >= 0.60 &
module_trait_q[, trait_name] < 0.05
]
cat("\nCandidate Stress Modules identified:\n")
print(candidate_MEs)
if (length(candidate_MEs) == 0) {
message("⚠ No module passed the exploratory screening threshold (|r| >= 0.6, q < 0.05).")
}3. Publication-Ready Labeled Heatmap & Multi-Factor Linear Modeling
Visualize the correlation matrix using WGCNA's labeledHeatmap(), displaying the Pearson $r$ coefficient, FDR-adjusted $q$-value, and sample count $n$ directly inside each cell:
# Construct multi-line annotation text matrix for each cell
text_matrix <- matrix(
paste0(
signif(module_trait_cor, 2), "\n",
"q=", signif(module_trait_q, 2), "\n",
"n=", module_trait_n
),
nrow = nrow(module_trait_cor),
dimnames = dimnames(module_trait_cor)
)
# Generate publication-grade heatmap
labeledHeatmap(
Matrix = module_trait_cor,
xLabels = colnames(trait_numeric),
yLabels = rownames(module_trait_cor),
ySymbols = rownames(module_trait_cor),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = text_matrix,
setStdMargins = FALSE,
cex.text = 0.68,
zlim = c(-1, 1),
main = "Module–Trait Correlations with BH-Adjusted q-Values"
)In agricultural plant trials with treatment-by-time dynamics (e.g., Drought vs Control measured at 0h, 6h, 24h, 72h), simple bivariate correlation conflates time progression with stress response. Fit general linear models directly on the module eigengene:
# Model module eigengene dynamics over experimental factors
lm_fit <- lm(MEs$MEblue ~ treatment * timepoint, data = sample_meta)
anova_table <- anova(lm_fit)
print(anova_table)Common Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.