【第5回】ゼロから始めるゲノム解析（R編）

Transcript

1. Predictive Modeling with Supervised Machine Learning Naoki Kimoto (@_kimoton) 【第５回】ゼロから始めるゲノム解析

2. アンケート結果 + α

3. • クロマチン構造と発現量との関係性を研究 するために、新規の定量的モデルを考案し た。 • 測定された遺伝子発現量は、開発したモデ ルによってPearsonの相関係数0.79～0.88 の精度で予測された。 • 遺伝子の90.44％を「オン/オフ」のカテゴリ

   に正しく分類した。このモデルの精度は、 AUC=0.95だった。 • H3K9acやH3K4me3などのヒストン修飾は、 遺伝子のオン/オフを特定するために重要 だが、H3K79me2およびH3K36me3などの ヒトン修飾は、遺伝子量の予測にとってより 重要。 応用例①エピジェネティックな修飾から遺伝子発現を予測 Dong, Greven, Kundaje, et al. 2012. “Modeling gene expression using chromatin features in various cellular contexts.” Genome Biol. 13 (9): R53.

4. • 特定のDNA調節領域でのヒストン修飾は、 遺伝子の活性化、不活性化、およびポージ ングに関連している。 • ChromaGenSVMと呼ばれる遺伝的アルゴリ ズムの最適化と組み合わせたサポートベク ターマシンを使用したクロマチン状態検出 の手法を開発した。 •

   ChromaGenSVMは、ENCODE projectで 実験的に確認されたHeLa細胞のエンハン サーの88％を予測した。 • エンハンサーがヒストンのメチル化とアセチ ル化の特定の組み合わせによって最もよく 予測されることを示した。 応用例②エンハンサーまたは他の調節領域の予測 Fernandez, and Miranda-Saavedra. 2012. “Genome-wide enhancer prediction from epigenetic signatures using genetic algorithm-optimized support vector machines.” Nucleic Acids Res. 40 (10): e77.

5. 機械学習は、大別すると教師あり学習と教師なし学習（＋強化学習） に分けられる。 教師あり学習 機械学習 教師あり学習 教師なし学習 分類問題 回帰問題 クラスタリング 入力データの正解ラベルなど、

   出力に相当するデータがある場合 入力データのみに基づいて データの特徴を学習 本日扱う範囲

6. 機械学習 vs 統計 機械学習 統計 機械学習と統計はいずれもデータから示唆を見出すという点で共通している。 最適化と パフォーマンスを 重視 統計的特性の発見と

   不確実性の推定 を重視 データを学習し 示唆を出す

7. どのような手法を使う場合でも、教師あり学習の基本的な流れは下記のようになる。Rで は、caretというパッケージを使用することでこの一連の流れを扱うことができる（同様の パッケージとして、mlr等もある）。 教師あり学習におけるステップ データの前処理 外れ値の除去・変 数選択・正規化・ データ変換・欠損 値の処理 データの分割

   学習データ・評価 データに分割 モデルの学習 データセットを用 いて損失関数を 最小化 モデルの評価 モデルの パフォーマンスを 評価 モデルの選択 様々なパラメータ を試し、最適なモ デルを選択

8. The Cancer Genome Atlasの神経膠芽腫腫瘍サンプルの遺伝子発現データを用いて、 CIMPサブタイプを持っているかを予測するモデルを構築 Q. ゲノムデータから疾病タイプの予測 ゲノムシーケンスとバイオインフォマティクスを使用し て、がんの原因となる遺伝子変異をカタログ化するプ ロジェクトで、2005年に開始された。

   WHO グレード４の悪性腫瘍。 脳の神経細胞を支える神経膠細胞が腫瘍化 したもの。 CpG Island Methylator Phenotypeの略。 遺伝子の転写を抑制するプロモーター領域の癌特異的 CpGアイランド高メチル化

9. The Cancer Genome Atlasの神経膠芽腫腫瘍サンプルの遺伝子発現データを用いて、 CIMPサブタイプを持っているかを予測するモデルを構築 入力データ # 遺伝子発現量のファイルを読み込む fileLGGexp=system.file("extdata", "LGGrnaseq.rds"

   , package="compGenomRData" ) # 遺伝子発現量 gexp=readRDS(fileLGGexp) head(gexp[,1:5]) ## TCGA-CS-4941 TCGA-CS-4944 TCGA-CS-5393 TCGA-CS-5394 TCGA-CS-5395 ## A1BG 72.2326 24.7132 46.3789 37.9659 19.5162 ## A1CF 0.0000 0.0000 0.0000 0.0000 0.0000 ## A2BP1 524.4997 105.4092 323.5828 19.7390 299.5375 ## A2LD1 144.0856 18.0154 29.0942 7.5945 202.1231 ## A2ML1 521.3941 159.3746 164.6157 63.5664 953.4106 ## A2M 17944.7205 10894.9590 16480.1130 9217.7919 10801.8461 # 患者のアノテーション情報のファイルを読み込む fileLGGann=system.file("extdata", "patient2LGGsubtypes.rds" , package="compGenomRData" ) # 患者のアノテーション情報 patient=readRDS(fileLGGann) head(patient) ## subtype ## TCGA-FG-8185 CIMP ## TCGA-DB-5276 CIMP ## TCGA-P5-A77X CIMP ## TCGA-IK-8125 CIMP ## TCGA-DU-A5TR CIMP ## TCGA-E1-5311 CIMP

10. 実験における測定誤差やデータ取得の仮定で生じる潜在的な問題により、データスケー ルに違いが生じていないか箱ひげ図により確認する。 理想的には、各腫瘍サンプルは遺伝子発現値の同様の分布を持っている。 データ前処理（正規化） boxplot(gexp[,1:50],outline=FALSE,col="cornflowerblue" ) サンプルあたりの遺伝子発現値は同じス ケールとなっており、正規化が既に行われ たと考えられる 最初の50個の腫瘍

    サンプルを確認 一方で、裾が長い分布に見えるため、 対数変換により扱いやすい分布に変換でき ると考えられる（次項）

11. 対数変換により扱いやすい分布に変換する。 データ前処理（変換） par(mfrow=c(1,2)) # 5番目の患者について、遺伝子発現量の分布を可視化する hist(gexp[,5],xlab="gene expression" ,main="", border="blue4", col="cornflowerblue"

    ) # 5番目の患者について、対数変換した遺伝子発現量の分布を可視化する hist(log10(gexp+1)[,5], xlab="gene expression log scale" , main="", border="blue4",col="cornflowerblue" ) 対数変換 gexp=log10(gexp+1) ※発現量が0の場合でも扱えるように 1を足す

12. 変動の少ない説明変数は目的変数を予測する上で重要でないと考えられる。 計算速度を低下させる要因となるため、前処理で除去することを検討する。 データ前処理（外れ値の除去） # データセットを転置する tgexp <- t(gexp) library(caret) #

    85% の値が等しい分散がほぼゼロの列を削除する # この関数はフィルターを作成するのみで、適用はまだ nzv=preProcess(tgexp,method="nzv",uniqueCut = 15) # "predict" 関数を用いてフィルターを適用する # nzv_tgexp 変数にフィルタリング後のデータフレームを代入する nzv_tgexp=predict(nzv,tgexp)

13. 変動の大きさについて、任意のカットオフ値を選択することもできる。 下記では、上位1000個の変数予測子を選択する。 データ前処理（外れ値の除去） par(mfrow=c(1,2)) # 5番目の患者について、フィルタリング前の遺伝子発現量の分布を可視化する hist(tgexp[,5],xlab="gene expression" ,main="", border="blue4",

    col="cornflowerblue" ) # 発現量上位1000個の遺伝子にフィルタリング SDs=apply(tgexp,2,sd ) topPreds=order(SDs,decreasing = TRUE)[1:1000] tgexp=tgexp[,topPreds] # 5番目の患者について、フィルタリング後の遺伝子発現量の分布を可視化する hist(tgexp[,5],xlab="gene expression" ,main="", border="blue4", col="cornflowerblue" ) 外れ値除去

14. 変動の大きさについて、任意のカットオフ値を選択することもできる。 標準化により、全ての予測変数の平均は0、標準偏差は1になる。 この処理は通常、一部の計算の数値安定性を向上させるために使用される。 データ前処理（標準化） par(mfrow=c(1,2)) # 5番目の患者について、標準化前の遺伝子発現量の分布を可視化する hist(tgexp[,5],xlab="gene expression" ,main="",

    border="blue4", col="cornflowerblue" ) # 遺伝子発現量を標準化 library(caret) processCenter= preProcess(tgexp, method = c("center")) tgexp=predict(processCenter,tgexp) # 5番目の患者について、標準化後の遺伝子発現量の分布を可視化する hist(tgexp[,5],xlab="gene expression" ,main="", border="blue4", col="cornflowerblue" ) 標準化

15. 一部の手法では予測変数間の相関を適切に処理できるため、外れ値の除去を行わなく ても問題ない。一方で予測変数が少ないほど、モデルの学習を高速化することができ る。 データ前処理（相関の高い変数の除去） # 高い相関を持つ変数を除去するためのフィルターを作成する # 2つの変数が相関を持つときいずれか 1つのみが除去される corrFilt=preProcess(tgexp,

    method = "corr",cutoff = 0.9) tgexp=predict(corrFilt,tgexp) M.before <- cor(tgexp) highlyCorDescr <- findCorrelation(M.before , cutoff = 0.9) tgexp <- tgexp[,-highlyCorDescr] M.after <- cor(tgexp)

16. 実験中の技術的な問題により、特定の遺伝子やゲノム位置の値がない場合（欠落値） が起こり得る。中央値や平均値、最近傍で補完することが一般的。 データ前処理（欠損値除去） # 欠損値が存在するデータを生成 missing_tgexp=tgexp missing_tgexp[ 1,1]=NA anyNA(missing_tgexp) #

    欠損値が存在するか確認 ## [1] TRUE # 欠損値を抽出 gexpnoNA=missing_tgexp[ , colSums(is.na(missing_tgexp)) == 0] # 中央値で補完 library(caret) mImpute=preProcess(missing_tgexp, method="medianImpute" ) imputedGexp=predict(mImpute,missing_tgexp) # 最近傍で補完 library(RANN) knnImpute=preProcess(missing_tgexp, method="knnImpute") knnimputedGexp= predict(knnImpute,missing_tgexp)

17. 学習データとテストデータにデータセットを分割する。 テスト用に独立したデータセットを分割することで、モデルのパフォーマンスを適切に評 価・調整することができる。 データの分割 学習データ テストデータ データセット 分割

18. 一般的なデータ分割手法。 単純にデータセットを学習データとテストデータに分割する。 データの分割 - ホールドアウト法 tgexp=merge(patient,tgexp, by="row.names") # push sample

    ids back to the row names rownames(tgexp)=tgexp[, 1] tgexp=tgexp[,-1] set.seed(3031) # シード値の設定 # データの70%について、インデックスを取得 intrain <- createDataPartition(y = tgexp[,1], p= 0.7)[[1]] # 学習データとテストデータを分割 training <- tgexp[intrain,] testing <- tgexp[-intrain,] 他の手法 • 交差検証 (Cross Validation) • ブートストラップ法

19. 最も簡単な推定手法として、類似したサンプルのラベルと同一のラベルを割り当てること が考えられる。概念的にはクラスタリングに近い。 Rではcaret::knn3()メソッドでk最近傍法を実行することができる。 k最近傍法によるCIMPサブタイプの推定 library(caret) knnFit=knn3(x=training[,-1], # 学習データ y=training[,1], #

    学習データの正解ラベル k=5) # テストデータに対して推定を実行 trainPred=predict(knnFit,training[, -1]) 考慮すべき最近傍数（ k）の 選定が重要

20. モデルを評価する際は、いくつかの指標を用いることができる。 二値分類で使用される指標は、以下の表の通り。 モデルの評価 True Positive (TP) False Positive (FP) False

    Negative (FN) True Negative (TN) Sensitivity Specificity Precision = Predicted as CIMP Predicted as noCIMP Actual CIMP Actual noCIMP ※正例、負例に偏りがある場合、これ らの指標はうまく機能しない。 例) 90個のCIMPのnoCIMPがある場合、すべて のサンプルをCIMPとして分類するとPrecisionは 0.9となる。

21. 不均衡なデータではSpecificityやPrecisionの代わりに、 下記のような指標が用いられる。 モデルの評価 Balanced Accuracy カッパ係数 平均正解率。バランスの取れた精度の評価指標。 Brodersen, K.H.; Ong,

    C.S.; Stephan, K.E.; Buhmann, J.M. (2010). The balanced accuracy and its posterior distribution. Proceedings of the 20th International Conference on Pattern Recognition, 3121-24. 多クラス分類問題にも適用可能。ランダムな分類モデルと比較 した際にどれだけパフォーマンスが良いかを示す指標 Pr(a): 観察された一致率 Pr(e): 期待一致率 McHugh ML. Interrater reliability: the kappa statistic. Biochem Med (Zagreb). 2012;22(3):276-82. Landis, J.R.; Koch, G.G. (1977). “The measurement of observer agreement for categorical data”. Biometrics 33 (1): 159–174

22. Rでモデルの各指標を確認するには、 caret :: confusionMatrix（）関数を用いる モデルの評価 # 学習データに対して、 k=5のk-NN 実行結果を取得 knnFit=knn3(x=training[,-1],

    # training set y=training[,1], # training set class labels k=5) # 学習データに対して推測を実行 trainPred=predict(knnFit,training[, -1],type="class") # 実際の正解ラベルと推測されたラベルを比較して精度指標を取得 confusionMatrix(data=training[,1],reference=trainPred) ## Confusion Matrix and Statistics ## ## Reference ## Prediction CIMP noCIMP ## CIMP 65 0 ## noCIMP 2 63 ## ## Accuracy : 0.9846 ## 95% CI : (0.9455, 0.9981) ## No Information Rate : 0.5154 ## P-Value [Acc > NIR] : <2e-16 ## ## Kappa : 0.9692 ## ## Mcnemar's Test P-Value : 0.4795 ## ## Sensitivity : 0.9701 ## Specificity : 1.0000 ## Pos Pred Value : 1.0000 ## Neg Pred Value : 0.9692 ## Prevalence : 0.5154 ## Detection Rate : 0.5000 ## Detection Prevalence : 0.5000 ## Balanced Accuracy : 0.9851 ## ## 'Positive' Class : CIMP ##

23. テストデータについても同様に 精度指標を取得する。 モデルの評価 # 学習データに対して推測を実行 testPred=predict(knnFit,testing[, -1],type="class") # 実際の正解ラベルと推測されたラベルを比較して精度指標を取得 confusionMatrix(data=testing[,1],reference=testPred)

    ## Confusion Matrix and Statistics ## ## Reference ## Prediction CIMP noCIMP ## CIMP 27 0 ## noCIMP 2 25 ## ## Accuracy : 0.963 ## 95% CI : (0.8725, 0.9955) ## No Information Rate : 0.537 ## P-Value [Acc > NIR] : 2.924e-12 ## ## Kappa : 0.9259 ## ## Mcnemar's Test P-Value : 0.4795 ## ## Sensitivity : 0.9310 ## Specificity : 1.0000 ## Pos Pred Value : 1.0000 ## Neg Pred Value : 0.9259 ## Prevalence : 0.5370 ## Detection Rate : 0.5000 ## Detection Prevalence : 0.5000 ## Balanced Accuracy : 0.9655 ## ## 'Positive' Class : CIMP

24. kの値が異なるモデルを試すことで、最適なkの値を探索することができる。 学習データの予測誤差がどのように変化するか確認する。 モデルのチューニング set.seed(101) k=1:12 # k の値の範囲 trainErr=c() #

    学習データの予測誤差を格納するためのベクトルを定義 for( i in k){ knnFit=knn3(x=training[,-1], # 学習データ y=training[,1], # 学習データの正解ラベル k=i) # 学習データに対して予測を実行 class.res=predict(knnFit,training[, -1],type="class") # 学習データの予測誤差 trainErr[i]= 1-confusionMatrix(training[,1],class.res)$overall[1] } # k ごとの学習データの予測誤差をプロット plot(k,trainErr,type="p",col="#CC0000",pch=20) # 平滑線を描画 lines(loess.smooth(x=k, trainErr,degree=2),col="#CC0000") k の値が大きいほど予測誤差が大きく なっていることが分かる。

25. テストデータについても、kの値で予測誤差がどのように変化するか確認する。 モデルのチューニング trainErr=c() # テストデータの予測誤差を格納するためのベクトルを定義 for( i in k){ knnFit=knn3(x=training[,-1],

    # 学習データ y=training[,1], # 学習データの正解ラベル k=i) # テストデータに対して予測を実行 class.res=predict(knnFit,testing[,-1],type="class") # 学習データの予測誤差 testErr[i]=1-confusionMatrix(testing[,1],class.res)$overall[1] } # k ごとの学習データの予測誤差をプロット plot(k,trainErr,type="p",col="#CC0000",pch=20) # 平滑線を描画 lines(loess.smooth(x=k, trainErr,degree=2),col="#CC0000") # k ごとの学習データの予測誤差をプロット points(k,testErr,col="#00CC66",pch=19) lines(loess.smooth(x=k,testErr,degree=2), col="#00CC66") # legendを追加 legend("bottomright",fill=c("#CC0000","#00CC66"), legend=c("training","test"),bty="n") k の値が小さい場合は過学習を起こしている ことが分かる。 ⇒ k = 6 程度が最適

26. 予測誤差は分散（Variance）とバイアス（Bias）の和となっている。 分散が大きいモデルは複雑でより学習データにフィットしやすい柔軟性を持つ モデルの複雑さとバイアスのトレードオフ k =2 と k=12のモデルを比較すると、k=2のモデルはより分 散が大きいモデルとなっていることが分かる。

27. 予測誤差の推定は、適切な複雑さのモデルを見つけるためのモデル調整の中心とな る。データ分割手法について再考する。 データ分割戦略再考 学習、検証、テストデータへの 分割 新たに検証データというデータ分割を行い、パラメータ調整には 検証データから算出される予測誤差を使用する。 パラメータ調整の際に検証データに過学習してしまうことを避け ることができる。データが多い場合に有用。 交差検証（Cross

    Validation） 学習データ・テストデータの分割だけでなく、 モデルチューニングにおいても利用可能。

28. 交差検証を用いた学習にはcaret::train()関数を用いる 交差検証によるモデルのチューニング set.seed(17) # trainControl が学習を制御する。 # 10-fold の交差検証を行う。 trctrl

    <- trainControl(method = "cv",number=10) # k-NN モデルの学習を行う knn_fit <- train(subtype ~., data = training, method = "knn", trControl= trctrl, tuneGrid = data.frame(k=1:12)) # 交差検証精度において最良の k を出力 knn_fit$bestTune ## k ## 4 # k ごとの予測誤差をプロット plot(x=1:12,1-knn_fit$results[, 2],pch=19, ylab="prediction error" ,xlab="k") lines(loess.smooth(x=1:12,1-knn_fit$results[, 2],degree=2), col="#CC0000" ) k = 4 で予測誤差が最小となっているこ とが分かる。

29. caret::train()関数を用いると、学習ロジックの変更が簡単に行える。 ブートストラップ法によるモデルのチューニング set.seed(17) # trainControl が学習を制御する。 # 20 回のブートストラップリサンプリングを行う。 trctrl

    <- trainControl(method = "boot",number=20,returnResamp= "all") # k-NN モデルの学習を行う knn_fit <- train(subtype ~., data = training, method = "knn", trControl= trctrl, tuneGrid = data.frame(k=1:12)) # 交差検証精度において最良の k を出力 knn_fit$bestTune ## k ## 4 # k ごとの予測誤差をプロット plot(x=1:12,1-knn_fit$results[, 2],pch=19, ylab="prediction error" ,xlab="k") lines(loess.smooth(x=1:12,1-knn_fit$results[, 2],degree=2), col="#CC0000" ) k = 4 で予測誤差が最小となっているこ とが分かる。 変更部分

30. 機械学習モデルによる予測の他の目的は、どの変数が予測に重要かを知ることにあ る。モデルに依存する重要度の測定法と、モデルに依存しない方法がある。 変数重要度 モデルに依存しない 統計的検定を関連付ける、若しくは適用することで、説明変数と 目的変数の関係性の強さを定量化する。 変数ごとに独立に考慮する。 モデルに依存 学習データ・テストデータモデルチューニングにおいても 通常、予測子間の相関関係を重要度の計算に組み込むことができ

    る。

31. RではDALEX::explain（）メソッドを使用することで、変数重要度を算出することができる。 変数重要度 library(DALEX) set.seed(102) # do permutation drop-out explainer_knn<- DALEX::explain(knn_fit,

    label="knn", data =training[,-1], y =as.numeric(training[,1])) viknn=feature_importance(explainer_knn, n_sample=50,type="dif ference") plot(viknn)

32. 不均衡データへの対処法 ゲノムデータに機械学習を応用する上で共通のハードルは、大規模な不均衡データ問題である。 これらの対処には下記のような手法が用いられる。 例）深刻なクラスの不均衡を伴う別の例は、エンハンサー予測が挙げられる。 ヒトゲノムは30億塩基対の長さであり、そのほとんどがエンハンサーアノテーションと重複しないため、 「エンハンサーではない」というネガティブセットがトレーニングを圧倒する。 クラスを均衡にする サンプリング 学習データを作成する際に、同じサイズのクラスを生成するようにサンプリングを行う。多 数クラスのサンプルからダウンサンプリング、若しくは少数クラスのサンプルからアップサ

    ンプリングすることで対処することが可能。 重みづけ いくつかの手法では、少数クラスに対して不均衡度による重みづけを行うことでこの問題 に対処することができる。ロジスティックベースの手法とブースティング手法は、重みづけ を行うことのできるアルゴリズムの一例。 異なる分類スコアカッ トオフの選択 過剰な真陽性または偽陽性を最小化するような予測スコアのカットオフ値を選択する。2ク ラス分類問題の古典的な予測カットオフは0.5だが、この値を調節することでSensitivity と Specificityを最適化することができる。ROC曲線が有用。

33. 学習データ内の相関構造を持った複雑なモデルは、統計的推論が困難となる上に、学習データと将 来のテストデータの相関構造が同じでない場合過学習に繋がる危険性がある。下記のような対処法 を検討する必要がある。 マルチコ（多重共線性）の扱い PCA 学習データに対してPCAを実行すると、共線性が除かれた新規の変数が作成される。一 方で、作成される変数は元の変数の線形結合となっており、変数重要度などの解釈が困 難となる。 変数選択 変数をフィルタリングすることで、共線性の強い変数を除くことができる。

    変数の線形結合は直接相関していないが、相関はあるような場合に対処することができ ない。 その他 正則化は過学習を抑えるための手法であり、共線性の影響を弱めることもできる。また、 決定木ベースの方法は、アルゴリズムの特性上共線性の影響を受けにくい。

34. Q＆A モデルの性能を調べる上で、バリアンス、バイアスを調べるのは必須か。 訓練データへの当てはまりと過学習抑制を同時に考慮する際に、モデルの評価として使 用することが可能ですが、必須という訳ではありません。 基本的には、学習データとテストデータの学習曲線を見るだけで十分な場合が殆どかと 思います。 mlxtendパッケージのbias-variance-decomposition関数等を用いることで、ブートストラッ プ法にて得られるバリアンスとバイアスの推定値を算出することができます。 参考： https://machinelearningmastery.com/calculate-the-bias-variance-trade-off/

    A. Q.

35. Q＆A 公開データがあったら教えてほしい Awesome Public Datasetsリポジトリに生物系公開データベースへのリンクがまとめられ ているので参考にしてみて下さい。 https://github.com/awesomedata/awesome-public-datasets#biology A. Q.