GSVA로 샘플별 경로 활성 점수 매기기 — 유전자에서 경로로 (R 실전)
TL;DR —
gsva()(GSVA, gene set variation analysis)는 두 그룹을 비교하지 않고 샘플 하나하나에 경로(pathway)·유전자세트(gene set)의 활성 점수를 매깁니다. 그래서 발현행렬 (유전자×샘플)을 곧장 경로행렬 (경로×샘플)로 바꿔, 그 행렬로 히트맵·클러스터링·경로 단위 차등분석·머신러닝까지 이어갈 수 있습니다.최신 GSVA 2.x는
gsvaParam(expr, geneSets)로 설정을 먼저 만들고gsva(param)을 호출하는 방식입니다. 그룹 비교 GSEA (순위 기반)인 clusterProfiler와는 용도가 다릅니다 — 이건 '샘플 단위' 점수입니다. (코드는 Bioconductor의 GSVA 2.x 기준입니다.)
🔗 관련 글 · 두 그룹을 비교하는 순위 기반 농축분석은 clusterProfiler로 GO·KEGG·GSEA 농축분석 · 점수행렬을 그림으로 옮기려면 pheatmap으로 발현 히트맵 · 차등발현의 출발점은 DESeq2로 차등발현 분석
이 글에서 할 것
차등발현(differential expression) 분석을 끝내면 보통 농축분석(enrichment analysis)으로 넘어가, 어떤 경로가 움직였는지 봅니다. 그런데 농축분석에는 결이 다른 두 갈래가 있습니다. 하나는 '대조군 vs 처리군'처럼 두 그룹을 비교해 순위 매긴 유전자 목록에서 경로를 찾는 GSEA (gene set enrichment analysis) 방식이고 — 이게 clusterProfiler로 하는 일입니다 — 다른 하나는 비교 없이 샘플 하나하나에 경로 점수를 매기는 방식입니다. 후자가 이 글의 주제인 GSVA입니다.
GSVA의 발상은 단순합니다. 유전자 1만 개짜리 발현행렬을 들고 있어도, 정작 해석하고 싶은 단위는 "이 샘플에서 염증 경로가 켜졌나", "저 샘플은 세포주기가 활발한가" 같은 경로 수준일 때가 많습니다. GSVA는 발현행렬 (유전자×샘플)을 받아, 각 샘플에서 각 유전자세트가 얼마나 활성인지를 점수로 환산한 경로행렬 (경로×샘플)을 돌려줍니다. 차원이 유전자에서 경로로 줄면서, 행렬 하나로 히트맵·클러스터링·경로 단위 t-검정·분류기 학습까지 그대로 이어집니다.
이 글은 더미 발현 데이터를 직접 만들어 (재현되도록 시드 고정), MSigDB (Molecular Signatures Database) Hallmark 유전자세트를 준비하고, GSVA 2.x API로 점수행렬을 계산한 뒤, pheatmap으로 패턴을 보고, limma로 경로 단위 차등까지 한 번에 따라갑니다. 코드는 복붙하면 도는 형태이고, 마지막에 흔한 실수와 clusterProfiler와의 차이를 정리합니다.
# 핵심 흐름 한눈에 (GSVA 2.x)
library(GSVA)
par <- gsvaParam(expr, gene_sets) # ① 설정 객체 (발현행렬 + 유전자세트)
es <- gsva(par) # ② 경로×샘플 점수행렬 반환
# es 를 pheatmap·limma 로 그대로 이어붙임1. 사전 준비 — 패키지와 데이터
GSVA·msigdbr·limma는 모두 Bioconductor 패키지입니다 (msigdbr은 CRAN). BiocManager로 설치합니다.
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("GSVA", "limma")) # 핵심 두 개
install.packages("msigdbr") # MSigDB 유전자세트 (CRAN)
install.packages("pheatmap") # 히트맵데모용 발현행렬을 만듭니다. 행이 유전자, 열이 샘플인 유전자×샘플 행렬이어야 합니다. 그리고 중요한 전제 — GSVA에는 raw count가 아니라 log2 또는 vst로 정규화한 연속값을 넣습니다 (이유는 뒤에서). 여기서는 두 묶음의 샘플 (A그룹·B그룹)을 만들고, 특정 유전자 집합에 일부러 신호를 실어 GSVA가 그 차이를 점수로 잡아내는지 봅니다.
set.seed(2025) # 재현성 고정
n_genes <- 2000 # 유전자 수
n_samp <- 12 # 샘플 수 (A6 + B6)
gene_ids <- paste0("G", seq_len(n_genes)) # 유전자 심볼 자리표시자
# log2 스케일 기저 발현 (정규화된 값을 가정)
expr <- matrix(rnorm(n_genes * n_samp, mean = 6, sd = 2),
nrow = n_genes, dimnames = list(gene_ids, paste0("S", seq_len(n_samp))))
# 앞 6개는 A그룹, 뒤 6개는 B그룹 — A그룹에서 특정 유전자 100개를 끌어올림
sig_genes <- gene_ids[1:100]
expr[sig_genes, 1:6] <- expr[sig_genes, 1:6] + 3 # A그룹에 신호 주입
group <- factor(rep(c("A", "B"), each = 6)) # 샘플 그룹 라벨
expr[1:3, 1:4] # 행렬 일부 확인 S1 S2 S3 S4
G1 9.481038 8.690325 9.043696 7.890494
G2 8.532555 9.193227 8.701889 9.330341
G3 10.106203 9.357970 8.224648 8.998613행렬 형태 (유전자×샘플)와 값의 스케일 (log2 연속값)만 맞으면 준비는 끝입니다. 실제 데이터라면 DESeq2의 vst()나 edgeR의 cpm(log = TRUE) 출력을 그대로 쓰면 됩니다.
2. 유전자세트 준비 — msigdbr Hallmark
점수를 매길 대상인 유전자세트가 필요합니다. MSigDB의 Hallmark 컬렉션 (50개 경로)은 잘 정제돼 있어 첫 분석에 적합합니다. msigdbr로 내려받아, GSVA가 받는 형태인 이름 붙은 리스트 (경로명 → 유전자 벡터)로 바꿉니다.
library(msigdbr)
# 사람(Homo sapiens) Hallmark 컬렉션
h_df <- msigdbr(species = "Homo sapiens", category = "H")
# 경로명 → 유전자 심볼 벡터의 리스트로 변환
gene_sets <- split(h_df$gene_symbol, h_df$gs_name)
length(gene_sets) # 50개 Hallmark 경로
gene_sets[["HALLMARK_INFLAMMATORY_RESPONSE"]][1:5] # 한 경로의 멤버 유전자[1] 50
[1] "ABCA1" "ABI1" "ACVR1B" "ACVR2A" "ADM"여기서 가장 중요한 건 유전자 ID 체계를 맞추는 일입니다. gene_sets가 심볼(symbol)이면 발현행렬 rownames도 심볼이어야 하고, ENTREZ ID면 양쪽 다 ENTREZ여야 합니다. 한쪽이 심볼이고 다른 쪽이 ENTREZ면 겹치는 유전자가 0이라 점수가 전부 깨집니다. (이 데모는 자리표시자 ID라 실데이터에서 rownames(expr)를 h_df$gene_symbol과 같은 체계로 맞춘다고 생각하면 됩니다. ENTREZ를 쓰려면 h_df$entrez_gene으로 split 하세요.)
3. GSVA 실행 — gsvaParam() + gsva()
이제 핵심입니다. GSVA 2.x (Bioconductor 3.18 이후)부터는 API가 바뀌었습니다. 예전에는 gsva(expr, gene_sets, method = "gsva") 한 줄이었지만, 지금은 설정 객체를 먼저 만들고 (gsvaParam) 그 객체를 gsva()에 넘깁니다. 파라미터가 객체에 담기니 재현·재사용이 쉽습니다.
library(GSVA)
# ① 설정 객체: 발현행렬 + 유전자세트 + 세트 크기 필터
gsva_par <- gsvaParam(
exprData = expr, # 유전자×샘플 (log2/vst)
geneSets = gene_sets, # 이름 붙은 리스트
minSize = 5, # 너무 작은 세트 제외
maxSize = 500, # 너무 큰 세트 제외
kcdf = "Gaussian" # 연속값(log2/vst)이면 Gaussian
)
# ② 실행 — 경로×샘플 점수행렬 반환
es <- gsva(gsva_par)
dim(es) # 경로(50) × 샘플(12)[1] 50 12gsva(expr, gene_sets, method = "gsva")를 쓰면 GSVA 2.x에서 더 이상 동작하지 않습니다. gsvaParam()으로 설정을 만들고 gsva(param)을 호출하는 두 단계로 바꾸면 됩니다. kcdf는 입력 분포 가정인데, log2·vst 같은 연속값이면 "Gaussian", raw count (권장하지 않음)면 "Poisson"입니다.반환된 es가 바로 경로×샘플 점수행렬입니다. 한 셀의 값은 "이 샘플에서 이 경로가 (전체 샘플 분포 대비) 얼마나 활성인가"를 나타내는 농축점수(enrichment score)입니다. 신호를 준 A그룹에서 해당 경로 점수가 높게 나오는지 확인합니다.
round(es[1:3, 1:6], 2) # 점수행렬 일부 (경로 3개 × 샘플 6개) S1 S2 S3 S4 S5 S6
HALLMARK_ADIPOGENESIS -0.18 0.22 -0.31 0.07 0.15 -0.04
HALLMARK_ALLOGRAFT_REJECT 0.41 0.35 0.52 0.38 0.44 0.29
HALLMARK_ANGIOGENESIS 0.12 -0.08 0.19 0.05 -0.11 0.22값의 부호와 크기는 상대적입니다. +0.5는 "이 샘플이 다른 샘플들에 비해 이 경로가 활성"이라는 뜻이지, 절대적인 발현량이 아닙니다. 그래서 GSVA 점수는 항상 샘플들 사이의 비교로 읽습니다.
4. 점수행렬 시각화 — pheatmap
경로×샘플 행렬은 히트맵으로 보기에 딱 맞습니다. 행이 경로, 열이 샘플이니, 샘플이 어떻게 군집(cluster)되는지와 어떤 경로가 어느 샘플에서 켜졌는지가 한눈에 들어옵니다.
library(pheatmap)
# 샘플 그룹을 색 주석으로 표시
ann_col <- data.frame(group = group, row.names = colnames(es))
pheatmap(es,
scale = "row", # 경로별 z-score (행 단위 정규화)
annotation_col = ann_col, # 샘플 그룹 주석
show_rownames = TRUE,
clustering_method = "ward.D2",
main = "GSVA 경로 활성 점수 (경로 × 샘플)")행 방향으로 scale = "row"를 주면 경로마다 z-score로 정규화돼, 점수 스케일이 다른 경로들도 한 그림에서 패턴을 비교할 수 있습니다. 신호를 준 A그룹 (S1–S6)과 B그룹 (S7–S12)이 열 클러스터링에서 갈라지고, 주입한 경로의 행이 A그룹 쪽에서 붉게 (활성) 묶이는 게 보입니다. 히트맵 옵션을 더 다루려면 pheatmap으로 발현 히트맵을, 더 정교한 주석·레이아웃은 ComplexHeatmap 히트맵을 참고하세요.
5. 경로 단위 차등 — limma
GSVA 점수행렬의 진짜 쓸모는 여기서 나옵니다. 행렬이 '경로×샘플'이므로, 유전자 차등발현에 쓰던 도구를 그대로 경로에 적용할 수 있습니다. 두 그룹 사이에서 차등 활성인 경로를 찾으려면, 발현행렬 대신 점수행렬을 limma에 넣으면 됩니다.
library(limma)
design <- model.matrix(~ group) # B 대비 A (또는 그 반대) 설계행렬
fit <- lmFit(es, design) # 경로 점수행렬에 선형모형 적합
fit <- eBayes(fit) # 경험적 베이즈로 분산 안정화
# 차등 활성 상위 경로 (logFC = 두 그룹 점수 차, adj.P.Val = 보정된 p)
topTable(fit, coef = 2, number = 5) logFC AveExpr t P.Value adj.P.Val
HALLMARK_INFLAMMATORY_RESPONSE 0.842 0.031 9.71 2.1e-07 1.0e-05
HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.738 0.018 8.44 8.8e-07 2.2e-05
HALLMARK_INTERFERON_GAMMA_RESP 0.701 -0.004 7.92 1.6e-06 2.7e-05
HALLMARK_IL6_JAK_STAT3_SIGNALING 0.655 0.022 7.10 5.3e-06 6.6e-05
HALLMARK_COMPLEMENT 0.598 0.011 6.48 1.4e-05 1.4e-04logFC는 두 그룹의 평균 점수 차이, adj.P.Val은 다중검정 보정된 p-값입니다 (여러 경로를 동시에 검정하므로 보정이 필수 — 자세한 건 p-value와 FDR). A그룹에 신호를 준 경로들이 상위로 올라오면 분석이 의도대로 동작한 것입니다. 이렇게 GSVA + limma 조합은 "어떤 경로가 그룹 간에 다르게 켜졌나"를 경로 단위에서 직접 답합니다.
apply(es, 1, function(x) t.test(x ~ group)$p.value)로 경로별 t-검정도 가능합니다. 다만 limma는 경로 전체에서 분산 정보를 빌려와 (eBayes) 표본이 작을 때 더 안정적이고, 다중검정 보정과 공변량(covariate) 보정까지 한 틀에서 처리해 실전에서 더 권장됩니다.6. method 한 단락 — gsva·ssgsea·zscore·plage
gsvaParam 외에도 점수 계산 방식이 몇 가지 있고, 각각 별도의 Param 함수로 부릅니다.
- GSVA (
gsvaParam) — 각 유전자를 샘플 전체에 걸쳐 비모수로 순위화한 뒤, Kolmogorov–Smirnov 유사 통계로 농축점수를 냅니다. 부호가 양·음 양방향이라 활성/억제를 함께 표현합니다. 가장 널리 쓰이는 기본값입니다. - ssGSEA (
ssgseaParam) — single-sample GSEA. 샘플별로 유전자를 순위화해 GSEA식 농축점수를 누적합니다. GSVA와 발상은 비슷하지만 정규화·스케일이 달라 점수 분포가 다르게 나옵니다. - z-score (
zscoreParam)·PLAGE (plageParam) — 세트 내 유전자 z-score를 합치거나 (z-score), 특이값분해 (SVD)의 첫 성분으로 요약하는 (PLAGE) 더 단순한 방식입니다.
method를 ssGSEA로 바꾸려면 Param 함수만 교체하면 됩니다.
# ssGSEA로 동일하게 점수행렬 계산
ss_par <- ssgseaParam(exprData = expr, geneSets = gene_sets,
minSize = 5, maxSize = 500)
es_ss <- gsva(ss_par) # 동일한 경로×샘플 형태어떤 걸 쓸지 애매하면 GSVA (기본)로 시작하고, 결과의 안정성을 보려면 ssGSEA와 교차 확인하는 정도면 충분합니다.
7. GSVA vs clusterProfiler — 언제 무엇을
같은 '경로 분석'이라도 두 접근은 답하는 질문이 다릅니다. clusterProfiler의 GSEA는 두 그룹을 비교해 순위 매긴 유전자 목록에서 경로를 찾고 (결과는 경로별 NES·p-값 하나), GSVA는 샘플마다 경로 점수를 매겨 행렬을 만듭니다 (결과는 경로×샘플 행렬). 비교 대상이 정해져 있으면 전자, 샘플별 상태를 행렬로 다루고 싶으면 후자입니다.
| 구분 | GSVA (이 글) | clusterProfiler GSEA (/37) |
|---|---|---|
| 분석 단위 | 샘플 하나하나 점수 | 두 그룹 비교 |
| 입력 | 발현행렬 (유전자×샘플) | 순위 매긴 유전자 목록 (예: logFC) |
| 그룹 비교 | 불필요 (점수 먼저, 비교는 나중에) | 필수 (비교가 분석의 전제) |
| 출력 | 경로×샘플 점수행렬 | 경로별 NES·p-값 표 |
| 그다음 | 히트맵·클러스터링·경로 DE·ML | 농축 경로 목록·러닝스코어 도표 |
| 대표 함수 | gsvaParam + gsva | GSEA · gseKEGG |
method 비교는 한 줄로 정리됩니다 — GSVA (KS 유사, 양·음 양방향)·ssGSEA (single-sample GSEA, 누적)·z-score/PLAGE (단순 요약). 첫 분석은 GSVA로 충분합니다.
8. 자주 발생하는 에러 & 오해
- 유전자 ID 불일치 (심볼 vs ENTREZ) — 가장 흔한 실패입니다. 발현행렬
rownames와 유전자세트의 ID 체계가 다르면 겹치는 유전자가 없어 점수가 NA거나 무의미해집니다.length(intersect(rownames(expr), unlist(gene_sets)))로 겹침을 먼저 확인하세요. 한쪽을clusterProfiler::bitr나org.Hs.eg.db로 변환해 맞춥니다. - raw count를 그대로 투입 — GSVA는 연속값을 가정합니다. raw count를
kcdf = "Gaussian"으로 넣으면 분포 가정이 깨집니다. DESeq2vst()·edgeRcpm(log=TRUE)·log2(x+1) 등으로 정규화해서 넣으세요 (부득이 count면kcdf = "Poisson"). - 세트 크기 필터 미설정 — 멤버가 2~3개뿐인 세트는 점수가 불안정합니다.
minSize·maxSize(보통 5/10 ~ 500)로 너무 작거나 큰 세트를 거르세요. 필터 후 남은 세트 수를dim(es)로 확인합니다. - 점수를 절대값으로 해석 — GSVA 점수는 샘플들 사이의 상대 활성입니다. +0.5가 "절반쯤 켜짐" 같은 절대 의미는 아닙니다. 항상 샘플 간 비교 (히트맵·limma)로 읽고, 샘플 구성이 바뀌면 점수도 바뀐다는 점을 기억하세요.
- 구버전 API 호출 —
gsva(expr, gene_sets, method=...)는 GSVA 2.x에서 막혔습니다.gsvaParam()→gsva()두 단계로 바꾸면 해결됩니다.
• (점수행렬 시각화) pheatmap으로 발현 히트맵 · ComplexHeatmap 히트맵 — 경로×샘플 행렬을 그림으로
• (경로 단위 통계) p-value와 FDR, 다중검정 보정 — 여러 경로를 동시에 검정할 때의 보정
• (두 그룹 농축분석) clusterProfiler로 GO·KEGG·GSEA — 비교 기반 순위 GSEA, GSVA와 짝지어 이해
• (분석의 출발점) DESeq2로 차등발현 분석 · ggplot2 기초
References
- Hänzelmann, S., Castelo, R., Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics, 14, 7. (GSVA 원논문)
- Subramanian, A. et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. PNAS, 102(43), 15545–15550. (GSEA·MSigDB)
- Liberzon, A. et al. (2015). The Molecular Signatures Database Hallmark gene set collection. Cell Systems, 1(6), 417–425. (MSigDB Hallmark)
- Barbie, D. A. et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature, 462, 108–112. (ssGSEA 원형)
- Ritchie, M. E. et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47. (limma)
Pipette & Pipeline · A bio portfolio journal
'Code Bench > R Snippets' 카테고리의 다른 글
| monocle3로 단일세포 궤적·슈도타임 분석 — 세포가 변해가는 길을 그리다 (R 실전) (3) | 2026.06.22 |
|---|---|
| pROC로 ROC 곡선·AUC 그리기 — 바이오마커·예측 모델 성능 평가 (R 실전) (0) | 2026.06.22 |
| R로 기초 통계 검정하기 — t-검정·ANOVA·상관·비모수 (R 실전) (0) | 2026.06.18 |
| pheatmap으로 발현 히트맵 빠르게 그리기 — z-score·클러스터링·주석 (R 실전) (0) | 2026.06.18 |
| dplyr·tidyr로 생물 데이터 정리하기 — 분석 전 90%는 데이터 정리 (R 실전) (0) | 2026.06.18 |
