読者です 読者をやめる 読者になる 読者になる

Note of Pediatric Surgery

腸内細菌、R、ときどき小児外科

メタボローム解析: パッケージ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. 参考

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 に近いほど予測的なモデルであるといえる
参考

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が出て先に進めませんでした。同様の質問がネットにもありましたが、解決策は出ていないようです。私としては判別モデルの作成までで用が住んでしまったので、この辺にしておきます。もし解決策がわかったら教えてください。