李添添，王俊杰

藥物發(fā)現(xiàn)是從化學數(shù)據(jù)庫中為特定的疾病尋找新藥并在臨床試驗中驗證其有效性和安全性的過程。這個過程通常需要10 年以上，期間消耗大量的時間和人力成本，但其結果往往不盡人意［1］。事實上，根據(jù)湯森路透生命科學事業(yè)部一份報告，2008-2010 年，108 種新的或重新利用的藥物中，有51%由于療效不佳而被宣告為失?。?］。這份報告提出以下兩點在藥物研發(fā)過程中十分重要：一是選擇新的及更合適的藥物靶點；二是在藥物研發(fā)的最初階段，盡可能篩選出療效良好的藥物。因此，預測藥物和靶點之間的相互作用具有重大意義。

然而長期以來，“藥物-靶點”相互作用的預測被認為是一個簡單的二分類問題［3-4］，即僅預測藥物與靶點之間是否存在相互作用，而很少對它們關系的親和力值進行評價。親和力值可以提供藥物與靶點相互作用的強度信息，能夠對候選藥物做出更為全面的評價［5］。目前在“藥物-靶點”親和力的預測任務中，Kronecker 正則化最小二乘［6］（Kronecker regularized least squares，KronRLS）是一種基于相似度的方法，即采用不同類型的藥物相似度和蛋白質(zhì)相似度評分矩陣作為特征，將“藥物-靶點”親和力的預測問題表述為一個回歸或秩預測問題；SimBoost［7］是一種新穎的使用梯度增強回歸樹的非線性方法，該方法同樣使用相似的矩陣和構造特征，其訓練數(shù)據(jù)的定義類似于KronRLS 方法。這兩種方法均是基于特征工程的傳統(tǒng)機器學習方法，但其預測結果的準確率仍不盡人意。得益于深度學習在圖像處理和語音識別的成功應用［8］，深度學習方法也被廣泛應用于生物信息學領域，如基因組學研究［9］和藥物發(fā)現(xiàn)［10］。深度學習的主要優(yōu)勢在于通過在每一層的神經(jīng)網(wǎng)絡中進行非線性轉換，可以更好地表示原始數(shù)據(jù)，從而有助于學習數(shù)據(jù)中隱藏的模式［11］。?ztürk Hakime 等［12］在2018年首次提出使用深度藥物-靶點親和力（deep drug-target binding affinity，DeepDTA）的方法預測“藥物-靶點”親和力。該方法使用藥物化學結構的一維表示作為藥物的輸入數(shù)據(jù)，氨基酸序列用于表示靶蛋白的輸入數(shù)據(jù)。但是在該方法中，氨基酸序列使用的獨熱編碼方式僅獨立地描述了每一種氨基酸，并沒有考慮肽鏈的上下游信息，也無法突出哪些氨基酸對靶蛋白有重要的修飾作用。因此，本文將改良以上DeepDTA 方法，構建一種準確率更高的基于深度神經(jīng)網(wǎng)絡的“藥物-靶點”親和力預測方法。

1 研究方法

研究方法概述如下：先對擬定藥物進行獨熱編碼，再通過雙向長短時記憶網(wǎng)絡［13-14］（bidirectional long short-term memory，biLSTM）預訓練語言模型對蛋白質(zhì)（氨基酸序列）進行編碼，隨后將藥物的獨熱編碼和蛋白質(zhì)的編碼通過預測網(wǎng)絡模塊進行深度學習，得出二者的相互作用分數(shù)，最后將預測結果在Davis 激酶結合親和力數(shù)據(jù)集［15］和KIBA大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集［7］上進行驗證。方法框架見圖1。

圖1 研究方法框架

1.1 藥物（化合物）的獨熱編碼機制

從PubChem BioAssay 數(shù)據(jù)庫中收集200 萬個結構多樣的化合物的簡化分子線性輸入規(guī)范（simplified molecular input line entry system，SMILES）序列，篩選出64 個描述符，每一個描述符對應特定的整數(shù)作為SMLIES 的獨熱編碼，如字符“C”對應整數(shù)1，“N”對應整數(shù)3，“O”對應整數(shù)5，“=”對應整數(shù)63，則SMILES“CN=C=O”的獨熱編碼為向量［1,3,63,1,63,5］。

1.2 雙向LSTM 預訓練氨基酸編碼機制

使用預先訓練的多層BiLSTM 獲得氨基酸序列的向量表征。首先，對一條氨基酸序列（r1,r2,…,rN），biLSTM 語言模型分別使用M個堆疊的LSTM 網(wǎng)絡從前向和后向2 個方向計算氨基酸出現(xiàn)的概率，2 個方向的LSTM 分別基于前向和后向語言模型的上下文輸出中間嵌入向量（即隱藏狀態(tài)向量），其中j=1,…,M。再對每一個氨基酸ri使用M層的雙向語言模型計算出2M+1 個嵌入向量E(ri)={hij|j=0,…,M}。然后通過聚合不同層的表示獲得其上下游的信息表示。因此，一條氨基酸序列（r1,r2,…,rN）經(jīng)過雙向語言模型編碼后表示為一組等長的向量Eco(S)=［Eco(r1),Eco(r2),…,Eco(rN)］。設定M=2，即使用2 組雙向LSTM 編碼氨基酸序列，將其中每個LSTM 的隱藏單元設定為32。BiLSTM 預訓練模型的網(wǎng)絡結構如圖2 所示。為了充分預訓練BiLSTM 模型，從STRING 數(shù)據(jù)庫［16］中收集了66 235 條氨基酸序列。在預訓練結束后將BiLSTM 的權重凍結，在下游任務預測“藥物-靶點”親和力時就不會改變其權重。

圖2 BiLSTM 預訓練模型網(wǎng)絡結構

由于使用獨熱編碼的SMILES序列和使用雙向語言模型編碼的氨基酸序列的長度不同，為了創(chuàng)建一個有效的表示形式，設定化合物的SMILES 的最大長度為100，氨基酸序列的最大長度為1 200，超過最大長度的化合物的SMILES序列和氨基酸序列將會被強制截斷為最大長度。

1.3 預測網(wǎng)絡模塊

預測網(wǎng)絡模塊包含 2 個卷積神經(jīng)網(wǎng)絡（convolution neural network，CNN）模塊及4 個全連接（fully connected，F(xiàn)C）層。2 個CNN 模塊分別用于提取蛋白質(zhì)和化合物的特征，4 個FC 層用于根據(jù)CNN 提取的特征預測蛋白質(zhì)和化合物之間的親和力。

1.3.1 CNN 模塊

每個CNN 模塊包含3 個堆疊的一維卷積層，每層卷積使用前一層的輸出作為其輸入。為了避免梯度消失問題，在每個卷積層上附加一個校正線性單元（rectified linear units，ReLU）。

1.3.2 FC 層

在每個CNN 模塊后分別使用2 個全局最大池化層，分別為蛋白質(zhì)和化合物生成高水平的特征向量，然后將2 個全局最大池化層的輸出連接到4個FC 層。在FC 層中，前2 層各包含1 024 個神經(jīng)元節(jié)點，第3 層包含512 個神經(jīng)元節(jié)點，第4層僅包含1 個神經(jīng)元節(jié)點。為了防止出現(xiàn)過擬合問題，在每個FC 層后面添加3 個速率為0.1 的隨機失活（dropout）。計算蛋白質(zhì)與化合物的相互作用分數(shù)，并通過Sigmoid 激活函數(shù)功能將相互作用分數(shù)調(diào)整成0～1 的數(shù)值。

為了訓練給定的神經(jīng)網(wǎng)絡，使用平方根均誤差目標函數(shù)作為損失函數(shù)，使用自適應矩估計算法優(yōu)化網(wǎng)絡參數(shù)［17］，默認學習率為0.01。

1.4 數(shù)據(jù)集介紹

Davis 激酶結合親和力數(shù)據(jù)集包含了激酶蛋白家族和相關抑制劑的選擇性分析及其各自的解離常數(shù)值，里面含有442 種蛋白質(zhì)及68 種化合物。KIBA 大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集起源于一種叫做KIBA 的方法，它將不同來源的激酶抑制劑生物活性結合起來。KIBA 數(shù)據(jù)集最初有467 個目標和52 498 種藥物，經(jīng)過濾后，該數(shù)據(jù)集僅包含具有至少10 種相互作用的藥物和靶點，總共產(chǎn)生229 種獨特的蛋白質(zhì)和2 111 種獨特的藥物。

1.5 實驗環(huán)境

模型訓練實驗環(huán)境中的硬件設施主要為GeForce GTX2080Ti 型 GPU；軟件設施主要為Ubuntu16.04 操作系統(tǒng)及Tensorflow 深度學習框架，其中Keras 的版本為2.2.5，算法實現(xiàn)語言采用Python 3.8。

1.6 模型評價

采用5 折交叉驗證法分別在Davis 激酶結合親和力數(shù)據(jù)集和KIBA大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集上評估本文構建的“藥物-靶點”親和力預測方法的性能。將每個數(shù)據(jù)集中的“藥物-蛋白質(zhì)”對平均分成5 份，選擇其中的4 份作為親和力已知的“藥物-蛋白質(zhì)”對輪流訓練本文提出的模型，將另外1 份作為親和力未知的“藥物-蛋白質(zhì)”對用于預測親和力，根據(jù)預測親和力和真實親和力計算均方誤差和一致性指數(shù)，將其作為評價結果，平均5 次的評價結果為最終的評價結果，然后將該結果與使用KronRLS、SimBoost 和DeepDTA 算法的預測結果進行比較。

1.6.1 均方誤差

使用均方誤差（mean squared error，MSE）衡量預測的“藥物-靶點”親和力值和真實值之間的差距，計算公式如下。

式中，P表示預測值，Y表示真實值，N表示所有樣本的個數(shù)。

1.6.2 一致性指數(shù)

使用一致性指數(shù)（concordance index，CI）衡量“藥物-靶點”親和力預測的性能［18］，計算公式如下。

式中，bi表示大親和力δi的預測值，bj表示小親和力δj的預測值，Z表示一個歸一化常數(shù)，h(x)表示階躍函數(shù)。

2 結果

相較于KronRLS、SimBoost、DeepDTA，本文所構建的方法在Davis 激酶結合親和力數(shù)據(jù)集和KIBA 大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集上均獲得了最高的CI 值和最低的MSE 值，見表1。

表1 4 種算法在2 個數(shù)據(jù)集上的實驗結果

3 討論

對于藥物研發(fā)，找出候選藥物與靶點之間的相互作用強度是至關重要的，而識別“藥物-靶點”的相互作用已成為早期藥物發(fā)現(xiàn)階段的關鍵步驟。但是采用藥物試驗的方法進行藥物篩選既昂貴又耗時。因此，迫切需要構建出能夠以最小的錯誤率來識別潛在的“藥物-靶點”相互作用的方法［19］。

根據(jù)已有研究［20］可知，藥物和蛋白靶點呈三維結構，二者的結合是一個相對復雜的過程。本方法雖然不能完全反映出藥物和靶點結合的復雜性，但三維結構不易獲得，且已有相關文獻論證了使用深度學習如DeepDTA 預測親和力的有效性。本文在藥物與靶點結合的三維結構不易獲取的情況下，借助深度學習強大的非線性建模能力，僅使用蛋白質(zhì)的氨基酸序列和藥物的一維化學結構來預測“藥物-靶點”親和力，具有更強的適用性，其與靜態(tài)氨基酸編碼方式（如DeepDTA）的不同之處在于，預訓練語言模型可以結合相鄰氨基酸的信息動態(tài)對氨基酸序列進行編碼，自動提取更為精細的氨基酸水平特征，讓這些特征可以在不同的氨基酸序列上下游之間有所區(qū)別。此外，為了預測靶蛋白和藥物之間的親和力，本文設計了2 個獨立的CNN 模塊，從原始化合物序列和經(jīng)過預訓練語言模型編碼的氨基酸序列中學習藥物和蛋白質(zhì)的特征，并將這些特征傳送到一個全連接的網(wǎng)絡中來預測親和力。

本文比較了KronRLS、SimBoost、DeepDTA 算法及本方法在 Davis 激酶結合親和力數(shù)據(jù)集和KIBA 大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集上的MSE值和CI 值。其中MSE 值越低、CI 值越高，說明方法預測結果越準確。在Davis 激酶結合親和力數(shù)據(jù)集中，SimBoost 和KronRLS 方法的性能類似；但是在KIBA 大規(guī)模激酶抑制劑生物活性數(shù)據(jù)集中，SimBoost 方法的CI 值高于KronRLS 方法的CI 值。KronRLS 是基于正則化最小二乘法的一種預測方法，其利用藥物和靶點的相似矩陣來獲得模型的參數(shù)值，在預測“藥物-靶點”親和力時僅依賴于藥物和靶點的相似性，無法對復雜的藥物和靶點的相互作用進行很好的預測。SimBoost 是一種基于特征工程的方法，需要專家來定義蛋白質(zhì)和化合物的相關特征。DeepDTA 算法使用深度學習來挖掘蛋白質(zhì)和藥物的特征［21］，因此其CI 值要高于KronRLS、SimBoos 這兩種傳統(tǒng)的方法。雖然DeepDTA 使用深度神經(jīng)網(wǎng)絡來模擬藥物和靶點復雜的相互作用過程，但是其獨熱編碼方法無法充分表達蛋白質(zhì)的氨基酸序列信息。而本方法在2 個數(shù)據(jù)集上均獲得了最高的CI 值和最低的MSE 值，說明本文使用的雙向語言模型學習氨基酸序列信息較DeepDTA 僅使用獨熱編碼的方式表達的信息更為準確，其預測能力要優(yōu)于DeepDTA，且無須專業(yè)人員來定義蛋白質(zhì)和化合物的相關特征，節(jié)約了人力資源及學習成本。

4 結語

預測藥物與靶點之間的親和力不僅可以提供更大的信息量，而且更具挑戰(zhàn)性［22］。本文在藥物與靶點結合的三維結構不易獲取的情況下，借助深度學習強大的非線性建模能力，僅使用蛋白質(zhì)的氨基酸序列和藥物的一維化學結構來預測“藥物-靶點”親和力，且其預測結果的準確率高于KronRLS、SimBoost 和DeepDTA 方法。然而本文未考慮到藥物的分子圖結構信息，下一步研究將嘗試將藥物的分子圖結構應用于“藥物-靶點”親和力的預測模型中，同時補充研究模型的可解釋性，以期獲得更滿意的預測結果。