0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@bacbac(Enoki)

セルフ夏期講習4日目：RNA-seqデータをPCAで可視化し、癌種を分類する

0
Posted at

はじめに

セルフ夏期講習の4日目として、RNA-seqデータを用いた癌種推定に取り組んだ。

今回使用したのは、UCI Machine Learning RepositoryのGene Expression Cancer RNA-Seqである。801検体について20,531遺伝子の発現量が記録され、目的変数は5種類の癌種となっている。

8月13日分の目標は、データの構造を理解し、PCAによる可視化と学習・テストデータへの分割を行うことだった。実際には8月14日に予定していたロジスティック回帰と交差検証まで進めることができた。

今回の進捗

  • RNA-seqデータと癌種ラベルの読み込み
  • 値の範囲と欠損値の確認
  • 癌種ごとの検体数を確認
  • 標準化とPCAによる可視化
  • 癌種比を維持した学習・テストデータへの分割
  • 標準化＋ロジスティック回帰によるベースライン作成
  • 5分割交差検証
  • 混同行列と誤分類検体の確認

データの構造

遺伝子発現量と癌種ラベルを、それぞれPandasのDataFrameとして読み込んだ。

import pandas as pd

df_data = pd.read_csv("data.csv", index_col=0)
df_label = pd.read_csv("labels.csv", index_col=0)

X = df_data
y = df_label["Class"]

Xは説明変数であり、1行が1検体、1列が1遺伝子に対応する。yは目的変数で、各検体の癌種が入っている。

癌種ごとの検体数は以下の通りだった。

癌種 検体数
BRCA 300
KIRC 146
LUAD 141
PRAD 136
COAD 78

BRCAが最も多く、COADが最も少ない。後でデータを分割するときには、この構成比をなるべく維持する必要がある。

PCAで遺伝子発現量の分布を見る

20,531個の遺伝子をそのまま図にすることはできない。そこで、StandardScalerで遺伝子ごとに標準化した後、PCAで2～3個の主成分へ圧縮した。

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

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

pca = PCA(n_components=3)
X_pca = pca.fit_transform(X_scaled)

df_pca = pd.DataFrame(
    X_pca,
    columns=["PC1", "PC2", "PC3"],
    index=X.index,
)
df_pca["Class"] = y

PCAは、データのばらつきをなるべく保ちながら、多数の特徴量を少数の軸へ要約する方法である。癌種ラベルはPCAの計算には使用せず、散布図を癌種ごとに色分けするためだけに使用した。

import matplotlib.pyplot as plt
import seaborn as sns

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

sns.scatterplot(
    data=df_pca,
    x="PC1",
    y="PC3",
    hue="Class",
    palette="tab10",
    alpha=0.7,
    s=50,
    ax=axes[0],
)

sns.scatterplot(
    data=df_pca,
    x="PC2",
    y="PC3",
    hue="Class",
    palette="tab10",
    alpha=0.7,
    s=50,
    ax=axes[1],
)

axes[0].set_title("PC1 vs PC3")
axes[1].set_title("PC2 vs PC3")
plt.tight_layout()
plt.show()

PC1とPC2の図では、KIRCが他の癌種から比較的離れた位置に集まっていた。一方、ほかの癌種には重なりが見られた。PC1とPC3の図では、LUADが広い範囲に分布し、ほかの癌種を覆うように重なっていた。

ただし、この結果だけから「癌種を分類できない」と判断することはできない。PCAは癌種を分ける方向ではなく、データ全体の分散が大きい方向を求める教師なし学習だからである。また、図に表示しているのは20,531次元の情報を要約した最初の2～3主成分だけである。

学習用とテスト用に分割する

分類モデルの評価に使用するため、全データを学習用80%、テスト用20%に分割した。

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X,
    y,
    test_size=0.2,
    random_state=42,
    stratify=y,
)

stratify=yを指定すると、癌種の構成比が学習データとテストデータでなるべく同じになる。COADのように検体数が少ないクラスがある場合には特に重要である。

各変数の意味は以下の通りである。

  • X_train：モデルの学習に使う遺伝子発現量
  • X_test：学習後のモデルを評価する遺伝子発現量
  • y_train：X_trainに対応する正解の癌種
  • y_test：X_testに対応する正解の癌種

ロジスティック回帰による癌種分類

少し余裕があったため、翌日に予定していたベースラインモデルの作成にも進んだ。

今回のデータは、約640件の学習検体に対して20,531個の特徴量がある。過学習とデータリークに注意し、標準化とL2正則化ロジスティック回帰をPipelineにまとめた。

from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

model = Pipeline([
    ("scaler", StandardScaler()),
    ("classifier", LogisticRegression(
        C=1.0,
        max_iter=3000,
        random_state=42,
    )),
])

model.fit(X_train, y_train)
y_test_pred = model.predict(X_test)

Pipelineを使うことで、StandardScalerは学習データだけから平均と標準偏差を求める。テストデータを含む全データで先に標準化すると、テストデータの情報が前処理へ混入するため注意が必要である。

ロジスティック回帰は、各遺伝子の発現量に重みを付けて癌種ごとの線形スコアを作り、正解癌種に高い確率を割り当てられるように係数を学習する。名前に「回帰」と付いているが、今回は5癌種の分類に使用している。

最初のホールドアウト評価では、訓練Accuracy、テストAccuracy、テストBalanced Accuracyがすべて1.0だった。

5分割交差検証

1回の分割だけでは、偶然分類しやすいテストデータになった可能性がある。そこで、癌種比を保ちながら検証担当のデータを5回入れ替えるStratifiedKFoldを使用した。

from sklearn.model_selection import StratifiedKFold, cross_validate

cv = StratifiedKFold(
    n_splits=5,
    shuffle=True,
    random_state=42,
)

scores = cross_validate(
    model,
    X,
    y,
    cv=cv,
    scoring={
        "accuracy": "accuracy",
        "balanced_accuracy": "balanced_accuracy",
        "macro_f1": "f1_macro",
    },
    n_jobs=-1,
    return_train_score=True,
)

結果は以下の通りだった。

評価指標 結果
各foldの訓練Accuracy 1.0000, 1.0000, 1.0000, 1.0000, 1.0000
各foldの検証Accuracy 1.0000, 1.0000, 0.9875, 0.99375, 0.9875
平均検証Accuracy 0.99375
検証Accuracyの標準偏差 0.00559
平均Balanced Accuracy 0.99407
平均Macro-F1 0.99381

検証Accuracyは98.75～100%であり、データの分け方が変わっても安定して高い性能が得られた。訓練Accuracyはすべて1.0だったが、検証性能との差は小さい。少なくとも今回のデータセット内では、顕著な過学習は確認されなかった。

今回はデータ全体に対する交差検証としてXとyを渡した。今後、最初に確保したX_testとy_testを完全に未使用の最終評価データとして扱う場合は、モデル選択の交差検証にはX_trainとy_trainだけを渡し、設定を決めた後にテストデータを一度だけ評価する。

PCA図では癌種同士に重なりがあったにもかかわらず、ロジスティック回帰では高い精度が得られた。これは、PCA図が2～3主成分だけを表示しているのに対し、ロジスティック回帰は20,531遺伝子すべてと正解ラベルを使って分類境界を学習しているためである。

今回学んだこと

  • PCAは分類器ではなく、高次元データの構造を確認するための方法である
  • PCA上で分布が重なっても、教師あり学習で分類できる可能性がある
  • 評価用データの情報を漏らさないため、標準化はデータ分割後に学習する
  • Pipelineを使うと、交差検証の各foldでも前処理を学習データだけに適用できる
  • 1回のテスト結果だけでなく、交差検証で分割に対する安定性を確認する
  • 特徴量数がサンプル数より多くても、正則化と適切な評価により高い汎化性能が得られる場合がある

感想

今回は、PCAの散布図から癌種の分離を考えるところから始めた。当初は、点が重なっているため癌種推定は難しいのではないかと考えた。しかし、PCAは癌種ラベルを見ずに分散の大きい方向を求めていること、ロジスティック回帰は全遺伝子と正解ラベルを使うことを理解すると、可視化と分類結果の違いを整理できた。

また、20,531遺伝子という特徴量数を見て過学習を心配したが、特徴量数だけで判断するのではなく、訓練性能と検証性能の差、交差検証のばらつき、正則化の有無を確認する必要があると分かった。

次回

次回は、同じ交差検証条件で以下の2つを比較する。

  1. StandardScaler → ロジスティック回帰
  2. StandardScaler → PCA → ロジスティック回帰

PCAによって特徴量数と計算量を減らしたとき、分類性能をどの程度維持できるか確認する。

加えて、言語処理100本ノック2025の第1章「準備運動」に着手する。RNA-seqの実践課題と並行しながら、自然言語処理の基礎も進めていく。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?