마이크로어레이 분석의 확장: 배치 효과 보정, 교차 플랫폼 통합, 메타분석, GSVA/ssGSEA, WGCNA를 위한 실전 템플릿을 제공합니다.

---

A) 배치 효과 보정 (ComBat / removeBatchEffect / SVA)

정규화(log2) 후, 생물학적 설계(group)를 보존하면서 batch를 보정합니다.

A-1) 배치 확인 (PCA)

Copy
# expr: (gene × sample) log2 행렬, meta: data.frame(sample, group, batch)
mat <- expr
pc  <- prcomp(t(mat), scale.=TRUE)
plot(pc$x[,1], pc$x[,2], pch=19, col=as.integer(meta$batch)+1,
     xlab="PC1", ylab="PC2"); legend("topright", legend=levels(meta$batch),
     col=2:(length(levels(meta$batch))+1), pch=19)

A-2) ComBat (sva 패키지)

생물학적 변수(group)는 모델로 포함하여 보존합니다.

Copy
if (!requireNamespace("sva", quietly=TRUE)) BiocManager::install("sva")
library(sva)

batch  <- factor(meta$batch)
group  <- factor(meta$group)
mod    <- model.matrix(~ group)          # 생물학적 효과 유지
expr_cb <- ComBat(dat=as.matrix(expr), batch=batch, mod=mod, par.prior=TRUE)

# PCA 다시 확인
pc2 <- prcomp(t(expr_cb), scale.=TRUE)
plot(pc2$x[,1], pc2$x[,2], pch=19, col=as.integer(group)+1,
     xlab="PC1", ylab="PC2")

A-3) removeBatchEffect (시각화용 간단 보정)

분석 모델 적합은 원본 expr + 디자인에 batch 포함을 권장. removeBatchEffect는 주로 시각화 용도.

Copy
library(limma)
expr_vis <- removeBatchEffect(expr, batch = batch, design = model.matrix(~ group))
# 히트맵/PCA 등 시각화에 expr_vis 사용

A-4) SVA (숨은 요인)

Copy
library(sva); library(limma)
design <- model.matrix(~ group)
# null 모델(생물학 효과 제외)과 full 모델 지정
design0 <- model.matrix(~ 1)
svobj <- sva(as.matrix(expr), mod=design, mod0=design0)   # sv: surrogate variables
design_sv <- cbind(design, svobj$sv)
fit  <- lmFit(expr, design_sv)
# ... 평소대로 contrasts.fit/eBayes 진행

B) 교차 플랫폼 통합 (Affy/Agilent/Illumina)

플랫폼별 표현행렬을 gene symbol 기준으로 정렬 후 병합하고, study를 배치로 처리합니다.

Copy
# expr_list: 리스트 형태로 각 스터디(혹은 플랫폼) 표현행렬(log2, gene × sample)
# 공통 유전자 교집합을 취해 column-bind
genes <- Reduce(intersect, lapply(expr_list, rownames))
mats  <- lapply(expr_list, function(x) x[genes, , drop=FALSE])
X     <- do.call(cbind, mats)

# 간단한 스케일링/정렬(선택)
Xn <- limma::normalizeBetweenArrays(X, method="quantile")

# 메타 정보
study <- factor(unlist(mapply(function(m, nm) rep(nm, ncol(m)),
                              expr_list, names(expr_list))))
group <- factor(meta$group[colnames(Xn)])  # 미리 준비된 그룹 라벨

# 배치 포함한 limma 모델
design <- model.matrix(~ 0 + group + study)
colnames(design)[1:length(levels(group))] <- levels(group)
fit <- lmFit(Xn, design)
# 이후 contrasts.fit/eBayes ...

C) 메타분석 (metafor: 효과크기 결합, 랜덤효과)

각 스터디에서 limma로 효과크기(logFC)와 t 통계를 구한 뒤, SE = |logFC / t|로 근사해 결합합니다.

Copy
if (!requireNamespace("metafor", quietly=TRUE)) install.packages("metafor")
library(metafor); library(dplyr); library(purrr)

# studies: 리스트(각각 topTable 결과 with rownames=gene, logFC, t, P.Value 등)
combine_gene <- function(gene){
  per <- lapply(studies, function(tt){
    if (!gene %in% rownames(tt)) return(NULL)
    lf <- tt[gene, "logFC"]; tval <- tt[gene, "t"]
    if (is.na(lf) || is.na(tval) || tval==0) return(NULL)
    se <- abs(lf / tval)
    data.frame(yi=lf, sei=se)
  }) %>% bind_rows()
  if (nrow(per) < 2) return(NULL)
  r <- rma(yi=per$yi, sei=per$sei, method="REML")
  data.frame(gene=gene, eff=r$b[[1]], se=sqrt(r$tau2 + mean(per$sei^2)),
             z=r$zval, p=r$pval)
}

genes_all <- Reduce(union, lapply(studies, rownames))
meta_res <- map_dfr(genes_all, combine_gene)
meta_res$padj <- p.adjust(meta_res$p, method="BH")
meta_res <- meta_res[order(meta_res$padj), ]
head(meta_res)
write.csv(meta_res, "meta_random_effects.csv", row.names=FALSE)

D) GSVA / ssGSEA (샘플별 경로 점수)

발현행렬과 유전자집합(리스트)을 입력하여 샘플별 경로 스코어를 계산합니다.

Copy
if (!requireNamespace("GSVA", quietly=TRUE)) BiocManager::install("GSVA")
if (!requireNamespace("msigdbr", quietly=TRUE)) install.packages("msigdbr")
library(GSVA); library(msigdbr); library(dplyr)

# expr: (gene × sample) log2
msig <- msigdbr(species="Homo sapiens", category="H") %>%
        dplyr::select(gs_name, gene_symbol)
gsets <- split(msig$gene_symbol, msig$gs_name)

# GSVA (기본) 또는 ssGSEA(method="ssgsea")
gsva_mat <- gsva(as.matrix(expr), gsets, method="gsva", kcdf="Gaussian", parallel.sz=1)
# ssGSEA: gsva(..., method="ssgsea")

write.csv(gsva_mat, "GSVA_scores.csv")

E) WGCNA (공발현 모듈)

유전자 간 상관 기반 네트워크로 모듈을 정의하고, 모듈-표현형 상관을 확인합니다.

Copy
if (!requireNamespace("WGCNA", quietly=TRUE)) install.packages("WGCNA")
library(WGCNA)
options(stringsAsFactors=FALSE); allowWGCNAThreads()

datExpr <- t(expr)                   # WGCNA는 sample × gene 형태 권장
gsg <- goodSamplesGenes(datExpr, verbose=3); if (!gsg$allOK) datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes]

# 소프트 임계치 선택
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector=powers, verbose=5)
power <- sft$powerEstimate %||% 6      # 추정 실패 시 기본값 6

# 모듈 탐색
net <- blockwiseModules(datExpr, power=power,
                        TOMType="unsigned", minModuleSize=30,
                        reassignThreshold=0, mergeCutHeight=0.25,
                        numericLabels=TRUE, pamRespectsDendro=FALSE,
                        saveTOMs=FALSE, verbose=3)
moduleColors <- labels2colors(net$colors)
MEs <- net$MEs

# 모듈-표현형 상관
trait <- model.matrix(~ 0 + meta$group); colnames(trait) <- levels(meta$group)
modTraitCor <- cor(MEs, trait, use="p"); modTraitP <- corPvalueStudent(modTraitCor, nrow(datExpr))
write.csv(modTraitCor, "WGCNA_module_trait_cor.csv")

F) 재현성 & 세션 정보

Copy
set.seed(123)
sessionInfo()

요약

• 배치는 ComBat/SVA로 처리하고, 시각화에는 removeBatchEffect 사용.
• 교차 플랫폼은 gene 기준 병합 + study를 배치로 모델링(2-color는 별도 분석 권장).
• 메타분석은 스터디별 효과/SE를 구해 랜덤효과로 결합(metafor).
• GSVA/ssGSEA로 샘플별 경로 활성도를 정량화, WGCNA로 공발현 모듈을 탐색.

← Back: Functional Enrichment Finish: Back to Microarray Tutorial E-learning Home