血中・尿中成分によって膵がん罹患を予測する機械学習モデルの構築

Last updated at 2024-08-01Posted at 2024-07-22

1. 本記事の趣旨

2024年2月4日から2024年8月5日までの期間でSAMURAI TERAKOYAのデータサイエンスオーダーメイドコースにおいてデータサイエンス、特に機械学習のうち、数値データに基づいた予測モデル構築を中心に学習した。この記事では、最終課題として選択した「Urinary biomarkers for pancreatic cancer」への取り組みを示す。

2. 課題の概要

膵がんは全てのがん種の中で最も予後の悪いがんであり、5年生存率は10%を下回る（最新がん統計）。その予後の悪さの原因のひとつに「発見の難しさ」が挙げられる。

Kaggleで公開されているデータセット「Urinary biomarkers for pancreatic cancer」の目的は、血液もしくは尿に含まれる下記のタンパク質群の濃度に基づいて、膵がん患者と非膵がん患者を区別するモデルを構築することである。

血中タンパク質

CA19-9

尿中タンパク質

クレアチニン
LYVE1
REG1B
TFF1
REG1A

このデータセットはKaggleでのホストが研究者から譲渡されたものを公開しているようだ。

研究者らはこのデータを用いて過去に論文を発表している。

A combination of urinary biomarker panel and PancRISK score for earlier detection of pancreatic cancer: A case–control study

SAMURAI TERAKOYAの最終課題として、これまでの学習に基づき、本データを使用して、膵がん患者と非膵がん患者を区別するモデルを構築することを目的とした。

3. データセットの紹介

3.1. 被験者数

分類	例数
健常人	183名
膵炎患者	208名
膵がん患者	199名
合計	590名

3.2. 特徴量

特徴量名	特徴量の意味	特徴量の備考
sample_id	検体ID
patient_cohort	データが取得されたタイミング	Cohort1、Cohort2
sample_origin	データが取得された施設	BPTB、LIV、ESP、UCL
age	年齢	連続変数
sex	性別	M＝男性、F＝女性
diagnosis	診断	1＝健常人、2＝膵炎患者、3＝膵がん患者
stage	膵がんのステージ	I、IA、IB、II、IIA、IIB、III、IV
benign_sample_diagnosis	膵炎疾患の種類	種類多数のため割愛
plasma_CA19_9	血中CA19-9濃度	連続変数
creatinine	尿中クレアチニン濃度	連続変数
LYVE1	尿中LYVE1濃度	連続変数
REG1B	尿中REG1B濃度	連続変数
TFF1	尿中TFF1濃度	連続変数
REG1A	尿中REG1A濃度	連続変数

3.3. ダウンロード可能なデータファイル

当データセットには下記2つのファイルが公開されている。

Debernardi et al 2020 data.csv
解析用ファイル
Debernardi et al 2020 documentation.csv
特徴量の説明ファイル

4. モデル構築環境

Google Colaboratory

5. モデル構築

5.1. 使用したライブラリ

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import lightgbm as lgb
import optuna
from sklearn.model_selection import train_test_split, 
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import classification_report, confusion_matrix, ConfusionMatrixDisplay, roc_curve, auc

5.2. データの読み込み

上記2つのデータセットをローカルにダウンロードした。ファイル名にスペースが含まれているが、スペースはエラーの原因になる可能性があることから、スペースをアンダースコア(_)に書き換えた。これらのファイルをGoogle driveにアップロードした後、データを読み込んだ。

まず、解析用ファイルを示す。

df = pd.read_csv('/content/drive/MyDrive/pdac_urine/Debernardi_2020_data.csv')
df.head()

つぎに、特徴量の説明ファイルを示す。

doc = pd.read_csv('/content/drive/MyDrive/pdac_urine/Debernardi_2020_doc.csv')
doc

5.3. データの概要の確認

df.shape

(590, 14)

Kaggleサイトでの説明通り590名分のデータを含まれていた。

df.info()

「stage」「benign_sample_diagnosis」「plasma_CA19_9」「REG1A」は欠損値を含むようだ。

欠損値を確認する。

missing_values = df.isnull().sum()
missing_features = missing_values[missing_values > 0]
print(missing_features)

df.describe()

「plasma_CA19_9」、「LYVE1」、「REG1B」、「TFF1」、「REG1A」について、第3四分位数と最大値との間に大きな差が見られる。このことは外れ値が存在することを示唆している。

5.4. 学習用データとテスト用データへの分割

ダウンロード可能なデータにはテストデータは含まれていないため、解析用ファイルの一部をテスト用データとして切り分けておく。

当データセットはデータ数が少ない。学習用データにデータの8割を使用するために、test_size = 0.2とした。

また、目的変数であるdiagnosisが学習用データとテスト用データへ均等に振り分けられるように引数にstratifyパラメータを設定した。

train_df, test_df = train_test_split(df, test_size=0.2, stratify=df['diagnosis'], random_state=42)

print('train_df:', train_df.shape)
print('test_df:', test_df.shape)

train_df: (472, 14)
test_df: (118, 14)

5.5. データの前処理

5.5.1. 特徴量の削除

まず、「sample_id」はモデル構築に使用しないので、削除する。

次に、本データを用いた予測モデル構築について、下記の注意書きが記されている。

原文：
The dataset includes information on stage of pancreatic cancer, and diagnosis for non-cancerous patients, but remember—these won't be available to a predictive model. The goal, after all, is to predict the presence of disease before it's diagnosed, not after!

日本語訳：
このデータセットには、膵がんのステージおよび非がん患者の診断に関する情報が含まれているが、これらの情報は予測モデルには利用できない。なぜならば、最終的な目標は、診断後ではなく、診断前に疾患の存在を予測することであるからだ。

この記載に沿って、「stage」と「benign_sample_diagnosis」をデータフレームから削除する。

データセットには「patient_cohort」という特徴量が含まれる。これは「データがいつ取られたか」に相当するが、予測結果はデータが取られたタイミングに依存してはならないため、「patient_cohort」も削除する。

データセットには「sample_origin」という特徴量が含まれる。これは「データがどの施設で取られたか」に相当するが、予測結果はデータが取られた施設に依存してはならないため、「sample_origin」も削除する。

columns_to_drop = ['sample_id', 'stage', 'benign_sample_diagnosis', 'patient_cohort', 'sample_origin']

train_df = train_df.drop(columns=columns_to_drop)
test_df = test_df.drop(columns=columns_to_drop)

print('train_df:', train_df.shape)
print('test_df:', test_df.shape)

train_df: (472, 9)
test_df: (118, 9)

5.5.2. 「diagnosis」の変換

本データを用いた予測モデル構築について、下記の注意書きが記されている。

原文：
The goal in this dataset is predicting diagnosis, and more specifically, differentiating between 3 (pancreatic cancer) versus 2 (non-cancerous pancreas condition) and 1 (healthy).

日本語訳：
このデータセットの目標は診断の予測であり、具体的には3（膵がん）と2（非がん性膵臓疾患）および1（健康）を区別することである。

この記載に沿って、「diagnosis」の「1」と「2」を「control」とし、「3」を「pdac」とする。（pdac= pancreatic ductal carcinoma）

# 置換のためのディクショナリー
replacement_dict = {1: 'control', 2: 'control', 3: 'pdac'}

# train_dfの「diagnosis」カラムを置換
train_df['diagnosis'] = train_df['diagnosis'].replace(replacement_dict)

# test_dfの「diagnosis」カラムを置換
test_df['diagnosis'] = test_df['diagnosis'].replace(replacement_dict)

5.5.3. ラベルエンコーディング

データフレームに含まれるカテゴリー変数（sexおよびdiagnosis）に対してラベルエンコーディングを実施する。

# データ型がobjectのカラムを取得
object_cols_train = train_df.select_dtypes(include='object').columns
object_cols_test = test_df.select_dtypes(include='object').columns

# ラベルエンコーダーのインスタンスを作成
label_encoder = LabelEncoder()

# トレーニングデータのラベルエンコーディング
for col in object_cols_train:
    train_df[col] = label_encoder.fit_transform(train_df[col])

# テストデータのラベルエンコーディング
for col in object_cols_test:
    test_df[col] = label_encoder.fit_transform(test_df[col])

5.5.4. データの確認

前処理を実施したデータの一部を表示する。

train_df.head()

test_df.head()

5.5.5. 外れ値の処理

train_dfを対象に、連続変数の特徴量についてstripplotを作成する。

# 特徴量の指定
columns = ['age', 'sex', 'plasma_CA19_9', 'creatinine', 'LYVE1', 'REG1B', 'TFF1', 'REG1A']

# 図と軸を作成
fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(10, 6))  # 4列2行のグリッドを作成
axes = axes.flatten()  # 2次元配列を1次元配列に変換

# 連続変数の列ごとにストリッププロットを描画
for i, column in enumerate(columns):
    sns.stripplot(x='diagnosis', y=column, data=train_df, jitter=True, ax=axes[i])
    axes[i].set_xlabel('Diagnosis')
    axes[i].set_ylabel(column)

plt.tight_layout()
plt.show()

「plasma_CA19_9」、「TFF1」、「REG1A」について、他と大きく離れた値が含まれている。

続いて、test_dfを対象に、連続変数の特徴量について、同様のコードでstripplotを作成する。

「REG1A」について、他と大きく離れた値が含まれている。

5.5.5.1. 「plasma_CA19_9」の解釈

「plasma_CA19_9」は血漿中CA19-9濃度を指す。CA19-9は主にがん細胞が産生する糖鎖である。この検査値は膵がんのマーカーのひとつで、この数値の上昇は膵がんの罹患や悪性進展に関連している。

この検査値の正常値は0～37（単位：U/mL）であるため、train_df中の30000を超える値は、正常値と比較して高い。過去に、50000以上のCA19-9を示す患者が報告されている（人間ドック (Ningen Dock), 2015, 30 巻, 1 号, p. 22-29）。

人間ドック (Ningen Dock), 2015, 30 巻, 1 号, p. 22-29より転載

今回のデータセットのデータ数は590と少数である。データ数を10倍、50倍した時に、一定の割合で万単位のCA19-9を示す患者が出現することが予想される。

したがって、本データセットに含まれる30000を超えるCA19-9は欠損値として扱わないこととした。

5.5.5.2. 「TFF1」や「REG1A」の解釈

これらは尿中のTFF1ならびにREG1Aタンパク質の濃度を指す。しかし、これらの因子が膵がん患者の尿中において10000（単位：pg/mL）前後の値を示す例を見つけることはできなかった。CA19-9と同様にがん細胞が産生するタンパク質であり、高レベルの値を示す可能性は十分にあることから、これら2つの因子についても外れ値として扱わないこととした。

5.5.5.3. 特徴量の相関性

特徴量間の相関性を検討するために、散布図を作成した。

%matplotlib inline

# カスタムパレットの作成
custom_palette = {0: 'blue', 1: 'red'}

sns.pairplot(train_df, hue='diagnosis', palette=custom_palette)

続いて、相関係数ヒートマップを作成した。

# 相関係数を計算
correlation_matrix = train_df.corr()

# ヒートマップを作成
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Matrix Heatmap')
plt.show()

下記の特徴量間の相関係数が高値（0.7以上）を示した。

組み合わせ	相関係数
TFF1 - REG1B	0.71

散布図と相関係数ヒートマップから、TFF1とREG1Bとの間に正の相関が認められた。

これらのうちどちらかを分析データから除外することを検討する必要があるかもしれない。

5.6. 適した予測モデル構築アルゴリズムの探索

欠損値や標準化、正規化などの前処理は実施していないが、ここで一旦、pycaretを使用して適した予測モデル構築アルゴリズムを探索した。

# ライブラリをインストール
pip install pycaret

# ライブラリのインポート
import pandas as pd
from pycaret.classification import *

# PyCaretのセットアップ
# ターゲットカラム名：'diagnosis'
clf = setup(train_df, target='diagnosis')

# モデルの比較
best_model = compare_models()

# ベストモデルの作成
model = create_model(best_model)

# モデルの評価
evaluate_model(model)

# モデルの保存
save_model(model, 'best_model')

Gradient Boosting ClassifierやLightGBM、XGBoostといった勾配ブースティング決定木系のアルゴリズムが上位を占めた。

上記3つのアルゴリズムのうち、Gradient Boosting Classifierは欠損値の処理が必要なアルゴリズムであるため、欠損値の処理方法が予測精度に影響を与える可能性がある。一方、lightGBMやXGBoostは欠損値を自動的に処理する機能を備えている。

そこで、予測モデル構築アルゴリズムとして、Gradient Boosting Classifierの次にF1 scoreが高値を示したlightGBMを選択した。

5.7. 学習用データを訓練用データと評価用データに分割

# 特徴量とターゲット変数を定義
X = train_df.drop(columns=['diagnosis'])
y = train_df['diagnosis']

# データを学習用と評価用に分割
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, stratify=df['diagnosis'], random_state=42)

# X_train, X_val, y_train, y_valの行列数を表示
print(X_train.shape, X_val.shape, y_train.shape, y_val.shape)

(377, 8) (95, 8) (377,) (95,)

5.8. LightGBMによる予測モデル構築

まず、ハイパーパラメーターを設定せずにモデルを構築する。

# LightGBM用データセットの作成
train_data = lgb.Dataset(X_train, label=y_train)
val_data = lgb.Dataset(X_val, label=y_val, reference=train_data)

# ハイパーパラメータの設定
params = {
    'boosting_type': 'gbdt',
    'objective': 'binary',
    'metric': ['binary_logloss', 'auc'],
}

# モデルの訓練
num_round = 100
bst = lgb.train(params, train_data, num_round, valid_sets=[val_data])

# テストデータでの予測
y_pred = bst.predict(X_val, num_iteration=bst.best_iteration)

# 二値分類なので、0.5を閾値にしてクラスを割り当てる
y_pred_binary = [1 if x >= 0.5 else 0 for x in y_pred]

Classification reportを作成する。

report = classification_report(y_val, y_pred_binary)
print(report)

非がん患者に対するF1 score = 0.91で、膵がん患者に対するF1 score = 0.79を示した。

Confusion matrixを作成する。

# Confusion Matrixを作成
cm = confusion_matrix(y_val, y_pred_binary)

# Confusion Matrixをプロット
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap=plt.cm.Blues)
plt.title("Confusion Matrix")
plt.show()

非がん患者66名のうち5名（7.58%）を膵がんと判断し、膵がん患者29名中7名（24.1%）を健康と判断している。

ROC曲線を作成する。

# ROC曲線を作成
fpr, tpr, _ = roc_curve(y_val, y_pred_binary)
roc_auc = auc(fpr, tpr)

# ROC曲線をプロット
plt.figure()
plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve')
plt.legend(loc="lower right")
plt.show()

ROC曲線のAUC = 0.84を示した。

以上の結果から、モデルについて改善の余地があると考えられた。

5.9. Optunaによるハイパーパラメーターチューニング

import lightgbm as lgb
import optuna
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split

def objective(trial):
    # ハイパーパラメータの設定
    param = {
        'objective': 'binary',
        'metric': 'binary_logloss',
        'boosting_type': 'gbdt',
        'num_leaves': trial.suggest_int('num_leaves', 20, 50),
        'learning_rate': trial.suggest_loguniform('learning_rate', 0.01, 0.1),
        'n_estimators': trial.suggest_int('n_estimators', 100, 300),
        'max_depth': trial.suggest_int('max_depth', 3, 8),
        'min_child_samples': trial.suggest_int('min_child_samples', 20, 50),
        'subsample': trial.suggest_uniform('subsample', 0.5, 1.0),
        'colsample_bytree': trial.suggest_uniform('colsample_bytree', 0.5, 1.0)
    }

    # モデルの学習
    gbm = lgb.LGBMClassifier(**param)
    gbm.fit(X_train, y_train, eval_set=[(X_val, y_val)])
    
    # モデルの予測
    y_pred = gbm.predict(X_val)
    
    # 評価指標としてaccuracyを使用
    accuracy = accuracy_score(y_val, y_pred)
    
    return accuracy

# Optunaのstudyを作成し、最適化を実行
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=100)

# 最適なハイパーパラメータの出力
print(f"Best trial parameters: {study.best_trial.params}")

Best trial parameters: {'num_leaves': 38, 'learning_rate': 0.07920690662206804, 'n_estimators': 173, 'max_depth': 7, 'min_child_samples': 39, 'subsample': 0.7399343402768092, 'colsample_bytree': 0.8944502338643244}

5.10. チューニングしたハイパーパラメーターによる予測モデルの再構築

これらのハイパーパラメーターを加えて、再度予測モデルを構築する。

# ハイパーパラメータの設定
params = {
    'boosting_type': 'gbdt',
    'objective': 'binary',
    'metric': ['binary_logloss', 'auc'], 
    'num_leaves': 38, 
    'learning_rate': 0.07920690662206804, 
    'n_estimators': 173, 
    'max_depth': 7, 
    'min_child_samples': 39, 
    'subsample': 0.7399343402768092, 
    'colsample_bytree': 0.8944502338643244
}

# モデルの訓練
num_round = 100
bst = lgb.train(params, train_data, num_round, valid_sets=[val_data])

# テストデータでの予測
y_pred = bst.predict(X_val, num_iteration=bst.best_iteration)

# 二値分類なので、0.5を閾値にしてクラスを割り当てる
y_pred_binary = [1 if x >= 0.5 else 0 for x in y_pred]

Classification report、Confusion matrix、ROC曲線は以下の通りである。

ハイパーパラメーターチューニング実施前後のF1 scoreは下記のとおりである。

	非がん患者	膵がん患者
前	0.91	0.79
後	0.93	0.84

デフォルトのLightGBMと比較して、ハイパーパラメーターチューニングしたLightGBMでは予測精度の改善が見られた。

5.11. 特徴量の重要度解析

特徴量の重要度評価に一般的に使用される「分割に使用された回数（split）」と「情報ゲイン（gain）」について解析した。

lgb.plot_importance(bst, max_num_features=10, importance_type='split')
plt.title("Feature Importance (split)")
plt.show()

lgb.plot_importance(bst, max_num_features=10, importance_type='gain')
plt.title("Feature Importance (gain)")
plt.show()

これらの結果は、特徴量の中で「plasma_CA19_9」と「LYVE1」が膵がん患者と非膵がん患者を区別する重要な因子であることを示唆している。

5.12. 特徴量の選択

「5.5.5.3. 特徴量の相関性」において、TFF1とREG1Bとの間に強い正の相関が見られた。

多重共線性の観点から、TFF1とREG1Bのどちらかをモデル構築から除外することを検討する。

今回モデル構築アルゴリズムとしてLightGBMを選択した。このアルゴリズムは決定木系であるため多重共線性の影響をうけにくいと考えられているが、特徴量除外が今回構築した予測モデルの精度に与える影響を検討した。

5.12.1. TFF1を除外した場合

Classification report、Confusion matrix、ROC曲線は以下の通りである。

膵がん患者に対するF1 scoreが0.84から0.82に低下した。

5.12.2. REG1Bを除外した場合

Classification report、Confusion matrix、ROC曲線は以下の通りである。

膵がん患者に対するF1 scoreが0.84から0.81に低下し、ROC曲線のAUCが0.87から0.86に低下した。

これらの結果から、TFF1とREG1Bの双方をモデル構築の特徴量として使用することとした。

6. テストデータを用いた予測

これまでに構築したモデルを使用して、冒頭で切り分けておいたテスト用データを対象に予測を実施する。

# test_dfのdiagnosisをy_testに、それ以外をX_testに分割
X_test = test_df.drop(columns=['diagnosis'])
y_test = test_df['diagnosis']
print(X_test.shape, y_test.shape)

(118, 8) (118,)

# テストデータでの予測
y_pred_test = bst.predict(X_test, num_iteration=bst.best_iteration)

# 二値分類なので、0.5を閾値にしてクラスを割り当てる
y_pred_test_binary = [1 if x >= 0.5 else 0 for x in y_pred_test]

Classification report、Confusion matrix、ROC曲線は以下の通りである。

非膵がん患者を対象としたF1 score = 0.94、膵がん患者を対象としたF1 score = 0.89を示し、Accuracy score = 0.92を示した。

非膵がん患者78名中74名（94.8%）を非膵がんと判別し、膵がん患者40名中35名（87.5%）を膵がんと判断した。

ROC曲線のAUC = 0.91を示した。

7. 結果の要約

今回、年齢や性別、血漿中CA19-9、尿中Creatinine、LYVE1、REG1B、TFF1、REG1Aの8つの特徴量から膵がんと非膵がんを区別するモデルの構築に取り組み、LightGBMをアルゴリズムとして選択した予測モデルを構築した。
非膵がんを約95％、膵がんを約88%の確率での正解が期待できる予測モデルを構築した。
8つの特徴量の中で、特に血漿中CA19-9と尿中LYVE1が膵がんと非膵がんを区別する重要な因子であると考えられた。

8. 考察

血漿中CA19-9は膵がんマーカーとして臨床使用されている。今回の検討からもその重要性が再確認された。本データセットの血漿中CA19-9は多くの欠損値を含んでいた。データが提供されていたら、どのようなモデルが構築されていただろうか。尿中成分も加味した膵がん予測においても血漿中CA19-9データはしっかり収集すべきである。
血漿中CA19-9と共に尿中LYVE1が重要な特徴量であることが示唆された。LYVE1（lymphatic vessel endothelial hyaluronan receptor 1）は、322アミノ酸残基からなるI型膜内在性糖タンパク質であり、主にリンパ管内皮細胞に発現している。これまでにLYVE1が膵がんの発症や悪性進展に関与するか否かは報告されていない。今後の機能解明が期待される。
尿中REG1BとTFF1の組み合わせについて、発現レベルの相関性が観察された。これらの病態生理学的意義の解明は、膵がんの発症機序解明につながるかもしれない。

9. モデル構築の限界

データ数が590と少なくことが予測精度に影響を与えている可能性がある。より多症例での検討が要求される。

10. データサイエンスオーダーメイドコースを終えて

学習を開始した6か月前は「データサイエンス」や「機械学習」「Python」などの単語は知っていたが、そのつながりや実際の解析方法については無知であった。

学習を終えた今、データサイエンスのレベルはまだまだ未熟ではあるが、データサイエンス全体のイメージが湧き、解析したいことをコードとして表現し、可視化できるようになった。

SAMURAI TERAKOYAの教材はわかりやすい表現で記載されており、またステップバイステップで画像解説を付いており、スムーズに学習を進めることができた。

成長できた最大の要因は素晴らしい講師に出会えたことである。先生は現在フリーのデータサイエンティストとして活躍なさっており、機械学習の基礎はもちろんのこと、実務での解析手順など教材には書かれていないことを多く教えてくださった。今後は、より実務を意識して学習を進めていきたい。

私は現在、基礎医学研究に従事していることもあり、受講期間は特に数値データに基づいた予測モデル構築を学習したが、データサイエンスのニーズはこの分野に留まらない。今後は画像認識や物体検出、自然言語処理も学習していきたい。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up