代谢组学数据通常是OPLS-DA或者PLS-DA来代替PCA

正常情况下我们的转录组表达量矩阵数据分析之前，会先出三张图，作为质量控制。我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

• 左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
• 中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
• 右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。

但是也有例外，比如我们的转录组表达量矩阵来自于一个没有分组的队列研究，TCGA数据库里面就有很多好几百个样品的单一癌症，实际上我们需要对这个单一癌症队列的异质性进行探索，算法非常多。如果我们有生物学背景，就可以根据指定生物学功能的基因集合，去对这个队列进行分群， 比如：

• 基于基因集的样品队列分组之层次聚类
• 基于基因集的样品队列分组之PCA，
• 基于基因集的样品队列分组之gsea等打分

实际上，这个过程在代谢组学数据里面就是OPLS-DA代替PCA，因为代谢组学矩阵即使我们有很明确的分组信息，它全局PCA通常是没办法像转录组表达量矩阵那样的成为一个三张图，详见：在生信技能树的教程：《你确定你的差异基因找对了吗？》

代谢组学数据多变量统计分析方法大致可以分为两类：

• 一类为非监督的学习方法，即在不给定样本标签的情况下对训练样本进行学习，如PCA、非线性映射（NLM) 等；
• 另一类为有监督的学习方法，即在给定样本标签的情况下对训练样本进行学习，如偏最小二乘判别分析（PLS-DA) 、基于正交信号校正的偏最小二乘判别分析（OPLS-DA) 、人工神经网络（ANN) 、支持向量机（SVM) 等。其中，PCA、PLS-DA和OPLS-DA是目前代谢组学领域中使用最为普遍的多变量统计分析方法。

我们仍然是选择前面的基于基因集的样品队列分组之层次聚类作为例子，可以看到：

(file = 'step1-output.Rdata')
cg=c('CD3D','CD3G CD247','IFNG','IL2RG','IRF1','IRF4','LCK','OAS2,STAT1')
cg
cg=cg[cg %in% rownames(dat)]
library("FactoMineR") #画主成分分析图需要加载这两个包
library("factoextra")  
dat.pca <- PCA(as.data.frame(t(dat[cg,]))  ) 
head(dat.pca$ind$coord[,1:2])

group_list=ifelse(dat.pca$ind$coord[,1] < 0 ,'neg','pos')

dat.pca <- PCA(as.data.frame(t(dat ))  )
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind =  group_list, # color by groups 
             addEllipses = T,
             legend.title = "Groups"
)

我们使用局部基因（免疫相关基因），无论是通过层次聚类，还是PCA分析，还是gsva等打分函数，都是可以把样品分成比较清晰的两个组别，它们两个组确实是在局部基因（免疫相关基因）上面的特征不一样，但是这样的局部基因（免疫相关基因）的分组，在全局水平上面是无法分开的！

在全局水平上面是无法分开

这个其实就是很多人咨询的一个问题， 为什么他自己的分组，在我们的生信技能树的教程：《你确定你的差异基因找对了吗？》发现没办法很好的区分开来，让他有点担心自己的数据分析能否继续。

其实在绝大部分代谢组数据里面，我们的分组，都是不太可能在全局PCA里面区分开来，所以有基于正交信号校正的偏最小二乘判别分析（OPLS-DA) 来代替PCA，有点类似于我们前面的使用局部基因（免疫相关基因）后的PCA，它天然就可以把我们的生物学分组很好的区分开来。

首先在bioconductor有一个包mixOmics可以做 偏最小二乘判别分析（PLS-DA) 分析，代码如下所示：

# 我们的 group_list 分组，就是这个样品队列的免疫高低
library(mixOmics)
pls_results <-plsda(as.data.frame(t(dat)), 
                     group_list ,
                     ncomp = 2)  
plotIndiv(pls_results , comp = c(1,2),
          group = group_list, style = 'ggplot2',ellipse = TRUE, 
          size.xlabel = 20, size.ylabel = 20, 
          size.axis = 25, pch = 15, cex = 5)

出图如下所示：

两个分组经过（PLS-DA) 分析

另外，在bioconductor有一个包可以做：The ropls R package implements the PCA, PLS(-DA) and OPLS(-DA) approaches，我这里就不再赘述。

总结一下

一个肿瘤队列，它全部的基因一般来说有两万多个，这样的矩阵充分的展现了肿瘤病人的异质性，这个队列可以有多种分组方式。但是如果我们仅仅是感兴趣免疫高低，那么可以强行把肿瘤病人分成免疫高低，类似的可以分成TP53基因突变与否的两个组。这样的分组，虽然说没有全局的差异，但是局部的差异仍然是有的。

如果是全局的PCA无法合理的展现你的分组， 就可以提取局部基因进行可视化，或者借助 偏最小二乘判别分析（PLS-DA) 、基于正交信号校正的偏最小二乘判别分析（OPLS-DA) 等方法帮你挑选合适的基因去可视化你的分组。

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。