基于貝葉斯統(tǒng)計方法的多重液相色譜-質(zhì)譜試驗數(shù)據(jù)匹配研究

崔 健

(中國石油大學(xué)勝利學(xué)院 基礎(chǔ)科學(xué)學(xué)院,山東 東營 257061)

液相色譜-質(zhì)譜聯(lián)用儀(LC-MS)是由液相色譜儀與質(zhì)譜儀結(jié)合而構(gòu)成的分析儀器,它結(jié)合了液相色譜儀有效分離熱不穩(wěn)性、高沸點化合物的分離能力與質(zhì)譜儀很強的組分鑒定能力,是一種分離分析復(fù)雜有機(jī)混合物的有效手段[1],是發(fā)現(xiàn)并分析生物標(biāo)志物中復(fù)雜肽信號的關(guān)鍵技術(shù)[2]。為了得到肽鏈更準(zhǔn)確信息,部分試驗采用二級質(zhì)譜聯(lián)用(MS/MS),通過碰撞誘導(dǎo)解離給出化合物的碎片離子等結(jié)構(gòu)信息,能量越大打成的碎片越多。由低級別離子對肽鏈成分進(jìn)行進(jìn)一步分析,可以降低對質(zhì)譜的要求,能夠獲取到肽鏈組成、準(zhǔn)確的電荷數(shù)目及時間等信息。在實際操作中,為了提高肽鏈檢測覆蓋率及量化準(zhǔn)確度,經(jīng)常采取對相同樣本的多次重復(fù)試驗[3],理論上同種肽鏈在不同次試驗中應(yīng)該出現(xiàn)在相同LC時間與M/Z位置,試驗譜圖應(yīng)該是一致的[4]。但是,由于試驗誤差不可避免,譜圖普遍存在時間偏移的情況,因此需要對多個譜圖進(jìn)行校準(zhǔn)[5]。目前,比較通用的軟件如Quil[6]、proteinquant[7]、msinspect[8]、OpenMS[9]和superhirn[10]等對于重復(fù)試驗數(shù)據(jù)校準(zhǔn)基本為整體時間譜圖校準(zhǔn)。對于復(fù)雜譜圖,例如較小的時間窗口中產(chǎn)生多個LC峰的情況,這樣時間修正就會存在修正錯誤的問題。因此,本次研究采用二級質(zhì)譜聯(lián)用(MS/MS)獲取的肽鏈信息作為訓(xùn)練序列,通過Warping函數(shù)來進(jìn)行時間校準(zhǔn),并聯(lián)合使用貝葉斯統(tǒng)計方法對Warping函數(shù)進(jìn)行提升,對任意峰對給出相關(guān)信號概率及非相關(guān)信號概率,并驗證有效性。最后,將多個數(shù)據(jù)的肽鏈信號對通過該方法進(jìn)行校準(zhǔn)匹配,并驗證覆蓋率。

1 數(shù)據(jù)分析

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

處理數(shù)據(jù)由RCMI Proteomics and Protein Biomarkers Cores試驗室產(chǎn)生,經(jīng)過LTQ OrbitrapVelos儀器處理的一組TAGE腫瘤樣本。LC-MS試驗是將蛋白質(zhì)切割成肽鏈,并使用試劑利用肽鏈斥水性不同的特性,將其沖入到質(zhì)譜儀中。斥水性不同導(dǎo)致肽鏈進(jìn)入質(zhì)譜儀形成譜圖的時間就不一致,形成了不同肽鏈時間上的區(qū)分。進(jìn)入質(zhì)譜儀中的肽鏈將隨機(jī)帶上電荷,根據(jù)不同肽鏈大小、質(zhì)量、帶電荷不同的特性,形成的質(zhì)量與電荷比(M/Z)值不同,形成質(zhì)荷比維度的區(qū)分。同一種肽鏈由于具備相同斥水性及質(zhì)荷比,因此理論上將出現(xiàn)在譜圖中的一個特定位置上,此類譜圖為Level 1數(shù)據(jù)(圖1)。由于儀器操作等影響,為了提高精度,一般將進(jìn)一步進(jìn)行MS/MS試驗,即從Level1譜圖中隨機(jī)選取位置,將肽鏈進(jìn)行成分分析確定肽鏈組成,稱為Level 2數(shù)據(jù)。

圖1 數(shù)據(jù)1譜圖

本次研究從多組試驗中選取了兩組數(shù)據(jù)(數(shù)據(jù)1與數(shù)據(jù)2)進(jìn)行分析,由MS/MS檢測到肽信號信息(圖2)。

圖2 MS/MS檢測結(jié)果韋恩圖

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

數(shù)據(jù)處理分為“數(shù)據(jù)預(yù)處理”、“訓(xùn)練與測試數(shù)據(jù)集生成”、“Warping函數(shù)及貝葉斯統(tǒng)計方法建?！薄ⅰ澳Ｐ万炞C及全集校準(zhǔn)”四個部分。

1.2.1 數(shù)據(jù)預(yù)處理

在數(shù)據(jù)預(yù)處理中,根據(jù)MS/MS信息表,生成肽信號合集,并計算肽鏈荷質(zhì)比(M/Z值),以肽鏈M/Z值為中心,前后20×10-6寬度,計算LC譜圖,生成肽鏈的全時間段XICs(圖3)。然后,在全時段XICs上進(jìn)行區(qū)間檢測。

圖3 肽鏈“CSTSSLLEACTFR”全時段XICs

1.2.2 訓(xùn)練與測試數(shù)據(jù)集生成

生成訓(xùn)練和測試數(shù)據(jù)集的原則是肽鏈信號必須具備可驗證的真實值(ground truth)。由MS/MS檢測到的肽信號即為真實信號,其具有的M/Z值與時間值即為真實可靠的值。因此,選擇圖2交集部分(共700個肽鏈)作為訓(xùn)練與測試數(shù)據(jù)集。在訓(xùn)練測試數(shù)據(jù)集中,首先進(jìn)行區(qū)間檢測預(yù)處理,能夠檢測到的區(qū)間,即具備較好的峰值。區(qū)間包含MS/MS時間點即為可用肽鏈。經(jīng)過區(qū)間檢測預(yù)處理,共599個肽鏈可用。采用2-折交叉驗證,即隨機(jī)選取一半作為訓(xùn)練,一半作為測試。訓(xùn)練序列用于warping函數(shù)與貝葉斯統(tǒng)計模型的生成,測試序列用于測試模型產(chǎn)生的肽信號匹配結(jié)果的準(zhǔn)確性(以MS/MS檢測值作為真實數(shù)據(jù))。

1.2.3 Warping函數(shù)及貝葉斯統(tǒng)計方法建模

對于生成的訓(xùn)練序列,選取由MS/MS確定的相關(guān)信號時間對生成Warping函數(shù)。在使用Warping函數(shù)擬合之前剔除偏差較大的奇異點。采取每個訓(xùn)練序列肽信號的數(shù)據(jù)1時間減去數(shù)據(jù)2時間,然后計算方差,采用平均值±3倍標(biāo)準(zhǔn)差作為標(biāo)準(zhǔn),將外部的點去掉,直方圖如圖4所示。

圖4 時間差直方圖

然后采用四階多項式作為Warping函數(shù)擬合,生成多項式參數(shù)，擬合后如圖5所示。

圖5 Warping函數(shù)擬合結(jié)果

下一步將計算由MS/MS檢測結(jié)果確定為相關(guān)信號的時間對、確定為非相關(guān)信號的時間對,分別距離Warping函數(shù)的時間差值,如圖6所示。

分別對相關(guān)信號時間差值和非相關(guān)信號時間差值建模。以相關(guān)信號時間差值建模為例,根據(jù)測試數(shù)據(jù)中的時間差,得到n個樣本t(t1,t2,t3,…,tn),其中ti是測試序列中第i對信號的時間差值。通過對樣本的直方圖觀察,基本符合正態(tài)分布特征,一般情況下正態(tài)分布的概率密度函數(shù)中包含的兩個參數(shù)μ和σ由樣本值進(jìn)行最大似然估計。使用貝葉斯統(tǒng)計方法進(jìn)行建模,首先設(shè)定先驗信息(Prior),即將參數(shù)μ和σ看作為兩個隨機(jī)變量,其服從以下分布特征:

p(μ|σ2)～N(μ0,σ2/κ0),

圖6 相關(guān)信號與非相關(guān)信號時間差直方圖

根據(jù)貝葉斯公式,參數(shù)μ和σ的聯(lián)合分布為

p(μ,σ2)=p(μ|σ2)p(σ2),

即為

簡化一下即可得:

下一步進(jìn)行后驗信息(posterior)計算,

p(μ,σ2|t)=p(t|μ,σ2)p(μ,σ2),其中p(μ,σ2)為先驗信息已經(jīng)計算獲得。

而p(t|μ,σ2)～N(μ,σ2)是μ和σ的正態(tài)分布。這樣計算p(t|μ,σ2)得:

p(μ,σ2|t)∝σ-1(σ2)-(1+(ν0+n)/2)×

令

得

現(xiàn)在已知先驗信息,后驗信息如下:

(1)

p(μ|σ2,t)～N(μn,σ2/κn),

(2)

(3)

p(x|μ,σ2,t)～N(μ,σ2).

(4)

式中,x為任意時間差,為隨機(jī)變量;t為已經(jīng)獲得的樣本值,那么下一步將計算p(x|t)的值,這樣就無須估算正態(tài)分布中μ和σ的值,直接由樣本值t計算任意時間差變量x的分布。

將(1)、(2)、(4)代入以上公式得:

(5)

(6)

1.2.4 模型驗證及全集校準(zhǔn)

模型的建立及測試是在訓(xùn)練與測試數(shù)據(jù)集上進(jìn)行的,即圖2的區(qū)域C部分。隨機(jī)選取一半作為訓(xùn)練序列建立以上模型,另一半數(shù)據(jù)驗證模型的有效性,并與MS/MS檢測的真實數(shù)據(jù)比對獲取模型準(zhǔn)確度。在驗證模型有效后,對圖2中區(qū)域A和區(qū)域B中的肽鏈信號通過模型進(jìn)行匹配。區(qū)域A中信號為由MS/MS檢測到的僅在數(shù)據(jù)1中有真實值的肽信號,通過模型匹配得到其在數(shù)據(jù)2中的匹配信號;同理,區(qū)域B中亦是如此。這樣就完成全集的校準(zhǔn)匹配,得到整體的匹配覆蓋率。

2 結(jié)果分析與問題討論

得到的結(jié)果主要有兩部分,一是通過貝葉斯統(tǒng)計方法改進(jìn)Warping函數(shù)校準(zhǔn)匹配有效性結(jié)果;二是全集最終校準(zhǔn)匹配結(jié)果。

2.1 模型有效性結(jié)果

本次研究進(jìn)行了10次測試,每次從訓(xùn)練與測試序列中隨機(jī)選取300個進(jìn)行Warping函數(shù)擬合,然后計算時間差,并用貝葉斯統(tǒng)計方法訓(xùn)練建立模型。另外,299個作為模型測試,一是單獨使用Warping函數(shù),判斷測試序列中時間距離Warping曲線最近的區(qū)間為匹配校準(zhǔn)區(qū)間;二是使用改進(jìn)的Warping函數(shù)與貝葉斯統(tǒng)計方法建立的模型,如果相關(guān)信號模型給出的概率大于非相關(guān)信號的模型概率,則判斷為匹配。以上兩種結(jié)果均與MS/MS時間點真實值進(jìn)行比對,計算準(zhǔn)確度如表1所示。

表1 測試結(jié)果對比

由表1看出,用Warping的測試結(jié)果準(zhǔn)確性均值為86.81%,通過Warping函數(shù)聯(lián)合使用貝葉斯統(tǒng)計方法建模準(zhǔn)確率均值達(dá)到93.08%,提高了6.27個百分點。

2.2 數(shù)據(jù)全集的校準(zhǔn)匹配

由MS/MS檢測到的數(shù)據(jù)1與數(shù)據(jù)2的肽鏈共4 247個,分布如圖2所示。交集共700個,通過區(qū)間檢測的信號共599個,在此基礎(chǔ)上采取Warping函數(shù)聯(lián)合貝葉斯統(tǒng)計方法建模進(jìn)行匹配。區(qū)域A中1 944個,區(qū)域B中1 603個,共3 547個肽鏈。根據(jù)模型給出的匹配與非匹配的概率大小,共能實現(xiàn)3 185個肽鏈在另外數(shù)據(jù)中的區(qū)間匹配,覆蓋率達(dá)到89.8%。

2.3 問題討論

基于以上研究結(jié)果,可以看出,交集的700個肽鏈中只有599個能被檢測到信號區(qū)間,檢測到區(qū)間的概率大約為85%。這說明部分被MS/MS檢測到的肽鏈信號非常弱,無法在level1數(shù)據(jù)中被檢測出來。這是由于區(qū)間檢測不準(zhǔn)確造成的,本文在區(qū)間檢測中采用的是簡單的低于最高峰值20%即在區(qū)間外的辦法,很多情況下并不有效,檢測不到真實的肽信號區(qū)間,因此,下一步將重點研究準(zhǔn)確的區(qū)間檢測算法。

3 結(jié)束語

采用Warping函數(shù)聯(lián)合貝葉斯統(tǒng)計方法建模對多次重復(fù)的液相色譜-質(zhì)譜數(shù)據(jù)進(jìn)行時間校準(zhǔn),根據(jù)MS/MS檢測值選取訓(xùn)練序列進(jìn)行時間擬合。通過測試序列驗證,單獨使用Warping的測試結(jié)果準(zhǔn)確性均值為86.81%;通過Warping函數(shù)聯(lián)合使用貝葉斯統(tǒng)計方法建模準(zhǔn)確率均值達(dá)到93.08%。同時,完成兩個譜圖的匹配校準(zhǔn),覆蓋率超過89%。對下一步進(jìn)行肽鏈量化提供了非常有意義的算法支撐。