ORA(과대표집 분석)와 GSEA(순위기반)를 이용해 GO/KEGG/MSigDB 경로를 해석합니다. (DEG 결과 또는 전체 유전자 순위를 입력)
처음 한 번만 설치.
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager")
BiocManager::install(c(
"clusterProfiler", "org.Hs.eg.db", "enrichplot",
"DOSE", "msigdbr", "fgsea", "pathview"
))
# 시각화
install.packages(c("ggplot2", "dplyr", "tibble", "readr"))
ORA는 DEG 목록(예: adj.P<0.05 & |logFC|≥1)과 배경(universe)이 필요합니다. GSEA는 전체 유전자 순위(예: t 통계량 또는 logFC×-log10P)가 필요합니다.
library(dplyr); library(tibble); library(AnnotationDbi); library(org.Hs.eg.db)
# 1) DEG 테이블 불러오기 (Step 03 출력)
deg <- read.csv("DEG_with_annotation.csv") # 열: PROBEID, SYMBOL, logFC, adj.P.Val 등
# 2) DEG 필터
deg_sig <- deg %>% filter(adj.P.Val < 0.05, abs(logFC) >= 1)
# 3) SYMBOL → ENTREZ (clusterProfiler는 ENTREZ 사용 권장)
map <- AnnotationDbi::select(org.Hs.eg.db,
keys = unique(deg$SYMBOL),
columns = c("ENTREZID","SYMBOL"),
keytype = "SYMBOL") %>% distinct(SYMBOL, .keep_all = TRUE)
deg_sig_entrez <- map$ENTREZID[match(deg_sig$SYMBOL, map$SYMBOL)] %>% na.omit()
universe_entrez <- map$ENTREZID[match(deg$SYMBOL, map$SYMBOL)] %>% na.omit()
length(deg_sig_entrez); length(universe_entrez)
과대표집 검정. DEG set과 universe를 명확히 지정하는 것이 중요합니다.
library(clusterProfiler); library(enrichplot); library(org.Hs.eg.db)
library(ggplot2)
# GO Biological Process
ego <- enrichGO(gene = deg_sig_entrez,
universe = universe_entrez,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
# KEGG
ekegg <- enrichKEGG(gene = deg_sig_entrez,
universe = universe_entrez,
organism = "hsa",
pvalueCutoff = 0.05)
# 결과 저장
write.csv(as.data.frame(ego), "ORA_GO_BP.csv", row.names=FALSE)
write.csv(as.data.frame(ekegg), "ORA_KEGG.csv", row.names=FALSE)
# 시각화 예시
dotplot(ego, showCategory=15) + ggtitle("GO BP ORA")
barplot(ekegg, showCategory=15) + ggtitle("KEGG ORA")
# 네트워크/개념도 예시(복수 유전자 연결)
cnetplot(ego, categorySize="pvalue", showCategory=8)
msigdbr로 H(Hallmark), C2(Canonical pathways) 등 선택 가능(인간: species = "Homo sapiens").
library(msigdbr); library(dplyr)
m_df <- msigdbr(species = "Homo sapiens", category = "H") # Hallmark
msig_h <- m_df %>% select(gs_name, entrez_gene)
# ORA용 리스트로 변환
msig_list <- split(msig_h$entrez_gene, msig_h$gs_name)
length(msig_list); head(names(msig_list))
전체 유전자 순위 벡터(named numeric) 필요. 이름=ENTREZ, 값=순위점수(예: limma t 통계).
library(fgsea); library(dplyr)
# 1) 전체 순위 만들기 (limma 결과 'tt' 사용 예시)
tt <- read.csv("DEG_limma_results.csv") # coef 1 결과 등
# SYMBOL→ENTREZ 매핑
map <- AnnotationDbi::select(org.Hs.eg.db, keys=tt$PROBEID,
columns=c("ENTREZID","SYMBOL"), keytype="PROBEID")
tt2 <- tt %>% mutate(PROBEID = rownames(tt)) %>% left_join(map, by="PROBEID")
# 순위점수: t 또는 logFC * -log10(P.Value) 등
tt2 <- tt2 %>% filter(!is.na(ENTREZID)) %>% distinct(ENTREZID, .keep_all=TRUE)
ranks <- tt2$t; names(ranks) <- tt2$ENTREZID
ranks <- sort(ranks, decreasing=TRUE)
# 2) fgsea 실행 (MSigDB Hallmark 예)
fg <- fgsea(pathways = msig_list, stats = ranks, nperm = 10000)
fg <- fg %>% arrange(padj)
write.csv(fg, "GSEA_fgsea_Hallmark.csv", row.names=FALSE)
# 3) 시각화
plotEnrichment(msig_list[["HALLMARK_INTERFERON_GAMMA_RESPONSE"]], ranks)
fgsea::plotGseaTable(msig_list[head(fg$pathway, 5)], ranks, fg, gseaParam=1)
pathview로 KEGG 경로에 발현 변화를 색으로 표시.
library(pathview)
# fold change 이름=ENTREZ
fc <- tt2$logFC; names(fc) <- tt2$ENTREZID
# 예: hsa04630 (JAK-STAT)
pathview(gene.data = fc, pathway.id = "hsa04630", species = "hsa", out.suffix="example")