基于深度學(xué)習(xí)的“藥物-靶點(diǎn)”親和力預(yù)測方法設(shè)計

李添添，王俊杰

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

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

1 研究方法

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

圖1 研究方法框架

1.1 藥物（化合物）的獨(dú)熱編碼機(jī)制

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

1.2 雙向LSTM 預(yù)訓(xùn)練氨基酸編碼機(jī)制

使用預(yù)先訓(xùn)練的多層BiLSTM 獲得氨基酸序列的向量表征。首先，對一條氨基酸序列（r1,r2,…,rN），biLSTM 語言模型分別使用M個堆疊的LSTM 網(wǎng)絡(luò)從前向和后向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)］。設(shè)定M=2，即使用2 組雙向LSTM 編碼氨基酸序列，將其中每個LSTM 的隱藏單元設(shè)定為32。BiLSTM 預(yù)訓(xùn)練模型的網(wǎng)絡(luò)結(jié)構(gòu)如圖2 所示。為了充分預(yù)訓(xùn)練BiLSTM 模型，從STRING 數(shù)據(jù)庫［16］中收集了66 235 條氨基酸序列。在預(yù)訓(xùn)練結(jié)束后將BiLSTM 的權(quán)重凍結(jié)，在下游任務(wù)預(yù)測“藥物-靶點(diǎn)”親和力時就不會改變其權(quán)重。

圖2 BiLSTM 預(yù)訓(xùn)練模型網(wǎng)絡(luò)結(jié)構(gòu)

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

1.3 預(yù)測網(wǎng)絡(luò)模塊

預(yù)測網(wǎng)絡(luò)模塊包含 2 個卷積神經(jīng)網(wǎng)絡(luò)（convolution neural network，CNN）模塊及4 個全連接（fully connected，F(xiàn)C）層。2 個CNN 模塊分別用于提取蛋白質(zhì)和化合物的特征，4 個FC 層用于根據(jù)CNN 提取的特征預(yù)測蛋白質(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é)點(diǎn)，第3 層包含512 個神經(jīng)元節(jié)點(diǎn)，第4層僅包含1 個神經(jīng)元節(jié)點(diǎn)。為了防止出現(xiàn)過擬合問題，在每個FC 層后面添加3 個速率為0.1 的隨機(jī)失活（dropout）。計算蛋白質(zhì)與化合物的相互作用分?jǐn)?shù)，并通過Sigmoid 激活函數(shù)功能將相互作用分?jǐn)?shù)調(diào)整成0～1 的數(shù)值。

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

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

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

1.5 實(shí)驗(yàn)環(huán)境

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

1.6 模型評價

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

1.6.1 均方誤差

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

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

1.6.2 一致性指數(shù)

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

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

2 結(jié)果

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

表1 4 種算法在2 個數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果

3 討論

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

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

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

4 結(jié)語

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