數據不依賴獲取的質譜數據的深度學習分析方法

何情祖，鐘傳奇，李 翔，帥建偉,3*，韓家淮,3*

(1.廈門大學物理科學與技術學院，福建廈門361005；2.廈門大學生命科學學院，福建廈門361102；3.廈門大學健康醫(yī)療大數據國家研究院，福建廈門361102)

蛋白質譜儀的采集策略可以分為數據依賴獲取(data-dependent acquisition,DDA)和數據不依賴獲取(data-independent acquisition，DIA)兩種.鳥槍(shotgun)法是一種典型的DDA采集模式[1].在鳥槍法實驗中，質譜儀自動選擇強度排名前N名的肽離子進行破碎，并記錄所得的碎片離子質譜圖.由于鳥槍法只選擇強度排名靠前的肽段進行打碎，具有一定隨機性，這使得鳥槍法檢測到的肽段重復性較差.

為了解決DDA檢測到的肽段重復性差的問題，2012年瑞士Gillet等[2]與AB-SCIEX公司聯合開發(fā)出DIA相應的采集模式——全碎片離子順序窗口化獲取(sequential windowed acquisition of all theoretical fragment ions，SWATH).SWATH是一種典型的DIA模式，與傳統的質譜數據采集方法相比，SWATH模式通過高速掃描每個荷質比窗口內所有的肽離子,并將其裂解成子離子，從而獲得肽段和碎片完整的信息.相比于DDA，DIA模式具有通量高、特異性好、靈敏度高、定量結果重復性高等特點.但是DIA的質譜數據極其復雜，不同肽段的信號混合在同一張譜圖內，加大了數據分析的難度.

為了分析DIA質譜數據，相關算法的開發(fā)一直是研究熱點.2014年Rost等[3]提出一種文庫依賴的SWATH質譜數據分析方法：OpenSWATH.該方法先使用DDA質譜數據進行數據庫搜索，建立肽段的目標文庫，然后根據目標文庫對DIA質譜數據進行定量分析.2015年Tsou等[4]提出一款開源軟件DIA-Umpire，直接分析DIA模式下的質譜數據，省去了DDA實驗.在使用質譜技術研究生物學問題時，一般會制備多組樣品，并且研究者主要關注不同樣品間有差異的蛋白.同年，廈門大學韓家淮課題組[5]提出一種DIA的質譜數據分析方法.該方法通過計算不同組實驗樣品之間肽段產生的母離子和子離子信號的相關性，從而生成虛擬譜圖并通過搜索引擎來鑒定肽段.同時Wang等[6]提出另一種用于DIA質譜數據的匹配分析方法MSPLIT-DIA(mixture-spectrum partitioning using libraries of identified tandem mass spectra-DIA).

此外，近兩年也出現了一些基于深度學習的質譜數據分析算法.2019年Tran等[7]提出了DeepNovo-DIA算法，將深度學習和從頭測序(de-novo sequencing)法結合起來，對DIA質譜數據直接進行肽段氨基酸序列的鑒定.同年，Gessulat等[8]提出了名為Prosit的基于深度學習的算法，直接從蛋白質數據庫中理論預測肽段的駐留時間和子離子強度，同樣摒棄了DDA實驗.2020年Demichev等[9]提出DIA-NN算法，使用深度學習替代OpenSWATH工作流中打分步驟，以期發(fā)現更多的肽段和蛋白.同年，Yang等[10]提出了DeepDIA算法，同樣使用深度學習從蛋白質數據庫中預測肽段的駐留時間和子離子強度，其性能要強于同類的Prosit算法.

為了能夠直接處理DIA質譜數據，生成虛擬譜圖，不需要DDA數據進行建庫，且輸出文件能夠與現有的肽段定性定量工作流兼容，本研究提出了基于深度學習的DIA質譜數據分析算法：Ultra-DIA.

1 Ultra-DIA算法分析流程

1.1 深度變分自動編碼器

圖1 3種自動編碼器結構示意圖Fig.1Schematic diagram of three kinds of autoencoder

自動編碼器是深度學習中一種重要的無監(jiān)督學習模型，主要用于數據的降維和特征抽取，它的本質是利用非線性神經網絡生成的低維特征來代替高維輸入.常用的自動編碼器有深度變分自動編碼器(deep variational autoencoder，DVAE)[11-12]、堆棧自動編碼器(stacked autoencoder,SAE)[13]和降噪自動編碼器(denoising autoencoder,DAE)[14]等.DVAE、SAE和DAE的結構分別如圖1所示，其編碼器和解碼器都為多層全連接網絡.其中，SAE和DAE常用的損失函數為均方差損失函數(mean square error loss, MSELoss).相較于SAE的簡單結構，DAE通過隨機擦除輸入數據的值來學習魯棒性更好的特征，而DVAE則進一步考慮高斯隨機采樣的操作，結合貝葉斯定理和KL散度(Kullback-Leibler divergence)來約束低維特征分布，從而提升了特征的表達能力和樣本重構的表現.Ultra-DIA采用DVAE作為深度學習模型.

自動編碼器實現的映射為：

φ：X→Z,

如圖1所示，DVAE主要分成以下3大模塊：編碼器、高斯隨機采樣器和解碼器.

假設存在Z的分布，但無法確定它的具體形式，只能期望在給定X的情況下推出Z分布，即p(Z|X).根據貝葉斯定理，可得：

由于p(X)表示輸入數據X的未知的概率分布，所以p(Z|X)無法直接計算，那么可以使用一個可解的分布q(Z|X)去近似p(Z|X).而KL散度D可以用來描述兩個分布的近似程度，其形式如下：

結合貝葉斯定理和KL散度，最終可以得到DVAE的損失L：

L=-D(q(Z|X)‖p(Z))+

EZ～q(X|Z)[logp(X|Z)]，

其中E表示期望.假設q(Z|X)服從均值為μ，方差為σ2的高斯分布Ν(μ,σ2)，而p(Z)服從標準正態(tài)分布Ν(0,1)，那么L的第一項為：

由于對Z的“采樣”操作不可求導，故DVAE引入重參數化技巧(reparameterization trick)，從正態(tài)分布ε～Ν(0,1)中獲取隨機數ε，使得采樣的操作不參與梯度下降，只讓采樣的結果參與即可.重參數化操作為：

1.2 Ultra-DIA分析流程

圖2 Ultra-DIA分析流程Fig.2Analysis workflow of Ultra-DIA

如圖2所示，Ultra-DIA的分析流程分為5步.

1) 制備樣品并通過AB SCIEX公司的TOF-5600質譜儀采集SWATH數據.由于原始數據為閉源的wiff格式，無法直接讀取數據內容，所以需要使用質譜數據分析軟件ProteoWizard包含的msconvert[15]轉換工具將wiff文件轉化為開源的mzXML格式文件.

2) 對DIA質譜數據進行預處理，將噪音信號通過信噪比閾值進行過濾，再對過濾后的一級質譜(MS1)和二級質譜(MS2)峰型信號進行歸一化，消除數量級對后續(xù)分析的影響.

3) 使用DVAE提取母離子和子離子的特征，進而對不同子離子進行分類.經過深度學習的識別后，可以得到母離子和子離子的組合和對應的駐留時間，即可生成虛擬譜圖.

4) 基于虛擬譜圖，使用Comet[16-17]和X!Tandem[18]等搜索引擎對這些譜圖進行數據庫搜索，得到肽段的定性結果.

5) 最后使用OpenSWATH工作流對肽段的定性結果進行定量分析，最終得到肽段和蛋白的定性和定量結果.

2 結果與分析

2.1 數據情況

為了驗證算法的性能，本研究制備了一個小鼠細胞的SWATH數據，梯度時間為30 min.這個數據包含100個MS1窗口，MS1質荷比范圍為400～1 200，MS2的質荷比范圍為100～1 800.在掃描時間上，MS1的掃描時間為250 ms，MS2的掃描時間為33 ms，一個循環(huán)的總時間約為3.6 s.

質譜數據包含3個維度的信息：質荷比、駐留時間和強度.質譜儀會采集不同時刻的譜圖信號，將這些譜圖按照時間排列就可以得到具有不同質荷比的離子的信號隨著時間的變化曲線，稱為色譜曲線.

在SWATH實驗中，混合蛋白質樣品被酶解后進入色譜柱進行初步分離，隨后再依次進入質譜儀.因此，肽段的信號會在駐留時間方向上具有色譜峰的特征.色譜峰是判斷肽段是否存在的重要依據.肽段離子的真實色譜峰呈現出高斯曲線的特征，但其一般不對稱且在峰開始或者結束位置會出現干擾峰.一般而言，高信噪比的色譜峰表明此時記錄到了離子信號，而低信噪比的色譜峰則是背景離子產生的干擾信號.在數據預處理階段要針對色譜峰進行過濾，保留高信噪比的離子.

不同質荷比的母離子產生的色譜峰可能互不相同，如圖3(a)所示，圖中包含了5 000個母離子在1 056～2 816 s內的提取離子流(XIC)色譜曲線.不同顏色的曲線表示不同荷質比的母離子.每條曲線都可能包含多個色譜峰，每個峰都代表不同的肽段離子.在同一條色譜曲線上的不同峰表示具有相同質荷比的不同肽段.從圖中可以看出，SWATH數據的MS1信號非常復雜，直接通過母離子來推斷肽段序列非常困難，所以需要MS2的信息來共同鑒定肽段.而在圖3(b)中，這是第50個MS1窗口內由母離子打碎產生的子離子的XIC色譜曲線.可以看出，子離子色譜曲線與母離子色譜曲線具有相似的特征，并且信噪比要高于母離子色譜曲線.

圖3 典型的母離子(a)和子離子(b)的色譜圖Fig.3Typical chromatograms of precursors (a) and fragments (b)

2.2 DVAE的訓練

深度學習的訓練是使用神經網絡的重要一步，在設計好誤差函數后，就可以通過優(yōu)化器不斷減小誤差，從而達到誤差最小.本研究對比了常見的6種深度學習優(yōu)化算法，分別為AdaDelta、Adagrad、自適應矩估計(adaptive moment estimation，Adam)、Adamax、Nadam和隨機梯度下降(stochastic gradient descent，SGD)算法[19].

如圖4所示，除SGD算法外，其余優(yōu)化算法的訓練誤差曲線在訓練約25次之后開始平穩(wěn)下降，并且降幅越來越小，最后趨于穩(wěn)定.由于Adam算法收斂速度快，本研究最終選取它作為模型的優(yōu)化器.對于Adam算法，參數設置為學習率0.001，β1=0.900，β2=0.999.

圖4 訓練損失曲線Fig.4Training loss curve

2.3 識別結果

2.3.1 肽段和蛋白的定性結果分析

Ultra-DIA的輸出文件經過搜索引擎Comet和X!Tandem進行數據庫搜索，再使用PeptideProphet、iProphet、ProteinProphet這3個軟件對數據庫搜索結果進行驗證打分，最后使用Mayu.perl腳本計算肽段水平的1%錯誤發(fā)現率(FDR)的得分，用這個得分過濾肽段和蛋白即得到定性分析結果.

如圖5所示，Ultra-DIA分析SWATH數據后找到的肽段和蛋白數量均大于主流軟件DIA-Umpire.在表示肽段的韋恩圖中，交集部分的肽段覆蓋了DIA-Umpire總肽段數的78.7%；而在表示蛋白的韋恩圖中，交集部分的蛋白對DIA-Umpire總蛋白數的占比達到91.7%.這說明在定性分析階段，Ultra-DIA能夠重現出主流軟件的結果，其找到的肽段和蛋白可信度高.此外，在肽段水平上，Ultra-DIA單獨發(fā)現的數量約是DIA-Umpire單獨發(fā)現的3.6倍，在蛋白水平上達到約8倍.可見Ultra-DIA能夠發(fā)現大量DIA-Umpire忽略的肽段和蛋白.

圖5 肽段(a)和蛋白(b)的定性分析結果Fig.5Identification analysis results of peptides (a) and proteins (b)

此外，Ultra-DIA總共發(fā)現了4 119個肽段和1 572 個蛋白，DIA-Umpire則僅為2 662個肽段和1 016 個蛋白.在肽段和蛋白的總數上，Ultra-DIA均比DIA-Umpire增加了54.7%，可見本文算法能夠找到更多的肽段和蛋白.

2.3.2 肽段和蛋白的定量結果分析

定性分析后的肽段和蛋白存在著假陽性，仍然需要進一步過濾，而過濾方式則是通過檢測肽段能否定量來進行.對于SWATH數據，能夠定量的肽段和蛋白才是研究人員所關注的.

本研究將定性分析后生成的譜圖庫和SWATH數據輸入給OpenSWATH軟件.OpenSWATH通過比較每個肽段對應的一組子離子色譜峰的相似度，對肽段進行評分；緊接著再通過PyProphet軟件對OpenSWATH輸出的評分結果進行綜合，最終過濾子離子相似度低的肽段，并計算出通過篩選的蛋白數量.通過篩選得到的肽段和蛋白即是可定量的肽段和蛋白.

如圖6所示，Ultra-DIA在肽段和蛋白的定量層次上，對DIA-Umpire的覆蓋率分別為79.4%和94.0%，相較于定性結果都有所上升.這說明在定量層次上，Ultra-DIA能夠重現出DIA-Umpire的更多結果.同時，Ultra-DIA獨自發(fā)現的肽段和蛋白數量也遠大于DIA-Umpire，這些被DIA-Umpire忽略的蛋白中可能隱藏著某些新蛋白.Ultra-DIA共發(fā)現3 830個肽段和1 349 個蛋白可以定量，而DIA-Umpire只發(fā)現2 373個肽段和820個蛋白可以定量.Ultra-DIA在肽段總數上比DIA-Umpire多61.4%，在蛋白總數上多64.5%.

圖6 肽段(a)和蛋白(b)的定量分析結果Fig.6Quantification analysis results of peptides (a) and proteins (b)

對比定性分析結果，Ultra-DIA有93.0%的肽段通過了定量篩選，而DIA-Umpire只有89.1%的肽段通過了定量篩選.同時，Ultra-DIA發(fā)現的蛋白有85.5% 可以定量，而DIA-Umpire發(fā)現的蛋白只有80.7%可以定量.這說明Ultra-DIA不僅找到了更多的肽段和蛋白，且其假陽性率更低，可信度更高.

圖7展示了肽段的電荷分布，Ultra-DIA和DIA-Umpire找到的肽段的電荷具有相似的分布，都是二價最多，三價次之，四價最少.這樣的分布符合儀器的設置和工作原理，說明Ultra-DIA找到的肽段可信度高.

圖7 肽段電荷分布Fig.7Peptide charge distribution

2.3.3 蛋白定量結果的豐度比較

在蛋白組學的研究中，樣品中的低濃度蛋白一直是關注的重點.許多關鍵蛋白或者新蛋白的濃度都很低，其質譜信號往往被噪聲淹沒，強度甚至接近噪聲.因此，發(fā)現的蛋白濃度越低說明算法的性能越優(yōu)越.

OpenSWATH的定量結果與蛋白的真實濃度成正比.本研究將OpenSWATH對蛋白信號的定量結果取以10為底的對數，繪制蛋白信號強度的分布直方圖，用來比較兩種算法發(fā)現的蛋白濃度分布.

如圖8所示，DIA-Umpire發(fā)現的蛋白以較高濃度居多，而Ultra-DIA則可以發(fā)現更低濃度的蛋白.這說明在分析相同數據的情況下，Ultra-DIA發(fā)現關鍵蛋白或者新蛋白的能力要強于DIA-Umpire.

圖8 蛋白信號強度分布Fig.8Protein signal intensity distribution

2.3.4 3種自動編碼器之間的比較

本研究對比了DVAE、SAE和DAE這3種自動編碼器在測試數據上的肽段鑒定結果.對于SAE和DAE，還對比了MSELoss和二分類交叉熵損失函數(binary cross entropy loss, BCELoss)對結果的影響.

如圖9所示，DVAE鑒定到的肽段數量最多，其次是采用BCELoss的DAE，因此在Ultra-DIA中本研究采用DVAE作為深度學習模型.此外，采用BCELoss的SAE和DAE鑒定到的肽段均多于采用MSELoss的SAE和DAE，這說明采用BCELoss提取到的特征更好.

圖9 不同自動編碼器的肽段鑒定結果Fig.9Identification results of peptides by different autoencoders

3 結 論

本研究針對DIA數據，提出了基于深度學習的算法Ultra-DIA.該算法先對DIA數據進行預處理，然后使用深度DVAE分析DIA數據后生成虛擬譜圖，再對虛擬譜圖使用Comet和X!Tandem等搜索引擎進行數據庫搜索，得到肽段和蛋白的定性分析結果和譜圖庫，最后使用OpenSWATH分析譜圖庫并對肽段和蛋白進行定量.

在測試數據集上，Ultra-DIA的性能要優(yōu)于DIA-Umpire，但仍有許多可以改進的地方：1) 從質譜儀工作原理出發(fā)，對數據進行進一步的預處理；2) 神經網絡的結構仍未達到最優(yōu)，還可以設計更多不同類型的網絡，通過對比最終找到的肽段和蛋白來確定哪種網絡結構最合適；3) 神經網絡的超參數仍未找到最優(yōu)，需要不斷測試；4) 測試數據的梯度時間較短，可以考慮使用不同梯度時間的數據來綜合檢驗算法性能.