WGCNA로 유전자 공발현 네트워크 분석하기 — 모듈부터 허브 유전자까지 (R 실전)
TL;DR — 차등발현(differential expression) 분석이 유전자를 하나씩 본다면, WGCNA (weighted gene co-expression network analysis)는 함께 켜지고 꺼지는 유전자 무리를 모듈(module)로 묶어 시스템 수준에서 봅니다. 워크플로는 한 줄로 요약됩니다 — vst 발현행렬 → soft-thresholding power 선택 → TOM (topological overlap matrix) → 모듈 검출 → 모듈 고유유전자(module eigengene) → 형질(trait) 연관 → 허브(hub) 유전자.
이 글은 DESeq2의 vst (variance-stabilizing transformation) 출력을 입력으로, 복붙하면 그대로 도는 WGCNA 코드를 따라갑니다 —
goodSamplesGenes→pickSoftThreshold→blockwiseModules→moduleEigengenes→ module–trait 상관 →signedKME(허브) →exportNetworkToCytoscape. (코드는 WGCNA 공식 튜토리얼 기준입니다.)
🔗 관련 글 · 같이 보면 좋은 RNA-seq 분석: DESeq2로 차등발현 분석 · p값과 FDR (다중검정보정) · edgeR 차등발현 분석
0. 왜 WGCNA인가 — 1분 직관
차등발현 분석은 "처리군 대비 어떤 유전자가 올라갔나/내려갔나"를 유전자별로 묻습니다. 좋은 질문이지만, 수백 개 DEG 목록을 받아 들고 "그래서 이게 무슨 생물학이냐"에서 막히기 쉽습니다.
- WGCNA는 질문을 바꿉니다 — "어떤 유전자들이 샘플 전체에서 발현 패턴을 공유하는가?" 같이 오르내리는 유전자는 같은 경로(pathway)·같은 조절자(regulator) 아래 있을 가능성이 큽니다.
- 이렇게 묶인 무리가 모듈(module)입니다. 각 모듈을 하나의 대표 패턴(고유유전자)으로 요약하면, 수만 개 유전자 문제가 수십 개 모듈 문제로 줄어듭니다.
- 모듈을 형질(trait, 예: 처리 여부·시간·표현형)과 상관 내면, "이 형질과 연관된 유전자 집단"을 통째로 잡을 수 있습니다. 그 모듈 안에서 가장 중심에 있는 유전자가 허브(hub) — 후속 검증의 1순위 후보입니다.
1. 준비물 — 설치
WGCNA는 CRAN에 있지만 의존 패키지 일부(GO.db, impute, preprocessCore)가 Bioconductor에 있어, Bioconductor로 한 번에 까는 편이 깔끔합니다. 입력 데이터를 만들 DESeq2도 함께 설치합니다.
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("WGCNA", "GO.db", "impute", "preprocessCore", "DESeq2")) # 최초 1회
library(WGCNA)
library(DESeq2)
options(stringsAsFactors = FALSE) # WGCNA 권장 전역 설정2. 입력 데이터 — 왜 vst 발현행렬인가
WGCNA는 상관(correlation) 기반이라 분산이 안정화된 연속값이 필요합니다. 원시 카운트(raw count)를 그대로 넣으면 안 됩니다 — 평균이 큰 유전자에 분산이 쏠려(과대산포) 상관이 왜곡됩니다. DESeq2의 vst (또는 log2(x+1), rlog)로 변환한 값을 씁니다.
또 하나 중요한 제약 — WGCNA는 샘플 수가 충분해야 합니다(상관을 안정적으로 추정하려면 보통 ≥15~20개, 많을수록 좋음). 샘플이 5~6개면 모듈 자체가 불안정하니 차등발현 쪽이 낫습니다.
# 예시: DESeq2 카운트 → vst. 실제로는 본인의 RNA-seq counts/coldata로 대체
dds <- DESeqDataSetFromMatrix(countData = counts, # 유전자 × 샘플 정수 행렬
colData = coldata, # 샘플 메타데이터(형질 포함)
design = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ] # 거의 안 읽힌 유전자 제거
vsd <- vst(dds, blind = TRUE) # 분산안정화 변환
datExpr0 <- t(assay(vsd)) # ★ WGCNA는 '샘플 × 유전자'로 전치
dim(datExpr0) # 예: 24 samples × 18000 genes# 저분산 유전자 필터 — 변화 없는 유전자는 모듈 형성에 무의미하고 메모리만 잡아먹음
vars <- apply(datExpr0, 2, var)
datExpr0 <- datExpr0[, vars > quantile(vars, 0.25)] # 분산 하위 25% 컷(상황에 맞게 조정)유전자 × 샘플이지만, WGCNA는 샘플 × 유전자입니다. t()로 전치하는 걸 빠뜨리면 이후 모든 단계가 어긋납니다.3. goodSamplesGenes — 품질 점검(QC)
분석 전에 결측·무분산 유전자와 이상(outlier) 샘플을 솎아 냅니다. goodSamplesGenes()가 자동 점검을 해 주고, 샘플 군집(clustering)으로 눈에 띄는 이상치를 봅니다.
gsg <- goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK # TRUE면 통과
if (!gsg$allOK) { # 문제가 있으면 불량 행/열 제거
datExpr0 <- datExpr0[gsg$goodSamples, gsg$goodGenes]
}
# 샘플 군집으로 이상치(outlier) 확인 — 혼자 멀리 떨어진 가지가 있으면 의심
sampleTree <- hclust(dist(datExpr0), method = "average")
plot(sampleTree, main = "Sample clustering — outlier check", cex = 0.7)
abline(h = 120, col = "red") # 이 높이 위로 튀는 샘플을 잘라낼 수 있음
datExpr <- datExpr0 # 통과한 행렬을 최종 datExpr로 Flagging genes and samples with too many missing values...
..step 1
allOK: TRUE군집도에서 한 샘플만 유난히 멀리 떨어져 있으면 cutreeStatic()으로 그 가지를 제거합니다. 이상치 하나가 전체 상관 구조를 흔들 수 있어, soft power 선택 전에 처리하는 게 중요합니다.
4. pickSoftThreshold — soft-thresholding power 선택
WGCNA의 핵심 아이디어 — 상관을 그대로 쓰지 않고 거듭제곱(power, β)으로 올려 약한 상관은 더 약하게, 강한 상관은 살립니다. 이렇게 하면 네트워크가 척도 없는 위상(scale-free topology)에 가까워집니다(소수 허브가 많은 연결을 갖는, 실제 생물학적 네트워크의 특성). 적절한 β를 고르는 게 pickSoftThreshold()입니다.
powers <- c(1:10, seq(12, 20, by = 2)) # 후보 power들
sft <- pickSoftThreshold(datExpr,
powerVector = powers,
networkType = "signed", # 부호 있는 네트워크(아래 설명)
verbose = 5)
sft$powerEstimate # R^2 >= 0.85 를 처음 넘는 power
head(sft$fitIndices) # power별 적합 지표 Power SFT.R.sq slope mean.k. median.k. max.k.
1 1 0.020 1.85 4200 4100 6800
...
6 12 0.870 -1.42 180 120 640 <- R^2 0.85 돌파- SFT.R.sq ( R² )이 0.85~0.90을 넘는 가장 낮은 power를 고릅니다. R²가 높을수록 척도 없는 위상에 잘 맞는다는 뜻입니다.
- 동시에 mean.k. (평균 연결도, k )가 너무 낮아지지 않아야 합니다 — 너무 높은 power는 네트워크를 과하게 희소하게 만듭니다.
pickSoftThreshold가 적당한 power를 못 찾으면(어떤 power도 R²=0.85 미달), 표본 수 기준 권장값을 쓰는 게 공식 FAQ의 조언입니다.
| 샘플 수 | unsigned / signed hybrid | signed |
|---|---|---|
| < 20 | 9 | 18 |
| 20–30 | 8 | 16 |
| 30–40 | 7 | 14 |
| > 40 | 6 | 12 |
5. blockwiseModules — 네트워크 + 모듈 한 번에
이제 인접행렬(adjacency) → TOM → 군집 → 모듈 검출을 한 번에 처리합니다. 한 줄짜리 래퍼 blockwiseModules()가 표준입니다. (단계를 직접 보고 싶으면 아래 "단계별" 박스 참고.)
net <- blockwiseModules(
datExpr,
power = 12, # 4단계에서 고른 soft power
networkType = "signed", # pickSoftThreshold와 동일하게
TOMType = "signed",
minModuleSize = 30, # 모듈 최소 유전자 수
mergeCutHeight = 0.25, # 비슷한 모듈 병합 기준(고유유전자 상관 0.75 이상)
numericLabels = TRUE, # 모듈을 숫자로(나중에 색으로 변환)
saveTOMs = FALSE,
maxBlockSize = 20000, # 메모리 한도 — 유전자 수보다 크게
verbose = 3
)
table(net$colors) # 모듈별 유전자 수 (0 = 어디에도 안 묶인 회색)
moduleColors <- labels2colors(net$colors) # 숫자 라벨 → 색 이름 0 1 2 3 4 5 6
410 3200 2100 1500 980 620 410 <- 모듈 0(grey)은 미할당 유전자각 유전자는 하나의 모듈(색)에 배정됩니다. 모듈 0(grey)은 어떤 모듈에도 들지 못한 유전자라 보통 해석에서 제외합니다. 모듈 덴드로그램은 plotDendroAndColors()로 봅니다.
plotDendroAndColors(net$dendrograms[[1]],
moduleColors[net$blockGenes[[1]]],
"Module", dendroLabels = FALSE,
main = "Gene dendrogram and module colors")blockwiseModules 한 줄은 사실 이 과정입니다.
adj <- adjacency(datExpr, power = 12, type = "signed") # 상관^power
TOM <- TOMsimilarity(adj, TOMType = "signed") # 위상중첩(TOM)
diss <- 1 - TOM # 거리 = 1 - TOM
tree <- hclust(as.dist(diss), method = "average") # 계층 군집
mods <- cutreeDynamic(tree, distM = diss, # 동적 가지치기
deepSplit = 2, minClusterSize = 30)
merged <- mergeCloseModules(datExpr, labels2colors(mods), # 비슷한 모듈 병합
cutHeight = 0.25)TOM은 단순 상관을 넘어 "두 유전자가 이웃을 얼마나 공유하는가"까지 반영해 모듈 경계를 또렷하게 만듭니다.
6. moduleEigengenes — 모듈 고유유전자(ME)
각 모듈을 하나의 대표 패턴으로 요약한 게 모듈 고유유전자(module eigengene, ME) — 그 모듈 발현행렬의 제1주성분(first principal component)입니다. 이제 "유전자 수만 개"가 아니라 "ME 수십 개"로 형질과 연관을 따집니다.
MEs0 <- moduleEigengenes(datExpr, colors = moduleColors)$eigengenes
MEs <- orderMEs(MEs0) # 비슷한 ME끼리 정렬(히트맵 보기 좋게)
head(MEs[, 1:4]) # 샘플 × 모듈ME 행렬 MEturquoise MEblue MEbrown MEgreen
Sample1 0.18 -0.22 0.05 0.31
Sample2 -0.24 0.19 -0.11 -0.08
...ME는 샘플마다 그 모듈이 얼마나 켜졌는지를 나타내는 한 개의 값입니다. 모듈 전체 발현을 한 숫자로 압축했으니, 형질과의 상관을 바로 계산할 수 있습니다.
7. module–trait 상관 — 형질과 모듈 연관
WGCNA의 하이라이트 — 각 ME를 형질(trait)과 상관 내어, "어떤 모듈이 어떤 표현형과 묶이는가"를 한 장의 히트맵으로 봅니다. 형질은 coldata의 수치형 변수(처리 0/1, 시간, 측정값 등)입니다.
traits <- as.data.frame(colData(vsd)[, c("condition_num", "time", "score")])
nSamples <- nrow(datExpr)
moduleTraitCor <- cor(MEs, traits, use = "p") # ME ↔ 형질 상관
moduleTraitPval <- corPvalueStudent(moduleTraitCor, nSamples) # 상관의 p값
# 셀에 "상관(p값)" 텍스트를 얹은 라벨 행렬
textMatrix <- paste0(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPval, 1), ")")
dim(textMatrix) <- dim(moduleTraitCor)
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(traits), yLabels = names(MEs),
ySymbols = names(MEs), colorLabels = FALSE,
colors = blueWhiteRed(50), textMatrix = textMatrix,
setStdMargins = FALSE, cex.text = 0.7,
main = "Module–trait relationships") condition_num time score
MEturquoise 0.82 0.21 0.74 <- 처리·score와 강한 양의 상관
(3e-06) (0.32) (4e-05)
MEblue -0.79 -0.15 -0.68
(5e-06) (0.48) (3e-04)빨강(양)·파랑(음)이 진하고 p값이 작은 셀이 형질과 강하게 연관된 모듈입니다. 위 예라면 turquoise 모듈이 처리·score와 강한 양의 상관 — 이 모듈을 파고들 가치가 있습니다.
8. signedKME — 허브(hub) 유전자 찾기
관심 모듈을 정했으면(예: turquoise), 그 안에서 가장 중심에 있는 유전자 = 허브를 뽑습니다. 기준은 모듈 소속도(module membership, kME) — 각 유전자가 자기 모듈 ME와 얼마나 상관되는가입니다. signedKME()로 계산합니다.
kME <- signedKME(datExpr, MEs) # 모든 유전자 × 모든 모듈 kME
# 열 이름 예: kMEturquoise, kMEblue ...
# turquoise 모듈 유전자만 골라 kME 높은 순으로
inMod <- moduleColors == "turquoise"
kmeTur <- kME[inMod, "kMEturquoise"]
hubs <- sort(kmeTur, decreasing = TRUE)
head(hubs, 10) # 상위 10개 허브 후보 GeneA GeneB GeneC GeneD GeneE
0.97 0.95 0.94 0.93 0.92 <- kME 높을수록 모듈 중심(허브)- kME (module membership)가 1에 가까울수록 그 유전자는 모듈을 대표합니다 — 허브 후보입니다.
- 한 줄로 각 모듈 최고 허브만 뽑으려면
chooseTopHubInEachModule(datExpr, moduleColors)도 있습니다. - 진짜 후보를 좁히려면 kME 높음 + 형질과의 유전자 유의성(gene significance) 높음을 함께 보세요(
cor(datExpr, trait)).
9. exportNetworkToCytoscape — 네트워크 내보내기
허브와 모듈 구조를 Cytoscape에서 시각화하려면, 모듈 내부 TOM을 엣지(edge)·노드(node) 표로 내보냅니다.
# 관심 모듈 유전자들의 TOM만 다시 계산(메모리 절약)
TOM <- TOMsimilarityFromExpr(datExpr, power = 12, networkType = "signed")
modGenes <- colnames(datExpr)[inMod]
modTOM <- TOM[inMod, inMod]
dimnames(modTOM) <- list(modGenes, modGenes)
exportNetworkToCytoscape(
modTOM,
edgeFile = "edges_turquoise.txt",
nodeFile = "nodes_turquoise.txt",
weighted = TRUE,
threshold = 0.10, # 약한 엣지는 잘라 그래프 단순화
nodeNames = modGenes,
nodeAttr = moduleColors[inMod]
)내보낸 edges_*.txt/nodes_*.txt를 Cytoscape에서 열면, 허브가 중앙에 모인 모듈 네트워크를 그릴 수 있습니다.
10. WGCNA는 언제 쓰나 — 공발현 모듈 vs 차등발현
둘은 경쟁이 아니라 보완입니다. 차등발현으로 "무엇이 변했나"를, WGCNA로 "그 변화가 어떤 시스템에 속하나"를 봅니다.
| 항목 | WGCNA (공발현 모듈) | 차등발현(DESeq2 / edgeR) |
|---|---|---|
| 보는 단위 | 유전자 무리(모듈) | 개별 유전자 |
| 핵심 질문 | 함께 움직이는 유전자는? | 그룹 간 변한 유전자는? |
| 그룹 라벨 | 불필요(비지도, 형질은 사후 연결) | 필수(처리 vs 대조) |
| 입력 | vst / log2 발현행렬 | 원시 카운트 |
| 샘플 수 | 많이 필요(≥15~20) | 적어도 됨(그룹당 ≥3) |
| 출력 | 모듈·ME·허브·형질 연관 | DEG 목록(logFC·FDR) |
| 강점 | 시스템·경로 수준 해석, 허브 발굴 | 명확한 효과 크기·유의성 |
| 한계 | 표본 적으면 불안정, 인과 아님 | 개별 유전자라 맥락 약함 |
- WGCNA를 쓸 때 — 샘플이 많고(시계열·코호트·다조건), 표현형과 연관된 유전자 모듈/허브를 찾고 싶을 때.
- 차등발현을 쓸 때 — 명확한 2~몇 그룹 비교에서 개별 DEG와 효과 크기가 목적일 때.
- 함께 쓰면 최고 — DEG로 후보를 잡고, WGCNA 모듈로 그 DEG들이 어떤 공발현 시스템에 속하는지 해석.
11. 자주 발생하는 에러 & 해결
- 모듈이 한두 개로 뭉치거나 회색(grey) 천지 — 표본이 너무 적거나(<15), soft power가 부적절. 표본을 늘리거나 FAQ 권장 power로.
pickSoftThreshold가 power를 못 찾음(R² 안 올라감) — 배치 효과(batch effect)나 강한 이상치 신호가 흔한 원인. 군집도로 확인하고 ComBat 등으로 보정 후 재시도. 그래도 안 되면 표본 수 기준 권장 power 사용.- 원시 카운트를 넣음 — 상관이 평균에 끌려 모듈이 깨집니다. 반드시 vst / rlog / log2 변환값을.
Error: cannot allocate vector of size ...— 유전자가 너무 많음. 저분산 유전자 필터링을 강화하거나maxBlockSize를 메모리에 맞춰 조정(블록 분할).- 모듈–형질 상관이 다 약함 — 형질이 수치형이 아니거나(요인을
as.numeric으로), 모듈이 과병합됨(mergeCutHeight낮추기). - 방향(전치) 실수 —
t()를 빠뜨려유전자 × 샘플로 넣으면 전부 어긋남.dim(datExpr)가 샘플 × 유전자인지 확인.
12. 확장 — 모듈 → clusterProfiler 농축분석
모듈은 "함께 움직이는 유전자 목록"이라, 그대로 기능 농축분석(enrichment)의 입력이 됩니다. 관심 모듈 유전자를 clusterProfiler에 넘기면 어떤 GO·KEGG 경로가 모이는지 보입니다.
library(clusterProfiler)
library(org.Hs.eg.db)
modGenes <- colnames(datExpr)[moduleColors == "turquoise"] # 모듈 유전자(SYMBOL 가정)
ego <- enrichGO(gene = modGenes, OrgDb = org.Hs.eg.db,
keyType = "SYMBOL", ont = "BP", # 생물학적 과정
pAdjustMethod = "BH", qvalueCutoff = 0.05)
head(ego) # 모듈에 모인 경로여기까지 오면 흐름이 완성됩니다 — 공발현 모듈 → 형질 연관 → 허브 + 경로 농축. 모듈 하나하나가 하나의 가설이 됩니다.
• (입력 만들기) DESeq2로 차등발현 분석 — vst 발현행렬의 출처 · edgeR 차등발현 분석
• (QC·시각화) PCA로 샘플 QC — 이상치·배치 점검 · ComplexHeatmap 히트맵
• (다음 단계) clusterProfiler 농축분석 — 모듈 → GO/KEGG 경로 해석
References
- Langfelder, P., Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 9, 559.
- 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. (soft-thresholding·scale-free·TOM)
- WGCNA Tutorials (Langfelder & Horvath, UCLA). goodSamplesGenes·pickSoftThreshold·blockwiseModules·moduleEigengenes·labeledHeatmap·signedKME·exportNetworkToCytoscape.
- WGCNA FAQ (Langfelder). 표본 수 기준 권장 soft power 표 · signed vs unsigned 네트워크 · scale-free 적합 실패 시 대처.
- Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 (vst 변환). Genome Biology, 15, 550.
- Yu, G., et al. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5), 284–287.
Pipette & Pipeline · A bio portfolio journal
'Code Bench > R Snippets' 카테고리의 다른 글
| limma-voom으로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2·edgeR과 무엇이 다른가 (R 실전) (0) | 2026.06.17 |
|---|---|
| GEOquery로 공개 발현 데이터 받아서 분석하기 — GEO에서 내 PC로, 재분석까지 (R 실전) (0) | 2026.06.17 |
| edgeR로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2와 무엇이 다른가 (R 실전) (0) | 2026.06.15 |
| 주성분분석(PCA)으로 RNA-seq 샘플 QC하기 — vst·plotPCA·prcomp·배치효과 (R 실전) (0) | 2026.06.11 |
| ComplexHeatmap으로 발현 히트맵 그리기 — z-score부터 주석·분할까지 (R 실전) (0) | 2026.06.10 |
