mixOmics

2020-04-27 本文已影响0人 nnlrl

在分析生物学问题时，我们会想着从多个角度入手分析样本。当我们拥有多种不同的组学数据的时候，组学的独立分析忽略了不同特征之间的关系，可能会遗漏关键的信息。mixOmics为我们提供了一种思路，我们可以将不同的组学数据在进行预处理后共同分析，并将研究问题主要放在特征的选择上，它可以应付多组学大样本量的分析，通过将特征投影到低维的子空间中捕获关键的变异源以便进行后续研究。mixOmics采用了有监督判别分析，在降维的同时识别最具判别性的因素，并将其用于预测新的样本。

在mixOmics中主要包括了三种分析模式：

在应对单类组学数据时，可以使用PCA，sPCA进行降维以及特征筛选，而面对有监督的数据类型时，可以转而使用PLSDA，sPLSDA在考虑样本分类信息的情况下将样本投射到低维空间中并进行分析
在应对多种组学数据时，可以使用DIABLO将不同数据联合起来，例如基因组，mRNA数据，miRNA数据，蛋白质表达谱以及代谢表达谱的联合分析
多个组学的相同分析类型MINT，比如对不同的mRNA数据进行联合分析

同时mixOmics也包括了多种可视化手段用以展示结果，针对样本的plotIndiv, plotArrow, plotDiablo, 针对特征的plotVar, plotLoadings, circosPlot, cim, network以及用于模型准确性验证的auroc, plot.tune, plot.pref

PLSDA和sPLSDA

可用于定量和定性

当我们正在分析一个数据集（如转录组学数据），我想将我的样本分类为已知的组，并预测新样本的类别。此外，我有兴趣找出导致这种变异的关键变量。
PLS-DA(Partial Least Squares Discriminant Analysis)偏最小二乘法判别分析，是多变量数据分析技术中的判别分析法，通过对主成分进行适当的旋转，可以有效的对组间观察值进行区分。自分析过程中，我们需要指定样本矩阵X以及样本分组信息Y

对于这组方法，需要选择三个参数：

保留的主成分数量
为sPLS-DA在每个主成分上选择的变量数keepX
模型的评估
作者建议使用n次重复的K折交叉验证选用不同数量的主成分模型进行结果评估，尽量选择准确率最高的模型作为后续分析

library(mixOmics)
data(srbct)
X <- srbct$gene
Y <- srbct$class 
summary(Y) ## class summary

## EWS  BL  NB RMS 
##  23   8  12  20

MyResult.plsda2 <- plsda(X,Y, ncomp=10)
set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat
MyPerf.plsda <- perf(MyResult.plsda2, validation = "Mfold", folds = 3, 
                  progressBar = FALSE, nrepeat = 10) # we suggest nrepeat = 50

plot(MyPerf.plsda, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")

关于第1项，该图输出分类错误率或平衡分类错误率，当每组样本数不平衡时，根据三个预测距离的标准差。这里我们可以看到，对于BER和最大距离，对于ncomp=3，似乎可以获得最佳性能（即低错误率）。
关于第2项，我们现在使用tune.splsda来评估在每个组件上选择的变量的最佳数量。我们首先建立一个keepX值的网格，每个组件一次评估一个组件。与上述类似，我们运行3折CV重复10次，最大距离预测定义如上。

list.keepX <- c(5:10,  seq(20, 100, 10))

set.seed(30) # for reproducbility in this vignette, otherwise increase nrepeat
tune.splsda.srbct <- tune.splsda(X, Y, ncomp = 3, # we suggest to push ncomp a bit more, e.g. 4
                                 validation = 'Mfold',
                                 folds = 3, dist = 'max.dist', progressBar = FALSE,
                                 measure = "BER", test.keepX = list.keepX,
                                 nrepeat = 10)   # we suggest nrepeat = 50

error <- tune.splsda.srbct$error.rate
ncomp <- tune.splsda.srbct$choice.ncomp$ncomp

select.keepX <- tune.splsda.srbct$choice.keepX[1:ncomp]

随后便可以使用选择好的参数构建模型

MyResult.splsda<- splsda(X, Y, ncomp = ncomp, keepX = select.keepX)

#样本
plotIndiv(MyResult.splsda, ind.names = FALSE, legend=TRUE,
          ellipse = TRUE, title="sPLS-DA - final result")
#由于我们选择了前三个主成分，我们也可以绘制3D图片展示所有主成分，注意安装rgl包
#plotIndiv(MyResult.splsda, style="3d")

#特征
plotVar(MyResult.splsda)                            # 3 Plot the variables

#添加背景
background <- background.predict(MyResult.splsda, comp.predicted=2,
                                 dist = "max.dist") 
plotIndiv(MyResult.splsda, comp = 1:2, group = srbct$class,
          ind.names = FALSE, title = "Maximum distance",
          legend = TRUE,  background = background)

#AUC
auc.plsda <- auroc(MyResult.splsda)

当我们构建了一个判别模型并通过AUC验证了模型的准确性后，我们还需要知道对结果贡献最大的特征，挑选出来并用于后续的分析

#在指定主成分中选择贡献最大的特征
selectVar(MyResult.splsda2, comp=1)$value

plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean')

DIABLO

DIABLO利用广义典型分析使用相同分析中多个组学的数据进行整合分析。这是mixOmics进行整合分析主要使用的方法，同时也有更多的函数支持。

#导入示例数据，包括mrna，mirna以及蛋白质谱
library(mixOmics)
data(breast.TCGA)
# extract training data and name each data frame
X <- list(mRNA = breast.TCGA$data.train$mrna, 
          miRNA = breast.TCGA$data.train$mirna, 
          protein = breast.TCGA$data.train$protein)
Y <- breast.TCGA$data.train$subtype
summary(Y)

在DIABLO训练过程中，我们需要注意3个参数：

设计矩阵
主成分以及keepX与单因素的PLSDA相同

list.keepX <- list(mRNA = c(16, 17), miRNA = c(18,5), protein = c(5, 5))

MyResult.diablo <- block.splsda(X, Y, keepX=list.keepX)
plotIndiv(MyResult.diablo) ## sample plot

plotVar(MyResult.diablo) ## variable plot

我们也可以进行一协调整

#根据亚型分组
plotIndiv(MyResult.diablo, 
          ind.names = FALSE, 
          legend=TRUE, cex=c(1,2,3),
          title = 'BRCA with DIABLO')

#某些数据集中可以省略标签以提高可读性。例如她，我们只显示蛋白质的名称
plotVar(MyResult.diablo, var.names = c(FALSE, FALSE, TRUE),
        legend=TRUE, pch=c(16,16,1))

#不同数据之间两两比较绘图，展示分类效果以及数据之间的相关性
plotDiablo(MyResult.diablo, ncomp = 1)

#circos图表示不同类型变量之间的相关性，cutoff限制展示的数量
circosPlot(MyResult.diablo, cutoff=0.7)

#cimDiablo使用层次聚类并按特征的表达水平绘制热图，行标签指定样本，列标签指定主成分中贡献最大的特征
cimDiablo(MyResult.diablo, color.blocks = c('darkorchid', 'brown1', 'lightgreen'), comp = 1, margin=c(8,20), legend.position = "right")

# 还可以使用network绘制网络图
#network(MyResult.diablo, blocks = c(1,2,3),color.node = c('darkorchid', 'brown1', 'lightgreen'), cutoff = 0.6, save = 'jpeg', name.save = 'DIABLOnetwork')

#AUC用于评估模型的准确性，同时还可以使用混淆矩阵
Myauc.diablo <- auroc(MyResult.diablo, roc.block = "miRNA", roc.comp = 2)
confusion.mat <- get.confusion_matrix(
                truth = breast.TCGA$data.test$subtype, 
                predicted = Mypredict.diablo$MajorityVote$centroids.dist[,2])
kable(confusion.mat)

在构建完模型后，我们同样也可以使用selectVar, plotLoadings来进行特征的筛选

mixOmics

PLSDA和sPLSDA

DIABLO

猜你喜欢

热点阅读