CN111272925A - 一种基于特征值确认的全二维气相色谱调制峰多层次限制归类方法 - Google Patents

Abstract

本发明涉及一种全二维气相色谱调制峰智能归类的方法，所述第二维(²D)调制峰来源于全二维气相色谱(GC×GC)‑质谱(MS)联用仪器，其中包括三个层次的限制。可用于：复杂体系的GC×GC‑MS三维(3D)数阵中共流区域的定位，便于数学分离。所述方案通过：1)将原始数据转换为二维(2D)切片的平铺矩阵；2)背景扣除和峰顶(宽)计算；3)以²D峰顶点、峰起点或峰终点时间差为第一个限制；4)构筑²D峰顶子矩阵，以第一维(¹D)剖面极值数为第二个限制；5)基于移动窗口和²D调制峰质谱信息的特征值计算为第三次限制；6)秩图计算与共流峰解析。所述方案在合适参数设置条件下，运行快速，结果可靠。

Description

一种基于特征值确认的全二维气相色谱调制峰多层次限制归 类方法

技术领域

本发明属于全二维气相色谱分析技术领域，具体地，涉及一种自动归类第二维(²D)调制峰的方法。

背景技术

全二维气相色谱(GC×GC)是一种新型分离分析技术，它将第一维(¹D)色谱柱的流出物质重进样到第二维(²D)色谱柱中进行再次分离和检测，解决了传统一维气相色谱在分离复杂样品时峰容量不足的问题。在这个分离系统中，²D调制峰矩阵在¹D方向上按一定调制周期叠加形成GC×GC三维(3D)数阵，研究者能从独立的3D峰中获取准确的定性和定量信息。然而，极端复杂样品的GC×GC数阵通常包含几百种甚至更多化合物的信息，依然存在少部分共流3D峰。针对这个问题，化学计量学家们通过目标子数阵截取和一系列的二阶校正方法从黑色体系中获取准确的化学信息。Tauler团队对不同样品GC×GC-TOFMS数据采用小波压缩和结构分析后，运用PARAFAC和MCR-ALS工具解析共流峰。He等提出2种数学分离方案用于GC×GC共流峰解析：1)，将3D子数阵转换为平铺的二维(2D)子矩阵，再使用矩阵分辨算法解析，如直接推导式演进特征投影法(HELP),选择性离子分析法(SIA)等；2)，用交替三线性分解法(ATLD)、平行因子分析法(PARAFAC2)、多元曲线分辨-交替最小二乘(MCR-ALS)等解析多个2D子矩阵叠加的目标3D子数阵。

然而，目标3D子数阵截取需要专业知识储备和人工判别，很难普及和自动化完成。更重要的是，相邻调制周期间的同一或多个组分的2D调制峰(峰簇)相互搀杂，使得目标子数阵截取区域不清晰。因此在进行二阶校正之前，一种子矩阵调制峰归类算法需要被应用在复杂体系的GC×GC数阵。这种预处理方案应能智能地将相邻2D切片的调制峰归类成不同共性的数据集，再进一步对有共流现象的数据组进行多元校正与多元分辨。在这个方案中，峰检测对于2D切片调制峰归类至关重要。目前，基于连续小波的算法被广泛应用于峰检测，如Haar小波。峰检测之后，更需要根据化学信息判断多个峰是否为(包含)相同组分，最终实现调制峰归类。van Stee等使用基本的色谱定律来计算2D峰的理论形状(来自同一分析物的一组峰)，以确定每个化合物的峰预计出现的区域。Peters等利用两条规律实现GC×GC–FID峰的聚类：²D保留时间的差异、¹D剖面的极值个数。但是，²D保留时间差异的阈值设置以及化合物的相似保留行为，使这种方法在高度重叠的GC×GC数阵中不能被较好地应用。因此，质谱信息成为GC×GC 3D数阵中²D调制峰归类的重要参数。Zheng等通过待校峰顶与参考峰顶的质谱相关系数来确定其对应的漂移候选点数为最佳的漂移点，从而实现各谱图的峰校正。然而，不同来源样品间的重叠峰或低信噪比峰的质谱相似度差异较大，无法达到研究者的预期目标。Zeng等提出了交替移动窗口因子分析法(AMWFA)，通过交替搜索发现和有效利用隐含在两组数据中的信息。Yang等提出基于质谱信息的子窗口因子分析法(SFA-MS)，它通过计算两组峰之间是否含有共同组分来最终实现峰对准。在本发明中，质谱信息和子窗口因子分析法(SFA)也可以被利用进行不同子矩阵的共有组分识别，最终实现GC×GC数阵的²D调制峰归类。

发明内容

本发明的目的是突破GC×GC仪器分析极端复杂样品依然存在的瓶颈，提供一种自动归类²D调制峰的方法。通过本发明，能够实现三维共流峰区域的快速定位(或二维共流峰的准确识别)，便于进行下一步的多元校正与多元分辨。

本发明是²D调制峰自动归类方法，适合于GC×GC-MS数据，也适用于LC×LC-MS数据。如图1所示，主要涉及三层次：²D保留时间限制、¹D剖面峰形限制、基于移动窗口和“2D切片”调制峰质谱信息的特征值计算。

本发明的²D调制峰自动归类方法如下：

(1)数据转换：将Shimadzu GC×GC-qMS原始qgd格式转换为mat格式，转换的新格式在Matlab软件中运行为2D切片的平铺矩阵。

在3D数阵中，²D维度被看做I轴，其中I(1，2，…i)是每个调制周期中的质谱扫描次数；J(i，2，…，j)是为每次扫描记录的预定义m/z范围；每个调制产生一个数据矩阵X(I，J)，X则为调制周期内的信号矩阵。假如有k个调制周期，则有X₁(I，J)，X₂(I，J)，...，X_k(I，J)，由于调制周期固定且扫描预设置m/z范围相同，因而每个X(i，j)在I轴上和J轴上等距，将其叠加起来能形成GC×GC-qMS的3D数组X(I，J，K)。在这个3D数阵中，I维(1,2，…，i)代表每个调制中的质谱扫描次数，J维(1，2，…，j)代表每次扫描的预定义m/z范围，K维(1,2，...，k)代表调制周期的数量。如图2所示，GC×GC-MS 3D数阵结构上又可以看作为一个三线性成分模型，X_aug(I，J，K)按K轴可切割成k个2D矩阵X(I，J)；再将这些X₁(I，J)，X₂(I，J)，...，X_k(I，J)按行排列，能形成一个大的2D矩阵X_aug(I×K，J)，我们称之为“2D切片平铺集”。

(2)“2D切片平铺集”的背景扣除与峰检测：在GC×GC-MS原始数据中，背景噪声通常会干扰色谱峰宽和峰高的计算，从而影响调制峰的识别。本发明采用了自适应迭代重加权惩罚最小二乘算法(airPLS)进行了“2D切片平铺集”的背景扣除，也可采用其它方案排除基线漂移。

背景扣除之后，以Haar函数为小波母函数的连续小波变换(Haar CWT)可增强求导过程中的信噪比，能很好地估算出峰的位置(宽度)。

连续小波变换可以看作对信号和缩放及平移的小波母函数乘积进行积分，设小波母函数ψ(t)伸缩因子为a，平移因子为τ，并记平移伸缩后的函数为ψa，r(t)，则：

Haar函数是最简单的小波母函数，它可以用下面公式描述：

Haar CWT可以用于计算色谱信号近似的n阶导，通过求导的方式可以得到每个色谱峰的起始点与结束点，即为峰宽；通过以Haar函数为母函数的小波变换检测出局部最大值，即为峰顶点；最后用信噪比剔除一些信号不够强的小峰。在本发明中，峰簇在峰检测过程中不好确认起始和结束点，被视为一个峰。

本发明采用了Haar CWT方案，不限制采用其它方案进行峰检测。

(3)²D调制峰保留时间限制：将3D数阵转换为“2D切片平铺集”后，同类/不同类调制峰在“2D切片平铺集”中相互交错的现象会使人工识别耗时耗力。因此在二阶校正前，应该使用一种自动识别分类方法，并确认3D目标子数阵重新构筑的范围。在本发明中，约束调制峰归类的第一个限制为²D峰顶点保留时间差。理论上，各“2D切片”中的调制峰顶点如相同(图3a)，则归为相同峰组；若不同则归为下一个峰组。然而在实际复杂样品测试中，重叠峰簇存在多个峰顶点，且不同¹D流出点(对应各“2D切片”)的成分浓度呈现规律性变化，这就会使不同切片中²D峰顶点检测结果易产生错误(图3b-e)。因此，我们除计算峰顶点的保留时间差之外，还增加了峰起点或峰终点的时间差计算。图3b、3c是峰终点时间差相同的顶点情况；图3d、3e是峰起点时间差相同的顶点情况。通过顶点、起点、终点的三种时间差限制，将²D轴保留时间差相同的调制峰归类到一个峰组。在实际的GC×GC测试中，仪器自身或者其它不理想操作会使相邻调制峰在²D轴上产生保留时间漂移现象。因此，本专利将这些时间差设置了一个阈值，小于这个阈值则归类到一个峰组；超过这个阈值，则不归类，如图3f，3g所示。

(4)¹D剖面色谱峰形限制：在(3)条件限制后，我们把相邻2D切片内同一²D保留时间的信号峰归为一类。但是，¹D色谱柱在相邻时间内洗脱出很多理化性质相似的组分，而某些组分又无法在²D色谱柱内展现出不同的保留行为。因此，经过(3)条件限制后归类的调制峰并不一定全属于同一组分，仍需要进行二次限制归类。在本发明中，我们通过单模态准则探究了二次分类法：在¹D方向，如只有一个化合物流出，理论色谱峰剖面应仅显示一个最大峰值，信号强度呈现出由低到高再到低的峰状；同理，一个化合物经过调制周期后相邻²D信号最大响应值，也应该在¹D方向上呈现由低到高再到低的峰状。

本发明计算一次归类好的调制峰子矩阵的峰顶点，形成峰顶子矩阵，然后通过求导的方式来获取峰顶点子矩阵的极值。除了起点，每一个极小值点也可为一个峰/峰的分界点，即得到不同的感兴趣组分的归类。在每一个峰组中，只能允许存在一个极大值峰(峰顶点)。如图4所示，a-i为一次限制分类后得到的一个峰组，内含九个峰簇；在二次限制分类中，将其峰顶点收集到一个矩阵内，求导，极小值(峰/峰交界)为d，g，则该峰组分为a-c、d-f、g-i三个峰组。

(5)基于移动窗口和调制峰质谱信息的特征值计算：经过两次限制后归类的一类色谱峰(峰簇)，可以通过各“2D切片”调制峰质谱信息相关性再次确认各²D调制峰间是否为共有组分。尽管GC×GC系统有很高的峰容量，但各“2D切片”的²D调制峰依然出现共流现象。在这些峰簇中，重叠或包埋现象会使各采样点的质谱信息存在差异，因此在进行各调制峰质谱信息相关性研究时，不能采用质谱相似度作为评价方法。如图5所示，SFA中的特征值在本发明中被用来判断两矩阵的质谱拟合程度，在有干扰峰的影响下也能很好地识别出是否含有共有组分。

在这个方案中，两个子矩阵窗口被选择用于一次的相关性比较，Y是基矩阵，用来确定感兴趣组分的化学信息；X是目标矩阵，用于搜索感兴趣组分。因此，选择一个具有代表性化学信息的基矩阵对峰组内各调制峰的划分具有重要的意义。如果一个峰组内基矩阵Y的调制峰是其余组中的干扰组分，就易引起错误的计算结果。在经过(4)中根据¹D剖面峰形限制之后，响应最大的调制点所容纳的调制峰(组分信息)最具有代表性意义。在这个剖面，¹D理论峰起点和¹D理论峰终点之间容纳的调制矩阵被视为进行两两相关性比较的数据范围。在响应最大的调制矩阵中，Haar小波被使用来计算最高调制峰(峰簇)起点和终点，该段数据被用作包含感兴趣组分信息的基矩阵Y。¹D组内其余n个调制点容纳的矩阵也通过Haar小波确认其窗口大小后作为目标矩阵Xi(i＝1,2,3…n)，用来搜索是否含有基矩阵Y中所包含的化学信息。紧接着，基矩阵Y和目标矩阵Xi进行两两比较，获取特征值(阈值0.85，用户也可以自行输入)，若特征值接近1，则说明该目标矩阵Xi包含基矩阵Y中的化学信息，会将这些调制峰(峰簇)归类为一个峰组；若没有接近阈值的特征值，则该目标矩阵Xi被视为不含有Y中的化学信息，会被移除基矩阵Y所在的峰组；被移除的矩阵重新组成一个峰组，并与其余特征值不含有接近阈值的Xi矩阵进行特征值(阈值0.85)计算，依次循环，直到该Y值所在的峰组内的除Y外的n个峰确认完为止，则可以获得含有同一感兴趣组分的新的一个或多个峰组。

首先通过对目标矩阵X与基矩阵Y进行SVD分解，得到相应的两组正交基矩阵与奇异值矩阵：

其中列正交阵U和行正交阵V分别表示得分矩阵(Score matrix)和荷载矩阵(Loading martix)，S为奇异值矩阵，其对角线表示奇异值，最大奇异值所对应的荷载矩阵的列向量表示其主要化学信息，次大奇异值对应的向量则表示次要化学信息，若有明显小于前几个的奇异值，则有理由相信其是噪声所引起的，E为误差矩阵。

为了确定基矩阵中的感兴趣成分，我们假设X,Y中组分数分别为m，n，共有组分数为c，设代表其组分的荷载矩阵中的向量集分别为E＝{e₁，e₂，…，e_n}，F＝{f₁，f₂，…，f_m}，X,Y中的共有组分S则可以由E和F的线性组合表示，如方程(6)所示：

其中，符号a_k与b_k分别为基矩阵E和F的线性组合系数，由于噪声干扰等误差，Ea_k和Fb_k不会完全相等。我们通过一个拟合的目标函数评价Ea_k和Fb_k的拟合程度：

我们设a_k ^TE^TFb_k＝d_k为特征值，通过应用E^TE＝I,F^TF＝I和s_k＝Ea_k＝Fb_k,求解a_k ^TE^TFb_k＝d_k的特征值，函数如等式(7)所示。

由方程(8)可知，如果向量S_k确实是矩阵X和Y中某一共有组分的波谱，则系数向量a_k与b_k就分别对应乘积矩阵E^TFF^TE和F^TEE^TF中特征值为1的特征向量。

这样，若特征值接近1，则函数f(a_k，b_k)越小，也就是说，Ea_k与Fb_k的拟合度越高，其对应的S_k就越可能是共有组分；相反，特征值越接近0，则函数f(a_k，b_k)越接近2，代表Ea_k与Fb_k的拟合度越低。其对应的S_k就越不可能是共有组分。通过特征值d_k可以用来判断X和Y子矩阵之间质谱信息相关性，并得到相关目标组分的共有组分秩图。通过共有组分秩图，可以清晰地知道目标矩阵中是否含有基矩阵中包含的化学信息。在本文中矩阵的每行代表一个光谱数据，矩阵的列数代表色谱峰的点数。因为S_k为X和Y的共有组分，由于基矩阵中只需要一个感兴趣组分，我们选择基矩阵的荷载矩阵内也只需要选取最大奇异值对应的向量，然后通过目标矩阵中的子窗口移动搜索是否含有基矩阵中的感兴趣组分的共有组分。

(6)如图6所示，²D调制峰自动归类后，通过秩图计算确认是否为共流峰。如为共流峰，可选择矩阵解析，如HELP、SIA法等；或重构三线性成分模型，如ATLD、PARAFAC2等。

从而，我们能快速且准确地锁定3D数阵中的共流峰区域，进一步开展多元校正与多元分辨。

相较于现有技术，本发明提供的GC×GC数阵中²D调制峰智能归类方法(ICMP)具有以下有益效果：

1.以模拟数据为例，三种限制条件结合的ICMP能较好地将各2D切片的²D调制峰自动归类，减少仅用前两种条件所带来的假阳性和假阴性风险。

2.基于移动窗口和调制峰质谱信息的特征值计算，适用于复杂体系中的调制峰(簇)(往往具有共流现象)。

3.以真实复杂的GC×GC-MS数阵为例，ICMP能在合适的参数设置下，自动、快速、准确地归类各²D调制峰，并识别3D共流峰区域，便于下一步的多元分辩与多元解析。

附图说明

图1 ICMP流程

图2全二维气相色谱-质谱联用的数据结构及其转化

图3“2D切片”中²D调制峰保留时间限制

图4“2D切片”的¹D剖面峰形限制

图5基于移动窗口和“2D切片”调制峰质谱信息的特征值计算

图6秩图计算与共流峰解析

图7模拟数据的ICMP运行结果

图8 ICMP操作界面

图9植物提取物GC×GC-qMS数阵的ICMP处理结果

图10SFA验证结果(a为2D切片平铺集色谱图和基矩阵Y质谱图，b为目标矩阵X1-Y特征值图、共有组分秩图和X1质谱图；c为目标矩阵X2-Y特征值图、共有组分秩图和X2质谱图；d为目标矩阵X3-Y特征值图、共有组分秩图和X3质谱图)

图11不同参数h设置下一段真实数据的归类结果(a、b和c分别为h＝5、h＝3、h＝0)

图12 ICMP-ATLD处理结果

图13 ICMP-SIA处理结果

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的实验方法，通常按照常规条件或按照制造厂商所建议的条件。除非另外说明，否则所有的百分数、比率、比例或份数按重量计。

除非另行定义，文中所使用的所有专业与科学用语与本领域熟练人员所熟悉的意义相同。此外，任何与所记载内容相似或均等的方法及材料皆可应用于本发明方法中。文中所述的较佳实施方法与材料仅作示范之用。

实施例1模拟3D数阵中各²D调制峰的智能归类

模拟3D数阵由真实的GC×GC-qMS数据通过剪切和拼接而成，用于解释ICMP算法的工作原理。首先，将模拟3D数阵转换为“2D切片平铺集”，再使用小波法检测峰(峰簇)信息。首先，在图7(a)曲线部分，小波法能在不同的信噪比情况下很好地估算出峰的顶点和峰的宽度(峰起点至峰终点)。为区分不同的峰(峰簇)，ICMP根据²D色谱保留时间约束进行峰(峰簇)分类。如图7(b)所示，²D坐标在138dot-166dot之间的峰组和²D坐标在173dot-190dot之间的峰组被划分。然而，相似结构的化合物往往在²D色谱柱上也具有相似的保留行为，不利于归类。因此，根据¹D色谱行为规律再次限制，一个峰组内的极值只能允许存在一个最大值；相反，ICMP就将具有相同²D保留时间的峰(峰簇)分割成不同的峰组。如图7(c)，根据¹D剖面峰形限制，²D坐标在173dot-190dot之间的峰组被分割成了两个峰组。通过两次限制后，仍然需要对比质谱信息来确认峰组内的峰是否含有同一感兴趣组分。在本发明中，研究者通过SFA确认峰组内各²D调制峰是否包含一个共有成分。若特征值大于阈值，则含有共有感兴趣组分，确认该峰组归类正确；若特征值小于该阈值，则该峰不含感兴趣组分，被排除出该峰组。如图7(c)中峰组蕴涵质谱信息的目标矩阵X与基矩阵Y进行了比较，导出特征值图均远小于阈值，故而将其再次分割成两个峰组(图7d)。

实施例2植物提取物GC×GC-qMS数阵的各²D调制峰归类结果

崖柏精油GC×GC-qMS分析：

GC条件1：Rtx-5Sil MS(30m×0.25mm×0.25μm)+BPX-50(2.5m×0.1mm×0.1μm)；柱温温度：100℃，进样口温度：280℃；分流比20:1；进样量：1.0μL；升温程序：100℃(保持2min)，以2℃/min升温至280℃，并保持5min。

GC条件2：Rtx-wax(30m×0.25mm×0.25μm)+BPX-5(2.5m×0.1mm×0.1μm)；柱温温度：50℃，进样口温度：250℃；分流比20:1；进样量：1.0μL；升温程序：50℃(保持2min)，以2℃/min升温至250℃，并保持40min。

MS条件1：载气：氦气；载气控制方式：恒线速度；离子化方式：EI；离子源温度：200℃；色谱接口温度：300℃；采集方式：全扫描模式(SCAN)；扫描速度：20000amu/sec；质量数范围：m/z 45～339amu；调制周期：4sec；热喷温度：300℃；热喷持续时间：350msec；冷喷流量：6L/min；采样频率：50Hz。

MS条件2：载气：氦气；载气控制方式：恒线速度；离子化方式：EI；离子源温度：200℃；色谱接口温度：250℃；采集方式：全扫描模式(SCAN)；扫描速度：20000amu/sec；质量数范围：m/z 45～339amu；调制周期：4sec；热喷温度：300℃；热喷持续时间：350msec；冷喷流量：6L/min；采样频率：50Hz。

ICMP操作步骤和参数：

1.将原始数据转换为mat格式，切割感兴趣的子矩阵，并命名为A1.mat，本发明给出两个截取的数据段‘53501-85300dot’，‘69001-75300dot’用来测试ICMP；

2.如图8，打开’MAIN’主窗口，并通过点击上方菜单栏选项进行操作；

3.使用’FILE’的’load_file’导入数据；

4.使用’FILE’中的’TIC’得到TIC图；

5.通过’BSB’进行背景扣除，以下为对应函数：

[xbc,xb]＝airPLS(TIC,lambda,2,0.05)；

[xs,xln]＝wavelet_denoising(xbc,l,6,'soft')；

6.在’PEAK detection’输入参数，点击’plot’按钮即可自动执行峰检测操作；

peaks0＝peak_detection(xs,2,xln,phi)；

peaks＝peak_clustering(xs,peaks0,1)；

[PEAKS]＝peak_grouping(peaks,mt)

％mt为调制周期，单位为秒

％phi为信噪比阈值

％pesks为检测到的峰值

7.通过ICMP进行自动峰聚类算法

1)通过’TDP’按钮将2D切片(按调制周期划分)的²D峰按²D保留时间差归类成TDP峰组；

[PEAKS1,TDP]＝peak_TDP(PEAKS,peaks,k,h1,h2,xs)；

％k为每个归类组跨越的最大调制周期数，一般为10个周期以上

％h1为两峰顶点在²D轴上允许的保留时间差，一般在5个扫描点以内

％h2为两峰起点在²D轴上允许的保留时间差，一般在10个扫描点以内

2)通过’CP’按钮将TDP峰组按¹D剖面色谱峰形限制再归类成CP峰组；

[CP]＝peak_CP(TDP,peaks,xs)；

％TDP为按²D保留时间差归类的峰组

％peaks为检测到的峰值

％xs为总离子色谱图数据

3)’SFA’按钮：将CP峰组按SFA检测是否含有相同组分，再进一步归类；

[CPSFA,a]＝peak_CPSFA(CP,peaks,peaks0,Y,xs,value)

％CP为按¹D剖面色谱峰形限制再归类的峰组

％peaks为检测到的峰值

％Y为原始数据

％xs为原始总流离子图数据

％value为SFA特征值选择阈值

4)导出SFA确认归类之后的2D图，点击对应的’plot’将不同归类峰组按不同颜色表示出来，并标以序号，在对应的坐标轴内显示。

CPplot(CP,xs)

8.通过’FSMWEFA’按钮绘制目标峰组秩图，并导入到坐标轴内，可判别是否为共流峰；

fun＝fsmwefa(Y,PG,nwin,CPSFA)

PG为选择的峰组序号

niwn为秩的个数

9.1)可直接对目标²D共流峰采用矩阵法解析；

2)对目标²D共流峰组叠加，ATLD半自动化解析重叠峰。

[AA]＝ATLD_plot(CP,xs,Y,PG,Kn)

[A,B,C,SSR,Tole,N_iter,AAA]＝CP_ATLD(CP,xs,Y,i,Kn,num,epsilon,Tol)

％PG为需要处理的峰组序号

％Kn为峰组内可能的峰数量

％num为可能的成分数

植物二级代谢物种类繁多、结构复杂，使用一维色谱难以获取植物提取物中全面的化学组分信息。为解决峰容量不足的问题，GC×GC技术被广泛用于开展植物提取物的整体化学组分研究。本发明使用崖柏提取物为测试样品，Rtx-5Sil MS(弱极性)/BPX-50(中等极性)和Rtx-wax(极性)/BPX-5(弱极性)分别设置为两次分离分析的¹D色谱柱/²D色谱柱，调谐时间均为6s。两组方法均能检测出大量萜烯类化合物，但两组数据中依然含有共流峰，无法准确定性。为继续深入挖掘，化学计量学家们可通过截取目标子数阵和多元校正多元分辨方法获取隐藏信息。然而在这个进程中，3D数阵转换的“2D切片平铺集”包含有同一化学组分的多个调制峰，且多个化学成分相互交错，这使得GC×GC的目标子数阵截取困难。在本发明中，ICMP算法通过归类“2D切片平铺集”的各²D调制峰，使每一个峰(峰簇)组都有同一感兴趣组分(群)，从而使3D子数阵截取界限清晰。

图9a)是崖柏提取物进行GC×GC分析的局部3D轮廓图，该数阵实质是²D矩阵在¹D方向上按一定调制周期叠加而成。GC×GC是一种新型分离分析技术，它将¹D色谱柱的流出物质不断地进样到²D色谱柱中进行再次分离和检测，共有700个调制周期(即¹D色谱时间范围为0-700dots)。在²D色谱再次分离过程中，流出物在6s内连续接受300次扫描，即²D色谱时间为0-300dots。为了检验ICMP对GC×GC-MS数阵各²D调谐峰的归类情况，我们截取了崖柏提取物数阵在¹D 231-251dots的子数阵进行详细的分析。如图9b)所示，该子数阵3D图能清晰地表观GC×GC-MS的“2D切片叠加集”数据结构，该段区域有多个相同的或不同的重叠峰堆叠，通过人眼难以识别其界限，给后续的“二阶校正”的人工切割目标子数阵造成了困难。为此，我们使用ICMP算法归类各²D调制峰，并识别组间界限，即确认重新构筑3D目标子数阵的范围。首先，我们将“2D切片叠加集”转换为“2D切片平铺集”(图9c)，并使用Haar小波检测22个调制峰或峰簇。为了归类这些信号，ICMP计算各调制峰(峰簇)的峰顶与峰宽，继而计算出每个调制峰(峰簇)在²D色谱保留时间的时间差。如图9d所示，通过该时间差约束峰(峰簇)分类，22个调制峰(峰簇)被归类成6个峰组。面对相似结构化合物可能具有相似²D峰位置的情况，ICMP通过¹D色谱行为规律，再次将相似²D保留时间的不同峰(峰簇)分割成不同的峰组。图9e所示，22个调制峰的6个峰组被再归类成8个峰组，分别是原第2个和第4个峰组均被检测出含有两个最大峰值，故而在最小峰值位置分离成不同的峰组。

最后，为防止²D保留时间差异及¹D剖面峰形判断错误，ICMP通过对比峰组内各调制峰(峰簇)的质谱信息来进行最后的约束。在这个过程中，SFA的特征值被用来比较基矩阵Y与目标矩阵X之间的各子窗口的质谱相关性，即判断是否含有共有感兴趣组分。首先，峰组内响应最大的调制点所容纳的调制峰被选为基矩阵，其它峰组内调制峰为目标矩阵；当两矩阵之间的特征值远小于1时，表示两峰的质谱信息不拟合，则该目标矩阵所在的峰被排除该峰组。我们选取图9e中的第3个峰组来举例说明。该段数据对应GC×GC-qMS数据240-243dots片段的4个调制峰(峰簇)，最大响应信号的调制峰所包含的数据作为基矩阵Y，与之相邻的峰簇所含数据分别作为目标矩阵X1、X2、X3，如图10a所示。我们提取基矩阵Y中的最大峰顶点的质谱信息，并将其在NIST14数据库中进行相似度匹配得到其组分为：Thujpsadiene。图10b1、10b2、10b3和10c1、10c2、10c3分别为基矩阵Y与目标矩阵X1与X2通过SFA确认之后得到的特征值图、共有组分秩图、X1与X2峰顶点的质谱信息。从图10观察到特征值dk(0.99)接近于1，其差异在仪器噪声允许范围内，表明X1包含的化学信息与Y所包含的化学信息为共有组分；而共有组分秩图表明X1矩阵的第5-9列数据与Y矩阵有共有组分；质谱信息图也验证了X1与Y具有同一组分。同理，目标矩阵X2与基矩阵Y的特征值接近于1，两者质谱信息高度相似，也具有同一感兴趣组分。而图10d1则为目标矩阵X3与基矩阵Y通过SFA确认之后得到的特征值图，特征值皆小于0.85；图10d2为共有组分秩图，表明两矩阵没有共有组分；图10d3中的X3质谱也与Y峰顶的质谱大为不同，所有结果均意味着X3与Y中没有共有组分。因此，在对8个峰组进行子窗口分析后，第3个峰组的第一个调制峰簇被确认为不含有该峰组内的感兴趣组分，故而排除该峰组，分割成两个峰组。其它峰组也都进行了ICMP特征值计算和分析，综合得到11个峰组，结果如图9f所示。在“2D切片平铺集”中，每一个2D切片的²D轴都相同(300个扫描点)，研究者能将所有2D切片再叠加在一起形成三线性成分模型(图9g)。

为了确定ICMP算法的最佳参数设置，本发明对不同参数的效果进行了对比。在前面论述中，小波算法已被证明能够检测各个信号峰，其检测峰数随信噪比(phi)的变化而变化。在本样本中，当phi设置为6时ICMP给出了合理的检测峰数。当phi越低，检测到的峰数越多；而当phi更低时，则会将噪声检测为信号峰。此外，在检测¹D色谱低强度峰中，设定正确的阈值是很有必要的。因此，针对¹D保留时间的信号强度，应为峰检测设置正确的phi值。对于²D保留时间差限制，其参数由两个阈值决定，分别是²D保留时间差，¹D保留时间长度。在这个限定中，分组效果主要由²D保留时间差所决定：²D保留时间差阈值设置越大，合并到一个峰组中的单个峰就越多，假阳性就越多；由于存在保留行为偏差，²D保留时间差阈值设置过小，就会将本该属于该峰组的峰值排除在外，造成假阴性。如图10，研究者对比了不同阈值设置下片段(231-251dots)的ICMP归类结果。在这个俯视图中，不同的散点代表不同峰组内的调制峰(峰簇)。当参数h＝6时，22个调制峰(峰簇)被归类成6个峰组(图11a)。当参数h＝3时，22个调制峰被(峰簇)归类成11个峰组(图11b)；其中，10、11峰组在¹D色谱行为约束中被归类为两个峰组，属于正确的归类；但是，1、2峰组在后面两次约束中都没有将其分开，说明应归类成一个峰组；同理，3、4峰组也应该被归类为一个峰组，所以该阈值设置错误。当参数h＝1时，由于²D保留时间轻微漂移，22个调制峰(峰簇)被归类成18个峰组(图11c)，属于同一峰组的调制峰(峰簇)被错误划分为不同峰组。综上研究，参数h在本数据中应被设置为6。因此，针对不同的数据选择不同的阈值：以归类峰准确性为目标，可以将阈值设置偏小；若数据量大，并考虑到后续还有第二次约束与特征值确认，可将阈值设置偏大。¹D色谱时间长度阈值(k)是ICMP算法允许归类的最大峰值个数，表示一个组分通过¹D色谱柱所经历的最大调制次数，一般设置为10左右，太小可能归类不完全，太大运行时间长，本发明设置为k＝14。

当感兴趣组分与干扰化合物共流时，通常要尽可能地先采用物理分离，再使用‘数学分离’来减少不必要的时间和成本消耗。针对GC×GC 3D共流峰，二阶校正是常用手段之一。在三线性分解前，重构目标3D子数阵是共流峰解析的先决条件。使用ICMP算法能有效地解决相邻调制周期间的同一或多个组分的²D色谱峰(峰簇)相互搀杂难题。从而，研究者准确地截取目标子数阵，并重构三线性成分模型。我们以三次限制归类后的第2个峰组为例，通过秩图确认该调制峰簇为重叠峰，并将这个235-239dots片段重构成一个3D子数阵。如图12所示，ATLD算法可以从这个复杂数阵中获取不同成分的定性与定量信息。其中，1：β-Himachalene，2：α-Chamigrene。对于单个²D调制峰，SIA算法能获得满意解析结果。以第238dot调制峰为例，通过秩图可判别该峰为重叠峰。基于最小熵原理可以找到两组分的选择性离子信息(图13)，再依据m/z 119和136获取两组分的纯色谱曲线，通过最小二乘法计算两组分的纯质谱。结合保留指数和NIST14数据库相似度匹配，两组分分别为1：β-Himachalene，2：α-Chamigrene。

Claims (6)

1.一种全二维气相色谱调制峰智能归类的方法，其特征在于：

所述第二维(²D)调制峰来源于全二维气相色谱(GC×GC)-质谱(MS)联用仪器，其中包括对²D保留时间限制、¹D剖面峰形限制、基于移动窗口和“二维(2D)切片”调制峰质谱信息的特征值计算等三层次限制。²D调制峰智能归类对快速且准确地锁定三维共流峰区域，便于下一步的解析具有重要的意义。

GC×GC数阵中²D调制峰智能归类方法(ICMP)步骤如下：

(1)数据转换：将GC×GC-MS原始qgd或其它格式文件转换为mat格式，转换的新格式在Matlab软件中运行为2D切片的平铺矩阵；

(2)“2D切片平铺集”的背景扣除与峰检测：采用自适应迭代重加权惩罚最小二乘算法(airPLS)或其它方法进行背景扣除，以Haar函数为母函数的小波变换或其它方法检测峰顶和峰宽；

(3)“2D切片”中²D调制峰保留时间限制：基于实际测试的复杂性，除了以²D峰顶点时间差为第一个限制，也增加了峰起点或峰终点的时间差限制；

(4)“2D切片”的¹D剖面峰形限制：通过单模态准则探究二次分类法，计算(3)归类好的调制峰子矩阵的峰顶点，形成峰顶子矩阵，然后通过求导的方式来获取峰顶子矩阵的极值，除了起点，每一个极小值点也可为一个峰/峰的分界点，同时每一个峰组中只能允许存在一个极大值(峰顶)；

(5)基于移动窗口和“2D切片”调制峰质谱信息的特征值计算：子窗口因子分析法(SFA)中特征值计算被用来判断两矩阵的质谱拟合程度，在有干扰峰的影响下也能很好地识别出是否含有共有组分，可作为智能归类的第三次限制；

(6)秩图计算：通过FSWMEFA获取秩图，直观地判断该峰组是否有共流峰；

(7)共流峰解析：使用多元校正和多元分辨法解析共流峰。

2.根据权利要求1所述的GC×GC数阵中²D调制峰智能归类的方法，其特征在于，²D保留时间限制、¹D剖面峰形限制、基于移动窗口和“2D切片”调制峰质谱信息的特征值计算等三层次限制的结合。

3.根据权利要求1所述的GC×GC数阵中²D调制峰智能归类的方法，适用于复杂体系中的调制峰(簇)(往往具有共流现象)。

4.根据权利要求1所述的三层次限制方案，可与多种背景扣除与峰检测方法组合。

5.根据权利要求1所述的GC×GC数阵中²D调制峰智能归类的方法，可智能归类GC×GC-MS仪器的²D调制峰，也适用于LC×LC-MS仪器。

6.权利要求1所述的²D调制峰智能归类方法所获结果，可用于：中药、石油和生物样品等复杂体系GC×GC-MS测试中三维共流峰区域的快速定位(或二维共流峰的准确识别)，便于多种方案的数学分离。

Images