基于圖卷積神經網絡的串聯質譜從頭測序

牟長寧，王海鵬，周丕宇，侯鑫行

（山東理工大學計算機科學與技術學院，山東淄博 255000）

（*通信作者電子郵箱hpwang@sdut.edu.cn）

0 引言

基于串聯質譜的蛋白質組學是生命科學研究的重要組成部分，近年來在探索細胞機制、疾病進程以及基因型和表型之間關系等研究上取得了巨大進展［1］?；诖撡|譜的蛋白質測序主流的方法是蛋白質數據庫搜索，常用工具有Mascot［2］、Comet［3］、MaxQuant［4］、pFind［5-6］等。該方法需要參考已有數據庫檢索候選肽序列，因此在未知生物蛋白、單克隆抗體測序等研究上失去優勢。另一種鑒定方法是從頭測序，該方法能夠直接從串聯質譜中推斷出氨基酸序列，無需數據庫作為參考，在鑒定未知生物肽序列上具有不可替代的作用。過去20 年間從頭測序方法進步顯著，應用較為廣泛的方案是基于圖論的思想，將質譜轉化為譜峰關系圖（spectrum graph），譜峰作為譜圖中的頂點，如果譜峰與譜峰之間的距離等于一個或者兩個氨基酸殘基的分子量，則兩個峰之間用一條邊相連；通過搜索圖中起始點到結束點的最優路徑得到產生這個質譜的候選肽序列。代表性工作包括：2003 年Ma 等［7］發表的PEAKS，通過預處理步驟（圖譜噪聲過濾和圖譜峰聚合）創建譜圖并用動態規劃算法來生成候選肽序列；2005年Frank等［8］發表了針對碰撞誘導裂解（Collision-Induced Dissociation，CID）質譜的PepNovo 算法，提出了一種基于概率網絡模型的候選肽序列評分方法；2010 年Chi 等［9］發表的pNovo，使用帶剪枝的深度優先搜索有效提升了在高能碰撞裂解（Higher-energy Collision Dissociation，HCD）質譜數據上的從頭測序性能；隨后同一團隊，在pNovo 基礎上開發了同時使用HCD 和電子轉運裂解（Electron Transfer Dissociation，ETD）數據的從頭測序方法pNovo+［10］，以及針對翻譯后修飾肽鑒定的OpenpNovo［11］，并在2019年發表了pNovo3［12］，將理論質譜預測用于候選肽重排序。另一類從頭測序方法則是基于機器學習和深度學習技術。2005年，NovoHMM算法［13］提出使用隱馬爾可夫模型解決從頭測序問題；2015 年，Novor［14］使用決策樹模型分別為碎片離子和氨基酸殘基進行打分，結合動態規劃推導肽序列；基于深度學習的DeepNovo［15］，通過基于卷積神經網絡（Convolutional Neural Network，CNN）的Ion-CNN 和Spectrum-CNN，以及長短期記憶（Long Short-Term Memory，LSTM）網絡模型融合的方式對肽序列進行預測。隨著從頭測序方法的改進，測序精度不斷得到提升，然而由于質譜儀中肽不完全碎裂等因素，導致質譜中碎片離子的覆蓋率較低，重要b離子或y離子峰丟失，大量噪聲干擾峰難以通過約束條件徹底清除，諸多因素致使從頭測序的精度仍然較低，嚴重制約了從頭測序在蛋白質組數據分析中的應用。因此提升肽段從頭測序準確性，對蛋白質組學研究具有重要意義。

在蛋白質組學中，深度學習方法已經應用到了預測肽段保留時間、理論質譜預測、翻譯后修飾、從頭測序、蛋白質結構預測等多個任務中［16-17］。深度學習的蓬勃發展，為質譜數據分析不斷提供新的方案啟示。本文在經典的譜峰關系圖方法基礎上，提出了一種基于圖卷積神經網絡（Graph Convolutional neural Network，GCN）的從頭測序方法denovo-GCN。該方法直接使用質譜數據作為輸入，簡化中間數據約束處理過程，在譜峰關系圖結構上按照碎裂位點為每個譜峰構造特征表示。通過在大規模數據上的訓練優化，能夠有效提升從頭測序的準確性。

1 denovo-GCN模型

1.1 圖卷積神經網絡與從頭測序

圖卷積神經網絡以其在圖數據上的強大建模能力，在知識圖譜、社交網絡等眾多領域得到了應用［18］。Kipf 等［19］對ChebNet［20］進行了簡化，提出了一種更加簡單的模型GCN，它相當于對一階切比雪夫卷積的再近似，降低了計算復雜度，并且可以通過堆疊多個GCN 擴大圖卷積神經網絡的感受野，實用性大大增強。GCN模型結構表述為式（1）：

其中：=A+I，A是圖的鄰接矩陣，包含了節點之間的連接信息；I是單位矩陣，加上I后得到的包含了自身節點和鄰接點的信息是頂點的度矩陣是激活函數；H(l) ∈Rn×m是第l層的激活矩陣；H(0)=X，X是由各節點特征向量xi組成的特征矩陣。隨后注意機制、序列模型等也用于圖中節點權重的計算，圖卷積神經網絡呈現出多樣化的發展。

從頭測序過程可以類比為語言翻譯或者圖像描述，最終目的是得到一個映射原始數據的序列表示。不同之處在于后者的原始數據是規則歐氏空間數據，而質譜數據是一組譜峰質荷比及其強度的數據對組成的集合。在基于圖論的從頭測序中，譜峰關系是由譜峰之間的距離來計算，形成譜峰連接圖。這種質譜數據圖結構化的表示方法與針對圖結構數據的圖卷積神經網絡十分契合。譜峰節點的特征則可以通過枚舉碎裂位點產生的離子與各譜峰的距離關系表示，借助圖卷積神經網絡的訓練學習能力將符合條件的譜峰點與干擾峰進行區分，預測當前位置的氨基酸身份，逐步實現氨基酸序列的推理。

1.2 圖的構建

在質譜數據上，使用圖卷積神經網絡的首要任務是構建譜峰連接圖。質譜數據中的關鍵信息包括母離子的質荷比、肽所帶電荷、譜峰。譜峰是碎片離子質荷比及其強度組成的數據對，將譜峰強度值按照同一質譜中最大強度值歸一化得到相對強度，相對強度最大值為1。單個質譜可以直觀表示為質荷比和相對強度的柱狀圖，x軸代表質荷比，y軸代表強度。若譜峰與譜峰之間的距離與一個氨基酸殘基的分子量的差值在設定誤差范圍內，則兩個譜峰之間建立一條邊。在構建譜圖前，需向原始譜圖添加序列端點的譜峰，分別為一個電荷（M（proton））、一個水分子量（M（H2O））、1 電荷肽的分子量（M（peptide））、肽失去一個水的分子量（M（peptide）-M（H2O））四個譜峰點，相對強度皆設置為1。設S=為譜（npeaks為譜峰數量），SA為峰與峰之間的差值矩陣，MASS_AA=（n=23，代表20 氨基酸殘基和3 種修飾后的氨基酸殘基）為氨基酸殘基質量集合，計算鄰接矩陣的過程用式（2）～（5）表示：

由式（2）計算譜峰差值矩陣絕對值與每個氨基酸殘基的誤差矩陣，如果誤差在給定ε內則將相應元素標記為1，若超出范圍則標記為0，然后將所有矩陣相加得到當前譜的鄰接矩陣；加入相同維度的單位矩陣作為節點自身的信息，避免構圖時譜峰為孤立峰，即譜中沒有相鄰位點產生的同類型離子譜峰，導致不存在邊與之相連造成信息丟失；再計算度矩陣并對鄰接矩陣進行歸一化。

將質譜數據處理成圖結構化數據是denovo-GCN 與DeepNovo處理質譜數據的不同之處，在DeepNovo中將串聯質譜數據對應成規則的歐氏數據，質荷比維度的數據按照質量精度0.01 Da（Dalton）進行擴展：假設譜中的最大質荷比為1 500.00 Da，整個譜離散化為150 000 個刻度，再將每個譜峰相對強度填入離散化后的刻度位置，卷積提取特征。而在denovo-GCN 中，譜峰之間的關系直接計算確定，不需要通過深度學習模型來學習這種關鍵信息。

1.3 譜峰特征構建

denovo-GCN 的另一個關鍵在于為質譜中的每一個譜峰建立特征。由于串聯質譜數據的特殊性，很難在只使用一組離子質荷比和譜峰強度數據條件下推斷出序列信息，因此必須利用肽碎裂產生的離子類型設計特征。肽段在HCD 模式下常見的碎片離子類型有b、y、b2+、y2+、b-H2O、y-H2O、b-NH3、y-NH3、a、a2+、a-H2O、a-NH3等［21］，在計算得到b離子或者y離子質荷比后便可根據母離子質荷比計算同一斷裂位點的其他離子質荷比。在模型中，設定了26種符號標記分別代表20種氨基酸殘基、3 種修飾后的氨基酸殘基、3 種特殊的標記（start、end、pad）。特征可以看作是當前碎裂位點產生的離子與譜峰的距離差值，構建過程如式（6）～（8）：

設ntoken為設定標記的個數，nions為使用的離子類型的種類，計算得到的理論質荷比矩陣為Mt大小為(1，ntoken×nions)，將其按第一維度復制得到Mt'(npeaks，ntoken×nions（)npeaks為譜峰數量）；將當前譜峰矩陣Mo(npeaks，1)，按第二維度復制得到Mo'大小同樣為(npeaks，ntoken×nions)，由式（6）計算譜峰與理論離子的誤差矩陣E，然后通過指數運算將誤差值縮放到區間（0，1）內，⊕代表將譜峰的相對強度Intensity(npeaks，1)拼接到E，形成了最終的特征矩陣F。

1.4 denovo-GCN的模型構建

denovo-GCN 的模型如圖1 所示：由質譜數據分別計算譜圖鄰接矩陣和初始特征矩陣。使用GCN 對質譜數據進行特征提取，按照譜峰的維度加和并使用Leaky ReLU 激活函數進行激活，再使用全連接層輸出，得到氨基酸類型的概率，輸出當前條件下的氨基酸身份。

圖1 denovo-GCN模型Fig.1 denovo-GCN model

新預測的氨基酸加入到序列后，計算下一個位點的特征矩陣，直至出現結束標記或者達到設定的序列最大長度。模型各層的參數大小設置如表1 所示，其中ntoken為設定標記的數目，nions為使用的離子類型的數目。訓練時標注肽序列中的每一個氨基酸作為標記，依次進行批訓練，初始學習率為0.001，根據模型訓練評價自適應調整學習率，最低學習率設置為10-5。由于肽序列中氨基酸出現的頻率差別很大，特別是帶有修飾的氨基酸殘基占比更少，因此在訓練時使用了Focal Loss 函數計算損失，該函數最初用于解決目標檢測中類別不平衡問題［22］。

表1 denovo-GCN模型中各層的參數Tab.1 Parameters of each layer in denovo-GCN

1.5 評價指標

通常從肽水平和氨基酸水平上評價從頭測序結果［12-15］。肽水平召回率和精確率分別為完全預測正確的肽序列占測試數據中所有肽序列的比例和接受的測序結果中肽序列總數的比例，氨基酸水平召回率和精確率分別為預測正確的氨基酸總數分別占測試數據中氨基酸總數的比例和接受的測序結果中氨基酸總數的比例。在氨基酸水平上，從N 端或C 端開始對應位置預測的氨基酸與標注一致則為正確，對于分子量相同的亮氨酸（Leucine，L）和異亮氨酸（Isoleucine，I），在同一位置時認為預測正確。

2 實驗與結果分析

2.1 數據集和模型結構優化

本文在ProteomeTools1（ID：PXD004732）數據集［23］上進行了模型的訓練和測試，確定了模型的結構、離子類型組合和采用的譜峰數量。該數據集來自人工合成蛋白質數據集，從proteomeXchange 蛋白質數據庫中獲得，根據MaxQuant搜索結果以得分score≥100、PIF≥0.7（Precursor Intensity Fraction）過濾選取高質量的肽譜匹配數據，最終得到204 996 條標注數據，并在實驗時以8∶1∶1 的比例隨機劃分訓練集、測試集、驗證集，集合劃分時相互不存在交集。實驗中構建譜峰關系圖時使用的質量誤差ε為0.02 Da。

不同層數的GCN 模型效果根據具體應用會有所差異。本節實驗設置最大譜峰數量為500，離子類型為12種，GCN的hidden size為256，實驗結果如表2 所示：實驗中采用2 層GCN的模型比使用1 層和3 層的模型在肽水平的召回率分別高出2.5個百分點和1.2個百分點，比直接使用全連接網絡高出了2.9個百分點，4層的模型與3層的模型效果基本一致；各組模型氨基酸的召回率在91.19%至92.19%。在氨基酸水平召回率相近的條件下，GCN模型明顯提高了肽水平的召回率，并在使用2 層GCN 結構時獲得最高召回率。因此，后續實驗皆采用2層的GCN結構。

表2 不同GCN層數模型的召回率對比 單位：%Tab.2 Comparison of recall by different GCN layers’models unit：%

2.2 碎片離子類型的選擇

肽段在高能碰撞裂解（HCD）碎裂模式下，主要產生b/y離子及帶二電荷的常規離子，也會產生常規離子失去水分子和失去氨分子的中性丟失離子，以及a型離子。為了測試不同離子類型組合對模型的影響，以b/y離子組合為基礎，測試了加入不同離子類型后的表現，該部分實驗譜峰數量設置為500，實驗結果如表3所示。在加入2電荷的b/y離子后肽召回率比只使用1電荷b/y離子時提升了16.0個百分點，氨基酸水平提升了7.3 個百分點。b、y、b2+、y2+在測序中起著關鍵作用，這與HCD 譜中關鍵離子為b/y離子的特性是一致的。當模型中繼續加入b/y離子的中性丟失離子（b-H2O、y-H2O、b-NH3、y-NH3）時，肽的召回率比使用4種常規離子增加了3.7個百分點，氨基酸水平增加了1.3 個百分點；在加入a型離子及其中性丟失離子（a、a2+、a-H2O、a-NH3）后模型肽水平召回率再次提升了2.1 個百分點。當離子從4 種增加到12 種時，氨基酸水平的召回率只提升了1.9個百分點，但肽的召回率提升了5.7個百分點。這說明，額外增加的8 種離子提供了更多測序信息。當譜中沒有出現某一碎裂位點的常規離子，但存在對應中性丟失的離子峰時，同樣可以為該處氨基酸身份的鑒定提供依據。因此豐富的離子類型組合可以提升測序的準確度。

表3 不同離子類型組合的召回率對比 單位：%Tab.3 Comparison of recall by different combinations of ion types unit：%

2.3 譜峰數量的影響

除離子類型組合會影響模型，每個譜采用的譜峰數量也會對模型產生影響。質譜中存在大量低豐度的離子峰和噪聲峰，基于圖論等其他從頭測序方法中會先對實驗譜消除一部分同位素峰和相對強度過低的峰。在denovo-GCN 中采用簡便的方式，保留相對強度在給定排名內的譜峰。為了驗證譜峰數量的影響，實驗以每個譜選取64 個峰為起始，每次實驗遞增64 個峰，最大峰數為640，實驗結果如圖2 所示。首先統計測序時使用的譜峰數量（used peaks）占全部數據的譜峰數量（total peaks）的變化曲線。當選取256 個譜峰進行實驗時，實驗中用到的譜峰數量占總數據的70.62%，此時譜中的關鍵峰基本納入到了考慮范圍內；選取譜峰數量為384 時占比達到89.63%；選取譜峰數量為512 時占比達到97.39%，接近全部數據。在譜峰數超過256 個時，肽召回率均值為77.84%，模型的準確率趨于穩定。當使用384 個譜峰時，基本將大部分譜峰納入到測序中，且使用384 個譜峰時訓練時間比使用512 個譜峰時減少了1/3，若考慮使用全部譜峰時可選擇512個譜峰。

圖2 肽水平的召回率隨譜峰數量的變化曲線Fig.2 Curve of peptide-level recall varying with number of spectral peaks

2.4 不同測序方法在ProteomeTools1數據上的對比

在確定了模型結構、離子類型組合、譜峰數量后在ProteomeTools1 數據集上對denovo-GCN（12 種離子類型，384個譜峰）、DeepNovo（version 0.0.1）、pNovo（version 3.1.3）、Novor（DeNovoGUI version 1.9.6）進行了測試。上述工具給出了預測肽序列的得分，將最終結果按照得分從小到大排序，給定分數t，計算肽水平的精確率（得分至少為t的實際正確肽數量/得分至少為t的肽數量）和召回率（得分至少為t的實際正確肽數量/數據中總的肽數量），畫出肽水平上的精確率-召回率（Precision-Recall，PR）曲線如圖3所示。

從圖3 可看出，denovo-GCN 的曲線明顯高于DeepNovo、Novor 的曲線，召回率在區間［0，0.5］內與pNovo 的曲線有重合的部分，召回率超過0.5時明顯高于pNovo。再分別計算各PR曲線下的面積，denovo-GCN 為0.731 8，DeepNovo 為0.613 8，pNovo為0.619 2，Novor為0.518 1。denovo-GCN 在同一數據上的測序性能要優于其他三種工具。

2.5 不同物種數據的交叉對比

在實際應用中，從頭測序更多的是解決未知物種蛋白的測序。因此，為了進一步驗證denovo-GCN 的測序表現，本節采用了DeepNovo 中的9 個HCD 數據集，進行物種間的交叉對比實驗，數據信息如表4所示。

表4 9個HCD數據集信息Tab.4 Information of 9 HCD datasets

每次使用其中的8 個數據集混合劃分訓練集、驗證集進行模型訓練，集合之間不存在肽序列交集，未參與模型訓練的1 個物種數據作為測試集。用相同的數據分別訓練DeepNovo和denovo-GCN（12 種離子類型，384 個譜峰），Novor 和pNovo直接使用其提供的軟件進行測序，測試結果如圖4所示。

圖4 denovo-GCN、Novor、pNovo、DeepNovo在9個HCD數據集上的實驗結果對比Fig.4 Experimental result comparison of denovo-GCN，Novor，pNovo，DeepNovo on 9 HCD datasets

圖4（a）是不同工具間氨基酸水平的召回率對比，denovo-GCN 比Novor 高出6.2～32.7 個百分點，比pNovo 高出7.6～14.9 個百分點，比DeepNovo 高出4.3～9.9 個百分點。圖4（b）在不同工具上氨基酸水平的精確率對比，denovo-GCN 比Novor 高出3.8～31.1 個百分點，比DeepNovo 高出4.1～10.0 個百 分 點，而pNovo 在H.sapiens 數 據、M.musculus 數 據、Candidatus 數據上比denovo-GCN 的精確率高出6.1 個百分點、3.7 個百分點、2.4 個百分點，其余數據上denovo-GCN 比pNovo 高出2.2～4.9 個百分點。圖4（c）在肽水平上不同工具的召回率對比，denovo-GCN 的肽的召回率比Novor 的高出9.8～21.1 個百分點，比pNovo 高出4.0～13.0 個百分點，比DeepNovo 高出2.1～10.7 個百分點。綜上實驗結果denovo-GCN相較于Novor、pNovo、DeepNovo，能夠測得更多的氨基酸，并且能夠轉化成更多正確的肽序列，測序能力超過了其他三種工具。相較于DeepNovo的模型結構，denovo-GCN模型更為精簡，使用圖來表達譜峰之間關系并結合圖卷積神經網絡的方式比CNN和LSTM模型在串聯質譜測序上更具優勢。

對于表4的9個物種的測試數據在以ProteomeTools1數據訓練的denovo-GCN、DeepNovo 模型上分別進行測試，并使用ProteomeToolsV2（ID：PXD010595）［24］人工合成肽的數據作為相似物種進行對比參照，得到的結果如表5 所示。在相似物種上兩個模型表現都要好于非同類物種的表現，而非同類物種上由于蛋白質差異，測序效果存在一定差距。這兩部分實驗結果表明denovo-GCN 的測序能力優于DeepNovo、pNovo、Novor。

表5 ProteomeTools1數據集訓練的模型在9個HCD數據集上的實驗結果 單位：%Tab.5 Experimental result on 9 HCD datasets by the models trained on ProteomeTools1 dataset unit：%

2.6 denovo-GCN與pNovo預測結果及序列分析

pNovo 是基于圖論的從頭測序的代表，在幾個測試數據上氨基酸召回率雖然低于DeepNovo，但肽的召回率卻與之接近。為了查看預測序列中出現的錯誤肽序列，在ProteomeTools1 測試數據上pNovo、denovo-GCN 測序結果和數據庫搜索結果之間的關系如圖5 所示：兩者有12 661 條數據測序結果相同，同時互有無法給出對方正確測序結果的數據，但denovo-GCN較pNovo多鑒定出了1 451條。

圖5 pNovo、denovo-GCN、數據庫搜索結果的文氏圖Fig.5 Venn diagram of pNovo，denovo-GCN，database search results

對兩者測序均為錯誤的結果進行分析，總結了測序時出現頻率較高的3 種錯誤類型，示例如表6 所示：1）當串聯質譜中的低質量區域，存在較多的亞胺離子和內部離子，而關鍵性的低質量常規離子峰與之不易區分甚至缺失，在構圖時會出現多條互相連接的邊，氨基酸位置難以確定；2）氨基酸殘基存在單個氨基酸分子量等于兩個小質量氨基酸之和或者兩種不同氨基酸分子量之和兩兩相等的情況，譜中兩端缺失了關鍵的b/y離子；3）在長序列譜或低質量譜中，離子峰更為復雜，譜峰可以對應多種氨基酸序列的組合，在測序時較難得出正確氨基酸組合。這也能夠解釋denovo-GCN 在不同物種數據實驗中氨基酸的召回率能夠達到60%以上，而肽序列的正確率卻在25%～48%。解決上述問題最直接的方法是提升質譜儀輸出數據的質量，而當前質譜數據條件下，解決上述問題的思路主要有兩個：1）算法模型輸出多個候選肽序列并進行再次打分，找出更優的序列；2）不斷探索創新測序算法，從而提高肽序列的正確率。

表6 pNovo與denovo-GCN結果中典型的序列錯誤示例Tab.6 Examples of typical sequence errors in pNovo and denovo-GCN results

3 結語

denovo-GCN 將質譜數據轉化為圖結構數據，根據肽碎裂產生的離子類型對每個譜峰點進行特征設計，將圖卷積神經網絡引入到從頭測序任務中，提升了串聯質譜測序的準確率，超過了基于圖論的從頭測序方法Novor、pNovo，以及基于CNN和LSTM 模型的DeepNovo。實驗結果表明充分利用肽碎片離子類型，選擇適當譜峰數量作為參數可以取得較為理想的效果。雖然denovo-GCN 實驗中同數據集上可以達到數據庫搜索結果70%的肽召回率，并且在不同物種測序上也好于其他工具，但不同物種數據的測試結果并未超過數據庫結果的50%。denovo-GCN 的測序效果會受到訓練數據的影響，可以通過擴大訓練數據種類來消除部分影響。提升從頭測序的準確性，仍是一項值得持續研究的課題，而另一方面，如何測定序列中修飾后的氨基酸類型也需要進一步研究。