李大舟,陳思思,高 巍,于錦濤
(沈陽化工大學 計算機科學與技術學院,遼寧 沈陽 110142)
近年來,藥物發現的技術和水平在不斷進步,促進了生物制劑技術和生物制藥開發的不斷發展。藥物發現是人類發現潛在新型藥物的過程,一般是通過將化合物庫、天然物質或提取物的合成小分子在完整細胞或整個生物體上進行表型篩選,從而識別在過程中具有理想治療效果的物質[1]。由于藥物發現的進步,制成的藥劑使得許多的疾病得以預防和治療。然而,由于目前醫學水平的限制,仍有許多疾病無法得以攻克,并且不斷有新型的病癥出現,所以進行新型藥物的研究和開發的需求十分迫切。
藥物中包含了特定的化合物分子,人體內的大部分化學反應都有蛋白質的參與,因此,掌握化合物-蛋白質相互作用(compound-protein Interaction,CPIs)在藥物發現上有著重要的作用,研究人員可以通過CPI識別篩選出有效的化合物,并且可以了解藥物產生副作用的原因。然而,通過生物實驗的方法來確定CPI十分耗時且費用高昂[2]。人類已知的蛋白質類型和化合物類型眾多,若通過生物實驗的方法來一一驗證它們之間是否存在相互作用,這幾乎很難完成的。因此,人們提出通過計算預測方法輔助CPI的研究,讓計算機來分析數據并進行預測,進而提高藥物發現的速度。
隨著人工智能的快速發展,機器學習(machine learning)已經應用于生活中的不同領域。使用傳統機器學習識別CPI的研究在不斷進步。2004年,Bredel和Jacoby[3]提出了一種從化學基因組學角度開發的預測方法,在統一的模型中同時考慮化合物和蛋白質的信息。在此之后,各種基于此想法的CPI預測模型不斷被提出。例如,在2008年,Jacob和Vert[4]利用化學結構和蛋白質家族之間的張量積作為特征,應用成對核的支持向量機來預測CPI。在2009年,Bleakley和Yamanishi[5]提出二部局部模型(BLM),利用化學結構和蛋白質的氨基酸序列之間的相似性度量,應用具有已知相互作用的支持向量機來預測CPI。為了降低化學基因組學空間的維度,在2012年,Cheng[6]提出使用特征選擇技術,使用選擇后的特征訓練支持向量機。在2013年,Tabei和Yamanishi[7]提出使用哈希算法改進線性支持向量機的預測性能,一次獲得化合物-蛋白質對的指紋。
傳統的機器學習往往由多個獨立的模塊組成,需要多個處理步驟,并且每一步的結果會影響下一步驟的好壞,而端到端的深度學習模型可以自動學習特征,且擁有學習海量數據的能力和強大的擬合能力,只需在輸入端輸入原始數據,模型自動在中間層提取數據的特征,最后在輸出端得到預測結果。在2016年,Kipf等人[8]提出圖神經網絡(graph convolutional network,GCN),該網絡能夠處理具有廣義拓撲圖結構的數據,目前主要應用于圖分類[9]、文本分類[10]、推薦系統[11]、疾病預測[12]等。在2018年,?ztürk等人[13]提出DeepDTA模型,利用卷積神經網絡(convolutional neural network,CNN)提取化合物和蛋白質的特征,然后將兩個特征向量拼接起來,經過全連接層輸出CPI二分類結果。在2019年,?ztürk等人[14]提出WideDTA模型,該模型類似于DeepDTA模型,不同之處是利用了兩個額外的特征以改善模型的性能,兩個特征分別是配體最大公共結構(LMCS)和蛋白質基序和結構域(PDM)。同年,Tsubaki等人[15]和Nguyen[16]分別提出CPI-GNN模型和GraphDTA模型,分別使用圖神經網絡(graph neural network,GNN)和圖卷積網絡(graph convolutional network,GCN)學習化合物分子圖的表示。在2019年,Schwaller等人[17]提出Transformer可用于化學反應預測,但是,仍局限于seq2seq任務。2019年,Yang等人[18]提出了XLNET模型,其基于自回歸(autoregressive,AR)語言模型實現了新的雙向編碼,考慮到在訓練過程中屏蔽的單詞與未屏蔽的單詞之間的關系。受XLNET在兩個序列之間獲得特征的強大能力的啟發,該文提出基于改進Attention Mask編解碼器模型,將化合物和蛋白質當作兩種類型的序列輸入到該模型中,最終得到化合物和蛋白質是否相互作用的預測結果。
該文提出的基于改進Attention Mask編解碼器的化合物與蛋白質預測模型的主體結構如圖1所示。

圖1 模型的基本框架
首先對原始化合物數據進行處理,得到原子矩陣和鄰接矩陣,然后根據關系矩陣得到化合物的分子表示矩陣;同時對原始蛋白質數據的氨基酸序列進行處理,使用Item2vec技術得到蛋白質嵌入矩陣;將蛋白質嵌入矩陣輸入到編碼器,最后將化合物和經編碼器處理的蛋白質表示矩陣輸入到解碼器中,得到相互作用的向量,經過線性變換,最終得到化合物和蛋白質是否相互作用的預測結果。
1.1.1 蛋白質數據處理
蛋白質是構成細胞的基本有機物,氨基酸脫水縮合組成多肽鏈,多肽鏈經過折疊后組成的具有空間結構的物質就是蛋白質。蛋白質序列可以通過其物理性質或其氨基酸序列進行編碼表示[19]。文中蛋白質原始數據表現形式為氨基酸首字母縮寫字符串,根據生物信息學信息可知,蛋白質序列由20種基本氨基酸組成,部分氨基酸英文名、中文名稱和首字母縮寫實例如人免疫球蛋白,其蛋白質氨基酸序列表示為“MEFGLSWVFLVAILEGVQCEVQLVESGGGLVQPGGSLRL SCAASGFTFSSHWMTWVRQTPGKRLEWVANVKQD GSARYYADSVRGRFTISRDNAKNSLYLQMDSLRADD TAVYYCARSTGIDYWGQGTLVTVSS”。
Item2vec是由Barkan[20]提出的一種用于學習和描述復雜句法和語義單詞關系的分布式向量表示技術,借鑒于Word2vec[21]的skip-gram with sampling (SGNS)的思路,將其運用于基于物品的協同過濾(item-based CF)上。Item2vec把原來蛋白質數據的高維稀疏的表示方式映射到低維稠密的向量空間中,然后用這個低維向量來表示該蛋白質,對于大量的蛋白質序列數據,可以通過Item2vec學習蛋白質序列的嵌入式表示,大大簡化下游建模。
基于前人的工作,該文將UniProt[22]中的人類蛋白質氨基酸序列進行預處理,作為一個語料庫,然后使用Item2vec訓練語料,設置蛋白質嵌入向量維度為128維,經過20輪迭代訓練了蛋白質嵌入模型。例如人免疫球蛋白中氨基酸序列長度為132,將其帶入訓練后的蛋白質嵌入模型中,通過嵌入算法將每一個氨基酸轉換為向量,對應一個長度為128的向量,最終人免疫球蛋白表示為大小為(132,128)的矩陣形式。
PCA是由Pearson[23]提出的一種統計方法,主要思想是將原始數據沿最大方差方向投影,得到原始數據的低維特征表示,從而實現數據的降維。通過PCA方法得到蛋白質嵌入向量實現,實現蛋白質特征維度的轉變,化合物的特征維度變換同理。以人免疫球蛋白為例,輸入的表示矩陣大小為(132,128),經PCA處理后的表示矩陣大小為(132,64)。
1.1.2 化合物數據處理
化合物是由兩種或兩種以上的元素組成的純凈物。簡化分子線性輸入規范(simplified molecular-input line-entry system,SMILES)是一種用于輸入和表示分子的線性符號,使用ASCII字符串來描述分子結構。文中化合物原始數據表現形式為SMILES字符串,例如吩噻嗪,其SMILES格式為C1=CC=C2SC3C=CC=CC3NC2=C1。依據化學特性劃分原子特征,原子特征列表如表1所示,每種原子的特征可以使用34維的向量表示。

表1 原子特征列表

續表1
RDKit是開源化學信息學與機器學習的工具包,支持機器學習方面的分子描述符的產生。該文通過使用RDKit封裝的函數對SMILES格式的化合物數據進行讀取和處理,得到化合物的原子矩陣和帶自環的鄰接矩陣,然后利用關系矩陣,得到分子的矩陣表示。
該文提出的模型沿用了經典的Encoder-Decoder結構,使用到并行化計算的自注意力機制,極大地縮短了訓練時間。該整體架構如圖2所示,其中編碼器部分主要由多頭自注意力層和前饋神經網絡層組成,解碼器部分主要由Attention Mask層、編碼器-解碼器注意力層和前饋神經網絡層組成。

圖2 編-解碼器架構
1.2.1 編碼器
編碼器的結構如圖2左側虛線框內所示,由3個編碼器塊堆疊而成,每一個編碼器塊都由兩個子層組成,并且每一個子層之間都使用了殘差連接和層歸一化操作。
編碼器的第一個子層是多頭自注意力層。自注意力的本質是通過當前詞來引入上下文的信息,以此增強對當前詞的表示。首先根據輸入的化合物的原子序列,通過線性變換得到Q、K、V的向量表示,然后根據公式1計算注意力值。
(1)

多頭自注意力層是包括了多個按比例縮放的自注意力層,可以在不改變參數量的情況下增強注意力的表現力,擴展模型專注不同位置的能力。多頭自注意力是對Q、K、V進行分組計算注意力值,如公式2所示,然后拼接所有注意力頭,計算過程如公式3所示。
(2)
MultiHead(Q,K,V)=Concat(head1,…,headn)Wo
(3)

編碼器的第二個子層是前饋神經網絡層。前饋神經網絡層(feed forward layer,FFL)是由兩層全連接神經網絡組成的,選擇ReLU作為激活函數,如公式4所示。該網絡層對注意力的輸出進行空間交換,增加了模型的表現能力。
FFN=Max(0,X*W4+b4)W5+b5
(4)
式中,X表示經多頭自注意力層的輸出矩陣,W4和W5表示權重矩陣,b4和b5表示網絡的偏置。
由于網絡不斷加深,數據的分布也在不斷地發生變化,同時可能會帶來梯度消失或爆炸等問題。加入殘差連接可以從一定程度上緩解因為梯度爆炸導致的網絡退化問題,而加入層歸一化可以保證數據的穩定分布,同時可以加速模型的優化速度。殘差連接和層歸一化操作如公式5所示。
Output=LN(X+(SubLayer(X)))
(5)
式中,X表示每個子層的輸出,SubLayer()表示子層本身的輸出,LN表示Layer Normalization,Layer Normalization的計算公式如下:
(6)
式中,μ、σ分別表示均值和方差,α表示縮放參數,β表示平移參數。
1.2.2 解碼器
編碼器的結構如圖2右側虛線框內所示,由3個解碼器塊堆疊而成,每一個解碼器塊都由三個子層組成,與編碼層一樣,每一個子層之間同樣使用了殘差連接和層歸一化操作。
解碼器的第一個子層是改進的Attention Mask層。傳統的自回歸模型的缺點是不能同時利用上文或者下文的信息,而傳統的自編碼模型的缺點是會導致預訓練階段和微調階段出現不一致的問題。改進的Attention Mask層部分避免了二者的缺點,在傳統的自回歸模型的模式下,引入全排列語言模型(permutation language modeling,PLM),保持當前詞的位置不變,對文本中的其他詞進行重新編排,使得當前中心詞的下文也有可能出現在中心詞的上文中,然后將句尾的一定量的詞進行遮掩,使用自回歸方式預測被遮掩的詞。全排列語言模型的優化目標最大似然化概率如公式7所示。
(7)
式中,T表示序列長度,ZT表示所有可能出現的排列序列,zt表示第t個元素。例如存在一個長度為T的序列,從序列的所有可能的排列序列中隨機采樣一個,然后通過計算來分解聯合概率成條件概率,并加權求和得到預測當前詞概率最大的參數θ,由此捕獲雙向的語境。具體的PLM操作是通過雙流自注意力機制實現的,雙流自注意力機制由內容流注意力機制和查詢流注意力機制組成,同時引入了兩個隱狀態,分別是內容隱狀態hzt和查詢隱狀態gzt。雙流注意力機制的計算過程如公式8和公式9所示。
(8)
(9)
式8中,上標m表示層數,Q值、K值和V值分別代表注意力機制中的查詢向量、鍵向量和值向量,zt表示z∈ZT的前t-1個元素。

例如原本輸入的句子是“1,2,3,4”,若經過PLM操作后的排列序列為“3,2,4,1”,表明在預測“2”的時候,可以看到上文的“3”的信息;當預測“4”的時候,可以看到上文“3”和“2”的信息,并以此類推。內容流和查詢流掩碼矩陣如圖3(c)右圖所示,通過掩碼矩陣,將句子改成隨機的排列組合,實現同時利用上下文信息預測當前詞。

圖3 雙流自注意力機制計算過程圖示
解碼器的第二個子層是編碼器-解碼器注意力層,它的輸入由兩部分構成,分別是掩碼多頭注意力層的輸出Q和編碼器的輸出K、V,通過注意力機制增強對當前詞的表示,并提取編碼器和解碼器間的交互信息。
解碼器的第三個子層是前饋神經網絡層,其工作原理與編碼器中的前饋神經網絡層一樣。該子層的輸入為編碼器-解碼器注意力層的輸出。最后該層的輸出是化合物和蛋白質相互作用的特征向量,將其經過softmax函數,最終得到化合物和蛋白質是否相互作用的概率。
本實驗在windows10系統下進行,使用Intel@i5-8265U作為計算單元,內存為8 GB。模型使用Pytorch框架進行搭建,版本為1.6.0+cu101。構建模型所用的代碼使用到RDKit庫。
文中用于訓練Item2vec模型的蛋白質數據來自于UniProt蛋白質數據庫[22]。選取UniProt蛋白質數據庫中Swiss-Prot子庫里的所有人類蛋白質序列作為一個語料庫,源數據格式如表3所示,總計20 413條,提取蛋白質的氨基酸序列數據,使用該數據對Item2vec模型進行預訓練,學習蛋白質的嵌入式表示。
文中化合物和蛋白質數據主要來源于Lifan[24]構建的GPCR標簽反轉數據集,據實驗驗證,標記反轉實驗可以有效地評估隱藏的配體偏差對模型的影響,降低基于化學基因組的化合物和蛋白質相互作用任務的常見風險。GPCR數據集主要有化合物信息、蛋白質信息和表示是否相互作用的布爾值,數據集包含了356種蛋白質和5 459種化合物的15 343種作用對。
對于GPCR組,隨機選擇500個配體,并將所有涉及這些配體的CPI負樣本匯集在一起。另外,選擇了500個配體,并將所有相關正樣本匯集在一起。在實驗設計后,最終建立了1 537個相互作用的GPCR測試集,剩余的數據集被用來確定超參數。
實驗中,采用二分類交叉熵損失函數、ROC曲線下面積(AUC)以及精度-召回率曲線(PRC)作為模型的評估指標。
二分類交叉熵(Binary Cross Entropy):是多分類softmax_cross_entropy的一種特殊情況,當只有兩類標簽時,即0或者1,使用邏輯回歸的損失函數,如公式10所示。

(10)

ROC(Receiver Operating Characteristic):以假正例率(FPR)為X軸、真正例率(TPR)為Y軸繪制的反映模型敏感性和精確性的趨勢走向的曲線。
AUC(Area Under Curve):ROC曲線下的面積。若分類器的性能越好,則AUC值越接近1。
PRC(Precision Recall Curve):以查全率(Recall)為X軸、查準率(Precision)為Y軸繪制的圖,可以對分類器的整體效果進行綜合評價。該評估指標引入“平衡點”(BEP)概念,當查全率等于查準率時取的值越大時,表明該分類器的性能越好。
該模型的編碼器和解碼器的層數各為3層,多頭注意力頭數為8個,經PCA處理后的蛋白質表示和原子表示的維度為64,編碼器和解碼器完全連接的前饋神經網絡層中隱藏單元數量為512,Dropout為0.2,學習率為1e-4,批尺寸大小為64。
該文使用GPCR測試集在模型上進行訓練,采用接收機工作特性曲線下面積(AUC)、準確召回率曲線(PRC)作為模型的評估指標。
從圖4中可以看出,隨著迭代次數的增加,模型的Loss值在逐漸變小,且愈加接近飽和,在迭代50輪前,模型的AUC值和PRC值的變化明顯,隨著迭代次增加,模型訓練愈加接近飽和,評估指標趨于平緩和穩定,模型的最優AUC值和PRC值分別為0.865和0.883。

圖4 模型訓練Loss變化和AUC值、PRC值變化
該文對模型進行調參試驗,采用控制變量法進行調參。實驗設置如表2所示,實驗結果如圖5所示。

表2 模型對比實驗設置

圖5 不同batchsize下Loss值、PRC值和AUC值變化
從圖5中可以看出,在同樣的迭代次數下,batchsize為64時,模型的Loss值相對于另外兩個更低,batchsize為96時效果較差。在batchsize為64時,模型的PRC和AUC值優于另外兩種情況,batchsize為96時效果較差。
該文對模型網絡結構也進行了對比實驗。實驗設置如表3所示,實驗結果如圖6所示。

表3 結構對比實驗設置

圖6 拓寬網絡和加深網絡的PRC和AUC變化
從圖6可以看出,在相同的迭代次數下,拓寬網絡的PRC和AUC相較于原始網絡在一開始處于較為落后的趨勢,后來逐漸接近;在相同的迭代次數下,加深網絡的PRC和AUC一直處于落后的趨勢。在其他設置不變的情況下,原始網絡的參數設定的PRC和AUC達到最優的情況。
該文選擇了經典的機器學習模型和兩種流行的行業常用模型與該文提出的模型進行對比,實驗結果如表4所示。在GPCR測試集上,該文提出的模型在AUC和PRC方面均優于對比的模型,在數據集上取得了較好的性能,AUC值和PRC值均有提升,表明該模型具有更強的學習蛋白質和化合物之間相互作用的能力。

表4 GPCR測試集性能
該文嘗試將深度學習技術應用于CPI預測的研究中,將該任務轉換成標簽二分類的問題進行解決。在使用傳統的編解碼器模型的基礎上,在解碼器中使用到改進的Attention Mask層,以此來處理蛋白質和化合物二分類任務。在AUC和PRC指標測評下,與其他模型相比,該文改進的模型在實驗上擁有更好的性能表現。
實驗結果表明,該模型可以學習期望的CPI特征,性能更穩定且準確率更高。如果將該模型應用于實際的藥物發現研究中,可以為藥物靶標選擇提供一定的參考價值,加快藥物發現的進程。同時深度學習不要求具備生物學和藥理學等專業知識,就可以得到數據背后的隱藏信息,且對于數據量特別大的數據具有明顯的優勢。然而,該模型構造了一個注意力矩陣,需求與輸入呈平方關系,因此,對內存和算力的需求非常高。