RNA-seq ヒートマップ — 正規化、クラスタリング、可視化を統合する

このトピックを終えると

教科書で学んだ正規化、階層的クラスタリング、Seabornヒートマップを組み合わせて、RNA-seqカウント行列を受け取り、自動的に正規化し、発現パターンが似ているサンプルと遺伝子をまとめて、論文に使えるヒートマップを作成するツールを自分で作成できるようになります。そして、あの美しい図が実は3つの層が重なり合ったものであることを理解します。

この記事は教育用の一般的な例です。RNA-seqは、今日の遺伝子発現研究における標準的な方法であるため、題材として選びました。

「サンプル順序を変えるだけで別の図になるのはなぜか？」—正規化・並び替えの落とし穴

RNA-seq実験が終わると、このようなカウント行列が得られます。

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 は1万単位の値。ヒートマップ全体がこの2つの遺伝子だけが目立つ真っ赤な線で埋め尽くされます。
サンプル列の順序 を実験時間順に配置したため、腫瘍と正常サンプルが混ざり合い、パターンが見えにくくなります。

これらの問題に直面し始めたときから、あなたは3つの異なる処理が必要であることを理解します。

正規化 — 表現範囲を遺伝子・サンプル間で比較可能にする
並び替え（クラスタリング） — 類似したものを隣に並べることで、パターンが目立つようにする
色マッピング — 値を視覚的に区別できる色に変換する

この3つの層がうまく組み合わさることで、論文に掲載できる図が得られます。これから、この3つの層を一つずつ作成し、組み立てていきます。

まずは完成品を見てみましょう（ブラックボックスを先に実行する）

これから作成するツールは、以下のように使用します。

python

plot_rnaseq_heatmap(
    counts_df,
    top_n_variable=50,        # 分散が大きい上位50個の遺伝子のみ
    cluster_samples=True,
    cluster_genes=True,
    save_path="heatmap.png",
)

この結果には、以下のような情報が含まれます。

text

=== ヒートマップの説明 ===
1. 腫瘍サンプル3個が自動的に左側にグループ化される
2. 正常サンプル3個が自動的に右側にグループ化される
3. 腫瘍で特に発現量の高い遺伝子グループが上側にグループ化される
4. 正常で発現量の高い遺伝子グループが下側にグループ化される
5. 色 = zスコア（サンプル間の分散を基準とした標準化された値）
6. 上部のデンドログラム = サンプルの階層構造
7. 左側のデンドログラム = 遺伝子の階層構造

ファイル保存：heatmap.png（300dpi、論文への挿入が可能）

同じデータを使用すれば、毎回同じ図が出力されます。再現性はファイル名に集約されています。

このツールはどのような部品で構成されているか（部品分解図）

text

RNA-seq ヒートマップパイプライン
   ┌──────────────────────────────────────────────────┐
   │  [入力] counts行列の読み込み ────── 部品：pandas       │  ← 完成品として提供（ツール）
   │              │                                      │
   │              ▼                                      │
   │  [1段階] 正規化（CPM・zスコア）                       │
   │        部品：normalization                            │  ← 自分で作成する ★
   │              │                                      │
   │              ▼                                      │
   │  [2段階] 偏差の大きい遺伝子を選択                            │
   │        部品：pandas変数計算                          │  ← 完成品として提供（ツール）
   │              │                                      │
   │              ▼                                      │
   │  [3段階] サンプル・遺伝子の階層的クラスタリング                    │
   │        部品：hierarchical clustering                   │  ← 自分で作成する ★
   │              │                                      │
   │              ▼                                      │
   │  [出力] seaborn clustermap（ヒートマップ＋デンドログラム）         │  ← 自分で組み立てる ★
   └──────────────────────────────────────────────────┘

部品	どこで学んだか	このツールで行うこと
pandas	`pandas-basics`	counts行列の取り扱い
正規化	`normalization`	遺伝子・サンプルのスケールを統一
階層的クラスタリング	`clustering-hierarchical`	似たものをまとめて並べる
seaborn heatmap	`seaborn-heatmap`	ヒートマップ＋デンドログラムを組み立てる

📌 これらの概念を初めて見る場合は（上部のリンクから）

正規化 — スケールの統一の原理

階層的クラスタリング — デンドログラムの成長規則

seaborn heatmap — 色と値のマッピング

自分で作成する新しい概念は、正規化・階層的クラスタリング・ヒートマップの組み立ての3つ。pandasはツールなので、完成品として提供します。たったの3つです — 認知の限界線内です。

ステップ1：データの準備（提供済みの完成版）

実際のRNA-seqカウント行列を模倣したデータを作成しましょう。実際には、このデータは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):
        # ハウスキーピング遺伝子：2つの条件で類似
        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"]
# ハウスキーピング遺伝子は2つの条件で類似
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正規化 ★ (正規化：サンプル間の比較)

✍️ 自分で記入する箇所。 部品 = 正規化。サンプルごとに総リード数が異なるため、絶対的なカウント数は比較できない。相対的な割合に変換する。

最初の正規化の理由は、サンプル間の総リード数の違いである。あるサンプルはシーケンシングがうまくいき、総3,000万リードを得たが、別のサンプルは2,000万リードしか得られなかった場合、同じ遺伝子のカウント数をそのまま比較することはできない。

CPM（Counts Per Million）: 総リード数を100万で正規化する。

text

CPM(gene, sample) = counts(gene, sample) / total_counts(sample) × 1,000,000

python

def to_cpm(counts):
    """サンプルごとの総カウント数で割り、100万でスケーリングする。"""
    library_sizes = counts.sum(axis=0)     # 各サンプルの総カウント数
    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は腫瘍において実際に高い値を示す
assert tumor_avg["BRCA1"] > normal_avg["BRCA1"]

これでサンプル間の比較が可能になった。しかし、遺伝子間の比較はまだできない。GAPDHのCPMは依然としてBRCA1のCPMよりも100倍大きい。これを解決するのがzスコアである。

🔎 CPMはなぜこれほど単純なのか（ドロー — 正規化） これは、「標準的なサイズにスケーリングする」という最も単純な形の正規化である。実際には、TPM、TMM、DESeq2のmedian-of-ratiosなど、より洗練された方法を使用する。概念的な根底は同じである。比較可能にする。

🤔 自己説明プロンプト CPMは、raw countsをそのまま使用するよりも、なぜ安全なのでしょうか？ 2つのサンプルの総リード数が5,000万と1,000万の場合、raw countsで遺伝子Aの発現をどのように比較できますか？（ヒント：できません。）

ステップ3：Zスコア正規化 ★（正規化：遺伝子間の比較）

✍️ 自分で埋めるセクション。 部品 = 正規化。各遺伝子を自身の平均と標準偏差で標準化し、遺伝子間のスケールの違いをなくします。

各遺伝子行について：

text

z = (value - mean) / std

これにより、各遺伝子の発現が自身の平均から何標準偏差離れているか が表現されます。GAPDHとBRCA1を同じスケールで比較できるようになります。

python

def to_zscore(cpm):
    """遺伝子ごとのzスコア。ログ変換後に標準化（実務慣行）。"""
    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の腫瘍サンプルはすべて正、正常サンプルはすべて負である必要がある（設計上）
assert (zscore.loc["BRCA1", ["tumor1","tumor2","tumor3"]] > 0).all()
assert (zscore.loc["BRCA1", ["normal1","normal2","normal3"]] < 0).all()

これで、すべての遺伝子が同じスケール（平均0、標準偏差1）になります。ヒートマップの色が遺伝子ごとに自由に広がります。この標準化が、ヒートマップの情報密度の決定的な鍵となります。

🤔 自己説明プロンプト Zスコアを適用する前に、なぜ log2(cpm + 1) を最初に適用したのでしょうか？生のCPMにZスコアを直接適用するとどうなるでしょうか？（ヒント：発現はログスケールで自然であり、大きな値が少数例外となり、標準偏差を支配する可能性があります。）

ステップ4：変動の大きい遺伝子のみを選択（完全版）

ヒートマップが情報を伝えるためには、変化のない遺伝子を排除する必要があります。条件に関係なく常に同じ値を示す遺伝子は、白い線を描くだけでノイズになります。

python

def select_top_variable(zscore, top_n=50):
    """zスコアの変動幅が最も大きい上位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)

✍️ 自分で埋める部分。 コンポーネント = 階層的クラスタリング。2つのベクトルの距離を計算し、近いものを順番にまとめて木構造を作成する。

発想：各サンプル（または遺伝子）がベクトルであるため、ベクトル間の距離を測定できる。最も近い2つをまとめてグループにし、そのグループを1つのベクトルとして扱い、さらにまとめていく。結果はデンドログラム — 階層的な木構造になる。

この方法は、scipyに完全な形で実装されている。

🔎 階層的クラスタリングとは（図解 — clustering-hierarchical） 「各点をそれぞれ別のグループとして開始 → 最も近い2つのグループを繰り返し結合 → 1つ残るまで。」結果の木構造を切り分けることで、必要な数のクラスタが得られる。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())
# 腫瘍が3つ、正常細胞が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")]
# 腫瘍が連続した位置にあること
assert max(tumor_positions) - min(tumor_positions) == 2

腫瘍の3つが自動的に横に並んだ。 腫瘍/正常細胞というラベルをクラスタリングに与えなくても、データ自体がそのグループを明らかにしている。これがクラスタリングの価値である — ラベルなしで構造を発見する。

🤔 自己説明プロンプト 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)

この1つの関数が、私たちが作成した正規化、クラスタリング、カラーマッピングをすべて組み立てます。各層の仕組みを理解していれば、パラメータを自信を持って調整できます。

部品を一つに — 完成されたダッシュボードクラス

全体のパイプラインを一つのクラスにまとめます。

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 発現ヒートマップ の縮小版です。実務で使用されるツール（Rのpheatmap、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スコア正規化後には、この2つの方法が似てきますが、代替案を知っておけば状況に合わせて選択できます。
サンプルサイズが小さい場合の危険性: サンプルが3+3のように非常に少ない場合、クラスタリングはノイズに影響を受けやすくなります。論文にこのようなヒートマップを入れる場合は、サンプルサイズが大きくないという事実を結果の解釈に反映する必要があります。

核心: 「正規化は比較可能性を、クラスタリングは構造を、ヒートマップは視覚化を創出します。」これらの3つの層がそれぞれ自身の役割を理解していれば、元のデータから発見までの経路が再現可能になります。皆さんが作り上げたものが、その経路の具体的なものです。

次のステップへ（下部のリンク）

正規化のさまざまな方法の詳細 → 正規化の原理
なぜ階層的クラスタリングはラベルなしで構造を見つけるのか → 階層的クラスタリング
先ほど作成したqPCRの結果をヒートマップにする → qPCR分析の自動化

実際にやってみよう（独立問題）

TPM拡張: 遺伝子長の情報（長さのCSVファイル）を追加で受け取り、CPMの代わりにTPMで正規化するto_tpm(counts, gene_lengths)関数を作成してください。
アノテーション列: clustermapにcol_colorsパラメータとして、腫瘍/正常条件の色分けバーを追加してください。（ヒント：条件ごとの色分けのシリーズを作成して渡す）
選択された遺伝子のみ: ユーザーが関心のある遺伝子のリストを直接指定できるオプションを追加してください（top_nの代わりに）。
挑戦 — 相関ヒートマップ: 遺伝子-遺伝子の相関行列（ピアソンの相関係数）を計算し、それをclustermapで表示します。共発現パターンを見つけるために使用します。

まとめ

「RNA-seqの結果から、ストーリー性のあるヒートマップを作成する」という問題を、以下の3つの要素に分解して解決しました。

正規化により、遺伝子とサンプル間の比較が可能になりました。
階層的クラスタリングにより、ラベルなしで腫瘍/正常グループの構造を自動的に検出しました。
seaborn clustermapにより、その結果を論文に掲載できる視覚的な形式にまとめました。

論文で見た美しいヒートマップの仕組みが、今や明らかになりました。それは、3つの要素がそれぞれの役割を正確に果たした結果にすぎません。皆さんが作成したパイプラインを使って、自分のデータから自分のストーリーを発見することができます。

この記事は一般的な教育例です。実際のRNA-seq解析（DESeq2、edgeR、Bioconductorなど）には、統計的検定、バッチ効果の補正、GOエンリッチメントなどが追加されます。詳細なバージョンは、この基本構造の上に構築するか、検証済みのライブラリに任せることができます。

BioPlayground

RNA-seq ヒートマップ：正規化、クラスタリング、可視化を統合

RNA-seq ヒートマップ — 正規化、クラスタリング、可視化を統合する

このトピックを終えると

「サンプル順序を変えるだけで別の図になるのはなぜか？」—正規化・並び替えの落とし穴

まずは完成品を見てみましょう（ブラックボックスを先に実行する）

このツールはどのような部品で構成されているか（部品分解図）

ステップ1：データの準備（提供済みの完成版）

2段階目 — CPM正規化 ★ (正規化：サンプル間の比較)

ステップ3：Zスコア正規化 ★（正規化：遺伝子間の比較）

ステップ4：変動の大きい遺伝子のみを選択（完全版）

5つのステップで作成 — 階層的クラスタリング ★ (clustering-hierarchical)

6つのステップで作成 — seabornのclustermapを組み立てる ★

部品を一つに — 完成されたダッシュボードクラス

別の道もある（マルチパス反省）

次のステップへ（下部のリンク）

実際にやってみよう（独立問題）

まとめ

💬 質問・コメント