点点鼠标就能完成的NMR代谢组学实战
本文有大量的外链，推荐阅读原文
代谢组学数据分析简介
代谢组学几乎完全不涉及生物信息学最核心的序列比对，包括武汉中科院数物所波谱国重实验室等主流科研机构都是利用化学计量学和多元统计分析方法，对通过核磁共振（NMR）、液相色谱质谱联用（LC-MS）、气相色谱质谱联用（GC-MS）等手段获取的小分子代谢物（相对分子质量MW小于1000）的定量数据进行处理。除此之外只有一些商业公司在做外包服务，因此相比基因组、转录组、蛋白质组显得非常小众。
目前在代谢组学中运用较多的多元统计分析方法，包括
主成分分析（Principal component analysis, PCA）
层次聚类分析（Hierarchical clustering analysis, HCA）
非线性映射（Nonlinear innuendo analysis, NLM）
分级聚类法（hierarchical cluster analysis, HCA）等非监督分类方法
以及K最邻近法（K-nearest neighbor classification method, K-NN）
偏最小二乘法-判别分析（Partial least squares projections to latent structures - discriminant analysis, PLS-DA）
正交偏最小二乘法-判别分析（Orthogonal projections to latent structures - discriminant analysis, OPLS-DA）
人工神经网络（Artificial Neural Network, ANN
支持向量机（Support Vector Machine, SVM）等监督分类方法。
篇幅所限，有关这些方法的具体内容和意义请参考相关文献。
代谢组学数据分析最常用的软件SIMCA是1987年由瑞典Umetrics公司研究开发的一款多元变量统计分析软件。官方网站随便一个课程视频都是一、两千美刀，其他学习资源目前也都限于各生物公司的培训资料，同样价格不菲。这里我利用SIMCA官方的教程和数据操作一遍完整的流程以飨读者。
下载并安装软件
官方网站的下载页面提供了30天免费试用版SIMCA，当前更新到第15版；下载解压得到.exe格式文件，安装后即可在Windows系统上运行。由于之前为了给转录组入门实战数据腾位置，我删掉了苹果机上的Windows双系统，这次我将临时在Parallels Desktop 13试用版虚拟机上运行Windows 10系统来进行操作。
教程及数据
官网下载页虽然也提供了官方教材Analysing_NMR_metabolomics_data_using_OPLS-DA.pdf，但是和数据NMR_METABOLOMICS.xls一样都非常不便下载，不过好在网上随便搜索一下也能找到很多流传的资源。方便起见，我上传到了我的GitHub里面，点击文件标题名就可以直接下载。
数据信息
该研究以杨树的文库中编码功能未知的转录因子MYB的基因PttMYB76为例，比较MYB76杨树与野生型杨树生长过程的代谢组学特征差异。实验材料包括来自MYB76杨树（编号为Ai, Bi, Ci）和野生型杨树（编号为Di, Ei, Fi）两组共6棵树的45个样品，其中i代表生长方向，即1代表顶部，8代表底部。加上12组重复观测值，数据集中总共有57组观测值，655种变量。由于这里我们主要目的是做数据分析，因此对实验细节不加赘述，有需要者可以参考原文：
S. Wiklund et.al A new metabonomic strategy for analysing the growth process of the poplar tree. Plant Biotechnology Journal 2005 3 pp 353-362
数据导入
打开菜单栏中的File，点击New选项中的Regular project弹出浏览界面，选择点击下载好的数据集NMR_METABOLOMICS.xls文件导入数据。
这里还需要按照软件要求调整数据格式，打开菜单栏中的Edit，点击Transpose将表格翻转90度；
在打开菜单栏中的Home，依次选择第一行后点击Primary variable ID，选择第二列后点击Secondary observation ID；最后点击右上角的Finish import，保存后即完成数据导入。
主成分分析
打开新建的模型M1，点击模型界面左上角的Workset…，选择Scale界面，点击左下角Find Primary ID…右侧的三角形，选择第一行的Select all（或者Ctrl+A）全选变量名；然后点击右侧Type项中下拉按钮，选择Par选项后，点击下方的Set all，对所有变量值进行Pareto标度化。
再选择Observations界面，长按Ctrl键，选出A、B、C三棵MYB76杨树对应观测值，在Set class栏中输入2，点击Set将MYB76杨树观测值编为2组；同理将D、E、F三棵野生型杨树对应观测值编为1组；再在下方的Model typr栏中选择PCA-X，点击确定完成设置。这里要求处理组样本需要具有较高的类别编号，对照样本需要具有较低的类别编号，这样做会有利于后面的直接解释。
选择新建的M1模型，点击上方工具栏中的Autofit自动拟合最佳模型。
生成的模型产生的8个主成分，这里我们看到第8主个成分增加后，模型预测性Q2没有明显增加，说明这个主成分对于模型的贡献以及非常微小。
为了便于分析，这里我们只选取前3个主成分，可以通过点击Autofit右侧的remove删去多余的主成分。
如果误删成分，可以点击remove上方的add恢复。（我没有点错！我只是为了给大家试错，对，没错，就是这样！）
最后的结果应该是这样：
点击上方工具栏中的DModX可以选择生成折线图或直方图查看离群值。
点击上方工具栏中的Scores可以选择Scatter生成以第一主成分t1和第二主成分t2分别为横纵坐标的t1-t2得分图。
也可以选择3D生成分别以三个主成分为坐标的三维得分图。
这里两组样本数据并没有得到很好的区分，我们需要在看一下t2-t3得分图。
实现这个得分图需要点击工具栏中的Scores选项下最右侧的图标，或者右键点击之前生成的得分图，点击第二个图标，进入设置界面。
、
在Scores页面中Selected栏下，X-axis对应的Comp选择2，Series对应的Comp选择3，点击确定即可生成t2-t3得分图.
这里我们对比两个得分图可以明显看出，后者对两组样本数据的区分效果优于前者。
此外，如果仔细观察t1-t2得分图，还可以发现生长方向i值有明显地沿某一方向依次变化的趋势，表现了杨树生长过程的代谢组学特征的规律性变化。
正交偏最小二乘法判别分析
右键点击模型M1，选择New as Model 1生成新的模型M2。
选择Observations界面，长按Ctrl键，选出最后12组重复观测值，右键选择Exclude；再在下方的Model typr栏中选择OPLS-DA，点击确定完成设置。
再次点击Autofit自动拟合最佳模型，得到一个1+4的OPLS模型。
点击上方工具栏中的Scores，选择Scatter即可生成OPLS-DA得分图。
主成分分析和正交偏最小二乘法判别分析的对比
显然正交偏最小二乘法判别分析的分类效果更佳！
解释正交偏最小二乘法判别分析模型
点击上方工具栏中的Loadings，选择Line生成图片后，右键点击图片，点击Properties，在Loadings页面中选择Selected栏下，X-axis对应的Item选择Var ID，Series对应的Item选择pq。
点击确定生成载荷图。
点击上方菜单栏中的Plot/List，点击工具栏中的Scatter，在Data series页面中Select data type项选择Observations and loadings，X-axis对应的Item选择p，Series对应的Item选择p(corr)。
点击确定生成S-plot。
标记点为可能的生物标志物。
诊断正交偏最小二乘法判别分析模型
模型概要
R2X(cum)是模型所解释的X预测和正交变异之和，即0.157+0.613=0.769，也可以被解释为X总变异的76.9%；
R2Y(cum) 是模型所解释Y的总变异, 此处为0.977；
Q2(cum)衡量模型预测好坏，此处为0.941；
预测变异 = X变异中与Y相关的部分
A代表于X和Y之间的相关主成分的数量。如果仅使用一个响应向量，则A总是1；
R2X衡量X变异中与Y相关部分，这里是0.157；
正交变异 = X变异中与Y无关的部分
A在此对应于不相关（正交）主成分的数量。每个正交分量都被表示出来并且可以被单独解释；
R2X是衡量X变异中与Y不相关的部分。每个主成分的信息都被分别表示出来；
首行中粗体的R2X（cum）是X变异中与Y无关部分的总和，此处为0.613；
噪声 = 1- 0.157 – 0.613 = 0.23 = 23%
_补充_
响应排序检验
点击上方菜单栏中的Analyze，点击工具栏中的Permutations，在Permutation Plot页面中Number of permutations to项选择999。
点击确定生成响应排序检验示意图，通常R2Y和Q2与直线的斜率越接近水平，模型越有可能过拟合。
变量权重重要性排序
点击工具栏中的VIP，选择VIP predictive。
即直接生成变量权重重要性排序图。
右键点击图片，选择Create栏中List，可生成变量权重重要性排序表，可以选择VIP值大于1的作为感兴趣的物质进行代谢通路分析。
许国旺. 代谢组学:方法与应用[M]. 科学出版社, 2008.
贾伟. 医学代谢组学[M]. 上海科学技术出版社, 2011.
想更多了解此文？