彭澤佳，張曉龍+

(1.武漢科技大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院，湖北 武漢 430065； 2.武漢科技大學(xué) 大數(shù)據(jù)科學(xué)與工程研究院，湖北 武漢 430065； 3.武漢科技大學(xué) 智能信息處理與實時工業(yè)系統(tǒng)湖北省重點實驗室，湖北 武漢 430065)

0 引 言

智能計算預(yù)測藥物靶點能為藥物研制降低近43%的實驗成本[1]。藥物靶點親和力預(yù)測是判斷藥靶相互作用的關(guān)鍵。藥物靶點親和力預(yù)測領(lǐng)域中，分子對接[2,3]方式無法作用于大規(guī)模數(shù)據(jù)集，基于相似度計算的KronRLS[4]和SimBoost[5]算法依賴特征工程導(dǎo)致相似空間中藥物靶點特征表達(dá)受到限制。為了克服這些缺陷，DeepDTA[6]和Wide-DTA[7]利用1D卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network，CNN)自動學(xué)習(xí)提取藥物靶點有效特征；DeepGS[8]運用雙向門控循環(huán)單元(bi-directional gated recurrent unit，BiGRU)和圖注意力網(wǎng)絡(luò)(graph attention network，GAT)分別提取藥物SMILES序列和藥物指紋(fingerprint(fp))的有效特征，利用1DCNN提取靶點蛋白序列的有效特征。利用深度學(xué)習(xí)模型提取特征是具有一定的優(yōu)勢，但是，這些單模型所能提取到的特征種類會受到限制。

近年來，深度學(xué)習(xí)混合模型在自然語言處理領(lǐng)域應(yīng)用廣泛。1DCNN與雙向長短時記憶(bi-directional long short term memory，BiLSTM)神經(jīng)網(wǎng)絡(luò)所形成的混合模型能夠有效地提取文本數(shù)據(jù)特征。靶點蛋白的線性序列能決定蛋白質(zhì)的三維結(jié)構(gòu)[9]，藥物的高維結(jié)構(gòu)也能從SMILES(simplified molecular input line entry system)序列中被獲取。因此，本文針對藥物靶點親和力預(yù)測單模型提取特征種類受限問題的主要研究內(nèi)容如下：

(1)利用BiLSTM代替?zhèn)鹘y(tǒng)循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network，RNN)和長短時記憶(long short term memory，LSTM)神經(jīng)網(wǎng)絡(luò)，提取藥物靶點序列局部特征的上下文聯(lián)系。

(2)利用CNN和BiLSTM的混合模型構(gòu)建藥物靶點特征提取器，并在策略設(shè)計中加入特征并行提取和特征融合技術(shù)。利用1DCNN提取局部特征，利用BiLSTM提取上下文聯(lián)系，并行提取的局部特征和上下文特征拼接融合為以向量表示的用于藥物靶點預(yù)測的特征信息。

1 全局特征提取策略

本文的全局特征提取策略以深度學(xué)習(xí)算法為特征計算基底，針對不同類型的特征構(gòu)建不同的特征提取器并行提取，分為藥物靶點序列特征向量化、局部特征提取、上下文特征提取、特征融合等4個步驟。

1.1 藥物靶點序列特征向量化

序列分割是藥物靶點序列數(shù)據(jù)特征向量化的第一步。本文采取分子分割方式，即以化學(xué)分子和氨基酸分子為分割單元拆分藥靶序列數(shù)據(jù)。該分割方式習(xí)得的詞向量能夠通過空間嵌入原理反映分子之間的分布特性和相互作用關(guān)系。以空格為間隔，將蛋白質(zhì)序列分割成單個氨基酸分子，將藥物SMILES序列分割成單個化學(xué)分子，分割原則為直接拆分，不改變分子的相對位置和排列順序。具體示例靶點蛋白序列片段“MTVKTEA”如下

MTVKTEA→M T V K T E A

整數(shù)編碼是藥物靶點序列數(shù)據(jù)特征向量化的第二步，致力于建立分割單元與嵌入空間中特征向量的映射聯(lián)系。給分割單元集中每個獨立的分割單元以唯一的整數(shù)表示，整數(shù)0不標(biāo)記任何的分割單元，分割單元集中藥靶分子的選取與DeepDTA[6]保持一致。定長處理分割后序列使其在詞嵌入學(xué)習(xí)過程中獲取有效特征表示。整數(shù)0用于輸入序列不足固定長度時的零填充處理，并映射到詞嵌入矩陣第0行行向量(零向量)。具體示例SMILES序列片段“COC1=C”如下

[C O C 1 = C]→[42 48 42 35 40 42]

詞嵌入式特征向量化是藥物靶點序列數(shù)據(jù)特征向量化的最后一步，即利用Keras框架中的Embedding層將分割編碼后的序列數(shù)據(jù)轉(zhuǎn)化成實數(shù)矩陣，并隨著后續(xù)預(yù)測任務(wù)學(xué)習(xí)更新詞向量，使其具有語義信息。序列數(shù)據(jù)在Embedding層中先轉(zhuǎn)化為獨熱(one-hot)稀疏矩陣，通過線性變換轉(zhuǎn)化為一個低維稠密矩陣，該矩陣即為輸入序列的特征向量矩陣。Embedding層詞向量訓(xùn)練是基于神經(jīng)網(wǎng)絡(luò)的分布式詞嵌入訓(xùn)練方法之一，其本質(zhì)上即為學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)中的一層全連接層參數(shù)，具體原理如圖1所示，V為由分割單元所組成的詞匯表大小，N為人為設(shè)定的詞嵌入維度，全連接層參數(shù)矩陣WV×N={wij} 是詞嵌入矩陣，該矩陣的行向量即為相應(yīng)行索引下標(biāo)所對應(yīng)整數(shù)編碼的分割單元的特征向量，該矩陣是所有分割單元的特征向量矩陣。

圖1 基于神經(jīng)網(wǎng)絡(luò)的詞嵌入原理

1.2 局部特征提取

藥物靶點線性序列由不同功能的子序片段以一定次序組合而成[10,11]。因此，本文選用1D卷積神經(jīng)網(wǎng)絡(luò)(CNN)提取藥物靶點線性序列局部特征。

假設(shè)藥物靶點序列W={w(1),w(2),…,w(i)}(1≤i≤l)，l為輸入序列W中分割單元個數(shù)。首先，將序列W利用Embedding層特征向量化為特征矩陣Sj={e(w(1)),e(w(2)),…,e(w(i))}(1≤i≤l)， 其中e(w(i)) 為對應(yīng)分割單元w(i) 的n維詞嵌入特征向量，e(w(i))∈n，Sj∈l×n； 將Sj作為卷積層輸入，用t個大小為r×n濾波器對特征矩陣Sj執(zhí)行步長s=1的卷積操作，提取Sj的局部特征值fi(1≤i≤l-r+1)， 得到局部特征向量C， 計算如下

(1)

C={f1,f2,…fl-r+1}

(2)

di=max(C)或者di=ave(V)

(3)

最后，將池化后得到的所有特征值di(1≤i≤t) 輸入全連接層(神經(jīng)元個數(shù)為n)進(jìn)行特征綜合提取，最終得到有效的局部特征向量V

V={d1,d2,…,dn}

(4)

1.3 上下文特征提取

大自然中存在著多種復(fù)雜的類似于人類溝通語言的“生命語言”，其完成著細(xì)胞間大部分重要信息的傳遞。其中，蛋白質(zhì)序列形式上雖然為生物序列但組成結(jié)構(gòu)上具有語言特性，存在著上下文依賴聯(lián)系，序列上單個殘基與其之前和之后的所有殘基均相關(guān)，殘基之間存在著復(fù)雜的長距離依賴信息[10]，根據(jù)殘基的上文信息和下文信息，能更準(zhǔn)確地習(xí)得其語義?？捎糜谒幬镒R別的SMILES字符串中蘊藏著豐富的信息，其化學(xué)分子之間和化學(xué)片段之間均存在著類似于文本上下文的語義信息，并可通過詞嵌入等算法習(xí)得其特定的特征模式[11]。

1D卷積神經(jīng)網(wǎng)絡(luò)(CNN)能夠運用局部感受野捕獲有效的局部特征，但是無法捕獲局部特征之間的上下文聯(lián)系。傳統(tǒng)的RNN可以有效挖掘時序方向上的上下文語義聯(lián)系，但其對很久之前信息的感知能力會隨著輸入序列長度的增加而下降，進(jìn)而產(chǎn)生梯度爆炸和長時依賴問題[12]；長短時記憶(LSTM)神經(jīng)網(wǎng)絡(luò)[13]能夠克服RNN的缺陷，但LSTM只能捕獲藥靶序列的前向依賴聯(lián)系，無法提取序列的后向依賴信息；而從LSTM改進(jìn)而來的BiLSTM能夠結(jié)合前向與后向信息。藥物靶點序列片段的語義特性既與其之前的片段信息有關(guān)，也與其之后的片段信息密切聯(lián)系，因此，利用BiLSTM代替RNN和LSTM，既引入了下文信息，也有效避免了傳統(tǒng)循環(huán)神經(jīng)網(wǎng)絡(luò)梯度消失等問題。

利用BiLSTM進(jìn)行局部特征的上下文特征提取。假設(shè)由卷積層提取的藥物靶點局部特征矩陣為Sk={s(w(1)),s(w(2)),…,s(w(i))}(1≤i≤l-r+1)， 其中，Sk∈(l-r+1)×t，t為濾波器數(shù)量。BiLSTM模型，如圖2所示，是由兩個上下疊加的LSTM網(wǎng)絡(luò)組成。將Sk每個時間步上的局部特征向量s(w(i)) 作為每個時刻t的BiLSTM的輸入，同一時刻t存在著兩個方向相反的LSTM門。LSTM門中利用記憶存儲單元捕獲之前長時間段的歷史信息，并利用各種門機(jī)制讓信息選擇式傳遞，在LSTM門中特征選擇方式為：

(5)

(6)

(7)

(8)

(9)

(5)最后獲得t時刻的上文特征信息輸出ht

(10)

其中，Wf、Wi、Wc、Wo分別代表遺忘門的權(quán)重矩陣、輸入門的權(quán)重矩陣、當(dāng)前輸入單元狀態(tài)的權(quán)重矩陣和輸出門的權(quán)重矩陣；bf、bi、bc、bo分別代表遺忘門偏置項、輸入門偏置項、當(dāng)前輸入單元偏置項和輸出門偏置項。

圖2 融合模型原理

(11)

(12)

(13)

BiLSTM在提取上下文特征信息過程中，利用前向LSTM模型提取子序片段上文特征信息，如式(11)所示；利用后向LSTM模型提取子序片段下文特征信息，如式(12)所示；最終，以concat形式拼接獲得包含上下文特征信息的特征矩陣，如式(13)所示。

1.4 特征融合

將局部特征向量和上下文特征向量使用keras的concatenate()方法進(jìn)行融合拼接，如圖2所示，得到藥物或者靶點最終的特征向量。特征融合代碼邏輯如下：

output=keras.layers.concatenate([cnn_output, cnn_bilstm_output],axis=-1)

1.5 全局特征提取策略描述及算法步驟

本文的全局特征提取策略，先將藥物靶點序列特征向量化，再利用局部特征提取器和下文特征提取器并行提取藥靶序列的局部特征和局部特征之間的上下文特征，最后將這兩種互補特征進(jìn)行融合，該策略適用于藥物和靶點蛋白的特征提取。

本文的局部特征提取器由輸入層、詞嵌入層、CNN層、特征存儲融合層、全連接層以及全局最大池化層組成，如圖2左分支所示。首先，由詞嵌入層將分割編碼后的藥物靶點序列特征向量化，詞嵌入維度為128維；后接四層卷積層，使用不同窗口大小的濾波器各128個并行提取藥物靶點序列不同子序片段特征；再將輸入特征矩陣和每層卷積操作所獲取的局部特征矩陣存儲進(jìn)由keras中concatenate()方法搭建的特征融合層，實現(xiàn)不同類型特征的存儲與融合；再經(jīng)Dense操作在每個時間步上對融合特征向量進(jìn)行特征的綜合再提取，可使所提取的特征具有更高層的生物含義；最后，使用全局最大池化層提取每個特征維度上的關(guān)鍵信息，通過取局部特征區(qū)域的最大值實現(xiàn)特征壓縮，得到最終有效的藥物靶點局部特征向量。

本文的上下文特征提取器由輸入層、詞嵌入層、卷積池化層、BiLSTM層、全連接層以及全局平均池化層組成，如圖2右分支所示。首先，由詞嵌入層將分割編碼后的藥物靶點序列特征向量化，詞嵌入維度為128維；后接四層卷積層和一層最大池化層，使用相同窗口大小、不同數(shù)量的濾波器提取藥靶序列局部特征，并下采樣提取關(guān)鍵局部特征；再利用BiLSTM層捕獲局部特征之間的上下文聯(lián)系；再經(jīng)全連接層綜合提取每個時間步上的上下文特征；最后，在特征維度上利用全局平均池化取特征均值，得到最終有效的藥物靶點上下文特征向量。

藥物靶點序列全局特征提取流程如圖3所示，流程解釋如下：

步驟1 從數(shù)據(jù)集中加載原始藥物靶點序列數(shù)據(jù)。

步驟2 藥物靶點序列特征向量化，具體操作見1.1節(jié)。

①以分子為分割單元分割藥物靶點序列數(shù)據(jù)；②整數(shù)編碼分割后序列；③利用Embedding層轉(zhuǎn)化藥靶序列為特征向量矩陣。

步驟3 局部特征和上下文特征并行提取。

A)局部特征提?。孩?DCNN并行提取不同局部片段特征；②concat存儲融合詞嵌入和各類局部特征；③Dense(帶TimeDistributed)進(jìn)行特征的綜合再提??；④GlobalMaxPooling1D提取特征維度上的最大值；

B)上下文特征提?。孩?DCNN+MaxPooling1D提取局部特征；②BiLSTM提取局部特征的上下文特征；③Dense(帶TimeDistributed)進(jìn)行特征的綜合再提??；④GlobalAveragePooling1D提取特征維度上的平均值；

步驟4 concat特征融合局部特征和上下文特征，得到藥物或者靶點的全局特征向量。

步驟5 concat特征融合藥物全局特征向量和靶點全局特征向量，得到藥物靶點對特征向量。

圖3 藥物靶點結(jié)合親和力預(yù)測實驗流程

2 藥物靶點結(jié)合親和力預(yù)測

本文的預(yù)測模型沿用DeepDTA預(yù)測模型部分[6]。將藥物靶點對特征向量輸入由全連接層和dropout機(jī)制組成的預(yù)測模型中得到藥物靶點結(jié)合親和力預(yù)測值，使用均方誤差(MSE)作為損失函數(shù)，如式(14)所示，其中V為真實值向量，Y為相應(yīng)的實際輸出的預(yù)測值向量，N為采樣數(shù)

(14)

3 實驗與分析

3.1 實驗環(huán)境

實驗環(huán)境如下：電腦端的操作系統(tǒng)為windows10，CPU為i5-7500，64位操作系統(tǒng)，8 GB內(nèi)存；服務(wù)器端的GPU資源為8個16 GB Tesla V100 SXM2，高性能服務(wù)器，125 GB運存，個人用戶目錄為3 TB內(nèi)存。

3.2 實驗數(shù)據(jù)

本文使用藥物靶點親和力預(yù)測領(lǐng)域的基準(zhǔn)數(shù)據(jù)集KIBA[14]和Davis[15]評估模型性能，其基本信息和實驗用度見表1。Davis數(shù)據(jù)集包含442條蛋白質(zhì)氨基酸序列和68條藥物SMILES序列，其中69%的藥物靶點相互作用對的結(jié)合親和力值為10 000 nM(pKd=5)預(yù)示著弱相互作用和無相互作用，同時意味著數(shù)據(jù)集中負(fù)樣本數(shù)偏多。所以，為了穩(wěn)定數(shù)值分布，He等[5]將原始的Kd值映射為對數(shù)空間中的pKd值，轉(zhuǎn)化公式如下

(15)

KIBA數(shù)據(jù)集起源于一種命名為KIBA的方法。KIBA方法利用優(yōu)化Ki、Kd和IC50所包含的統(tǒng)計信息一致性來構(gòu)建KIBA值。KIBA數(shù)據(jù)集由229條獨立蛋白質(zhì)氨基酸序列和2111條獨立藥物SMILES序列組成。本文使用經(jīng)過過濾的KIBA數(shù)據(jù)集版本，其中每個藥物和靶點至少發(fā)生10次相互作用[5]。

表1 實驗數(shù)據(jù)集分布

3.3 評價指標(biāo)

本文通過計算Concordance Index(CI)和Mean Squared Error(MSE)指標(biāo)來衡量所提模型對藥靶結(jié)合親和力預(yù)測任務(wù)的性能。CI評估輸出連續(xù)值模型的排序性能

(16)

該指標(biāo)衡量兩個隨機(jī)藥靶對的預(yù)測結(jié)合親和力值是否與它們的真實結(jié)合親和力值有著相同的順序。其中bx是兩者中較大結(jié)合親和力δx的預(yù)測值，by是兩者中較小結(jié)合親和力δy的預(yù)測值，Z為歸一化常量， h(m) 是一個分段函數(shù)[4]

(17)

MSE衡量預(yù)測值向量和真實值向量之間的不同，其原理已經(jīng)解釋于2節(jié)。

3.4 實驗參數(shù)

本文實驗參數(shù)的選取與藥物靶點序列特性相關(guān)，參數(shù)改變會直接影響模型性能。在參數(shù)設(shè)置時采取5折交叉驗證評估模型性能，測試結(jié)果取測試集上的5次平均測試性能。Davis數(shù)據(jù)集在模型構(gòu)造時引入dropout機(jī)制減緩訓(xùn)練時可能存在的過擬合情況。KIBA數(shù)據(jù)集在卷積操作時引入

表2 局部特征提取器參數(shù)設(shè)定

表3 全局特征提取器參數(shù)設(shè)定

表4 實驗超參數(shù)設(shè)定

BatchNormalization機(jī)制加速收斂。具體參數(shù)設(shè)定見表2、表3、表4。

在設(shè)定局部特征提取器參數(shù)時，藥物SMILES序列滑動窗口大小比較[2,4,6,8,12,16]，靶點氨基酸序列滑動窗口大小比較[3,5,7,11,15]，滑動窗口數(shù)量分別比較[64,96,128,256]，dropout分別比較[0.1,0.3,0.5]；在設(shè)定全局特征提取器參數(shù)時，藥物SMILES序列滑動窗口大小比較[4,6,8,12]，靶點氨基酸序列滑動窗口大小比較[5,7,8,11]，過濾器數(shù)量逐層遞加對比實驗，池化窗口大小對比了[2,4]，下采樣因子對比了[1,2]，padding方式對比了[“SAME”,”VALID”]，BiLSTM隱藏大小對比了[64,128,256,512]，Dense隱藏大小隨著LSTM層的隱藏數(shù)對比了[100,200,400,800,1000]；在設(shè)定超參數(shù)時，學(xué)習(xí)率對比了[0.01,0.001,0.0001]，因為KIBA數(shù)據(jù)集的大小近乎為Davis的4倍，故KIBA數(shù)據(jù)集的epoch對比實驗了[100,150,200,300]，選取200時的特征提取效果最佳，Davis的epoch對比實驗了[50,100,150,200]，選取100時的特征提取效果最佳；batch大小對比實驗了[64,128,256,512]。參數(shù)對比實驗時均采取控制變量法。

3.5 實驗結(jié)果及分析

首先，本文將基于深度學(xué)習(xí)的單模型和混合模型進(jìn)行對比實驗，實驗結(jié)果見表5和表6。單CNN模型結(jié)果選自DeepDTA[6]，單BiLSTM模型和CNN-BiLSTM串行混合模型均抽取自本文策略所構(gòu)模型中相應(yīng)部分。

表5 KIBA數(shù)據(jù)集上單模型與混合模型結(jié)果對比

表6 Davis數(shù)據(jù)集上單模型與混合模型結(jié)果對比

如表5所示，在KIBA數(shù)據(jù)集上，本文方法比單CNN模型在平均CI指標(biāo)和平均MSE值上分別提升了2.1%和降低了2.8%，比單BiLSTM模型在平均CI指標(biāo)和平均MSE值上分別提升了1.3%和降低了1.7%，比CNN-BiLSTM串行混合模型在平均CI指標(biāo)和平均MSE值上分別提升了0.4%和降低了0.6%。如表6所示，Davis數(shù)據(jù)集上的本文方法性能均優(yōu)于單CNN模型、單BiLSTM模型和CNN-BiLSTM模型性能。如表5和表6共同所示，本文的深度并行全局特征提取策略，利用CNN提取藥靶序列局部特征，利用CNN-BiLSTM提取上下文特征，將兩種并行提取的互補特征進(jìn)行融合，取得了比單模型和串行混合模型都好的預(yù)測性能，同時能夠適用于不同的藥物靶點數(shù)據(jù)集。

其次，本文所提的全局特征提取策略還與藥物靶點親和力預(yù)測中的經(jīng)典以及前沿模型結(jié)果進(jìn)行對比，實驗結(jié)果見表7和表8。KronRLS[4]算法和SimBoost[5]算法引入特征工程構(gòu)建的相似度矩陣預(yù)測藥物靶點親和力，有效解決了傳統(tǒng)分子對接[2,3]難以作用于大規(guī)模數(shù)據(jù)集的問題，但是藥物靶點在相似空間中的特征表達(dá)受到限制，實驗結(jié)果顯示本文方法及其它深度學(xué)習(xí)方法取得了更好的預(yù)測結(jié)果；DeepDTA[6]模型利用1D卷積神經(jīng)網(wǎng)絡(luò)(CNN)自動學(xué)習(xí)提取藥物靶點特征，有效解決了人工提取特征的缺陷，提高了預(yù)測性能，但是CNN僅能提取局部特征而忽略了特征之間的上下文聯(lián)系，實驗結(jié)果顯示本文方法取得了更好的預(yù)測結(jié)果；DeepGS[8]利用GAT和BiGRU分別提取藥物指紋(fp)和SMILES序列特征，利用CNN提取靶點蛋白特征，有效增強(qiáng)了藥物特征表示卻忽視了靶點蛋白的特征表示，實驗結(jié)果顯示本文方法取得了更好的預(yù)測結(jié)果。

表7 KIBA數(shù)據(jù)集的結(jié)果對比

表8 Davis數(shù)據(jù)集的結(jié)果對比

4 結(jié)束語

本文為有效解決單模型提取特征種類受限問題，提出了一種深度并行全局特征提取策略。該策略利用CNN和特征存儲融合層構(gòu)建局部特征提取器，實現(xiàn)了局部特征有效地提取、存儲與壓縮，并緩解了數(shù)據(jù)集差異對特征提取效率的影響；利用CNN-BiLSTM串行混合模型構(gòu)建上下文特征提取器，提取局部特征之間的上下文聯(lián)系；最后，特征融合這兩種互補特征得到藥物靶點的全局特征向量。本文致力于設(shè)計全局特征提取策略，在利用詞嵌入特征向量化藥靶序列時只采用了基本的詞嵌入方式(Embedding)，后期會研究不同的詞嵌入生成方式對藥物靶點親和力預(yù)測的影響。