高沖，秦玉芳，陳明

1.上海海洋大學(xué)信息學(xué)院，上海201306；2.農(nóng)業(yè)農(nóng)村部漁業(yè)信息重點實驗室，上海201306

前言

惡性腫瘤（癌癥）是嚴(yán)重影響人類健康的疾病之一［1］。雖然傳統(tǒng)治療癌癥的方法（放療、化療）有明顯的治療效果，但是大量研究表明腫瘤具有異質(zhì)性［2］，患有相同癌癥的病人使用相同的治療方法卻有不同的療效?；诖?，個性化醫(yī)療應(yīng)運而生，它關(guān)注每一位患者的特異性特征，其中測量患者對藥物的反應(yīng)是一個關(guān)鍵問題［3-4］。

隨著高通量基因組學(xué)技術(shù)的發(fā)展，藥物基因組學(xué)成為測量患者對藥物反應(yīng)的一個重要方法［5］。研究者通常通過基因或蛋白質(zhì)表達譜等分子圖譜來測量細胞對藥物的反應(yīng)，進而建立相應(yīng)的計算模型預(yù)測藥物反應(yīng)［6］。Gillet 等［7］發(fā)現(xiàn)在細胞系模型和臨床具有相關(guān)性的前提下，這些計算模型能識別決定藥物反應(yīng)的分子因素，并對患者群體進行相應(yīng)的個性化藥物治療。許多研究機構(gòu)開發(fā)了諸如癌癥細胞系百科全書（Cancer Cell Line Encyclopedia, CCLE）和腫瘤藥物敏感性基因組學(xué)（Genomics of Drug Sensitivity in Cancer,GDSC）等包含基因表達數(shù)據(jù)和拷貝數(shù)變異等基因組學(xué)數(shù)據(jù)以及藥物反應(yīng)值在內(nèi)的大型數(shù)據(jù)庫，這些大型數(shù)據(jù)集為識別新的藥物靶點和藥物反應(yīng)標(biāo)記物提供了更多的可能性［8］；同時，這也為開發(fā)藥物反應(yīng)計算模型提供了依據(jù)，如Papillon-Cavanagh 等［9］利用CCLE 和癌癥基因組計劃（Cancer Genome Project,CGP）數(shù)據(jù)集建立預(yù)測藥物反應(yīng)的線性模型，發(fā)現(xiàn)基因組預(yù)測因子能夠驗證對特定藥物的反應(yīng)；Masica等［10］利用CCLE數(shù)據(jù)集構(gòu)建多變量組合改變組織（MOCA）模型來識別藥物反應(yīng)的組合生物標(biāo)志物；Menden 等［11］利用GDSC 數(shù)據(jù)集和機器學(xué)習(xí)算法建立基于細胞系的基因組特征和藥物化學(xué)特性的藥物敏感性預(yù)測模型，并通過對實驗結(jié)果和已有事實的對比驗證該模型的有效性。

近年來，許多研究者根據(jù)基因水平特征建立抗癌藥物敏感性預(yù)測模型［6］。如Costello等［12］把基因表達譜或拷貝數(shù)變異等基因組學(xué)數(shù)據(jù)用于預(yù)測抗癌藥物反應(yīng)，發(fā)現(xiàn)基于基因表達數(shù)據(jù)建立的抗癌藥物敏感性預(yù)測模型具有很好的預(yù)測性能；Geeleher等［13］采用嶺回歸算法建立抗癌藥物反應(yīng)預(yù)測模型，同時使用獨立數(shù)據(jù)集驗證了該模型的有效性。這些方法大多基于基因表達數(shù)據(jù)等基因水平特征，在獨立研究中的重復(fù)性有限，這對生物學(xué)解釋提出了挑戰(zhàn)［14］。有研究表明考慮基因間相互作用行為比僅僅關(guān)注單個基因行為在預(yù)測藥物反應(yīng)上具有更好的預(yù)測效果［15］。通路數(shù)據(jù)庫是系統(tǒng)分析基因功能，聯(lián)系基因組信息和功能信息的知識庫。通路作為基因功能集合能夠提高預(yù)測模型的預(yù)測能力和解釋能力［16］。Wang 等［17］把通路數(shù)據(jù)和基因表達譜應(yīng)用到藥物敏感性預(yù)測，研究表明在CCLE 數(shù)據(jù)集的24 種藥物中，基于通路的模型較基于基因的模型具有更好的預(yù)測性能，并且基于通路的模型能識別更多藥物相關(guān)的基因或通路，具有更好的生物學(xué)解釋；然而該方法僅僅把通路作為基因集合，沒有考慮通路中基因互相作用關(guān)系。

針對以上問題，本研究提出一種整合通路網(wǎng)絡(luò)中高連接度基因和基因表達數(shù)據(jù)推斷通路活性，建立抗癌藥物敏感性預(yù)測模型，簡記為PHG（Pathway Hub Gene）。首先利用通路數(shù)據(jù)和STRING數(shù)據(jù)庫得到每個通路的基因相互作用網(wǎng)絡(luò)表，從該網(wǎng)絡(luò)表中選擇高連接度基因；然后分別計算每一個通路的活性向量；最后合并所有通路的活性向量，得到通路活性特征矩陣，以此作為抗癌藥物敏感性預(yù)測模型的輸入。10折交叉驗證的實驗結(jié)果表明，在17-AAG等大多數(shù)抗癌藥物上，并不是通路中所有基因都對藥物敏感性預(yù)測有幫助，考慮通路中的關(guān)鍵基因較通路全部基因構(gòu)建預(yù)測模型具有更好的預(yù)測效果，同時驗證了基于通路的模型較基于基因的模型能給出更好的生物學(xué)解釋。

1 數(shù)據(jù)集及預(yù)處理

本研究的基因表達和藥物IC50 值數(shù)據(jù)來自于CCLE 數(shù)據(jù)庫，下載地址為 https://portals.broadinstitute.org/ccle/data；同時為了獨立檢驗，也下載了GDSC 數(shù)據(jù)庫中的基因表達和藥物IC50 數(shù)據(jù)，下載地址為https://www.cancerrxgene.org/。為消除實驗技術(shù)和實驗平臺所導(dǎo)致的基因表達量誤差，采用以基因為中心的RMA 標(biāo)準(zhǔn)化算法對基因表達譜進行標(biāo)準(zhǔn)化處理。經(jīng)過標(biāo)準(zhǔn)化后，CCLE 基因表達譜共有18 900 個基因和1 036 個細胞系樣本，GDSC 基因表達譜中共有9 920個基因和697個細胞系樣本。

本研究使用IC50 值衡量藥物敏感性，類似于Wang 等［17］的做法，對藥物反應(yīng)IC50 值做log 變換。由于基因表達數(shù)據(jù)中的一些細胞系樣本在藥物反應(yīng)數(shù)據(jù)里不存在，所以本研究選取在基因表達數(shù)據(jù)和藥物反應(yīng)數(shù)據(jù)中同時存在的細胞系進行分析。例如，對于藥物AEW541 來說，NCIH2196_LUNG 細胞系存在基因表達譜中，但在AEW541 藥物反應(yīng)數(shù)據(jù)中沒有該細胞系，所以在做AEW541 藥物的敏感性預(yù)測時需去除該細胞系。

本研究使用的通路數(shù)據(jù)來自京都基因和基因組數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes,KEGG）中的通路數(shù)據(jù)庫，在該數(shù)據(jù)庫中下載每個通路的基因集，最終的通路數(shù)據(jù)集包括389 個通路，共有14 097 個基因。通路中基因間相互作用關(guān)系表可從STRING 數(shù)據(jù)庫中獲得，下載地址為https://www.string-db.org/cgi/download，STRING數(shù)據(jù)庫包含5 090個物種、24 584 628 種蛋白和3 123 056 667 個相互作用關(guān)系［18］，本研究下載的數(shù)據(jù)來自數(shù)據(jù)庫最新版本（Version 11.0）。

2 模型

2.1 基于通路高連接基因的通路活性推斷

為推斷通路活性，本研究不僅考慮通路中每個基因的表達水平，還關(guān)注了通路中基因間相互作用關(guān)系，基因相互作用關(guān)系在預(yù)測藥物敏感性具有更好的魯棒性［15］。首先從STRING 數(shù)據(jù)庫中得到每個通路（基因集）中基因間互相作用網(wǎng)絡(luò)表，表中的（Gi,Gj）表示基因Gi和基因Gj在通路中是相互連接的；接著根據(jù)通路互相作用網(wǎng)絡(luò)表計算通路活性向量。

計算每個基因在通路互相作用網(wǎng)絡(luò)表中的度，由于通路網(wǎng)絡(luò)中高連接度的Hub 基因?qū)φ麄€通路的功能起著更關(guān)鍵的作用［19］，所以從該網(wǎng)絡(luò)表中選擇高連接度Hub 基因來進行分析。將基因的度降序排序，選擇排名在前10%的基因作為Hub 基因，圖1中的Gh1,Gh2,…,Ghk為通路p1的Hub 基因；計算Hub基因表達值的平均值作為活性值，然后合并每個細胞系樣本中的活性值得到該通路的活性向量?；钚灾涤嬎愎饺缦拢?/p>

其中,hk表示通路P1中Hub基因的數(shù)量；gij表示基因i在細胞系樣本j中表達值；p1j表示細胞系樣本j在通路p1中的活性值。

按照上面的方法，計算所有通路的活性向量，合并得到通路活性矩陣（列為細胞系樣本，行為通路）。假定有l(wèi) 個通路，分別記為P1，P2，…,Pl，按照上述方法計算所有通路的活性向量后得到通路活性矩陣(pij)，其中i為通路，j為細胞系樣本。

總的來說，可將基因表達譜和通路中關(guān)鍵基因信息分析整合得到通路活性矩陣，以此來預(yù)測癌癥藥物敏感性?；谕分懈哌B接度基因模型的流程如圖1所示。

圖1 利用通路中高連接度基因表達推斷通路活性Fig.1 Using the expression of genes with high connectivity in the pathway to infer pathway activity

2.2 建立藥物反應(yīng)預(yù)測模型

將得到的通路活性特征矩陣作為預(yù)測模型的輸入，藥物敏感性水平作為模型輸出，根據(jù)均方誤差（Mean Square Error, MSE）來調(diào)試優(yōu)化模型的參數(shù)，并進行訓(xùn)練與預(yù)測。本研究采用機器學(xué)習(xí)中的彈性網(wǎng)作為預(yù)測算法。

彈性網(wǎng)是一種使用L1 和L2 范數(shù)作為先驗正則項訓(xùn)練的線性回歸模型［20］。這種組合可以學(xué)習(xí)到類似于Lasso 的一個稀疏模型，同時還保留嶺回歸的正則化屬性，既能實現(xiàn)重要特征變量的選擇，又能處理強相關(guān)性特征數(shù)據(jù)，具有較好的群組效應(yīng)，結(jié)合了嶺回歸和Lasso 回歸的優(yōu)點。因此，彈性網(wǎng)尤其適用于有多個特征彼此相關(guān)的場合。在基于通路/基因的預(yù)測模型中，作為特征的通路/基因相互之間實際上都是有聯(lián)系的。因此，本研究選用彈性網(wǎng)回歸算法來構(gòu)建預(yù)測模型，并使用R語言中g(shù)lmnet包實現(xiàn)彈性網(wǎng)算法。調(diào)整和優(yōu)化模型主要通過網(wǎng)格搜索，在1 000個參數(shù)中尋找最優(yōu)參數(shù)，其中α:[0.1,1 ]設(shè)置10 個參數(shù)，λ:[exp-5,exp5]設(shè)置100 個參數(shù)，使用10 折交叉驗證選取最優(yōu)參數(shù)。

3 結(jié)果與討論

3.1 通路之間的重疊性

本研究使用Jaccard 指數(shù)來評價兩個通路之間的重疊性。通過對通路間重疊性的研究，分析通路是否具有特異性，是否對實驗產(chǎn)生較大的誤差。Jaccard指數(shù)計算公式如下：

其中，P1∩P2表示同時存在于通路P1和通路P2的基因；P1∪P2表示存在于通路P1或P2的基因。由式（1）可以發(fā)現(xiàn)，當(dāng)兩個通路完全不同時，即兩個通路沒有相同的基因，則Jaccard 指數(shù)為0，當(dāng)兩個通路的基因集完全相同時，則Jaccard 指數(shù)為1。因此，所有通路對的Jaccard 指數(shù)在0 到1 變化不等。計算所有通路對的Jaccard 指數(shù)，結(jié)果顯示約30%通路對的重疊性小于0.6，大多數(shù)通路的Jaccard 指數(shù)小于0.2，這說明通路之間的重疊性較低，降低了因通路之間的重疊過高而引起的模型誤差。

3.2 基于通路活性的藥物敏感性預(yù)測性能

比較分析文獻［17］中的方法（DiffRank），本研究提出基于通路中所有基因推斷通路活性的方法，即PAG（All Gene of Pathway）。為了把基于通路模型和基于基因模型進行對比，還提出基于基因模型的方法AG（All Gene）。

PAG 方法和PHG 方法的不同在于PHG 方法在推斷通路活性時使用的是通路中高連接度的關(guān)鍵基因，而PAG 方法使用通路中所有基因來計算活性值，進而得到通路活性矩陣，以此作為預(yù)測模型的輸入。此外，基于基因模型的AG 方法是直接使用基因表達矩陣作為藥物敏感性預(yù)測模型的輸入，而不考慮通路信息，基因模型中的細胞系為樣本，基因表達值為特征。

本研究使用彈性網(wǎng)算法訓(xùn)練通路活性矩陣，10折交叉驗證選擇最優(yōu)參數(shù)，并使用最優(yōu)參數(shù)下的MSE 作為預(yù)測模型性能的評價標(biāo)準(zhǔn)。圖2 給出了基于CCLE 數(shù)據(jù)集中24 種藥物在4 種模型下進行藥物敏感性預(yù)測的結(jié)果。PHG 方法在17-AAG 等12 種藥物上具有最好的預(yù)測效果，在AZD6244 等6 種藥物上的預(yù)測效果是次好的；PAG 方法在Irinotecan 等4種藥物上具有最好的預(yù)測效果，在17-AAG等11種藥物上具有次好的效果。通過PHG 和PAG 對比分析，發(fā)現(xiàn)并不是通路中所有基因都會對藥物敏感性預(yù)測有幫助，只選取通路中連接緊密的基因進行預(yù)測可能更具有魯棒性。AG 方法在AZD6244 等7 種藥物上具有最好的效果，在Erlotinib 等6 種藥物上的預(yù)測效果是次好的。對比基于通路模型和基于基因模型可以發(fā)現(xiàn)基于通路模型有較好的預(yù)測性能?？偟膩碚f，對于一些藥物，使用基于通路中高連接度基因的計算分析方法取得了最好的預(yù)測效果，更有利于藥物敏感性預(yù)測。

圖2 不同模型對CCLE中24種藥物的預(yù)測性能Fig.2 Predictive performance of different models for 24 kinds of drugs in CCLE

3.3 模型再現(xiàn)性

本研究中的模型再現(xiàn)性是指在一個數(shù)據(jù)集上訓(xùn)練數(shù)據(jù)，在另一個數(shù)據(jù)集上測試數(shù)據(jù)，然后再交換數(shù)據(jù)集重新訓(xùn)練和測試。與CCLE 數(shù)據(jù)集相比，GDSC數(shù)據(jù)集中基因表達矩陣和通路數(shù)據(jù)推斷通路活性矩陣的特征數(shù)量較少。為了實驗的有效性，本研究從基于CCLE 基因表達譜推斷通路活性矩陣中隨機抽取和GDSC相同數(shù)量的特征，以便訓(xùn)練和預(yù)測。

對于給定的抗癌藥物，隨機選擇50 次相等數(shù)量的特征數(shù)據(jù)輸入到藥物敏感性模型，然后計算MSE，把50 次MSE 的平均值作為驗證模型再現(xiàn)性預(yù)測性能。共計算了24 種藥物在CCLE 數(shù)據(jù)集和GDSC 數(shù)據(jù)集上的預(yù)測性能，表1 中列出了Paclitaxel 等4 種藥物在CCLE 數(shù)據(jù)集和GDSC 數(shù)據(jù)集上的預(yù)測性能。在基于藥物Paclitaxel 敏感性預(yù)測的模型再現(xiàn)性中，當(dāng)以GDSC 數(shù)據(jù)作為訓(xùn)練集，CCLE 數(shù)據(jù)作為測試集時，MSE 為5.91；當(dāng)以CCLE 數(shù)據(jù)為訓(xùn)練集，GDSC 數(shù)據(jù)為測試集時，MSE 為5.52，這表明PHG 方法在藥物Paclitaxel 上的藥物敏感性預(yù)測具有較好的模型再現(xiàn)性。

表1 PHG方法在4種藥物的模型再現(xiàn)性Tab.1 Model reproducibility of PHG method in 4 kinds of drugs

另外，對于藥物17-AAG和PD0332991，以GDSC數(shù)據(jù)作為訓(xùn)練集訓(xùn)練模型，同時用此模型測試CCLE數(shù)據(jù)，發(fā)現(xiàn)具有較低的MSE，即較好的預(yù)測性能，然而當(dāng)以CCLE 數(shù)據(jù)訓(xùn)練模型，再以GDSC 數(shù)據(jù)測試模型，則有較高的誤差，這表明PHG 方法在這兩種藥物上使用基于GDSC 基因表達譜作為訓(xùn)練集時會得到較好的模型，具有較好的預(yù)測性能。相反，對于藥物PHA665752，以CCLE 基因表達數(shù)據(jù)作為訓(xùn)練集構(gòu)建藥物敏感性預(yù)測模型則會得到較好的預(yù)測性能。

3.4 識別藥物相關(guān)聯(lián)基因的通路

本研究把通路數(shù)據(jù)和基因表達譜整合得到通路活性評分，并以此構(gòu)建預(yù)測模型，進一步識別癌癥標(biāo)記物，從而給出生物學(xué)解釋。當(dāng)利用通路中高連接度基因數(shù)據(jù)和彈性網(wǎng)算法建立預(yù)測模型時，彈性網(wǎng)中非零系數(shù)對應(yīng)的特征是預(yù)測細胞對藥物反應(yīng)的重要數(shù)據(jù)［1］。因此，本研究采用了彈性網(wǎng)算法中非零系數(shù)統(tǒng)計與抗癌藥物相關(guān)聯(lián)基因的通路數(shù)量。在24種藥物中，19 種藥物包含靶向基因的通路都能識別出來（表2）。

表2 藥物相關(guān)基因的通路數(shù)量Tab.2 Number of pathways for drug-related genes

例如，對于藥物L(fēng)apatinib，使用PHG 方法能識別彈性網(wǎng)中非零系數(shù)對應(yīng)的MicroRNAs in cancer、Breast cancer 和 EGFR tyrosine kinase inhibitor resistance 等9 個特征通路，其中MicroRNAs in cancer通路包含ABCB1、EGFR 和ERBB2 等靶向基因，Breast cancer 通路包含EGFR 和ERBB2 等靶向基因。總的來說，基于通路高連接度基因的藥物敏感性預(yù)測模型能夠識別藥物相關(guān)聯(lián)基因的通路，具有更好的生物學(xué)解釋能力。

4 結(jié)論

本研究提出一種基于通路中高連接度基因的抗癌藥物敏感性預(yù)測方法（PHG）；對基因表達譜、通路數(shù)據(jù)和藥物敏感性IC50 值進行綜合分析，綜合考慮不同因素的作用，提取高連接度基因集合，然后計算通路活性矩陣，進而通過機器學(xué)習(xí)技術(shù)進行抗癌藥物敏感性預(yù)測分析，并把識別的標(biāo)記與已有研究進行對比分析，驗證基因/通路與藥物之間的聯(lián)系。實驗表明，基于通路中高連接度基因模型相比其他通路或基因模型有更好的預(yù)測性能。通路中并不是所有的基因都對藥物敏感性預(yù)測起到促進作用，而是一些關(guān)鍵基因更為重要。本研究提出的計算方法為通路活性預(yù)測模型的發(fā)展提供了參考。