Chapter 7: Intramodular Hub Gene Prioritization & Module Membership (kME) Analysis
Chapter 7: Intramodular Hub Gene Prioritization & Module Membership (kME) Analysis
Quantifying module membership (kME) as correlation between gene expression and module eigengenes; calculating Student t-test p-values and BH-adjusted q-values; ranking central hubs; and annotating with Arabidopsis TAIR locus descriptions.
1. The Concept of Module Membership (kME)
Within a detected co-expression module, not all transcripts exert equal topological influence. Intramodular hub genes occupy the geometric center of the module cluster, exhibiting the highest number of co-expression connections to fellow module members.
While intramodular connectivity ($k_{in}$) can be calculated by summing adjacency weights ($k_{in, i} = \sum_{j \in \text{module}} a_{ij}$), Horvath and Dong demonstrated that Module Membership ($k_{ME}$) provides an equivalent, computationally superior, and statistically interpretable continuous metric:
Where $x_i$ is the expression profile of gene $i$ across all samples, and $ME_q$ is the eigengene of module $q$. If $k_{ME, i} \approx +1.0$, the individual gene's transcriptional kinetics mirror the overarching trajectory of the entire module. A high $k_{ME}$ isolates master coordinators and potential transcription factor drivers.
2. Computing Signed kME, Statistical Significance & Hub Ranking
WGCNA's signedKME() calculates module membership across all network genes simultaneously. To identify top candidates, isolate the module displaying the strongest trait correlation, compute Student's t-test p-values, and apply Benjamini-Hochberg FDR correction:
# Select candidate module with highest significant drought correlation
if (length(candidate_MEs) == 0) {
stop("No candidate module passed screening. Select a module based on documented trait correlation.")
}
target_ME <- candidate_MEs[1] # e.g., "MEblue"
stopifnot(target_ME %in% colnames(MEs))
# Extract module color (strip leading 'ME')
target_color <- sub("^ME", "", target_ME)
stopifnot(target_color %in% unique(module_colors))
# Isolate all genes assigned to target module
target_genes <- names(module_colors)[module_colors == target_color]
cat("Target Module:", target_color, "| Total Genes:", length(target_genes), "\n")
# Calculate signed kME across all genes
kme_res <- signedKME(
expr_mat,
datME = MEs,
corFnc = "cor",
corOptions = "use = 'p'"
)
# Extract kME vector for our target module
kme_col <- paste0("kME", target_color)
stopifnot(kme_col %in% colnames(kme_res))
kme <- kme_res[target_genes, kme_col]
# Compute Student's asymptotic p-value and BH FDR q-value
kme_p <- corPvalueStudent(kme, nSamples = nrow(expr_mat))
hub_df <- data.frame(
gene = names(kme),
kME = as.numeric(kme),
p_kME = as.numeric(kme_p),
q_kME = p.adjust(kme_p, method = "BH"),
stringsAsFactors = FALSE
)
# Rank genes by descending absolute kME (highest intramodular hubness first)
hub_df <- hub_df[order(-abs(hub_df$kME), hub_df$q_kME), , drop = FALSE]
# Annotate with Arabidopsis TAIR Gene Names and Functional Descriptions
hub_df$symbol <- mapIds(
org.At.tair.db,
keys = hub_df$gene,
keytype = "TAIR",
column = "SYMBOL",
multiVals = "first"
)
hub_df$description <- mapIds(
org.At.tair.db,
keys = hub_df$gene,
keytype = "TAIR",
column = "GENENAME",
multiVals = "first"
)
# Display Top 15 Intramodular Hub Regulators
cat("\n=== Top 15 Intramodular Hub Regulators in Module", target_color, "===\n")
print(head(hub_df[, c("gene", "symbol", "kME", "q_kME", "description")], 15))3. Prioritizing Master Regulators in Plant Stress Cascades
When inspecting top hubs in plant drought and abiotic stress modules, look for convergence between high intramodular connectivity ($|k_{ME}| \ge 0.80$) and direct Gene Significance ($|GS_{\text{drought}}| \ge 0.60$):
Key Plant Transcription Factor Families Operating as Network Hubs
- DREB / CBF (AP2/ERF Family): Dehydration-responsive element-binding factors (e.g., DREB2A /
AT5G05410) controlling ABA-independent osmotic stress responses. - bZIP / AREB / ABF Family: ABA-responsive element binding factors (e.g., AREB1/ABF2 /
AT1G45249) orchestrating stomatal closure and LEA protein accumulation. - WRKY Superfamily: Stress-responsive factors (e.g., WRKY40 /
AT1G80840) modulating the trade-off between growth, immune signaling, and abiotic drought survival. - NAC Domain Factors: Regulators of senescence, drought-induced leaf shedding, and reactive oxygen species (ROS) scavenging enzymes.
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.