近年、腸内フローラは、健康や病気に大きな影響を与えることが分かっており、その構成や機能を解析することで、新たな治療法や診断法の開発につながる可能性が示唆されています。
ここでは、腸管性下痢症(以下、EDD)患者と健康者の腸内マイクロバイオームデータを活用して、機械学習の手法を用いた予測モデルの構築例を見つけたので、再現してみました。
使用言語
R
データセットのダウンロード
以下のサイトのSession1をワーキングディレクトリの直下にダウンロードしてください。Session1ディレクトリ直下のinputディレクトリにメタデータとOTUテーブルが格納されています。
https://github.com/glbio-mlmb/MLMB_materials
Rパッケージのインストールと読み込み
#caretとstringrはCRANからインストール
install.packages("caret")
install.packages("stringr")
#randomForestはBioconductorからインストール
install.packages("BiocManager")
BiocManager::install("randomForest")
#pROCはGitHubからインストール
install.packages("devtools")
# または
install.packages("remotes")
devtools::install_github("xrobin/pROC")
# または
remotes::install_github("xrobin/pROC")
# パッケージの読み込み
library(caret)
library(randomForest)
library(pROC)
library(stringr)
予測モデルの構築と評価
以下で予測モデルの構築と評価を行っていきます。
参考元で実行されたコードの再現と各コードの詳細を記します。
現在のディレクトリパスの取得
current_dir <- dirname(as.character(rstudioapi::getSourceEditorContext()$path))
コードの詳細
rstudioapi::getSourceEditorContext()は、RStudioのソースエディターの状態に関する情報を返す。
$pathは、その情報の中からスクリプトのファイルパスだけを抽出するための演算
dirname()は、ファイルパスからディレクトリ名だけを取り出す関数
メタデータの読み込み
metadata <- read.table(paste0(current_dir,"input/edd_singh.metadata.txt"),sep="\t",header=T,row.names=1, stringsAsFactors=TRUE)
内容は、EDD患者と健康体のメタ情報、Stattusが患者か健康体かを示し、患者のPathogenが原因となったバクテリア（このデータは、16Sアンプリコンシークエンス解析なので）を示していると考えられる。
メタデータからのデータ抽出や編集
#メタデータは、304行 5列で構成されている。
dim(metadata)
[1] 304 5
# EDD患者の人数は、222人、健康体の数は82人
table(metadata$DiseaseState)
EDD H
222 82
# メタデータ中のtimepoint1のデータ数を見てみる
metadata <- metadata[metadata$Time.Point == 1, ]
dim(metadata)
283 5
# timepoint1のデータの内、EDD患者は201人 健康体は82人
table(metadata$DiseaseState)
EDD H
201 82
OTUテーブルの読み込み
otu <- read.table(paste0(current_dir,"input/edd_singh_otu_table"),sep="\t",header=T,row.names=1, stringsAsFactors=TRUE)
OTUテーブルの内容は以下の画像の通り、 数字が、各種微生物の存在量であると考えられる。
メタデータとOTUテーブルの共通しているサンプル数を見てみる。
common_samples <- intersect(rownames(metadata), colnames(otu))
length(common_samples)
282
コードの詳細
intersect()は、二つのオブジェクトの共通の値を取得する関数。
rownames()は、データフレームの行名を取得する関数。
colnames()は、データフレームの列名を取得する関数。
length()は、ベクトルやリストなどのオブジェクトの長さを取得する関数。
したがって、このコードは、metadataの行名とotuの列名の共通部分をintersect()で取得し、その長さをlength()で計算して、common_samplesという変数に代入するという処理を行っています。
出力結果から、common_samplesの長さは282であることがわかります。つまり、metadataとotuには282個の共通のサンプルがあるということです。
common_samplesの中身を確認してみる
otu <- otu[,common_samples]
dim(otu)
3698 282
metadata <- metadata[common_samples,]
dim(metadata)
282 5
all(rownames(metadata) == colnames(otu))
#TRUE
コードの詳細
このコードは、otuとmetadataという二つのデータフレームを、common_samplesというベクトルに含まれるサンプル名に基づいて部分的に抽出し、その次元を確認し、最後に行名と列名が一致しているかどうかを検査しています。具体的には、以下のような処理を行っています。
otu <- otu[,common_samples] この行では、otuというデータフレームの列を、common_samplesというベクトルに一致するものだけに絞り込んでいます。common_samplesには、otuとmetadataの両方に存在するサンプル名が入っていると思われます。この操作により、otuの列数はcommon_samplesの要素数と同じになります。
dim(otu) この行では、otuの次元を表示しています。[1] 3698 282という出力から、otuには3698行と282列があることがわかります。行はOTU（operational taxonomic unit）を、列はサンプルを表していると思われます。
metadata <- metadata[common_samples,] この行では、metadataというデータフレームの行を、common_samplesというベクトルに一致するものだけに絞り込んでいます。この操作により、metadataの行数はcommon_samplesの要素数と同じになります。
dim(metadata) この行では、metadataの次元を表示しています。[1] 282 5という出力から、metadataには282行と5列があることがわかります。行はサンプルを、列はメタデータを表していると思われます。
all(rownames(metadata) == colnames(otu)) この行では、metadataの行名とotuの列名がすべて一致しているかどうかを調べています。TRUEという出力から、両者は同じ順序で同じサンプル名を持っていることがわかります。これは、otuとmetadataを結合したり分析したりする際に重要な前提条件です。
データの前処理
OTUテーブルからのノイズ除去
#remove OTUs with fewer than 10 reads
## rowSums -- gives sum for each row of dataframe
otu <- otu[which(rowSums(otu)>=10),]
dim(otu)
1542 282
#remove OTUs which were present in fewer than 1% of samples
otu <- otu[which(rowSums(otu>0) >= ncol(otu)*.01),]
dim(otu)
1319 282
コードの詳細
このコードは、otuというデータフレームからノイズとなるOTU（operational taxonomic unit）を除去するためのものです。具体的には、以下のような処理を行っています。
otu <- otu[which(rowSums(otu)>=10),] この行では、otuの各行の合計値が10以上のものだけを抽出しています。rowSums ()は、データフレームの各行の和を計算する関数。
which ()は、論理ベクトルの中でTRUEとなる要素の位置を返す関数。この操作により、otuの行数は1542に減ります。
dim(otu) この行では、otuの次元を表示している。[1] 1542 282という出力から、otuには1542行と282列があることがわかります。行はOTUを、列はサンプルを表していると思われます。
otu <- otu[which(rowSums(otu>0) >= ncol(otu)*.01),] この行では、otuの各行の0より大きい値の数が、列数の1%以上のものだけを抽出しています。rowSums (otu>0)は、otuの各行の0より大きい値の個数を計算する関数です。
ncol ()は、データフレームの列数を返す関数。この操作により、otuの行数は1319に減ります。
dim(otu) この行では、otuの次元を再び表示。[1] 1319 282という出力から、otuには1319行と282列があることがわかります。
微生物の分類単位のデータの集約
## collapse OTUs to genus level by summing their respective abundance counts
## Split taxa names by ";" into 8 parts (taxonomic levels)
tax_table <- str_split_fixed(rownames(otu),";",n=8)
## Append genus name (6th column in genus) to the otu table
otu <- data.frame(genus=tax_table[,6],otu)
## Filter out otus where genus is not characterized.
otu <- otu[!(otu$genus=="g__"),]
dim(otu)
944 283
## summarize otu table by genus label
otu <- aggregate(otu[,-1], by=list(otu$genus), sum)
rownames(otu) <- otu[,1]
otu <- otu[,-1]
dim(otu)
122 282
上記のコードが実行された後のOTUテーブルの内容
1列目が、属名になっている。
コードの詳細
このコードは、OTU（operational taxonomic unit）と呼ばれる微生物の分類単位のデータを集約するためのものです。具体的には、以下のような処理を行っています。
str_split_fixed関数を使って、rownames(otu)に含まれるタクソン名（分類階級ごとに;で区切られた文字列）を8つの部分に分割し、tax_tableという行列に格納します。
data.frame関数を使って、tax_tableの6列目（属名を表すg__の後に続く部分）をgenusという列名でotuというデータフレームに追加します。
otuのgenus列で、g__となっているもの（属名が特定されていないもの）を除外します。
dim関数を使って、otuの行数と列数を確認します。この時点では、944行（OTU数）と283列（genus列と282個のサンプル列）になっています。
aggregate関数を使って、otuのgenus列でグループ化し、同じ属に属するOTUの各サンプルでのカウント数（豊富度）を合計します。otuのgenus列をby引数に、カウント数の列をx引数に、合計関数sumをFUN引数に指定します。
rownames関数を使って、otuの行名をgenus列の値に置き換えます。otuのgenus列は不要になるので、削除します。
再びdim関数を使って、otuの行数と列数を確認します。この時点では、122行（属数）と282列（サンプル数）になっています。
正規化
#calculate relative abundance of each genus by dividing its value by the total reads per sample
otu <- sweep(otu,2,colSums(otu),"/")
正規化後のOTUテーブル
コードの詳細
このコードは、otuデータフレームの各属の相対存在量を計算するためのものです。具体的には、以下のような処理を行っています。
colSums(otu)は、otuの各列（サンプル）の合計値（総リード数）を計算する関数。
sweep(otu,2,colSums(otu),"/")は、otuの各要素を、対応する列の合計値で割る関数。2は列方向を指定する引数。"/"は割り算を指定する引数。
otu <- sweep(otu,2,colSums(otu),"/")は、上記の関数の結果をotuという変数に代入する操作です。
これにより、otuの各要素は、その属のサンプルでの相対存在量（割合）に変換されます。
機械学習のための前準備
## Prepare for training
x <- data.matrix(otu)
## transpose to make rows are samples and feature (i.e. genera) as columns.
x <- t(x)
変数 x の内容は以下の通り、
コードの詳細
このコードは、otuデータフレームを機械学習のために準備するためのものです。具体的には、以下のような処理を行っています。
data.matrix(otu)は、otuデータフレームを数値行列に変換する関数。この関数は、データフレーム内のすべての変数を数値モードに変換し、それらを行列の列として結合します。因子や順序付き因子は、その内部コードに置き換えられます。
x <- data.matrix(otu)は、上記の関数の結果をxという変数に代入する操作です。この時点で、xは数値行列になります。行はOTU（operational taxonomic unit）を、列はサンプルを表しています。
t(x)は、xという行列の転置を計算する関数。この関数は、行列の行と列を入れ替えます。つまり、行列の要素x[i,j]は、転置後の要素x[j,i]になります。
x <- t(x)は、上記の関数の結果をxという変数に再代入する操作。この時点で、xは転置された数値行列になります。行はサンプルを、列はOTUを表しています。
DiseaseState列に含まれる因子のレベル確認
## relevel disease-state to make H as control
levels(metadata$DiseaseState)
"EDD" "H"
コードの詳細
このコードは、metadataというデータフレームのDiseaseStateという列に含まれる因子のレベルを確認するためのものです。1
levels関数2は、因子のレベルを取得または設定する関数。因子とは、カテゴリー変数を表すデータ構造で、レベルと呼ばれる値を持ちます。
metadata$DiseaseStateは、metadataというデータフレームのDiseaseStateという列を指定する方法。この列には、疾患の状態を表す因子が入っていると思われます。
levels(metadata$DiseaseState)は、metadataのDiseaseState列のレベルを取得する関数です。この関数の結果は、"EDD" "H"というベクトルになります。これは、DiseaseState列には、EDDとHという二つのレベルがあることを意味します。
訓練とテスト（Random Forest）
set.seed(1000)
# Split dataset into 80% training, and 20% test
train_index <- createDataPartition(metadata$DiseaseState, ## outcome
p = 0.8, ## percentage of training samples
list = FALSE ## show subsamples as matrix, not list
# times = 10 ## This will create 10 different 80% subsamples
)
View(train_index)
x.train <- x[train_index,]
y.train <- metadata$DiseaseState[train_index]
x.test <- x[-train_index,]
y.test <- metadata$DiseaseState[-train_index]
train_control <- trainControl(
method = "cv",
number = 5, ## also try 10
summaryFunction=twoClassSummary, # computes area under the ROC curve
classProbs = TRUE ## required for scoring models using ROC
)
set.seed(1000)
rf_train <- train( x = x.train, y = as.factor(y.train),
method='rf',
metric="ROC", ## default accuracy
trControl = train_control)
rf_train
# mtry ROC Sens Spec
# 2 0.9541896 0.95625 0.6659341
# 62 0.9469437 0.91250 0.7560440
# 122 0.9442995 0.91250 0.6659341
## mtry: Number of variables randomly sampled as candidates at each split.
rf_train$resample
# ROC Sens Spec Resample
# 1 0.9375000 0.96875 0.6923077 Fold1
# 2 0.9776786 0.96875 0.7142857 Fold2
# 3 0.9519231 0.93750 0.6923077 Fold5
# 4 0.9399038 0.93750 0.6153846 Fold4
# 5 0.9639423 0.96875 0.6153846 Fold3
## We can modify tuning grid using tuneGrid param in train().
rf_test <- predict(rf_train, x.test)
rf_test
[1] EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD H H
[21] EDD EDD EDD EDD EDD EDD EDD EDD EDD EDD H H H EDD EDD EDD EDD H EDD EDD
[41] H H EDD H H H H H H H EDD EDD H EDD EDD EDD
Levels: EDD H
# compare predicted outcome and true outcome
conf_matrix <- confusionMatrix(rf_test, y.test)
conf_matrix$table
# Reference
# Prediction EDD H
# EDD 39 5
# H 1 11
コードの詳細
このコードは、Rでランダムフォレストという機械学習の手法を用いて、二値分類（EDDかどうか）の問題を解くためのものです。ランダムフォレストとは、複数の決定木を組み合わせて予測を行うアルゴリズム。コードの内容は以下のように説明できます。
set.seed(1000) は、乱数の生成に用いるシード値を1000に設定することで、再現性を保つためのものです。
createDataPartition() は、caretパッケージの関数で、データを訓練用とテスト用に分割するためのもの。この例では、metadataというデータフレームのDiseaseStateという列を目的変数として、その値に応じてデータを層化抽出し、訓練用に80%、テスト用に20%の割合で分けています。list = FALSEとすることで、分割後のデータのインデックスを行列として返します。train_indexというオブジェクトには、訓練用データのインデックスが格納されます。
View(train_index) は、train_indexというオブジェクトの内容を表形式で表示するためのもの。
x.train <- x[train_index,] は、xデータフレームの訓練用データのインデックスに対応する行を抽出して、x.trainというオブジェクトに代入することで、説明変数の訓練用データを作成するためのものです。同様に、y.train <- metadata$DiseaseState[train_index] は、目的変数の訓練用データを作成するためのものです
。
x.test <- x[-train_index,] は、xデータフレームの訓練用データのインデックスに対応しない行を抽出して、x.testというオブジェクトに代入することで、説明変数のテスト用データを作成するためのものです。同様に、y.test <- metadata$DiseaseState[-train_index] は、目的変数のテスト用データを作成するためのものです。
train_control <- trainControl() は、caretパッケージの関数で、モデルの訓練に用いるパラメータを設定するためのもの。この例では、method = "cv"とすることで、交差検証という手法を用いて、モデルの性能を評価することを指定しています。交差検証とは、データをk個のグループに分けて、そのうちの一つをテスト用に、残りを訓練用に使ってモデルを構築し、テスト用データでモデルの精度を測定するということを、k回繰り返して、平均した精度を求める手法。 number = 5とすることで、交差検証のグループ数を5に設定しています。summaryFunction = twoClassSummaryとすることで、二値分類の問題に適した評価指標を計算する関数を指定して。 この関数は、ROC曲線下面積（AUC）、感度、特異度という指標を返します。classProbs = TRUEとすることで、各クラスに属する確率を計算することを指定しています。これは、ROC曲線下面積を求めるために必要です。
set.seed(1000) は、再び乱数のシード値を設定することで、再現性を保つためのものです。
rf_train <- train() は、caretパッケージの関数で、モデルの訓練を行うためのもの。この例では、x = x.train, y = as.factor(y.train)とすることで、説明変数と目的変数の訓練用データを指定しています。
as.factor()は、目的変数を因子型というデータ型に変換するためのものです。因子型とは、カテゴリー変数を表すデータ型。method = 'rf’とすることで、ランダムフォレストを用いることを指定しています。
metric = "ROC"とすることで、モデルの最適なパラメータを選択するための評価指標をROC曲線下面積に設定しています。
trControl = train_controlとすることで、先ほど設定したパラメータをモデルの訓練に用いることを指定しています。
rf_trainというオブジェクトには、訓練されたモデルが格納。
rf_train は、訓練されたモデルの内容を表示するためのもの。この例では、ランダムフォレストのパラメータの一つであるmtryという値を3通り（2, 62, 122）に変えて、それぞれのモデルのROC曲線下面積、感度、特異度を交差検証で求めています。mtryとは、各決定木を作る際に、ランダムに選択される説明変数の数を表すパラメータ。この例では、mtry = 2のときに最も高いROC曲線下面積（0.9541896）を得ているので、この値が最適なパラメータとして選択されます。
※mtryは、ランダムフォレストのパラメータの一つで、各分岐点でランダムに選択される変数の数を表します。この値が大きいほど、モデルの複雑さが増し、過学習のリスクが高まります。逆に小さいほど、モデルの汎化性能が向上しますが、バイアスが大きくなる可能性があります。適切な値を見つけるためには、交差検証やグリッドサーチなどの手法を用いることが一般的です。
rf_train$resampleは、train関数で作成したランダムフォレストのモデルの性能を評価するために、5分割交差検証を行った結果を表示。ROC曲線下の面積を表し、1に近いほど良いモデルと言えます。Sensは感度、Specは特異度と呼ばれる指標で、それぞれ正例を正しく判定する割合と負例を正しく判定する割合を表します。これらの値も1に近いほど良いモデルと言えます。Resampleは、交差検証で使われたデータの分割を表します。Fold1からFold5までの5つの分割があります。
tuneGridは、train関数の引数の一つで、モデルのパラメータを探索する範囲や方法を指定。例えば、mtryの値を1から10まで変化させて、最も良いモデルを選択する場合は、tuneGrid = data.frame(mtry = 1:10)というように設定します。この引数を省略すると、train関数はデフォルトの値や方法でパラメータを探索します。
rf_testは、predict関数で作成したランダムフォレストのモデルに、テストデータ（x.test）を入力して、予測したクラス（EDDかどうか）を表す。
conf_matrixは、confusionMatrix関数で作成した混同行列。混同行列は、予測したクラスと実際のクラスとの一致度を表す表です。表の行が予測したクラス、列が実際のクラスを表します。対角線上のセルは、予測と実際が一致した数を表し、それ以外のセルは、予測と実際が不一致した数を表します。 例えば、EDDと予測してEDDだったものは39個、EDDと予測してHだったものは5個、Hと予測してEDDだったものは1個、Hと予測してHだったものは11個です。混同行列からは、正解率や適合率、再現率などの指標を計算することができます。
モデルの評価
## Compute precision
conf_matrix$byClass
# Sensitivity Specificity Pos Pred Value Neg Pred Value Precision Recall
# 0.9750000 0.6875000 0.8863636 0.9166667 0.8863636 0.9750000
# F1 Prevalence Detection Rate Detection Prevalence Balanced Accuracy
# 0.9285714 0.7142857 0.6964286 0.7857143 0.8312500
## Can also spit out probability instead of predicted class
rf_test <- predict(rf_train, x.test, type = "prob")
rf_test
rf_test <- rf_test[,1]
re_testの結果
EDD H ER0004 0.714 0.286 ER0014 0.916 0.084 ER0017 0.846 0.154 ER0018 0.662 0.338 ER0019 0.870 0.130 ER0020 0.832 0.168 ER0022 0.976 0.024 ER0036 0.892 0.108 ER0046 0.712 0.288 ER0051 0.990 0.010 ER0058 0.942 0.058 ER0073 0.994 0.006 ER0078 0.668 0.332 ER0097 0.998 0.002 ER0099 0.912 0.088 ER0100 0.556 0.444 ER0111 0.650 0.350 ER0118 0.920 0.080 ER0123 0.454 0.546 ER0127 0.468 0.532 ER0134 0.880 0.120 ER0138 0.548 0.452 ER0140 0.916 0.084 ER0151 0.690 0.310 ER0157 0.954 0.046 ER0160 0.948 0.052 ER0164 0.892 0.108 ER0169 0.724 0.276 ER0171 0.972 0.028 ER0172 0.846 0.154 ER0177 0.314 0.686 ER0178 0.410 0.590 ER0191 0.436 0.564 ER0195 0.604 0.396 ER0199 0.962 0.038 ER0202 0.734 0.266 ER0219 0.546 0.454 ER0222 0.400 0.600 ER0232 0.718 0.282 ER0233 0.964 0.036 ER0240 0.240 0.760 ER0243 0.250 0.750 ER0267 0.534 0.466 ER0270 0.204 0.796 ER0271 0.168 0.832 ER0281 0.324 0.676 ER0303 0.474 0.526 ER0304 0.228 0.772 ER0323 0.342 0.658 ER0327 0.466 0.534 ER0348 0.996 0.004 ER0369 0.588 0.412 ER0384 0.496 0.504 ER0392 0.944 0.056 ER0410 0.786 0.214 ER0495 0.732 0.268コードの詳細
このコードは、ランダムフォレストを使って、EDDとHという二つのクラスに分類する問題を解いた後、そのモデルの性能を評価するためのものです。以下に詳しく説明します。
conf_matrix$byClassは、confusionMatrix関数で作成した混同行列のオブジェクトから、各クラスに対する評価指標を取り出すためのもの。この例では、Sensitivity（感度）、Specificity（特異度）、Pos Pred Value（陽性的中率）、Neg Pred Value（陰性的中率）、Precision（適合率）、Recall（再現率）、F1（F値）、Prevalence（母集団における陽性率）、Detection Rate（検出率）、Detection Prevalence（予測における陽性率）、Balanced Accuracy（平衡正解率）という指標が表示されます。
rf_test <- predict(rf_train, x.test, type = “prob”)は、predict関数で作成したランダムフォレストのモデルに、テストデータ（x.test）を入力して、各クラスに属する確率を出力するためのもの。type = "prob"とすることで、確率の出力を指定しています。rf_testというオブジェクトには、各サンプルに対するEDDとHの確率が格納されます。
rf_test <- rf_test[,1]は、rf_testというオブジェクトから、EDDの確率のみを取り出して、再代入するためのものです。この例では、EDDの確率がrf_testの1列目にあるので、rf_test[,1]としています。これにより、rf_testはEDDの確率のベクトルになります。
モデル評価(グラフの描画)
## ROC curve
rf <- roc(y.test,rf_test) ## pROC package
auc <- rf$auc
auc
# Area under the curve: 0.9633
## Plot ROC curve
plot(rf, col="blue",legacy.axes = TRUE)
コードの詳細
このコードは、モデルの性能を評価するためにROC曲線とAUCという指標を計算して描画するためのものです。以下に詳しく説明します。
ROC曲線とは、分類モデルの性能を表すグラフ。ROC曲線は、横軸に1-特異度（偽陽性率）、縦軸に感度（真陽性率）をとり、カットオフ値を変化させたときのモデルの挙動を示します。ROC曲線が左上隅に近いほど、モデルの性能が高いと言えます。
AUCとは、ROC曲線下面積の略で、ROC曲線の形状を一つの数値で表す指標。AUCは0から1の範囲の値をとり、1に近いほどモデルの性能が高いと言えます。
rf <- roc(y.test,rf_test)は、pROCパッケージのroc関数を使って、テストデータの真のクラス（y.test）とモデルの予測したEDDの確率（rf_test）から、ROC曲線のオブジェクト（rf）を作成するためのもの。
auc <- rf$aucは、ROC曲線のオブジェクト（rf）から、AUCの値（auc）を取り出すためのものです。
aucは、AUCの値（auc）を表示するためのものです。この例では、AUCは0.9633となっています。
plot(rf, col=“blue”,legacy.axes = TRUE)は、plot関数を使って、ROC曲線のオブジェクト（rf）をグラフに描画するためのものです。colは、曲線の色を指定する引数で、この例では青色（blue）にしています。legacy.axesは、x軸の目盛りを1-特異度に合わせて左が0、右が1となるように表示する引数で、この例ではTRUEにしています。
plotの結果
出典
コードおよびデータの提供元
https://github.com/glbio-mlmb/MLMB_materials
コードの元となった論文
https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2019.00579/full
モデル構築に活用したデータの提供元
https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2019.00579/full#B56