CN113049664B - 一种基于质谱代谢组学的通路分析建模方法 - Google Patents

Abstract

本发明提供一种质谱代谢组学的代谢通路分析建模方法，包括如下步骤：S1：采集正常生物和患病生物的生物样本，对采集样本得到的质谱转化成mzML文件，对质谱文件进行中心化、去噪以及对齐处理，得到样本的代谢物的二维数据矩阵；S2：对样本的代谢物的二维数据矩阵以及样本类别矢量矩阵，做中心化和Unite Variance归一化处理；S3：获取代谢物‑通路映射关系，并通路权重系数进行优化；S4：进行通路排序，调整惩罚因子并确定通路被选频率，利用通路的被选频率对通路进行排序；本发明提供的方法，通过将代谢物分组为相互之间有重叠的通路，将偏最小二乘法和Group Lasso相结合建立回归模型，并引入通路权重系数和惩罚因子，实现基于“分组稀疏”的多元回归。

Description

一种基于质谱代谢组学的通路分析建模方法

技术领域

本发明涉及数据分析领域，特别是指一种基于质谱代谢组学的通路分析建模方法。

背景技术

现代高通量的质谱分析技术为我们提供生物体丰富的分子信息，使得我们可以系统研究机体的内源性分子在外界刺激、病理生理学改变以及基因突变等条件下的变化。基于质谱的代谢组学作为综合分析机体代谢的现代分析技术，在学术研究领域被广泛采用，尤其是生物医学和植物科学等领域。

数据分析是质谱代谢组学研究中的关键步骤，近二十年来人们提出了各种的数据分析策略和工具，用于解读数据中生物学差异，揭示其潜在的生物生理病理机制。传统的数据分析策略集中于识别在两个感兴趣的条件之间表现出差异的个体代谢物，然而其结果无法关联上生物环境(如，代谢物通路)，需要结合过表达分析(Over-representationanalysis，ORA)将得到的差异代谢物映射到预先定义的生物通路中，从而进一步富集得到扰动通路。这些基于代谢物的通路分析方法忽略了代谢物浓度，结果受到差异代谢物识别结果影响很大。与ORA不同，定量富集分析(Quantitative enrichment analysis，QEA)无需进行差异代谢物识别，直接利用Global-test方法统计来自同一通路的一组代谢物在两类样本间是否存在显著差异。从系统生物学的角度来看，许多的代谢物参与多条代谢通路，即通路间存在重叠，通路之间并不是独立的，它们彼此关联构成整个代谢网络。此外，通路在某些特性上是不等价的，例如不同通路，其包含代谢物的数量、“通路-谢物”对应关系等均存在差异。这些对通路分析提出了三个方面的挑战。第一，如何考虑代谢物和通路之间的相互关系；第二，对于通路间重叠代谢物，如何分配其在各通路中的权重；第三，如何消除各通路间的差异，使得不同通路之间具有可比性。

针对第一个问题，我们可以通过整合所有代谢物来构建一个综合模型来实现，代谢物和通路之间的相互作用被表征成多个通路同时建模来解决。近期研究人员利用多块偏最小二乘(multi-block partial least square，MB-PLS)分析方法对的所有的“组”数据(通路数据)进行建模，使用通路重要性投影(pathway importance in projection,PIP)参量来评估各通路的重要性，并证明了这类基于通路水平(“组”数据)的方法较传统的基于代谢物水平的方法可靠性更高。对于第二问题，通常有两种策略，一种是从数据集中移除所有重叠的代谢物；另一种是将重叠的代谢物复制/均摊到其涉及的各条通路中。然而，这两种策略在筛选扰动通道时都会导致较高的假阳性率(FDR)。对于最后一个问题，有人建议在多条通路联合建模时引入通路内代谢物数量的平方根成做为通路的权重系数，以减少通路之间的不等价性。然而，代谢物数量并不是导致通路间不可比的唯一原因，选择合适的通路权重系数仍是一个挑战。

从网络角度来看，许多的代谢物参与多条代谢通路，即存在通路重叠问题。虽然这类重叠代谢物的代谢水平是由涉及的多条通路共同作用的，但当这类代谢物的代谢水平出现扰动时，扰动通常只是由一条(或部分)通路引起。现有的通路分析方法往往将重叠代谢物的扰动得分均摊在其涉及的各条通路中，进而增加了分析结果的假阳率；其次，由于各通路包含代谢物数量不同、通路间的存在交叠等情况，在多通路联合建模时，通路权重系数难以准确估计；另外，由于代谢组学研究样本量有限，分析结果往往存在一定的随机性。

发明内容

本发明的目的在于克服现有技术的不足，针对当前质谱代谢组学通路分析中的通路重叠问题，提出了一种基于重叠通路的通路分析方法，该方法整合所有代谢物构建出一个综合模型，将代谢物分组为相互之间有重叠的通路，采用权重系数优化算法来分配各个通路的权重，有效地解决了上述背景技术中的问题。

本发明采用如下技术方案：

一种质谱代谢组学的通路分析方法，其特征在于，包括如下步骤：

S1：对采集得到的质谱转化成mzML文件，对质谱文件进行中心化、去噪以及对齐处理，得到样本的代谢物的二维数据矩阵；

S2：对样本的代谢物的二维数据矩阵以及样本类别矢量矩阵，做中心化和单位方差归一化处理；

S3：根据京都基因和基因组百科全书据库(Kyoto Encyclopedia of Genes andGenomes,KEGG)，建立“代谢物-通路”间的一对多关系；构建ogPLS模型，对通路权重系数进行优化；

S4：进行通路排序，调整惩罚因子并确定通路被选择的频率，利用通路的选择频率对通路进行排序。

具体地，所述步骤S3对代谢物-通路映射关系的通路权重系数进行优化，具体包括：

S31：初始化通路权重系数ω_g：

S32：初始化通路选择计数器c_g：c_g＝0,g＝1,2,…,G，通路选择的经验分布为：

S33：通路权重系数的调整：计算通路选择的经验分布Π^ω和零分布

之间的距离：

若D(Π^null|Π^ω)<∈，停止步骤S33，其中∈≥0表示收敛精度；

从g＝1到g＝G，循环

S34：循环迭代：转步骤S32。

具体地，所述步骤S4进行通路排序，调整惩罚因子并确定通路被选择的频率，利用通路的选择频率对通路进行排序，具体包括：

S41：确定惩罚因子λ值：

初始化：令λ＝1.0，K＝min(N,M)，β>0表示学习率，S＝φ表示所选通路集；

ogPLS建模：用ogPLS方法计算

若

则S＝{S,g}；

调整惩罚因子λ：如果|S|≠L，L为所要选择的通路数量，则更新λ＝λ+β·(|S|-L)；回到ogPLS建模步骤继续ogPLS建模直至满足|S|＝L；

S42：确定通路被选择的频率：

Bootstrapping建模，初始化通路选择计数器c_g＝0,g＝1,2,…,G；从n＝1到n＝1000循环：对数据进行50％的随机采样构建一个ogPLS模型M_n；若通路

被模型M_n选择，则c_g＝c_g+1；

S43：利用通路的选择频率对通路进行排序：

利用通路选择的频率对通路进行排序F^λ＝(c₁,c₂,…,c_G)/N，通路被选择的频率越高，该通路的排序越靠前。

具体地，构建所述ogPLS模型的方法，具体包括：

S421：初始化通路特定基向量

N为样本数量，M为N个样本中代谢物类别的数量；

S422：计算ogPLS模型的基向量u_k：从k＝1到K循环：

从g＝1到G计算

其中

Z^(g)表示组

的收缩矩阵，

是组

中检测到的代谢物，i＝1,2,…,M是在整个数据集中检测到的代谢物；

计算

若‖u_k‖₂≠0，则计算：

t_k＝Xu_k，

否则‖u_k‖₂＝0，则K＝k；

S423：建立回归模型：

其中，U＝(u₁，u₂，…，u_K)，X表示N个样本M个代谢物的二维数据矩阵，y表示原始响应。

由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：

(1)本发明的生物样本是正常生物和患病生物的生物样本(可以是血液、尿液、粪便等)，此类生物样本是常用的检验样品，其收集和处理过程简单方便，花费低，便于推广应用。

(2)现有技术中的方法大多假设通路与通路间是正交(不重叠)的，但在生物学中代谢物不只参与一个通路，通路间是交叠的，本发明将代谢物分组为交叠的通路，在生物学上具有可解释性。

(3)本发明将PLS和Group Lasso相结合建立回归模型，既恰当地处理了代谢组学数据中的多重共线性，又实现了多元回归的分组稀疏，有效地降低了通路识别结果的假阳性率。

(4)本发明通过引入通路权重系数实现去偏目的，解决了在改变通路的识别中，某些通路天生具有比其他通路高的被选择的可能性而使通路选择结果产生偏差的问题。

以下结合附图和具体实施例对本发明作进一步详述。

附图说明

图1为本发明实施例提供的各通路选择频率的示意图；

图2为本发明与其他两种方法(PIP和Globaltest)的累积代谢物覆盖率的对比结果；

图3为本发明与其他两种方法(PIP和Globaltest)的累积代谢物重叠率的对比结果。

具体实施方式

以下通过具体实施方式对本发明作进一步的描述。

本发明解决其技术问题所采用的技术方案是：通过将代谢物分组为相互之间有重叠的通路，将PLS和Group Lasso结合建立回归模型，并引入通路权重系数和惩罚因子，实现基于“分组稀疏”的多元回归。

本发明包括以下述步骤：

S1：对采集得到的质谱转化成mzML文件，对质谱文件进行中心化、去噪以及对齐处理，得到样本的代谢物的二维数据矩阵；采集实验样本，对实验样本处理得到健康志愿者样本和患病志愿者样本，对样本进行LC-MS/MS实验得到质谱。

S11：首先样本收集与检测：

所有志愿者(总样本量为N)在采血前一晚禁食，采血后进行结肠镜检查。根据活组织切片检查的结果将志愿者分为健康者(样本量为N₁)和疾病患者(样本量为N₂)两类。每个血液样本凝结45分钟后，以2000转/分钟速度离心10分钟，离心后其取上清液储存于-80℃冰箱中直至实验。

对待检测的代谢组学生物样本进行靶向LC-MS/MS实验,质谱实验在Agilent1260LC(Agilent Technologies,Santa Clara,CA)AB Sciex QTrap 5500MS(AB Sciex,Toronto,Canada)仪器上进行，在正、负模式下各检测一次，使用MultiQuant 2.1软件(ABSciex，Toronto，Canada)对谱峰进行积分处理。

S12：谱图预处理：

对采集得到的质谱先转化成mzML文件，然后用R读取文件，使用XCMS v3.0.2软件对质谱文件进行中心化、去噪、对齐等处理，得一个N个样本M个代谢物的二维数据矩阵X，矩阵的每一行代表一个样本信息，每一列代表一个代谢物信息，矩阵大小为N×M。

S2：对样本的代谢物的二维数据矩阵以及样本类别矢量矩阵，做中心化和UniteVarience归一化处理；

待处理的样本数据表示为X，样本类别矢量表示为Y，大小为N×1。分别对X和Y做中心化和Unite Varience归一化处理。

S3：获取代谢物-通路映射关系，并通路权重系数进行优化；

S31：获取代谢物-通路映射关系

京都基因与基因组百科全书(KEGG；http://www.genome.jp/kegg/)是一个可公开访问的数据库,其中包含了一系列人工绘制的代谢通路图，从该数据库中下载所研究数据集包含的M个代谢物参与的所有通路。

为使结果更加可靠和可解释，排除掉那些检测到的代谢物少于三种的代谢通路。

S32：通路权重系数ω_g优化：

(1)初始化ω_g：

(2)通路选择的经验分布：

初始化通路选择计数器c_g：c_g＝0,g＝1,2,…,G；

从n＝1到n＝N循环做：

a.对原始响应y随机置换得到一个新的响应y′，然后在数据(y′,X)上构建一个ogPLS模型M_n；

b.调谐λ使得仅有一个通路被选择，当且仅当

时，唯一的通路

被选择；

c.如果通路

被模型M_n选择，c_g＝c_g+1。

那么，通路选择的经验分布为：

(3)通路权重系数的调整：

a.计算通路选择的经验分布Π^ω和零分布

之间的距离：

b.∈≥0表示收敛精度，如果D(Π^null|∏^ω)<∈，停止通路权重系数调整步骤。

c.0<α<1表示学习率，从g＝1到g＝G，循环做：

(4)循环迭代：回到(2)继续循环。

S4：进行通路排序，调整惩罚因子并确定通路被选择的频率，利用通路的选择频率对通路进行排序。

S41：确定λ值：

a初始化：令λ＝1.0，K＝min(N,M)，β>0表示学习率(β>0)，S＝φ表示所选通路集；

b ogPLS建模：

用ogPLS方法计算

如果

则S＝{S,g}。

c调整惩罚因子λ：

如果|S|≠L(L为所要选择的通路数量)，则更新λ＝λ+β·(|S|-L)；

回到步骤b继续ogPLS建模知道满足|S|＝L。

S42：Bootstrapping建模

a.初始化通路选择计数器c_g＝0,g＝1,2,…,G

b.从n＝1到n＝1000循环做：

对数据进行50％的随机采样构建一个ogPLS模型M_n；

如果通路

被模型M_n选择，则c_g＝c_g+1。

S43：利用通路选择的频率对通路进行排序F^λ＝(c₁,c₂,…,c_G)/N，通路被选择的频率越高，该通路的排序越靠前。

所述步骤中的Overlapping-group PLS(ogPLS)模型的方法为：

u_k表示一个M×1的ogPLS模型的基向量；

表示通路特定基向量，满足

如果代谢物

是一个M×1的稀疏向量，它的大多数元素都为0，且这些0元素在优化过程中保持不变，为了突出这一点，这里给组

定义了一个收缩矩阵：

其中，

是组

中检测到的代谢物，i＝1,2,…,M是在整个数据集X中检测到的代谢物。然后我们可以定义一个收缩的通路特定基向量：

ogPLS的具体步骤如下：

(1)初始化通路特定基向量

(2)计算u_k：

从k＝1到K循环：

a从g＝1到G计算：

b计算

c如果‖u_k‖₂≠0，则计算：

t_k＝Xu_k

否则‖u_k‖₂＝0，则K＝k，回到(3)

(3)建立回归模型：

其中，U＝(u₁，u₂，…，u_K)，X表示N个样本M个代谢物的二维数据矩阵，y表示原始响应。

下面通过具体实例进行说明：

1.研究对象

本发明使用在印第安纳州大学医学院采集的158个血清样本，分别来自66个结直肠癌患者和92个健康志愿者，所有结直肠癌患者均是第一次被诊断出，采血之前未进行任何手术治疗、化学治疗和放射治疗。

2.试剂和仪器设备

实验试剂包括重水、三氯甲烷、甲醇、磷酸二氢钠、磷酸氢二钾、乙酸铵、醋酸。

实验耗材包括5mm核磁管、移液器、1000μL和200μL枪头、5mL EP管、记号笔、乳胶手套和口罩。

实验设备包括移液器、低温高速离心机、安捷伦1260LC和AB Sciex QTrap5500MS。

3.样本收集和预处理

采集66个结直肠癌患者和92个健康志愿者的血浆，采集到的血液静置45分钟后以2000转/分钟的转速离心10分钟，离心后其取上清液储存于-80℃冰箱中直至实验。

4.LC/MS谱测定

50μL血清样本与150μL甲醇涡旋混合，并静置于-20℃冰箱20分钟；取出样本以14000转/分钟的速度离心10分钟，取上清液；对于剩下的蛋白质小球，再次加入300μL甲醇、涡旋混合10分钟、14000转/分钟离心10分钟，之后取上清液；两轮中提取的上清液混合后进行干燥；干燥后的固体复溶于500μL的99.8％乙酸铵溶液(5mM乙酸铵，40％水/60％乙腈)+0.2％醋酸(5.13μM L-酪氨酸-13C2/22.54μM L-乳酸钠-13C3)，其中两种同位素标记的标准物被用于监控实验系统稳定性。LC-MS/MS实验在Agilent 1260LC(AgilentTechnologies,Santa Clara,CA)和AB Sciex QTrap 5500MS(AB Sciex,Toronto,Canada)仪器上进行。每个样本正、负模式下各检测一次。靶标质谱数据在质谱多反应监测(multiple reaction monitoring，MRM)模式下获取。在正、负模式下分别检测到59对和99个母离子-子离子对(总共158对)。最后，利用MultiQuant 2.1软件(AB Sciex，Toronto，Canada)对谱峰进行积分处理。最终113个代谢物被识别和量化，2019年3月16日发布的KEGG数据库共下载了81条人类(智人)代谢通路(包含1498个代谢物)。为使结果更加可靠和可解释，本发明排除了检测到的代谢物少于三种的代谢通路。经过数据预处理，从CRC数据集中只纳入30个代谢通路和81个代谢产物进行进一步分析。

5.数据处理

对CRC数据集的健康和癌症样本的数据矩阵X＝(x_ij)_158×81做中心化和UniteVarience归一化处理。用y＝(y_i)_158×1表示响应变量，根据样本所属的类别设为：

6.获取代谢物-通路映射关系

(1)2013年3月16日，从KEGG数据库下载了81条人类代谢物通路信息(包含1498个代谢物)，

(2)为使结果更加可靠和可解释，本发明排除了那些包含的检测到的代谢物少于三种的代谢通路。经过数据预处理，从CRC数据集中只纳入30个代谢通路和81个代谢产物进行进一步分析。

7.通路权重系数优化

用M^h～c表示健康和癌症样本数据(y,X)的ogPLS模型，使用所述的通路权重优化算法对M^h～c的通路权重系数进行优化和设置。

8.通路排序

由于在具体实施时已验证分组排序算法对参数λ是鲁棒的，此处选择λ＝20来对各个通路进行排序。

9.方法验证与对比

(1)背景技术中的PIP方法和Global-Test方法是根据通路重要程度或者PIP值来对通路进行进行排序的。图1列出了本发明和上述两种方法用于健康类和癌症类通路分析的结果。

通路通过颜色的深浅进行了排序，通路的排序越靠前，对应的颜色越深。通路的独立性I_g被定义为：

其中，

是代谢物m参与的通路数，

是通路g的代谢物数，0<I_g≤1和I_g＝1表示与其他通路不重叠的通路。在统计意义上通路的独立性越高，通路分析结果的可靠性越高。

从图1中可以看出通过ogPLS选择出的的排序靠前的通路，它们的独立性较高。PIP和Global-Test的结果更加一致，而ogPLS的结果与其他两种方法有所不同。例如，在PIP和Global-Test方法中，D-谷氨酰胺和D-谷氨酸代谢通路和精氨酸生物合成通路的排名较前，而在ogPLS中则较后。进一步的独立性分析表明，两种通路的独立性分别为0.24和0.16。这说明这两条通路与其他通路有高度重叠，且这两条通路不可靠。

(2)两种度量，定义排序后前r个通路的累积代谢物覆盖率#(r)和累积代谢物重叠率O(r)来描述三种方法的排序行为：

其中，r是排序的顺序，

是排序为i的通路，#(r)是在排名前r的通路中代谢物的覆盖率(从排名1到排名r)，O(r)是在排名前r的通路中代谢物的重叠率。

从图2可以看出，横坐标为选择的通路，纵坐标为选择的代谢物，与PIP和Global-Test相比，ogPLS倾向于选择代谢物覆盖率更高的通路。

从图3可以看出，横坐标为选择的通路，纵坐标为选择的代谢物的重叠率，对于前几个显著通路，PIP和Global-Test方法比ogPLS方法更容易选择代谢物重叠率高的通路。例如，PIP和Global-Test方法的Ο(r＝5)比ogPLS方法更高，在ogPLS方法中Group Lasso惩罚可以有效减少重叠通路的竞争力。

上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。

Claims (2)

1.一种质谱代谢组学的通路分析建模方法，其特征在于，包括如下步骤：

S1：采集正常生物和患病生物的生物样本，对采集样本得到的质谱转化成mzML文件，对质谱文件进行中心化、去噪以及对齐处理，得到样本的代谢物的二维数据矩阵；

S2：对样本的代谢物的二维数据矩阵以及样本类别矢量矩阵，做中心化和单位方差归一化处理；

S3：获取代谢物-通路映射关系，并通路权重系数进行优化；

S4：进行通路排序，调整惩罚因子并确定通路被选择的频率，利用通路的选择频率对通路进行排序；

所述步骤S3对代谢物-通路映射关系的通路权重系数进行优化，具体包括：

S31：初始化通路权重系数ω_g：

其中，

表示的是该通路组内包括的代谢物的数目；

S32：初始化通路选择计数器c_g：c_g＝0,g＝1,2,…,G，通路选择的经验分布为：

S33：通路权重系数的调整：计算通路选择的经验分布Π^ω和零分布

之间的距离：

若D(Π^null|∏^ω)＜∈，停止优化过程，其中∈≥0且∈表示收敛精度；

从g＝1到g＝G，循环

S34：循环迭代：转步骤S32；

所述步骤S4进行通路排序，调整惩罚因子并确定通路被选择的频率，利用通路的选择频率对通路进行排序，具体包括：

S41：确定惩罚因子λ值：

初始化：令λ＝1.0，K＝min(N,M)，β＞0且β表示学习率，S＝φ且s表示所选通路集；

ogPLS建模：用ogPLS方法计算

若

则S＝{S,g}；

调整惩罚因子λ：如果|S|≠L，L为所要选择的通路数量，则更新λ＝λ+β·(|S|-L)；回到ogPLS建模步骤继续ogPLS建模直至满足|S|＝L；

S42：确定通路被选择的频率：

自助法建模，初始化通路选择计数器c_g＝0,g＝1,2,…,G；从n＝1到n＝1000循环：对数据进行50％的随机采样构建一个ogPLS模型M_n；若通路

被模型M_n选择，则c_g＝c_g+1；

S43：利用通路的选择频率对通路进行排序：

利用通路选择的频率对通路进行排序F^λ＝(c₁,c₂,…,c_G)/T，通路被选择的频率越高，该通路的排序越靠前。

2.根据权利要求1所述的一种质谱代谢组学的通路分析建模方法，其特征在于，构建所述ogPLS模型的方法，具体包括：

S421：初始化通路特定基向量

K＝min(N,M)，N为样本数量，M为N个样本中代谢物类别的数量；

S422：计算ogPLS模型的基向量u_k：从k＝1到K循环：

从g＝1到G计算

其中

Z^(g)表示组

的收缩矩阵，

是组

中检测到的代谢物，i＝1,2,…,M是在整个数据集中检测到的代谢物；

计算

若||u_k||₂≠0，则计算：

t_k＝Xu_k，

否则||u_k||₂＝0，则K＝k；

S423：建立回归模型：

其中X表示N个样本M个代谢物的二维数据矩阵，y表示原始响应，U＝(u₁，u₂，…，u_K)。

Images