発現クラスタリング — 相関関係と並び替えによる遺伝子共発現グループの探索
このトピックを終えると
教科書で学んだ numpy、相関、ソート を組み合わせて、複数のサンプルにわたる遺伝子発現データから、一緒に動く（共発現）グループを検出する階層的クラスタリングを自分で実装できるようになります。scikit-learnやscipy.clusterが裏で何をしているのか、コードを通じて理解します。
この記事は、教育用一般的な例 です。実戦的なクラスタリングには、pandas/scanpy/scikit-learnを使用します。ここでは、これらのツールの内部概念を正確に扱います。
「この5,000個の遺伝子のうち、どの遺伝子が一緒に動くのか？」—目で見ることはできない
あなたにはRNA-seqデータがあります。5,000個の遺伝子×30サンプルの発現行列です。そして、以下の質問を投げかけます。
- どの遺伝子が一緒にオン・オフになるのか？
- この共発現グループは、どのような生物学的機能に対応するのか？
5,000個の遺伝子の30次元ベクトルを目で比較することは不可能です。自動的に関係性を見つけ、整理する必要があります。
単純なアプローチ: 遺伝子ごとに発現プロファイルをプロットします。5,000個のグラフ。結局、関係性は見えません。
現実的なアプローチ: 2つの軸があります。
- 相関関係: 2つの遺伝子の発現プロファイルがどれだけ一緒に動くかを定量化します。ピアソンの相関係数。
- 階層的クラスタリング: 最も類似したペアから始めて、グループに統合します。結果は、系統樹のようなものになります。
これらを組み合わせると、一緒に動く遺伝子グループが、自然な階層構造として現れます。
ブラックボックスからコンポーネントへ
コンポーネント 1: 発現マトリックスと相関
import numpy as np
def load_expression_matrix(path: str) -> tuple[np.ndarray, list[str], list[str]]: """ 戻り値: (matrix, gene_names, sample_names) matrix の形状: (n_genes, n_samples) """ with open(path) as f: header = f.readline().strip().split("\t") sample_names = header[1:] gene_names = [] rows = [] for line in f: parts = line.strip().split("\t") gene_names.append(parts[0]) rows.append([float(x) for x in parts[1:]]) return np.array(rows), gene_names, sample_names
expr, genes, samples = load_expression_matrix("expression.tsv")print(f"遺伝子 {len(genes)}, サンプル {len(samples)}")print(f"行列の形状: {expr.shape}")次に、相関マトリックス。各要素は2つの遺伝子の間のピアソンの相関です。
def compute_correlation_matrix(expr: np.ndarray) -> np.ndarray: """ expr の形状: (n_genes, n_samples) 戻り値: (n_genes, n_genes) 相関マトリックス """ return np.corrcoef(expr)
corr = compute_correlation_matrix(expr)print(f"corr の形状: {corr.shape}") # (5000, 5000)print(f"対角成分: {corr[0, 0]}") # 1.0 (自分自身)np.corrcoef は、各行を1つの変数として扱います。つまり、expr が (遺伝子, サンプル) の形状の場合、結果は遺伝子-遺伝子の相関になります。
コンポーネント 2: 距離マトリックス
クラスタリングでは、類似度 (similarity) ではなく、距離 (distance) を使用します。相関は -1 ～ +1 の類似度であるため、次の変換で距離にします。
def correlation_to_distance(corr: np.ndarray) -> np.ndarray: """ 距離 = 1 - 相関。 正の相関が強い → 距離が小さい。 負の相関 (逆方向に動く) → 距離が大きい。 """ return 1.0 - corr
dist = correlation_to_distance(corr)注意: 別の定義もあります。sqrt(2 * (1 - corr)) は、ユークリッド距離と一貫性のある定義です。目的に応じて選択してください。
コンポーネント 3: 階層的クラスタリング (Single Linkage)
凝集的なアプローチ: 最初は、各遺伝子が個別のクラスタです。距離が最も近い2つのクラスタを結合することを繰り返します。1つのクラスタが残るまで。
2つのクラスタ間の距離を定義する必要があります。Single linkage は、最も単純な定義です。2つのクラスタの要素のペアの中で最小の距離です。
def hierarchical_clustering_single(dist: np.ndarray) -> list[tuple[int, int, float]]: """ 戻り値: 結合イベントのリスト [(cluster_a, cluster_b, merge_distance), ...] 各結合後、新しいクラスタ ID は元の数 + イベントのインデックスで割り当てられます。 """ n = dist.shape[0] active_clusters = {i: [i] for i in range(n)} cluster_distances = {(i, j): dist[i, j] for i in range(n) for j in range(i + 1, n)} events: list[tuple[int, int, float]] = [] next_id = n while len(active_clusters) > 1: best_pair = min(cluster_distances, key=cluster_distances.get) i, j = best_pair merge_dist = cluster_distances[best_pair] events.append((i, j, merge_dist)) new_cluster = active_clusters[i] + active_clusters[j] del active_clusters[i] del active_clusters[j] active_clusters[next_id] = new_cluster # 新しいクラスタと残りのクラスタの距離を計算 (single linkage = min) new_distances = {} for existing_id in active_clusters: if existing_id == next_id: continue existing = active_clusters[existing_id] min_d = min( dist[a, b] for a in new_cluster for b in existing ) new_distances[(min(existing_id, next_id), max(existing_id, next_id))] = min_d cluster_distances = { k: v for k, v in cluster_distances.items() if i not in k and j not in k } cluster_distances.update(new_distances) next_id += 1 return events時間計算量: 単純な実装では O(n³)。5,000個の遺伝子では処理しきれません。実際には O(n² log n) のアルゴリズムを使用します。このチュートリアルは概念の理解が目的であるため、単純な実装を使用します。
コンポーネント 4: 結果からクラスタを抽出
結合イベントのリストから、特定の数のクラスタを抽出します。
def cut_dendrogram(events: list[tuple[int, int, float]], n_leaves: int, num_clusters: int): parent = list(range(n_leaves + len(events))) def find(x): while parent[x] != x: parent[x] = parent[parent[x]] x = parent[x] return x def union(x, y, new_id): rx, ry = find(x), find(y) parent[rx] = new_id parent[ry] = new_id # 最後の (num_clusters - 1) 個の結合を無視し、それ以前まで union events_to_apply = events[:-num_clusters + 1] if num_clusters > 1 else events next_id = n_leaves for a, b, _ in events_to_apply: union(a, b, next_id) next_id += 1 clusters = {} for leaf in range(n_leaves): root = find(leaf) clusters.setdefault(root, []).append(leaf) return list(clusters.values())
cluster_lists = cut_dendrogram(events, n_leaves=len(genes), num_clusters=10)for i, members in enumerate(cluster_lists): print(f"クラスタ {i}: {len(members)} 遺伝子")Union-Find データ構造を使用すると、結合・検索は O(α(n)) ≈ O(1) で行われます。大規模なデータセットでも役立ちます。
ソーティングとソートされたヒートマップ
クラスタリングの視覚的な成果は、再配置されたヒートマップに現れます。各クラスタの遺伝子が隣接する行になるようにソートします。
def cluster_order(cluster_lists: list[list[int]]) -> list[int]: """各クラスタ内で元のインデックス順序を維持しながら、クラスタごとにソートします。""" result = [] for members in cluster_lists: result.extend(sorted(members)) return result
new_order = cluster_order(cluster_lists)reordered_expr = expr[new_order]この再配置された行列をヒートマップとして描画すると、各クラスタが明確なブロックとして表示されます。
import matplotlib.pyplot as plt
def plot_heatmap(matrix: np.ndarray, ax=None) -> None: if ax is None: _, ax = plt.subplots(figsize=(6, 8)) im = ax.imshow(matrix, cmap="RdBu_r", aspect="auto") plt.colorbar(im, ax=ax)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))plot_heatmap(expr, ax1)ax1.set_title("Original")plot_heatmap(reordered_expr, ax2)ax2.set_title("Clustered")plt.tight_layout()plt.show()この再配置前と後の比較は、クラスタリングの力を視覚的に確認できる決定的なものです。
フェーディング — あなたが埋めるべき3つの空白
空白1：完全連結 / 平均連結
シングル連結の代わりに、完全連結（2つのクラスタ要素ペアのうち、最大の距離）または平均連結（平均距離）を実装する。
def hierarchical_clustering_complete(dist: np.ndarray) -> list[tuple[int, int, float]]: # シングル連結のコードとほとんど同じ。minをmaxに置き換える。 # ... # クラスタ間の距離を計算する部分： # TODO: min(dist[a, b] for ...) を max(dist[a, b] for ...) にする。 passヒント: 関数を1つパラメータ化する。linkage: Callable[[list[float]], float] = min or max or (lambda xs: sum(xs) / len(xs)).
空白2：クラスタの代表遺伝子
各クラスタで、他のメンバーとの平均相関が最も高い遺伝子＝代表遺伝子。
def find_hub_genes( corr: np.ndarray, clusters: list[list[int]], gene_names: list[str]) -> list[tuple[str, float]]: """各クラスタの代表遺伝子とその平均相関。""" hubs = [] for members in clusters: # TODO: 各メンバーについて、他のメンバーとの平均相関を計算する。 # 最も高い値を持つ遺伝子をhubとして選択する。 pass return hubsヒント:
best_score = -1best_gene = Nonefor m in members: score = np.mean([corr[m, other] for other in members if other != m]) if score > best_score: best_score = score best_gene = gene_names[m]空白3：機能的注釈（外部ツールとの連携）
クラスタの遺伝子リストをKEGG pathwayやGO termにマッピングする。ここでは、CSVファイルからマッピングをロードする。
def enrich_clusters( clusters: list[list[str]], gene_to_pathway_csv: str) -> dict: """ 各クラスタで最も多く出現する上位3つのpathway。 """ # TODO: CSVファイルをロードする（gene、pathwayの列）。 # 各クラスタでpathwayの出現頻度を計算する。 # 各クラスタの上位3つのpathwayを返す。 passヒント: from collections import Counter; counter = Counter(); for g in members: counter.update(gene_pathways.get(g, [])).
考察 — 実際の共発現ツールとの違い
scipy.cluster.hierarchy: 実際のツールは、このモジュールの linkage 関数を使用します。C言語で記述された O(n² log n) アルゴリズムです。あなたの O(n³) の単純な実装と比較して、100倍以上高速です。
WGCNA: 遺伝子共発現ネットワークの標準ツールです。符号付き隣接行列、ソフトスレッショルディング、モジュール保存など、高度な概念を使用します。あなたのツールは、これの最小限の骨格です。
scanpy: シングルセルRNA-seqの標準スタックです。クラスタリングは通常、Leidenアルゴリズム — グラフベースのアプローチ — を使用し、あなたが扱った階層的クラスタリングとは異なる系統です。
バッチ効果: 実際の発現データは、異なる実験バッチから取得したデータを結合する必要があります。バッチ効果が、実際の共発現シグナルを汚染します。Combat、Harmonyなどのツールがこの問題を解決します。
動的ツリーカット: あなたの cut_dendrogram は、クラスタ数を指定します。実際のWGCNAは、デンドログラムの形状を見て、動的にカットするアルゴリズムを使用します（Dynamic Tree Cut）。
拡張プロジェクト
1. scipyによる再実装: 独自のクラスタリングを scipy.cluster.hierarchy.linkage で置き換え、パフォーマンスを比較する。
2. デンドログラムの可視化: scipy.cluster.hierarchy.dendrogram を使用してツリー構造を可視化する。独自の events リストを scipy 形式に変換する。
3. GO タームエンリッチメント: gseapy を使用して、各クラスタの GO タームエンリッチメントを実際に実行する。
4. Streamlit ダッシュボード: ユーザーがクラスタ数をスライダーで調整すると、ヒートマップと代表遺伝子がリアルタイムで更新される。
このセクションの構成要素
- [F] numpy: 発現行列、相関行列の操作。
np.corrcoefの実用的な活用。 - [F] 相関関係: Pearson の相関係数を、2 つの遺伝子の類似度を定量的に評価するための基準として使用。
- [F] ソート: クラスタの再配置によるヒートマップのブロック化。Union-Find を使用して結果を抽出。
- [W] matplotlib: ヒートマップの可視化（完成したスクリプトを提供）。
[F] = ユーザーが自分で実装 / [W] = 完成したコードとして提供。