BioPlayground

🧬
목록으로

RNA-seq 히트맵 — 정규화·클러스터링·시각화를 하나로 엮기

지저분한 raw counts에서 시작해 CPM 정규화·z-score·계층적 클러스터링·seaborn 히트맵까지 한 번에. 논문 그림을 스스로 만들 수 있게 됩니다.

중급
|
100분
|
검증 완료 (2026-07)
RNA-seqcounts matrixCPMz-scorehierarchical clusteringdifferential expression샘플 계층 구조발현 히트맵
진행률0/8 (0%)

RNA-seq 히트맵 — 정규화·클러스터링·시각화를 하나로 엮기

이 토픽을 마치면

교과서에서 배운 정규화, 계층적 클러스터링, seaborn 히트맵 을 엮어서, RNA-seq counts 행렬을 받아 자동으로 정규화하고, 발현 패턴이 비슷한 샘플과 유전자를 옆으로 모아, 논문에 넣을 수 있는 히트맵을 만드는 도구를 직접 만들 수 있습니다. 그 예쁜 그림 하나가 사실 세 개의 층이 겹친 것임을 이해하게 됩니다.

이 글은 교육용 일반 예제 입니다. RNA-seq는 오늘 발현 연구에서 표준 방법이라 소재로 삼았습니다.

"샘플 순서만 바꾸면 다른 그림이 나오는데?" — 정규화·정렬의 함정

RNA-seq 실험이 끝나면 이런 counts 행렬이 나옵니다.

text
gene       tumor1  tumor2  tumor3  normal1  normal2  normal3
BRCA1      1520    1830    1650    340      420      380
TP53       850     920     780     650      710      680
MYC        4200    5100    4700    1200     1350     1280
GAPDH      18000   19500   17800   17200    17900    18300
ACTB       15000   16800   15500   14800    15600    15200
...

이 표를 그냥 히트맵으로 색칠하면 어떻게 될까요?

  • BRCA1 은 몇백~몇천 대의 값. 눈에 잘 안 보이는 진한 파란 셀들.
  • GAPDH · ACTB 는 만 단위 값. 히트맵 전체가 이 두 유전자만 눈에 띄는 새빨간 줄로 변합니다.
  • 샘플 열 순서 를 실험 시간순으로 놓았더니, tumor와 normal이 마구 섞여 있어서 패턴이 안 보입니다.

이 문제들을 만나기 시작하는 순간부터 여러분은 세 개의 다른 처리 가 필요함을 깨닫습니다.

  1. 정규화 — 표현 범위를 유전자·샘플 사이에서 비교 가능하게
  2. 정렬(클러스터링) — 비슷한 것끼리 옆에 놓아야 패턴이 눈에 보임
  3. 색상 매핑 — 값을 시각적으로 구별 가능한 색으로 변환

이 세 층이 잘 맞물려야 논문에 넣을 수 있는 그림이 나옵니다. 지금부터 이 세 층을 하나씩 만들어 조립합니다.

먼저 완성품을 봅시다 (블랙박스 먼저 돌리기)

우리가 만들 도구는 이렇게 씁니다.

python
plot_rnaseq_heatmap(
    counts_df,
    top_n_variable=50,        # 편차 큰 상위 50개 유전자만
    cluster_samples=True,
    cluster_genes=True,
    save_path="heatmap.png",
)

결과는 이런 이야기를 담습니다.

text
=== 히트맵 이야기 ===
1. tumor 샘플 3개가 자동으로 왼쪽에 모임
2. normal 샘플 3개가 자동으로 오른쪽에 모임
3. tumor에서 특히 강한 유전자 그룹이 위쪽으로 모임
4. normal에서 강한 유전자 그룹이 아래쪽으로 모임
5. 색상 = z-score (샘플 간 편차 기준의 표준화 값)
6. 상단 덴드로그램 = 샘플 계층 구조
7. 좌측 덴드로그램 = 유전자 계층 구조

파일 저장: heatmap.png (300dpi, 논문 삽입 가능)

같은 데이터로 매번 동일한 그림이 나옵니다. 재현성이 파일 이름 하나에 담깁니다.

이 도구는 어떤 부품으로 조립되는가 (부품 분해도)

text
RNA-seq 히트맵 파이프라인
   ┌──────────────────────────────────────────────────┐
   │  [입력] counts 행렬 로딩 ────── 부품: pandas       │  ← 완성 제공 (도구)
   │              │                                      │
   │              ▼                                      │
   │  [1단계] 정규화 (CPM · z-score)                       │
   │        부품: normalization                            │  ← 직접 만든다 ★
   │              │                                      │
   │              ▼                                      │
   │  [2단계] 편차 큰 유전자 선택                            │
   │        부품: pandas 변수 계산                          │  ← 완성 제공 (도구)
   │              │                                      │
   │              ▼                                      │
   │  [3단계] 샘플·유전자 계층적 클러스터링                    │
   │        부품: hierarchical clustering                   │  ← 직접 만든다 ★
   │              │                                      │
   │              ▼                                      │
   │  [출력] seaborn clustermap (히트맵 + 덴드로그램)         │  ← 직접 조립한다 ★
   └──────────────────────────────────────────────────┘
부품어디서 배웠나이 도구에서 하는 일
pandaspandas-basicscounts 행렬 다루기
정규화normalization유전자·샘플 스케일 통일
계층적 클러스터링clustering-hierarchical비슷한 것끼리 옆에 모으기
seaborn heatmapseaborn-heatmap히트맵 + 덴드로그램 조립

📌 이 개념들 처음 본다면 (상단 진입 링크)

직접 만드는 새 개념은 정규화 · 계층적 클러스터링 · 히트맵 조립 3개. pandas는 도구라 완성본으로 드립니다. 딱 3개죠 — 인지 한계선 안입니다.

만들기 1단계 — 데이터 준비 (완성 제공)

실제 RNA-seq counts 행렬을 흉내낸 데이터를 만듭시다. 실무에서는 이 데이터가 salmon·featureCounts·HTSeq 등에서 CSV로 나옵니다.

python
import numpy as np
import pandas as pd


rng = np.random.default_rng(42)


def make_rnaseq_counts():
    genes = ["BRCA1", "TP53", "MYC", "KRAS", "EGFR", "PIK3CA", "PTEN", "APC",
             "GAPDH", "ACTB", "B2M", "HPRT1", "TBP"]  # 마지막 5개는 하우스키핑
    tumor_samples = [f"tumor{i+1}" for i in range(3)]
    normal_samples = [f"normal{i+1}" for i in range(3)]
    samples = tumor_samples + normal_samples


    counts = np.zeros((len(genes), len(samples)), dtype=int)
    for gi, gene in enumerate(genes):
        # 하우스키핑: 두 조건에서 비슷
        base_tumor  = 1000 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        base_normal = 300 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        if gene in {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}:
            base_tumor = base_normal = 15000  # 하우스키핑
        for si, sample in enumerate(samples):
            base = base_tumor if sample.startswith("tumor") else base_normal
            counts[gi, si] = max(0, int(base * rng.lognormal(0, 0.15)))
    return pd.DataFrame(counts, index=genes, columns=samples)


counts = make_rnaseq_counts()
print(counts)


# 검증
assert counts.shape == (13, 6)
assert (counts >= 0).all().all()
# tumor의 BRCA1이 normal보다 높아야 함 (설계상)
assert counts.loc["BRCA1", "tumor1"] > counts.loc["BRCA1", "normal1"]
# 하우스키핑은 두 조건이 비슷
housekeeping_tumor = counts.loc["GAPDH", ["tumor1","tumor2","tumor3"]].mean()
housekeeping_normal = counts.loc["GAPDH", ["normal1","normal2","normal3"]].mean()
assert abs(housekeeping_tumor / housekeeping_normal - 1) < 0.3

만들기 2단계 — CPM 정규화 ★ (normalization: 샘플 간 비교)

✍️ 직접 채우는 구간. 부품 = 정규화. 샘플마다 총 read 수가 다르므로 절대 counts는 비교 불가. 상대적 비율로 바꾼다.

첫 번째 정규화 이유는 샘플 간 총 read 수의 차이 입니다. 어떤 샘플은 시퀀싱이 잘 되어 총 3천만 read를 얻었고, 어떤 샘플은 2천만 read만 얻었다면, 같은 유전자의 counts를 그대로 비교할 수 없습니다.

CPM (Counts Per Million): 총 read 수를 백만으로 정규화.

text
CPM(gene, sample) = counts(gene, sample) / total_counts(sample) × 1,000,000
python
def to_cpm(counts):
    """샘플별 총 counts로 나눠 백만으로 스케일링."""
    library_sizes = counts.sum(axis=0)     # 각 샘플의 총 counts
    cpm = counts.div(library_sizes, axis=1) * 1_000_000
    return cpm


cpm = to_cpm(counts)
print(cpm.round(1).head())


# 검증: 각 샘플의 CPM 합이 정확히 1,000,000
assert np.allclose(cpm.sum(axis=0).values, 1_000_000)
# 비교: raw counts로는 못했던 비교가 CPM으로 가능
tumor_avg = cpm[["tumor1","tumor2","tumor3"]].mean(axis=1)
normal_avg = cpm[["normal1","normal2","normal3"]].mean(axis=1)
# BRCA1이 tumor에서 실제로 더 높음
assert tumor_avg["BRCA1"] > normal_avg["BRCA1"]

이제 샘플 간 비교가 가능해졌습니다. 하지만 유전자 간 비교 는 아직 안 됩니다. GAPDH의 CPM은 여전히 BRCA1의 CPM보다 100배 큽니다. 이걸 해결하는 게 z-score입니다.

🔎 CPM은 왜 이렇게 단순한가 (드로어 — normalization) "표준 크기로 스케일링" 이라는 정규화의 가장 단순한 형태입니다. 실무에서는 TPM · TMM · DESeq2의 median-of-ratios 등 더 정교한 방법을 씁니다. 개념적 뿌리는 같습니다 — 비교 가능하게 만들기.

🤔 자기설명 프롬프트 CPM은 왜 raw counts를 그대로 쓰는 것보다 안전할까요? 두 샘플의 총 read 수가 5천만과 1천만일 때, raw counts로 유전자 A의 발현을 어떻게 비교할 수 있을까요? (힌트: 못합니다.)

만들기 3단계 — Z-score 정규화 ★ (normalization: 유전자 간 비교)

✍️ 직접 채우는 구간. 부품 = 정규화. 각 유전자를 자기 평균·편차로 표준화해 유전자 간 스케일 차이 없앰.

각 유전자 행에 대해:

text
z = (value - mean) / std

이렇게 하면 각 유전자의 발현이 자기 평균에서 몇 배의 표준편차만큼 떨어져 있는가 로 표현됩니다. GAPDH와 BRCA1이 같은 스케일에서 비교 가능해집니다.

python
def to_zscore(cpm):
    """유전자별 z-score. log 변환 후 표준화 (실무 관행)."""
    log_cpm = np.log2(cpm + 1)                        # +1은 log(0) 방지
    means = log_cpm.mean(axis=1)
    stds = log_cpm.std(axis=1, ddof=0)
    # 편차 0(모든 값 같음)인 유전자는 안전하게 0으로
    stds = stds.replace(0, 1)
    return log_cpm.sub(means, axis=0).div(stds, axis=0)


zscore = to_zscore(cpm)
print(zscore.round(2))


# 검증: 각 유전자의 평균 ≈ 0, 편차 ≈ 1
row_means = zscore.mean(axis=1)
row_stds = zscore.std(axis=1, ddof=0)
assert np.allclose(row_means.values, 0, atol=1e-9)
# 편차 있는 유전자만 검사 (편차 0으로 나눈 것 제외)
active_genes = cpm.std(axis=1) > 0
assert np.allclose(row_stds[active_genes].values, 1, atol=0.05)
# BRCA1의 tumor 샘플들은 양수, normal은 음수여야 함 (설계상)
assert (zscore.loc["BRCA1", ["tumor1","tumor2","tumor3"]] > 0).all()
assert (zscore.loc["BRCA1", ["normal1","normal2","normal3"]] < 0).all()

이제 모든 유전자가 같은 스케일 (평균 0, 편차 1)에 있습니다. 히트맵 색상이 유전자별로 자유롭게 스팬됩니다. 이 표준화가 히트맵 정보 밀도의 결정적 열쇠입니다.

🤔 자기설명 프롬프트 Z-score 전에 왜 log2(cpm + 1)을 먼저 씌웠을까요? Raw CPM에서 그냥 z-score를 하면 어떻게 될까요? (힌트: 발현이 로그 스케일로 자연스럽고, 큰 값이 소수 이상치가 편차를 지배할 수 있습니다.)

만들기 4단계 — 편차 큰 유전자만 선택 (완성 제공)

히트맵이 정보를 담으려면 변화 없는 유전자 는 없애야 합니다. 조건과 무관하게 항상 같은 값인 유전자는 흰색 줄만 그리고 노이즈만 됩니다.

python
def select_top_variable(zscore, top_n=50):
    """z-score 변동 폭 상위 top_n 유전자."""
    # 여기서는 편차 대신 range(max-min)을 씀 — 강한 위/아래 편차 유전자 우선
    variability = zscore.max(axis=1) - zscore.min(axis=1)
    top_genes = variability.sort_values(ascending=False).head(top_n).index
    return zscore.loc[top_genes]


top = select_top_variable(zscore, top_n=8)
print(top.index.tolist())


# 검증: 8개 유전자만 남음
assert len(top) == 8
# 하우스키핑(변화 적음)이 상위 8에서 배제되어야
housekeeping = {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}
selected_housekeeping = set(top.index) & housekeeping
# 대부분 배제 (일부는 노이즈로 포함될 수 있음)
assert len(selected_housekeeping) <= 2

만들기 5단계 — 계층적 클러스터링 ★ (clustering-hierarchical)

✍️ 직접 채우는 구간. 부품 = 계층적 클러스터링. 두 벡터의 거리를 계산해 가까운 것끼리 순차적으로 묶어 나무를 만든다.

발상: 각 샘플(또는 유전자)이 벡터이므로 벡터 간 거리 를 잴 수 있습니다. 가장 가까운 두 개를 묶어 그룹으로 만들고, 그 그룹을 하나의 벡터로 취급해 계속 묶어 나갑니다. 결과가 덴드로그램(dendrogram) — 위계적 나무입니다.

이 방법은 scipy에 완결된 형태로 있습니다.

🔎 계층적 클러스터링이란 (드로어 — clustering-hierarchical) "각 점을 자기 그룹으로 시작 → 가장 가까운 두 그룹을 반복 병합 → 하나 남을 때까지." 결과 나무를 자르면 원하는 개수의 클러스터가 나옵니다. k-means와 달리 클러스터 개수를 미리 정할 필요가 없습니다.

python
from scipy.cluster.hierarchy import linkage, dendrogram


def cluster_order(data, axis="rows", method="average", metric="euclidean"):
    """계층적 클러스터링으로 재정렬된 인덱스 반환."""
    matrix = data.values if axis == "rows" else data.values.T
    Z = linkage(matrix, method=method, metric=metric)
    # leaves_list는 나무 잎(원본)의 최적 정렬 순서
    from scipy.cluster.hierarchy import leaves_list
    order = leaves_list(Z)
    if axis == "rows":
        return data.index[order]
    return data.columns[order]


# 샘플과 유전자 각각 클러스터 순서로 재정렬
sample_order = cluster_order(top, axis="cols")
gene_order   = cluster_order(top, axis="rows")


reordered = top.loc[gene_order, sample_order]
print("샘플 순서:", list(sample_order))
print("유전자 순서:", list(gene_order))


# 검증: 재정렬된 표의 값 총합은 원본과 같음
assert np.allclose(reordered.values.sum(), top.values.sum())
# tumor 3개와 normal 3개가 각각 뭉쳐야 함 (설계상)
sample_names = list(sample_order)
tumor_positions = [i for i, s in enumerate(sample_names) if s.startswith("tumor")]
normal_positions = [i for i, s in enumerate(sample_names) if s.startswith("normal")]
# tumor들이 연속된 위치에 있어야
assert max(tumor_positions) - min(tumor_positions) == 2

tumor 3개가 자동으로 옆에 모였습니다. 우리가 tumor/normal이라고 라벨을 클러스터링에 안 알려줬는데도 데이터 자체가 그 그룹을 드러냅니다. 이게 클러스터링의 가치입니다 — 라벨 없이 구조를 발견합니다.

🤔 자기설명 프롬프트 method="average", metric="euclidean" 을 골랐습니다. method="ward" 로 바꾸면 어떻게 될까요? metric="correlation" 이라면? 실제로 바꿔서 결과를 비교해보세요.

만들기 6단계 — seaborn clustermap 조립 ★

seaborn의 clustermap 이 위의 클러스터링과 히트맵을 한 번에 그려줍니다. 우리가 앞에서 손으로 조립한 원리를 알고 있으니, 이 함수의 파라미터가 무엇을 하는지 정확히 이해할 수 있습니다.

python
import matplotlib
matplotlib.use("Agg")   # Pyodide/서버 호환
import matplotlib.pyplot as plt
import seaborn as sns


def plot_rnaseq_heatmap(counts_df, top_n=50, save_path=None):
    """전체 파이프라인: 정규화 → 상위 변동 → clustermap."""
    cpm = to_cpm(counts_df)
    z = to_zscore(cpm)
    top = select_top_variable(z, top_n=top_n)


    g = sns.clustermap(
        top,
        cmap="RdBu_r",           # 빨강-파랑 (양음 명확)
        center=0,                # 0 (평균)을 중앙에
        vmin=-2, vmax=2,         # 색상 범위 고정
        method="average",
        metric="euclidean",
        figsize=(10, max(5, top_n * 0.15)),
        cbar_kws={"label": "z-score (log2 CPM)"},
    )
    if save_path:
        g.savefig(save_path, dpi=300, bbox_inches="tight")
    return g


# 실행 (Colab/Jupyter)
g = plot_rnaseq_heatmap(counts, top_n=8, save_path=None)
plt.close("all")


# 검증: 반환된 grid가 seaborn clustermap 객체
from seaborn.matrix import ClusterGrid
assert isinstance(g, ClusterGrid)
# 재정렬된 데이터가 예상 형태
assert g.data2d.shape == (8, 6)

이 한 함수가 우리가 만든 정규화·클러스터링·색상 매핑을 전부 조립합니다. 각 층의 원리를 알고 있으면 파라미터를 자신 있게 튜닝할 수 있습니다.

부품을 하나로 — 완성된 대시보드 클래스

전체 파이프라인을 하나의 클래스로 묶습니다.

python
class RNAseqDashboard:
    def __init__(self, counts_df: pd.DataFrame):
        self.counts = counts_df
        self.cpm = to_cpm(counts_df)
        self.zscore = to_zscore(self.cpm)


    def summary(self) -> pd.DataFrame:
        return pd.DataFrame({
            "n_genes": [len(self.counts)],
            "n_samples": [self.counts.shape[1]],
            "total_reads_mean": [self.counts.sum(axis=0).mean()],
            "top_variable_gene": [
                (self.zscore.max(axis=1) - self.zscore.min(axis=1)).idxmax()
            ],
        })


    def plot(self, top_n=50, save_path=None):
        return plot_rnaseq_heatmap(self.counts, top_n=top_n, save_path=save_path)


dash = RNAseqDashboard(counts)
print(dash.summary())


# 검증
assert dash.summary()["n_genes"].iloc[0] == 13
assert dash.summary()["n_samples"].iloc[0] == 6
# 가장 편차 큰 유전자는 진짜 조건 간 차이가 큰 것이어야 함 (BRCA1 · MYC · KRAS 중 하나)
assert dash.summary()["top_variable_gene"].iloc[0] in {"BRCA1", "MYC", "KRAS"}

이 도구가 오늘 논문에서 자주 보는 RNA-seq 발현 히트맵 의 축소판입니다. 실무 도구(pheatmap in R, ComplexHeatmap 등)는 여기에 색상 어노테이션, 커스텀 거리 함수, 통계 검정 오버레이 등을 추가한 것뿐입니다.

다른 길도 있다 (멀티패스 성찰)

  • TPM vs CPM: CPM은 샘플 크기만 정규화합니다. TPM (Transcripts Per Million) 은 유전자 길이도 함께 정규화해서 유전자 간 비교에도 더 정확합니다. 언제 무엇: 샘플 간 비교만 = CPM / 유전자 간 절대 비교 = TPM.
  • DESeq2 정규화: 차등 발현 분석(differential expression) 을 목적으로 하면 CPM/TPM 대신 DESeq2의 median-of-ratios 나 TMM (edgeR) 를 씁니다. 이들이 통계 검정 가정에 더 잘 맞습니다.
  • k-means 대안: 계층적 클러스터링은 O(n²) 메모리라 유전자가 수만 개면 무리입니다. 그때는 k-means 나 UMAP + HDBSCAN 으로 대체합니다. 라벨 없는 구조 발견이라는 목적은 같습니다.
  • 거리 함수 선택: euclidean 은 값의 절대 크기에 민감합니다. correlation 은 패턴 모양(오르내림)에만 집중합니다. z-score 정규화 후엔 두 방법이 비슷해지지만, 대안을 알고 있으면 상황에 맞게 선택할 수 있습니다.
  • 표본 크기가 작을 때 위험: 샘플이 3+3처럼 매우 적으면 클러스터링이 노이즈에 흔들립니다. 논문에 이런 히트맵을 넣을 때는 표본이 크지 않다는 사실을 결과 해석에 반영 해야 합니다.

핵심: "정규화가 비교 가능성을, 클러스터링이 구조를, 히트맵이 눈을 만듭니다." 세 층이 각각 자기 역할을 알고 있으면 원본 데이터에서 발견까지의 경로가 재현 가능합니다. 여러분이 방금 만든 것이 그 경로의 실체입니다.

다음 단계로 (하단 출구 링크)

직접 해보기 (독립 문제)

  1. TPM 확장: 유전자 길이 정보(길이 CSV)를 추가로 받아 CPM 대신 TPM으로 정규화하는 to_tpm(counts, gene_lengths) 를 짜세요.
  2. 어노테이션 컬럼: clustermap 에 col_colors 파라미터로 tumor/normal 조건 색상 막대를 추가하세요. (힌트: 조건별 색상 시리즈를 만들어 전달)
  3. 선택된 유전자만: 사용자가 관심 유전자 리스트를 직접 지정하는 옵션을 추가하세요 (top_n 대신).
  4. 도전 — 상관 히트맵: 유전자-유전자 상관 행렬(Pearson correlation)을 계산해서 그것을 clustermap하세요. co-expression 패턴 발견에 씁니다.

정리

우리는 "RNA-seq 결과에서 이야기가 있는 히트맵을 만드는" 문제를, 세 개의 부품으로 정복했습니다.

  • 정규화 가 유전자·샘플 사이의 비교를 가능하게 만들었습니다.
  • 계층적 클러스터링 이 라벨 없이 tumor/normal 그룹 구조를 자동 발견했습니다.
  • seaborn clustermap 이 그 결과를 논문에 넣을 수 있는 시각적 형태로 조립했습니다.

논문에서 보던 그 예쁜 히트맵의 마법이 이제 벗겨졌습니다. 세 개의 층이 정확히 자기 몫을 한 결과일 뿐입니다. 여러분이 만든 파이프라인으로 여러분의 데이터에서 여러분의 이야기를 발견할 수 있습니다.

이 글은 일반 교육 예제 입니다. 실무 RNA-seq 분석(DESeq2, edgeR, Bioconductor 등)은 여기에 통계 검정, batch effect 보정, GO enrichment 등이 추가된 것입니다. 그 상세한 버전은 이 뼈대 위에 여러분이 붙이거나, 검증된 라이브러리에 맡기면 됩니다.