ComplexHeatmap으로 발현 히트맵 그리기 — z-score부터 주석·분할까지 (R 실전)

TL;DR — 논문에서 가장 자주 보는 발현 시각화가 히트맵입니다. R에서는 ComplexHeatmap이 사실상 표준 — 주석 바, 패널 분할, 여러 히트맵 결합까지 됩니다. 이 글은 시뮬레이션 행렬(복붙하면 그대로 돌아감)로 기본 히트맵 → z-score → 색·군집 → 주석 → 분할까지 쌓아 올립니다.

그리고 입문자가 100% 걸리는 첫 함정 — 행(유전자) z-score는 t(scale(t(mat)))입니다. scale()은 열 단위로 표준화하기 때문에, 유전자를 표준화하려면 전치(t)로 감쌌다 풀어줘야 합니다. (코드는 R 4.5.2·ComplexHeatmap 2.26.1에서 실제 실행해 검증했습니다.)

🔗 관련 글  ·  R 시각화 3부작: Volcano plot으로 DEG 시각화  ·  DESeq2로 차등발현 분석  ·  clusterProfiler 농축분석

그림 1. z-score 발현 히트맵 — 열은 Ctrl/Treat로 분할, 위 주석 바로 그룹 표시, 군집이 두 그룹을 가른다.

0. 준비물 — 설치

ComplexHeatmap은 Bioconductor 패키지입니다. 색 함수용 circlize도 함께 받습니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("ComplexHeatmap", "circlize"))   # 최초 1회
library(ComplexHeatmap)
library(circlize)
packageVersion("ComplexHeatmap")   # 2.26.1 (이 글 기준)

1. 예제 데이터 — 재현 가능한 시뮬레이션 행렬

다운로드 없이 따라 하도록, 유전자 45개 × 샘플 10개(대조군 5·처리군 5)를 만듭니다. 그중 15개 유전자를 처리군에서 일부러 바꿔 두면, 나중에 군집이 두 그룹을 갈라 줍니다.

set.seed(1)
n_genes <- 45; n_each <- 5
samples <- c(paste0("Ctrl", 1:n_each), paste0("Treat", 1:n_each))
group   <- factor(rep(c("Ctrl", "Treat"), each = n_each), levels = c("Ctrl", "Treat"))

mat <- matrix(rnorm(n_genes * length(samples), mean = 8, sd = 1.5),
              nrow = n_genes, ncol = length(samples))
treat <- which(group == "Treat")
mat[1:8,  treat] <- mat[1:8,  treat] + 3   # 8개 유전자 처리군에서 상향
mat[9:15, treat] <- mat[9:15, treat] - 3   # 7개 유전자 처리군에서 하향
rownames(mat) <- paste0("gene", 1:n_genes)
colnames(mat) <- samples

dim(mat)        # 45 10
is.matrix(mat)  # TRUE  ← data.frame 아니라 numeric matrix여야 함

2. 가장 단순한 히트맵

Heatmap() 한 줄이면 됩니다. 기본값이 행·열 모두 자동 군집입니다.

Heatmap(mat)

그런데 이대로면 절대 발현량이 큰 유전자가 화면을 지배합니다. '8 근처 vs 11 근처'의 차이는 묻히고, 그냥 밝은 줄·어두운 줄만 보이죠. 그래서 거의 항상 z-score를 씁니다.

3. ⚠️ 핵심 — 행(유전자) z-score는 t(scale(t(mat)))

그림 2. scale()은 열 단위 표준화 — 유전자(행)를 z-score 하려면 전치로 감쌌다 풀어준다.

발현 히트맵의 관례는 유전자마다 표준화(z-score)해서, '절대량'이 아니라 '샘플 간 상대 패턴'을 보는 것입니다. 여기서 첫 함정 —

# scale()은 '열' 단위로 표준화한다. 유전자(행)를 표준화하려면 전치(t)로 감쌌다 풀어준다.
mat_z <- t(scale(t(mat)))

max(abs(rowMeans(mat_z)))   # ~0   (각 행 평균이 0)
range(apply(mat_z, 1, sd))  # 1 1  (각 행 표준편차가 1)

scale(mat)만 쓰면 열(샘플)이 표준화돼 엉뚱한 축이 정규화됩니다. 유전자 패턴을 보려면 반드시 t(scale(t(mat))) — 전치해서 행을 열로 보낸 뒤 표준화하고 다시 전치해 되돌립니다.

⚠️ 분산이 0인 유전자(모든 샘플 값이 같음)는 z-score에서 NaN이 됩니다. 미리 걸러 주세요 — mat <- mat[matrixStats::rowVars(mat) > 0, ].

4. 색과 군집 다듬기

색은 circlize::colorRamp2로 발산형(diverging) 스케일을 만듭니다. z-score엔 음수·0·양수를 가르는 3색이 잘 맞습니다.

col_fun <- colorRamp2(c(-2, 0, 2), c("#2E9FDF", "white", "#E76F51"))

Heatmap(mat_z, name = "z-score", col = col_fun,
        clustering_distance_rows = "pearson",   # 상관 거리(패턴 유사도)
        clustering_method_rows   = "ward.D2",
        show_row_names = TRUE,
        row_names_gp = gpar(fontsize = 7),
        column_title = "Samples", row_title = "Genes")

• colorRamp2(c(-2, 0, 2), ...) — ±2를 넘는 값은 완전 채도로 '클램프'됩니다. z 범위에 맞춰 끊는 점을 정하세요(극단값이 다 뭉개지지 않게).
• clustering_distance_rows = "pearson" — 상관(패턴) 기준 군집. 발현 패턴 비교엔 유클리드보다 자주 쓰입니다.

5. 주석 — 샘플 그룹을 색 바로

HeatmapAnnotation으로 히트맵 위에 그룹 바를 올립니다.

ha <- HeatmapAnnotation(
  group = group,
  col   = list(group = c(Ctrl = "#2E9FDF", Treat = "#E76F51"))
)

Heatmap(mat_z, name = "z-score", col = col_fun, top_annotation = ha,
        clustering_distance_rows = "pearson", clustering_method_rows = "ward.D2",
        row_names_gp = gpar(fontsize = 7))

이제 열 위에 파랑(Ctrl)·주황(Treat) 바가 생겨, 군집이 두 그룹을 제대로 가르는지 한눈에 보입니다.

6. 분할 — 그룹·클러스터별 패널

column_split·row_split으로 히트맵을 패널로 쪼갭니다. 비교가 훨씬 또렷해집니다.

ht <- Heatmap(mat_z, name = "z-score", col = col_fun, top_annotation = ha,
              column_split = group,   # 열을 Ctrl / Treat로 분할
              row_split    = 2,       # 행 덴드로그램을 2덩이로 자름
              clustering_distance_rows = "pearson",
              clustering_method_rows   = "ward.D2",
              row_names_gp = gpar(fontsize = 7))
draw(ht)   # 함수·반복문·png() 안에서는 draw()를 명시해야 그려진다

💡 draw(ht)는 레이아웃된 객체를 돌려줘, row_order(ht)·column_order(ht)로 군집 순서를 꺼낼 수 있습니다(분할하면 리스트로 반환).

7. 실제 DESeq2 결과에 꽂기

시뮬레이션 대신 진짜 데이터는 DESeq2 결과를 이렇게 넣습니다.

# library(DESeq2); library(matrixStats)
vsd <- vst(dds, blind = FALSE)                       # 분산 안정화 변환값
topGenes <- head(order(rowVars(assay(vsd)),          # 변동 큰 유전자 30개
                       decreasing = TRUE), 30)
mat   <- assay(vsd)[topGenes, ]                      # [유전자 × 샘플] 행렬
mat_z <- t(scale(t(mat)))                            # 똑같은 z-score 관용구
Heatmap(mat_z, name = "z-score", col = col_fun,
        top_annotation = HeatmapAnnotation(condition = colData(vsd)$condition))

핵심은 — 원시 카운트가 아니라 vst()/rlog() 변환값을 쓰고, 변동 큰(또는 차등발현) 유전자만 골라 z-score한다는 것입니다.

8. 흔한 실수 & 버전 주의

그림 3. 히트맵 4대 실수 — 행 z-score 잊기·data.frame 입력·분산 0 방치·draw() 누락.

• ❌ scale(mat)로 끝내기 — 열이 표준화됨. 유전자 z-score는 t(scale(t(mat))). (이 글의 1번 함정.)
• ❌ data.frame을 그대로 넣기 — Heatmap()은 numeric matrix가 필요. 유전자 ID를 '열'로 두면 전체가 문자형으로 깨짐 → as.matrix() + ID는 rownames()로.
• ❌ 분산 0 유전자 방치 — z-score에서 NaN → 미리 rowVars > 0으로 거르기.
• ❌ 함수·반복문 안에서 draw() 빠뜨림 — 최상위 Heatmap(...)은 자동으로 그려지지만, 함수·png() 안에선 draw(ht) 명시 필요.
• 🔗 여러 히트맵 결합 — ht1 + ht2는 가로(행 공유), ht1 %v% ht2는 세로(열 공유).
• ⚠️ ComplexHeatmap vs pheatmap — 한 장 빠르게는 pheatmap이 간결하지만, 주석·분할·결합은 ComplexHeatmap만 됩니다.

9. 전체 스크립트 & 핵심 정리

library(ComplexHeatmap); library(circlize)

# 1) 데이터
set.seed(1)
mat <- matrix(rnorm(45 * 10, 8, 1.5), 45, 10)
g <- factor(rep(c("Ctrl","Treat"), each = 5), levels = c("Ctrl","Treat"))
tr <- which(g == "Treat"); mat[1:8, tr] <- mat[1:8, tr] + 3; mat[9:15, tr] <- mat[9:15, tr] - 3
rownames(mat) <- paste0("gene", 1:45)
colnames(mat) <- c(paste0("Ctrl", 1:5), paste0("Treat", 1:5))

# 2) 행 z-score  ← 핵심 관용구
mat_z <- t(scale(t(mat)))

# 3) 색 + 주석 + 분할
col_fun <- colorRamp2(c(-2, 0, 2), c("#2E9FDF", "white", "#E76F51"))
ha <- HeatmapAnnotation(group = g, col = list(group = c(Ctrl = "#2E9FDF", Treat = "#E76F51")))
draw(Heatmap(mat_z, name = "z-score", col = col_fun, top_annotation = ha,
             column_split = g, row_split = 2,
             clustering_distance_rows = "pearson", clustering_method_rows = "ward.D2",
             row_names_gp = gpar(fontsize = 7)))

• ✅ 히트맵 표준 도구 = ComplexHeatmap — 주석·분할·결합이 강점.
• ⚠️ 행 z-score = t(scale(t(mat))) (scale은 열 단위) — 1번 함정.
• ✅ 흐름 — 행렬(numeric) → z-score → colorRamp2 색 → 군집(pearson/ward.D2) → 주석 → 분할 → draw().
• ✅ 실제 데이터는 vst() 변환 + 변동 큰 유전자를 골라 z-score.

10. 다음 학습 추천

• 📊 DESeq2로 차등발현(DEG) 분석 — 히트맵에 넣을 발현 행렬을 만드는 곳 (실전)
• 🌋 Volcano plot으로 DEG 시각화 — 같은 DEG를 다른 그림으로 (비교)
• 🧬 clusterProfiler 농축분석 — DEG에 생물학적 의미 입히기 (응용)

References

1. Gu, Z., Eils, R., Schlesner, M. (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 32(18), 2847–2849.
2. Gu, Z. (2022). Complex heatmap visualization. iMeta, 1(3), e43.
3. Gu, Z., et al. (2014). circlize implements and enhances circular visualization in R. Bioinformatics, 30(19), 2811–2812.
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.

Pipette & Pipeline · A bio portfolio journal — bridging research, data, and industry.