基于GC-MS的复杂体系解析是化学计量学的重要内容。GC-MS数据解析可由如下模型刻画[1]:
其中X为气相色谱质谱联用色谱仪产生的量测矩阵, 其中每一列表示质谱中的质荷比(m/z)记录, 每一行表示色谱在不同保留时间的响应强度。组分矩阵S的每一列为化学纯组分的质谱数据, 权重矩阵C中的每一列为S中相应位置的化学纯组分在不同保留时间的色谱响应值。E为系统误差矩阵。当量测矩阵X已知, C与S均未知的情况下, 解析问题为典型的黑色体系问题, 一般仅可采用非监督算法(如基于主成分分析的方法)尝试解决[2]。
对于GC-MS数据, 充分利用质谱数据库(如NIST数据库)的信息可显著降低系统的不确定性。其方法主要分为两类:第一类, 首先进行色谱峰识别和重叠峰解析, 再利用质谱数据库进行质谱检索, 实现定性功能; 第二类, 首先根据X提供的信息, 在数据库中检索一定数量的相关纯质谱, 作为矩阵S的估计, 然后使用合适的回归算法(如非负最小二乘回归)估计矩阵C, 实现定性定量分析。第一类方法的主要挑战在于色谱重叠峰解析, 其经典方法包括多元分辨技术[3]、小波分析[4]和神经网络[5]等, 但该类方法基于色谱峰形判别, 要求色谱峰具有一定的分辨率, 解析严重重叠峰时往往失效。第二类方法的关键步骤为质谱检索和回归算法。
关于质谱检索, 可追溯至20世纪70年代提出的概率基础匹配(PBM)算法[6], 后来许多文献基于该算法进行改进, 并将其广泛应用于商业软件, 文献[7]对此类算法进行了综合比较, 并认为NIST MS Search中的算法综合性能最优。然而, 此类检索算法假设待测质谱纯度较高, 对混合质谱(色谱峰重叠时可能出现)无法保证有效检出各纯组分[8]。为在混合质谱假设下进行质谱检索, 文献[9]中提出一种参考谱加权存在指数方法, 该方法在许多后续GC-MS数据解析算法中作为检索工具[8, 10-13]。然而, 实际计算表明, 当质谱基线水平较高时, 该算法所选参考谱数目较多, 增加后续分析的不确定性。另外, 使用单一指标遍历大规模数据库, 检索的时间复杂度颇高。
关于回归算法, 许多文献[8-11]推荐非负最小二乘法解决该回归问题, 采用方差分析或主成分分析估计待测谱中所含的纯组分。但非负最小二乘法为实现数学上的最佳拟合往往将相关性较低的质谱进行强行拼凑, “过拟合”现象突出, 经常出现错误的解析结果。
鉴于上述传统算法的不足, 本文提出基于稀疏模型的GC-MS数据解析算法, 在严重重叠峰解析中取得了较好效果, 算法的主要特点如下:在质谱检索方面, 为降低检索的时间复杂度, 提出分步检索方案。首先利用质谱碎片规律, 结合索引技术进行快速粗筛; 然后, 为进一步降低所选参考谱集的规模, 提出更为精细的强峰高概率出峰准则和耐挤压性准则进行参考谱剔除。在解决非负最小二乘法的“过拟合”问题方面, 提出采用易于提取测量矩阵X“主要结构”的稀疏优化模型。
如图 1所示, 通过仪器获得GC-MS数据后, 首先通过有效的峰选择算法[12, 13]或用户交互的办法确定色谱峰, 将峰起点与终点之间的数据作为量测矩阵X, 将峰顶位置所对应的质谱作为待测混合谱实施质谱检索以获得组分矩阵S, 确定X与S后, 通过稀疏模型解算权重矩阵C, 完成解析过程。其中蓝色字体标定的部分为算法的非平凡部分, 是应当重点考察的内容。
主要质谱筛除步骤将质谱数据库中的所有质谱作为参考谱, 考察其与待测混合质谱的相关性。进行质谱分析前, 先对待测质谱与参考谱均进行规整化处理。规整化时将质谱中最大峰的强度缩放至1000, 其余各峰按比例缩放。另外, 为简化处理, 将非整数的m/z按四舍五入法则设定为整数。
分子离子峰和基峰是非常重要的质谱碎片特征, 可作为筛除标准。为简化分子离子峰的认定, 本文考虑最右端质量数, 即质谱图最右端峰簇中相对丰度最大的峰所对应的m/z, 以标准质谱数据库(如NIST谱库)为基础, 预先建立最右端质量数索引和基峰索引, 以期加快质谱检索速度。最右端质量数索引将所有参考谱按最右端质量数分类, 以最右端质量数作为类标, 每个类标存储所有对应参考谱在质谱数据库中的位置。基峰索引存储数据库中所有参考谱的基峰位置, 存储形式为键值对, 以参考谱位置为键, 以基峰位置为值。下文详述具体筛除步骤或准则。
考察待测混合质谱中的任一有效m/z, 通过查询最右端质量数索引, 可得相应m/z对应的所有参考谱。然后合并所有有效m/z的对应参考谱列表, 即得所需候选质谱集。
基于上一步所得候选质谱集, 考察其中的每个参考谱, 通过查询基峰索引获得其基峰位置(m/z)。考虑待测混合质谱中相应位置的出峰强度, 若该强度低于某阈值T(默认T=300), 则将所考虑参考谱从候选质谱集中剔除。
设定相对丰度超过一定阈值的峰为强峰, 对参考谱中的任意强峰, 若混合谱中相应位置的出峰相对丰度与该峰的比值低于阈值Q=T/1000, 则标记为异常, 若标记为异常的强峰数目超过2, 则剔除所考察的参考谱。反过来, 考察混合谱中的强峰(强度高于T), 若参考谱相应位置的出峰强度与待测谱的强度之比低于Q, 则标记为异常, 若标记为异常的峰数目超过混合谱中强峰总数的一半, 则剔除所考察的参考谱。
将参考谱与混合谱对齐后, 对参考谱各个m/z进行同比例缩放, 直至参考谱在任何m/z处的出峰都低于待测混合质谱相应位置的出峰。将参考谱的任一有效m/z在压缩后与压缩前的丰度比定义为挤压比例, 若挤压比例小于阈值Q, 则剔除所考察参考谱。
实施强峰高概率出峰准则和耐挤压性准则时, 应当忽略相对丰度小于2%的m/z, 因为它们有可能是背景或噪声, 纳入计算可能导致误剔除。强峰高概率出峰准则中, 混合谱中的强峰在参考谱中的存在性要求较弱, 其目的在于防止在重叠峰情形下排除符合要求的参考谱。经历以上过程后, 将所得参考谱集合中的所有质谱按列组装为矩阵S。
给定X和S, 若采用非负最小二乘法估计式(1)中的C, 存在“过拟合”问题。为此, 可做适当的正则化处理, 本文采取的正则化方法拟对非负最小二乘模型做一定程度的稀疏惩罚, 如下式所示:
其中‖.‖2为矩阵的2-范数, 即所有分量的平方和之平方根; ‖.‖1为矩阵的1-范数, 即所有分量的绝对值之和, 采用1-范数作为正则化项可形成稀疏结果[14], 易于勾勒数据的主要结构, 有效降低噪声敏感性和无关数据参与拟合的可能性。式(2)中, 超参数λ控制模型的稀疏程度, λ=0时模型退化为非负最小二乘模型。关于稀疏模型超参数λ, 本文设置其默认值为10, 相对丰度约为103数量级的质谱图, 稀疏惩罚程度并不高。实验结果表明, 轻微的稀疏惩罚便有助于提取待测谱的主要结构。不断增加λ的值可能更有利于抽取质谱框架性结构, 但拟合误差亦将同步提升, 同样容易导致定性错误。超参数的自适应选择方法是一项颇具挑战性的待研究内容, 从应用的角度看, 用户交互与可视化选取仍然不失为目前的最佳方案[12]。
求得稀疏优化问题式(2)的最优解C后, 可将矩阵C第i列的和作为组分i的定量估计, 亦即:
事实上, 该定量估计PAi类似于求组分i的响应强度沿保留时间的积分, 或峰面积(peak area)。j为求和的哑变量。将各组分按峰面积排名, 峰面积排名靠前的(一般取1~2种)作为可能的定性估计。
本文使用Python编程语言实现上述算法。基础数据处理使用Numpy与Pandas函数库, 稀疏优化相关部分的实现调用Scikit-Learn库中的Lasso模型, 可视化采用Matplotlib函数库。算法默认参数设置为:质谱检索阶段的阈值T=300;稀疏模型超参数默认情况下设置为λ=10。采用NIST 11质谱数据库, 共含212961张参考谱。
标准品:丁酸乙酯(ethyl butyrate, 纯度≥99%)、肉桂酸甲酯(methyl cinnamate, 纯度≥98%)、γ-十二内酯(γ-dodecalactone, 纯度≥98%)、肉桂酸正丙酯(n-propyl cinnamate, 纯度≥98%)、愈创木酚(guaiacol, 纯度≥99%)、乙基麦芽酚(ethylmaltol, 纯度≥98%)均购自Admas Reagent公司(中国); 己酸乙酯(ethyl caproate, 纯度≥99%)、正戊醇(1-pentanol, 纯度≥99%)、葵酸乙酯(ethyl decanoate, 纯度≥99%)、5-庚基二氢-2(3H)-呋喃酮(又名γ-十一内酯, 5-heptyldihydro-2(3H)-furanone, 纯度≥97%)、吲哚(indole, 纯度≥99%)、3-乙酰基吡啶(1-(3-pyridinyl)-ethanone, 纯度≥98%)、四甲基吡嗪(tetramethylpyrazine, 纯度≥98%)、甲基吡嗪(methylpyrazine, 纯度≥99%)、6-甲基-5-庚烯-2-酮(6-methyl-5-hepten-2-one, 纯度≥97.5%)和丁香酚(eugenol, 纯度≥99%)均购自比利时Acros Organics公司; 3-庚烯-2-酮(3-hepten-2-one, 纯度≥96%)、辛酸乙酯(ethyl caprylate, 纯度≥98%)、庚酸乙酯(ethyl heptanate, 纯度≥97%)、5-乙基二氢-2(3H)-呋喃酮(又名γ-己内酯, 5-ethyldihydro-2(3H)-furanone)、壬酸乙酯(ethyl nonanoate, 纯度≥95%)、正己醇(1-hexanol, 纯度≥98%)、正庚醇(1-heptanol, 纯度≥98%)、正辛醇(1-octanol, 纯度≥98%)、壬醇(1-nonanol, 纯度≥92%)、5-丙基二氢-2(3H)-呋喃酮(又名γ-庚内酯, dihydro-5-propyl-2(3H)-furanone, 纯度≥98%)、5, 6, 7, 8-四氢喹喔啉(5, 6, 7, 8-tetrahydroquinoxaline, 纯度≥98%)、4-乙基愈创木酚(4-ethylguaiacol, 纯度≥97%)、4-乙基苯酚(4-ethylphenol, 纯度≥97%)、2-乙酰基吡咯(2-acetylpyrrole, 纯度≥98%)、2-羟基-3-乙基-环戊-2-烯-1-酮(2-hydroxy-3-ethyl-2-cyclopenten-1-one, 纯度≥97%)和2-羟基-3-甲基-环戊-2-烯-1-酮(2-hydroxy-3-methyl-2-cyclopenten-1-one, 纯度≥98%)均购自日本东京化成公司; 乙基吡嗪(ethylpyrazine, 纯度≥98%)、乙酸苯乙酯(2-phenylethyl acetate, 纯度≥99%)和麦芽酚(maltol, 纯度≥97%)均购自上海国药集团; 三甲基吡嗪(trimethylpyrazine, 纯度≥98%)购自百灵威公司; 烟碱(nicotine, 纯度≥98%)为自制样品。
取上述37种常见烟用香料的标准品各5 mg于100 mL的容量瓶中, 加入乙醇至刻度线, 摇匀, 再用乙醇定容, 即配制成各标准品浓度均为5×10-5 g/mL的标准溶液。取1 mL标准溶液样品于进样小瓶中, 进样分析。
Aglient 7890气相色谱仪, 配5975型质谱检测器(美国Agilent公司); GERSTEL三合一(固相微萃取、静态顶空、溶液进样)自动进样器(德国GERSTEL公司); CP323S-OCE天平(感量0.0001 g, 德国Sartorious公司); 移液器(德国Eppendorf公司); 无水乙醇(色谱纯, 美国Dikma公司)。
色谱柱:DB-WAX毛细管柱(60.0 m×250 μm×0.25 μm)购自美国Aglient公司; 载气为高纯氦气, 流速为1.0 mL/min, 分流比为10:1;进样口温度为230 ℃; 传输线温度为250 ℃; 电离能量为70 eV; 离子源温度为230 ℃; 四极杆温度为150 ℃; 扫描范围为35~450 amu。
将上述香料混合溶液在气相色谱-质谱联用仪上进样分析, 共设计两种实验方案并获得两组不同色谱条件下的数据。
数据1(D1):起始炉温为50 ℃, 保持1 min, 再以3 ℃/min的速率升高到240 ℃, 保持2 min, 总运行时间为66.33 min。
数据2(D2):起始炉温为80 ℃, 保持5 min, 再以40 ℃/min的速率升高到220 ℃, 保持5 min, 总运行时间为13.5 min。
D1是常规实验条件下获得的数据。D2与D1相比, 由于进行了快速升温, 总运行时间约为D1的1/5, D2色谱共流出峰现象较为严重, 对后续数据处理与分析提出了挑战。安捷伦工作站MS Search与本文算法对D1的分析都能得到令人满意的结果, 但使用MS Search分析D2时, 若干纯组分的检出出现问题, 如表 1所示, “not found”表示MS Search检索排名在10名以后, “low matching rate”表示未进入前5名。从表 1可看出, 本文算法由于采用整个色谱峰数据以及稀疏优化方法, 有效地避免了上述情况的发生, 并且有效地处理了若干重叠峰情形。
质谱检索的效果应当从检索性能、检索正确率和剩余参考谱数量等方面考察。本文提出的最右端质量数索引和基峰索引旨在解决性能问题。若未使用索引技术进行粗筛, 经Numpy充分优化后, 使用WREI (weighted reference existing index)算法[9]在普通台式机上单次检索的平均执行时间约为15 s; 使用索引技术进行粗筛后, 单次检索平均计算时间约为0.5 s, 达到实时水平。
检索完成后的剩余参考谱数目影响后续分析质量。一般而言, 剩余参考谱数目较少且所得参考谱集合包含期望参考谱的方案为更佳方案。图 2为对数据D2中所有色谱峰分析完成后, 对剩余参考谱数目的统计分析。其中纵坐标为检索后剩余参考谱数目的对数值(以10为底); 横坐标为检索步骤, 步骤1和2分别对应最右端质量数符合准则和基峰符合准则, 步骤3对本文算法而言对应强峰高概率出峰准则和耐挤压性准则, 对WREI而言对应经历步骤1与步骤2后再进行参考谱加权存在指数检索算法。WREI算法参数采用文献[9]中的默认值, 即参数C=20, 检索阈值为90%。图中红色与蓝色虚线所围区域分别给定了WREI和本文算法的剩余参考谱数目的范围, 红色与蓝色实线分别给定了WREI算法与本文算法的算术平均值。
从图 2可见, 经历最右端质量数准则和基峰准则后, 参考谱数目已从2.13×105降至104量级, 经历步骤3后, 两种方法均可将参考谱数目控制在103量级。然而, WREI存在若干接近104量级的样本。从图 1左下角的内嵌Voilin图可见, WREI算法的剩余参考谱数目的平均值大于102, 而本文算法的剩余参考谱数目平均值为101左右, 且多数样本在平均数以下。
另外, WREI在默认参数下, 尚有正己醇(1-hexanol)等6种参考谱未检出, 本文算法仅麦芽酚(maltol)未检出。本文检索算法在性能和精度方面均可为后续分析提供较满意的参考谱集。
稀疏约束的目的在于一定程度上克服了非负最小二乘法的“过拟合”效应。为验证其效果, 考察了数据D2中壬醇色谱峰的解析, 如图 3与图 4所示。图 3为非负最小二乘法(即非稀疏情形, 相当于λ=0)解析结果。其中, 图 3e为实际色谱数据, 其中上方红色曲线为总离子流图(TIC), 下方的曲线簇为各m/z随保留时间的变化曲线; 图 3d为各组分的解算强度随保留时间的变化曲线, 峰面积排名前两位的质谱图分别为图 3a和图 3b, 待测质谱为图 3c。由图 3可见, 峰面积排名前两位的纯组分与待测谱比较, 均未得到较好拟合, 一般难以推断该待测质谱对应组分为壬醇, 不符合预期。实际上, 在输出结果中壬醇的峰面积排名已至第五。
图 4为稀疏模型(λ=10)解析结果。此时壬醇峰面积排名第一, 由图 4d可知壬醇的峰面积相较其他组分有明显优势。通过观察质谱图(见图 4a和图 4c)可知, 壬醇标准质谱与待测谱吻合较好。可见, 以壬醇作为定性估计较为合理, 与实际情况吻合。
稀疏惩罚的主要作用是使得最优解变得稀疏, 以提取质谱的“主要结构”, 降低相关性较低质谱强行参与拼凑的可能性。若有单个质谱与待测质谱吻合较好, 原则上推荐单个质谱作为定性结果。倘若即便实施各种程度的稀疏惩罚, 单个质谱始终未与待测质谱较好地吻合, 需考虑重叠出峰的可能性。
重叠峰解析一直是GC-MS复杂体系解析所面临的挑战, 经典算法一般基于色谱峰形进行重叠峰解析。其中, 经典方法为切线法和均线法等解析几何方法[15], 由于这类方法误差较大, 近年许多学者关注高斯峰拟合等主流数值计算方法[16-18]。基于色谱峰形的方法对严重重叠峰情形将失效[19, 20]。
D2中共流出峰现象较严重, 如图 5所示。图 5a为实际色谱数据, 视觉上可观察到单个峰包, 使用固定窗口因子分析法[21]跟踪5个主奇异值的曲线, 如图 5c所示, 可见明显高于基线的曲线仅有一条, 意味着该峰为单峰或为重叠比较严重的重叠峰。另外, 经测试各种程度的稀疏惩罚均未得到与待测谱吻合较好的单一质谱, 亦可判断其为严重重叠峰。给予一定程度的稀疏惩罚(λ=10)得到的结果如图 5b所示, 其中峰面积前两位的参考谱相比其他参考谱呈现较强优势, 恰好得到预期结果2-羟基-3-甲基-环戊-2-烯-1-酮与5-丙基二氢-2(3H)-呋喃酮。若使用非负最小二乘法, 则出现明显的“过拟合”现象, 无法得到预期结果。
为验证上述分析结果, 考察各组分的质谱图, 如图 6所示。图 6a与图 6b分别为峰面积排名第一和第二的质谱图, 图 6d为以上两张质谱图按峰面积加权求和并进行最大值对齐后的混合谱, 图 6c为待测质谱。观察4张质谱图可知, 排名第一和第二组分的质谱图都无法单独与待测谱进行较好的匹配, 但按解算结果加权求和后的混合谱与待测谱匹配较好, 并与实际结果一致。
在数据D2中, 于保留时间为7.46 min处有另一组严重重叠峰。与上述情况类似, 若使用非负最小二乘法, 依然无法得到满意的解析结果。使用稀疏模型解析, 则可有效克服“过拟合问题”, 所得纯组分为5-庚基二氢-2(3H)-呋喃酮与γ-十二内酯(见表 1), 其峰面积比为1.55:1.44。二者皆为所配标准溶液中的物质, 与预期相符。
综上所述, 本文算法在严重重叠峰解析方面较有效。
本文提出一种GC-MS数据解析算法, 该算法包括一种高效的分步检索技术以及基于该检索结果的稀疏模型解析, 以得到定性定量分析结果。实验表明, 该方法具有较好的精度和性能, 且在严重重叠峰解析中表现出良好效果。
2018,
Vol. 36 Issue (8): 772-779
PDF
2018,
Vol. 36 Issue (8): 772-779
PDF