高沖，秦玉芳，陳明

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

前言

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

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

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

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

1 數據集及預處理

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

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

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

2 模型

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

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

計算每個基因在通路互相作用網絡表中的度，由于通路網絡中高連接度的Hub 基因對整個通路的功能起著更關鍵的作用［19］，所以從該網絡表中選擇高連接度Hub 基因來進行分析。將基因的度降序排序，選擇排名在前10%的基因作為Hub 基因，圖1中的Gh1,Gh2,…,Ghk為通路p1的Hub 基因；計算Hub基因表達值的平均值作為活性值，然后合并每個細胞系樣本中的活性值得到該通路的活性向量。活性值計算公式如下：

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

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

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

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

2.2 建立藥物反應預測模型

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

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

3 結果與討論

3.1 通路之間的重疊性

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

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

3.2 基于通路活性的藥物敏感性預測性能

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

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

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

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

3.3 模型再現(xiàn)性

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

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

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

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

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

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

表2 藥物相關基因的通路數量Tab.2 Number of pathways for drug-related genes

例如，對于藥物Lapatinib，使用PHG 方法能識別彈性網中非零系數對應的MicroRNAs in cancer、Breast cancer 和 EGFR tyrosine kinase inhibitor resistance 等9 個特征通路，其中MicroRNAs in cancer通路包含ABCB1、EGFR 和ERBB2 等靶向基因，Breast cancer 通路包含EGFR 和ERBB2 等靶向基因?？偟膩碚f，基于通路高連接度基因的藥物敏感性預測模型能夠識別藥物相關聯(lián)基因的通路，具有更好的生物學解釋能力。

4 結論

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