mixOmicsの概要

mixOmicsとは

次世代シーケンサー（NGS）などの普及により、網羅的にデータを取得する技術「オミクス」が今や生物学研究の主流になっています。しかし、サンプル数（N）に対して非常に多くの変数（P）を扱う事になるオミクス解析では、適切な多変量解析の手法を選択して利用することが重要です。また、１種類のオミクスだけはで説明できない現象も多く、今後はゲノムやトランスクリプトーム、プロテオーム、メタボロームといった複数のオミクスの情報を統合することで、複雑系ネットワークの解析を進めていくことが期待されています。そうなると、ますます多変量解析の統計手法の重要性が増してくると言えます。

出典：http://mixomics.org

そこで登場するのが、今回紹介するmixOmicsです。RのパッケージであるmixOmicsは、マルチオミクスデータの統合解析のために開発されました。mixOmicsのwebサイト（http://mixomics.org：英語）では、mixOmicsの特徴が以下のようにまとめられているように、特に変数選択（バイオマーカーの抽出）に適したパッケージになっています。例えば、トランスクリプトームやメタボロームのデータの中から、病気に関係していそうな変数（遺伝子や代謝物）をまとめて取り出したいような場合に役に立つのです。

mixOmics offers a wide range of multivariate methods for the exploration and integration of biological datasets with a particular focus on variable selection

DIABLOとMINT

出典：http://mixomics.org

mixOmicsのパッケージには複数の手法が含まれていますが、その中でも２つのフレームワークがオミクスの統合解析用に用意されています。そのフレームワークとは、DAIBLO（Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omits studies）とMINT（Multivariate INTegrative method）です。

いずれも複数のオミクスデータを統合するための手法ですが、違いは統合の仕方がN-integrationかP-integationにあります。

出典：http://mixomics.org/mixdiablo/

N-integrationとは、同じサンプルに対して異なるオミクスを利用して得られた複数のデータセットを統合する方法です。いずれのデータセットでもサンプル数（N）が等しいのに対して、ゲノムやプロテオームなどデータセットごとに変数の数（P）は異なっているのが特長です。

一方で、P-integationは、同じ変数を測定した異なるデータセットを統合する方法です。独立した異なる研究を統合するために用います。統合するいずれのデータセットでも、変数が同じ（Pが等しい）のに対して、Nは等しくなくても良いのが特長です。

ここでは、異なる種類のオミクスデータを統合して解析するために、前者のDIABLOに注目をして流れを紹介したいと思います。

解析の流れ

パッケージのインストール

さっそく、Rでパッケージのインストールを行います。パッケージはBioconductorに公開されているので、BiocManagerパッケージがインストールされていない場合は先にこちらをインストールします。

## BiocManagerがインストールされていない場合
# install.packages("BiocManager")

## mixOmicsのインストール
BiocManager::install('mixOmics')

インストールが正しく行えた場合、パッケージの読み込みを行うと以下のようにパッケージのバージョン（6.16.0）が表示されます。

> library(mixOmics)

Loaded mixOmics 6.16.0
Thank you for using mixOmics!
Tutorials: http://mixomics.org
Bookdown vignette: https://mixomicsteam.github.io/Bookdown
Questions, issues: Follow the prompts at http://mixomics.org/contact-us
Cite us:  citation('mixOmics')

読み込みを行った時にエラーが表示される場合などは、Rのバージョンが古い可能性もあるので確認してみてください（https://cran.r-project.org)

データの読み込み

では、DIABLOを行っていきます。

まずは、先ほどインストールしたパッケージを読み込みます。

library(mixOmics)

次に、解析で使うデータを読み込みます。ここでは、TCGA（The Cancer Genome Atlas）と呼ばれる、米国がん研究所と米国ヒトゲノム研究所の共同プロジェクトによるオープンソースのデータを利用します。mixOmicsのパッケージの中には、このTCGAの公開データから作った、乳がんの患者に関するマルチオミクスデータが用意されています。このデータは、150サンプルのトレーニングデータと70サンプルのテストデータに分かれています。それぞれには、3つのデータセット（mRNA、miRNA、タンパク質）と各サンプルのカテゴリ（Basal, Her2 and LumAのどれか）に関するデータが含まれています。今回行うのは、これらのデータセットを統合して、乳がんのサブタイプ(サンプルのカテゴリ)を判別できるような変数を見つけ出してくるということです。

data('breast.TCGA')

データを読み込んだら、トレーニングデータの中の３つのオミクスデータセットを統合したリストを作ります。

data = list(mRNA = breast.TCGA$data.train$mrna, 
            miRNA = breast.TCGA$data.train$mirna, 
            proteomics = breast.TCGA$data.train$protein)

また、各サンプルのカテゴリ情報をYとしておきます。

Y = breast.TCGA$data.train$subtype

これで、データの準備は完了です。

パラメータの設定（チューニング）

モデル構築にはblock.splsda関数を用います。その際に、データとサンプルのカテゴリ情報以外に、引数として、ncomp、list.keepX、designを渡す必要があります。

sgccda.res = block.splsda(X = data, Y = Y, ncomp = ncomp, 
                          keepX = list.keepX, design = design)

① Design matrix

潜在変数間の相関関係とカテゴリの判別性能のバランスを決めるのが、Design matrixと呼ばれるハイパーパラメータです。3つのデータセットを統合する今回の例では、3×3の対称な行列になります。このDesign matrixの各要素の値を大きく（0.5〜1）するとデータセット間の相関が最大化され、小さく（0〜0.5）するとYの識別性が最大化されます。データ間の相関に関する事前知識を使って値を選択したり、解析の目的に応じて値を設定したりしますが、ここでは全てのデータセット間で0.1と設定しておきます。

design = matrix(0.1, ncol = length(data), nrow = length(data), 
                dimnames = list(names(data), names(data)))
diag(design) = 0

design

以下のような行列が設定されます。

##            mRNA miRNA proteomics
## mRNA        0.0   0.1        0.1
## miRNA       0.1   0.0        0.1
## proteomics  0.1   0.1        0.0

② Number of components

次に、潜在変数の成分数を設定する必要があります。任意に設定することもできますが、pref関数を用いたクロスバリデーションによるパラメータのチューニング法を紹介します。以下では、10分割のクロスバリデーションを10回繰り返しています（folds = 10, nrepeat = 10）。

sgccda.res = block.splsda(X = data, Y = Y, ncomp = 5, 
                           design = design)

set.seed(123) 
perf.diablo = perf(sgccda.res, validation = 'Mfold', folds = 10, nrepeat = 10)

plot(perf.diablo) 

ncomp = perf.diablo$choice.ncomp$WeightedVote["Overall.BER", "centroids.dist"]

④ keepX

DIABLOでは、それぞれのデータセットから変数選択を行います。それぞれのデータセットから選択する変数の数をkeepXに格納します。ここでも、チューミング用の関数tune.block.splsdaを用いたクロスバリデーションにより値を決定します。先程と同様に、10分割のクロスバリデーションを10回繰り返します。test.keepXに候補となる値を入れます。PCのCPU性能によっては、候補の値が多すぎたり、クロスバリデーションの計算回数が多すぎたりすると計算時間がかなり多くなります。適宜調整してくいださい。

test.keepX = list (mRNA = c(5,10,15),
                   miRNA = c(5,10,15),
                   proteomics = c(5,10,15))

BPPARAM <- BiocParallel::MulticoreParam(workers = parallel::detectCores()-1)

tune.TCGA = tune.block.splsda(X = data, Y = Y, ncomp = ncomp, 
                              test.keepX = test.keepX, design = design,
                              validation = 'Mfold', folds = 10, nrepeat = 10,
                              BPPARAM = BPPARAM, dist = "centroids.dist")


list.keepX = tune.TCGA$choice.keepX
list.keepX

以上でパラメータの設定は完了です。チューニング②と③を飛ばす場合はこちら。

sgccda.res = block.splsda(X = data, Y = Y, ncomp = 2, 
                          keepX = list(mRNA=c(10,10),miRNA=c(10,10),proteomics=c(10,10)), design = design)

モデル構築と結果の出力

設定したパラメータを用いてモデルを構築します。モデル構築にはblock.splsda関数を用います。

sgccda.res = block.splsda(X = data, Y = Y, ncomp = ncomp, 
                          keepX = list.keepX, design = design)

これでsgccda.resに結果が格納されました。

選択された変数はselectVar関数で取得できます。

selectVar(sgccda.res, comp = 1)

ここから、結果をグラフに出力していきます。mixOmicsには沢山のグラフ出力用の関数が用意されています。順に紹介していきます。

まずはplotDAIBLO関数で、第一成分間の相関を確認します。

plotDiablo(sgccda.res, ncomp = 1)

plotIndiv関数を用いると、サンプルを潜在変数空間に射影することができます。

plotIndiv(sgccda.res, ind.names = FALSE, legend = TRUE, ellipse = TRUE)

plotVar関数で選択された変数を潜在変数空間に射影することもできます。

plotVar(sgccda.res, var.names = FALSE, style = 'graphics', legend = TRUE, 
        pch = c(16, 17, 15), cex = c(2,2,2), col = c('darkorchid', 'brown1', 'lightgreen'))

circosPlot関数は、選択された変数間の相関を可視化します。また、line = TRUEに設定しておくことで、各変数のカテゴリごとの発現レベルを表示することもできます。cutoffでは相関の足切りの値を指定します。

circosPlot(sgccda.res, cutoff = 0.7, line = TRUE, size.labels = 1.5)

network関数では、選択された変数のネットワークのグラフを出力することができます。

network(sgccda.res, blocks = c(1,2,3), cutoff = 0.4)

plotLoadings関数では、各成分に対する選択変数の負荷量（loading）を表示することができます。横軸は負荷量の大きさを表し、色はどのカテゴリで最も発現しているかを表しています。

plotLoadings(sgccda.res, comp = 1, contrib = 'max', method = 'median')

cimDIABLO関数は、クラスターイメージマップ（ヒートマップ）を出力します。変数が列、サンプルが行に相当します。凡例が重なってしまう場合は、size.legendやmargin引数を適当に指定してあげます。

cimDiablo(sgccda.res, size.legend = 1, margin=c(5,15))

モデルの評価

次に構築したモデルの精度の評価を行います。再度pref関数を用いて、10分割のクロスバリデーションを10回繰り返して誤差を計算します。

perf.diablo = perf(sgccda.res, validation = 'Mfold', M = 10, nrepeat = 10, dist = 'centroids.dist')
perf.diablo$MajorityVote.error.rate

## $centroids.dist
##                  comp1      comp2
## Basal       0.02000000 0.03777778
## Her2        0.22000000 0.13666667
## LumA        0.05866667 0.01866667
## Overall.ER  0.07933333 0.04800000
## Overall.BER 0.09955556 0.06437037

auroc関数を用いると、ROC曲線を出力することもできます。

auc.splsda = auroc(sgccda.res, roc.block = "miRNA", roc.comp = 2)

外部データのカテゴリを予測

Predict関数を使って、テストデータのカテゴリを予測することができます。

data.test.TCGA = list(mRNA = breast.TCGA$data.test$mrna, 
                      miRNA = breast.TCGA$data.test$mirna)

predict.diablo = predict(sgccda.res, newdata = data.test.TCGA)

混同行列を出力することで、予測されたカテゴリが実際のカテゴリに一致していたかを確認します。

confusion.mat = get.confusion_matrix(truth = breast.TCGA$data.test$subtype, 
                     predicted = predict.diablo$WeightedVote$centroids.dist[,2])
confusion.mat

##      predicted.as.Basal predicted.as.Her2 predicted.as.LumA
## Basal                 20                 1                 0
## Her2                   0                14                 0
## LumA                   0                 1                34

get.BER(confusion.mat)

## [1] 0.06825397

以上で、DIABLOの主な機能を紹介しました。

参考

公式webサイト：http://mixomics.org
mixOmics vignette：https://mixomicsteam.github.io/Bookdown/index.html