BioPlayground

一覧へ

機械学習で遺伝子発現データを分析する

scikit-learnで遺伝子発現データをクラスタリングし、PCAで可視化し、腫瘍/正常を分類する方法。

中級
|
120分
|
検証済み (2026-06)
クラスタリングPCA分類遺伝子発現scikit-learnK-means
進捗0/12 (0%)

機械学習で遺伝子発現データを分析する

このトピックを終えたら

scikit-learnで遺伝子発現データをクラスタリング（K-means）し、PCAで次元を圧縮し、腫瘍/正常サンプルを分類できるようになります。

機械学習とは？

実験で経験を積むと、ゲル写真を見ただけで「このバンドパターンはmutationだ」と判断できるようになりますよね？ 機械学習はコンピュータにこの「パターン認識」能力をデータで学習させることです。

バイオでよく使うML：

  • 教師なし学習: ラベルなしでデータのパターンを見つける（クラスタリング、PCA）
  • 教師あり学習: 正解のあるデータで分類モデルを学習（腫瘍 vs 正常）

データ準備：仮想遺伝子発現データ

実際のRNA-seqデータを模した仮想データを作ります。100サンプル（正常50 + 腫瘍50）で20遺伝子の発現量を測定した状況です。

python
import numpy as np
import pandas as pd


np.random.seed(42)


n_samples = 100
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]
sample_labels = ["Normal"] * 50 + ["Tumor"] * 50


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0   # Gene_01~05 上方制御
tumor_expr[:, 15:] -= 2.0  # Gene_16~20 下方制御


expression = np.vstack([normal_expr, tumor_expr])
df = pd.DataFrame(expression, columns=gene_names)
df["Label"] = sample_labels


print(f"データサイズ: {df.shape}")
print(f"サンプル: Normal {sum(df['Label']=='Normal')}, Tumor {sum(df['Label']=='Tumor')}")
print(f"\n最初の5行:")
print(df.head().to_string())


assert df.shape == (100, 21)
assert sum(df["Label"] == "Normal") == 50
assert sum(df["Label"] == "Tumor") == 50

PCA：高次元データを2Dに圧縮

20遺伝子の発現量を人が見れる2次元に圧縮します。PCA（主成分分析）はデータの分散が最も大きい方向を見つけます。

python
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA


np.random.seed(42)
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0


expression = np.vstack([normal_expr, tumor_expr])
labels = ["Normal"] * 50 + ["Tumor"] * 50


# 1. 標準化（平均0、分散1）
scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)


# 2. PCAを適用
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)


print(f"元の次元: {expression.shape[1]}個の遺伝子")
print(f"圧縮後: {X_pca.shape[1]}個の主成分")
print(f"説明分散: PC1={pca.explained_variance_ratio_[0]:.1%}, PC2={pca.explained_variance_ratio_[1]:.1%}")
print(f"合計: {sum(pca.explained_variance_ratio_):.1%}")


assert X_pca.shape == (100, 2)
assert pca.explained_variance_ratio_[0] > pca.explained_variance_ratio_[1]
python
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA


np.random.seed(42)
n_genes = 20


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0


expression = np.vstack([normal_expr, tumor_expr])
labels = np.array(["Normal"] * 50 + ["Tumor"] * 50)


scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)


fig, ax = plt.subplots(figsize=(8, 6))
for label, color in [("Normal", "#2E86AB"), ("Tumor", "#E74C3C")]:
    mask = labels == label
    ax.scatter(X_pca[mask, 0], X_pca[mask, 1], c=color, s=40, alpha=0.7, label=label, edgecolors="white", linewidth=0.5)


ax.set_xlabel(f"PC1 ({pca.explained_variance_ratio_[0]:.1%})", fontsize=12)
ax.set_ylabel(f"PC2 ({pca.explained_variance_ratio_[1]:.1%})", fontsize=12)
ax.set_title("PCA — Normal vs Tumor", fontsize=14, fontweight="bold")
ax.legend(fontsize=11)
ax.grid(True, alpha=0.2)
fig.tight_layout()
fig.savefig("pca_plot.png", dpi=150)
plt.close(fig)


print("pca_plot.png 保存完了")
assert X_pca.shape == (100, 2)

K-meansクラスタリング：教師なし学習

ラベルなしでデータだけからグループを見つけます。K-meansは最も直感的なクラスタリングアルゴリズムです。

python
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans


np.random.seed(42)
n_genes = 20


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0


expression = np.vstack([normal_expr, tumor_expr])
true_labels = [0] * 50 + [1] * 50


scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)


# K-means: 2つのクラスタに分類
kmeans = KMeans(n_clusters=2, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)


# クラスタが実際のラベルとどれくらい一致するか確認
from sklearn.metrics import adjusted_rand_score
ari = adjusted_rand_score(true_labels, clusters)
print(f"クラスタ 0: {sum(clusters == 0)}サンプル")
print(f"クラスタ 1: {sum(clusters == 1)}サンプル")
print(f"Adjusted Rand Index: {ari:.3f} (1.0 = 完全一致)")


assert len(set(clusters)) == 2
assert ari > 0.5
print("K-meansクラスタリングが正常/腫瘍の区別に成功")

分類：教師あり学習

今度は正解（Normal/Tumor）がわかっているデータでモデルを学習し、新しいサンプルを分類します。

python
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, classification_report


np.random.seed(42)
n_genes = 20


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0


X = np.vstack([normal_expr, tumor_expr])
y = np.array([0] * 50 + [1] * 50)  # 0=Normal, 1=Tumor


# 1. データ分割（80%学習、20%テスト）
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)


# 2. 標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)


# 3. Random Forest学習
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train_scaled, y_train)


# 4. 予測と評価
y_pred = clf.predict(X_test_scaled)
accuracy = accuracy_score(y_test, y_pred)


print(f"精度: {accuracy:.1%}")
print(f"\n学習データ: {len(X_train)}個、テストデータ: {len(X_test)}個")
print(f"\n分類レポート:")
print(classification_report(y_test, y_pred, target_names=["Normal", "Tumor"]))


assert accuracy > 0.8
assert len(X_train) == 80
assert len(X_test) == 20

Feature Importance：どの遺伝子が重要か？

Random Forestは分類に最も貢献した遺伝子（feature）を教えてくれます。

python
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier


np.random.seed(42)
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]


normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0


X = np.vstack([normal_expr, tumor_expr])
y = np.array([0] * 50 + [1] * 50)


scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_scaled, y)


importances = clf.feature_importances_
sorted_idx = np.argsort(importances)[-10:]  # 上位10個


fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(range(len(sorted_idx)), importances[sorted_idx], color="#2E86AB")
ax.set_yticks(range(len(sorted_idx)))
ax.set_yticklabels([gene_names[i] for i in sorted_idx], fontsize=10)
ax.set_xlabel("Feature Importance", fontsize=12)
ax.set_title("Top 10 Important Genes", fontsize=14, fontweight="bold")
fig.tight_layout()
fig.savefig("feature_importance.png", dpi=150)
plt.close(fig)


top_gene = gene_names[sorted_idx[-1]]
print(f"最も重要な遺伝子: {top_gene}")
print(f"上位5遺伝子: {[gene_names[i] for i in sorted_idx[-5:]]}")


assert len(importances) == 20
assert sum(importances) > 0.99

やってみよう（Faded Example）

空欄を埋めて、scikit-learnのfit-predictパターンを完成させてください。

穴埋め問題python
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(data)
kmeans = KMeans(n_clusters=)
clusters = kmeans.(X_scaled)

よくあるエラーと解決法

Q: ValueError: could not convert string to float

文字列の列（例：遺伝子名、ラベル）がデータに含まれています。df.select_dtypes(include=[np.number])で数値列のみを選択してください。

Q: 精度が50%しかありません

データがランダムと変わらないということです。標準化（StandardScaler）をしたか、特徴量（遺伝子）が十分に異なるパターンを持っているか確認してください。

Q: ConvergenceWarning: Number of distinct clusters found

K-meansが収束しませんでした。n_init=10、max_iter=300を増やすか、データのスケーリングを確認してください。

おめでとうございます！DevBench共通 + A-1トラックのMVP 5トピックをすべて完了しました。

💬 質問・コメント

0件のコメント

ログインせずに投稿できます。ゲスト投稿は投稿者自身で編集・削除できません。

0/2000

読み込み中...