馬媛媛

摘要：文章是對(duì)多次重復(fù)液相質(zhì)譜（LC-MS）實(shí)驗(yàn)得到的蛋白質(zhì)肽鏈生物數(shù)據(jù)進(jìn)行時(shí)間校準(zhǔn)建模分析，目的是校準(zhǔn)匹配多次實(shí)驗(yàn)中的肽鏈信號(hào)，為蛋白質(zhì)量化提供準(zhǔn)確信息。本項(xiàng)目研究一種基于時(shí)間特征的LC-MS生物實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)校準(zhǔn)算法，解決了現(xiàn)在生物數(shù)據(jù)處理中的實(shí)際問題，具有現(xiàn)實(shí)意義。

Abstract： This paper is a modeling analysis of large biological data of protein peptide chains obtained by liquid chromatography-mass spectrometry （LC-MS） experiments. The purpose of this study was to align the signals of the same peptide in different datasets， in order to provide the accurate quantification information. This project， which provides a method for the alignment based on the time feature to solve the actual problem of biological data processing， has practical significance.

關(guān)鍵詞：蛋白質(zhì)肽鏈；時(shí)間特征；統(tǒng)計(jì)建模

Key words： protein peptide chain；time feature；statistical modeling

中圖分類號(hào)：C37 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1006-4311（2018）21-0194-03

0 引言

隨著生物實(shí)驗(yàn)技術(shù)的高速發(fā)展，生命科學(xué)研究獲得大量生物實(shí)驗(yàn)數(shù)據(jù)，主要包括基因組學(xué)、蛋白質(zhì)組學(xué)等生物學(xué)大數(shù)據(jù)，這些數(shù)據(jù)均具有4V的特性：①數(shù)據(jù)量大（Volume）：目前基因組學(xué)中只需花費(fèi)幾千美元幾個(gè)小時(shí)即可完成一個(gè)人基因組的解析，大量的物種得以測序解析，數(shù)據(jù)成爆炸性增長。②數(shù)據(jù)多樣化（Variety）：生物信息學(xué)中兩大分類：基因組學(xué)和蛋白質(zhì)組學(xué)中，實(shí)驗(yàn)儀器種類繁多，產(chǎn)生的數(shù)據(jù)格式也各不相同。同時(shí)，利用不同的生物信息分析軟件或分析流程處理得到的結(jié)果也是千差萬別。③有價(jià)值（Value）：隨著生物信息學(xué)的發(fā)展，越來越多有價(jià)值的信息從生物大數(shù)據(jù)中挖掘出來，這些價(jià)值不僅體現(xiàn)在其在生物科研領(lǐng)域，而且已應(yīng)用于健康和醫(yī)學(xué)等領(lǐng)域。④高速（Velocity）：不僅體現(xiàn)在數(shù)據(jù)采集量急劇增長，也表現(xiàn)在數(shù)據(jù)的多樣化和價(jià)值性上。

本項(xiàng)目處理的數(shù)據(jù)是由液相質(zhì)譜（LC-MS）實(shí)驗(yàn)得到的蛋白質(zhì)肽鏈生物大數(shù)據(jù)。實(shí)驗(yàn)的目的是為了確定蛋白質(zhì)中生物標(biāo)志物（biomarker）。這些標(biāo)志物是可以標(biāo)記系統(tǒng)、器官、組織、細(xì)胞及亞細(xì)胞結(jié)構(gòu)或功能的改變或可能發(fā)生的改變的生化指標(biāo)，可用于疾病診斷、判斷疾病分期或者用來評(píng)價(jià)新藥或新療法在目標(biāo)人群中的安全性及有效性。目前尋找和發(fā)現(xiàn)有價(jià)值的生物標(biāo)志物已經(jīng)成為科學(xué)研究的一個(gè)重要熱點(diǎn)。本項(xiàng)目重點(diǎn)研究多次重復(fù)LC-MS實(shí)驗(yàn)中的肽鏈信號(hào)匹配算法，目的是為尋找生物標(biāo)志物提供準(zhǔn)確的量化信息[1][2]。

1 LC-MS實(shí)驗(yàn)簡介

在LC-MS實(shí)驗(yàn)中，首先，將蛋白質(zhì)切割成肽鏈，并放入容器中；其次，進(jìn)入實(shí)驗(yàn)的LC部分，用化學(xué)試劑將容器中肽鏈沖入到質(zhì)譜儀，由于不同肽鏈具有不同的斥水性，因此進(jìn)入質(zhì)譜儀的時(shí)間便有所不同，形成三維譜圖中的時(shí)間軸（Time）；再次，進(jìn)入實(shí)驗(yàn)的MS部分，肽鏈進(jìn)入質(zhì)譜儀后將隨機(jī)粘上電荷，根據(jù)質(zhì)量和電荷比的不同打到檢測板上的位置就不同，這樣就形成三維譜圖中的質(zhì)荷比軸（M/Z）；最后，相同時(shí)間，相同質(zhì)荷比的位置上基本由同一種肽鏈組成，個(gè)數(shù)越多強(qiáng)度越大形成三維譜圖中的強(qiáng)度軸（Intensity）。經(jīng)過一次MS處理的數(shù)據(jù)稱為Level1數(shù)據(jù)，數(shù)據(jù)粗糙，但是全面。生物實(shí)驗(yàn)中經(jīng)過在LC-MS實(shí)驗(yàn)后再進(jìn)行一次MS實(shí)驗(yàn)，得到Level2數(shù)據(jù)。Level2數(shù)據(jù)是從Level1數(shù)據(jù)中隨機(jī)抽取生成的，點(diǎn)數(shù)多但是覆蓋率不足，量化準(zhǔn)確性不高，基本用來和肽鏈庫進(jìn)行比對(duì)，確定Level1數(shù)據(jù)中的部分肽鏈組成，進(jìn)而確定Level1數(shù)據(jù)中的蛋白質(zhì)組成。

2 欲解決的問題

2.1 問題的提出

目前，在LC-MS生物大數(shù)據(jù)處理中的重要任務(wù)就是對(duì)各種檢測到的肽鏈進(jìn)行量化分析，面臨的一個(gè)重要問題就是對(duì)相同樣本的多次重復(fù)實(shí)驗(yàn)數(shù)據(jù)中肽鏈產(chǎn)生的信號(hào)進(jìn)行校準(zhǔn)識(shí)別，這項(xiàng)工作對(duì)于減少多次重復(fù)實(shí)驗(yàn)產(chǎn)生量化誤差，提高量化準(zhǔn)確性是至關(guān)重要的。但是進(jìn)行相同樣本的多次重復(fù)實(shí)驗(yàn)生物大數(shù)據(jù)的在特征檢測、校準(zhǔn)識(shí)別、量化分析時(shí)，理論上在重復(fù)實(shí)驗(yàn)條件完全一致的情況下，同一種肽鏈在不同重復(fù)實(shí)驗(yàn)數(shù)據(jù)中的相對(duì)應(yīng)位置（相同時(shí)間值，相同M/Z值）應(yīng)該產(chǎn)生相同的特征峰值。實(shí)際中，由于各種誤差因素的存在，重復(fù)實(shí)驗(yàn)數(shù)據(jù)的時(shí)間軸也會(huì)產(chǎn)生較大差異，這樣就無法對(duì)同一肽鏈在多組數(shù)據(jù)中進(jìn)行相關(guān)峰值識(shí)別校準(zhǔn)，進(jìn)一步說就無法量化分析。這就需要我們對(duì)多次重復(fù)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行校準(zhǔn)。

2.2 方法思路

我們要處理重復(fù)實(shí)驗(yàn)數(shù)據(jù)1和數(shù)據(jù)2，通過和ms2實(shí)驗(yàn)產(chǎn)生的Level2數(shù)據(jù)比對(duì)，如圖1所示，數(shù)據(jù)1與數(shù)據(jù)2在ms2中共同檢測的肽鏈共有700個(gè)（即為交集部分，區(qū)域B），這些肽鏈都能分別在數(shù)據(jù)1和數(shù)據(jù)2中找到相應(yīng)的信號(hào)區(qū)間。但是，通過ms2實(shí)驗(yàn)的檢測，數(shù)據(jù)1中依然有1944個(gè)肽鏈沒有數(shù)據(jù)2中找到（區(qū)域A部分），數(shù)據(jù)2中依然有1603個(gè)肽鏈沒有在數(shù)據(jù)1中找到（區(qū)域C部分）。那么我們將通過數(shù)據(jù)1與數(shù)據(jù)2中共同檢測到的部分建立數(shù)學(xué)模型，然后對(duì)于區(qū)域A中的在數(shù)據(jù)2中找到匹配的肽信號(hào)區(qū)間，對(duì)于區(qū)域B中的在數(shù)據(jù)1中找到匹配的肽信號(hào)區(qū)間。

2.3 數(shù)據(jù)處理流程及算法思想（圖2）

2.3.1 步驟1：數(shù)據(jù)的預(yù)處理

①由實(shí)驗(yàn)室獲取的生物大數(shù)據(jù)為mzxml格式的大數(shù)據(jù)文件，我們使用matlab中的mzxmlread函數(shù)讀取初始的mzxml文件，將mzxml實(shí)驗(yàn)數(shù)據(jù)讀出生成原始的level1數(shù)據(jù)、level2數(shù)據(jù)、原始峰值信息、level1的時(shí)間信息，并保存，同時(shí)生成實(shí)驗(yàn)數(shù)據(jù)三維譜圖如圖3。

②讀入ms2實(shí)驗(yàn)數(shù)據(jù)生成的肽鏈信息總表和數(shù)據(jù)的level1信息，根據(jù)ms2數(shù)據(jù)總表中的肽鏈的質(zhì)量值（mass）和電荷（charge state）計(jì)算出質(zhì)荷比（m/z值），按照肽鏈m/z值前后選取20ppm寬度計(jì)算LC譜圖，以獲取該肽鏈在兩組數(shù)據(jù)中的XICs（用來確定肽鏈可能產(chǎn)生的LC峰）[3]。然后對(duì)XICs做區(qū)間檢測，我們使用肽鏈主峰位置檢測到高強(qiáng)度峰區(qū)域在背景噪聲的標(biāo)準(zhǔn)偏差的三倍作為閾值，高于閾值的區(qū)間被認(rèn)為是候選LC峰區(qū)間。

③對(duì)于給定的肽，如果分別在數(shù)據(jù)1和數(shù)據(jù)2對(duì)應(yīng)的XICs中，檢測出n、m個(gè)候選LC峰區(qū)間，則會(huì)有n*m候選LC峰區(qū)間對(duì)。然而，只有一對(duì)是給定的肽在兩個(gè)重復(fù)實(shí)驗(yàn)數(shù)據(jù)中對(duì)應(yīng)產(chǎn)生的[4]。我們再處理XICs選取候選區(qū)間過程中，首先區(qū)間中形成的波峰的點(diǎn)的個(gè)數(shù)要多于6個(gè)，然后按照每一個(gè)區(qū)間的最高信號(hào)值由高到低選出前10個(gè)區(qū)間，并保存區(qū)間時(shí)間的起始和結(jié)束位置。

2.3.2 步驟2：生成訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)

對(duì)于訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)我們應(yīng)該選取有ms2時(shí)間點(diǎn)并檢測到包含時(shí)間點(diǎn)區(qū)間的肽鏈。這樣，我們再測試模型的時(shí)候才有真實(shí)值做比對(duì)，才能檢測模型的準(zhǔn)確性。那么對(duì)數(shù)據(jù)1和數(shù)據(jù)2中在ms2檢測后重復(fù)的部分，即同時(shí)在數(shù)據(jù)1和數(shù)據(jù)2中檢測到ms2時(shí)間點(diǎn)的肽鏈，共700個(gè)。以這700個(gè)肽鏈為基礎(chǔ)，我們首先對(duì)區(qū)間檢測的結(jié)果和ms2時(shí)間進(jìn)行比對(duì)，選出區(qū)間包含ms2時(shí)間的肽鏈，經(jīng)過篩選有599個(gè)肽鏈符合條件。這樣我們隨機(jī)選取400個(gè)作為訓(xùn)練序列訓(xùn)練統(tǒng)計(jì)模型，剩下的199個(gè)作為測試序列測試模型準(zhǔn)確性，重復(fù)5次，準(zhǔn)確性取平均值。

2.3.3 步驟3：建立訓(xùn)練模型

我們了解到產(chǎn)生時(shí)間偏移是隨機(jī)的，且有直方圖可以觀測到基本符合正態(tài)分布。因此，如果出現(xiàn)未確定的兩個(gè)區(qū)間的時(shí)間差Δt，我們需要根據(jù)已經(jīng)得到的相關(guān)時(shí)間差樣本t1計(jì)算概率p（Δt|t1），我們需要根據(jù)已經(jīng)得到的非相關(guān)時(shí)間差樣本t2計(jì)算概率p（Δt|t2）。我們可以根據(jù)相關(guān)概率p（Δt|t1）/p（Δt|t2）比值是否大于1來判斷是否為相關(guān)區(qū)間。我們知道正態(tài)分布的公式為：

2.3.4 步驟4：訓(xùn)練模型并測試，得出模型匹配成功率

根據(jù)步驟3中用400個(gè)訓(xùn)練序列訓(xùn)練出來的匹配模型，我們將199個(gè)測試序列輸入模型輸出匹配結(jié)果。同時(shí)由于測試序列具有ms2的檢測結(jié)果，因此，我們將測試結(jié)果與ms2結(jié)果做比對(duì)，所謂匹配成功即為模型輸出數(shù)據(jù)1與數(shù)據(jù)2中的匹配區(qū)間對(duì)均能覆蓋數(shù)據(jù)1中該肽鏈的ms2時(shí)間和數(shù)據(jù)2中該肽鏈的ms2時(shí)間。通過對(duì)學(xué)習(xí)模型進(jìn)行測試，得到模型檢測的準(zhǔn)確率。重復(fù)進(jìn)行5次隨機(jī)選取訓(xùn)練和測試，平均準(zhǔn)確率結(jié)果作為算法的準(zhǔn)確率。

2.3.5 步驟5：對(duì)所有數(shù)據(jù)應(yīng)用模型進(jìn)行匹配校準(zhǔn)

在對(duì)模型進(jìn)行測試之后，我們將模型應(yīng)用于圖1中數(shù)據(jù)1的A區(qū)域和數(shù)據(jù)2的C區(qū)域，應(yīng)用過程是：以數(shù)據(jù)1的A區(qū)域中某一肽鏈為例，由于該肽鏈在數(shù)據(jù)1中被ms2檢測到，有m/z和時(shí)間信息等，但是并沒有在數(shù)據(jù)2中被ms2檢測到。因此，我們先分別在數(shù)據(jù)1與數(shù)據(jù)2中處理生成對(duì)應(yīng)的XICs，然后進(jìn)行區(qū)間檢測，那么在數(shù)據(jù)1中我們根據(jù)ms2檢測到的時(shí)間信息確定出準(zhǔn)確區(qū)間，同時(shí)數(shù)據(jù)2中該肽鏈檢測到的區(qū)間均為候選區(qū)間。我們將數(shù)據(jù)2中的候選區(qū)間與數(shù)據(jù)1中的準(zhǔn)確區(qū)間的時(shí)間差作為模型輸入，判斷相關(guān)性的概率值和非相關(guān)性的概率值的比值，如果大于1那么我們認(rèn)為找到了該肽鏈在數(shù)據(jù)2中的區(qū)間。如果多個(gè)區(qū)間概率比值大于1，那么相關(guān)性概率值最大的就是我們要選擇的區(qū)間。這樣我們將可以將數(shù)據(jù)1區(qū)域A和數(shù)據(jù)2區(qū)域C中的肽鏈分別在數(shù)據(jù)2中和數(shù)據(jù)1中找到了匹配區(qū)間。

3 數(shù)據(jù)處理結(jié)果

我們對(duì)算法的驗(yàn)證是通過對(duì)交集肽信號(hào)進(jìn)行在有真實(shí)值條件下測試準(zhǔn)確率和對(duì)待校準(zhǔn)集合進(jìn)行無真實(shí)值情況下完成匹配。

3.1 基于ms2檢測結(jié)果的模型測試結(jié)果

我們進(jìn)行了5次測試，每次都是隨機(jī)選取400作為訓(xùn)練、199作為測試，對(duì)測試結(jié)果以MS2時(shí)間點(diǎn)為真實(shí)值進(jìn)行比對(duì)，得到區(qū)間匹配準(zhǔn)確度結(jié)果如表1。

3.2 無ms2檢測結(jié)果的待校準(zhǔn)集合匹配結(jié)果

數(shù)據(jù)1與數(shù)據(jù)2并集共4247個(gè)肽鏈，數(shù)據(jù)1中待校準(zhǔn)匹配的個(gè)數(shù)為1944個(gè)，數(shù)據(jù)2中待校準(zhǔn)匹配的個(gè)數(shù)為1603個(gè)。經(jīng)過我們對(duì)數(shù)據(jù)1和數(shù)據(jù)2非交集中共3547個(gè)肽鏈信號(hào)進(jìn)行算法匹配，最后得到區(qū)間結(jié)果的共3098對(duì)，校準(zhǔn)匹配的覆蓋率達(dá)到87.34%。

這樣實(shí)現(xiàn)了對(duì)兩組數(shù)據(jù)大部分的肽鏈的匹配校準(zhǔn)工作，且匹配成功的概率為96.32%，而且達(dá)到了比較高的匹配覆蓋率。

參考文獻(xiàn)：

[1]寧康，陳挺.生物醫(yī)學(xué)大數(shù)據(jù)的現(xiàn)狀與展望[J].北京：科學(xué)通報(bào)，2015，Z1.

[2]胡瑞峰，邢小燕，孫桂波，孫曉波.大數(shù)據(jù)時(shí)代下生物信息技術(shù)在生物醫(yī)藥領(lǐng)域的應(yīng)用前景[J].北京：藥學(xué)學(xué)報(bào)，2014，11.

[3]Smith R， Ventura D， Prince J T. LC-MS alignment in theory and practice： a comprehensive algorithmic review.[J]. Briefings in Bioinformatics， 2015， 16（1）：104.

[4]Bielow C， Mastrobuoni G， Kempa S. Proteomics Quality Control： Quality Control Software for MaxQuant Results[J]. Journal of Proteome Research， 2015， 15（3）.