郭帥旗,閆效鶯
(西安石油大學 計算機學院,西安 710065)
近年來,卷積神經(jīng)網(wǎng)絡(luò)(Convolution Neural Network,CNN)在圖像處理、機器視覺和自然語言處理等領(lǐng)域均取得了較好的應用。但隨著大量非結(jié)構(gòu)化圖數(shù)據(jù),如交通網(wǎng)絡(luò)、社交網(wǎng)絡(luò)和生物信息學領(lǐng)域中蛋白質(zhì)作用關(guān)系網(wǎng)絡(luò)(PPI)等的出現(xiàn),圖卷積神經(jīng)網(wǎng)絡(luò)(Graph Convolution Network,GCN)應運而生,并迅速發(fā)展,已成為圖表示學習的重要方法之一。
隨著健康理念的發(fā)展,精準醫(yī)療已成為疾病診療的方向,如何在分子水平上精確測定個體病人對藥物治療的響應情況,是精準醫(yī)療的基礎(chǔ)和關(guān)鍵。然而,對個體病人進行大量藥物臨床試驗是不可行的,因此需要建立模型,預測藥物對個體疾病的敏感性響應關(guān)系。隨著高通量測序技術(shù)的發(fā)展,藥物敏感性數(shù)據(jù)庫,如NCI-60 數(shù)據(jù)[1]、癌癥細胞系百科全書(CCLE)[2]和癌癥藥物敏感性基因組學數(shù)據(jù)庫(GDSC)[3]等相繼發(fā)布,這些數(shù)據(jù)庫中整合收錄了大量細胞系與藥物之間的敏感性關(guān)系數(shù)據(jù)?;诖?,近些年國內(nèi)外學者提出了許多藥物-細胞系響應的預測方法。這些方法大致可分為基于回歸模型[4]、基于網(wǎng)絡(luò)推斷[5]、基于矩陣分解[6]和基于深度學習[7]的預測方法等等。如:Iorio 等人[4]通過彈性網(wǎng)絡(luò)和LASSO 構(gòu)建基因表達值與響應值之間的回歸預測模型,但其缺點是忽略了藥物的相關(guān)信息;Yang 等人[5]利用網(wǎng)絡(luò)表示學習提取特征,并采用SVM 預測藥物響應關(guān)系(NRL2DRP);Stanfield 等人[8]提出基于異構(gòu)網(wǎng)絡(luò)的帶重啟隨機游走算法預測藥物敏感性響應(HRWR);Yan 等人[6]提出一種具有可解釋性的三矩陣分解方法,預測藥物敏感性響應(TMF);Li 等人[9]提出一種應用堆疊的深度自動編碼器方法預測藥物敏感性響應(DeepDSC)。
雖然上述算法已極大推動了藥物-細胞系作用關(guān)系預測的研究,但預測精度仍有很大的提升空間。特別是近幾年圖神經(jīng)網(wǎng)絡(luò)的出現(xiàn),已成功應用于藥物-靶蛋白,藥物-藥物作用關(guān)系預測等生物信息學相關(guān)問題研究[10-11]之中??紤]藥物-細胞系敏感關(guān)系中網(wǎng)絡(luò)節(jié)點的異質(zhì)性,本文對適用于同構(gòu)網(wǎng)絡(luò)的GCN 算法進行改進,提出了一種新的基于異構(gòu)圖卷積網(wǎng)絡(luò)和深度神經(jīng)網(wǎng)絡(luò)的藥物-細胞系響應預測方法(HGCNDCP)。該方法首先分別計算藥物相似性和細胞系相似性,并融合藥物相似性特征、細胞系相似性特征以及藥物-細胞系響應關(guān)系,構(gòu)建異構(gòu)網(wǎng)絡(luò);然后在異構(gòu)網(wǎng)絡(luò)上使用圖卷積操作,通過不斷聚合鄰居節(jié)點特征,可同時捕獲異構(gòu)圖網(wǎng)絡(luò)的拓撲結(jié)構(gòu)特征和節(jié)點特征,得到藥物和細胞系兩類對象的特征表示數(shù)據(jù),最后使用深度神經(jīng)網(wǎng)絡(luò)(DNN)預測藥物-細胞系響應關(guān)系,并在GDSC 數(shù)據(jù)集中對算法進行驗證。
GCN 是CNN 算法在圖數(shù)據(jù)領(lǐng)域應用的產(chǎn)物。GCN 模型可用于捕獲非歐氏空間中存在的復雜網(wǎng)絡(luò)關(guān)系及對象(或?qū)嶓w)間的各種依賴關(guān)系。該模型通過不斷聚合鄰居節(jié)點信息,來更新自身節(jié)點特征,可同時獲捕圖結(jié)構(gòu)拓撲特征和節(jié)點特征。因此,GCN 可從原始圖數(shù)據(jù)和節(jié)點特征中,更好地進行特征表示與學習。基于同構(gòu)圖的GCN 模型定義如下:
給定圖G =(V,E),其中V是n個節(jié)點的集合,E是節(jié)點之間邊的集合;對應的鄰接矩陣記作A,其元素aij代表節(jié)點vi與vj之間的連接關(guān)系,節(jié)點特征矩陣記為X∈Rn×p;一 階GCN 模型定義為為歸一化的鄰接矩陣,Hl、Wl分別為l 層的節(jié)點表示和映射權(quán)重矩陣。
1.2.1 構(gòu)建藥物相似性矩陣
藥物特征描述符主要包括化學結(jié)構(gòu)描述符、分子指紋等。本文采用Pubchem 數(shù)據(jù)庫中記錄的藥物分子指紋描述符[6],將藥物表示為881 維的子結(jié)構(gòu)特征。因此,藥物特征矩陣可表示為Fd∈RN×p,p =881,N為藥物數(shù)目。其中,藥物di的特征記為=[si1,…,sil,…,sip]。若第l個子結(jié)構(gòu)特征在藥物di中存在sil =1,否則sil =0。由于具有相似化學結(jié)構(gòu)的藥物在細胞系中表現(xiàn)出相似的反應,在此使用Jaccard 系數(shù)計算藥物結(jié)構(gòu)相似性,其公式如下:
1.2.2 構(gòu)建細胞系相似性矩陣
對于細胞系的基因表達譜數(shù)據(jù)來說,每個細胞系均包含16 383個基因的表達值,因此細胞系的特征矩陣可表示為Fc∈RM×q、q =16 383,M為細胞系數(shù)目,其中細胞系ci的特征記為fci =[gi1,…,gil,…,giq],gil表示第l個基因在細胞系ci中的表達值。由于具有相似基因表達譜的細胞系會表現(xiàn)出相似的藥物反應。本文使用皮爾遜相關(guān)系數(shù)計算細胞系之間的相似性,公式如下:
1.2.3 藥物-細胞系響應關(guān)系網(wǎng)絡(luò)
已知的藥物-細胞系響應關(guān)系可以表示為二分圖G =(V,E),其中V ={Vc,Vd}表示藥物和細胞系兩類節(jié)點,E表示藥物與細胞系之間已知的IC50 響應值。本文使用Iorio 等人[12]提供的閾值,將已觀察響應值劃分為敏感性和耐藥性兩類。其中,敏感響應16 804 個,耐藥響應125 647 個,未知響應33 595個,由此構(gòu)建鄰接矩陣為Acd。
由藥物相似性矩陣、細胞系相似性矩陣和藥物-細胞系響應關(guān)系構(gòu)建形成的異構(gòu)網(wǎng)絡(luò)模型如圖1 所示。
圖1 藥物-細胞系響應異構(gòu)網(wǎng)絡(luò)Fig.1 Drug-cell line heterogeneous network model
基于異構(gòu)圖卷積的藥物-細胞系響應預測模型如圖2 所示,算法實現(xiàn)步驟如下:
圖2 HGCNDCP 網(wǎng)絡(luò)模型圖Fig.2 The flowchart of HGCNDCP pipeline
Input:藥物相似性網(wǎng)絡(luò)Sd、細胞系相似性網(wǎng)絡(luò)Sc、藥物-細胞系二分圖網(wǎng)絡(luò)Acd(邊權(quán)重“1”和“0”,分別表示敏感性和耐藥性響應類別);
Output:藥物-細胞系響應關(guān)系預測得分。
Step 1對Sd、Sc以及Acd按如下方式重構(gòu),得到異構(gòu)網(wǎng)絡(luò)鄰接矩陣A和特征矩陣S。
Step 2矩陣歸一化
Step 3基于異構(gòu)網(wǎng)絡(luò)的圖卷積操作,得到藥物和細胞系嵌入特征表示為F':
Step 4特征向量聚合,將藥物嵌入特征和細胞系嵌入特征拼接形成藥物-細胞系對的特征X∈RK×P,K為樣本數(shù),P為樣本特征的維度,見公式(6)。
Step 5構(gòu)建預測器。使用深度神經(jīng)網(wǎng)絡(luò)(DNN)作為HGCNDCP 的預測器。
其中,Zout和Zk(k =0,…,l)是DNN 模型中對應權(quán)重Wout、Wk和偏置bout、bk的隱層神經(jīng)元,Z0=X。y∈RK×t為K個藥物-細胞系樣本對的預測值。
本文使用開源的GDSC[3]作為基準數(shù)據(jù)集,網(wǎng)址為http:/ /www.cancerrxgene.org/,其中包括256 個藥物、1 001 個癌癥細胞系,以及藥物和細胞系的對數(shù)變換半抑制濃度值IC50[13]。該值代表要使50%的細胞生長受到抑制所需的藥物濃度,是藥物-細胞系響應的測量值。考慮到實驗的具體進行,需要對GDSC 基準數(shù)據(jù)集進行篩選、清洗等預處理。經(jīng)預處理后,本文得到183 種同時具有化學結(jié)構(gòu)特征和藥物反應數(shù)據(jù)的藥物,962 種同時具有基因組特征和藥物反應數(shù)據(jù)的細胞系。在這些藥物與細胞系之間,藥物-細胞系響應總共有176 046 個,其中敏感響應16 804 個,耐藥響應125 647,未知響應33 595個。
使用PyCharm 集成開發(fā)環(huán)境,Pytorch 1.7.1 作為框架。采用文獻[5]的驗證方法,基于上述預處理后的數(shù)據(jù)集,采用5-CV 交叉驗證方法進行實驗,即將數(shù)據(jù)隨機分成大致相等的5 份,每一份輪流作為測試樣本,其余4 份做訓練集。對測試集中每個藥物-細胞系對樣本進行預測,并將預測結(jié)果與實際標簽進行對比。使用ROC曲線下面積AUC表征模型預測性能。AUC值越大,表示算法性能越好。
本文算法包括基于異構(gòu)GCN 的特征表示和DNN 預測器兩部分,其中第一部分的損失函數(shù)采用二元加權(quán)交叉熵,第二部分采用二元交叉熵,見公式(10)、公式(11)。
公式(10)中,p(aij)為aij的真實標簽,q(aij)=是由異構(gòu)GCN 生成的兩類節(jié)點嵌入特征向量內(nèi)積計算出的預測概率,Wpos為負樣本與正樣本數(shù)目比的權(quán)重;公式(11)中,s(aij)是DNN預測的得分值,其值越大表示該樣本呈現(xiàn)敏感性的概率越大,T為類別數(shù)量,β為權(quán)值衰減系數(shù)。式子前一項旨在對所有類別計算平均損失,后一項旨在為權(quán)值矩陣和偏置矩陣提供L2 范數(shù)約束。兩部分均使用Adam 優(yōu)化器[14]、ReLU 激活函數(shù)和批量歸一化[15]處理。通過尋找損失函數(shù)的最小值和最佳精度對參數(shù)進行網(wǎng)格搜索。其中,基于異構(gòu)GCN 的特征表示部分,分別設(shè)置嵌入特征數(shù)ne ={5,25,50,75,100},并根據(jù)式(10)計算訓練損失。由圖3 可見,不同數(shù)量的嵌入特征,訓練過程相似,在200 輪訓練之前,損失快速下降,在1 000 輪之后,陸續(xù)趨于收斂狀態(tài)。其中潛在因子數(shù)ne為75和100 時,損失值相差無幾,誤差可縮小至10-3數(shù)量級。本文設(shè)定ne =100,學習率lr =0.01。DNN 預測器的各層維度分別為[200,128,96,64,2],lr =3.25e-5,衰減系數(shù)β =1e-5。
圖3 隱層節(jié)點數(shù)對異構(gòu)圖GCN 特征提取的影響Fig.3 Influence of the number of latent factors on the heterogeneous GCN feature extractor
采用5 -CV 交叉驗證方法,將本文算法HGCNDCP 與HNMDRP[16]、HRWR[8]、NRL2DRP[5]和TMF[6]算法進行比較。數(shù)據(jù)集中敏感性數(shù)據(jù)為正樣本,耐藥性為負樣本。ROC曲線下面積AUC結(jié)果詳見表1。
表1 算法性能比較結(jié)果Tab.1 Performance comparison
由表1 可見,HGCNDCP 的AUC值比HNMDRP提高了20.97%、比HRWR[8]算法提高了10.14%、比NRL2DRP 提高了16.06%、比TMF 提高了12.32%,證明HGCNDCP 具有更優(yōu)的預測性能。
為了評估不同k-CV 對模型性能的影響,本文分別進行2-CV、5-CV 和10-CV 交叉驗證,其對應AUC結(jié)果如圖4 所示。
圖4 不同交叉驗證方法下預測性能Fig.4 Prediction performance under different cross validation methods
結(jié)果表明,預測精度隨著訓練數(shù)據(jù)集的增多而增加,10-CV 驗證的訓練數(shù)據(jù)集大于2-CV 驗證和5-CV 驗證,其AUC值也高于兩者。
為了更高效地預測藥物-細胞系之間的敏感性響應關(guān)系,本文在圖卷積神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上提出了基于異構(gòu)圖卷積網(wǎng)絡(luò)的藥物-細胞系響應預測方法(HGCNDCP)。研究表明:
(1)使用藥物結(jié)構(gòu)分子指紋特征數(shù)據(jù)、細胞系的基因表達譜數(shù)據(jù)和藥物-細胞系作用關(guān)系數(shù)據(jù),對學習藥物、細胞系的特征提取有重要影響。
(2)通過構(gòu)建異構(gòu)網(wǎng)絡(luò),使用圖卷積神經(jīng)網(wǎng)絡(luò),能夠有效地聚合鄰居特征信息,得到較好的藥物和細胞系的表征。
(3)通過使用GDSC 數(shù)據(jù)集,并與其它算法的一系列實驗比對,HGCNDCP 具有較高的預測精度,能夠較好地預測藥物細胞系響應,從而為藥物細胞系響應預測提供有效的思路和方法。