メタボローム解析: パッケージroplsのtutorialをやってみた
1. はじめに
OPLS-DAを行う際の強力なツールとして、パッケージroplsの使い方をまとめていきます。OPLD-DAだけでなく、PCA、PLS(-DA),OPLSも行うことができ、R2とQ2によるquality metrics、permutation diagnosis、VIP valuesの計算、外れ値を計測するためのscoreとorthogonal distanceの計算もでき、様々なグラフも描画することができます。OPLS-DAに関しては、前回記事を参照にして下さい。
pediatricsurgery.hatenadiary.jp
2. 参考
- Tutorial
- Analysis of the Human Adult Urinary Metabolome Variations with Age, Body Mass Index, and Gender by Implementing a Comprehensive Workflow for Univariate and OPLS Statistical Analyses
- Vignette
3. Rでのワークフロー
3-1. {ropls}をインストールする
## try http:// if https:// URLs are not supported source("https://bioconductor.org/biocLite.R") biocLite("ropls") library(ropls)
3-2. データセットsacurine
data(sacurine) attach(sacurine) # sacurineはlistになっていて、下記3つのmatrixが含まれている names(sacurine) [1] "dataMatrix" "sampleMetadata" "variableMetadata"
3-3. str: Printed summary of an R object
次にstrF()でデータの中身を見てみます。
strF(dataMatrix) dim class mode typeof size NAs min mean median max 183 x 109 matrix numeric double 0.2 Mb 0 -0.3 4.2 4.3 6 (2-methoxyethoxy)propanoic acid isomer (gamma)Glu-Leu/Ile ... Valerylglycine isomer 2 Xanthosine HU_011 3.019766011 3.888479324 ... 3.889078716 4.075879575 HU_014 3.81433889 4.277148905 ... 4.181765852 4.195761901 ... ... ... ... ... ... HU_208 3.748127215 4.523763202 ... 4.634338821 4.487781609 HU_209 4.208859398 4.675880567 ... 4.47194762 4.222953354 strF(sampleMetadata) age bmi gender numeric numeric factor nRow nCol size NAs 183 3 0 Mb 0 age bmi gender HU_011 29 19.75 M HU_014 59 22.64 F ... ... ... ... HU_208 27 18.61 F HU_209 17.5 21.48 F strF(variableMetadata) msiLevel hmdb chemicalClass numeric character character nRow nCol size NAs 109 3 0 Mb 0 msiLevel hmdb chemicalClass (2-methoxyethoxy)propanoic acid isomer 2 Organi (gamma)Glu-Leu/Ile 2 AA-pep ... ... ... ... Valerylglycine isomer 2 2 AA-pep:AcyGly Xanthosine 1 HMDB00299 Nucleo
3-4. Principal Component Analysis (PCA)
sacurine.pca <- opls(dataMatrix) PCA 183 samples x 109 variables standard scaling of predictors R2X(cum) pre ort Total 0.501 8 0
- (左上) Scree plot: 第3因子くらいまでで大部分を説明できる
- (右上) Outlier
- (左下) Score plot
- (右下) Loading plot
これに群分け情報を追加してプロットしてみる。
genderFc <- sampleMetadata[, "gender"] plot(sacurine.pca, typeVc = "x-score", parAsColFcVn = genderFc, # parameter tobe used as colors ( Factor of character type or Vector of numeric type ) parEllipsesL = TRUE) # ellipseを追加するか
3-5. PLS ( - DA )
判別分析は教師あり学習なので、教師 ( 群分けデータ ) を供することが必要となるこの場合、genderFcで性別のデータを使用している。
sacurine.plsda <- opls(dataMatrix, genderFc) PLS-DA 183 samples x 109 variables standard scaling of predictors and response(s) R2X(cum) R2Y(cum) Q2(cum) RMSEE pre ort pR2Y pQ2 Total 0.275 0.73 0.584 0.262 3 0 0.05 0.05
- (左上) significance diagnostic
- (右上) inertia barplot: Inertia ( 分散 )
- (左下) outlier diagnostics
- (右下) X-score plot: 累積R2X, R2Y, Q2Yが記載される
参考: R2とQ2に関して
- PLS回帰直線を作成する際に算出されるPLSモデルの評価指標
R2
- 説明分散
- モデルの直線性を示す
- モデル構築に用いた学習データに対する当てはまりを表わす
- 値が 1 に近いほど精度の高いモデルであるといえる
Q2
- 予測的説明分散
- モデルの予測性能を示す
- モデルの予測性を表わし、値が 1 に近いほど予測的なモデルであるといえる
参考
- FT-NIR を用いたメタボリックフィンガープリンティングによる 緑茶製造における蒸熱工程管理の重要指標である「蒸し度」の推定モデルの構築
- RNAi 法におけるトランスフェクション試薬の 遺伝子発現抑制率予測モデルの構築
3-6. OPLS ( - DA )
sacurine.oplsda <- opls(dataMatrix, sampleMetadata[, "gender"], predI = 1, # 第何主成分まで使用するか、OPLSの時は1になる orthoI = NA) # number of components which are orthogonal, PLS: 0, OPLS: NA OPLS-DA 183 samples x 109 variables and 1 response standard scaling of predictors and response(s) R2X(cum) R2Y(cum) Q2(cum) RMSEE pre ort pR2Y pQ2 Total 0.275 0.73 0.602 0.262 1 2 0.05 0.05
とまあこんな感じになります。この後、もう1つのデータセットを使用してvalidateしていくのですが
trainVi <- ropls::getSubsetVi(sacurine.oplsda)
のところでerrorが出て先に進めませんでした。同様の質問がネットにもありましたが、解決策は出ていないようです。私としては判別モデルの作成までで用が住んでしまったので、この辺にしておきます。もし解決策がわかったら教えてください。