伍毅子, 魏維偉, 楊華武*, 陳增萍, 鄒 有, 李燕春, 庹蘇行, 尹雙鳳, 鐘科軍

(1. 湖南中煙工業(yè)有限責(zé)任公司, 湖南 長沙 410014; 2. 湖南大學(xué)化學(xué)化工學(xué)院, 湖南 長沙 410082; 3. 中南大學(xué)網(wǎng)絡(luò)教育學(xué)院, 湖南 長沙 410083)

基于GC-MS的復(fù)雜體系解析是化學(xué)計量學(xué)的重要內(nèi)容。GC-MS數(shù)據(jù)解析可由如下模型刻畫[1]:

X=CST+E

(1)

其中X為氣相色譜質(zhì)譜聯(lián)用色譜儀產(chǎn)生的量測矩陣,其中每一列表示質(zhì)譜中的質(zhì)荷比(m/z)記錄,每一行表示色譜在不同保留時間的響應(yīng)強(qiáng)度。組分矩陣S的每一列為化學(xué)純組分的質(zhì)譜數(shù)據(jù),權(quán)重矩陣C中的每一列為S中相應(yīng)位置的化學(xué)純組分在不同保留時間的色譜響應(yīng)值。E為系統(tǒng)誤差矩陣。當(dāng)量測矩陣X已知,C與S均未知的情況下,解析問題為典型的黑色體系問題,一般僅可采用非監(jiān)督算法(如基于主成分分析的方法)嘗試解決[2]。

對于GC-MS數(shù)據(jù),充分利用質(zhì)譜數(shù)據(jù)庫(如NIST數(shù)據(jù)庫)的信息可顯著降低系統(tǒng)的不確定性。其方法主要分為兩類:第一類,首先進(jìn)行色譜峰識別和重疊峰解析,再利用質(zhì)譜數(shù)據(jù)庫進(jìn)行質(zhì)譜檢索,實(shí)現(xiàn)定性功能;第二類,首先根據(jù)X提供的信息,在數(shù)據(jù)庫中檢索一定數(shù)量的相關(guān)純質(zhì)譜,作為矩陣S的估計,然后使用合適的回歸算法(如非負(fù)最小二乘回歸)估計矩陣C,實(shí)現(xiàn)定性定量分析。第一類方法的主要挑戰(zhàn)在于色譜重疊峰解析,其經(jīng)典方法包括多元分辨技術(shù)[3]、小波分析[4]和神經(jīng)網(wǎng)絡(luò)[5]等,但該類方法基于色譜峰形判別,要求色譜峰具有一定的分辨率,解析嚴(yán)重重疊峰時往往失效。第二類方法的關(guān)鍵步驟為質(zhì)譜檢索和回歸算法。

關(guān)于質(zhì)譜檢索,可追溯至20世紀(jì)70年代提出的概率基礎(chǔ)匹配(PBM)算法[6],后來許多文獻(xiàn)基于該算法進(jìn)行改進(jìn),并將其廣泛應(yīng)用于商業(yè)軟件,文獻(xiàn)[7]對此類算法進(jìn)行了綜合比較,并認(rèn)為NIST MS Search中的算法綜合性能最優(yōu)。然而,此類檢索算法假設(shè)待測質(zhì)譜純度較高,對混合質(zhì)譜(色譜峰重疊時可能出現(xiàn))無法保證有效檢出各純組分[8]。為在混合質(zhì)譜假設(shè)下進(jìn)行質(zhì)譜檢索,文獻(xiàn)[9]中提出一種參考譜加權(quán)存在指數(shù)方法,該方法在許多后續(xù)GC-MS數(shù)據(jù)解析算法中作為檢索工具[8,10-13]。然而,實(shí)際計算表明,當(dāng)質(zhì)譜基線水平較高時,該算法所選參考譜數(shù)目較多,增加后續(xù)分析的不確定性。另外,使用單一指標(biāo)遍歷大規(guī)模數(shù)據(jù)庫,檢索的時間復(fù)雜度頗高。

關(guān)于回歸算法,許多文獻(xiàn)[8-11]推薦非負(fù)最小二乘法解決該回歸問題,采用方差分析或主成分分析估計待測譜中所含的純組分。但非負(fù)最小二乘法為實(shí)現(xiàn)數(shù)學(xué)上的最佳擬合往往將相關(guān)性較低的質(zhì)譜進(jìn)行強(qiáng)行拼湊,“過擬合”現(xiàn)象突出,經(jīng)常出現(xiàn)錯誤的解析結(jié)果。

鑒于上述傳統(tǒng)算法的不足,本文提出基于稀疏模型的GC-MS數(shù)據(jù)解析算法,在嚴(yán)重重疊峰解析中取得了較好效果,算法的主要特點(diǎn)如下:在質(zhì)譜檢索方面,為降低檢索的時間復(fù)雜度,提出分步檢索方案。首先利用質(zhì)譜碎片規(guī)律,結(jié)合索引技術(shù)進(jìn)行快速粗篩;然后,為進(jìn)一步降低所選參考譜集的規(guī)模,提出更為精細(xì)的強(qiáng)峰高概率出峰準(zhǔn)則和耐擠壓性準(zhǔn)則進(jìn)行參考譜剔除。在解決非負(fù)最小二乘法的“過擬合”問題方面,提出采用易于提取測量矩陣X“主要結(jié)構(gòu)”的稀疏優(yōu)化模型。

1 算法描述

如圖1所示,通過儀器獲得GC-MS數(shù)據(jù)后,首先通過有效的峰選擇算法[12,13]或用戶交互的辦法確定色譜峰,將峰起點(diǎn)與終點(diǎn)之間的數(shù)據(jù)作為量測矩陣X,將峰頂位置所對應(yīng)的質(zhì)譜作為待測混合譜實(shí)施質(zhì)譜檢索以獲得組分矩陣S,確定X與S后,通過稀疏模型解算權(quán)重矩陣C,完成解析過程。其中藍(lán)色字體標(biāo)定的部分為算法的非平凡部分,是應(yīng)當(dāng)重點(diǎn)考察的內(nèi)容。

圖 1 算法流程總覽Fig. 1 Overview of the proposed algorithm X: the measuring matrix; S: the component matrix; C: the weighting matrix; λ: the parameter for adjusting sparsity.

1.1 質(zhì)譜檢索

主要質(zhì)譜篩除步驟將質(zhì)譜數(shù)據(jù)庫中的所有質(zhì)譜作為參考譜,考察其與待測混合質(zhì)譜的相關(guān)性。進(jìn)行質(zhì)譜分析前,先對待測質(zhì)譜與參考譜均進(jìn)行規(guī)整化處理。規(guī)整化時將質(zhì)譜中最大峰的強(qiáng)度縮放至1 000,其余各峰按比例縮放。另外,為簡化處理,將非整數(shù)的m/z按四舍五入法則設(shè)定為整數(shù)。

分子離子峰和基峰是非常重要的質(zhì)譜碎片特征,可作為篩除標(biāo)準(zhǔn)。為簡化分子離子峰的認(rèn)定,本文考慮最右端質(zhì)量數(shù),即質(zhì)譜圖最右端峰簇中相對豐度最大的峰所對應(yīng)的m/z,以標(biāo)準(zhǔn)質(zhì)譜數(shù)據(jù)庫(如NIST譜庫)為基礎(chǔ),預(yù)先建立最右端質(zhì)量數(shù)索引和基峰索引,以期加快質(zhì)譜檢索速度。最右端質(zhì)量數(shù)索引將所有參考譜按最右端質(zhì)量數(shù)分類,以最右端質(zhì)量數(shù)作為類標(biāo),每個類標(biāo)存儲所有對應(yīng)參考譜在質(zhì)譜數(shù)據(jù)庫中的位置?；逅饕鎯?shù)據(jù)庫中所有參考譜的基峰位置,存儲形式為鍵值對,以參考譜位置為鍵,以基峰位置為值。下文詳述具體篩除步驟或準(zhǔn)則。

1.1.1最右端質(zhì)量數(shù)符合準(zhǔn)則

考察待測混合質(zhì)譜中的任一有效m/z,通過查詢最右端質(zhì)量數(shù)索引,可得相應(yīng)m/z對應(yīng)的所有參考譜。然后合并所有有效m/z的對應(yīng)參考譜列表,即得所需候選質(zhì)譜集。

1.1.2基峰符合準(zhǔn)則

基于上一步所得候選質(zhì)譜集,考察其中的每個參考譜,通過查詢基峰索引獲得其基峰位置(m/z)?？紤]待測混合質(zhì)譜中相應(yīng)位置的出峰強(qiáng)度,若該強(qiáng)度低于某閾值T(默認(rèn)T=300),則將所考慮參考譜從候選質(zhì)譜集中剔除。

1.1.3強(qiáng)峰高概率出峰準(zhǔn)則

設(shè)定相對豐度超過一定閾值的峰為強(qiáng)峰,對參考譜中的任意強(qiáng)峰,若混合譜中相應(yīng)位置的出峰相對豐度與該峰的比值低于閾值Q=T/1 000,則標(biāo)記為異常,若標(biāo)記為異常的強(qiáng)峰數(shù)目超過2,則剔除所考察的參考譜。反過來,考察混合譜中的強(qiáng)峰(強(qiáng)度高于T),若參考譜相應(yīng)位置的出峰強(qiáng)度與待測譜的強(qiáng)度之比低于Q,則標(biāo)記為異常,若標(biāo)記為異常的峰數(shù)目超過混合譜中強(qiáng)峰總數(shù)的一半,則剔除所考察的參考譜。

1.1.4耐擠壓性準(zhǔn)則

將參考譜與混合譜對齊后,對參考譜各個m/z進(jìn)行同比例縮放,直至參考譜在任何m/z處的出峰都低于待測混合質(zhì)譜相應(yīng)位置的出峰。將參考譜的任一有效m/z在壓縮后與壓縮前的豐度比定義為擠壓比例,若擠壓比例小于閾值Q,則剔除所考察參考譜。

實(shí)施強(qiáng)峰高概率出峰準(zhǔn)則和耐擠壓性準(zhǔn)則時,應(yīng)當(dāng)忽略相對豐度小于2%的m/z,因為它們有可能是背景或噪聲,納入計算可能導(dǎo)致誤剔除。強(qiáng)峰高概率出峰準(zhǔn)則中,混合譜中的強(qiáng)峰在參考譜中的存在性要求較弱,其目的在于防止在重疊峰情形下排除符合要求的參考譜。經(jīng)歷以上過程后,將所得參考譜集合中的所有質(zhì)譜按列組裝為矩陣S。

1.2 稀疏模型

給定X和S,若采用非負(fù)最小二乘法估計式(1)中的C,存在“過擬合”問題。為此,可做適當(dāng)?shù)恼齽t化處理,本文采取的正則化方法擬對非負(fù)最小二乘模型做一定程度的稀疏懲罰,如下式所示:

min‖CST-X‖22+λ‖C‖1s.t.C≥0

(2)

其中‖.‖2為矩陣的2-范數(shù),即所有分量的平方和之平方根;‖.‖1為矩陣的1-范數(shù),即所有分量的絕對值之和,采用1-范數(shù)作為正則化項可形成稀疏結(jié)果[14],易于勾勒數(shù)據(jù)的主要結(jié)構(gòu),有效降低噪聲敏感性和無關(guān)數(shù)據(jù)參與擬合的可能性。式(2)中,超參數(shù)λ控制模型的稀疏程度,λ=0時模型退化為非負(fù)最小二乘模型。關(guān)于稀疏模型超參數(shù)λ,本文設(shè)置其默認(rèn)值為10,相對豐度約為103數(shù)量級的質(zhì)譜圖,稀疏懲罰程度并不高。實(shí)驗結(jié)果表明,輕微的稀疏懲罰便有助于提取待測譜的主要結(jié)構(gòu)。不斷增加λ的值可能更有利于抽取質(zhì)譜框架性結(jié)構(gòu),但擬合誤差亦將同步提升,同樣容易導(dǎo)致定性錯誤。超參數(shù)的自適應(yīng)選擇方法是一項頗具挑戰(zhàn)性的待研究內(nèi)容,從應(yīng)用的角度看,用戶交互與可視化選取仍然不失為目前的最佳方案[12]。

求得稀疏優(yōu)化問題式(2)的最優(yōu)解C后,可將矩陣C第i列的和作為組分i的定量估計,亦即:

PAi=∑jCji

(3)

事實(shí)上,該定量估計PAi類似于求組分i的響應(yīng)強(qiáng)度沿保留時間的積分,或峰面積(peak area)。j為求和的啞變量。將各組分按峰面積排名,峰面積排名靠前的(一般取1～2種)作為可能的定性估計。

1.3 算法實(shí)現(xiàn)

本文使用Python編程語言實(shí)現(xiàn)上述算法?；A(chǔ)數(shù)據(jù)處理使用Numpy與Pandas函數(shù)庫,稀疏優(yōu)化相關(guān)部分的實(shí)現(xiàn)調(diào)用Scikit-Learn庫中的Lasso模型,可視化采用Matplotlib函數(shù)庫。算法默認(rèn)參數(shù)設(shè)置為:質(zhì)譜檢索階段的閾值T=300;稀疏模型超參數(shù)默認(rèn)情況下設(shè)置為λ=10。采用NIST 11質(zhì)譜數(shù)據(jù)庫,共含212 961張參考譜。

2 實(shí)驗部分

2.1 香料混合溶液制備

標(biāo)準(zhǔn)品:丁酸乙酯(ethyl butyrate,純度≥99% )、肉桂酸甲酯(methyl cinnamate,純度≥98% )、γ-十二內(nèi)酯(γ-dodecalactone,純度≥98% )、肉桂酸正丙酯(n-propyl cinnamate,純度≥98% )、愈創(chuàng)木酚(guaiacol,純度≥99% )、乙基麥芽酚(ethylmaltol,純度≥98% )均購自Admas Reagent公司(中國);己酸乙酯(ethyl caproate,純度≥99% )、正戊醇(1-pentanol,純度≥99% )、葵酸乙酯(ethyl decanoate,純度≥99% )、5-庚基二氫-2(3H)-呋喃酮(又名γ-十一內(nèi)酯,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)-呋喃酮(又名γ-己內(nèi)酯,5-ethyldihydro-2(3H)-furanone)、壬酸乙酯(ethyl nonanoate,純度≥95% )、正己醇(1-hexanol,純度≥98% )、正庚醇(1-heptanol,純度≥98% )、正辛醇(1-octanol,純度≥98% )、壬醇(1-nonanol,純度≥92% )、5-丙基二氫-2(3H)-呋喃酮(又名γ-庚內(nèi)酯,dihydro-5-propyl-2(3H)-furanone,純度≥98% )、5,6,7,8-四氫喹喔啉(5,6,7,8-tetrahydroquinoxaline,純度≥98% )、4-乙基愈創(chuàng)木酚(4-ethylguaiacol,純度≥97% )、4-乙基苯酚(4-ethylphenol,純度≥97% )、2-乙?；量?2-acetylpyrrole,純度≥98% )、2-羥基-3-乙基-環(huán)戊-2-烯-1-酮(2-hydroxy-3-ethyl-2-cyclopenten-1-one,純度≥97% )和2-羥基-3-甲基-環(huán)戊-2-烯-1-酮(2-hydroxy-3-methyl-2-cyclopenten-1-one,純度≥98% )均購自日本東京化成公司;乙基吡嗪(ethylpyrazine,純度≥98% )、乙酸苯乙酯(2-phenylethyl acetate,純度≥99% )和麥芽酚(maltol,純度≥97% )均購自上海國藥集團(tuán);三甲基吡嗪(trimethylpyrazine,純度≥98% )購自百靈威公司;煙堿(nicotine,純度≥98% )為自制樣品。

取上述37種常見煙用香料的標(biāo)準(zhǔn)品各5 mg于100 mL的容量瓶中,加入乙醇至刻度線,搖勻,再用乙醇定容,即配制成各標(biāo)準(zhǔn)品濃度均為5×10-5g/mL的標(biāo)準(zhǔn)溶液。取1 mL標(biāo)準(zhǔn)溶液樣品于進(jìn)樣小瓶中,進(jìn)樣分析。

2.2 儀器、試劑及實(shí)驗條件

Aglient 7890氣相色譜儀,配5975型質(zhì)譜檢測器(美國Agilent公司); GERSTEL三合一(固相微萃取、靜態(tài)頂空、溶液進(jìn)樣)自動進(jìn)樣器(德國GERSTEL公司); CP323S-OCE天平(感量0.000 1 g,德國Sartorious公司);移液器(德國Eppendorf公司);無水乙醇(色譜純,美國Dikma公司)。

色譜柱:DB-WAX毛細(xì)管柱(60.0 m×250 μm×0.25 μm)購自美國Aglient公司;載氣為高純氦氣,流速為1.0 mL/min,分流比為10∶1;進(jìn)樣口溫度為230 ℃;傳輸線溫度為250 ℃;電離能量為70 eV;離子源溫度為230 ℃;四極桿溫度為150 ℃;掃描范圍為35～450 amu。

2.3 實(shí)驗設(shè)計

將上述香料混合溶液在氣相色譜-質(zhì)譜聯(lián)用儀上進(jìn)樣分析,共設(shè)計兩種實(shí)驗方案并獲得兩組不同色譜條件下的數(shù)據(jù)。

數(shù)據(jù)1(D1):起始爐溫為50 ℃,保持1 min,再以3 ℃/min的速率升高到240 ℃,保持2 min,總運(yùn)行時間為66.33 min。

數(shù)據(jù)2(D2):起始爐溫為80 ℃,保持5 min,再以40 ℃/min的速率升高到220 ℃,保持5 min,總運(yùn)行時間為13.5 min。

3 結(jié)果與討論

D1是常規(guī)實(shí)驗條件下獲得的數(shù)據(jù)。D2與D1相比,由于進(jìn)行了快速升溫,總運(yùn)行時間約為D1的1/5, D2色譜共流出峰現(xiàn)象較為嚴(yán)重,對后續(xù)數(shù)據(jù)處理與分析提出了挑戰(zhàn)。安捷倫工作站MS Search與本文算法對D1的分析都能得到令人滿意的結(jié)果,但使用MS Search分析D2時,若干純組分的檢出出現(xiàn)問題,如表1所示,“not found”表示MS Search檢索排名在10名以后,“l(fā)ow matching rate”表示未進(jìn)入前5名。從表1可看出,本文算法由于采用整個色譜峰數(shù)據(jù)以及稀疏優(yōu)化方法,有效地避免了上述情況的發(fā)生,并且有效地處理了若干重疊峰情形。

表 1 數(shù)據(jù)2(D2)中使用MS Search分析異常的純組分Table 1 Pure components with resolution failure using MS Search in data 2 (D2)

3.1 質(zhì)譜檢索效果分析

質(zhì)譜檢索的效果應(yīng)當(dāng)從檢索性能、檢索正確率和剩余參考譜數(shù)量等方面考察。本文提出的最右端質(zhì)量數(shù)索引和基峰索引旨在解決性能問題。若未使用索引技術(shù)進(jìn)行粗篩,經(jīng)Numpy充分優(yōu)化后,使用WREI (weighted reference existing index)算法[9]在普通臺式機(jī)上單次檢索的平均執(zhí)行時間約為15 s;使用索引技術(shù)進(jìn)行粗篩后,單次檢索平均計算時間約為0.5 s,達(dá)到實(shí)時水平。

檢索完成后的剩余參考譜數(shù)目影響后續(xù)分析質(zhì)量。一般而言,剩余參考譜數(shù)目較少且所得參考譜集合包含期望參考譜的方案為更佳方案。圖2為對數(shù)據(jù)D2中所有色譜峰分析完成后,對剩余參考譜數(shù)目的統(tǒng)計分析。其中縱坐標(biāo)為檢索后剩余參考譜數(shù)目的對數(shù)值(以10為底);橫坐標(biāo)為檢索步驟,步驟1和2分別對應(yīng)最右端質(zhì)量數(shù)符合準(zhǔn)則和基峰符合準(zhǔn)則,步驟3對本文算法而言對應(yīng)強(qiáng)峰高概率出峰準(zhǔn)則和耐擠壓性準(zhǔn)則,對WREI而言對應(yīng)經(jīng)歷步驟1與步驟2后再進(jìn)行參考譜加權(quán)存在指數(shù)檢索算法。WREI算法參數(shù)采用文獻(xiàn)[9]中的默認(rèn)值,即參數(shù)C=20,檢索閾值為90% 。圖中紅色與藍(lán)色虛線所圍區(qū)域分別給定了WREI和本文算法的剩余參考譜數(shù)目的范圍,紅色與藍(lán)色實(shí)線分別給定了WREI算法與本文算法的算術(shù)平均值。

從圖2可見,經(jīng)歷最右端質(zhì)量數(shù)準(zhǔn)則和基峰準(zhǔn)則后,參考譜數(shù)目已從2.13×105降至104量級,經(jīng)歷步驟3后,兩種方法均可將參考譜數(shù)目控制在103量級。然而,WREI存在若干接近104量級的樣本。從圖1左下角的內(nèi)嵌Voilin圖可見,WREI算法的剩余參考譜數(shù)目的平均值大于102,而本文算法的剩余參考譜數(shù)目平均值為101左右,且多數(shù)樣本在平均數(shù)以下。

另外,WREI在默認(rèn)參數(shù)下,尚有正己醇(1-hexanol)等6種參考譜未檢出,本文算法僅麥芽酚(maltol)未檢出。本文檢索算法在性能和精度方面均可為后續(xù)分析提供較滿意的參考譜集。

圖 2 本文算法與WREI算法的剩余參考譜數(shù)目對比Fig. 2 Comparison of the numbers of residual reference spectra between this method and weighted reference existing index (WREI) algorithm Insert: distribution of the number of the residual reference spectra (violin diagram for Step 3). N: the number of the residual reference spectra.

圖 3 非稀疏情形下數(shù)據(jù)D2中保留時間4.88 min處的 色譜峰解析Fig. 3 Analysis of the chromatographic peak at 4.88 min in data D2 in the case of non sparsity a. mass spectrum of pure 1; b. mass spectrum of pure 2; c. mass spectrum of the peak top in the mixture; d. regression coefficients of pures, each curve presents the variation of regression coefficients of some compound; e. original chromatographic data, containing TIC curve (red) and m/z curves.

圖 4 稀疏情形下數(shù)據(jù)D2中保留時間4.88 min處的 色譜峰解析Fig. 4 Analysis results of the chromatographic peak at 4.88 min in data D2 in the case of sparsity a. mass spectrum of pure 1; b. mass spectrum of pure 2; c. mass spectrum of the peak top in the mixture; d. regression coefficients of pures, each curve presents the variation of regression coefficients of some compound; e. original chromatographic data, containing TIC curve (red) and m/z curves.

3.2 稀疏模型的優(yōu)勢

稀疏約束的目的在于一定程度上克服了非負(fù)最小二乘法的“過擬合”效應(yīng)。為驗證其效果,考察了數(shù)據(jù)D2中壬醇色譜峰的解析,如圖3與圖4所示。圖3為非負(fù)最小二乘法(即非稀疏情形,相當(dāng)于λ=0)解析結(jié)果。其中,圖3e為實(shí)際色譜數(shù)據(jù),其中上方紅色曲線為總離子流圖(TIC),下方的曲線簇為各m/z隨保留時間的變化曲線;圖3d為各組分的解算強(qiáng)度隨保留時間的變化曲線,峰面積排名前兩位的質(zhì)譜圖分別為圖3a和圖3b,待測質(zhì)譜為圖3c。由圖3可見,峰面積排名前兩位的純組分與待測譜比較,均未得到較好擬合,一般難以推斷該待測質(zhì)譜對應(yīng)組分為壬醇,不符合預(yù)期。實(shí)際上,在輸出結(jié)果中壬醇的峰面積排名已至第五。

圖4為稀疏模型(λ=10)解析結(jié)果。此時壬醇峰面積排名第一,由圖4d可知壬醇的峰面積相較其他組分有明顯優(yōu)勢。通過觀察質(zhì)譜圖(見圖4a和圖4c)可知,壬醇標(biāo)準(zhǔn)質(zhì)譜與待測譜吻合較好。可見,以壬醇作為定性估計較為合理,與實(shí)際情況吻合。

稀疏懲罰的主要作用是使得最優(yōu)解變得稀疏,以提取質(zhì)譜的“主要結(jié)構(gòu)”,降低相關(guān)性較低質(zhì)譜強(qiáng)行參與拼湊的可能性。若有單個質(zhì)譜與待測質(zhì)譜吻合較好,原則上推薦單個質(zhì)譜作為定性結(jié)果。倘若即便實(shí)施各種程度的稀疏懲罰,單個質(zhì)譜始終未與待測質(zhì)譜較好地吻合,需考慮重疊出峰的可能性。

3.3 在嚴(yán)重重疊峰解析中的應(yīng)用

重疊峰解析一直是GC-MS復(fù)雜體系解析所面臨的挑戰(zhàn),經(jīng)典算法一般基于色譜峰形進(jìn)行重疊峰解析。其中,經(jīng)典方法為切線法和均線法等解析幾何方法[15],由于這類方法誤差較大,近年許多學(xué)者關(guān)注高斯峰擬合等主流數(shù)值計算方法[16-18]。基于色譜峰形的方法對嚴(yán)重重疊峰情形將失效[19,20]。

圖 5 數(shù)據(jù)D2中保留時間5.493 min處的嚴(yán)重重疊峰解析Fig. 5 Analysis of the severe overlapping peak at 5.493 min in data D2 a. original chromatographic data, contains TIC curve and m/z curves; b. regression coefficients of the pures, each curve presents the variation of regression coefficients of some compound; c. singular values of the original chromatographic matrix, contains 5 singular curves.

D2中共流出峰現(xiàn)象較嚴(yán)重,如圖5所示。圖5a為實(shí)際色譜數(shù)據(jù),視覺上可觀察到單個峰包,使用固定窗口因子分析法[21]跟蹤5個主奇異值的曲線,如圖5c所示,可見明顯高于基線的曲線僅有一條,意味著該峰為單峰或為重疊比較嚴(yán)重的重疊峰。另外,經(jīng)測試各種程度的稀疏懲罰均未得到與待測譜吻合較好的單一質(zhì)譜,亦可判斷其為嚴(yán)重重疊峰。給予一定程度的稀疏懲罰(λ=10)得到的結(jié)果如圖5b所示,其中峰面積前兩位的參考譜相比其他參考譜呈現(xiàn)較強(qiáng)優(yōu)勢,恰好得到預(yù)期結(jié)果2-羥基-3-甲基-環(huán)戊-2-烯-1-酮與5-丙基二氫-2(3H)-呋喃酮。若使用非負(fù)最小二乘法,則出現(xiàn)明顯的“過擬合”現(xiàn)象,無法得到預(yù)期結(jié)果。

圖 6 2-羥基-3-甲基-環(huán)戊-2-烯-1-酮與5-丙基二氫-2(3H)- 呋喃酮峰面積加權(quán)混合譜與待測譜對比Fig. 6 Comparison between the peak area weighted mixture of the two pures (2-hydroxy-3-methyl-2-cyclopenten-1-one and dihydro-5-propyl-2(3H)-furanone) and the mass spectrum to be resolved a. mass spectrum of pure 1; b. mass spectrum of pure 2; c. original chromatographic data; d. peak area weighted mixture mass spectrum of pure 1 and pure 2.

為驗證上述分析結(jié)果,考察各組分的質(zhì)譜圖,如圖6所示。圖6a與圖6b分別為峰面積排名第一和第二的質(zhì)譜圖,圖6d為以上兩張質(zhì)譜圖按峰面積加權(quán)求和并進(jìn)行最大值對齊后的混合譜,圖6c為待測質(zhì)譜。觀察4張質(zhì)譜圖可知,排名第一和第二組分的質(zhì)譜圖都無法單獨(dú)與待測譜進(jìn)行較好的匹配,但按解算結(jié)果加權(quán)求和后的混合譜與待測譜匹配較好,并與實(shí)際結(jié)果一致。

在數(shù)據(jù)D2中,于保留時間為7.46 min處有另一組嚴(yán)重重疊峰。與上述情況類似,若使用非負(fù)最小二乘法,依然無法得到滿意的解析結(jié)果。使用稀疏模型解析,則可有效克服“過擬合問題”,所得純組分為5-庚基二氫-2(3H)-呋喃酮與γ-十二內(nèi)酯(見表1),其峰面積比為1.55∶1.44。二者皆為所配標(biāo)準(zhǔn)溶液中的物質(zhì),與預(yù)期相符。

綜上所述,本文算法在嚴(yán)重重疊峰解析方面較有效。

4 結(jié)論

本文提出一種GC-MS數(shù)據(jù)解析算法,該算法包括一種高效的分步檢索技術(shù)以及基于該檢索結(jié)果的稀疏模型解析,以得到定性定量分析結(jié)果。實(shí)驗表明,該方法具有較好的精度和性能,且在嚴(yán)重重疊峰解析中表現(xiàn)出良好效果。