WGCNA für Omics-Module: Von Koexpression zu klinisch interpretierbaren ClusternWGCNA for Omics Modules: From Co-expression to Clinically Interpretable Clusters

Abstract

WGCNA (Weighted Gene Co-expression Network Analysis) identifiziert Module ko-exprimierter Gene in großen Expressionsdatensätzen und korreliert diese mit klinischen oder phänotypischen Merkmalen. Anders als diff. Expressionsanalysen, die einzelne Gene testen, betrachtet WGCNA das gesamte Transkriptom als Netzwerk – und findet Gruppen von Genen, die gemeinsam reguliert werden. Für Teams in der translationalen Forschung liefert WGCNA damit einen systembiologischen Blick auf Krankheitsmechanismen, der über Gen-Listen hinausgeht.

Typisches Projektszenario

Ein Onkologie-Labor analysiert RNA-Seq-Daten von 120 kolorektalen Karzinom-Patienten mit unterschiedlichen Tumorstadien (I–IV), Mikrosatelliten-Instabilitätsstatus (MSI-H vs. MSS) und Therapieansprechraten (Responder vs. Non-Responder). Die normalisierte Expressionsmatrix (15.000 variabelste Gene × 120 Proben) wird in WGCNA geladen. Ziel: Module identifizieren, die signifikant mit dem Therapieansprechen korrelieren, und deren Hub-Gene als potenzielle therapeutische Targets oder Biomarker-Kandidaten priorisieren.

Welches Problem löst WGCNA?

• Von Genlisten zu biologischen Modulen: Ein DE-Experiment liefert oft 1.000–5.000 signifikante Gene. WGCNA gruppiert diese in Module, die häufig biologischen Pathways entsprechen – zum Beispiel ein Modul für Immunantwort, eines für Zellzyklus, eines für Metabolismus.
• Trait-Assoziation: Jedes Modul wird über sein Eigengene (erste Hauptkomponente) mit klinischen Variablen korreliert. So findet man nicht nur, welche Gene verändert sind, sondern welche Gene gemeinsam mit dem Krankheitsphänotyp variieren.
• Hub-Gen-Identifikation: Innerhalb jedes Moduls werden die am stärksten vernetzten Gene (Hub-Gene) identifiziert – diese sind oft biologisch die relevantesten und gute Kandidaten für funktionelle Validierung.

Warum Teams WGCNA einsetzen

1. Hypothesenfreie Entdeckung: WGCNA findet biologisch relevante Gengruppen ohne Vorannahmen über Pathways.
2. Dimensionsreduktion: Statt 15.000 Gene einzeln zu testen, werden 20–50 Module als biologisch interpretierbare Einheiten analysiert – das reduziert das Multiple-Testing-Problem.
3. Netzwerk-Topologie: Die gewichtete Netzwerkstruktur enthült Information über die Stärke der Ko-Expression, die bei einfacher Korrelation verloren geht.
4. Scale-free Topologie: WGCNA erzwingt eine skalenfreie Netzwerk-Topologie, die der Architektur biologischer Netzwerke entspricht – wenige hoch vernetzte Hub-Gene, viele periphere Gene.
5. Reproduzierbarkeit: Module können zwischen Datensätzen verglichen werden (modulePreservation()), was die biologische Relevanz validiert.

Netzwerkkonstruktion

WGCNA konstruiert das Ko-Expressionsnetzwerk in mehreren Schritten, die den mathematischen Kern des Verfahrens ausmachen. Zunächst wird die paarweise Pearson-Korrelation aller Gene berechnet. Dann wird die Korrelationsmatrix in eine Adjazenzmatrix transformiert: a_ij = |cor(x_i, x_j)|^β, wobei β der Soft-Thresholding-Power ist. Der pickSoftThreshold()-Algorithmus wählt β so, dass das resultierende Netzwerk approximativ skalenfrei ist – typisch sind Werte von β = 6–14 für unsigned und β = 12–20 für signed Netzwerke. Ein R² > 0.8 für den Scale-Free-Fit gilt als ausreichend.

Das Soft-Thresholding hat gegenüber einem Hard-Threshold (Korrelation > Schwelle = 1, sonst = 0) den Vorteil, dass die kontinuierliche Natur der Ko-Expression erhalten bleibt. Gene mit starker Korrelation erhalten hohe Adjazenzwerte, Gene mit schwacher Korrelation niedrige – aber keine werden willkürlich auf 0 gesetzt.

Die Adjazenzmatrix wird in eine topologische Überlappungsmatrix (TOM) umgerechnet: TOM_ij = (∑_u a_iua_uj + a_ij) / (min(k_i, k_j) + 1 – a_ij). Die TOM misst nicht nur die direkte Korrelation, sondern berücksichtigt auch gemeinsame Nachbarn – Gene, die ähnliche Konnektivitätsprofile haben, erhalten hohe TOM-Werte, auch wenn ihre direkte Korrelation moderat ist. Die Module werden dann durch hierarchisches Clustering auf 1–TOM (Dissimilaritätsdistanz) identifiziert.

R-Code: WGCNA-Pipeline

library(WGCNA)
options(stringsAsFactors = FALSE)
allowWGCNAThreads(nThreads = 4)

# Daten vorbereiten
expr <- read.csv("crc_vst_top15k.csv", row.names = 1)
datExpr <- t(expr)  # Proben x Gene
traits <- read.csv("crc_clinical.csv", row.names = 1)

# Soft-Threshold bestimmen
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector = powers,
                          networkType = "signed hybrid")
beta <- sft$powerEstimate  # typisch: 12-14

# Netzwerk + Module in einem Schritt
net <- blockwiseModules(
  datExpr, power = beta,
  networkType = "signed hybrid",
  TOMType = "signed",
  minModuleSize = 30,
  mergeCutHeight = 0.25,
  numericLabels = TRUE,
  saveTOMs = FALSE,
  verbose = 3
)

# Modul-Trait-Korrelation
MEs <- net$MEs
moduleTraitCor <- cor(MEs, traits, use = "p")
moduleTraitPval <- corPvalueStudent(moduleTraitCor, nrow(datExpr))

# Visualisierung der Modul-Trait-Beziehungen
textMatrix <- paste0(signif(moduleTraitCor, 2), "
(",
                     signif(moduleTraitPval, 1), ")")
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = colnames(traits),
               yLabels = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               main = "Module-Trait Relationships")

Beispielausgabe

# pickSoftThreshold: power = 12, R^2 = 0.87
# blockwiseModules: 23 modules detected
# Module sizes: 31 - 1,247 genes
# Module 0 (grey, unassigned): 892 genes

# Top module-trait correlations (response):
# ME_turquoise:  r = -0.54, p = 1.2e-10  (immune)
# ME_blue:       r =  0.48, p = 3.4e-08  (cell cycle)
# ME_brown:      r = -0.41, p = 2.1e-06  (metabolism)
# ME_green:      r =  0.38, p = 8.7e-06  (EMT)

Diagnostische Plots

Vergleich mit Alternativen

Merkmal	WGCNA	MEGENA	CEMiTool
Netzwerktyp	Gewichtet, skalenfrei	Planar gefiltert	Gewichtet (WGCNA-Kern)
Moduldetektion	Hierarchisches Clustering + dynamischer Baumschnitt	Multi-Skalen-Clustering	Automatisch (WGCNA-basiert)
Skalierbarkeit	~20k Gene, ~500 Proben	~15k Gene, ~300 Proben	~20k Gene (automatisiert)
Trait-Korrelation	Eigengene-basiert	Hub-Score-basiert	Eigengene + ORA
Benutzerfreundlichkeit	Mittel (viele Parameter)	Komplex	Hoch (wenige Parameter)
Modulerhaltung	modulePreservation()	Extern	Integriert

Statistische Vertiefung

Die Wahl zwischen unsigned, signed und signed hybrid Netzwerken beeinflusst die biologische Interpretation fundamental. Unsigned Netzwerke verwenden |cor|^β und gruppieren Gene, die stark positiv oder negativ korreliert sind, in dasselbe Modul. Signed Netzwerke verwenden (0.5 + 0.5·cor)^β und trennen positiv und negativ korrelierte Gene in verschiedene Module. Signed hybrid verwendet |cor|^β für positive Korrelationen und 0 für negative – ein Kompromiss, der biologisch oft am sinnvollsten ist, weil ko-regulierte Gene typischerweise positiv korrelieren.

Die Modul-Eigengene sind die erste Hauptkomponente der Expressionsdaten aller Gene in einem Modul. Sie erklären typisch 40–70% der Modulvariation und dienen als Zusammenfassung des gesamten Moduls. Die Korrelation eines Gen-Expressionsprofils mit dem Eigengene heißt Module Membership (kME) – Gene mit kME > 0.8 sind Hub-Gene und treiben die Modul-Biologie. Die Gene Significance (GS) misst die Korrelation jedes Gens mit einem externen Trait. Gene mit hoher kME und hoher GS sind die vielversprechendsten Kandidaten für experimentelle Validierung.

Die modulePreservation()-Funktion testet, ob die identifizierten Module in einem unabhängigen Datensatz replizierbar sind. Das Z_summary-Maß kombiniert Dichte- und Konnektivitätsmetriken: Z_summary > 10 zeigt starke, Z_summary 2–10 moderate Erhaltung. Module, die nicht erhalten sind, können datensatzspezifische Artefakte oder Batch-Effekte widerspiegeln.

Zitationen

1. Langfelder P, Horvath S (2008). “WGCNA: an R package for weighted correlation network analysis.” BMC Bioinformatics, 9, 559.
2. Zhang B, Horvath S (2005). “A general framework for weighted gene co-expression network analysis.” Statistical Applications in Genetics and Molecular Biology, 4, Article 17.
3. Langfelder P, Horvath S (2012). “Fast R functions for robust correlations and hierarchical clustering.” Journal of Statistical Software, 46(11).

Fazit

WGCNA verwandelt Gen-Listen in biologisch interpretierbare Netzwerk-Module und ist damit das Mittel der Wahl für explorative Analysen großer Kohorten. Die Stärke liegt in der Fähigkeit, tausende Gene zu funktionellen Einheiten zusammenzufassen und diese direkt mit klinischen Phänotypen zu verbinden. Limitierungen: (1) Die Ergebnisse hängen sensitiv von der Parameterw ahl ab (β, minModuleSize, mergeCutHeight). (2) WGCNA setzt normalverteilte, homoskedastische Daten voraus – eine Varianzstabilisierung (VST/rlog) ist Pflicht. (3) Netzwerke sind korrelationsbasiert und implizieren keine Kausalität. Aber als erster Schritt zur Hypothesengenerierung in Omics-Studien ist WGCNA unersetzlich.

Dokumentation

• WGCNA-Homepage
• WGCNA-Tutorials

library(WGCNA)
options(stringsAsFactors = FALSE)
allowWGCNAThreads(nThreads = 4)

# Prepare data
expr <- read.csv("crc_vst_top15k.csv", row.names = 1)
datExpr <- t(expr)  # samples x genes
traits <- read.csv("crc_clinical.csv", row.names = 1)

# Determine soft threshold
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector = powers,
                          networkType = "signed hybrid")
beta <- sft$powerEstimate  # typically: 12-14

# Network + modules in one step
net <- blockwiseModules(
  datExpr, power = beta,
  networkType = "signed hybrid",
  TOMType = "signed",
  minModuleSize = 30,
  mergeCutHeight = 0.25,
  numericLabels = TRUE,
  saveTOMs = FALSE,
  verbose = 3
)

# Module-trait correlation
MEs <- net$MEs
moduleTraitCor <- cor(MEs, traits, use = "p")
moduleTraitPval <- corPvalueStudent(moduleTraitCor, nrow(datExpr))

# Visualization of module-trait relationships
textMatrix <- paste0(signif(moduleTraitCor, 2), "
(",
                     signif(moduleTraitPval, 1), ")")
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = colnames(traits),
               yLabels = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               main = "Module-Trait Relationships")

# pickSoftThreshold: power = 12, R^2 = 0.87
# blockwiseModules: 23 modules detected
# Module sizes: 31 - 1,247 genes
# Module 0 (grey, unassigned): 892 genes

# Top module-trait correlations (response):
# ME_turquoise:  r = -0.54, p = 1.2e-10  (immune)
# ME_blue:       r =  0.48, p = 3.4e-08  (cell cycle)
# ME_brown:      r = -0.41, p = 2.1e-06  (metabolism)
# ME_green:      r =  0.38, p = 8.7e-06  (EMT)