徐國保，陳媛曉，王 驥

（廣東海洋大學(xué)電子與信息工程學(xué)院，廣東湛江 524088）

（*通信作者電子郵箱xuguobao@126.com；zjouwangji@163.com）

0 引言

藥物與人類疾病息息相關(guān)，而藥物靶標(biāo)點(diǎn)的確認(rèn)是藥物研發(fā)工作的開始，因此，快速準(zhǔn)確的預(yù)測藥物-靶標(biāo)的相互作用是藥物研發(fā)的關(guān)鍵。然而，受到成本、通量等的影響，傳統(tǒng)的用于闡明藥物-靶標(biāo)關(guān)系的生物實(shí)驗(yàn)難以展開，很多潛在的藥物靶標(biāo)相互作用關(guān)系尚未被研發(fā)出來。傳統(tǒng)的計(jì)算方法主要有分子對(duì)接方法［1］和基于配體的方法［2］。然而，當(dāng)靶標(biāo)蛋白的三維結(jié)構(gòu)不明確時(shí)，分子對(duì)接方法的預(yù)測性能下降；當(dāng)只有少數(shù)已知配體與靶標(biāo)結(jié)合時(shí)，基于配體的方法預(yù)測結(jié)果往往較差。

過去十年，國內(nèi)外學(xué)者致力于利用機(jī)器學(xué)習(xí)方法來預(yù)測藥物-靶標(biāo)相互作用關(guān)系，這些機(jī)器學(xué)習(xí)方法主要可分為有監(jiān)督學(xué)習(xí)和半監(jiān)督學(xué)習(xí)兩大類：有監(jiān)督學(xué)習(xí)如二分圖局部模型方法（Bipartite Local Model，BLM）［3］、基于核的回歸模型方法［4］及基于分子特征方法［5］等；半監(jiān)督學(xué)習(xí)如基于拉普拉斯正則化的最小二乘法NetLapRLS［6］。文獻(xiàn)［6］中提出的NetLapRLS通過整合已知的化學(xué)結(jié)構(gòu)信息、基因序列數(shù)據(jù)和藥物-蛋白相互作用網(wǎng)絡(luò)對(duì)藥物-蛋白質(zhì)相互作用關(guān)系進(jìn)行了預(yù)測［7］。文獻(xiàn)［8］使用帶高斯相互作用屬性（Gaussian Interaction Profile，GIP）核的正則化最小二乘法（Regularized Least Squares，RLS）分類器，結(jié)合局部二分圖模型（Bipartite Local Model，BLM），提出基于鄰居相互作用譜的局部二分圖模型（Bipartite Local Model-Neighbor-based Interaction-profile Inferring，BLM-NII），此外，有監(jiān)督學(xué)習(xí)的方法還有二分圖模型［9］、基于網(wǎng)絡(luò)的推斷模型（Network-based Inference，NBI）［10］等。文獻(xiàn)［11］提出異質(zhì)網(wǎng)絡(luò)上的可重啟隨機(jī)游走（Networkbased Random Walk with Restart on the Heterogeneous network，NRWRH）方法，該方法結(jié)合已知的藥物-靶標(biāo)關(guān)聯(lián)關(guān)系、藥物相似網(wǎng)絡(luò)及靶標(biāo)相似網(wǎng)絡(luò)，構(gòu)建異質(zhì)網(wǎng)絡(luò)并在該網(wǎng)絡(luò)上執(zhí)行隨機(jī)游走算法。上述兩種方法都沒有考慮到靶標(biāo)信息未知的藥物。此外，文獻(xiàn)［12］中結(jié)合標(biāo)記數(shù)據(jù)（已知和未知）提出一種具有網(wǎng)絡(luò)一致性的藥物靶標(biāo)相互作用半監(jiān)督預(yù)測方法，但是該方法嚴(yán)重依賴于藥物和靶標(biāo)的相似性。

以上這些方法雖然提高了機(jī)器學(xué)習(xí)性能，但它們僅使用了少量的生物特征信息，在預(yù)測潛在的藥物靶標(biāo)關(guān)系中可能丟失一些重要的信息，如藥物靶標(biāo)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)信息等。

針對(duì)現(xiàn)有技術(shù)的不足，在前人工作的基礎(chǔ)上，本文利用圖卷積網(wǎng)絡(luò)（Graph Convolutional Network，GCN）結(jié)合自編碼技術(shù)來預(yù)測潛在的藥物-靶標(biāo)關(guān)系對(duì)。在考慮藥物與靶標(biāo)多種生物信息與拓?fù)渚W(wǎng)絡(luò)結(jié)構(gòu)的前提下，綜合運(yùn)用深度學(xué)習(xí)知識(shí)，設(shè)計(jì)合理的預(yù)測方案，提高藥物-靶標(biāo)相互作用預(yù)測準(zhǔn)確度。實(shí)驗(yàn)結(jié)果表明，該方法能夠有效預(yù)測藥物-靶標(biāo)相互作用關(guān)系，且具有較強(qiáng)的魯棒性。

1 基于圖卷積的藥物靶標(biāo)關(guān)聯(lián)預(yù)測算法

1.1 數(shù)據(jù)集構(gòu)建

本文的數(shù)據(jù)集來自文獻(xiàn)［13］，可以在https：//github.com/luoyunan/DTINet 中下載，該數(shù)據(jù)集包含1 923 個(gè)已知的藥物-靶標(biāo)對(duì)，1 512 種不同類型的靶標(biāo)蛋白質(zhì)和708 種不同類型的藥物化合物。本文的任務(wù)是要利用已知的藥物-靶標(biāo)關(guān)聯(lián)數(shù)據(jù)及藥物和靶標(biāo)特征數(shù)據(jù)，從未標(biāo)記的藥物靶標(biāo)中尋找潛在的藥物靶標(biāo)對(duì)。在現(xiàn)有的數(shù)據(jù)集中只有藥物相似矩陣和靶標(biāo)相似矩陣，若直接使用稠密相似矩陣作為圖，計(jì)算是非常耗時(shí)的。此外，稠密圖還會(huì)產(chǎn)生噪聲，影響模型性能，因此，需要利用已知數(shù)據(jù)進(jìn)行特征提取。本文特征提取的方法參照文獻(xiàn)［14］。

為了得到稀疏圖來避免耗時(shí)的計(jì)算，藥物G(u) ∈Rd×d和靶標(biāo)G(v) ∈Rt×t之間的相似性可表示為：

其中：Sim(i，j)為藥物相似矩陣和靶標(biāo)相似矩陣，h(x)是x的鄰居集合。本文分別取前10 個(gè)藥物相似矩陣鄰居節(jié)點(diǎn)和前50 個(gè)靶標(biāo)相似矩陣鄰居節(jié)點(diǎn)，得到G(u)和G(v)后，將其作為藥物和靶標(biāo)的特征輸入。

1.2 圖卷積

將已知藥物與靶標(biāo)之間的關(guān)聯(lián)關(guān)系表示成二分圖，則本文提到的預(yù)測任務(wù)可以定義為在這樣一個(gè)圖上進(jìn)行半監(jiān)督預(yù)測。

給定一個(gè)二分圖G={V，E}，其中V=(vd，vt)表示nd個(gè)藥物節(jié)點(diǎn)和nt個(gè)靶標(biāo)節(jié)點(diǎn)，Xd=分別表示藥物特征矩陣和靶標(biāo)特征矩陣。由于藥物和靶標(biāo)節(jié)點(diǎn)的特征維度都很高，傳統(tǒng)的相似性度量方法如歐幾里得距離無法取得很好的效果。為此，本文使用譜圖卷積有效利用圖拓?fù)浜凸?jié)點(diǎn)特征信息。

目前在圖數(shù)據(jù)上使用卷積濾波器的方法大致可分為兩類：空域圖卷積［15-16］和譜圖卷積［17-19］?？沼驁D卷積本質(zhì)是不斷聚合節(jié)點(diǎn)的鄰居信息，即直接將卷積操作定義在每個(gè)節(jié)點(diǎn)的連接關(guān)系上，文獻(xiàn)［20］中曾指出這種方法存在的問題。相對(duì)于空間卷積，譜圖卷積則將卷積網(wǎng)絡(luò)濾波器與圖信號(hào)同時(shí)變換到傅立葉域后進(jìn)行處理。

譜圖卷積可以定義為濾波器gθ=diag(θ)（θ∈RN）與信號(hào)x∈RN在傅里葉域的乘積：

其中UTx表示x的圖傅里葉變換，在這里gθ可以看作L特征向量的函數(shù)，即gθ(Λ)，Λ是特征值對(duì)角矩陣。

當(dāng)圖中節(jié)點(diǎn)數(shù)量多、節(jié)點(diǎn)關(guān)系復(fù)雜時(shí)，拉普拉斯矩陣L進(jìn)行特征分解需要很大的計(jì)算量，為了解決這個(gè)問題，采用切比雪夫多項(xiàng)式Tk(x)直到第k階的截?cái)嗾归_來近似gθ(Λ)。

文獻(xiàn)［17］通過限制k=1 并將L的最大特征值近似為2，進(jìn)一步簡化了譜圖卷積的定義：

1.3 模型概述

如前所述，藥物和靶標(biāo)之間的關(guān)聯(lián)預(yù)測問題可以當(dāng)作一個(gè)半監(jiān)督預(yù)測問題。很多基于圖卷積神經(jīng)網(wǎng)絡(luò)（GCN）的方法主要用于解決同質(zhì)網(wǎng)絡(luò)上的節(jié)點(diǎn)分類問題，為充分利用圖卷積，使其能夠解決異質(zhì)、二部、有屬性網(wǎng)絡(luò)的預(yù)測問題，文獻(xiàn)［21］首次結(jié)合圖卷積與自編碼技術(shù)，提出基于圖卷積的MicroRNA 和抗藥性關(guān)聯(lián)預(yù)測（Graph Convolution for association between MicroRNA and Drug Resistance，GCMDR）算法。為了確保模型有效訓(xùn)練，本文在文獻(xiàn)［21］的基礎(chǔ)上，引入了集成學(xué)習(xí)（Ensemble Learning，EL）中的堆疊思想，將兩個(gè)組件線性組合在一起，聯(lián)合訓(xùn)練。

給定鄰接矩陣M∈，其中Nd為藥物節(jié)點(diǎn)數(shù)量，Nt為蛋白質(zhì)節(jié)點(diǎn)數(shù)量，矩陣的值Mij表示藥物i與蛋白質(zhì)j是否已通過生物實(shí)驗(yàn)驗(yàn)證存在關(guān)聯(lián)：

模型的目標(biāo)是通過構(gòu)建基于圖卷積的編碼器[Fd，F(xiàn)t]=fen(v，ε，Xd，Xt)來學(xué)習(xí)藥物和靶標(biāo)的嵌入特征F，并且通過構(gòu)建解碼器M′=fde(Fd，F(xiàn)t)來預(yù)測新鏈接，式中Xd∈分別表示藥物和靶標(biāo)的輸入特征矩陣，分別表示藥物和靶標(biāo)學(xué)習(xí)到的特征矩陣。

為此，本文提出的模型由兩種不同類型的層組成：1）用于在藥物和靶標(biāo)相互作用網(wǎng)絡(luò)圖上整合其節(jié)點(diǎn)特征的編碼層；2）使用上一層學(xué)到的嵌入特征來預(yù)測全連接交互網(wǎng)絡(luò)的編碼層。算法結(jié)構(gòu)如圖1。

圖1 本文算法結(jié)構(gòu)Fig.1 Structure of the proposed algorithm

1.3.1 編碼層

編碼層的輸入包括藥物和靶標(biāo)的原始特征矩陣Xd、Xt和鄰接矩陣M。為了把藥物特征矩陣和靶標(biāo)特征矩陣整合成一個(gè)輸入矩陣，定義一個(gè)新的特征矩陣：

同樣，鄰接矩陣重新定義為：

然后，GCDTI 對(duì)矩陣X進(jìn)行行歸一化：Xrw=D-1X，其中表示輸入信號(hào)矩陣。根據(jù)式（6），可得到一個(gè)圖卷積矩陣G：

通過引入G的權(quán)重矩陣We和偏重矩陣Be構(gòu)建隱藏層，選擇ReLU函數(shù)作為激活函數(shù)，則編碼層的輸出F如下：

式中可訓(xùn)練權(quán)重矩陣We∈為傅里葉系數(shù)矩陣，它將矩陣G轉(zhuǎn)化為描述藥物靶標(biāo)節(jié)點(diǎn)與潛在因子之間的關(guān)聯(lián)的隱藏矩陣F，Ne表示潛在因子的數(shù)量，是手動(dòng)設(shè)置的。編碼層的輸出是輸入特征到隱藏空間的投影，其由兩部分組成，分別是藥物的嵌入特征矩陣Fd和靶標(biāo)的嵌入特征矩陣Ft。

1.3.2 解碼層

為了重構(gòu)藥物-靶標(biāo)關(guān)聯(lián)矩陣，構(gòu)建解碼器M′=fde(Fd，F(xiàn)t)如下：

式中權(quán)重矩陣Wd∈RL×L描述了隱藏層潛在因子的相似性。在本文模型中使用文獻(xiàn)［22］提出的初始化方法來隨機(jī)初始化矩陣We、Be、Wd。

顯然，輸出矩陣M′與輸入矩陣M的維度相同，M′中的值表示所有藥物-靶標(biāo)對(duì)的權(quán)值。所有在矩陣M中值為0 的藥物-靶標(biāo)對(duì)將由解碼器賦予一個(gè)預(yù)測值，預(yù)測得分高的藥物-靶標(biāo)對(duì)更有可能是關(guān)聯(lián)的。

此外，為了在半監(jiān)督學(xué)習(xí)下訓(xùn)練模型，本文使用了負(fù)抽樣方法。在每個(gè)訓(xùn)練階段，隨機(jī)選擇未標(biāo)記的藥物-靶標(biāo)對(duì)作為負(fù)樣本進(jìn)行訓(xùn)練。給定訓(xùn)練集，模型嘗試最小化以下?lián)p失函數(shù)：

1.4 交叉驗(yàn)證

為了評(píng)估模型性能，使用了k折交叉驗(yàn)證。將數(shù)據(jù)隨機(jī)分成k份，每一份輪流作測試樣本（假設(shè)關(guān)聯(lián)關(guān)系未知）。在性能評(píng)估中，將測試樣本和所有未標(biāo)記的藥物-靶標(biāo)對(duì)都視為候選樣本，若測試樣本在所有候選樣本中排名較前，則表明該模型具有良好的預(yù)測性能。

在本文中，k分別取2、5、10，由于數(shù)據(jù)集中的樣本數(shù)量有限，因此對(duì)數(shù)據(jù)集進(jìn)行10 次劃分，最后取10 次實(shí)驗(yàn)的平均值作為模型整體的性能指標(biāo)。

2 實(shí)驗(yàn)結(jié)果與分析

2.1 算法運(yùn)行環(huán)境

使用PyCharm 集成開發(fā)環(huán)境，TensorFlow 2.0.0 作為框架。

2.2 圖卷積有效性評(píng)估

在本節(jié)中就圖卷積整合原始輸入特征數(shù)據(jù)的能力進(jìn)行了評(píng)估，具體來說，將模型輸入原始特征數(shù)據(jù)與刪除輸入特征的情況進(jìn)行對(duì)比。為此，將式（9）中輸入特征矩陣A的每個(gè)值都替換為1，在這種情況下，由于所有節(jié)點(diǎn)的特征都相同，因此圖卷積就沒有意義了。本節(jié)使用5 折交叉驗(yàn)證法來做有無特征輸入的對(duì)比實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果表明，在沒有任何特征輸入的情況下，模型的受試者工作特性（Receiver Operating Characteristic，ROC）曲線下的面積（Area Under ROC Curve，AUC）為0.888 9，而有特征輸入時(shí)，模型的AUC 為0.920 1，明顯高于沒有特征輸入的情況。實(shí)驗(yàn)結(jié)果證明提出的基于圖卷積的模型能夠有效地整合輸入特征數(shù)據(jù)。

2.3 k折交叉驗(yàn)證的預(yù)測性能

為了評(píng)估本文提出的藥物靶標(biāo)關(guān)聯(lián)預(yù)測模型的性能，分別使用了2 折、5 折和10 折交叉驗(yàn)證，使用不同交叉驗(yàn)證的AUC平均值見表1。

表1 不同交叉驗(yàn)證方法的預(yù)測性能Tab.1 Prediction performance of different cross validation methods

從表中可以看出，預(yù)測精度隨著訓(xùn)練數(shù)據(jù)集的增加而增加，由于10 折交叉驗(yàn)證的訓(xùn)練數(shù)據(jù)集大于2 折交叉驗(yàn)證和5折交叉驗(yàn)證，因此其平均AUC最高，為0.924 6±0.004 8。

此外，圖2和圖3給出了不同交叉驗(yàn)證方法下的訓(xùn)練損失和訓(xùn)練誤差，訓(xùn)練損失和訓(xùn)練誤差分別由式（12）和式（12）的第一項(xiàng)計(jì)算而得。從圖中的曲線可以看出，不同交叉驗(yàn)證方法下GCDTI 的訓(xùn)練過程是相似的。在大多數(shù)實(shí)驗(yàn)中，訓(xùn)練損失和訓(xùn)練誤差可以分別在第300 次遍歷數(shù)據(jù)集和第250 次遍歷數(shù)據(jù)集之前收斂到下界，說明采用不同的交叉驗(yàn)證方法時(shí)，數(shù)據(jù)的差異對(duì)計(jì)算過程的影響很小，說明模型具有很強(qiáng)的魯棒性。

圖2 不同交叉驗(yàn)證方法下的訓(xùn)練損失Fig.2 Training loss under different cross validation methods

圖3 不同交叉驗(yàn)證方法下的訓(xùn)練誤差Fig.3 Training error under different cross validation methods

2.4 負(fù)樣本對(duì)模型性能的影響

由于在本文所收集的數(shù)據(jù)庫中只有正樣本，因此需要找到負(fù)樣本進(jìn)行半監(jiān)督訓(xùn)練來提高模型的預(yù)測性能。為此，對(duì)未標(biāo)記的藥物-靶標(biāo)對(duì)進(jìn)行采樣，以生成負(fù)樣本進(jìn)行訓(xùn)練。然而，負(fù)樣本的數(shù)量也會(huì)對(duì)模型的預(yù)測性能產(chǎn)生影響，大量的負(fù)樣本可以為訓(xùn)練提供數(shù)據(jù)資源，提高模型性能，但這同時(shí)也可能會(huì)造成訓(xùn)練數(shù)據(jù)不平衡問題。因此，負(fù)樣本數(shù)量的選擇對(duì)準(zhǔn)確預(yù)測GCDTI 模型是非常重要的。在每次采樣中，負(fù)樣本集固定為正樣本集的p倍。圖4 展示了不同負(fù)樣本集對(duì)預(yù)測性能的影響，從圖中可以看出，當(dāng)負(fù)樣本數(shù)是正樣本數(shù)的10倍時(shí)，該模型的預(yù)測性能最高，AUC 為0.920 2±0.011 1。當(dāng)p設(shè)為0 時(shí)，意味著沒有使用負(fù)樣本，只用正樣本進(jìn)行訓(xùn)練，此時(shí)AUC 遠(yuǎn)低于有負(fù)樣本的情況。AUC 從p=0 到p=10 的變化體現(xiàn)了負(fù)樣本對(duì)GCDTI的重要性和有效性。

圖4 不同負(fù)樣本數(shù)量時(shí)的預(yù)測性能Fig.4 Prediction performance when having different negative sample numbers

2.5 不同潛在因子數(shù)的比較

由于本文提出的模型是基于潛在因子模型構(gòu)建的，因此隱藏層的大小對(duì)其性能預(yù)測至關(guān)重要。本節(jié)討論潛在因子數(shù)量L對(duì)模型性能的影響，本節(jié)的實(shí)驗(yàn)結(jié)果基于10 折交叉驗(yàn)證法并且負(fù)樣本是正樣本的10 倍。從圖5 可看出，當(dāng)L在5～80取值時(shí)，AUC 的平均值呈單峰分布，且當(dāng)L=25時(shí)，模型性能達(dá)到最優(yōu)。25 這個(gè)數(shù)字可能反映了藥物與靶標(biāo)之間真實(shí)關(guān)聯(lián)的數(shù)量。

圖5 潛在因子數(shù)對(duì)模型的影響Fig.5 Influence of the number of latent factors on the model

2.6 與其他方法的比較

為了進(jìn)一步評(píng)估模型的預(yù)測性能，將實(shí)驗(yàn)結(jié)果與其他5種較為先進(jìn)的藥物靶標(biāo)關(guān)聯(lián)預(yù)測方法在同一數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果進(jìn)行比較，這5 種方法包括BLM-NII［8］、NetLapRLS［23］、異構(gòu)網(wǎng)絡(luò)模型（Heterogeneous Network Model，HNM）［24］、多相似度矩陣分解（multiple similarities Collaborative Matrix Factorization，CMF）模型［25］、藥物-靶標(biāo)相互作用預(yù)測的網(wǎng)絡(luò)集成方法（Network integration approach for Drug-Target Interaction prediction，DTINet）［13］。比較結(jié)果見圖6，從圖6可以看出，GCDTI的平均AUC最高，為0.9246±0.004 8，比DTINet 高1.13 個(gè)百分點(diǎn)。本節(jié)實(shí)驗(yàn)結(jié)果基于10 折交叉驗(yàn)證并且負(fù)樣本是正樣本的10 倍。這些實(shí)驗(yàn)表明利用端到端學(xué)習(xí)的模型架構(gòu)，當(dāng)需要預(yù)測大量的藥物和靶標(biāo)數(shù)據(jù)的關(guān)聯(lián)關(guān)系時(shí)，GCDTI有潛力成為一種可靠的預(yù)測方法。

圖6 不同預(yù)測方法的性能比較Fig.6 Performance comparison of different prediction methods

3 結(jié)語

為了更高效地識(shí)別潛在的藥物-靶標(biāo)關(guān)系對(duì)，本文利用圖卷積神經(jīng)網(wǎng)絡(luò)結(jié)合自編碼技術(shù)，提出GCDTI 模型。該方法通過輸入已知的藥物-靶標(biāo)關(guān)系對(duì)以及藥物和節(jié)點(diǎn)的特征信息，以端到端學(xué)習(xí)的方式提取藥物和靶標(biāo)的嵌入特征。一系列實(shí)驗(yàn)結(jié)果表明，從模型中學(xué)習(xí)到的低維嵌入特征能夠有效地表達(dá)藥物靶標(biāo)之間的相互作用關(guān)系。但是在當(dāng)前的藥物靶標(biāo)相互作用模型中，圖卷積只能輸入數(shù)值型數(shù)據(jù)，因此其他非數(shù)值型的特征仍不適用于當(dāng)前的模型，接下來的工作將進(jìn)一步研究解決這個(gè)問題的方案。