LC-MS下机数据的处理方法及处理装置与流程

lc-ms下机数据的处理方法及处理装置
技术领域
1.本发明涉及分析化学领域的数据处理技术，特别是涉及一种对样本经lc-ms(液相色谱质谱联用仪)分析得到的下机数据进行处理获得代谢组数据集的方法及装置。


背景技术：

2.代谢组学(metabolomics)是涉及代谢产物的化学过程的科学研究。具体而言，代谢组学是对特定的细胞过程遗留下的特殊化学指纹，即对它们的小分子代谢产物的整体研究。代谢组(metabolome)定义为在一个生物细胞、组织、器官或生物体中所有的代谢产物(化合物)的集合。
3.液相色谱质谱联用仪(lc-ms)因其卓越的灵敏度和广泛的覆盖范围成为代谢组学研究中常用的技术平台，通过它能快速而极为准确地鉴定生物小分子化合物。
4.lc-ms的高灵敏度导致对下机数据的处理相当复杂，主要的困难或陷阱包含以下几个方面：峰识别、峰积分(峰积分的目标是获得这些色谱峰的保留时间、高度和面积)、保留时间校正、加合离子去除、去除噪音和批次效应去除。目前的代谢组质谱数据处理方式多种多样，然而，对于以上提及困难的解决方法没有一致的标准，这就导致不同的处理方法获得的代谢组数据不同，降低了可比性与可信度。
5.色谱峰保留时间应具有高度可重复性，以便在样品之间准确选择峰。然而，给定代谢物(指同一化合物)的峰保留时间通常可能在相同检测条件不同检测批次中略有变化，可能是由色谱柱老化、样品过载、流动相ph值不稳定以及色谱柱温度和压力变化引起的。而且，即便在同一个样本中，不同的化合物对上述影响因素的反应不同，从而导致峰保留时间漂移的程度不同。所有这些在分析过程中都很难严格控制，尤其是在涉及数百或数千个样本的大规模代谢组学研究中。因此，保留时间的对齐是数据预处理中的一个重要步骤，以确保来自不同样品和不同检测批次中相同代谢物的峰对齐到相同的保留时间从而指向相同的代谢物或成分(化合物)。
6.lc-ms技术的高通量加上对大型实验的需求导致数据预处理，即跨样品的代谢物量化，成为主要瓶颈。目前市面上也有许多商用或非商用的数据处理软件，如compound discoverer,msdial等。但是，现有软件最大的限制是单次处理的最大样本量受限，单次只能处理数十个样本，且极其耗费计算资源。对于数百甚至数千样本量的大型项目，必须分批次处理质谱数据，耗费极大的人力与时间成本。


技术实现要素：

7.鉴于以上所述现有技术的缺点，本发明的目的在于提供一种高通量的代谢组lc-ms下机数据的处理方法及装置。
8.为实现上述目的及其他相关目的，本发明第一方面提供一种液相色谱质谱联用仪数据的处理方法，所述方法至少包括以下步骤：
9.s1：读取lc-ms的下机数据、制作上机序列集(正模式b1和负模式b2)及参数文件，
转换下机数据的一级全扫描图谱信息对每个数据进行解卷积并扣除背景基线，获得正模式数据集(a1)和负模式数据集(a2)。
10.s2：针对步骤s1获得的正模式数据集(a1)和负模式数据集(a2)，用至少一种方式校正每个样本间的保留时间漂移，再将每个样本内的色谱峰与参考数据库内的参考化合物特征进行匹配，将匹配上的化合物的色谱峰归为正模式靶向数据集(c1、c2、c3)和负模式靶向数据集(c4、c5、c6)；另将剩余未能匹配成功的色谱峰进行样本间的峰匹配，归为非靶向数据集(d1-d6)。
11.s3：对步骤s2不同方式下获得的正模式靶向数据集和负模式靶向数据集(c1和c4，c2和c5，c3和c6)的化合物分别进行匹配审查，按照只保留主要加合离子峰的原则扣除次要加合离子峰，并进行冲突化合物审查，对冲突化合物合并或过滤，针对不同方式分别建立靶向合并数据集(e1、e2和e3)。
12.s4：当进行了多种方式的保留时间校正时，对步骤s3获得的靶向合并数据集(e1、e2和e3)中不同方式下同一色谱峰在不同样本内的强度进行比较，保留最高样本检出率对应方式下该色谱峰的强度和化合物信息，建立平台靶向数据集(f)。
13.s5：将步骤s2和s3中获得的其中一种方式下产生的非靶向数据集和靶向合并数据集合并，计算两两色谱峰之间的相关性，保留低于相关性阈值的非靶向色谱峰，与步骤s4中得到的平台靶向数据集(f)合并，获得去冗余数据集(g)。
14.s6：对步骤s5获得的去冗余数据(g)进行缺失值填充处理后进行批次效应校正，再进行数据归一化，然后进行标准化转换，得到鉴定出化合物的、去除了批次效应的代谢组数据集(h)。
15.优选的是，在步骤s1中，批量将下机数据raw文件(文件后缀为“.raw”)转换为binary文件(文件后缀为“.bin”)，再将binary文件再转化为hdf5文件(文件后缀为“.h5”)。
16.优选的是，步骤s1中，以用户自定义的适配文件中的参数值，将下机数据解卷积后，去除数据的噪音和背景后，按响应阈值筛选每个时间窗格中的色谱峰。
17.优选的是，在步骤s2中，所述方式包括：
18.方式一：参照样本的校正与化合物识别；
19.方式二：内部校准的校正与化合物识别；
20.方式三：不进行保留时间校正的化合物识别。
21.所述方式一是按设定的参数调整参数对步骤s1的数据集(a1或a2)中每个样本的保留时间漂移进行校正，并依据参考数据库中已储存的参考化合物不同加合离子的保留时间和质荷比进行化合物识别。方式一设定的参数为：参照样本中各校准化合物编号，化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
22.方式一保留时间漂移校正方法为：根据检测序列中邻近的两个参照样本内校准化合物的保留时间逐一对在这两个参照样本之间的其他样本的保留时间按区段进行校正，即将整个保留时间轴划分为多个时间段，在每个时间段内，计算参照样本中各校准化合物的实际保留时间与预期保留时间之间的差值，再使用局部线性校正的方法对其他生物样本的色谱峰的保留时间进行线性缩放。
23.具体的方式一保留时间漂移校正方法例如：在整个检测序列中，生物样本的保留时间按最邻近的头尾两个参照样本的保留时间进行校正，计算每个校准化合物在头尾两个
参照样本内的保留时间差(rt
1-rt2＝δrt)，将头尾两个参照样本中间的多生物样本对应的保留时间进行线性缩放至该校准化合物在头尾两个参照样本内的保留时间的均值(rt1+0.5δrt)。
24.所述方式二是按设定的参数调整参数对步骤s1的数据集(a1或a2)中每个样本的保留时间漂移进行校正，并依据参考数据库中已储存的各参考化合物的不同加合离子的保留时间和质荷比进行化合物识别。方式二设定的参数为：化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
25.方式二保留时间漂移校正的方法为：不采用参照样本而是根据样本检出率和平均强度皆为最高的原则，在各保留时间轴上各时间段中筛选色谱峰作为校正其他生物样本的校正峰，再根据校正峰的保留时间轴逐一对其他样本的时间轴进行校正，用筛选出的校正峰的保留时间与质荷比信息匹配参考化合物，并将各匹配上参考化合物的校正峰的保留时间缺失值进行填充，再将各校正峰的实际保留时间替换为预期保留时间，最后对其他色谱峰的保留时间分段进行线性缩放。
26.具体的方式二保留时间漂移校正的方法例如：(a)筛选校正峰：按筛选校正峰的最低样本检出率(80％)，在各保留时间段(设定为0-1,1-2,2-5,5-8,8-10,10-11,11-13,13-16分钟)中用每个时间段对应的最低峰强度(分别为2+e07,5+e06,8+e06,2+e06,9+e05,5+e07,7+e07,9+e07)从参考化合物的色谱峰特征中筛选出备选校正峰；依据已存储于参考数据库中各参考化合物的不同加合离子的保留时间和质荷比对备选的校正峰进行化合物识别，筛选校正峰的化合物识别保留时间窗口为0.08分钟，筛选校正峰的化合物识别质量精度偏差为5ppm，筛选的备校正峰匹配上的参考化合物的预期保留时间即为校正峰的预期保留时间。对各校正峰保留时间的缺失值进行填充，填充值为该校正峰在所有检出样本中保留时间的中位值。(b)经过缺失值填充，每一个校正峰在所有样本中都有保留时间的值。将所有样本中的各校正峰一一对应后，将各校正峰的保留时间调整为预期保留时间，其余峰随之进行线性缩放完成校正。
27.所述方式三是指不对保留时间漂移进行校正，直接按设定的参数依据参考数据库中已储存的各参考化合物的不同加合离子的保留时间和质荷比对步骤s1的数据集(a1或a2)中每个样本进行化合物注释，产生靶向数据集(c3或c6)。设定的参数为：化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
28.优选的是，步骤s2中将色谱峰与参考化合物特征进行匹配的方法，包括以下步骤：1)对每个样本里的每个色谱峰，都与参考数据库中参考化合物不同加合离子的参考保留时间和参考质荷比进行匹配，其中，允许的质量精度偏差和保留时间窗口由用户自定义；2)计算每个加合离子的样本检出数、平均保留时间、平均质荷比和平均离子强度，并与参考保留时间、参考质荷比和参考离子强度进行对比。此步骤的目的是，一次过地对大批量的样本进行化合物鉴定，并返回代谢物谱易于判断的概括性特征，利于进行快速且精准的化合物鉴定。
29.优选的是，步骤s2中对剩余色谱峰进行样本间的峰匹配的方法，是将所有色谱峰在样本间进行两两对比，将在设定的质量精度偏差(如5ppm)以内和保留时间偏差(如0.03分钟)以内的峰作为同一化合物在不同样本间检测出的峰，归为同一个非靶向色谱峰。
30.优选的是，步骤s3中审查化合物匹配结果的方法，按参考化合物主要加合离子和
各次要加合离子之间的关系，与样本中与之匹配的色谱峰组合间的关系进行比较。如果关系相似，认为是正确的化合物匹配结果，按照只保留主要加合离子峰的原则扣除次要加合离子峰；如果关系不相似，则认为是错误的匹配结果，此色谱峰组合归入未能匹配成功的非靶向数据集。关系相似至少需要满足以下条件：1)匹配到了主要加合离子，该离子峰在所有加合离子峰中强度最高；2)匹配到了第二加合离子或者第三加合离子，该离子峰的保留时间与主要加合离子峰的保留时间偏差在0.01分钟以内，且该离子峰的强度与主要加合离子峰强度的比值和参考化合物该次要加合离子峰与参考化合物主要加合离子峰强度的比值近似(两个比值偏差在0.1以内)；3)主要加合离子的检出样本数要不低于次要加合模式的检出样本数。
31.优选的是，步骤s3中过滤冲突化合物时，在化合物匹配结果正确的色谱峰中，比较具有相同保留时间和质荷比的色谱峰，筛选并合并满足特征的冲突化合物，不满足特征的色谱峰被过滤。保留的冲突化合物以化合物名称、保留时间偏差和质量精度偏差以“|”符号隔开合并表示。保留的冲突化合物至少需要满足以下特征：1)质量精度偏差小于设定值(如5ppm)，2)冲突化合物的检出样本数不低于此冲突组中最高检出样本数的设定百分比(如90％)。
32.优选的是，在步骤s5中，对步骤s2中采用参照样本的校正方式校正得到的非靶向数据集(d1和d4)，步骤s3中得到的靶向合并数据集(e1)中两两色谱峰的峰强度进行相关性比对，小于设定的相关性系数阈值(例如0.9)的非靶向色谱峰保留并入步骤s4中得到的平台靶向数据集(f)。
33.优选的是，步骤s6中，进行缺失值填充处理的方法采用python语言的pandas包进行，按实验设计进行调整：1)针对样本检出率低于设定百分数(如50％)的色谱峰，采用仪器检出下限的定值插补的方式进行处理，对样本检出率高于设定百分数(如50％)的色谱峰，采用插补该色谱峰中位值的方式进行处理；或者2)按实验设计分组，分别对每一个分组按每组的仪器检出下限进行定值插补；或者3)采用定值插补的方式进行处理，插补的值为仪器检出下限的五分之一。
34.步骤s6中，进行批次效应校正的方法，按实验设计和效果进行调整，包括1)normae方法：采用python语言的normae脚本进行，输入步骤s1中的上机序列集的测量批次和测量顺序，使用默认参数，对色谱峰在不同批次的样本之间的强度进行校正，去除仪器响应波动的影响；或者2)总和方法：计算每个样本中所有色谱峰的强度总和，再计算所有样本强度总和的平均值和标准差，将每个样本强度总和减去该平均值后除以该标准差，获得缩放系数，最后将每个样本中的色谱峰强度都除以该样本的缩放系数；或者3)中位值方法：计算每个样本中所有色谱峰的强度中位值，再计算所有样本强度中位值的平均值和标准差，将每个样本强度中位值减去该平均值后除以该标准差，获得缩放系数，最后将每个样本中的色谱峰强度都除以该样本的缩放系数；或者4)分位数方法：将每个样本中的色谱峰强度按从高到低排序，计算排序第一的色谱峰强度的均值，并将所有样本中排序第一的色谱峰强度替换为该均值，同理对排序第二、三
……
的色谱峰强度进行替换。
35.步骤s6中，数据归一化的方法，按实验设计和效果进行调整，包括1)均值归一化：将每个样本中的色谱峰强度减去样本色谱峰强度的均值；或者2)z-score标准化：将每个样本中的色谱峰强度减去样本色谱峰强度的均值后，除以样本色谱峰强度的标准差；或者3)
离差标准化：将每个样本中的色谱峰强度减去样本色谱峰强度的均值后，除以样本色谱峰强度的强度范围(最大值减最小值)。
36.步骤s6中，标准化转换的方法选自log2标准化方法。
37.本发明第二方面提供一种代谢组质谱数据的分析方法，对本发明第一方面lc-ms下机数据的处理方法得到的代谢组数据集进行批次效应校正效果分析，包括以下步骤：
38.1)使用批次效应校正前的代谢组质谱数据，使用r语言的prcomp()函数，观察参照样本与待测样本代谢组在pc1和pc2上的空间位置；
39.2)使用批次效应校正后的代谢组数据，按照步骤1)同样的方法进行pca分析，对比两步骤产生的pca图。
40.本发明再一方面提供一种生物代谢组的液相色谱质谱联用仪数据的处理装置，执行前述对lc-ms下机数据处理方法的操作，所述装置包括：
41.获取模块，用于获取下机数据、上机序列集及参数文件，并对下机数据的一级全扫描图谱信息进行转换得到数据集；
42.保留时间校正及化合物识别模块，设定至少一种方式，用于校正样本间的保留时间漂移后，再将每个样本内的色谱峰与参考数据库内的参考化合物特征进行匹配，将匹配上的化合物的色谱峰归为靶向数据集；并将剩余未能匹配成功的色谱峰进行样本间的峰匹配，获得非靶向数据集；
43.化合物匹配及冲突化合物审查模块，用于审查每个样本内的色谱峰与数据库内的参考化合物特征进行匹配的关系是否正确，过滤或合并冲突化合物的色谱峰；
44.保留时间校正审查模块，用于比较不同保留时间校正及化合物识别方式下匹配上相同化合物的色谱峰；
45.冗余峰过滤模块，用于过滤相关性高的色谱峰；
46.批次效应校正模块，用于对所获得的去冗余数据进行缺失值处理，然后根据参数文件中的批次信息对缺失值处理后的数据进行批次效应校正，然后进行数据归一化，最后进行标准化转换，得到代谢组数据集。
47.本发明又一方面提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现前述对lc-ms下机数据的处理方法和前述代谢组质谱数据的分析方法。
48.如上所述，本发明的代谢组lc-ms下机数据处理方法及装置，具有以下有益效果：
49.本发明的代谢组lc-ms数据处理方法及装置针对不同化合物保留时间漂移程度不同的情况，设计有三种方式对保留时间进行校正并进而进行化合物鉴定。然后，针对鉴定出的化合物，审核其在不同方式下的样本检出率和强度并进行合并。这是因为大规模代谢组学研究的保留时间漂移范围可能很大，单一方式的保留时间校正可能没有办法完美校正所有样本的偏移，导致只能在部分样本中鉴定出的化合物。本发明同时采用三种保留时间校正与化合物识别方法，可以起到互相补充、互为印证的作用。与传统的单一保留时间校正与化合物识别方法相比，在化合物鉴定上具有显著的优越性。
50.本发明的代谢组lc-ms下机数据处理流程，其最大特征在于能够一次输入数百甚至数千个样本的lc-ms检测下机数据，获得一个鉴定出尽可能多化合物的、去除批次效应的、去冗余的、可直接用于后续统计学分析的代谢组数据集(其为由样本名称、色谱峰编号、
化合物名称和峰强度构成的矩阵)。本发明针对lc-ms数据处理的困难进行各个击破并实现全局最优，显著降低了单个样本的处理时间和人力资源，在保证准确性的同时，使得该方法相比传统方法具有更好的适应性和应用范围。
附图说明
51.图1为本发明一实施例中代谢组lc-ms下机数据处理方法的流程图。
52.图2为本发明一实施例中代谢组lc-ms下机数据处理装置的示意图。
53.图3为本发明一实施例中代谢组lc-ms下机数据处理方法流程中涉及的具体数据集展示流程图。
54.图4为参照样本的一级全扫描色谱图(a幅为色谱图，b幅为质谱图)，显示了一个参照样本的下机数据，以总离子流图形式展示。
55.图5所示为某一示例标准化合物的散点图，显示采用参照样本的校正方法时，所有参照样本中该校准化合物按上机顺序排列的保留时间，以展示保留时间漂移。
56.图6所示为某一示例化合物的散点图，显示采用内部校准的校正方法时，所有样本中该化合物按上机顺序排列的保留时间，以展示保留时间漂移。
57.图7为本发明一实施例中保留时间校正审查的示例散点图。
58.图8为本发明一实施例中批次效应校正分析的示例散点图。
具体实施方式
59.本技术涉及一种对样本经lc-ms(液相色谱质谱联用仪)分析得到的下机数据进行处理获得代谢组数据集的方法及装置。
60.液相色谱(lc)，可根据化合物极性在时间维度(x轴)上分离样品中复杂混合物的组分，各组分(各化合物)流经质谱检测器时响应的连续信号在色谱图上显示为特征峰，质谱仪(ms)测量各特征峰对应化合物的质荷比(y轴)及特征峰信号强度(z轴)。由保留时间(x轴)、质荷比(y轴)和信号强度(z轴)组合定义的特征峰被称作色谱峰，全部色谱峰的集合生成一个能包含有样品中各组分信息的三维矩阵。本技术中的下机数据，是指样本通过例如acquity uplc i-class液相色谱系统(沃特世)和例如q-exactive质谱系统(赛默飞世尔科技)联用的lc-ms检测产生的三维矩阵，以q-exactive质谱系统默认的raw数据格式呈现。
61.参数文件指包含下述步骤中所有参数设定值的r语言脚本。
62.代谢组数据集，是指下机数据经本技术处理后最终产出的与诸多代谢物质(化合物)相关的信息，其为由色谱峰编号、化合物名称和在样本中检出的峰强度构成的数据集。
63.被分离样品组分从进样开始到柱后出现该组分浓度极大值时的时间，也即从进样开始到出现某组分色谱峰的顶点时为止所经历的时间，称为此组分的保留时间，用rt表示，常以分钟(min)为时间单位。
64.参考数据库：包含多个参考化合物在与样本相同色谱、质谱、离子源条件下产生的参考化合物特征。
65.参考化合物：是包含精确已知浓度的物质的化学标准品，参考样本中的校准化合物属于参考化合物的一种。
66.总离子流色谱图(tic)是通过将属于同一扫描的所有质谱峰的强度相加而创建的
色谱图(图4就是tic)。
67.参考化合物特征为某个特定化合物在相同的lc-ms实验环境中产生的色谱峰，包含了常见的电喷雾离子源(esi)加合离子类型的母离子的碎片信息，除了正离子模式下常见的[m+h]
+
和负离子模式下常见的[m-h]-之外，还包含了[m+h-h2o]
+
，[m+na]
+
，[m+nh4]
+
，[m+h-nh3]
+
，[2m+h]
+
，[m-h-h2o]-，[2m-h]-和[m-2h]
2-等形式的加合离子数据。
[0068]
正模式与负模式指质谱系统中电喷雾电离(electrospray ionization,esi)的两种离子化模式。esi将液态样品转为气相，并将不带电的分子变为带正电或者负电的分子离子。
[0069]
特定的物质通过esi后会根据其本身的物理化学性质被电离成一种到多种不同的加合离子，其中占比最多的是主要加合离子，占比第二多的是次要加合离子。
[0070]
同一个色谱峰可能会被匹配上不同的化合物，这些化合物就被称为冲突化合物。
[0071]
质荷比指带电离子的质量与所带电荷之比值，以m/z表示。质量精度偏差是描述实测离子质荷比与准确离子质荷比之间差异的指标(以ppm表示)。
[0072]
主成分分析(principle component analysis,pca)是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些不相关变量称为主成分(principle component,pc)。pc1为第一主成分，是数据集中对方差贡献最大的特征，同理，pc2为第二主成分。
[0073]
以下通过特定的具体实施例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。
[0074]
需要说明的是，本实施例中所提供的图示仅以示意方式说明本发明的基本构想，虽图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。
[0075]
以一个用于人类血浆代谢组学研究的数据为例，说明本发明提供的一种代谢组lc-ms下机数据的处理系统。
[0076]
本实施例针对3300个样本进行lc-ms检测，样本包括生物样本(血浆)和参照样本。参照样本由数个化学标准品(包含校准化合物)混合而成，这些标准品的参考信息存在于参考数据库中。参照样本在每批血浆(例如，每批检测的生物样本数量为40个)检测开始、中间(例如，每间隔10个生物样本)和结束时上机检测。每个标准品在参照样本下机数据的色谱峰可用于评估该峰对称性，还可以评估色谱稳定性和质量精度，还可以用来进行保留时间校正。
[0077]
本实施例中，所有样本使用acquity uplc i-class液相色谱系统(沃特世)和q-exactive质谱系统(赛默飞世尔科技)采集信息，其中质谱仪为配有电喷雾电离源(esi)的静电场轨道阱(orbitrap)质谱。具体操作为：将待检批次的各样本玻璃进样小瓶中的3μl上清液注入经1.8μm，2.1mm
×
30mm的acquity uplc hss t3色谱柱分离(流动相a为含0.1％甲酸水溶液，流动相b为含0.1％甲酸的乙腈溶液，均为色谱纯)，随即由联用的q-exactive质谱系统进行检测，获得该批次各检测样品的下机数据，均为q-exactive质谱系
统默认的raw数据格式。
[0078]
下面参照图1至图3，说明针对lc-ms下机数据的处理方法和处理装置。
[0079]
对lc-ms检测的下机数据的处理具体包括如下步骤，是由图2所示的代谢组质谱数据处理装置的各个模块来执行的。本发明一实施例中数据处理方法流程(参见图1)中涉及的具体数据集如图3所示。
[0080]
【步骤s1】读取lc-ms的下机数据、制作上机序列集及参数文件，转换下机数据的一级全扫描图谱信息对每个数据进行解卷积并扣除背景基线，获得正模式数据集a1和负模式数据集a2。
[0081]
1)首先读取该批次样本检测产生的所有待处理的下机数据raw文件的一级全扫描图谱信息(参见图4所示)，批量将raw文件(文件后缀为“.raw”)转换为binary文件(文件后缀为“.bin”)，对每个数据进行解卷积并扣除背景基线并转化为hdf5文件(文件后缀为“.h5”)。正模式所有hdf5文件为数据集a1，负模式所有hdf5文件为数据集a2。
[0082]
参照样本的一级全扫描色谱图参见图4所示，该参照样本的a1数据集如下表：
[0083]
peak_idrtmzintensity14.280215102.0462352698.424.3136504102.0462240689.330.512868102.97001744925.0省略4048个
………………
405217.47795993.747190258.9
[0084]
注：表中，rt表示峰保留时间(x轴值)，mz表示质荷比(y轴值)，intensity表示峰强度(即信号强度，z轴值)。以下同。
[0085]
该参照样本的a2数据集如下表：
[0086]
peak_idrtmzintensity10.625573100.93361060720.323.237853100.93341880686.730.626843101.9338133571.5省略3158个
………………
31629.458827867.2982002522
[0087]
本实施例中某一生物样本的数据集a1和数据集a2如下表所示。
[0088]
peak_idrtmzintensity14.746851128.468595820.328.074239128.4695105615.538.496752129.0016132886.4省略12813个
………………
1281711.156785794.5988146583.4
[0089]
2)接着制作所有样本的正模式上机序列集b1和负模式上机序列集b2。举例来说，一个样本分别检测了正模式和负模式，则需要分别制作正模式和负模式两个上机序列表，本实施例3300个样本的示例数据集b1和b2见表1和表2(表中省略了中间上机样本3294个)。
上机序列表包含以下五个变量：
[0090]
file_name：该样本对应的hdf5文件名。
[0091]
group：样本对应的分组，生物样本标为“us”，参照样本标为“qm”。
[0092]
msrt_day：测量批次，按测量开始日期从早到晚排序，并从1开始编号，每个批次内的编号一致。
[0093]
msrt_order：测量顺序，按测量先后顺序从前到后排序，并从1开始顺序编号至所有测量批次的最后一个样本，每个批次内的编号不一致。
[0094]
sample_id：样本编号，一个样本在同一平台内的样本编号唯一，且在上机序列表中分别对应上正模式的hdf5文件名和负模式的hdf5文件名。
[0095]
表1：正模式上机序列表b1
[0096][0097][0098]
表2：负模式上机序列表b2
[0099]
file_namegroupmsrt_daymsrt_ordersample_idsmp001_pnqm11smp0001smp002_pnus12smp0002smp003_pnus13smp0003
……
(省略7个)us1
…………
smp0011_pnqm111smp0011smp0012_pnus112smp0012
……
(省略3287个)
……………………
smp3300_pnqm743300smp3300
[0100]
3)制作参数文件
[0101]
参数文件为r语言脚本，内含后续所有步骤所需的参数，用户可根据实际需求调整参数。以下步骤中列举的本实施例中使用到的参数设定值都可通过此文件进行调整。参数包括：
[0102]
参照样本中各校准化合物的编号：各校准化合物在参考数据库中的化合物编号。
[0103]
化合物识别的保留时间窗口：进行化合物识别时，允许的实测保留时间和参考数据库中该化合物的保留时间之间的最大偏差范围。
[0104]
化合物识别的质量精度偏差：进行化合物识别时，允许的实测质荷比和参考数据库中该化合物的质荷比之间的最大偏差范围。
[0105]
筛选校正峰的最低样本检出率：筛选校正峰时允许的最低样本检出率。
[0106]
筛选校正峰的化合物识别保留时间窗口：筛选校正峰的化合物时，进行该化合物识别允许的实测保留时间和参考数据库中参考化合物的保留时间之间的最大偏差范围。
[0107]
筛选校正峰的化合物识别质量精度偏差：筛选校正峰时，进行化合物识别允许的实测质荷比和参考数据库中参考化合物的质荷比之间的最大偏差范围。
[0108]
筛选校正峰的保留时间段：分保留时间段筛选校正峰。
[0109]
筛选校正峰的保留时间段对应最低峰强度：每个保留时间段内筛选校正峰时，允许的最低峰强度。
[0110]
【步骤s2】针对步骤s1获得的数据集a1和数据集a2，校正每个样本间的保留时间漂移，再将每个样本内的色谱峰与参考数据库内的参考化合物特征进行匹配，将匹配上的化合物的色谱峰归为靶向数据集(c1-c6)，并将剩余未能匹配成功的色谱峰进行样本间的峰匹配，归为非靶向数据集(d1-d6)。
[0111]
校正和识别方式至少包括以下三种：
[0112]
方式一：参照样本的校正与化合物识别
[0113]
作为步骤s2的一种方式，可以将前述数据集a1和b1输入保留时间校正及化合物识别模块，该模块中加载有参考数据库，参考化合物的信息(如表3所示)已存储于参考数据库中，其中参照样本的校准化合物给定化合物编号，按设定的参数调整参数对每个样本的保留时间漂移进行校正，并依据各参考化合物的不同加合离子的保留时间和质荷比进行化合物识别。设定的参数为：参照样本中各校准化合物编号，化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
[0114]
1)保留时间漂移校正的具体操作：在整个检测序列中，生物样本的保留时间按最邻近的头尾两个参照样本的保留时间进行校正。参照样本是同一样品的技术重复，因此所有参照样本中都有这十个校准化合物一一对应的色谱峰。一个校准化合物在头尾两个参照样本间的保留时间漂移可以反映中间十个生物学样本在此校准化合物保留时间附近色谱峰的保留时间漂移程度。假设保留时间的漂移是线性匀速的。因此，头尾两个参照样本内校准化合物一一对应后，计算每个校准化合物在头尾两个参照样本内的保留时间差(rt
obs-rt
ref
＝δrt)，将两个参照样本中间的十个生物样本对应的保留时间(rt
obs
+0.09δrt，rt
obs
+0.18δrt，
…
,rt
obs
+0.91δrt)进行线性缩放至该校准化合物在两个参照样本内的保留时间的均值(rt
obs
+0.5δrt)。
[0115]
本方式某一参照样本中各校准化合物的示例总离子流图见图4(图4为十种校准化合物的混和样本的正模式总离子流色谱图tic)，各校准化合物的色谱峰已用箭头标出，a幅色谱峰上数值表示保留时间，b幅色谱峰上数值表示质荷比。
[0116]
本方式中采用的校准化合物的预期保留时间与实测保留时间见图5(以其中一个校准化合物rfc0001为例)，图中横坐标表示该校准化合物上机检测顺序，纵坐标表示保留时间，参照样本中校准化合物的实测保留时间以点示意，预期保留时间以虚线表示。该图显示了该校准化合物rfc0001的保留时间漂移情况。
[0117]
2)对采用参照样本的校正方式进行保留时间校正后的a1和b1的数据进行化合物
识别，即与参考数据库中参考化合物的色谱峰特征进行匹配，产生正模式靶向数据集c1和正模式非靶向数据集d1。正模式下的靶向数据集c1中包含数据集a1中所有匹配上化合物的色谱峰(亦为离子峰)。将剩余未能匹配成功的色谱峰进行样本间的峰匹配，归为非靶向数据集d1，包含数据集a1中剩余的所有未匹配上化合物的色谱峰。下表展示一个生物样本经该步处理后的正模式靶向数据集c1，“compound_id”列显示了校准化合物的标号，“adduct”列显示加合方式；剩余的所有未匹配上化合物的色谱峰归集到正模式非靶向数据集d1中，d1表格中没有“compound_id”和“adduct”列。
[0118]
peak_idrtmzintensitycompound_idadduct1114.9890123.04471546892.6rfc0001[m+2h-hcooh]2+1124.9890273.0757371254.1rfc0001[m+h-h2o]+省略2416个
…………………………
1146811.1567794.5988146583.4rfc1963[m+h]+
[0119]
按照同样的方式，由前述数据集a2和b2，产生负模式靶向数据集c4和负模式非靶向数据集d4。
[0120]
样本间峰匹配的方法为，所有色谱峰在样本间进行两两对比，质量精度偏差5ppm以内和保留时间偏差为0.03分钟以内的峰，可以认为是同一个物质在不同样本间检测出的峰，归为同一个非靶向色谱峰。
[0121]
方式二：内部校准的校正与化合物识别
[0122]
作为步骤s2的另一种方式，可以将前述数据集a1和b1输入保留时间校正及化合物识别模块，该模块中加载有参考数据库，参考化合物的信息(如表3所示)已存储于参考数据库中，按设定的参数调整参数对每个样本的保留时间漂移进行校正，并依据各参考化合物的不同加合离子的保留时间和质荷比进行化合物识别。设定的参数为：筛选校正峰的最低样本检出率为80％，筛选校正峰的保留时间段设定为0-1,1-2,2-5,5-8,8-10,10-11,11-13,13-16分钟，每个时间段对应的最低峰强度为2+e07,5+e06,8+e06,2+e06,9+e05,5+e07,7+e07,9+e07，筛选校正峰的化合物识别保留时间窗口为0.08分钟，筛选校正峰的化合物识别质量精度偏差为5ppm，化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
[0123]
1)保留时间漂移校正的具体操作：(a)筛选校正峰：按筛选校正峰的最低样本检出率(80％)，在各保留时间段(设定为0-1,1-2,2-5,5-8,8-10,10-11,11-13,13-16分钟)中用每个时间段对应的最低峰强度(分别为2+e07,5+e06,8+e06,2+e06,9+e05,5+e07,7+e07,9+e07)从参考化合物的色谱峰特征(特征峰)中筛选出备选校正峰。备选的校正峰输入保留时间校正及化合物识别模块，该模块中加载有参考数据库，参考化合物的信息(如表3所示)已存储于参考数据库中，依据各参考化合物的不同加合离子的保留时间和质荷比进行化合物识别，筛选校正峰的化合物识别保留时间窗口为0.08分钟，筛选校正峰的化合物识别质量精度偏差为5ppm。通过筛选的备选校正峰即为校正峰，匹配上的参考化合物的预期保留时间即为校正峰的预期保留时间(图6中虚线表示)。对各校正峰保留时间的缺失值进行填充，填充值为该校正峰在所有检出样本中保留时间的中位值。(b)经过缺失值填充，每一个校正峰在所有样本中都有保留时间的值。将所有样本中的各校正峰一一对应后，将各校正峰的保留时间调整为预期保留时间，其余峰随之进行线性缩放，即达到了保留时间校正的目的。
[0124]
本方式筛选出的参考化合物的预期保留时间与实测保留时间见图6(以其中一个参考化合物anc4688为例)，图中横坐标表示上机检测顺序，纵坐标表示保留时间，参考化合物的实测保留时间以点示意，预期保留时间以虚线表示。该图显示了该参考化合物的保留时间漂移情况。
[0125]
2)对采用内部校准的校正方式进行保留时间校正后的a1和b1数据进行化合物注释，即与数据库内的参考化合物特征进行匹配，产生正模式靶向数据集c2，将剩余未能匹配成功的色谱峰进行样本间的峰匹配，归为正模式非靶向数据集d2(数据集形式同前)。按照同样的方式，由前述数据集a2和b2，产生负模式靶向数据集c5和负模式非靶向数据集d5。
[0126]
样本间峰匹配的方法为，所有色谱峰在样本间进行两两对比，质量精度偏差5ppm以内和保留时间偏差为0.03分钟以内的峰，可以认为是同一个物质在不同样本间检测出的峰，因此归为同一个非靶向色谱峰。
[0127]
方式三：不进行保留时间校正的化合物识别
[0128]
作为步骤s2的又一种方式，不对保留时间漂移进行校正，直接将前述数据集a1和b1输入化合物识别模块，该模块中加载有参考数据库，参考化合物的信息(如表3所示)已存储于参考数据库中，按设定的参数依据各参考化合物的不同加合离子的保留时间和质荷比进行化合物注释，产生正模式靶向数据集c3，将剩余未能匹配成功的色谱峰进行样本间的峰匹配，归为正模式非靶向数据集d3(数据集形式同前)。按照同样的方式，由前述数据集a2和b2，产生负模式靶向数据集c6和负模式非靶向数据集d6。设定的参数为：化合物识别的保留时间窗口为正负0.4分钟，化合物识别的质量精度偏差为10ppm。
[0129]
样本间峰匹配的方法为，所有色谱峰在样本间进行两两对比，质量精度偏差5ppm以内和保留时间偏差为0.03分钟以内的峰，可以认为是同一个物质在不同样本间检测出的峰，因此归为同一个非靶向色谱峰。
[0130]
化合物识别使用的参考化合物的信息示例(以一个校准化合物rfc0001为例)见表3。
[0131]
表3：参考数据库(以d-(+)-catechin(儿茶素)示例)
[0132]
modecompound_idcompound_nameadductmzrtintratioposrfc0001d-(+)-catechin[m+2h-hcooh]2+123.04474.9168.00+e081.0posrfc0001d-(+)-catechin[m+h-h2o]+273.07574.9161.92+e070.24posrfc0001d-(+)-catechin[m+h]+291.08654.9167.2+e060.09negrfc0001d-(+)-catechin[m-h]-289.07174.9181.00+e081.0negrfc0001d-(+)-catechin[2m-h]-579.15034.91800
[0133]
【步骤s3】对步骤s2不同方式下获得的正模式靶向数据集和负模式靶向数据集(c1和c4，c2和c5，c3和c6)的化合物分别进行匹配审查，按照只保留主要加合离子峰(即各加合离子峰中相对强度最高的峰)的原则扣除次要加合离子峰，并进行冲突化合物审查，对冲突化合物合并或过滤，针对不同方式分别建立靶向合并数据集(e1、e2和e3)。
[0134]
将上述步骤s2方式一中得到的正模式数据集c1和负模式数据集c4输入化合物匹配及冲突化合物审查模块，按参考化合物主要加合离子和各次要加合离子之间的关系，与样本中与之匹配的色谱峰组合间的关系进行比较。如果关系相似，认为是正确的化合物匹配结果并给匹配度评级，按照只保留主要加合离子峰的原则扣除次要加合离子峰；如果关
系不相似，则认为是错误的匹配结果，此色谱峰组合归入未能匹配成功的色谱峰矩阵。关系相似至少需要满足以下条件：1)匹配到了主要加合离子，该离子峰在所有加合离子峰中强度最高；2)匹配到了第二加合离子或者第三加合离子，该离子峰的保留时间与主要加合离子峰的保留时间偏差在0.01分钟以内，且该离子峰的强度与主要加合离子峰强度的比值和参考化合物该次要加合离子峰与参考化合物主要加合离子峰强度的比值近似(两个比值偏差在0.1以内)；3)主要加合离子的检出样本数要不低于次要加合模式的检出样本数。
[0135]
化合物匹配审查示例如下，参见表4。
[0136]
表4：化合物匹配审查(以儿茶素示例)
[0137]
peakmodeadductmean_rtmean_intobs_ratiocountrtintratio111pos[m+2h-hcooh]2+4.9891.00+e061.010004.9168.00+e081.0112pos[m+h-h2o]+4.9882.80+e050.285004.9161.92+e070.24113pos[m+h]+00004.9167.2+e060.09114neg[m-h]-4.9685.10+e051.05004.9161.00+e081.0115neg[2m-h]-4.9734.08+e050.82004.91600
[0138]
如此例所示，正模式和负模式数据集c1(离子峰111-113)和c4(离子峰114-115)中都匹配到了参考化合物——标号rfc0001的校准化合物儿茶素，审查认为其中正模式数据集c1的111号峰匹配满足了1)匹配到了主要加合离子[m+2h-hcooh]2+，该离子峰在所有加合离子峰(111-115号)中强度(mean_int)最高；2)匹配到了第二加合离子[m+h-h2o]+，该离子峰(112号峰)的平均保留时间(mean_rt)为4.988分钟，主要加合离子峰(111号峰)的保留时间(mean_rt)为4.989分钟，二者的差值为0.001分钟，在0.01分钟以内；且该112号离子峰的强度与111号主要加合离子峰强度的比值(obs_ratio)为0.28，参考化合物该次要加合离子峰与参考化合物主要加合离子峰强度的比值(ratio)为0.24的差值，二者差值为0.04，小于0.1；3)对应主要加合离子峰的检出样本数(count)为1000，高于对应第二加合离子峰的检出样本数(count)500。因此审查可以确认c1中的色谱峰111匹配的化合物是儿茶素。
[0139]
步骤s3中冲突化合物审查时，在匹配审查确认化合物匹配结果正确的色谱峰中，比较具有相同保留时间和质荷比的色谱峰并进行筛选，保留的冲突化合物至少需要满足以下特征：1)质量精度偏差小于5ppm，2)冲突化合物的检出样本数不低于此冲突组中最高检出样本数的90％。不满足该标准的色谱峰被过滤。保留的化合物以化合物名称、保留时间偏差和质量精度偏差以“|”符号隔开合并表示。
[0140]
化合物冲突过滤示例如下，参见表5。
[0141]
表5：化合物冲突审查(以儿茶素和4',5,5'-三羟基-3-甲氧基-2'-甲基-2-联苯羧酸示例)
[0142][0143]
如此例所示，表5中的色谱峰111对应化合物儿茶素和色谱峰222(也如表4示例方
式从c1数据集中确认得到)对应的化合物4',5,5'-trihydroxy-3-methoxy-2prime-methyl-2-biphenylcarboxylic acid(4',5,5'-三羟基-3-甲氧基-2'-甲基-2-联苯羧酸)有相同的保留时间(mean_rt)和质量精度(mean_mz)，且与各自参考保留时间偏差(rt_diff)都小于0.01分钟，质量精度偏差(mz_difff_ppm)都小于5ppm，检出样本数量(count)相等，都是[m+2h-hcooh]2+加合离子且都是主要加合离子。因此表5中的色谱峰111和色谱峰222对应的为冲突化合物，合并为同一化合物，化合物名称(compound)为“儿茶素|4',5,5'-三羟基-3-甲氧基-2'-甲基-2-联苯羧酸”，保留时间偏差(rt_diff)为
“‑
0.073|0.071”分钟，质量精度偏差(mz_difff_ppm)为
“‑
3.1|0.1”ppm。将该化合物注释信息归入靶向合并数据集e1中。
[0144]
按照同样的方式，由数据集c2和c5获得靶向合并数据集e2，由数据集c3和c6获得靶向合并数据集e3。
[0145]
【步骤s4】当进行了多种方式的保留时间校正时，将步骤s3获得的靶向合并数据集(e1、e2和e3)，对不同保留时间校正方式下同一色谱峰在不同样本内的强度进行合并，保留最高强度，建立平台靶向数据集(f)。
[0146]
由于步骤s2中采用了多种方式的保留时间校正，将步骤s3中得到的三个靶向合并数据集e1，e2和e3输入保留时间校正审查模块，在该模块评估不同保留时间校正方式下每个化合物注释的表现，每个注释取在不同保留时间校正方式下在每个样本(smp0001
……
)中的最高强度。最终获得平台靶向数据集f，数据集f形式为二维矩阵，“行”为每一色谱峰，“列”为该色谱峰注释上的化合物信息和该峰在每个样本中的强度，表6以色谱峰111为例显示其在数据集f中的形式。
[0147]
步骤s4示例见图7。此示例是某化合物在各检出样本中的保留时间(上图)和强度(下图)，按检测训练从左到右排列。参照样本的保留时间校正方法用交叉
“×”
表示，内部校准的保留时间校正方法用实心三角
“▲”
表示，不校正保留时间的方法用空心圆圈
“○”
表示。重叠的图形表示在不同的校正方法下找到的相同的色谱峰。此示例可以认为不校正保留时间的方式具有最高的样本检出率(出现该化合物的样本数量最多)，因此，在平台靶向数据集f中保留该化合物在不校正保留时间的方式下找到的化合物信息和样本强度。
[0148]
【步骤s5】不管执行了几种保留时间校正方法，都只使用其中一种校正方式产生的非靶向数据集和靶向数据集进行合并，来计算两两色谱峰之间的相关性并过滤冗余非靶向色谱峰。
[0149]
例如，将步骤s2和s3中获得的使用相同保留时间校正方式下产生的非靶向数据集(d1和d4)和靶向合并数据集(e1)合并，计算两两色谱峰之间的相关性，保留低于相关性阈值的非靶向色谱峰，与步骤s4中得到的平台靶向数据集(f)合并，获得去冗余数据集(g)。
[0150]
将步骤s2中采用内部校准的校正方式校正得到的非靶向数据集d1和d4，步骤s3中得到的靶向合并数据集e1一起输入冗余峰过滤模块，在该模块中对两两色谱峰的峰强度进行相关性比对，设定相关性系数阈值为0.9，保留相关性系数小于0.9的非靶向色谱峰，并与步骤s4中得到的平台靶向数据集f合并，获得去冗余数据集g，数据集g的形式与数据集f相同。
[0151]
【步骤s6】对步骤s5获得的去冗余数据集(g)进行缺失值填充处理后进行批次效应校正，然后进行标准化转换，得到鉴定出化合物的、去除了批次效应的代谢组数据集(h)。
[0152]
具体操作为：将去冗余数据集g输入批次效应校正模块，首先进行缺失值填充处理，然后进行批次效应校正，最后进行标准化转换，得到最终代谢组数据集h。其中：
[0153]
缺失值填充处理的方法包括采用python语言的pandas包进行，针对样本检出率低于50％的色谱峰，采用仪器检出下限的定值插补的方式进行处理，对样本检出率高于50％的色谱峰，采用插补该色谱峰中位值的方式进行处理，将色谱峰的缺失值填充完整。
[0154]
批次效应校正的方法采用python语言的normae
[0155]
(https://github.com/luyiyun/normae)脚本进行，输入如表1和表2中所示的测量批次和测量顺序，使用默认参数，对色谱峰在不同批次的样本之间的强度进行校正，去除仪器响应波动的影响。该步骤同时对参照样本和生物样本进行批次效应校正，以便对批次效应的校正效果进行评估。
[0156]
批次效应的校正效果评估也称批次效应校正效果分析。首先，对于去冗余数据集g进行缺失值填充处理后的数据，做pca分析，作为批次效应校正前的数据。具体而言，使用r语言的prcomp()函数，观察参照样本与待测生物样本代谢组在pc1和pc2上的空间位置。然后，对于进行了批次效应校正后的数据，再做pca分析，作为批次效应校正后的数据。同样，使用r语言的prcomp()函数，观察参照样本与待测生物样本代谢组在pc1和pc2上的空间位置。批次效应校正效果分析的示例见图8。图中灰色三角为待测样本，黑色方块为参照样本。可见批次效应校正前，参照样本和待测样本空间位置重叠，而批次效应校正后，两组的空间位置不重叠且聚集在各自的中心，认为批次效应校正的效果好。其余情况可认为效果不佳。
[0157]
数据归一化的方法选择均值归一化，即将每个样本中的色谱峰强度减去样本色谱峰强度的均值。
[0158]
标准化转换的方法选自log2标准化方法，即将样本中每个色谱峰的原始强度转化为以2为底的对数。将转化后的数据归入最终代谢数据集中。最终代谢组数据集h中某一化合物的示例见表7(以d-(+)-catechin儿茶素示例)。
[0159]
以上详细介绍了本发明的代谢组lc-ms下机数据的处理方法，对应的，本发明进一步提供一种生物代谢组的液相色谱质谱联用仪数据的处理装置，能够执行前述对lc-ms下机数据处理方法的操作，所述装置包括：
[0160]
获取模块，用于获取下机数据、上机序列集及参数文件，并对下机数据的一级全扫描图谱信息进行转换得到数据集；
[0161]
保留时间校正及化合物识别模块，设定至少一种方式，用于校正样本间的保留时间漂移后，再将每个样本内的色谱峰与参考数据库内的参考化合物特征进行匹配，将匹配上的化合物的色谱峰归为靶向数据集；并将剩余未能匹配成功的色谱峰进行样本间的峰匹配，获得非靶向数据集；
[0162]
化合物匹配及冲突化合物审查模块，用于审查每个样本内的色谱峰与数据库内的参考化合物特征进行匹配的关系是否正确，过滤或合并具有相同保留时间和质荷比的冲突化合物的峰；
[0163]
保留时间校正审查模块，用于确定不同保留时间校正及化合物识别方式下匹配上相同化合物的色谱峰；
[0164]
冗余峰过滤模块，用于过滤相关性高的色谱峰；
[0165]
批次效应校正模块，用于对所获得的去冗余数据进行缺失值处理，然后根据参数
文件中的批次信息对缺失值处理后的数据进行批次效应校正，然后进行数据归一化，最后进行标准化转换，得到代谢组数据集。
[0166]
本发明还依据前述lc-ms下机数据处理方法，提供一种计算机可读存储介质，其上存储有计算机程序，该程序被计算机的处理器执行时实现包括步骤s1至步骤s6的处理方法。
[0167]