TL;DR — Seurat로 군집화까지 끝내면 "어떤 세포가 몇 종류 있나"는 답했지만, 분화나 활성화처럼 세포가 어떤 순서로 변해가는지는 아직 모릅니다. 이걸 재구성하는 도구가 R 패키지 monocle3입니다. 세포 상태의 변화를 잇는 궤적(trajectory)을 학습하고, 시작점 기준으로 진행 정도를 매기는 슈도타임(pseudotime)을 계산합니다.

흐름은 cell_data_set (CDS) 만들기 → preprocess_cds (PCA) → reduce_dimension (UMAP) → cluster_cells → learn_graph (주그래프 학습) → order_cells (루트 지정 → 슈도타임) → plot_cells로 색칠하기입니다. 마지막에 슈도타임을 따라 변하는 유전자를 graph_test로 찾습니다.

군집(cluster)이 세포를 이산적인 그룹으로 나눈다면, 궤적은 그 사이의 연속적인 변화를 그립니다. 분화 단계처럼 경계가 흐릿하게 이어지는 과정엔 궤적이 더 맞습니다.

🔗 관련 글  ·  그 전 단계인 군집화는 Seurat로 단일세포 RNA-seq 분석  ·  궤적 유전자를 히트맵으로 보려면 ComplexHeatmap 히트맵  ·  출판용 그림 문법은 ggplot2 기초

이 글에서 할 것

Seurat 글에서 단일세포 RNA-seq (scRNA-seq) 데이터를 품질관리(QC)하고 군집으로 나눈 뒤 셀타입까지 붙였습니다. 그 다음 자연스럽게 떠오르는 질문이 있습니다. "이 세포들이 정지 상태로 흩어져 있는 게 아니라, 어떤 순서로 변해가는 중이라면 그 경로는 어떻게 생겼을까?" 줄기세포에서 성숙 세포로 가는 분화, 면역세포의 활성화처럼 연속적인 상태 변화를 다룰 때 던지는 질문입니다.

여기에 답하는 두 개념이 궤적(trajectory)과 슈도타임(pseudotime)입니다. 궤적은 세포 상태의 변화를 연결한 그래프입니다. 차원축소 공간에서 비슷한 세포끼리 이웃하도록 펼친 뒤, 그 위에 변화의 흐름을 따라가는 곡선(주그래프, principal graph)을 얹습니다. 슈도타임은 그 곡선 위에서 시작점으로부터 각 세포가 얼마나 진행했는지를 나타내는 값입니다. 이름에 '타임'이 들어가지만 실제 시간(분·시간·일)이 아니라 진행 정도의 순서라는 점이 핵심입니다. 세포를 한 시점에 한꺼번에 잡아도, 그 안에 분화의 여러 단계가 섞여 있으면 그 단계를 한 줄로 세울 수 있다는 발상입니다.

이 글은 monocle3 (Cao et al. 2019)로 작은 더미 데이터를 만들어 CDS를 구성하고, 전처리·차원축소·군집화를 거쳐 주그래프를 학습하고, 루트(시작 세포)를 지정해 슈도타임을 매기고, 슈도타임으로 색칠한 궤적을 그립니다. 마지막에 슈도타임을 따라 발현이 변하는 유전자를 graph_test로 찾고, 군집과 궤적의 차이, 흔한 실수를 정리합니다.

# 핵심 흐름 한눈에 (monocle3)
library(monocle3)
cds <- new_cell_data_set(expr, cell_metadata, gene_metadata)  # ① CDS 생성
cds <- preprocess_cds(cds, num_dim = 50)   # ② PCA
cds <- reduce_dimension(cds)               # ③ UMAP
cds <- cluster_cells(cds)                  # ④ 군집·파티션
cds <- learn_graph(cds)                    # ⑤ 주그래프 학습
cds <- order_cells(cds, root_pr_nodes = root)  # ⑥ 루트 → 슈도타임
plot_cells(cds, color_cells_by = "pseudotime") # ⑦ 색칠

1. 사전 준비 — 설치와 입력 데이터

monocle3은 CRAN이 아니라 GitHub에서 설치합니다. 의존성으로 Bioconductor 패키지 몇 개(BiocGenerics·SingleCellExperiment 등)가 먼저 필요합니다.

# Bioconductor 의존성 먼저
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("BiocGenerics", "DelayedArray", "DelayedMatrixStats",
                       "limma", "S4Vectors", "SingleCellExperiment",
                       "SummarizedExperiment", "batchelor"))

# monocle3 본체 (GitHub)
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle3")

monocle3의 기본 자료구조는 cell_data_set (CDS)입니다. 만드는 데 필요한 건 세 가지입니다. 발현 행렬(유전자 × 세포, 보통 희소행렬), 세포 정보를 담은 cell_metadata (행=세포), 유전자 정보를 담은 gene_metadata (행=유전자)입니다. 여기서는 설명용으로 분화의 세 단계(초기·중간·말기)를 흉내 낸 작은 더미 데이터를 만듭니다.

library(monocle3)
set.seed(42)                                  # 재현성 고정

n_genes <- 200
n_cells <- 600                                # 단계당 200세포 × 3단계

# 잠재 진행도(latent progression) 0→1 을 따라 세포를 배치
progress <- sort(runif(n_cells))              # 0~1 사이 진행도

# 발현 행렬: 일부 유전자는 진행도에 따라 증가/감소하도록 설계
expr <- matrix(rpois(n_genes * n_cells, lambda = 1),
               nrow = n_genes, ncol = n_cells)
# 앞쪽 30개 유전자는 진행도가 커질수록 발현↑ (분화 마커 흉내)
for (g in 1:30)
  expr[g, ] <- rpois(n_cells, lambda = 1 + 8 * progress)
# 다음 30개 유전자는 진행도가 커질수록 발현↓ (초기 마커 흉내)
for (g in 31:60)
  expr[g, ] <- rpois(n_cells, lambda = 1 + 8 * (1 - progress))

rownames(expr) <- paste0("gene", seq_len(n_genes))
colnames(expr) <- paste0("cell", seq_len(n_cells))

# 세포 메타데이터: 진행도 구간으로 3단계 라벨
cell_meta <- data.frame(
  row.names = colnames(expr),
  stage = cut(progress, breaks = c(-Inf, 0.33, 0.66, Inf),
              labels = c("초기", "중간", "말기"))
)

# 유전자 메타데이터: gene_short_name 열은 필수 (plot_cells가 참조)
gene_meta <- data.frame(
  row.names = rownames(expr),
  gene_short_name = rownames(expr)
)

cds <- new_cell_data_set(expression_data = expr,
                         cell_metadata   = cell_meta,
                         gene_metadata   = gene_meta)
cds

class: cell_data_set
dim: 200 600
assays(1): counts
rownames(200): gene1 gene2 ... gene199 gene200
colnames(600): cell1 cell2 ... cell599 cell600
colData names(2): stage Size_Factor

gene_metadata에 gene_short_name 열을 넣은 건 빠뜨리기 쉬운 필수 조건입니다. plot_cells가 유전자를 이름으로 찾을 때 이 열을 참조하므로, 없으면 나중에 에러가 납니다. 실제 분석이라면 발현 행렬은 Read10X로 읽은 카운트, cell_metadata는 샘플·조건·군집 정보가 됩니다.

2. 전처리와 차원축소 — preprocess_cds → reduce_dimension

CDS가 준비되면 주성분분석(PCA)으로 차원을 줄이고(preprocess_cds), 그 결과를 UMAP으로 2차원에 펼칩니다(reduce_dimension). Seurat의 RunPCA·RunUMAP에 해당하는 단계입니다.

# ② 전처리: 정규화 + PCA (num_dim = 사용할 주성분 수)
cds <- preprocess_cds(cds, num_dim = 50)

# (선택) 배치효과가 있으면 여기서 보정
# cds <- align_cds(cds, alignment_group = "batch")

# ③ 차원축소: UMAP (기본 reduction_method = "UMAP")
cds <- reduce_dimension(cds)

# 단계 라벨로 색칠해 UMAP 확인
plot_cells(cds, color_cells_by = "stage",
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)   # 아직 궤적 학습 전

num_dim은 PCA에서 몇 개의 주성분을 쓸지 정합니다. 너무 작으면 변화의 축을 놓치고, 너무 크면 잡음까지 끌어와 군집·궤적이 흔들립니다. 실제 데이터에서는 plot_pc_variance_explained(cds)로 분산 설명 곡선이 꺾이는 지점을 보고 정합니다. 배치(batch)가 섞여 있다면 align_cds로 보정한 뒤 차원축소를 해야 합니다. 보정을 건너뛰면 생물학적 변화가 아니라 배치 차이가 궤적으로 잡히는 함정에 빠집니다.

3. 군집화 — cluster_cells

cluster_cells는 UMAP 공간에서 세포를 군집(cluster)으로 나누고, 동시에 파티션(partition)이라는 더 큰 단위로도 묶습니다. 파티션은 서로 떨어진 세포 덩어리로, monocle3는 파티션마다 별도의 궤적을 학습합니다. 즉 연결되지 않은 세포 집단을 억지로 한 궤적에 잇지 않습니다.

# ④ 군집화 (resolution 으로 군집 해상도 조절)
cds <- cluster_cells(cds, resolution = 1e-2)

# 군집으로 색칠
plot_cells(cds, color_cells_by = "cluster",
           show_trajectory_graph = FALSE)

# 파티션 확인 (궤적이 그려지는 단위)
plot_cells(cds, color_cells_by = "partition",
           show_trajectory_graph = FALSE)

resolution을 키우면 군집이 잘게 쪼개지고, 줄이면 뭉칩니다. 군집 해상도가 너무 높으면 하나의 연속 분화가 여러 조각으로 끊겨 궤적이 부자연스러워지고, 너무 낮으면 서로 다른 경로가 한 덩어리로 뭉쳐 분기(branch)를 놓칩니다. 데이터에 맞춰 몇 번 조정하는 게 보통입니다.

4. 주그래프 학습 — learn_graph

이제 궤적의 뼈대를 학습합니다. learn_graph는 세포가 펼쳐진 공간 위로 변화의 흐름을 따라가는 주그래프(principal graph)를 적합합니다. 이 곡선이 "어디서 시작해 어디로 갈라지고 어디로 끝나는지"의 골격입니다.

# ⑤ 주그래프 학습
cds <- learn_graph(cds)

# 궤적(주그래프)을 단계 라벨 위에 겹쳐 그리기
plot_cells(cds,
           color_cells_by = "stage",
           label_groups_by_cluster = FALSE,
           label_leaves = TRUE,            # 가지 끝(잎) 표시
           label_branch_points = TRUE)     # 분기점 표시

그림에 세포 위로 검은 선이 얹히는데, 이게 주그래프입니다. 끝점(잎, leaf)은 분화의 종착 상태, 분기점(branch point)은 한 경로가 둘 이상으로 갈리는 지점입니다. 아직 이 그래프는 방향이 없습니다. 어느 끝이 시작이고 어느 끝이 마지막인지는 다음 단계에서 사람이 지정해야 합니다.

5. 루트 지정과 슈도타임 — order_cells

궤적에 방향을 주는 단계입니다. order_cells에 루트(root, 시작 지점)를 알려 주면, 그 지점에서부터 그래프를 따라간 거리로 각 세포의 슈도타임이 정해집니다. 루트를 지정하지 않으면 슈도타임의 기준점이 없어 값 자체가 의미를 잃습니다 — 가장 중요한 한 단계입니다.

# 루트 주그래프 노드를 프로그램으로 고르는 헬퍼
# (여기서는 'stage == 초기' 세포가 가장 많은 노드를 시작점으로)
get_root_node <- function(cds, start_stage = "초기") {
  cell_ids <- which(colData(cds)$stage == start_stage)
  closest <- igraph::V(principal_graph(cds)$UMAP)$name[
    as.numeric(names(which.max(table(
      principal_graph_aux(cds)$UMAP$pr_graph_cell_proj_closest_vertex[cell_ids, ]
    ))))
  ]
  closest
}

# ⑥ 루트 지정 → 슈도타임 계산
cds <- order_cells(cds, root_pr_nodes = get_root_node(cds))

# ⑦ 슈도타임으로 색칠 (연속 색 = 진행 정도)
plot_cells(cds,
           color_cells_by = "pseudotime",
           label_leaves = FALSE,
           label_branch_points = FALSE)

RStudio에서 인자 없이 order_cells(cds)만 실행하면 마우스로 루트를 클릭하는 창이 뜹니다. 탐색할 땐 편하지만, 재현 가능한 분석에는 위처럼 root_pr_nodes로 시작 노드를 코드로 못박는 편이 낫습니다. 슈도타임으로 색칠한 그림에서 색이 한쪽 끝(루트)에서 반대쪽으로 매끄럽게 번지면, 세포들이 그 방향으로 진행하는 연속 과정으로 잘 정렬된 것입니다. 각 세포의 슈도타임 값은 pseudotime(cds)로 꺼내 다른 분석에 쓸 수 있습니다.

6. 군집(이산) vs 궤적(연속), 그리고 슈도타임 유전자

같은 데이터라도 군집과 궤적은 보는 각도가 다릅니다. 둘은 경쟁이 아니라 보완 관계입니다.

궤적을 손에 쥐었으면 다음 질문은 "슈도타임을 따라 발현이 변하는 유전자는 무엇인가"입니다. graph_test가 주그래프 상에서 공간적으로 자기상관(autocorrelation)이 큰 유전자, 즉 궤적의 특정 구간에 발현이 몰리는 유전자를 찾아 줍니다. 통계량은 모런 지수(Moran's I)입니다.

# 주그래프를 따라 발현이 변하는 유전자 검정 (Moran's I)
deg_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)

# 유의한 유전자만: morans_I 가 클수록 궤적을 따라 강하게 변함
sig_genes <- subset(deg_res, q_value < 0.05)
sig_genes <- sig_genes[order(-sig_genes$morans_I), ]
head(sig_genes[, c("gene_short_name", "morans_I", "q_value")])

# 상위 유전자를 슈도타임 궤적 위에 발현으로 표시
plot_cells(cds,
           genes = head(rownames(sig_genes), 4),
           show_trajectory_graph = FALSE,
           label_cell_groups = FALSE,
           label_leaves = FALSE)

        gene_short_name   morans_I       q_value
gene7             gene7  0.6831042  3.10e-58
gene42           gene42  0.6552180  1.22e-54
gene15           gene15  0.6310977  8.04e-51
gene38           gene38  0.5994461  2.71e-47

여기서 모런 지수(Moran's I)는 −1에서 +1 사이이고, 값이 클수록 그 유전자가 궤적의 특정 구간에 발현이 집중된다는 뜻입니다. 더미 데이터에서 진행도에 따라 증가·감소하도록 심어 둔 앞쪽 유전자들이 상위에 올라옵니다. 실제 분석이라면 이 유전자 목록을 슈도타임 순으로 정렬해 ComplexHeatmap으로 그리면, 분화 단계마다 켜지고 꺼지는 유전자 프로그램이 한눈에 보입니다.

7. 자주 발생하는 에러 & 오해

• 루트를 지정하지 않음 → 슈도타임이 무의미 — order_cells에 root_pr_nodes(또는 root_cells)를 주지 않으면 기준점이 없어 슈도타임 값이 뒤죽박죽이 됩니다. 시작 세포는 알고리즘이 아니라 생물학(줄기·전구세포 마커, 정지 상태)으로 정해 코드로 못박으세요.
• monocle2 API와 혼동 — monocle3는 이전 버전(monocle2)과 함수가 다릅니다. newCellDataSet·differentialGeneTest·reduceDimension(method = "DDRTree")는 monocle2 문법이라 monocle3에서는 동작하지 않습니다. monocle3에서는 new_cell_data_set·graph_test·reduce_dimension을 씁니다.
• gene_short_name 열 누락 — gene_metadata에 gene_short_name 열이 없으면 plot_cells(genes = ...)가 유전자를 못 찾아 에러를 냅니다. CDS를 직접 만들 땐 이 열을 꼭 넣으세요.
• 배치효과를 보정하지 않음 — 여러 샘플·실험 배치가 섞였는데 align_cds로 보정하지 않으면, 생물학적 변화가 아니라 배치 차이가 궤적으로 잡힙니다. reduce_dimension 전에 align_cds(cds, alignment_group = "batch")로 정렬하세요.
• 군집 해상도를 방치 — cluster_cells의 resolution이 너무 높으면 연속 분화가 토막 나고, 너무 낮으면 분기를 놓칩니다. 궤적이 어색하면 해상도부터 조정해 보세요.
• 슈도타임을 실제 시간으로 오해 — 슈도타임은 시작점 기준의 상대적 진행 순서일 뿐, 분·시간 같은 물리적 시간이 아닙니다. 두 세포의 슈도타임 차이가 실제 경과 시간에 비례하지 않습니다.

References

1. Cao, J., Spielmann, M., Qiu, X. et al. (2019). The single-cell transcriptional landscape of mammalian organogenesis. Nature, 566, 496–502. (monocle3 원논문)
2. Trapnell, C., Cacchiarelli, D., Grimsby, J. et al. (2014). The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature Biotechnology, 32, 381–386. (슈도타임 개념)
3. Qiu, X., Mao, Q., Tang, Y. et al. (2017). Reversed graph embedding resolves complex single-cell trajectories. Nature Methods, 14, 979–982. (주그래프 학습)
4. Street, K. et al. (2018). Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics, 19, 477. (대안 궤적 도구)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 어떤 유전자나 점수가 질병군과 대조군을 얼마나 잘 가르는지 보려면 ROC (receiver operating characteristic) 곡선과 AUC (area under the curve)를 씁니다. R의 pROC 패키지로 roc(response, predictor) 한 줄이면 곡선·AUC·신뢰구간·최적 컷오프가 다 나오고, roc.test()로 두 바이오마커의 AUC까지 DeLong (DeLong et al. 1988) 검정으로 비교합니다.

ROC는 모든 컷오프에서 민감도(sensitivity)와 1−특이도(specificity)가 그리는 궤적이고, AUC는 0.5(무작위)에서 1.0(완벽) 사이입니다. 이 글은 더미 데이터로 점수 → ROC → AUC → 컷오프 → 두 모델 비교까지 복붙하면 도는 형태로 따라갑니다.

🔗 관련 글  ·  검정 자체가 처음이면 R로 기초 통계 검정하기  ·  p-값을 여러 번 볼 때의 보정은 p-value와 FDR, 다중검정 보정  ·  점수의 출발점인 차등발현은 DESeq2로 차등발현 분석

이 글에서 할 것

차등발현(differential expression)이나 어떤 점수 모델을 만들고 나면 곧장 따라오는 질문이 있습니다. "이 유전자(혹은 위험점수)가 환자군과 정상군을 정말 잘 구분하나?" 단순히 두 그룹의 평균이 다른 것(통계 검정의 영역)과, 그 차이로 개별 샘플을 분류할 수 있느냐는 다른 문제입니다. 후자를 평가하는 표준 도구가 ROC 곡선과 AUC입니다.

발상은 이렇습니다. 바이오마커 점수에 컷오프(cutoff)를 하나 정하면, 그보다 높은 샘플은 '양성'으로 부릅니다. 컷오프를 낮추면 진짜 환자를 더 많이 잡지만(민감도↑) 정상인도 양성으로 잘못 부르고(특이도↓), 높이면 그 반대입니다. 이 맞교환(trade-off)을 모든 컷오프에 대해 그린 게 ROC 곡선입니다. 가로축은 1−특이도(거짓양성률, FPR), 세로축은 민감도(참양성률, TPR)입니다. 곡선이 왼쪽 위 구석에 붙을수록 좋은 분류기이고, 그 아래 면적이 AUC입니다.

이 글은 pROC 패키지 (Robin et al. 2011)로 더미 데이터를 만들어 ROC를 그리고, AUC와 신뢰구간(confidence interval)을 구하고, Youden 지수로 최적 컷오프를 찾고, 두 바이오마커의 AUC를 DeLong 검정으로 비교하고, ggplot2로 출판용 곡선까지 그립니다. 마지막에 방향(direction) 뒤집힘 같은 흔한 실수와 ROC vs PR 곡선의 차이를 정리합니다.

# 핵심 흐름 한눈에 (pROC)
library(pROC)
roc_obj <- roc(response, predictor)   # ① 정답 라벨 + 점수 → ROC 객체
plot(roc_obj); auc(roc_obj)           # ② 곡선 + AUC
coords(roc_obj, "best",               # ③ Youden 최적 컷오프
       best.method = "youden")
roc.test(roc_a, roc_b)                # ④ 두 AUC 비교 (DeLong)

1. 사전 준비 — 패키지와 데이터

pROC와 ggplot2는 모두 CRAN 패키지입니다. install.packages()로 설치합니다.

install.packages("pROC")      # ROC·AUC·DeLong 검정
install.packages("ggplot2")   # ggroc로 출판용 곡선

데모용 데이터를 만듭니다. 핵심은 질병군/대조군 라벨과 그에 대한 바이오마커 점수 두 개입니다. 현실성을 위해 두 그룹의 점수 분포를 일부러 겹치게 둡니다 (완벽히 갈리면 AUC가 1.0이라 재미가 없습니다). 환자군 점수를 대조군보다 평균만 조금 높게 잡습니다.

set.seed(2025)                                   # 재현성 고정
n <- 200                                         # 샘플 수 (환자 100 + 정상 100)

# 정답 라벨: 0 = 정상(control), 1 = 질병(case)
status <- factor(rep(c("control", "case"), each = 100),
                 levels = c("control", "case"))  # 레벨 순서를 명시 (뒤에서 중요)

# 바이오마커 A 점수: 환자군 평균이 살짝 높지만 분포가 겹침
biomarker_A <- c(rnorm(100, mean = 5.0, sd = 1.5),   # control
                 rnorm(100, mean = 6.2, sd = 1.5))   # case

# 바이오마커 B 점수: A보다 두 그룹을 더 잘 가름 (평균 차이 큼)
biomarker_B <- c(rnorm(100, mean = 5.0, sd = 1.3),   # control
                 rnorm(100, mean = 7.0, sd = 1.3))   # case

dat <- data.frame(status, biomarker_A, biomarker_B)
head(dat, 3)                                     # 데이터 확인

   status biomarker_A biomarker_B
1 control    3.879049    4.574687
2 control    5.539645    5.045524
3 control    5.117417    4.929475

라벨이 두 수준(정상·질병)이고 점수가 연속값이면 준비는 끝입니다. 실제 분석이라면 status는 진단 결과, biomarker_A는 유전자 발현·위험점수·검사 수치 등이 됩니다.

2. ROC 객체 만들기 — roc()

pROC의 출발점은 roc() 하나입니다. 정답 라벨(response)과 점수(predictor)를 넣으면 ROC 곡선 계산에 필요한 모든 것이 든 객체가 나옵니다.

library(pROC)

# response = 정답 라벨, predictor = 바이오마커 점수
roc_A <- roc(response = dat$status,
             predictor = dat$biomarker_A,
             levels = c("control", "case"),   # 음성, 양성 순서 명시
             direction = "<")                  # case의 점수가 더 '높다'고 지정

roc_A                                          # 객체 요약 출력

Call:
roc.default(response = dat$status, predictor = dat$biomarker_A, levels = c("control", "case"), direction = "<")

Data: dat$biomarker_A in 100 controls (dat$status control) < 100 cases (dat$status case).
Area under the curve: 0.7138

levels와 direction을 명시한 건 의도적입니다. levels = c("control", "case")는 "앞이 음성, 뒤가 양성"을 못박고, direction = "<"는 "양성(case)의 점수가 음성보다 높다"는 뜻입니다. 이 둘을 지정하지 않으면 pROC가 데이터에서 자동으로 추정하는데, 점수가 거꾸로 든 경우(낮을수록 질병) AUC가 0.5 아래로 뒤집혀 나옵니다. 가장 흔한 실수라 뒤에서 따로 다룹니다.

3. 곡선 그리기 + AUC + 신뢰구간

ROC 객체를 plot()에 넘기면 곡선이 그려지고, auc()로 면적을, ci.auc()로 신뢰구간을 구합니다.

# 곡선 그리기 (대각선 = 무작위 분류 기준선)
plot(roc_A,
     col = "#2E7D6B", lwd = 3,
     print.auc = TRUE,                # 곡선 위에 AUC 표시
     auc.polygon = TRUE,              # AUC 면적 음영
     auc.polygon.col = "#2E7D6B22",
     legacy.axes = TRUE,              # x축을 1-특이도(FPR)로 (0→1 방향)
     main = "바이오마커 A의 ROC 곡선")

# AUC 값과 95% 신뢰구간
auc(roc_A)                            # 점추정
ci.auc(roc_A)                         # DeLong 방식 95% CI (기본)
ci.auc(roc_A, method = "bootstrap")   # 부트스트랩 CI (대안)

Area under the curve: 0.7138
95% CI: 0.6427-0.7849 (DeLong)

legacy.axes = TRUE를 준 건 가로축을 익숙한 1−특이도(0에서 1로 증가)로 그리기 위해서입니다. 이걸 빼면 pROC 기본값은 가로축을 특이도(1에서 0)로 그려, 곡선 모양은 같아도 축 방향이 거꾸로 보여 헷갈립니다. AUC 0.71의 95% 신뢰구간이 0.64–0.78이고 하한이 0.5보다 충분히 크니, "이 바이오마커는 무작위(0.5)보다 유의하게 잘 가른다"고 말할 수 있습니다.

4. 최적 컷오프 — Youden 지수

AUC는 곡선 전체의 요약일 뿐, 실제로 양성/음성을 가르려면 컷오프 하나를 정해야 합니다. 가장 널리 쓰는 기준이 Youden 지수로, 민감도 + 특이도 − 1을 최대로 만드는 지점(ROC 곡선이 대각선에서 가장 멀리 떨어진 점)을 고릅니다.

# Youden 지수로 최적 컷오프와 그때의 민감도·특이도
coords(roc_A, x = "best",
       best.method = "youden",
       ret = c("threshold", "sensitivity", "specificity"))

  threshold sensitivity specificity
1  5.612471        0.71        0.66

이 컷오프(약 5.61)에서 점수가 그보다 높으면 양성으로 부르며, 그때 환자의 71%를 잡고(민감도) 정상인의 66%를 정상으로 맞춥니다(특이도). 임상 맥락에 따라 기준은 달라집니다. 놓치면 안 되는 병이라면 민감도를 우선해 컷오프를 낮추고, 불필요한 추가검사를 줄이려면 특이도를 높입니다. 특정 민감도(예: 90%)에서의 컷오프가 필요하면 다음처럼 직접 지정합니다.

# 민감도 0.90을 보장하는 컷오프에서의 특이도
coords(roc_A, x = 0.90, input = "sensitivity",
       ret = c("threshold", "specificity"))

5. 두 바이오마커 비교 — roc.test() (DeLong)

바이오마커 A와 B 중 어느 쪽이 더 잘 가르는지 묻는다면, 두 ROC를 겹쳐 그리고 AUC를 통계적으로 비교합니다. 같은 샘플에서 잰 두 점수는 짝지어진(paired) 데이터이므로, pROC는 상관을 반영하는 DeLong 검정을 씁니다.

# 바이오마커 B의 ROC도 만들기
roc_B <- roc(dat$status, dat$biomarker_B,
             levels = c("control", "case"), direction = "<")

# 두 곡선 겹쳐 그리기
plot(roc_A, col = "#2E7D6B", lwd = 3, legacy.axes = TRUE,
     main = "바이오마커 A vs B")
plot(roc_B, col = "#E76F51", lwd = 3, add = TRUE)   # add = TRUE로 겹쳐 그림
legend("bottomright",
       legend = c(sprintf("A (AUC %.3f)", auc(roc_A)),
                  sprintf("B (AUC %.3f)", auc(roc_B))),
       col = c("#2E7D6B", "#E76F51"), lwd = 3)

# AUC 차이의 통계적 유의성 (DeLong, paired)
roc.test(roc_A, roc_B, method = "delong")

	DeLong's test for two correlated ROC curves

data:  roc_A and roc_B
Z = -3.42, p-value = 0.000625
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
     0.7138      0.8462

B의 AUC (0.85)가 A (0.71)보다 높고, DeLong 검정의 p = 0.0006으로 그 차이가 우연이 아님이 확인됩니다. 즉 바이오마커 B가 두 그룹을 유의하게 더 잘 가릅니다. 세 개 이상을 비교할 땐 쌍마다 roc.test()를 돌리고, 여러 번 검정하므로 p-값을 보정합니다 (자세한 건 p-value와 FDR).

6. AUC 해석 기준 + ROC vs PR 곡선

AUC를 한 숫자로 읽을 때의 대략적 기준은 다음과 같습니다. 다만 분야·표본·목적에 따라 달라지는 관행적 눈금일 뿐, 절대 기준은 아닙니다.

한 가지 주의할 점이 있습니다. ROC는 클래스 불균형에 둔감합니다. 정상이 99%, 환자가 1%인 데이터에서는 거짓양성률(FPR)의 분모가 워낙 커서, 모델이 양성 예측을 남발해도 ROC 곡선이 낙관적으로 보입니다. 이런 희귀 양성·불균형 상황에서는 정밀도(precision)와 재현율(recall)을 축으로 쓰는 PR (precision-recall) 곡선이 모델 성능을 더 정직하게 드러냅니다.

요약하면, 클래스가 균형 잡혔으면 ROC/AUC로 충분하고, 양성이 드물면 PR 곡선을 함께 보는 게 안전합니다. PR 곡선은 PRROC나 precrec 같은 패키지로 그립니다.

7. 자주 발생하는 에러 & 오해

• AUC가 0.5보다 작게 나옴 (방향 뒤집힘) — 가장 흔한 실수입니다. 점수가 "낮을수록 질병"인데 pROC가 "높을수록 질병"으로 잡으면 곡선이 대각선 아래로 내려가 AUC < 0.5가 됩니다. 해결은 roc()에서 direction과 levels를 명시하는 것입니다. levels = c("control", "case")로 음성·양성 순서를, direction = "<"(양성 점수가 높음) 또는 direction = ">"(양성 점수가 낮음)로 방향을 못박으세요. AUC가 0.5 아래면 거의 항상 direction을 반대로 둔 것입니다.
• 클래스 라벨이 문자/숫자로 섞임 — response는 두 수준이어야 합니다. 0/1, "control"/"case" 어느 쪽이든 되지만, levels로 어느 쪽이 양성인지 명시하지 않으면 pROC가 알파벳·숫자 순으로 추정해 의도와 어긋날 수 있습니다. factor(..., levels = ...)로 순서를 고정하는 습관을 들이세요.
• AUC만 보고 컷오프를 무시 — AUC는 곡선 전체의 요약일 뿐, 실제 진단에는 컷오프 하나가 필요합니다. coords(roc, "best")로 최적 컷오프와 그때의 민감도·특이도를 함께 보고하세요. 임상에서는 AUC보다 "민감도 90%에서 특이도 얼마"가 더 중요할 때가 많습니다.
• 불균형 데이터에서 ROC를 과신 — 양성이 희귀하면 ROC가 낙관적으로 보입니다. 양성 비율이 한쪽으로 크게 쏠렸다면 PR 곡선을 함께 보세요(위 표).
• 학습·평가를 같은 데이터로 — 모델을 만든 데이터에서 그대로 AUC를 재면 과적합(overfitting)으로 부풀려집니다. 가능하면 학습/검증 데이터를 나누거나 교차검증(cross-validation)으로 AUC를 추정하세요. 이 글의 더미 예시는 개념 설명용이라 같은 데이터를 썼습니다.

References

1. Robin, X. et al. (2011). pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 12, 77. (pROC 원논문)
2. DeLong, E. R., DeLong, D. M., Clarke-Pearson, D. L. (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics, 44(3), 837–845. (DeLong 검정)
3. Youden, W. J. (1950). Index for rating diagnostic tests. Cancer, 3(1), 32–35. (Youden 지수)
4. Saito, T., Rehmsmeier, M. (2015). The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets. PLoS ONE, 10(3), e0118432. (ROC vs PR)

Pipette & Pipeline · A bio portfolio journal

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 12

반환된 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-04

logFC는 두 그룹의 평균 점수 차이, adj.P.Val은 다중검정 보정된 p-값입니다 (여러 경로를 동시에 검정하므로 보정이 필수 — 자세한 건 p-value와 FDR). A그룹에 신호를 준 경로들이 상위로 올라오면 분석이 의도대로 동작한 것입니다. 이렇게 GSVA + limma 조합은 "어떤 경로가 그룹 간에 다르게 켜졌나"를 경로 단위에서 직접 답합니다.

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는 샘플마다 경로 점수를 매겨 행렬을 만듭니다 (결과는 경로×샘플 행렬). 비교 대상이 정해져 있으면 전자, 샘플별 상태를 행렬로 다루고 싶으면 후자입니다.

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"으로 넣으면 분포 가정이 깨집니다. DESeq2 vst()·edgeR cpm(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() 두 단계로 바꾸면 해결됩니다.

References

1. Hänzelmann, S., Castelo, R., Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics, 14, 7. (GSVA 원논문)
2. 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)
3. Liberzon, A. et al. (2015). The Molecular Signatures Database Hallmark gene set collection. Cell Systems, 1(6), 417–425. (MSigDB Hallmark)
4. Barbie, D. A. et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature, 462, 108–112. (ssGSEA 원형)
5. 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

TL;DR — 두 그룹 비교는 t.test(), 세 그룹 이상은 aov() + TukeyHSD(), 두 변수의 관계는 cor.test(), 정규성⁠(normality)을 못 믿겠으면 wilcox.test()·kruskal.test() 같은 비모수⁠(non-parametric) 검정으로 갑니다. 검정을 고르는 기준은 단 세 가지 — 그룹이 몇 개인가 (2 vs 3 이상), 짝지어졌는가 (대응 vs 독립), 정규분포를 따르는가 (모수 vs 비모수)입니다.

이 글은 더미 발현 데이터를 직접 만들어, 가정 점검 (shapiro.test·var.test) → 2그룹 t-검정 (독립·대응·Welch) → wilcox.test → 일원분산분석 (ANOVA) + 사후검정 → kruskal.test → cor.test 순서로 복붙하면 도는 코드를 따라갑니다. 효과크기⁠(effect size)와 다중비교 보정도 함께 짚습니다. 모두 base R 내장 함수라 따로 설치할 패키지가 없습니다. (코드는 R 4.5 기준, 출력은 실제 실행값입니다.)

🔗 관련 글  ·  검정 결과의 p-값을 어떻게 다룰지 (다중검정 보정)는 p-value와 FDR, 다중검정 보정  ·  결과를 그림으로 옮기려면 ggplot2 기초  ·  차등발현 같은 대규모 검정은 DESeq2로 차등발현 분석

이 글에서 할 것

실험 데이터를 손에 쥐면 가장 먼저 떠오르는 질문이 "이 차이가 통계적으로 의미가 있나"입니다. 처리군이 대조군보다 발현이 높아 보이는데, 그게 진짜 효과인지 우연인지 가르는 일이죠. R에는 이걸 위한 검정 함수가 base에 거의 다 들어 있습니다 — t-검정 (t-test), 분산분석 (ANOVA), 상관⁠(correlation), 비모수 검정까지 한 줄짜리 함수로 돌아갑니다.

문제는 함수를 모르는 게 아니라 어떤 검정을 언제 쓰는지를 헷갈리는 경우가 많다는 점입니다. 그래서 이 글은 검정을 고르는 흐름부터 잡고, 더미 데이터로 하나씩 실행해 결과를 읽습니다. 데모는 단순합니다. 대조군·처리군 발현값을 만들고 (재현 가능하도록 시드 고정), 가정을 점검한 뒤, 2그룹 비교 → 3그룹 분산분석 → 상관 분석을 차례로 돌립니다. 출력에 찍히는 t·F·r·p-값을 어떻게 해석하는지까지 정리합니다.

# base R만으로 — 설치할 패키지 없음
t.test(value ~ group, data = expr)      # 두 그룹 평균 비교
aov(expr ~ geno, data = df3)            # 세 그룹 이상 (분산분석)
cor.test(x, y, method = "pearson")      # 두 변수의 상관

0. 검정 고르기 — 1분 흐름

검정을 고르는 결정은 세 갈래 질문으로 끝납니다. 순서대로 따라가면 됩니다.

• ① 그룹이 몇 개인가 — 두 그룹이면 t-검정 계열, 세 그룹 이상이면 분산분석 (ANOVA) 계열입니다. 세 그룹을 두 개씩 짝지어 t-검정을 반복하면 안 됩니다 (p-값이 부풀려지는 다중비교 문제 — 뒤에서 다룹니다).
• ② 짝지어졌는가 — 같은 개체를 처리 전후로 잰 값처럼 대응⁠(paired) 표본이면 짝을 묶어 검정하고, 서로 다른 개체에서 잰 독립⁠(independent) 표본이면 그냥 두 집단을 비교합니다. 이 구분을 놓치면 검정력이 떨어지거나 가정이 깨집니다.
• ③ 정규분포를 따르는가 — 데이터가 정규성⁠(normality)을 만족하면 평균을 비교하는 모수 (parametric) 검정 (t-검정·ANOVA), 그렇지 않거나 표본이 작고 순위 자료면 비모수⁠(non-parametric) 검정 (wilcox.test·kruskal.test)으로 갑니다. 비모수는 평균 대신 순위⁠(rank)를 비교합니다.

1. 사전 준비 — 데이터 만들기

설치할 패키지가 없습니다. 검정 함수는 전부 base R (stats 패키지)에 들어 있으니 바로 시작합니다. 먼저 대조군·처리군 발현값을 만듭니다. 처리군 평균을 일부러 조금 높여, 검정이 차이를 잡아내는지 봅니다.

set.seed(2025)                                 # 재현성 고정
n <- 8                                         # 군당 표본 수
ctrl <- rnorm(n, mean = 10.0, sd = 1.2)        # 대조군 발현
trt  <- rnorm(n, mean = 12.2, sd = 1.2)        # 처리군 발현 (평균 높임)

expr <- data.frame(
  value = c(ctrl, trt),
  group = factor(rep(c("Ctrl", "Trt"), each = n))   # 그룹 라벨
)
head(expr, 4)                                  # 구조 확인

     value group
1 10.74491  Ctrl
2 10.04277  Ctrl
3 10.92779  Ctrl
4 11.52699  Ctrl

검정에 들어가는 데이터 형태는 보통 긴 형식⁠(long format)입니다. 측정값 한 열 (value)과 그룹 라벨 한 열 (group)로 두면, value ~ group 같은 공식⁠(formula)으로 함수에 그대로 넘길 수 있습니다. 데이터 정리가 필요하면 dplyr·tidyr 데이터 정리를 참고하세요.

2. 가정 점검 — 정규성·등분산

모수 검정 (t-검정·ANOVA)은 두 가지를 가정합니다 — 데이터가 정규분포를 따르고, 비교하는 집단의 분산이 비슷하다(등분산)는 것입니다. 검정을 돌리기 전에 이걸 확인합니다.

shapiro.test(ctrl)             # 정규성 검정 (귀무가설: 정규분포다)
shapiro.test(trt)
var.test(value ~ group, data = expr)   # 등분산 검정 (F-검정)

	Shapiro-Wilk normality test
data:  ctrl
W = 0.94679, p-value = 0.6789

	F test to compare two variances
data:  value by group
F = 0.28378, num df = 7, denom df = 7, p-value = 0.1185

shapiro.test() (Shapiro-Wilk 검정)는 귀무가설이 "데이터가 정규분포를 따른다"입니다. 그래서 p > 0.05면 정규성을 기각하지 못함 — 즉 정규분포로 봐도 무방하다는 뜻입니다 (여기선 p = 0.68이라 통과). var.test()는 두 집단의 분산이 같은지 보는 F-검정이고, 역시 p > 0.05면 등분산⁠(equal variance)으로 봅니다 (여기선 p = 0.12). 셋 이상의 등분산은 Levene 검정 (car::leveneTest)이나 Bartlett 검정 (bartlett.test)을 씁니다.

3. 2그룹 비교 — t.test()

가정을 통과했으니 두 그룹 평균을 비교합니다. R의 t.test()는 기본이 Welch t-검정으로, 두 집단의 분산이 다를 수 있다고 가정합니다 (var.equal = FALSE가 기본값). 그래서 등분산이 의심스러워도 안전한 게 기본 동작입니다.

t.test(value ~ group, data = expr)         # Welch t-검정 (등분산 가정 안 함, 기본)

	Welch Two Sample t-test
data:  value by group
t = -3.7418, df = 10.677, p-value = 0.003426
alternative hypothesis: true difference in means between group Ctrl and group Trt is not equal to 0
95 percent confidence interval:
 -2.5912808 -0.6673877
sample estimates:
mean in group Ctrl  mean in group Trt
          10.48409           12.11343

출력을 읽는 법은 이렇습니다. t = -3.74는 검정통계량 t-값 (두 평균 차이를 표준오차로 나눈 값), df = 10.677은 자유도 (Welch라 정수가 아닙니다), p-value = 0.0034는 두 평균이 같다는 귀무가설 아래 이 차이가 나올 확률입니다. p < 0.05이므로 두 군의 평균이 다르다고 봅니다. 95% 신뢰구간⁠(confidence interval, CI) [-2.59, -0.67]이 0을 포함하지 않는 것도 같은 결론입니다 — 부호가 음수인 건 Ctrl − Trt라 대조군이 더 낮다는 뜻입니다.

등분산이 확실하면 고전적인 Student t-검정도 쓸 수 있습니다 (df가 정수로 바뀜).

t.test(value ~ group, data = expr, var.equal = TRUE)   # Student t-검정 (등분산 가정)
# t = -3.7418, df = 14, p-value = 0.002189

같은 개체를 처리 전후로 잰 대응⁠(paired) 표본이면 paired = TRUE를 켭니다. 짝을 묶어 차이를 검정하므로 개체 간 변동을 제거해 검정력이 올라갑니다.

# before/after는 같은 샘플을 두 번 잰 값 (순서·길이가 짝으로 맞아야 함)
t.test(after, before, paired = TRUE)       # 대응 t-검정
# t = 7.1658, df = 9, p-value = 5.275e-05

효과크기⁠(effect size)도 함께 봅니다. p-값은 "차이가 있는가"만 말하지 "얼마나 큰가"는 말하지 않기 때문입니다. 두 그룹 비교의 표준 효과크기는 Cohen's d — 평균 차이를 합동 표준편차 (pooled SD)로 나눈 값입니다.

# Cohen's d = (평균 차) / 합동 표준편차
m1 <- mean(trt); m0 <- mean(ctrl)
s_p <- sqrt(((n-1)*var(trt) + (n-1)*var(ctrl)) / (2*n - 2))
(m1 - m0) / s_p                            # 1.87 → 큰 효과 (관행: 0.2 작음·0.5 중간·0.8 큼)

4. 비모수 2그룹 — wilcox.test()

정규성을 만족하지 못하거나 표본이 작고 이상치가 있으면 비모수 검정으로 갑니다. 두 그룹용은 wilcox.test() — 독립 표본이면 Mann-Whitney U 검정, 대응 표본 (paired = TRUE)이면 Wilcoxon 부호순위 검정입니다. 평균 대신 순위⁠(rank)를 비교하므로 분포 가정이 느슨합니다.

wilcox.test(value ~ group, data = expr)    # Mann-Whitney U (독립·비모수)
# paired = TRUE 면 Wilcoxon 부호순위 검정

	Wilcoxon rank sum exact test
data:  value by group
W = 5, p-value = 0.002953

출력의 W는 순위합 기반 검정통계량이고, 해석은 p-값으로 같습니다 (p < 0.05면 두 군의 분포 위치가 다름). 모수 t-검정과 결론이 대체로 일치하지만, 비모수는 평균이 아니라 분포의 위치 이동⁠(location shift)을 봅니다. 동점⁠(tie)이 있으면 정확 p-값 대신 정규근사를 쓴다는 경고가 뜰 수 있는데, 해석에는 지장이 없습니다.

5. 3그룹 이상 — aov() + TukeyHSD()

그룹이 셋 이상이면 t-검정을 짝마다 반복하는 대신 일원분산분석 (one-way ANOVA)으로 "적어도 한 군은 다른가"를 한 번에 검정합니다. R에서는 aov()로 모형을 적합하고 summary()로 표를 봅니다.

# 세 유전형: WT(야생형)·KO(녹아웃)·OE(과발현)
fit <- aov(expr ~ geno, data = df3)        # 일원분산분석
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)
geno         2  29.52  14.760   19.51 1.63e-05 ***
Residuals   21  15.89   0.757

분산분석표에서 핵심은 F value = 19.51과 Pr(>F) = 1.63e-05입니다. F-값은 군 간 분산을 군 내 분산으로 나눈 비로, 클수록 군 차이가 두드러집니다. p < 0.05이므로 "세 군의 평균이 모두 같다"는 귀무가설을 기각합니다. 다만 ANOVA는 어느 군이 다른지는 말해 주지 않습니다 — 그래서 사후검정⁠(post-hoc test)이 필요합니다.

TukeyHSD()는 모든 군 쌍을 비교하되, 다중비교에 따른 p-값 부풀림을 자체적으로 보정합니다 (Tukey의 HSD — honestly significant difference).

TukeyHSD(fit)                              # 모든 쌍 비교 + 다중비교 보정

$geno
            diff       lwr      upr     p adj
KO-WT 0.06274305 -1.033493 1.158979 0.9885979
OE-WT 2.38340203  1.287166 3.479638 0.0000558
OE-KO 2.32065899  1.224423 3.416895 0.0000779

표의 각 행이 한 쌍의 비교입니다. diff는 평균 차, lwr·upr는 그 차이의 95% 신뢰구간, p adj는 보정된 p-값입니다. 여기서는 OE-WT와 OE-KO가 p < 0.001로 유의하지만 KO-WT는 p = 0.99로 차이가 없습니다 — 즉 OE만 도드라지고 KO와 WT는 비슷하다는 뜻입니다. ANOVA의 전체 p-값만 봤다면 놓쳤을 구조를, 사후검정이 짚어 줍니다.

6. 비모수 3그룹 — kruskal.test()

세 군 이상인데 정규성을 못 믿겠으면 ANOVA의 비모수 대응인 kruskal.test() (Kruskal-Wallis 검정)를 씁니다. 역시 순위 기반으로 "적어도 한 군은 다른가"를 검정합니다.

kruskal.test(expr ~ geno, data = df3)      # Kruskal-Wallis (비모수 ANOVA)

	Kruskal-Wallis rank sum test
data:  expr by geno
Kruskal-Wallis chi-squared = 15.365, df = 2, p-value = 0.0004608

p < 0.05이므로 군 간 차이가 있습니다. 어느 쌍이 다른지 보려면 사후검정으로 Dunn 검정(dunn.test::dunn.test 또는 FSA::dunnTest)을 쓰고, 거기서 다중비교 보정을 적용합니다 — TukeyHSD가 ANOVA의 짝이라면, Dunn 검정이 Kruskal-Wallis의 짝입니다.

7. 상관 — cor.test()

두 연속 변수가 함께 움직이는지 (예: 유전자 A 발현 ↔ 유전자 B 발현)는 상관⁠(correlation)으로 봅니다. cor.test()는 상관계수 r과 그 p-값을 함께 줍니다.

cor.test(x, y, method = "pearson")         # 피어슨: 선형 관계 (정규성 가정)

	Pearson's product-moment correlation
data:  x and y
t = 9.3529, df = 38, p-value = 2.119e-11
cor
0.8349581

cor = 0.83이 피어슨 상관계수 r입니다 (−1~+1, 0이면 무관, 부호는 방향). p < 0.05이므로 두 변수의 상관이 0이 아니라고 봅니다. 피어슨 (Pearson)은 선형⁠(linear) 관계를 재고 정규성을 가정하므로, 관계가 곡선이거나 이상치가 있으면 순위 기반의 스피어만 (Spearman)으로 바꿉니다.

cor.test(x, y, method = "spearman")        # 스피어만: 단조⁠(순위) 관계, 비모수
# rho = 0.737, p-value = 3.419e-07

스피어만은 값 대신 순위로 상관을 재므로 (r 대신 ρ, rho), 단조⁠(monotonic)이기만 하면 곡선 관계도 잡고 이상치에 덜 휘둘립니다. 관계의 모양을 모르겠으면 둘 다 보고 비교하면 됩니다. 관계를 식으로 모형화하려면 회귀⁠(regression)로 넘어갑니다 — 한 줄이면 됩니다.

lm(y ~ x)                                  # 선형회귀 (기울기·절편 추정)
summary(lm(y ~ x))$coefficients            # 계수·표준오차·t·p

8. 검정 선택 표

9. 자주 발생하는 에러 & 오해

• 세 그룹을 t-검정으로 짝마다 비교 — 가장 흔한 실수입니다. 비교를 반복하면 위양성이 쌓이니, 3군 이상은 aov() + TukeyHSD() (또는 kruskal.test() + Dunn)로 한 번에 검정하고 사후 보정을 받으세요.
• 독립 표본에 paired = TRUE — 서로 다른 개체인데 대응 검정을 켜면 결과가 무의미합니다. 같은 대상의 전/후일 때만 대응입니다. 반대로 진짜 대응 자료를 독립으로 검정하면 검정력을 버립니다.
• 정규성·등분산을 점검하지 않고 t-검정 — 작은 표본에서 분포가 치우치면 t-검정의 p-값을 믿기 어렵습니다. shapiro.test·qqnorm으로 보고, 애매하면 wilcox.test로 가세요. R의 t.test는 기본이 Welch라 등분산을 가정하지 않는 점은 다행입니다.
• 다중비교 보정 누락 — 여러 유전자·여러 조건을 한꺼번에 검정하면 p-값을 그대로 쓰면 안 됩니다. p.adjust(pvals, method = "BH")로 보정하세요 (BH = Benjamini-Hochberg, FDR 제어). 예: p.adjust(c(0.001, 0.008, 0.02, 0.04, 0.31), "BH") → 0.005, 0.020, 0.033, 0.050, 0.310. 자세한 건 p-value와 FDR에서 다룹니다.
• 상관을 인과로 해석 — r이 크고 p가 작아도 "A가 B를 일으킨다"는 뜻은 아닙니다. 상관은 동반 변화일 뿐, 인과는 실험 설계로 따로 보여야 합니다.
• cor.test가 NA를 반환 — 입력에 NA가 있으면 결과가 NA로 나옵니다. use = "complete.obs"를 주거나 결측을 먼저 제거하세요.

References

1. R Core Team. (2024). R: A language and environment for statistical computing. (t.test·aov·cor.test·wilcox.test·kruskal.test·p.adjust 등 stats 패키지 기본 함수)
2. Student. (1908). The probable error of a mean. Biometrika, 6(1), 1–25. (t-검정의 원형)
3. Welch, B. L. (1947). The generalization of "Student's" problem when several different population variances are involved. Biometrika, 34(1–2), 28–35. (Welch t-검정 — t.test의 기본)
4. Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate. Journal of the Royal Statistical Society B, 57(1), 289–300. (p.adjust(method = "BH") — 다중검정 보정)
5. Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Erlbaum. (Cohen's d — 효과크기 관행)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 발현 히트맵을 가장 빠르게 그리는 패키지가 pheatmap (pretty heatmaps)입니다. 핵심은 단 한 줄, pheatmap(mat). 여기에 scale = "row" (유전자별 z-score, 행 표준화)만 더하면 발현이 큰 유전자에 가려져 있던 패턴이 한눈에 드러납니다. 행·열 계층 클러스터링⁠(hierarchical clustering), 샘플 그룹 색띠⁠(annotation), 색·컷·저장까지 인자 몇 개로 끝납니다.

이 글은 더미 발현행렬을 직접 만들어, 기본 히트맵 → 행 표준화 → 클러스터링 → 샘플 주석 → 색·cutree·저장 순서로 복붙하면 도는 코드를 따라갑니다. pheatmap과 ComplexHeatmap의 차이, 자주 터지는 실수도 정리합니다. (코드는 pheatmap 1.0.12 기준입니다.)

🔗 관련 글  ·  더 정밀한 히트맵이 필요하면: ComplexHeatmap으로 발현 히트맵 (고급)  ·  히트맵에 올릴 차등발현 결과는 DESeq2로 차등발현 분석  ·  샘플 군집을 먼저 보려면 PCA로 샘플 QC

이 글에서 그릴 것

차등발현 분석이 끝나면 보통 다음 그림이 발현 히트맵입니다. 상위 변동 유전자 수십 개를 골라, 어떤 유전자군이 어떤 샘플에서 함께 오르내리는지 한 장으로 보여 주죠. R에서 이걸 가장 빠르게 그리는 도구가 pheatmap입니다 — 함수 하나에 인자 몇 개만 얹으면 z-score 히트맵에 주석띠까지 완성됩니다.

데모는 단순합니다. 대조군 3개·처리군 3개, 유전자 30개짜리 더미 발현행렬을 만들고⁠(재현 가능하도록 시드 고정), scale = "row"로 유전자별 z-score를 매긴 뒤, 샘플 그룹을 색띠로 붙인 히트맵을 그립니다. 최종 결과물은 빨강⁠(평균보다 높음)·파랑⁠(낮음)으로 발현 패턴이 정리되고, 행·열이 클러스터링된 히트맵입니다.

library(pheatmap)
pheatmap(mat,
         scale = "row",                 # ★유전자별 z-score — 발현 패턴 비교의 핵심
         annotation_col = anno_col)      # 샘플 그룹 색띠

0. 왜 행 표준화⁠(scale = "row")인가 — 1분 직관

히트맵을 그릴 때 가장 흔한 실수가 스케일링을 건너뛰는 것입니다. 원시 발현값⁠(vst 또는 log 발현)을 그대로 그리면, 절대 발현이 큰 유전자 몇 개가 색 범위를 독차지합니다. 항존유전자⁠(housekeeping gene)는 늘 빨갛고, 정작 보고 싶은 저발현 전사인자의 미세한 변화는 죄다 같은 색으로 뭉개지죠.

행 표준화⁠(z-score)는 이 문제를 해결합니다. 유전자⁠(행)마다 평균을 빼고 표준편차로 나눠, 모든 유전자를 평균 0·표준편차 1의 같은 잣대로 맞춥니다. 그러면 절대 발현량은 사라지고 "이 샘플에서 이 유전자가 자기 평균보다 높은가 낮은가"라는 상대 패턴만 남습니다. 발현량의 크기가 아니라 패턴을 비교하는 게 히트맵의 목적이므로, 이게 핵심입니다.

pheatmap의 편한 점은 이 표준화를 scale = "row" 한 인자로 처리한다는 것입니다. (참고로 ComplexHeatmap에서는 t(scale(t(mat)))처럼 직접 z-score를 계산해 넘기는 경우가 많은데, pheatmap은 내부에서 알아서 합니다.)

1. 준비물 — 설치

pheatmap은 CRAN 패키지라 설치가 가볍습니다. 상위 변동 유전자를 고를 때 분산을 빠르게 계산하려면 matrixStats (rowVars)도 함께 있으면 편합니다.

install.packages("pheatmap")     # 최초 1회
install.packages("matrixStats")  # rowVars로 상위 변동 유전자 선택 (선택)

library(pheatmap)                # 히트맵 함수
library(matrixStats)             # rowVars (없으면 apply로 대체 가능)

pheatmap이 받는 입력은 유전자 × 샘플 수치행렬입니다. 행이 유전자, 열이 샘플이고, 값은 보통 vst (분산안정화 변환) 또는 log2 정규화 발현값입니다. 카운트 원값을 그대로 넣지는 않습니다 — 분포가 한쪽으로 쏠려 색이 무너지니, 먼저 정규화·변환을 거친 행렬을 씁니다.

2. 더미 데이터 만들기 + 상위 변동 유전자 고르기

실제로는 DESeq2의 vst()나 assay()에서 행렬을 가져오지만, 여기서는 재현 가능한 더미 행렬을 만듭니다. 대조군 3개·처리군 3개에, 일부 유전자는 군 간 차이가 나도록 설계합니다.

set.seed(42)                                  # 재현성 고정
genes   <- paste0("Gene", sprintf("%02d", 1:30))
samples <- c(paste0("Ctrl", 1:3), paste0("Trt", 1:3))

# 기본 발현 (log 스케일 가정): 평균 8, 약간의 노이즈
mat <- matrix(rnorm(30 * 6, mean = 8, sd = 1),
              nrow = 30, dimnames = list(genes, samples))

# 앞쪽 10개 유전자는 처리군에서 상향 (군 차이를 일부러 부여)
mat[1:10, 4:6] <- mat[1:10, 4:6] + 2.5

dim(mat)                                       # 30 × 6
mat[1:3, ]                                      # 값 확인

[1] 30  6
          Ctrl1    Ctrl2    Ctrl3    Trt1     Trt2     Trt3
Gene01  9.371    7.435    8.363    10.632   10.404   11.298
Gene02  7.435    8.116    7.773    10.111   10.231   10.512
Gene03  8.363    8.718    8.575    11.286   10.490   10.954

실전에서는 유전자가 수만 개라 전부 그릴 수 없습니다. 변동이 큰 상위 유전자⁠(보통 50~100개)만 골라야 패턴이 보입니다. 분산⁠(또는 표준편차) 기준으로 정렬해 위에서 잘라냅니다.

# 분산이 큰 순서로 상위 유전자 인덱스 (실전: top 50~100)
rv  <- rowVars(mat)                            # matrixStats — 유전자별 분산
top <- order(rv, decreasing = TRUE)[1:20]      # 데모라 20개만
mat_top <- mat[top, ]

# matrixStats가 없으면:
# rv <- apply(mat, 1, var)

3. 기본 히트맵 — pheatmap(mat)

먼저 아무 옵션 없이 그려 봅니다. pheatmap은 기본으로 행·열 모두 클러스터링하고, 값에 따라 색을 입힙니다.

pheatmap(mat_top)                              # 기본: 행·열 클러스터링 + 원시값 색

이 그림에는 문제가 있습니다. 스케일링을 안 했으니, 발현이 큰 유전자 행이 진하게 깔리고 작은 유전자는 흐릿합니다. 색이 절대 발현량을 따라가, 군 간 패턴이 잘 안 보이죠. 다음 단계에서 행 표준화를 켭니다.

4. 행 표준화 — scale = "row" (핵심)

scale = "row" 하나만 추가합니다. 이제 각 유전자⁠(행)가 z-score로 변환돼, 자기 평균 대비 높낮이만 색으로 표현됩니다.

pheatmap(mat_top,
         scale = "row")                        # ★유전자별 z-score (행 표준화)

이 한 줄로 그림이 완전히 달라집니다. 처리군에서 올린 앞쪽 유전자들이 일관되게 빨강, 대조군에서 파랑으로 뚜렷이 갈립니다. scale의 값은 세 가지 — "row"⁠(유전자별, 발현 히트맵의 표준), "column"⁠(샘플별, 드묾), "none"⁠(원값 그대로, 기본값)입니다. 발현 패턴을 볼 때는 거의 항상 "row"입니다.

5. 클러스터링 제어 — cluster_rows / cluster_cols

pheatmap은 기본으로 행·열을 계층 클러스터링해 덴드로그램⁠(dendrogram)을 그립니다. 거리·연결 방식을 바꾸거나, 클러스터링을 끌 수 있습니다.

pheatmap(mat_top,
         scale = "row",
         cluster_rows = TRUE,                  # 행⁠(유전자) 클러스터링
         cluster_cols = TRUE,                  # 열⁠(샘플) 클러스터링
         clustering_distance_rows = "correlation",  # 상관 거리 (발현 패턴엔 흔히 유용)
         clustering_distance_cols = "euclidean",    # 유클리드 거리 (기본)
         clustering_method = "complete")       # 연결법: complete·average·ward.D2 등

행 클러스터링은 함께 오르내리는 유전자⁠(공발현 유전자군)를 묶고, 열 클러스터링은 발현이 비슷한 샘플을 묶습니다. 발현 패턴의 모양을 보고 싶을 때는 거리를 "correlation"⁠(1 − 상관계수)으로 두는 경우가 많습니다 — 절대 크기가 아니라 패턴의 유사도로 묶이기 때문입니다. 샘플 순서를 원본대로 두고 싶으면 cluster_cols = FALSE로 끄면 됩니다.

6. 샘플 주석 — annotation_col + annotation_colors

히트맵 위에 샘플이 어느 그룹인지 색띠를 붙이면 해석이 훨씬 쉬워집니다. 행 이름이 샘플 ID와 일치하는 데이터프레임을 만들어 annotation_col에 넘깁니다. 유전자⁠(행) 쪽 주석은 annotation_row로 같은 방식입니다.

# 행 이름 = 샘플 ID (mat의 열 이름과 정확히 일치해야 함)
anno_col <- data.frame(
  Condition = rep(c("Control", "Treated"), each = 3),
  row.names = samples
)

# 주석 색 직접 지정 (선택)
anno_colors <- list(
  Condition = c(Control = "#6BAB97", Treated = "#E76F51")
)

pheatmap(mat_top,
         scale = "row",
         annotation_col    = anno_col,         # 열 위 색띠 (샘플 그룹)
         annotation_colors = anno_colors)      # 색띠 색 지정

annotation_col의 행 이름이 행렬 열 이름과 어긋나면 주석이 통째로 안 붙거나 NA로 나옵니다 — all(rownames(anno_col) == colnames(mat_top))로 먼저 확인하세요. 여러 변수를 넣으면⁠(예: Condition, Batch) 색띠가 여러 줄로 쌓입니다.

7. 색·컷·간격·저장 — 마감 손질

마지막으로 발표용 그림을 다듬습니다. 색 팔레트, 덴드로그램 자르기, 샘플 그룹 사이 간격, 그리고 파일 저장입니다.

# 파랑–흰–빨강 발산형 팔레트
my_color <- colorRampPalette(c("#2166AC", "white", "#B2182B"))(100)

pheatmap(mat_top,
         scale = "row",
         color = my_color,                     # 발산형 색 (z-score에 적합)
         breaks = seq(-2, 2, length.out = 101),# 색 경계: −2~2로 고정 (이상치 포화 방지)
         cutree_rows = 2,                       # 행 덴드로그램을 2개 군으로 절단
         cutree_cols = 2,                       # 열도 2개 군으로 절단
         show_rownames = TRUE,                  # 유전자 이름 표시 (많으면 FALSE 권장)
         show_colnames = TRUE,
         fontsize_row = 8,                      # 행 글자 크기
         main = "상위 변동 유전자 z-score 히트맵")

각 인자의 쓸모는 이렇습니다. color + breaks는 z-score의 색 범위를 −2~2로 고정해, 극단값 하나가 색을 다 먹어 버리는 일을 막습니다. cutree_rows·cutree_cols는 덴드로그램을 k개 군으로 잘라 사이에 흰 간격을 넣어 군집을 시각적으로 분리합니다⁠(클러스터링이 켜져 있어야 동작). show_rownames = FALSE는 유전자가 많아 이름이 겹칠 때 필수입니다. 군 사이를 수동으로 띄우려면 클러스터링을 끄고 gaps_col = c(3)처럼 위치를 직접 줄 수도 있습니다.

저장은 두 가지 길이 있습니다.

# 방법 A — filename 인자 (가장 간단, 확장자로 포맷 결정)
pheatmap(mat_top, scale = "row", annotation_col = anno_col,
         filename = "heatmap.png", width = 6, height = 7)

# 방법 B — ggsave (pheatmap은 grid 객체라 $gtable을 넘김)
p <- pheatmap(mat_top, scale = "row", silent = TRUE)
ggplot2::ggsave("heatmap.pdf", p$gtable, width = 6, height = 7)

filename을 주면 화면에 띄우지 않고 바로 파일로 떨어집니다⁠(확장자가 .png·.pdf·.tiff면 그 포맷). pheatmap의 반환값은 grid 그래픽 객체라, ggplot의 ggsave로 저장할 때는 $gtable을 넘깁니다.

8. 주요 인자 정리

9. pheatmap vs ComplexHeatmap — 언제 무엇을

둘 다 발현 히트맵을 그리지만 결이 다릅니다. pheatmap은 함수 하나로 z-score·클러스터링·주석을 끝내는 빠르고 쉬운 도구라, 탐색 단계나 표준적인 발현 히트맵에 적합합니다. ComplexHeatmap (Gu et al. 2016)은 여러 히트맵을 가로로 잇고, 행을 임의 기준으로 분할하고, 복잡한 주석⁠(막대그래프·점)을 얹는 정밀하고 유연한 도구라, 그림 한 장에 많은 정보를 담는 출판용 figure에 강합니다. 빠르게 패턴만 보려면 pheatmap, 정교한 조합 그림이 필요하면 ComplexHeatmap으로 넘어가면 됩니다⁠(ComplexHeatmap 글 참고).

10. 자주 발생하는 에러 & 해결

• 고발현 유전자만 새빨갛고 나머지는 패턴이 안 보임 — 스케일링을 안 한 경우입니다. scale = "row"를 켜 유전자별 z-score로 바꾸세요. 발현 히트맵의 1순위 점검 항목입니다.
• 유전자 이름이 겹쳐 까맣게 떡짐 — 행이 너무 많은데 show_rownames = TRUE라서입니다. show_rownames = FALSE로 끄거나, 상위 변동 유전자를 더 적게⁠(예: 50개) 고르세요.
• Error: 'x' must be numeric 또는 NA로 색이 깨짐 — 행렬에 문자열 열이 섞였거나 NA가 있습니다. mat <- as.matrix(mat)로 수치행렬화하고, mat <- mat[complete.cases(mat), ]로 NA 행을 제거하세요. 표준편차가 0인 유전자⁠(모든 값이 같음)도 scale = "row"에서 NaN을 만드니 함께 거릅니다.
• 주석 색띠가 안 붙거나 회색 NA로 나옴 — annotation_col의 행 이름이 행렬 열 이름과 다릅니다. all(rownames(anno_col) == colnames(mat))가 TRUE인지 확인하세요⁠(순서까지 일치해야 안전).
• cutree_rows가 무시됨 — cluster_rows = FALSE로 클러스터링을 끈 상태에서는 덴드로그램이 없어 자를 수 없습니다. cutree는 해당 축 클러스터링이 켜져 있을 때만 동작합니다.
• ggsave가 빈 파일을 만듦 — pheatmap 결과를 그대로 넘겨서입니다. 반환 객체의 $gtable을 넘기거나, 가장 간단하게는 filename = 인자를 쓰세요.

References

1. Kolde, R. (2019). pheatmap: Pretty Heatmaps. R package version 1.0.12. CRAN. (pheatmap 패키지·scale·annotation_col·cutree)
2. Gu, Z., Eils, R., Schlesner, M. (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 32(18), 2847–2849. (pheatmap 대비 정밀·복합 히트맵)
3. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. (vst — 히트맵 입력 변환)
4. R Core Team. (2024). R: A language and environment for statistical computing. (colorRampPalette·scale·hclust 기본 함수)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 발현행렬·메타데이터·DEG 표를 다룰 때 실제 시간의 대부분은 통계가 아니라 데이터 정리에 들어갑니다. tidyverse의 dplyr (행·열 조작)과 tidyr (형식 변환)이 그 일을 맡죠. dplyr의 다섯 동사 — filter (행)·select (열)·mutate (파생열)·arrange (정렬)·summarise+group_by (집계) — 와 tidyr의 pivot_longer/pivot_wider (wide↔long), 파이프 |>만 손에 익으면 대부분의 정리가 한 흐름으로 이어집니다.

이 글은 더미 발현 데이터를 코드로 만들어 결합·집계·필터·형식 변환까지 복붙하면 그대로 도는 예시로 따라갑니다. ggplot2 글과 함께 'tidyverse 기초'를 이루는 짝입니다. (코드는 dplyr 1.1+ / tidyr 1.3+ 기준으로 실제 실행해 확인했습니다.)

🔗 관련 글  ·  정리한 데이터를 그림·분석으로: ggplot2 기초  ·  DESeq2로 차등발현 분석  ·  GEOquery로 공개데이터 받기

이 글에서 할 것

분석 코드를 짜다 보면 정작 DESeq2나 ggplot2를 부르기 전에, 표 모양을 맞추는 데 대부분의 시간을 씁니다. 발현행렬은 유전자가 행·샘플이 열인데 그림을 그리려면 long (긴 형식)으로 바꿔야 하고, 메타데이터의 그룹 정보는 발현값과 다른 표에 들어 있어 결합⁠(join)해야 하죠. DEG 결과에서는 유의한 유전자만 골라 정렬해야 합니다.

이 글은 그 정리 과정을 dplyr·tidyr로 처리합니다. 유전자 × 샘플 발현행렬과 샘플 메타데이터를 코드로 만든 뒤 — ① 두 표를 결합하고 ② 그룹별 평균 발현을 구하고 ③ DEG 표에서 유의 유전자를 추리고 ④ wide 발현행렬을 long으로 바꿔 ggplot2에 넘길 준비까지, 한 단계씩 따라갑니다.

# 이 글이 반복해서 쓰는 한 패턴 — 데이터를 파이프로 흘려보낸다
expr |>
  pivot_longer(-gene, names_to = "sample", values_to = "expr") |>   # wide → long
  left_join(meta, by = "sample") |>                                 # 메타 결합
  summarise(mean_expr = mean(expr), .by = c(gene, group))           # 그룹별 평균

0. tidy data — 정리의 목표가 뭔가

dplyr·tidyr가 전제하는 데이터 모양이 tidy data (정돈된 데이터)입니다. 규칙은 세 가지뿐입니다 — 한 변수는 한 열, 한 관측은 한 행, 한 종류의 관측 단위는 한 표. 발현 데이터로 말하면 "어떤 유전자가 어떤 샘플에서 얼마로 발현됐다"가 한 관측이니, tidy 형태에서는 gene·sample·expr 세 열에 한 측정이 한 행으로 들어갑니다 (Wickham 2014).

그런데 우리가 흔히 받는 발현행렬은 샘플이 열로 펼쳐진 wide (넓은 형식)입니다. 사람이 보기엔 편하지만 ggplot2·통계 함수는 tidy long을 좋아하죠. 그래서 tidyr의 pivot_longer로 wide를 long으로 눕히는 작업이 거의 항상 첫 단계가 됩니다. 정리의 목표는 결국 "지금 표를 tidy 형태로 만드는 것"이라고 보면 됩니다.

1. 사전 준비 — 설치와 예시 데이터

dplyr·tidyr는 둘 다 tidyverse 묶음에 들어 있습니다. tidyverse를 깔면 dplyr (데이터 조작)·tidyr (형식 변환)·ggplot2 (시각화) 등이 한 번에 들어와 편합니다.

install.packages("tidyverse")   # dplyr + tidyr + ggplot2 ... 한 번에 (최초 1회)

library(dplyr)   # 행·열 조작: filter·select·mutate·arrange·summarise·group_by
library(tidyr)   # 형식 변환: pivot_longer·pivot_wider·separate·drop_na

이제 예시 데이터를 코드로 만듭니다. 유전자 6개 × 샘플 6개 발현행렬, 샘플별 그룹을 담은 메타데이터, 그리고 DEG 결과 표 — 세 개를 더미로 생성합니다. set.seed로 난수를 고정해 누가 돌려도 같은 값이 나오게 합니다.

set.seed(1)

# 발현행렬: 유전자⁠(행) × 샘플⁠(열), 첫 열은 gene 이름 (wide 형식)
genes <- paste0("Gene", 1:6)
expr <- data.frame(
  gene = genes,
  S1 = round(rnorm(6, 8, 2), 2), S2 = round(rnorm(6, 8, 2), 2),
  S3 = round(rnorm(6, 8, 2), 2), S4 = round(rnorm(6, 9, 2), 2),
  S5 = round(rnorm(6, 9, 2), 2), S6 = round(rnorm(6, 9, 2), 2)
)

# 샘플 메타데이터: 샘플마다 그룹·배치 (별도 표)
meta <- data.frame(
  sample = paste0("S", 1:6),
  group  = rep(c("control", "treated"), each = 3),
  batch  = rep(c("b1", "b2"), times = 3)
)

# DEG 결과 표: logFC와 보정 p값 (DESeq2/limma 출력을 흉내)
deg <- data.frame(
  gene       = genes,
  logFC      = c(2.4, -1.8, 0.3, 1.2, -0.5, 3.1),
  adj.P.Val  = c(0.001, 0.004, 0.62, 0.03, 0.41, 0.0008)
)

head(expr, 3)

   gene    S1    S2    S3    S4    S5    S6
1 Gene1  6.75  9.84  9.62 11.36  8.39  9.39
2 Gene2  8.37  8.16  9.74  8.88  9.78 11.93
3 Gene3  6.33  9.49  6.85  9.07  9.59  9.94

2. 파이프 — |> vs %>%

본격적으로 들어가기 전에 파이프 (pipe)부터 짚겠습니다. 파이프는 왼쪽 결과를 오른쪽 함수의 첫 인자로 넘기는 연산자라, f(g(x)) 같은 중첩을 x |> g() |> f()처럼 읽는 순서대로 펼 수 있습니다.

R에는 파이프가 둘 있습니다. base R에 내장된 |> (R 4.1+)와, magrittr 패키지의 %>% (tidyverse가 오래 써 온 것)입니다. 일상적인 정리 작업에서는 거의 똑같이 동작하니, 새로 배운다면 의존성이 없는 |>를 권합니다.

# 아래 셋은 같은 일을 한다 — head(arrange(deg, adj.P.Val), 3)
deg |> arrange(adj.P.Val) |> head(3)    # base 파이프 (권장)
deg %>% arrange(adj.P.Val) %>% head(3)   # magrittr 파이프 (tidyverse 전통)

3. dplyr 다섯 동사 — filter·select·mutate·arrange

dplyr는 데이터 조작을 다섯 개의 동사 (verb) 함수로 정리합니다. 전부 "데이터프레임을 받아 데이터프레임을 돌려준다"는 규칙이 같아, 파이프로 줄줄이 이을 수 있습니다.

deg |>
  filter(adj.P.Val < 0.05, abs(logFC) > 1) |>   # 행 고르기: 유의 + 효과크기 큰 DEG
  mutate(direction = if_else(logFC > 0, "up", "down")) |>   # 파생열 추가
  select(gene, logFC, direction) |>             # 열 고르기 (순서도 정렬)
  arrange(desc(abs(logFC)))                      # |logFC| 내림차순 정렬

   gene logFC direction
1 Gene6   3.1        up
2 Gene1   2.4        up
3 Gene2  -1.8      down
4 Gene4   1.2        up

한 줄씩 뜯어보면 — filter는 조건을 만족하는 행만 남기고 (쉼표로 나열하면 AND), mutate는 기존 열로 새 열을 계산하며, select는 필요한 열만 골라 순서까지 정합니다. arrange는 정렬인데 desc()로 감싸면 내림차순입니다. 여기서 adj.P.Val < 0.05 & abs(logFC) > 1은 DEG를 추릴 때 가장 흔히 쓰는 조건이죠.

여기에 더해 자주 쓰는 보조 동사도 있습니다.

deg |> count(adj.P.Val < 0.05)        # 조건별 개수 세기
meta |> distinct(group)               # 중복 제거한 고유값
deg |> slice_max(logFC, n = 2)        # logFC 상위 2개 행

4. group_by + summarise — 그룹별 집계

데이터를 그룹으로 묶어 요약 (aggregate)하는 건 정리에서 가장 자주 하는 일입니다. group_by로 묶고 summarise로 집계하면 "그룹마다 한 줄"로 줄어듭니다. 먼저 발현행렬을 long으로 바꿔야 그룹 집계가 자연스러운데, long 변환은 6절에서 자세히 보고 여기서는 결과만 씁니다.

# expr를 long으로 바꾼 뒤 메타 결합 (자세한 변환은 6절)
long <- expr |>
  pivot_longer(-gene, names_to = "sample", values_to = "expr") |>
  left_join(meta, by = "sample")

# 유전자 × 그룹별 평균·표준편차 발현
long |>
  group_by(gene, group) |>
  summarise(mean_expr = mean(expr),
            sd_expr   = sd(expr),
            n         = n(),       # 그룹 내 샘플 수
            .groups   = "drop")    # 집계 후 그룹 해제

# A tibble: 12 × 5
   gene  group   mean_expr sd_expr     n
   <chr> <chr>       <dbl>   <dbl> <int>
 1 Gene1 control      8.74   1.69      3
 2 Gene1 treated      9.71   1.49      3
 3 Gene2 control      8.76   0.860     3
 ...

summarise 안에서 mean·sd·n 여러 통계를 한 번에 계산했습니다. 집계가 끝나면 .groups = "drop"으로 그룹을 풀어 줘야 뒤 단계에서 의도치 않은 그룹 연산을 막을 수 있습니다 (예전 코드의 ungroup()과 같은 역할). dplyr 1.1+에서는 summarise(..., .by = c(gene, group))처럼 .by 인자로 그룹을 그 자리에서만 적용할 수도 있어, group_by + ungroup을 한 줄로 대체합니다.

여러 열에 같은 요약을 한꺼번에 적용할 때는 across가 깔끔합니다.

# 모든 숫자형 샘플 열을 z-점수로 정규화 (mutate + across)
expr |>
  mutate(across(where(is.numeric), ~ (.x - mean(.x)) / sd(.x)))   # 열마다 평균 0·분산 1

5. join — 발현값과 메타데이터 결합

발현값과 그룹 정보가 다른 표에 있을 때, 공통 열을 기준으로 붙이는 게 join (결합)입니다. dplyr의 left_join은 왼쪽 표를 다 남기고 오른쪽에서 짝을 찾아 붙이고, inner_join은 양쪽에 다 있는 행만 남깁니다. by로 기준 열을 지정합니다.

long <- expr |>
  pivot_longer(-gene, names_to = "sample", values_to = "expr")

long |> left_join(meta, by = "sample") |> head(4)   # sample 기준으로 group·batch 붙이기

# A tibble: 4 × 5
  gene  sample  expr group   batch
  <chr> <chr>  <dbl> <chr>   <chr>
1 Gene1 S1      6.75 control b1
2 Gene1 S2      9.84 control b2
3 Gene1 S3      9.62 control b1
4 Gene1 S4     11.4  treated b2

기준 열의 이름이 양쪽에서 다르면 by = c("sample" = "sample_id")처럼 짝을 명시합니다. join 뒤에는 행 수가 예상대로인지 꼭 확인하세요 — 오른쪽 표에 같은 키가 여러 번 있으면 행이 의도치 않게 불어날 수 있습니다.

6. pivot_longer / pivot_wider — wide ↔ long

tidyr의 핵심은 표의 형식을 바꾸는 두 함수입니다. pivot_longer는 여러 열에 펼쳐진 값을 한 열로 눕혀 long (긴 형식)으로 만들고, pivot_wider는 거꾸로 한 열을 여러 열로 펼쳐 wide (넓은 형식)로 되돌립니다.

# wide → long: 샘플 열⁠(S1~S6)을 sample·expr 두 열로 눕히기
long <- expr |>
  pivot_longer(
    cols      = -gene,          # gene만 빼고 나머지 열을 모은다
    names_to  = "sample",       # 옛 열 이름이 들어갈 새 열
    values_to = "expr")         # 옛 셀 값이 들어갈 새 열
head(long, 3)

# A tibble: 3 × 3
  gene  sample  expr
  <chr> <chr>  <dbl>
1 Gene1 S1      6.75
2 Gene1 S2      9.84
3 Gene1 S3      9.62

이제 36행짜리 long이 됐고, ggplot2·통계 함수에 바로 넘길 수 있는 tidy 형태입니다. 반대로 long을 다시 발현행렬 모양으로 펴려면 pivot_wider를 씁니다.

# long → wide: 원래 발현행렬 모양으로 복원
long |>
  pivot_wider(names_from = sample, values_from = expr)   # sample 값이 열 이름이 된다

names_from은 어떤 열의 값을 새 열 이름으로 쓸지, values_from은 어떤 열의 값을 셀에 채울지를 가리킵니다. pivot_longer의 names_to/values_to와 방향이 반대라고 기억하면 헷갈리지 않습니다.

문자열이 한 열에 뭉쳐 있을 때는 separate로 쪼개고 unite로 합칩니다.

# "Gene1_control" 같은 합성 ID를 두 열로 분리
deg |>
  mutate(label = paste(gene, c("up","dn","ns","up","ns","up"), sep = "_")) |>
  separate(label, into = c("g", "dir"), sep = "_")   # _ 기준으로 g·dir 분리

7. 결측값 — drop_na / replace_na

실제 데이터에는 NA (결측값)가 섞이기 마련입니다. tidyr의 drop_na는 NA가 있는 행을 버리고, replace_na는 NA를 특정 값으로 채웁니다.

# 일부러 결측을 넣은 표
long_na <- long
long_na$expr[c(2, 5)] <- NA

long_na |> drop_na(expr)                          # expr이 NA인 행 제거
long_na |> mutate(expr = replace_na(expr, 0))     # NA를 0으로 채우기

drop_na()를 열 지정 없이 부르면 아무 열에든 NA가 있는 행을 모두 버리니, 보통은 drop_na(expr)처럼 대상 열을 명시하는 게 안전합니다.

8. dplyr · tidyr 치트시트

자주 쓰는 함수를 한 표로 정리했습니다.

9. 자주 발생하는 에러 & 해결

• filter에서 =를 썼다 — filter(group = "treated")는 에러입니다. 비교는 반드시 ==를 씁니다 (filter(group == "treated")). =는 인자 대입, ==는 같은지 비교라는 점을 늘 구분하세요.
• filter가 NA 행을 말없이 버린다 — filter(adj.P.Val < 0.05)는 adj.P.Val이 NA인 행을 조건 불만족으로 보고 제외합니다. NA도 남기려면 filter(adj.P.Val < 0.05 | is.na(adj.P.Val))처럼 명시합니다.
• group_by 후 그룹을 안 풀었다 — summarise 뒤에도 그룹이 남아 있으면 이어지는 mutate가 그룹 단위로 계산돼 엉뚱한 결과가 납니다. .groups = "drop"이나 ungroup()으로 풀어 주세요.
• pivot에서 names_to/values_to를 빠뜨렸다 — pivot_longer는 새로 생길 두 열의 이름을 지정해야 합니다. 빠뜨리면 기본 이름 (name·value)이 붙어 뒤 코드와 어긋납니다.
• join 후 행이 불어났다 — 오른쪽 표의 키가 유일하지 않으면 (한 sample이 여러 번 등장) 행이 곱으로 늘어납니다. join 전에 meta |> distinct(sample) 또는 count(sample)로 키 중복을 확인하세요.

References

1. Wickham, H. (2014). Tidy Data. Journal of Statistical Software, 59(10), 1–23. (tidy data 3원칙)
2. Wickham, H., Çetinkaya-Rundel, M., Grolemund, G. (2023). R for Data Science (2nd ed.). O'Reilly. (dplyr·tidyr·파이프 표준 교재, r4ds.hadley.nz)
3. dplyr 1.1+ documentation (tidyverse). filter·select·mutate·arrange·summarise·group_by·across·join·.by.
4. tidyr 1.3+ documentation (tidyverse). pivot_longer·pivot_wider·separate·unite·drop_na·replace_na.

Pipette & Pipeline · A bio portfolio journal

TL;DR — salmon·kallisto는 빠른 pseudo/selective-alignment로 transcript 수준의 정량⁠(quantification)을 냅니다⁠(quant.sf / abundance.h5). 그런데 DESeq2·edgeR·limma 같은 차등발현⁠(differential expression, DE) 도구는 보통 gene 수준 카운트 행렬에서 출발하죠. 그 둘 사이를 잇는 게 tximport입니다 — transcript→gene 집계와 평균 transcript 길이 보정을 한 번에 처리하고, DESeqDataSetFromTximport() 한 줄로 DE 입력까지 연결합니다.

이 글은 salmon 정량 폴더에서 시작해 tx2gene 매핑 → tximport() → 반환 구조 확인 → DESeq2 입력까지, 복붙하면 도는 코드로 따라갑니다. kallisto·edgeR·limma 경로와 countsFromAbundance 선택 기준, 자주 터지는 에러도 정리합니다. (코드는 tximport·DESeq2 공식 vignette 기준입니다.)

🔗 관련 글  ·  카운트 행렬이 준비된 다음 단계: DESeq2로 차등발현 분석  ·  edgeR 차등발현  ·  GEOquery로 공개데이터 받기

이 글에서 할 것

이 블로그의 DE 글들⁠(DESeq2·edgeR·limma-voom)은 전부 카운트 행렬이 이미 있다는 데서 출발합니다. 그런데 그 카운트는 어디서 오나요? 요즘 RNA-seq 정량은 salmon·kallisto로 transcript 수준 추정치를 먼저 뽑고, 그걸 gene 수준으로 모으는 흐름이 표준입니다. 이 글은 바로 그 'DESeq2 글의 한 단계 앞'을 채웁니다.

데모는 단순합니다 — salmon이 샘플마다 만들어 둔 quant.sf 폴더가 있다고 가정하고, tximport로 읽어 gene 수준 카운트로 모은 뒤, DESeqDataSetFromTximport()에 그대로 넣습니다. 최종 결과물은 DESeq2가 바로 받을 수 있는 dds 객체입니다.

txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ condition)
# 이 한 흐름이 transcript 정량을 gene 수준 DE 입력으로 바꾼다

0. 왜 tximport인가 — 1분 직관

salmon (selective alignment)·kallisto (pseudo-alignment)는 리드를 게놈에 일일이 정렬하지 않고, transcriptome을 기준으로 어느 transcript에서 나왔을 확률이 높은지를 빠르게 추정합니다. 출력은 transcript마다의 추정 카운트와 TPM (transcripts per million)이죠. 문제는 두 가지입니다.

• 수준이 다르다. salmon·kallisto는 transcript 수준인데, DE 분석은 보통 gene 수준에서 합니다. 한 유전자의 여러 isoform을 합쳐야 합니다.
• 길이 보정이 필요하다. 같은 유전자라도 샘플마다 우세하게 발현된 isoform이 다르면 평균 transcript 길이가 달라집니다. 단순히 카운트만 더하면 이 차이가 묻혀 편향이 생깁니다.

tximport는 이 둘을 함께 처리합니다 — transcript 추정치를 gene 단위로 합산하고, 샘플별 평균 transcript 길이를 offset으로 넘겨 DESeq2·edgeR이 정규화에 반영하도록 합니다⁠(Soneson et al. 2015). 덕분에 단순 합산보다 정확한 gene 수준 추론이 됩니다.

1. 준비물 — 설치

Bioconductor 패키지 tximport가 핵심이고, 연결할 DE 패키지로 DESeq2⁠(edgeR·limma도 동일하게 가능)를 깝니다. 메타데이터를 자동으로 붙이고 싶으면 tximeta도 유용합니다. kallisto의 .h5를 읽으려면 rhdf5가 필요합니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("tximport", "DESeq2", "tximeta", "rhdf5"))   # 최초 1회

library(tximport)   # 정량 불러오기 + transcript→gene 집계
library(DESeq2)     # DESeqDataSetFromTximport로 바로 연결

여기서 salmon·kallisto 출력이 무엇인지 먼저 짚어 둡시다. salmon은 샘플 폴더마다 quant.sf⁠(transcript별 추정 카운트·TPM·유효 길이)를, kallisto는 abundance.h5·abundance.tsv·run_info.json을 남깁니다. tximport는 이 파일들을 가리키는 경로 벡터만 받으면 됩니다.

2. 파일 경로 벡터 만들기

tximport의 첫 입력은 각 샘플의 정량 파일 경로를 담은 벡터입니다. 핵심은 벡터에 이름⁠(names)을 붙이는 것 — 그 이름이 곧 카운트 행렬의 열 이름이 됩니다.

# 샘플 메타데이터⁠(보통 별도 .csv/.tsv로 관리)
dir     <- "/path/to/project"
coldata <- read.csv(file.path(dir, "samples.csv"))   # run, condition 등 컬럼
coldata$condition <- factor(coldata$condition, levels = c("control", "treated"))

# salmon 출력 구조: <dir>/salmon/<run>/quant.sf
files <- file.path(dir, "salmon", coldata$run, "quant.sf")
names(files) <- coldata$run     # 열 이름이 될 샘플 ID

all(file.exists(files))         # 모두 TRUE여야 진행

[1] TRUE

경로가 하나라도 FALSE면 거기서 멈추세요 — 폴더 구조나 압축 여부⁠(quant.sf vs quant.sf.gz)부터 확인하는 게 빠릅니다. salmon을 --gzip으로 돌렸다면 파일명이 quant.sf.gz입니다.

3. tx2gene — transcript ID를 gene ID로 잇는 지도

tximport가 transcript를 gene으로 모으려면 어떤 transcript가 어떤 gene에 속하는지 알려 주는 2열 표가 필요합니다 — 1열 transcript ID, 2열 gene ID. 이게 tx2gene입니다. GENCODE·Ensembl GTF에서 만들거나, TxDb·EnsDb 같은 주석 패키지에서 뽑습니다.

# 방법 A — TxDb 주석 패키지에서 (예: 사람 hg38)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v44.annotation.gtf")   # 또는 기존 TxDb 사용
k    <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")          # 1열 TXNAME, 2열 GENEID

head(tx2gene, 3)

              TXNAME            GENEID
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000488147.1 ENSG00000227232.5

# 방법 B — 이미 만들어 둔 매핑 표를 읽기 (GENCODE 동봉 파일 등)
tx2gene <- read.csv(file.path(dir, "tx2gene.gencode.v44.csv"))

4. tximport 실행 — salmon

이제 본론입니다. type = "salmon"으로 지정하고 tx2gene을 넘기면, tximport가 모든 샘플의 quant.sf를 읽어 gene 수준으로 집계합니다.

txi <- tximport(files,
                type      = "salmon",
                tx2gene   = tx2gene,
                ignoreTxVersion = TRUE)   # ID 버전 불일치 방지⁠(필요 시)

names(txi)

[1] "abundance"           "counts"              "length"
[4] "countsFromAbundance"

tximport 한 줄이 진행 상황을 찍으며⁠(reading in files 1 2 3 ...) transcript 추정치를 gene 단위로 모읍니다. 반환값은 리스트입니다.

kallisto일 때는?

도구만 바꾸면 됩니다. kallisto는 abundance.h5⁠(권장, rhdf5 필요) 또는 abundance.tsv를 읽습니다.

# kallisto 출력: <dir>/kallisto/<run>/abundance.h5
files_k <- file.path(dir, "kallisto", coldata$run, "abundance.h5")
names(files_k) <- coldata$run

txi_k <- tximport(files_k, type = "kallisto", tx2gene = tx2gene)
# .tsv를 쓰면 GENCODE의 'ENST|ENSG|...' 헤더 분리를 위해 ignoreAfterBar = TRUE

salmon·kallisto·RSEM 모두 같은 인터페이스로 읽힙니다 — type만 바뀌고 나머지는 동일합니다.

5. 반환 구조 — counts · abundance · length

tximport가 돌려주는 리스트의 세 핵심 원소를 이해하면, 뒤 단계가 전부 쉬워집니다.

dim(txi$counts)              # gene × sample 카운트 행렬
txi$counts[1:3, 1:2]         # DE 도구가 받을 추정 카운트
txi$length[1:3, 1:2]         # gene별 평균 transcript 길이⁠(샘플마다 다름)

[1] 38194     6
                  sample1 sample2
ENSG00000000003   742.30  689.11
ENSG00000000419   531.07  604.88
ENSG00000000457   198.44  221.50

• $counts — gene × sample 추정 카운트 행렬. DESeq2·edgeR·limma가 받는 그 행렬입니다. salmon의 확률적 할당 때문에 정수가 아닌 실수일 수 있습니다⁠(정상).
• $abundance — 같은 모양의 TPM 행렬. 라이브러리 크기·길이가 보정된 발현량이라, 빠른 탐색·시각화나 샘플 간 단순 비교에 씁니다⁠(DE 검정 입력으로 직접 쓰진 않습니다).
• $length — gene별 평균 transcript 길이 행렬. 샘플마다 우세 isoform이 달라 값이 다릅니다. 이 행렬이 길이 offset으로 DE 도구에 넘어가, isoform 전환에 따른 편향을 보정합니다.

6. DESeqDataSetFromTximport — DESeq2로 연결

가장 깔끔한 길입니다. txi 리스트를 통째로 넘기면, DESeq2가 $counts를 카운트로, $length를 길이 offset으로 자동 등록합니다. 별도 정규화 코드가 필요 없습니다.

# coldata의 행 순서가 txi$counts의 열 순서와 같아야 함
rownames(coldata) <- coldata$run
stopifnot(all(rownames(coldata) == colnames(txi$counts)))

dds <- DESeqDataSetFromTximport(txi,
                                colData = coldata,
                                design  = ~ condition)
dds <- DESeq(dds)
res <- results(dds)                  # 이후는 일반 DESeq2 워크플로와 동일

여기서부터는 DESeq2 글과 완전히 같습니다 — results()로 logFC·padj를 뽑고, 시각화하면 됩니다. tximport는 그 출발점인 dds를 만들어 주는 역할까지입니다.

edgeR·limma로 연결하려면

edgeR·limma도 길이 offset만 제대로 넘기면 됩니다. 최신 edgeR⁠(≥ 4.10)은 전용 함수가 있고, 그 이전 버전이라면 countsFromAbundance로 길이 보정된 카운트를 만들어 DGEList에 넣는 방식이 이식성이 좋습니다.

# (a) 최신 edgeR — 길이 offset을 알아서 등록
library(edgeR)
y <- DGEListFromTximport(txi)        # edgeR >= 4.10
y <- normLibSizes(y)

# (b) 이식성 있는 방식 — 길이 보정 카운트로 DGEList 구성
txi_ls <- tximport(files, type = "salmon", tx2gene = tx2gene,
                   countsFromAbundance = "lengthScaledTPM")
y <- DGEList(txi_ls$counts, group = coldata$condition)
y <- normLibSizes(y)                 # 이후 voom 또는 glmQLFit

lengthScaledTPM으로 만든 카운트는 길이·라이브러리 크기가 이미 반영돼 있어, 별도 offset 없이 edgeR·limma-voom에 바로 넣을 수 있습니다.

7. countsFromAbundance — 언제 무엇을 쓰나

countsFromAbundance는 "카운트를 TPM⁠(abundance)에서 다시 만들지" 정하는 옵션입니다. DE 도구마다 권장값이 다릅니다.

규칙은 단순합니다 — DESeq2로 갈 거면 "no" (기본)로 두고 offset을 자동으로 쓰게 하고, edgeR·limma에 카운트만 직접 넣을 거면 "lengthScaledTPM"으로 길이 보정을 카운트에 녹이세요. transcript 수준 차등사용⁠(DTU)을 볼 때만 "dtuScaledTPM" + txOut = TRUE입니다.

8. salmon vs 정렬 기반⁠(featureCounts) — 무엇이 다른가

카운트를 만드는 길은 크게 둘입니다. salmon·kallisto 같은 pseudo/selective-alignment 방식과, STAR·HISAT2로 게놈에 정렬한 뒤 featureCounts·htseq-count로 세는 정렬 기반 방식입니다. tximport는 전자의 출력을 받습니다.

둘 다 유효하지만, 요즘 표준 RNA-seq 정량은 속도와 isoform 인지 때문에 salmon·kallisto로 기우는 추세입니다. 그리고 그 출력을 DE로 잇는 다리가 정확히 tximport입니다. 게놈 좌표가 필요한 분석⁠(신규 transcript 발굴, 변이 연계 등)에서는 정렬 기반이 여전히 제 몫을 합니다.

9. 자주 발생하는 에러 & 해결

• tximport 후 카운트가 거의 0이거나 gene이 몇 개뿐 — tx2gene의 transcript ID와 salmon 인덱스 ID의 버전 불일치가 1순위입니다. head(rownames(read.delim(files[1])))로 실제 ID를 보고, ignoreTxVersion = TRUE를 주거나 tx2gene ID를 맞추세요.
• Error: all(file.exists(files)) is not TRUE — 경로가 틀렸거나 압축⁠(quant.sf.gz)입니다. files를 직접 출력해 폴더 구조와 파일명을 눈으로 확인하세요.
• tx2gene 관련 경고로 transcript가 누락 — GTF·TxDb 버전이 salmon 인덱스를 만든 주석과 다른 경우입니다. 정량과 같은 버전의 주석으로 tx2gene을 만드는 게 원칙입니다.
• kallisto .h5를 못 읽음 — rhdf5 미설치입니다. BiocManager::install("rhdf5") 후 다시 시도하거나, abundance.tsv로 읽으세요⁠(GENCODE면 ignoreAfterBar = TRUE).
• DESeqDataSetFromTximport에서 행/열 불일치 — colData의 행 순서가 txi$counts의 열 순서와 달라서입니다. rownames(coldata)와 colnames(txi$counts)가 같은지 stopifnot으로 먼저 확인하세요.
• 카운트가 정수가 아니라고 경고 — salmon의 확률적 할당 때문에 실수 카운트는 정상입니다. DESeqDataSetFromTximport는 이를 알아서 처리합니다⁠(직접 반올림하지 마세요).

References

1. Soneson, C., Love, M. I., Robinson, M. D. (2015). Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research, 4, 1521. (tximport 방법론·길이 offset)
2. Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods, 14(4), 417–419.
3. Bray, N. L., Pimentel, H., Melsted, P., Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification (kallisto). Nature Biotechnology, 34(5), 525–527.
4. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. (DESeqDataSetFromTximport)
5. Love, M. I., Soneson, C., et al. (2020). Tximeta: Reference sequence checksums for provenance identification in RNA-seq. PLoS Computational Biology, 16(2), e1007664.
6. tximport vignette (Bioconductor). tximport·tx2gene·countsFromAbundance·ignoreTxVersion·DGEListFromTximport.

Pipette & Pipeline · A bio portfolio journal

TL;DR — ggplot2 (그래픽 문법, grammar of graphics)는 ggplot(data, aes(x, y)) + geom_*()처럼 레이어를 +로 쌓아 그림을 그립니다. 데이터·미적매핑⁠(aesthetic mapping)·기하 도형⁠(geom)을 분리해서 생각하는 게 전부예요.

이 글은 발현·그룹 더미 데이터로 산점도·막대⁠(평균±오차)·박스플롯을 그리고, 색·소분할⁠(facet)·테마·저장까지 복붙하면 그대로 도는 코드로 따라갑니다. 이 블로그의 Volcano plot·히트맵·PCA 글이 쓰는 문법의 바닥을 채우는 글입니다. (코드는 ggplot2 4.0 / 3.5+ 기준으로 실제 실행해 확인했습니다.)

🔗 관련 글  ·  이 문법을 응용한 바이오 그래프: Volcano plot 시각화⁠(ggplot2 응용)  ·  ComplexHeatmap 히트맵  ·  PCA로 샘플 QC

이 글에서 그릴 것

최종 결과물은 산점도 + 막대그래프 + 박스플롯, 그리고 유전자별로 쪼갠 facet 패널입니다. 발현값⁠(log2 발현)과 그룹⁠(대조군·처리군) 더미 데이터를 코드로 만들어, ggplot2 문법을 한 겹씩 쌓아 보겠습니다. 핵심은 단 한 줄로 요약됩니다 — 데이터를 넣고⁠(ggplot), 무엇을 어디에 그릴지 매핑하고⁠(aes), 어떤 도형으로 그릴지 고른다⁠(geom_*). 나머지는 전부 이 위에 +로 얹는 장식입니다.

# ggplot2의 한 문장 — 이 패턴만 외우면 절반은 끝
ggplot(데이터, aes(x = 변수1, y = 변수2)) +   # 데이터 + 미적매핑
  geom_point()                                # 기하 도형⁠(여기선 점)

1. 사전 준비 — 설치와 예시 데이터

ggplot2는 단독 패키지로 써도 되고, 데이터 정리용 패키지 묶음인 tidyverse의 일부로 설치해도 됩니다. tidyverse를 깔면 ggplot2 (그래픽 문법)·dplyr (데이터 가공)·tidyr 등이 한 번에 들어와 편합니다.

install.packages("tidyverse")   # ggplot2 + dplyr + tidyr ... 한 번에 (최초 1회)
# 또는 그래프만 필요하면:
# install.packages("ggplot2")

library(ggplot2)   # 그래프 문법

예시는 바이오에서 가장 흔한 형태 — 샘플마다 유전자 발현값이 있고, 대조군·처리군 그룹이 붙은 긴 형식⁠(long format) 데이터프레임으로 만듭니다. set.seed()로 난수를 고정했으니 여러분 손에서도 똑같은 숫자가 나옵니다.

set.seed(42)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")

expr <- data.frame(
  sample   = rep(paste0("S", 1:12), each = 4),                       # 샘플 12개
  group    = rep(rep(c("Control", "Treated"), each = 4), times = 3), # 대조/처리
  gene     = rep(genes, times = 12),                                 # 유전자 4종
  log2expr = round(rnorm(48, mean = 8, sd = 1.5), 2),                # log2 발현값
  depth    = round(runif(48, 5, 30), 1)                              # 시퀀싱 깊이⁠(가짜)
)
expr$group <- factor(expr$group, levels = c("Control", "Treated"))   # 대조군을 기준으로

# 처리군에서 EGFR만 발현을 올려, '효과'를 일부러 심는다
expr$log2expr[expr$gene == "EGFR" & expr$group == "Treated"] <-
  expr$log2expr[expr$gene == "EGFR" & expr$group == "Treated"] + 2.5

head(expr, 4)

  sample   group  gene log2expr depth
1     S1 Control BRCA1    10.06  13.3
2     S1 Control  TP53     7.15  17.9
3     S1 Control  EGFR     8.54  23.6
4     S1 Control   MYC     8.95  20.5

무엇이 보이나 — 한 행이 한 관측⁠(샘플 × 유전자)입니다. 열은 sample·group·gene·log2expr·depth. ggplot2는 이렇게 한 행 = 한 점인 긴 형식을 가장 좋아합니다. 엑셀처럼 유전자가 열마다 흩어진 넓은 형식⁠(wide)이면 tidyr::pivot_longer()로 먼저 녹여야 합니다.

2. 첫 그래프 — ggplot() + aes() + geom_point()

가장 단순한 산점도부터. 깊이⁠(depth)와 발현⁠(log2expr)의 관계를 점으로 찍습니다. 세 조각만 기억하세요 — ggplot(데이터), aes(x, y), geom_point().

ggplot(expr, aes(x = depth, y = log2expr)) +   # 데이터 + x·y 매핑
  geom_point()                                  # 점으로 그린다

무엇이 보이나 — 빈 좌표계 위에 점 48개가 찍힙니다. aes() 안에 적은 변수는 데이터 열 이름이고, ggplot2가 알아서 축 범위·눈금·라벨을 잡아 줍니다. +로 다음 줄에 geom_*을 얹는 게 핵심이에요. +는 반드시 줄 끝에 두세요 — 줄 맨 앞에 두면 R은 앞 줄에서 문장이 끝난 줄 알고 에러를 냅니다.

3. 색·모양 매핑 — aes 안 vs 밖

이제 그룹을 색으로 구분해 봅시다. 그룹별로 색을 다르게 하려면 color를 aes() 안에 넣습니다. 데이터 열에 따라 변하는 건 전부 aes() 안, 모든 점에 똑같이 적용할 고정값은 aes() 밖입니다.

ggplot(expr, aes(x = depth, y = log2expr,
                 color = group, shape = group)) +  # 그룹 → 색·모양 (매핑)
  geom_point(size = 3, alpha = 0.8)                # 크기·투명도는 고정값 (밖)

무엇이 보이나 — 대조군·처리군이 다른 색·다른 모양으로 갈리고, 오른쪽에 범례가 자동으로 생깁니다. 반면 size = 3·alpha = 0.8은 aes() 밖에 있으니 모든 점에 똑같이 적용되고 범례에 안 뜹니다. 이 안/밖 구분이 ggplot2 초보가 가장 많이 헷갈리는 지점입니다.

4. 막대그래프 — 평균±오차 (geom_col / stat_summary)

막대는 두 갈래입니다. 값을 그대로 막대 높이로 쓰면 geom_col(), 행 개수를 세서 빈도 막대를 그리면 geom_bar()⁠(기본 stat = "count"). 바이오에서 자주 쓰는 그룹 평균 ± 표준오차 막대는 stat_summary() 한 줄로 깔끔하게 나옵니다 — 평균을 따로 계산할 필요가 없어요.

ggplot(expr, aes(x = gene, y = log2expr, fill = group)) +
  stat_summary(fun = mean, geom = "col",                 # 그룹 평균을 막대로
               position = position_dodge(0.8), width = 0.7) +
  stat_summary(fun.data = mean_se, geom = "errorbar",     # 평균±표준오차 막대
               position = position_dodge(0.8), width = 0.2) +
  labs(x = NULL, y = "평균 log2 발현", fill = "그룹")

무엇이 보이나 — 유전자 4종 × 그룹 2개의 평균 막대가 나란히⁠(position_dodge) 서고, 위에 오차막대⁠(errorbar)가 얹힙니다. EGFR만 처리군에서 막대가 확 올라간 게 보이죠 — 우리가 일부러 심은 효과입니다. fill은 막대·박스처럼 면을 칠하는 미적매핑이고, color는 선·점의 색입니다. 둘을 혼동하면 막대 테두리만 칠해지니 주의하세요.

5. 박스플롯 — 분포를 한눈에 (geom_boxplot)

막대는 평균만 보여 주지만, 박스플롯⁠(box plot)은 중앙값·사분위수·이상치까지 분포를 통째로 보여 줍니다. 발현 데이터엔 박스플롯이 더 정직한 선택일 때가 많아요. 점을 흩뿌리는 geom_jitter()를 겹치면 실제 데이터 개수까지 드러납니다.

ggplot(expr, aes(x = group, y = log2expr, fill = group)) +
  geom_boxplot(width = 0.6, alpha = 0.8) +
  geom_jitter(width = 0.12, size = 1.6, alpha = 0.5)   # 점을 살짝 흔들어 겹침 방지

무엇이 보이나 — 그룹별 박스⁠(상자·수염)와 그 위에 흩뿌린 점들이 함께 보입니다. 레이어는 쓴 순서대로 위에 쌓이므로, geom_boxplot을 먼저, geom_jitter를 나중에 둬야 점이 박스 위에 올라옵니다. 순서를 바꾸면 박스가 점을 덮어 버립니다.

6. 소분할⁠(facet) — 유전자별로 쪼개기

지금 박스플롯은 유전자를 뭉뚱그렸습니다. 유전자마다 따로 보고 싶을 때 facet_wrap()이 데이터를 자동으로 쪼개 작은 패널 여러 개로 깝니다 — 코드 한 줄로요.

ggplot(expr, aes(x = group, y = log2expr, fill = group)) +
  geom_boxplot(width = 0.6) +
  facet_wrap(~ gene, nrow = 1)        # gene 값마다 패널 하나씩, 한 줄로

무엇이 보이나 — BRCA1·TP53·EGFR·MYC 네 패널이 한 줄로 나란히 뜨고, EGFR 패널에서만 처리군 박스가 위로 올라간 게 또렷합니다. 변수 둘로 격자를 짜려면 facet_grid(행변수 ~ 열변수)를 씁니다. facet_wrap(~ gene)처럼 물결표⁠(~) 뒤에 쪼갤 변수를 적는다고 기억하세요.

7. 색·축·라벨·테마 — 외형 다듬기

마지막으로 발표용으로 다듬습니다. 색은 scale_fill_manual()로 직접 지정, 제목·축 라벨은 labs(), 전체 외형은 theme_*() + 세부 theme()로 손봅니다. 이 레이어들은 그림의 내용이 아니라 표현을 바꿉니다.

ggplot(expr, aes(x = group, y = log2expr, fill = group)) +
  geom_boxplot(width = 0.6) +
  facet_wrap(~ gene, nrow = 1) +
  scale_fill_manual(values = c(Control = "#6BAB97",        # 색 직접 지정
                               Treated = "#E76F51")) +
  labs(title = "유전자별 발현 비교",                        # 제목·축·범례 라벨
       x = NULL, y = expression(log[2]~expression),
       fill = "그룹") +
  theme_bw(base_size = 13) +                                # 흑백 테마, 글자 크기
  theme(legend.position = "top",                            # 범례를 위로
        strip.background = element_rect(fill = "grey95"))   # facet 라벨 배경

무엇이 보이나 — 손맛 나는 색, 위로 올라간 범례, 깔끔한 흰 배경의 발표용 그림이 됩니다. theme_bw()·theme_minimal()·theme_classic() 같은 완성형 테마를 먼저 고른 뒤, 세부는 theme()로 덮어쓰는 게 정석입니다. 축을 로그로 바꾸려면 scale_y_log10()을, 시간에 따른 추세는 geom_line(linewidth = 1)을 얹으면 됩니다. (ggplot2 3.4부터 선 굵기는 size가 아니라 linewidth입니다.)

8. 저장 — ggsave()

화면이 아니라 파일로 저장할 땐 ggsave() 한 줄. 마지막에 그린 그림을 잡아 저장하며, 크기⁠(인치)·해상도⁠(dpi)를 직접 지정합니다. 논문·발표용은 dpi = 300 이상을 권장합니다.

p <- ggplot(expr, aes(group, log2expr, fill = group)) +
  geom_boxplot() + facet_wrap(~ gene, nrow = 1)

ggsave("expr_boxplot.png", plot = p,
       width = 9, height = 4, dpi = 300)   # 9×4인치, 300dpi PNG
# 확장자만 바꾸면 형식이 바뀜: .pdf .svg .tiff ...

무엇이 보이나 — 작업 폴더에 expr_boxplot.png가 떨어집니다. 확장자만 .pdf·.svg로 바꾸면 벡터 그래픽으로 저장돼 무한 확대해도 안 깨집니다 — 논문 그림은 벡터가 정석입니다. plot =을 생략하면 가장 최근 그림이 저장됩니다.

geom 치트시트

자주 쓰는 기하 도형⁠(geom)과 통계 변환을 한 표로 정리합니다.

자주 발생하는 에러 & 해결

• color가 안 먹히고 막대 테두리만 칠해짐 — 면을 칠하는 건 fill, 선·점을 칠하는 건 color입니다. 막대·박스는 aes(fill = group)을 쓰세요.
• Error: unexpected '+' / 그림이 한 겹만 나옴 — +를 줄 맨 앞에 뒀기 때문입니다. +는 항상 앞 줄 끝에 붙여, R이 "다음 줄로 이어진다"고 알게 하세요.
• 그룹별로 색이 안 갈리고 전부 한 색 — color를 aes() 밖에 적었습니다. 데이터에 따라 변해야 하면 aes(color = group)처럼 안으로 옮기세요.
• 막대·축 순서가 알파벳순으로 멋대로 — ggplot2는 문자형을 알파벳순으로 정렬합니다. 원하는 순서는 factor(x, levels = c(...))로 요인 수준을 먼저 지정하세요⁠(이 글의 Control → Treated처럼).
• 색이 의도와 다르게 그라데이션으로 나옴 — 색 매핑 변수가 연속형⁠(숫자)이면 그라데이션, 범주형⁠(요인·문자)이면 구분색입니다. 범주로 쓰려면 factor()로 감싸세요.
• Error ... object 'depth' not found — aes() 안의 이름은 데이터프레임의 열 이름이어야 합니다. 철자·대소문자를 names(expr)로 확인하세요.

References

1. Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis (2nd ed.). Springer. (그래픽 문법의 표준 교재)
2. Wickham, H., et al. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686.
3. ggplot2 공식 문서 (tidyverse.org/ggplot2) — aes·geom_*·facet_wrap·scale_*·theme·ggsave 레퍼런스.
4. Wilkinson, L. (2005). The Grammar of Graphics (2nd ed.). Springer. (ggplot2가 따르는 그래픽 문법의 원전)
5. Healy, K. (2018). Data Visualization: A Practical Introduction. Princeton University Press. (ggplot2 실전 시각화 안내)

Pipette & Pipeline · A bio portfolio journal

TL;DR — RNA-seq 차등발현⁠(differentially expressed gene, DEG) 분석의 'DE 3대장' 중 마지막이 limma-voom입니다. DESeq2·edgeR이 카운트를 음이항⁠(negative binomial, NB) 모델로 직접 다룬다면, limma-voom은 voom (variance modeling at the observation level)으로 카운트를 로그-CPM 연속값 + 가중치⁠(precision weight)로 바꾼 뒤, 마이크로어레이용으로 만들어진 limma의 선형모형⁠(linear model) + 경험적 베이즈⁠(empirical Bayes)를 그대로 적용합니다.

이 글은 edgeR 글과 같은 pasilla 데이터⁠(초파리, 7 샘플)로 복붙하면 그대로 도는 limma-voom 코드를 따라갑니다 — DGEList → filterByExpr → calcNormFactors → voom → lmFit → eBayes → topTable. 마지막에 DESeq2 vs edgeR vs limma-voom을 모델·검정·강점·선택 기준으로 정리합니다. (코드는 limma/edgeR 공식 vignette 기준이며, 실제 실행해 출력을 옮겼습니다.)

🔗 관련 글  ·  같이 보면 좋은 RNA-seq 분석: DESeq2로 차등발현 분석  ·  edgeR 차등발현⁠(DESeq2 비교)  ·  p값과 FDR⁠(다중검정보정)

이 글에서 만들 것

최종 결과물은 voom 평균-분산 plot + DEG 표 + volcano plot입니다. 카운트 행렬을 넣으면 voom이 RNA-seq 특유의 평균-분산 관계를 한 장의 그림으로 보여 주고, 유전자마다 logFC·AveExpr·t·P.Value·adj.P.Val가 한 표로 나옵니다. 아래 미리보기처럼 topTable() 한 줄이면 상위 DEG가 정렬돼 나옵니다.

              logFC AveExpr     t  P.Value  adj.P.Val    B
FBgn0025111    2.92    6.25  26.8  2.5e-10   1.0e-06  14.2
FBgn0003360   -3.13    7.78 -24.9  5.0e-10   1.1e-06  13.6
FBgn0039155   -4.62    4.68 -27.9  1.7e-10   1.0e-06  13.3

0. 왜 limma-voom인가 — 1분 직관

RNA-seq 카운트는 정수이고 평균이 크면 분산도 커집니다⁠(과대산포, overdispersion). DESeq2·edgeR은 이 성질을 음이항⁠(NB) 분포로 정면 돌파합니다. limma-voom의 발상은 다릅니다 — "그 평균-분산 관계를 추정해서 가중치로 바꿔 버리면, 잘 만들어진 정규분포용 도구⁠(limma)를 그대로 쓸 수 있지 않을까?"

• limma는 원래 마이크로어레이⁠(연속 형광값)용으로 설계된, 선형모형 + 경험적 베이즈의 성숙한 프레임워크입니다. 복잡한 설계·블로킹·대비⁠(contrast)를 다루는 기능이 풍부합니다.
• 문제는 카운트입니다. 저발현 유전자는 분산이 크고 고발현 유전자는 작은데, 이걸 무시하고 선형모형을 돌리면 가정이 깨집니다.
• voom (variance modeling at the observation level)이 다리를 놓습니다 — 평균-분산 추세를 LOWESS로 적합해, 각 관측에 분산의 역수에 해당하는 가중치⁠(precision weight)를 부여합니다. 저발현·고분산 관측은 가볍게, 고발현·저분산 관측은 무겁게 보도록.

1. 준비물 — 설치

Bioconductor 패키지 limma와 edgeR을 씁니다. limma-voom은 카운트 입력·정규화에 edgeR의 DGEList·calcNormFactors를 그대로 쓰기 때문에 둘 다 필요합니다. 예제 데이터 pasilla, volcano용 ggplot2도 설치합니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("limma", "edgeR", "pasilla"))   # 최초 1회
install.packages("ggplot2")

library(edgeR)    # DGEList · filterByExpr · calcNormFactors 제공
library(limma)    # voom · lmFit · eBayes · topTable 제공
library(pasilla)
library(ggplot2)

2. 예제 데이터 — pasilla 카운트 불러오기

pasilla는 초파리 세포에서 pasilla 유전자를 RNAi로 억제한 RNA-seq로, 처리⁠(treated) 3개 + 비처리⁠(untreated) 4개 = 7 샘플입니다. edgeR 글과 똑같은 데이터를 쓰는 건 의도된 선택입니다 — 뒤에서 세 도구의 결과를 같은 데이터 위에서 비교하려는 것이죠.

# pasilla 패키지에 동봉된 유전자 카운트 행렬⁠(.tsv)을 읽는다
fn  <- system.file("extdata", "pasilla_gene_counts.tsv",
                   package = "pasilla", mustWork = TRUE)
cts <- as.matrix(read.csv(fn, sep = "\t", row.names = "gene_id"))

# 그룹⁠(조건) 벡터: 열 이름이 'treated'로 시작하면 처리군
group <- factor(ifelse(grepl("^treated", colnames(cts)), "treated", "untreated"),
                levels = c("untreated", "treated"))   # 비처리를 기준⁠(reference)으로
dim(cts); group

[1] 14599     7
[1] untreated untreated untreated untreated treated treated treated
Levels: untreated treated

levels의 첫 값이 기준⁠(reference)입니다. untreated를 앞에 두면 logFC가 "비처리 대비 처리" 방향이 됩니다 — 부호 해석이 직관적이죠.

3. DGEList — 카운트를 담는 그릇

limma-voom도 출발은 edgeR과 같습니다. 카운트와 그룹을 DGEList에 담습니다. voom이 이 객체를 입력으로 받아 로그-CPM과 가중치를 계산합니다.

dge <- DGEList(counts = cts, group = group)
dge$samples

               group lib.size norm.factors
untreated1 untreated 13972512            1
untreated2 untreated 21911438            1
...
treated3     treated 10343856            1

lib.size⁠(라이브러리 크기)가 자동 계산돼 들어가고, norm.factors는 아직 전부 1입니다 — TMM이 5단계에서 채웁니다. voom의 가중치는 이 라이브러리 크기·정규화 인자를 반영하므로 순서가 중요합니다.

4. filterByExpr — 저발현 유전자 거르기

거의 발현되지 않는 유전자는 평균-분산 추세를 망치고 검정력만 깎습니다. edgeR의 filterByExpr()로 "그룹 크기를 고려해 충분히 발현된" 유전자만 남깁니다⁠(수동 CPM 컷오프보다 권장). voom 이전에 거는 게 핵심입니다 — 저발현 유전자가 남으면 평균-분산 곡선이 왼쪽에서 치솟아 가중치 추정이 흔들립니다.

keep <- filterByExpr(dge, group = group)        # 그룹 설계를 반영한 자동 필터
dge  <- dge[keep, , keep.lib.sizes = FALSE]      # 필터 후 라이브러리 크기 재계산
table(keep)

keep
FALSE  TRUE
 6680  7919

keep.lib.sizes = FALSE가 중요합니다 — 유전자를 버렸으니 라이브러리 크기를 다시 합산해야 이후 CPM·가중치가 정확합니다. 14599개 중 7919개가 남았습니다.

5. calcNormFactors — TMM 정규화

샘플마다 시퀀싱 깊이가 다르고, 소수의 초고발현 유전자가 구성 편향⁠(composition bias)을 만듭니다. edgeR의 기본 정규화 TMM (trimmed mean of M-values)으로 샘플별 스케일 인자를 구합니다. voom은 이 norm.factors를 effective library size에 반영해 로그-CPM을 계산합니다.

dge <- calcNormFactors(dge)        # method = "TMM"가 기본값
round(dge$samples$norm.factors, 3)

[1] 1.025 1.006 0.966 0.948 1.084 0.984 0.993

TMM은 카운트를 바꾸지 않습니다. norm.factors라는 곱셈 보정값만 만듭니다. limma-voom·edgeR·DESeq2가 정규화 철학은 달라도⁠(TMM vs median-of-ratios) 출발점은 같습니다 — "라이브러리 크기 차이를 보정한 뒤 비교한다."

6. design — 실험 설계 행렬

검정에 앞서 설계 행렬⁠(design matrix)을 만듭니다. voom·lmFit 모두 이 행렬을 받습니다. 2그룹 비교는 ~ group 한 줄이면 됩니다. 배치⁠(batch)가 있으면 ~ batch + group처럼 보정 변수를 먼저 넣습니다.

design <- model.matrix(~ group)    # 1열=절편, 2열=grouptreated(우리가 볼 효과)
colnames(design)

[1] "(Intercept)"   "grouptreated"

관심 효과는 2번째 열 grouptreated입니다 — 뒤에서 coef = 2로 지목합니다. 셀 단위 평균이 필요하면 ~ 0 + group 후 makeContrasts()로 대비를 직접 짤 수도 있습니다⁠(9단계 박스).

7. voom — 평균-분산을 가중치로 (limma-voom의 심장)

이 단계가 핵심입니다. voom()은 두 가지를 동시에 합니다 — (1) 카운트를 로그2-CPM⁠(log2 counts per million) 연속값으로 변환하고, (2) 평균-분산 추세를 LOWESS로 적합해 각 관측에 가중치⁠(precision weight)를 부여합니다. plot = TRUE면 그 추세를 그림으로 보여 줍니다.

v <- voom(dge, design, plot = TRUE)   # 카운트 → 로그2-CPM + 가중치, 추세 plot
v                                      # EList 객체: $E(로그CPM), $weights(가중치)

An object of class "EList"
$E           : 7919 genes × 7 samples  (로그2-CPM 값)
$weights     : 7919 × 7                (관측별 precision weight)
$design      : (Intercept), grouptreated

v$E는 로그2-CPM 발현행렬, v$weights는 같은 크기의 가중치 행렬입니다. 저발현 유전자의 관측은 가중치가 작고⁠(예: 3~7), 고발현 유전자는 큽니다⁠(예: 40~50). lmFit이 이 가중치로 가중최소제곱⁠(weighted least squares)을 풉니다.

8. lmFit + eBayes — 가중 선형모형 + 경험적 베이즈

이제 limma 본연의 영역입니다. lmFit()이 유전자마다 가중 선형모형을 적합하고, eBayes()가 경험적 베이즈로 분산 추정을 안정화합니다 — 수만 개 유전자의 분산 정보를 빌려 와 각 유전자의 분산을 공통 추세 쪽으로 수축시키고, moderated t-통계량을 만듭니다.

fit <- lmFit(v, design)    # 가중치를 반영한 선형모형 적합
fit <- eBayes(fit)         # 경험적 베이즈 → moderated t-통계량

edgeR의 QL F-검정, DESeq2의 Wald 검정에 대응하는 것이 limma의 moderated t-검정입니다. 셋 다 "유전자별 분산을 그대로 믿지 말고, 전체 분포에서 빌려 와 안정화한다"는 경험적 베이즈 철학을 공유합니다 — 검정 통계량의 형태만 다릅니다.

9. topTable — 결과 표 뽑기

topTable()로 차등발현 결과를 정렬된 표로 받습니다. coef = 2는 설계 행렬의 2번째 열⁠(grouptreated)을 봅니다.

topTable(fit, coef = 2, number = 6)    # 상위 6개 DEG

              logFC AveExpr     t  P.Value  adj.P.Val    B
FBgn0025111    2.92    6.25  26.8  2.5e-10   1.0e-06  14.2
FBgn0003360   -3.13    7.78 -24.9  5.0e-10   1.1e-06  13.6
FBgn0029167   -2.19    7.87 -24.5  5.7e-10   1.1e-06  13.5
FBgn0039155   -4.62    4.68 -27.9  1.7e-10   1.0e-06  13.3
FBgn0035085   -2.57    5.23 -21.6  1.9e-09   3.0e-06  12.2
FBgn0051092    2.31    3.18  18.0  1.0e-08   1.4e-05  10.4

각 열의 뜻 — logFC⁠(log2 발현 변화, 음수=처리군에서 감소), AveExpr⁠(평균 로그2-CPM 발현), t⁠(moderated t-통계량), P.Value⁠(원시 p값), adj.P.Val⁠(BH로 보정한 q값, FDR), B⁠(log-odds, 차등발현일 사후 로그 승산).

유의 유전자 수를 세고 전체 표를 데이터프레임으로 받습니다.

summary(decideTests(fit))               # adj.P.Val < 0.05 기준 상향/하향
res <- topTable(fit, coef = 2, number = Inf, sort.by = "P")
res$gene <- rownames(res)
sum(res$adj.P.Val < 0.05)               # 유의 DEG 총 개수

       grouptreated
Down            321
NotSig         7249
Up              349

adj.P.Val < 0.05 기준 670개⁠(하향 321 + 상향 349)가 유의합니다. decideTests() 기본은 adjust.method = "BH", p.value = 0.05 — 컷오프를 바꾸려면 인자를 직접 넘기세요⁠(예: p.value = 0.01). FDR의 의미는 p값과 FDR 글에 자세히.

design2 <- model.matrix(~ 0 + group)
colnames(design2) <- c("untreated", "treated")
v2   <- voom(dge, design2)
fit2 <- lmFit(v2, design2)
cm   <- makeContrasts(treated - untreated, levels = design2)
fit2 <- eBayes(contrasts.fit(fit2, cm))
topTable(fit2, number = 4)        # ~ group 결과와 상위 유전자 동일

10. 시각화 — MA plot & volcano

MA plot은 limma의 plotMD() 한 줄, volcano plot은 topTable 결과로 ggplot2를 그립니다.

# (a) MA plot — x=AveExpr(평균발현), y=logFC, 유의 유전자 자동 강조
plotMD(fit, coef = 2, status = decideTests(fit)[, 2],
       main = "limma-voom MA plot (pasilla)")
abline(h = c(-1, 1), col = "grey40", lty = 2)   # |logFC|=1 보조선

# (b) volcano plot — x=logFC, y=-log10(P.Value)
res$sig <- ifelse(res$adj.P.Val < 0.05 & abs(res$logFC) > 1, "DEG", "n.s.")
ggplot(res, aes(logFC, -log10(P.Value), color = sig)) +
  geom_point(size = 1.4, alpha = 0.7) +
  scale_color_manual(values = c(DEG = "#C5573A", n.s. = "grey75")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  labs(x = "log2 fold change", y = expression(-log[10]~italic(P))) +
  theme_bw()

• MA plot — 발현 변화가 평균 발현에 따라 어떻게 분포하는지. 저발현 구간에서 점이 양옆으로 퍼지면 정상입니다.
• volcano plot — 위로 갈수록 유의, 좌우로 갈수록 효과 큼. 오른쪽 위·왼쪽 위 모서리가 신뢰할 만한 DEG입니다. (volcano 해석은 Volcano plot 글 참고)
• limma는 volcanoplot(fit, coef = 2)라는 내장 함수도 있어, ggplot 없이 한 줄로 그릴 수 있습니다.

11. DE 3대장 — DESeq2 vs edgeR vs limma-voom

같은 pasilla 데이터를 세 도구에 넣으면 상위 DEG 목록이 대체로 높게 일치합니다⁠(셋 다 경험적 베이즈로 분산을 안정화하니까). 차이는 모델·검정·잘 맞는 상황에 있습니다.

• 소표본·저카운트⁠(그룹당 2~3개) — DESeq2 / edgeR이 무난합니다. NB로 카운트를 직접 모델링하는 편이 극저카운트에서 안정적입니다.
• 샘플이 많거나 설계가 복잡할 때⁠(시계열·다요인·배치) — limma-voom이 유리합니다. 빠르고, limma의 풍부한 대비 기능과 블로킹⁠(duplicateCorrelation)을 그대로 씁니다.
• 샘플이 아주 적을 때⁠(가중치 추정이 부담) — voom 대신 limma-trend⁠(eBayes(lmFit(logCPM), trend = TRUE))가 더 간단하고 안정적일 수 있습니다.

12. 자주 발생하는 에러 & 해결

• could not find function "DGEList" — library(edgeR)를 안 했음. limma만으로는 부족합니다. edgeR을 먼저 로드하세요⁠(limma도 함께 딸려 옵니다).
• Error in lmFit ... design — model.matrix(~ group)을 빠뜨렸거나, 샘플 순서가 카운트 열과 어긋남. colnames(cts)와 group 순서를 맞추세요.
• voom plot 추세선이 우상향하거나 점이 두 덩어리 — 저발현 유전자 필터가 부족⁠(filterByExpr 누락)하거나 강한 배치 효과⁠(batch effect)의 신호. 필터를 먼저 적용하고, 배치는 설계 행렬에 변수로 넣으세요.
• 유의 유전자가 0개 — TMM·filterByExpr을 건너뛰었거나 coef 지목이 틀림. colnames(design)로 관심 열 위치를 확인⁠(보통 coef = 2)하세요.
• 대비⁠(contrast)가 안 먹힘 — ~ 0 + group 설계에서 makeContrasts(..., levels = design)의 levels 인자를 빠뜨림. 설계 행렬 자체를 넘겨야 합니다.
• logFC 부호가 반대 — factor의 levels 순서 문제. 기준⁠(reference)을 levels 맨 앞에 두세요.

13. 확장 — 복잡한 설계: 블로킹과 duplicateCorrelation

limma-voom이 빛나는 지점이 바로 여기입니다. 같은 개체에서 여러 번 측정했거나⁠(반복측정), 짝지은 설계⁠(paired)처럼 관측이 독립이 아닐 때, limma는 duplicateCorrelation()으로 그 상관 구조를 추정해 모형에 넣습니다 — edgeR·DESeq2로는 다루기 까다로운 경우입니다.

# 예: 환자⁠(patient)를 블록으로 묶은 짝지은 설계
v    <- voom(dge, design)
corfit <- duplicateCorrelation(v, design, block = patient)   # 블록 내 상관 추정
v    <- voom(dge, design, block = patient,                   # 가중치 재계산
             correlation = corfit$consensus)
fit  <- lmFit(v, design, block = patient, correlation = corfit$consensus)
fit  <- eBayes(fit)

배치를 단순 보정만 할 거라면 ~ batch + condition으로 충분합니다. 블록 내 상관까지 살려야 할 때 duplicateCorrelation을 씁니다⁠(공식 limma User's Guide의 권장 패턴).

References

1. Law, C. W., Chen, Y., Shi, W., Smyth, G. K. (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15, R29.
2. Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., Smyth, G. K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47.
3. limma User's Guide (Bioconductor). voom·lmFit·eBayes·topTable·makeContrasts·contrasts.fit·duplicateCorrelation·plotMD·volcanoplot.
4. Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments (moderated t-검정·경험적 베이즈). Statistical Applications in Genetics and Molecular Biology, 3, Article 3.
5. Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data (DGEList·filterByExpr·TMM). Bioinformatics, 26(1), 139–140.
6. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550.
7. Brooks, A. N., et al. (2011). Conservation of an RNA regulatory map between Drosophila and mammals (pasilla 데이터셋). Genome Research, 21, 193–202.

Pipette & Pipeline · A bio portfolio journal

TL;DR — 논문에 "데이터는 GEO에 GSE12345로 올렸다"는 한 줄이 보이면, 그 발현 데이터를 내 R 세션으로 그대로 가져와 재분석할 수 있습니다. GEOquery는 그 다리입니다 — getGEO("GSExxxxx") 한 줄이면 발현행렬⁠(expression matrix)·샘플 메타데이터⁠(metadata)·프로브 주석⁠(annotation)이 한 객체로 들어옵니다. 워크플로는 짧습니다 — GEO 검색 → getGEO → exprs / pData / fData → 그룹 라벨 정리 → 다운스트림⁠(limma·DESeq2).

단, 마이크로어레이⁠(microarray)와 RNA-seq는 GEO에서 오는 형태가 다릅니다 — 어레이는 정규화된 발현행렬이 바로 오지만, RNA-seq는 보통 원시 카운트⁠(raw count)가 보충 파일⁠(supplementary file)로 따로 옵니다. 이 글은 그 차이까지 포함해, 복붙하면 그대로 도는 GEOquery 코드를 따라갑니다. (코드는 Bioconductor GEOquery 비네트 기준입니다.)

🔗 관련 글  ·  받은 데이터를 바로 분석에 투입: DESeq2로 차등발현 분석  ·  edgeR 차등발현 분석  ·  RNA-seq 정규화란?

0. 왜 GEOquery인가 — 1분 직관

블로그의 RNA-seq 분석 글들⁠(DESeq2·edgeR·정규화·PCA)은 모두 "발현 데이터가 이미 손에 있다"는 전제로 시작합니다. 그런데 실전에서 가장 먼저 막히는 건 "그 데이터를 대체 어디서 가져오나"입니다. 답의 8할은 GEO (Gene Expression Omnibus) — NCBI가 운영하는 세계 최대 공개 발현 데이터 저장소입니다.

• 논문 Methods나 Data Availability에 적힌 GSE 번호 하나면, 그 연구의 발현행렬과 샘플 정보를 통째로 받을 수 있습니다. 내 실험 데이터가 없어도 남의 공개 데이터로 재분석·검증·연습이 가능합니다.
• 문제는 GEO 웹페이지에서 손으로 파일을 받아 R로 읽으면 포맷이 제각각이라 번거롭다는 점입니다. GEOquery는 이 과정을 함수 한 줄로 줄입니다 — 다운로드·파싱·객체화를 한 번에.
• 받은 결과는 R의 표준 발현 객체인 ExpressionSet (S4 클래스) 또는 SummarizedExperiment라, limma·DESeq2·edgeR 같은 분석 패키지에 바로 연결됩니다.

1. 준비물 — 설치와 GEO 개념

GEOquery는 Bioconductor 패키지입니다. 어레이 데이터를 다룰 limma, 카운트를 다룰 DESeq2도 함께 깔아 둡니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("GEOquery", "limma", "DESeq2"))  # 최초 1회

library(GEOquery)

먼저 GEO의 세 가지 식별자를 구분해야 코드가 헷갈리지 않습니다. 셋은 포함 관계입니다 — GSE 하나가 여러 GSM을 묶고, 각 GSM은 하나의 GPL 위에서 측정됩니다.

2. getGEO — 시리즈 데이터 받기

핵심 함수는 단 하나, getGEO()입니다. 재현 가능한 공개 예시로 GSE10072 (폐 선암종 종양 vs 정상, Affymetrix HG-U133A 어레이, 107 샘플)를 받아 보겠습니다. 마이크로어레이라 정규화된 발현행렬이 series matrix에 바로 들어 있습니다.

# GSE 시리즈 다운로드 — 처음엔 시간이 걸림⁠(파일을 받아 캐시에 저장)
gse <- getGEO("GSE10072", GSEMatrix = TRUE)   # GSEMatrix=TRUE가 기본⁠(빠른 series matrix)

# ★ getGEO는 '리스트'를 반환 — 한 GSE가 여러 플랫폼을 가질 수 있기 때문
length(gse)        # 플랫폼 개수⁠(여기선 1)
eset <- gse[[1]]   # 첫 번째 플랫폼의 ExpressionSet을 꺼냄
eset               # 객체 요약: assayData(발현), phenoData(메타), featureData(주석)

ExpressionSet (storageMode: lockedEnvironment)
assayData: 22283 features, 107 samples
  element names: exprs
phenoData: sampleNames: GSM254625 GSM254626 ... (107 total)
  varLabels: title geo_accession ... characteristics_ch1
featureData: featureNames: 1007_s_at 1053_at ... (22283 total)

요약을 보면 22283 프로브 × 107 샘플의 발현행렬, 샘플별 표현형⁠(phenotype) 표, 프로브 주석이 한 객체에 들어 있는 게 보입니다. 이게 ExpressionSet의 세 칸입니다.

3. exprs / pData / fData — 세 칸 꺼내기

ExpressionSet에서 실제로 쓸 세 덩어리를 뽑는 접근자⁠(accessor)는 외우기 쉽습니다 — 발현은 exprs, 표현형은 pData, 피처⁠(프로브) 주석은 fData입니다.

# ① 발현행렬 — 행=프로브, 열=샘플 (어레이는 보통 log2 정규화된 강도)
expr <- exprs(eset)
dim(expr)            # 22283 × 107
expr[1:4, 1:3]       # 좌상단만 살짝

# ② 샘플 메타데이터⁠(표현형) — 행=샘플, 열=변수
meta <- pData(phenoData(eset))   # pData(eset)으로 줄여도 동일
dim(meta)            # 107 × (수십 개 변수)
colnames(meta)[1:8]  # title, geo_accession, source_name_ch1, characteristics_ch1 ...

# ③ 피처⁠(프로브) 주석 — 프로브 ID ↔ 유전자 심볼 등
feat <- fData(eset)
colnames(feat)       # ID, "Gene Symbol", "Gene Title", ENTREZ_GENE_ID ...

> expr[1:4, 1:3]
          GSM254625 GSM254626 GSM254627
1007_s_at     9.84      9.71      10.02
1053_at       6.12      6.33       6.05
117_at        7.45      7.20       7.51
121_at        9.03      8.88       9.14

세 줄이 끝입니다. expr은 분석에 넣을 발현행렬, meta는 그룹을 가를 표현형 표, feat은 프로브를 유전자 이름으로 바꿀 사전입니다. 이제 "정상 vs 종양" 같은 비교 라벨만 정리하면 분석 준비가 끝납니다.

4. 그룹 라벨 정리 — 메타데이터에서 비교군 만들기

다운스트림 분석은 "어느 샘플이 어느 그룹인가"를 알아야 합니다. 그 정보는 meta의 characteristics_ch1나 source_name_ch1 같은 열에 자유 텍스트⁠(free text)로 들어 있어, 깔끔한 요인⁠(factor)으로 바꿔 줘야 합니다.

# 어떤 열에 그룹 정보가 있는지 눈으로 확인
meta$source_name_ch1[1:6]
# [1] "Adenocarcinoma of the Lung" "Adenocarcinoma of the Lung"
#     "Normal Lung Tissue" ...

# 텍스트에서 종양/정상 라벨을 뽑아 요인으로 — grepl로 패턴 매칭
group <- ifelse(grepl("Normal", meta$source_name_ch1), "Normal", "Tumor")
group <- factor(group, levels = c("Normal", "Tumor"))   # 대조군을 앞에
table(group)                                             # 그룹별 샘플 수 확인

group
Normal  Tumor
    49     58

5. 탐색 — 차원·결측·로그 변환 점검

분석에 넣기 전, 발현행렬이 멀쩡한지 빠르게 점검합니다. 핵심은 세 가지 — 차원이 맞나, 결측⁠(NA)이 있나, 이미 로그 변환됐나입니다.

dim(expr)                  # 프로브 × 샘플 — 라벨 길이⁠(107)와 열 수가 같아야
sum(is.na(expr))           # 결측값 개수 — 어레이는 종종 NA가 섞임
range(expr, na.rm = TRUE)  # 값의 범위 — log2면 대략 2~16, 원시강도면 수천~수만

# 이미 log 스케일인지 빠른 판정⁠(값의 최댓값이 크면 원시 강도일 가능성)
qx <- quantile(expr, c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
qx

   0%   25%   50%   75%  100%
 2.30  5.41  6.58  8.12 14.92      <- 범위가 2~15 → 이미 log2 정규화됨

위처럼 값이 대략 2~15 사이면 이미 log2로 정규화된 series matrix라 추가 변환이 필요 없습니다. 만약 최댓값이 수만 단위로 크면 원시 강도이니 log2(expr + 1)로 변환합니다. 결측이 많은 프로브는 expr <- expr[complete.cases(expr), ]로 솎아 냅니다.

# 원시 강도일 때만⁠(범위가 매우 클 때) 로그 변환
if (max(expr, na.rm = TRUE) > 100) {
  expr <- log2(expr + 1)
}

# 빠른 시각 점검 — 샘플 간 분포가 비슷한지⁠(boxplot) · 이상 샘플 탐지
boxplot(expr[, 1:10], las = 2, main = "Sample distributions (log2)")

이제 expr⁠(발현행렬) + group⁠(비교 라벨)이 준비됐습니다. 이 둘이 곧 limma⁠(어레이) 또는 DESeq2/edgeR⁠(카운트)의 입력입니다.

6. 마이크로어레이 vs RNA-seq — GEO에서 뭐가 다른가

여기가 GEOquery 실전의 가장 중요한 분기점입니다. getGEO()로 받았을 때 exprs()에 값이 바로 차 있으면 보통 마이크로어레이입니다. 반대로 exprs()가 비거나⁠(0행) 이상하면, 십중팔구 RNA-seq라 원시 카운트가 보충 파일에 따로 있습니다.

RNA-seq — getGEOSuppFiles로 원시 카운트 받기

RNA-seq 시리즈⁠(예: GSE164073)라면 보충 파일을 받아 직접 읽습니다.

# ① 먼저 어떤 보충 파일이 있는지 '목록만' 확인⁠(fetch_files=FALSE)
supp <- getGEOSuppFiles("GSE164073", fetch_files = FALSE)
supp[, c("fname")]          # 예: GSE164073_raw_counts.txt.gz

# ② 실제 다운로드 — 현재 작업폴더 아래 GSE164073/ 폴더로 받아짐
getGEOSuppFiles("GSE164073")

# ③ 받은 gzip 카운트 파일을 행렬로 읽기 — 첫 열을 유전자 ID(행이름)로
counts <- read.delim(gzfile("GSE164073/GSE164073_raw_counts.txt.gz"),
                     row.names = 1, check.names = FALSE)
counts <- as.matrix(counts)
dim(counts)                 # 유전자 × 샘플 정수 행렬 → DESeq2/edgeR 입력

7. 자주 발생하는 에러 & 해결

• exprs(eset)가 0행이거나 비어 있음 — 거의 항상 RNA-seq 시리즈입니다. series matrix엔 카운트가 없으니 getGEOSuppFiles()로 보충 파일을 받아 읽으세요⁠(6절).
• 다운로드 타임아웃 / 중간에 끊김 — 큰 series matrix는 기본 타임아웃⁠(60초)을 넘길 수 있습니다. options(timeout = 600)으로 늘리고 다시 시도하세요. 캐시 덕에 받다 만 파일은 이어받습니다.
• Error: ... VROOM_CONNECTION_SIZE — 시리즈가 아주 클 때 파서 버퍼가 부족하다는 신호입니다. Sys.setenv(VROOM_CONNECTION_SIZE = 500000) 후 getGEO()를 다시 부르면 해결됩니다.
• length(gse)가 2 이상⁠(여러 플랫폼) — 한 GSE가 두 칩에 걸친 경우입니다. 각 원소가 독립 ExpressionSet이니 gse[[1]], gse[[2]]를 따로 다루거나, 분석할 플랫폼⁠(GPL)만 고르세요.
• 프로브 ID는 있는데 유전자 심볼이 없음 — fData(eset)에 "Gene Symbol" 열이 없으면, 플랫폼 주석을 따로 받습니다 — getGEO("GPL96")의 Table에 프로브↔유전자 매핑이 있고, 어레이 기종이면 Bioconductor 주석 패키지⁠(예: hgu133a.db)로도 매핑됩니다.
• 같은 프로브가 여러 개⁠(중복 유전자) — 한 유전자에 프로브가 여러 개일 수 있습니다. 분석 전 발현 평균/최댓값으로 대표 프로브 하나만 남기거나⁠(limma::avereps), 프로브 수준 그대로 분석 후 마지막에 매핑하세요.

8. 다운스트림 연결 — limma / DESeq2 · edgeR

받은 데이터의 형태에 따라 분석 도구가 갈립니다. 정규화된 어레이 강도는 limma, 정수 카운트는 DESeq2/edgeR입니다.

# ── (A) 마이크로어레이 → limma ──  expr=log2 정규화 행렬, group=요인
library(limma)
design <- model.matrix(~ group)             # Normal(기준) vs Tumor
fit <- lmFit(expr, design)                  # 선형모형 적합
fit <- eBayes(fit)                          # 경험적 베이즈 보정
topTable(fit, coef = 2, number = 10)        # 상위 차등발현 프로브⁠(logFC·adj.P)

# ── (B) RNA-seq 카운트 → DESeq2 ──  counts=정수 행렬, coldata=샘플 정보
library(DESeq2)
coldata <- data.frame(group = group, row.names = colnames(counts))
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = coldata,
                              design    = ~ group)
dds <- DESeq(dds)
results(dds)                                # logFC·padj 차등발현 결과

여기까지 오면 GEO의 공개 데이터가 완전한 분석 파이프라인으로 들어왔습니다 — GEO 검색 → GEOquery 다운로드 → 발현행렬·메타 → 그룹 라벨 → limma/DESeq2. 이 블로그의 차등발현·정규화·농축분석 글들이 모두 이 입력에서 출발합니다.

보너스 — 데이터셋 고르기: GEO에서 무엇을 받을까

재분석할 GSE를 고를 땐 세 가지를 봅니다. ① 플랫폼 — 같은 종·같은 기술⁠(어레이 vs RNA-seq)을 골라야 비교가 깔끔합니다. ② 샘플 수 — 그룹당 최소 3개, 통계적 힘을 원하면 그룹당 6개 이상이 안전합니다. ③ 실험 설계 — 비교군이 명확한지⁠(종양 vs 정상, 처리 vs 대조), 메타데이터에 그룹 정보가 잘 적혀 있는지 GEO 페이지에서 미리 확인하세요. GEO 검색창에서 "lung cancer"[Title] AND "Homo sapiens"[Organism] AND "expression profiling by array"[DataSet Type] 같은 필터로 좁히면 원하는 시리즈를 빠르게 찾습니다.

References

1. Davis, S., Meltzer, P. S. (2007). GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 23(14), 1846–1847.
2. GEOquery Vignette (Davis, S., Bioconductor). getGEO·GSEMatrix·exprs·pData·fData·getGEOSuppFiles·getRNASeqData·hasRNASeqQuantifications.
3. Barrett, T., et al. (2013). NCBI GEO: archive for functional genomics data sets — update. Nucleic Acids Research, 41(D1), D991–D995. (GSE·GSM·GPL 구조)
4. Edgar, R., Domrachev, M., Lash, A. E. (2002). Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Research, 30(1), 207–210.
5. Ritchie, M. E., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47.
6. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550.

Pipette & Pipeline · A bio portfolio journal