一種新的基于多重液相色譜-質(zhì)譜實(shí)驗(yàn)肽信號(hào)峰形相似性的校準(zhǔn)算法

崔 健，董曉睿，商 凱，陳 強(qiáng)，祁 鑫，崔 浩

(1.中國石油大學(xué) 勝利學(xué)院 信息技術(shù)系，山東 東營 257016；2.中國石油大學(xué) 計(jì)算機(jī)與通信工程學(xué)院，山東 青島 266580)

液相色譜-質(zhì)譜(LC-MS)是發(fā)現(xiàn)并分析生物標(biāo)志物中復(fù)雜肽信號(hào)的關(guān)鍵技術(shù)，其中對實(shí)驗(yàn)譜中的肽信號(hào)進(jìn)行檢測和量化至關(guān)重要。理論上相同樣本多次重復(fù)實(shí)驗(yàn)得到的譜圖是一致的，即同種肽鏈應(yīng)在不同譜圖的相同位置(相同LC時(shí)間與m/z值)產(chǎn)生相同信號(hào)。但由于實(shí)驗(yàn)誤差，多次重復(fù)實(shí)驗(yàn)譜圖會(huì)存在較大差異，需對譜數(shù)據(jù)進(jìn)行校準(zhǔn)，減小誤差。

圖1 肽鏈“KVEDMMK”生成的XICsFig.1 XICs generated by “KVEDMMK”A.XICs of peptide “KVEDMMK”in data 1;B.XICs of peptide “KVEDMMK”in data 2

圖2 肽鏈“AGGPTTPLSPTR”的相關(guān)峰信號(hào)匹配Fig.2 Corresponding peak of “AGGPTTPLSPTR” A.LC peak between 2 800-3 400 of peptide “AGGPTTPLSPTR” in data 1;B.LC peak between 2 800-3 400 of peptide “AGGPTTPLSPTR” in data 2

根據(jù)液相色譜-二級(jí)質(zhì)譜(LC-MS/MS)實(shí)驗(yàn)標(biāo)識(shí)的LC峰時(shí)間位置，通常使用翹曲函數(shù)(Warping function)對時(shí)間特征進(jìn)行校準(zhǔn)。目前，采用翹曲函數(shù)對LC時(shí)間軸校準(zhǔn)，通常先計(jì)算肽鏈的m/z值，然后固定m/z值對整個(gè)時(shí)間譜圖進(jìn)行匹配。但由于時(shí)間差產(chǎn)生的隨機(jī)性，該方法并不能完全校準(zhǔn)。有研究者提出了基于翹曲函數(shù)的改進(jìn)算法，如2002年Nielsen等[1]提出的相關(guān)優(yōu)化翹曲函數(shù)算法(Correlation optimized warping,COW)；2004年Eilers[2]提出的參數(shù)時(shí)間翹曲函數(shù)算法(Parametric time warping,PTW)；2006年van Nederkassel等[3]提出的半?yún)?shù)時(shí)間翹曲函數(shù)算法(Semi-parametric time warping,STW)及Jaitly等[4]提出的液相質(zhì)譜數(shù)據(jù)翹曲函數(shù)法(lcmswarp)等。Voss等[5]提出了一種將相關(guān)特征峰對和整體時(shí)間校正相結(jié)合的算法，該法關(guān)注同時(shí)校正多重實(shí)驗(yàn)數(shù)據(jù)，但對數(shù)據(jù)的處理效果比OpenMS軟件[6]略差。此外，當(dāng)實(shí)驗(yàn)樣本比較復(fù)雜時(shí)，在一個(gè)m/z值相同的提取離子色譜圖(XICs)中會(huì)有多個(gè)LC峰出現(xiàn)在一個(gè)狹窄的LC時(shí)間窗口，極可能導(dǎo)致相應(yīng)特征峰的錯(cuò)誤識(shí)別。圖1顯示了肽鏈“KVEDMMK”在本文所處理的實(shí)驗(yàn)數(shù)據(jù)1和2中產(chǎn)生的XICs，其LC峰信號(hào)充滿噪聲，主峰附近也分布很多噪聲峰，即使用OpenMS或Msinspect[7]等軟件進(jìn)行多個(gè)數(shù)據(jù)集處理，也無法避免此類問題。

時(shí)間誤差的隨機(jī)性及噪聲會(huì)導(dǎo)致匹配結(jié)果準(zhǔn)確性降低。以肽鏈“AGGPTTPLSPTR”為例(圖2)，在數(shù)據(jù)1中檢測到多個(gè)肽鏈信號(hào)(Peak1～6)，與MS/MS檢測結(jié)果對比，確認(rèn)Peak3為真正的肽信號(hào)峰；數(shù)據(jù)2中由MS/MS檢測到的真實(shí)肽信號(hào)已標(biāo)出。以數(shù)據(jù)2為基準(zhǔn)匹配數(shù)據(jù)1中的真實(shí)信號(hào)，即在數(shù)據(jù)1中的6個(gè)信號(hào)中找到與數(shù)據(jù)2中真實(shí)信號(hào)相匹配的信號(hào)Peak3。而從時(shí)間間隔上顯示匹配結(jié)果為Peak2，而非Peak3。因此，除時(shí)間特征外，還需引入其他特征提高校準(zhǔn)匹配的準(zhǔn)確性。

目前，僅MS/MS均能識(shí)別出的肽鏈與LC-MS峰重疊的一小部分可以使用Quil[8]、Proteinquant[9]、Msinspect[10]、OpenMS[11-12]和Superhirn[13]等軟件進(jìn)行重復(fù)實(shí)驗(yàn)數(shù)據(jù)量化。MaxQuant軟件[14-16]可以大大提高量化范圍，是因?yàn)橛蒑S/MS檢測的肽鏈可以至少量化1次，但在所有數(shù)據(jù)集中可以同時(shí)量化的總肽數(shù)量有限，只能是多數(shù)據(jù)MS/MS肽鏈信號(hào)的交集。這導(dǎo)致多次重復(fù)實(shí)驗(yàn)數(shù)據(jù)量化的覆蓋率較低。

針對兩個(gè)重復(fù)LC-MS實(shí)驗(yàn)數(shù)據(jù)，本文采用普通的區(qū)間檢測方法，選取多次重復(fù)實(shí)驗(yàn)數(shù)據(jù)中均被MS/MS檢出的肽鏈信號(hào)作為訓(xùn)練數(shù)據(jù)集。以MS/MS檢測到的肽鏈的m/z值及LC時(shí)間值為真實(shí)值(Ground truth)，訓(xùn)練數(shù)據(jù)集中所有肽鏈在兩個(gè)數(shù)據(jù)中均具有真實(shí)值。從訓(xùn)練數(shù)據(jù)集中隨機(jī)選取部分肽鏈，基于其真實(shí)值建立統(tǒng)計(jì)學(xué)習(xí)模型，訓(xùn)練數(shù)據(jù)集中剩余部分作為測試序列進(jìn)行模型測試，以模型給出的最高分值作為匹配結(jié)果，再比對真實(shí)值，計(jì)算百分比作為檢測模型的準(zhǔn)確性(Accuracy)。最后，僅在1個(gè)數(shù)據(jù)中被MS/MS檢測到的肽鏈通過模型匹配其在另外數(shù)據(jù)中的相關(guān)區(qū)間(無MS/MS檢測結(jié)果)，提升校準(zhǔn)后肽信號(hào)的覆蓋率。

圖3 數(shù)據(jù)1和數(shù)據(jù)2的MS/MS檢測肽鏈信息文氏圖Fig.3 Venn diagram of MS/MS peptides information of data 1 and data 2

1 實(shí)驗(yàn)部分

1.1 數(shù)據(jù)來源

本文處理的數(shù)據(jù)由RCMI Proteomics and Protein Biomarkers Cores at UTSA實(shí)驗(yàn)室提供，為經(jīng)過LTQ Orbitrap Velos儀器處理的TAGE腫瘤數(shù)據(jù)(腫瘤樣本的多次重復(fù)實(shí)驗(yàn)數(shù)據(jù))，選取2組數(shù)據(jù)(數(shù)據(jù)1與數(shù)據(jù)2)進(jìn)行分析。每個(gè)數(shù)據(jù)分為Level1和Level2。數(shù)據(jù)中每個(gè)數(shù)據(jù)點(diǎn)包括3個(gè)坐標(biāo)時(shí)間值、質(zhì)荷比值、強(qiáng)度值。數(shù)據(jù)1中Level1的有用數(shù)據(jù)(強(qiáng)度值不為0)為11487個(gè)，Level2的有用數(shù)據(jù)(強(qiáng)度值不為0)為58636個(gè)。數(shù)據(jù)2中Level1的有用數(shù)據(jù)(強(qiáng)度值不為0)為11446個(gè)，Level2的有用數(shù)據(jù)(強(qiáng)度值不為0)為59573個(gè)。數(shù)據(jù)1中被MS/MS檢測到的肽鏈為1944個(gè)，數(shù)據(jù)2中被MS/MS檢測到的肽鏈為1603個(gè)，交集為700個(gè)，并集為2847個(gè)(圖3)。

1.2 數(shù)據(jù)處理

1.2.1數(shù)據(jù)預(yù)處理如圖4所示，根據(jù)實(shí)驗(yàn)1和實(shí)驗(yàn)2的MS/MS信息表，生成MS/MS肽信號(hào)合集，并計(jì)算相應(yīng)肽信號(hào)的質(zhì)荷比(m/z值)。然后在數(shù)據(jù)1和數(shù)據(jù)2的Level 1數(shù)據(jù)中，分別計(jì)算肽鏈m/z值(前后各取20×10-6寬度)下的LC譜圖，獲取全時(shí)間段的XICs；在全時(shí)段XICs上進(jìn)行區(qū)間檢測，在數(shù)據(jù)1與數(shù)據(jù)2中分別獲得的區(qū)間段均為候選信號(hào)區(qū)間。

1.2.2訓(xùn)練數(shù)據(jù)集生成具備可測試的真實(shí)值是選取訓(xùn)練數(shù)據(jù)集的前提。本研究以MS/MS檢測到的肽信號(hào)的m/z值與時(shí)間值為真實(shí)可靠的值。首先，選取圖3中的交集部分作為生成訓(xùn)練數(shù)據(jù)集的基礎(chǔ)(共700個(gè)肽鏈)，將預(yù)處理區(qū)間檢測后包含MS/MS時(shí)間點(diǎn)的肽鏈選作訓(xùn)練數(shù)據(jù)集合(共599個(gè)肽鏈)。再在訓(xùn)練數(shù)據(jù)集中隨機(jī)選取一半作為訓(xùn)練序列，產(chǎn)生時(shí)間差統(tǒng)計(jì)學(xué)習(xí)模型以及峰形相似性模型；剩余部分作為測試序列，測試模型匹配結(jié)果的準(zhǔn)確性(以MS/MS檢測值作為真實(shí)數(shù)據(jù)比對)。

1.2.3統(tǒng)計(jì)學(xué)習(xí)模型生成基于以下兩個(gè)假設(shè)建立模型:①同一肽鏈在重復(fù)實(shí)驗(yàn)中產(chǎn)生信號(hào)區(qū)間的位置(包括m/z與時(shí)間)理論上一致，不同種肽鏈產(chǎn)生的信號(hào)位置有差別；②同一肽鏈重復(fù)實(shí)驗(yàn)產(chǎn)生的信號(hào)形狀理論上一致，不同種肽鏈產(chǎn)生的信號(hào)形狀有差別。在訓(xùn)練數(shù)據(jù)集中選取訓(xùn)練肽信號(hào)k個(gè)，以圖2為例:肽鏈“AGGPTTPLSPTR”在數(shù)據(jù)1中的真實(shí)信號(hào)峰(Peak 3)與數(shù)據(jù)2中的真實(shí)信號(hào)峰為相關(guān)信號(hào)峰對，數(shù)據(jù)1中其他信號(hào)峰(除去Peak 3)與數(shù)據(jù)2中真實(shí)信號(hào)峰為非相關(guān)信號(hào)峰對，分別計(jì)算相關(guān)信號(hào)峰對的時(shí)間差與峰形相似性，以及非相關(guān)信號(hào)峰對的時(shí)間差與峰形相似性。時(shí)間差即為區(qū)間最高值的時(shí)間差值，峰形相似性為計(jì)算兩個(gè)信號(hào)的線性回歸決定系數(shù)r2的值。r2反映了兩個(gè)數(shù)列的相似程度，如數(shù)列x、y的r2值反映了數(shù)列y的變化可用數(shù)列x的變化來解釋的百分比，計(jì)算公式如下:

r2=SSreg/SStot=1-SSres/SStot

其中，SStot為總平方和，SSreg為回歸平方和，SSres為殘差平方和。r2結(jié)果在0～1之間，SStot在數(shù)據(jù)確定后始終為固定值。估計(jì)的準(zhǔn)確性越低，則SSres越大，r2越接近0；反之，則r2越接近1，即峰形越相似r2值越接近1。

時(shí)間差特征統(tǒng)計(jì)特性如圖5所示。時(shí)間差直方圖基本符合正態(tài)分布，相關(guān)信號(hào)峰對的時(shí)間差相對非相關(guān)信號(hào)峰對的時(shí)間差更集中(圖5A、B)。針對時(shí)間差樣本，采用最大似然估計(jì)生成兩個(gè)正態(tài)分布模型的參數(shù)。用相關(guān)峰對的時(shí)間差樣本t估計(jì)正態(tài)分布模型f(Δt|t)的參數(shù)μ和σ:

得到相關(guān)信號(hào)峰對的時(shí)間差模型f(Δt|μ,σ2):

同理使用非相關(guān)峰對時(shí)間差樣本得到非相關(guān)信號(hào)峰對的時(shí)間差模型。圖5C、D為相關(guān)和非相關(guān)信號(hào)峰對的正態(tài)分布時(shí)間差模型。

峰形相似性特征統(tǒng)計(jì)特性見圖 6，相關(guān)信號(hào)峰對的相似性集中在0.7以上(圖6A)，而非相關(guān)信號(hào)峰對的相似性比較分散(圖6B)。采用gamma分布進(jìn)行擬合:

gamma分布具有兩個(gè)參數(shù)k和θ，采用matlab中g(shù)amfit( )函數(shù)進(jìn)行g(shù)amma分布參數(shù)的極大似然估計(jì)得數(shù)值解，得到相關(guān)信號(hào)峰對的相似性模型和非相關(guān)信號(hào)峰對的相似性模型，峰形相似性模型的區(qū)分度較明顯(圖6C)。

1.2.4LC峰匹配校準(zhǔn)基于圖3的交集數(shù)據(jù)建立時(shí)間差與峰形相似性的統(tǒng)計(jì)模型，并測試模型有效性后，將統(tǒng)計(jì)模型用于圖3中去交集部分肽鏈的校準(zhǔn)匹配，即根據(jù)僅在1個(gè)數(shù)據(jù)中由MS/MS檢測到的肽鏈信號(hào)，匹配其在另一數(shù)據(jù)中的相關(guān)信號(hào)。最終由MS/MS檢測到的每一個(gè)肽鏈均能在數(shù)據(jù)1和2中匹配到相關(guān)信號(hào)峰對。

圖6 相關(guān)峰與非相關(guān)峰峰形相似性模型Fig.6 Peak shape similarity of corresponding and non-corresponding peak modelsA.histogram of corresponding peak similarity(相關(guān)峰對峰形相似性直方圖)；B.histogram of non-corresponding peak similarity(非相關(guān)峰對峰形相似性直方圖)；C.model of peak similarity(峰形相似性模型)

2 結(jié)果與討論

2.1 結(jié)果分析

2.1.1模型測試結(jié)果對模型進(jìn)行10次測試，每次均從訓(xùn)練數(shù)據(jù)集中隨機(jī)選取300個(gè)肽鏈信號(hào)作為訓(xùn)練，剩余299個(gè)肽鏈信號(hào)作為測試，將模型得到的匹配結(jié)果與MS/MS檢測值(真實(shí)值)進(jìn)行比對得到準(zhǔn)確度。對時(shí)間模型和峰形模型分別進(jìn)行單獨(dú)測試，然后按照峰形和時(shí)間模型不同的權(quán)重配比進(jìn)行綜合測試，得到最優(yōu)配比后再測試1次。

單獨(dú)使用時(shí)間模型進(jìn)行10次獨(dú)立測試，準(zhǔn)確率的平均測試值為96.3%，標(biāo)準(zhǔn)差為0.76。單獨(dú)使用峰形模型進(jìn)行10次獨(dú)立測試，準(zhǔn)確率的平均測試值為66.3%，標(biāo)準(zhǔn)差為1.79。

按照不同權(quán)重配比，綜合使用時(shí)間模型與峰形相似性模型，結(jié)果如表1所示?？傮w上時(shí)間模型比重越高，準(zhǔn)確性越高，這與單獨(dú)分析時(shí)間模型結(jié)果優(yōu)于峰形相似性模型的結(jié)果相符。但峰形相似性模型對整體準(zhǔn)確性有所提升，配比為8∶2時(shí)的準(zhǔn)確性達(dá)98.3%。

表1 模型在不同權(quán)重下測試結(jié)果Table 1 Testing result of different weights

采用8∶2的權(quán)重配比，綜合使用時(shí)間與峰形模型進(jìn)行10次測試，準(zhǔn)確率的平均值為97.8%，標(biāo)準(zhǔn)差為0.77。準(zhǔn)確率平均值提高了1.5%，對于并集2 847個(gè)肽鏈，提高大約為42個(gè)肽鏈的校準(zhǔn)。通過Wilcoxon rank sum test對兩種方法的結(jié)果進(jìn)行顯著性檢驗(yàn):其中方法1單純使用時(shí)間方法校準(zhǔn)，方法2采用時(shí)間峰形綜合方法校準(zhǔn)。Wilcoxon rank sum test將觀測值和零假設(shè)的中心位置之差絕對值的秩分別按照不同的符號(hào)相加作為其檢驗(yàn)統(tǒng)計(jì)量，檢驗(yàn)成對的觀測數(shù)據(jù)之差是否來自均值為0的總體(產(chǎn)生數(shù)據(jù)的總體是否具有相同的均值)。本文隨機(jī)分配訓(xùn)練集和測試集，共10次，方法1和方法2的結(jié)果均在具備相同訓(xùn)練集和測試集的情況下得到。計(jì)算觀測數(shù)據(jù)之差，共10個(gè)樣本，假設(shè)兩種方法的結(jié)果無顯著差異，即H0:兩種方法的檢測結(jié)果無顯著差異；H1:兩種方法的檢測結(jié)果有差異；在取顯著水平為0.05的條件下，通過Wilcoxon rank sum test得到P-value為0.001,h值為1，即在5%的顯著水平下拒絕H0，表明方法2相比方法1的結(jié)果有明顯改進(jìn)。

2.1.2數(shù)據(jù)并集的校準(zhǔn)匹配根據(jù)圖3，數(shù)據(jù)1與數(shù)據(jù)2并集共4 247個(gè)肽鏈，區(qū)域1中有1 944個(gè)(數(shù)據(jù)1中待校準(zhǔn)匹配的個(gè)數(shù))，區(qū)域2中有1 603個(gè)(數(shù)據(jù)2中待校準(zhǔn)匹配的個(gè)數(shù))，并集共3 547個(gè)肽鏈。通過算法匹配，最后得到匹配區(qū)間的肽鏈共3 226對，校準(zhǔn)匹配的覆蓋率達(dá)91.0%。

2.2 討 論

以上研究存在以下問題:

①區(qū)間檢測準(zhǔn)確性需提高。本文的區(qū)間檢測是以基礎(chǔ)峰值位置檢測到高強(qiáng)度峰區(qū)域內(nèi)背景噪聲標(biāo)準(zhǔn)偏差的3倍作為閾值，高于閾值的信號(hào)被納入?yún)^(qū)間，長度超過連續(xù)6個(gè)點(diǎn)的信號(hào)被認(rèn)為是候選LC峰區(qū)間。但數(shù)據(jù)1和數(shù)據(jù)2的MS/MS交集共700個(gè)肽鏈，僅檢測到599個(gè)包含MS/MS時(shí)間點(diǎn)區(qū)間的肽鏈，檢測覆蓋率為85%。

②數(shù)學(xué)模型區(qū)分度需提高。通過驗(yàn)證，時(shí)間差模型區(qū)分度好，但在噪聲較多的XICs中，仍受大量干擾信號(hào)影響。峰形相似性模型的引入雖有所改善，但區(qū)分度比時(shí)間差模型差。本文僅用線性回歸決定系數(shù)(r2)描述兩個(gè)信號(hào)相似性的值，模型的準(zhǔn)確性有待提高。

③雙模型的混合應(yīng)用。本文對兩個(gè)模型的綜合使用采取簡單的設(shè)置權(quán)值后相加的方法，今后將探索建立一個(gè)數(shù)學(xué)模型對兩個(gè)特征進(jìn)行統(tǒng)一。

3 結(jié) 論

本文通過采用統(tǒng)計(jì)學(xué)習(xí)的方法，利用多次重復(fù)的液相色譜-質(zhì)譜實(shí)驗(yàn)的譜圖中肽信號(hào)的時(shí)間差與峰形相似性兩個(gè)特征，選取訓(xùn)練數(shù)據(jù)集建立統(tǒng)計(jì)模型，并對模型有效性進(jìn)行驗(yàn)證，完成了對譜圖的校準(zhǔn)，并實(shí)現(xiàn)了多個(gè)肽信號(hào)對的匹配，準(zhǔn)確性達(dá)98%以上，覆蓋率達(dá)91.0%，為基于多次重復(fù)LC-MS實(shí)驗(yàn)數(shù)據(jù)的肽鏈量化提供了有意義的數(shù)據(jù)支撐。

致 謝:Michelle Zhang及RCMI Proteomics and Protein Biomarkers Cores at UTSA實(shí)驗(yàn)室為本研究提供了生物實(shí)驗(yàn)數(shù)據(jù)，并為論文寫作與修改提供巨大幫助。