基于集成學(xué)習(xí)的腫瘤藥物敏感性預(yù)測研究

黃鵬杰，林勇，張夢歡，呂琳，劉振浩，裴瀟倜，許林鋒,，謝鷺

1.上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093；2.上海生物信息技術(shù)研究中心，上海201203；3.中國科學(xué)院分子細胞科學(xué)卓越創(chuàng)新中心，上海200031

前言

腫瘤藥物敏感性預(yù)測研究是精準(zhǔn)醫(yī)學(xué)的重要分支。當(dāng)前，腫瘤精準(zhǔn)醫(yī)學(xué)不斷發(fā)展，基于生物標(biāo)志物和臨床病理特點的抗腫瘤靶向治療方法得到普遍應(yīng)用［1］。在合適的時間對不同癌癥患者采取個性化預(yù)防、治療和用藥指導(dǎo)，能使某些類型腫瘤的客觀緩解率和5年生存率均有較大程度提高，同時對癌癥病人的副作用比傳統(tǒng)方法小，這是精準(zhǔn)醫(yī)學(xué)在腫瘤領(lǐng)域應(yīng)用的成功典范［2］。

近年來，人工智能因為底層技術(shù)突破和海量數(shù)據(jù)驅(qū)動而不斷發(fā)展，而生物信息領(lǐng)域隨著第二代測序技術(shù)的發(fā)展產(chǎn)生大量的組學(xué)和藥物反應(yīng)數(shù)據(jù)，將兩者結(jié)合能極大輔助腫瘤精準(zhǔn)醫(yī)學(xué)的發(fā)展?；蚪M學(xué)研究不斷深入，大量研究表明腫瘤的產(chǎn)生、分化以及藥物治療和基因密切相關(guān)［3］，通過基因表達譜來預(yù)測腫瘤藥物的敏感性應(yīng)答成為當(dāng)下熱點。Geeleher等［4］通過乳腺癌的細胞系和藥物數(shù)據(jù)，建立了以嶺回歸（Ridge regression）算法為基礎(chǔ)的基因表達和藥物敏感性的回歸模型；Menden等［5］利用CGP數(shù)據(jù)集，通過神經(jīng)網(wǎng)絡(luò)算法來建立藥物應(yīng)答預(yù)測模型；基于隨機森林算法，Riddick 等［6］利用NCI-60 數(shù)據(jù)訓(xùn)練回歸模型來預(yù)測藥物敏感性；Gupta 等［7］研究了基于支持向量機和多種基因組類型的回歸方法。

這些方法不僅促進了癌癥藥物基因組學(xué)的發(fā)展，也為預(yù)測藥物敏感性提供了有效途徑，但當(dāng)前研究缺少對大批量藥物和細胞系基因數(shù)據(jù)的有效處理，同時模型所采用算法不多且精準(zhǔn)度還不夠高。本文工作較上述研究具有以下特點：（1）藥物種類多，且對細胞系和基因特征有效篩選用來建模；（2）采用集成學(xué)習(xí)算法，大部分藥物預(yù)測準(zhǔn)確率高，相比單模型有顯著提升；（3）對特定藥物進行功能分析，從生物學(xué)角度提供模型的可解釋性。

本文基于GDSC2 數(shù)據(jù)庫中809 個腫瘤細胞系樣本的RNA-seq 數(shù)據(jù)和183 種藥物與相關(guān)細胞系反應(yīng)的IC50 值，經(jīng)過樣本篩選和基因特征降維后，結(jié)合boosting集成和4種機器學(xué)習(xí)方法來構(gòu)建藥物敏感性分類預(yù)測模型。實驗結(jié)果表明，AdaBoost+SVM 集成方法相比單模型和其他集成模型，在整體藥物集上的ROC 曲線下面積（AUC）、精準(zhǔn)度（Precision）、準(zhǔn)確率（Accuracy）、查全率（Recall）、綜合評價指標(biāo)（F1）等指標(biāo)中具有更好效果。同時本文以藥物Nutlin-3A(-)為例進行特異性分析，通過特征排序篩選和富集分析，從生物學(xué)和統(tǒng)計學(xué)驗證了該藥物的作用靶向通路P53。綜上，本文的模型能對腫瘤病人的臨床精準(zhǔn)用藥起到指導(dǎo)和建議作用。

1 基于集成學(xué)習(xí)的腫瘤藥物敏感性預(yù)測方法

圖1 分類預(yù)測方法流程圖Fig.1 Flowchart of classification prediction method

1.1 數(shù)據(jù)集和預(yù)處理

本研究數(shù)據(jù)集來源于腫瘤藥敏多組學(xué)數(shù)據(jù)庫（Genomics of Drug Sensitivity in Cancer,GDSC），該數(shù)據(jù)庫有多于1 000 種人類癌癥細胞系的基因組特征和它們對不同抗癌藥物的敏感性實驗數(shù)據(jù)，是癌癥細胞藥物敏感性和藥物反應(yīng)分子標(biāo)志物信息的最大公共資源。GDSC 目前由GDSC1 和GDSC2 兩個數(shù)據(jù)庫組成，后者主要收錄2015年至今的測序和試驗結(jié)果，它使用代謝測定法(CellTitreGlo)來確定細胞活力，是基于改進的技術(shù)、設(shè)備和程序等得到的最新數(shù)據(jù)［8］。

本文的數(shù)據(jù)集為GDSC2，包含有809 個細胞系，198 個抗癌化合物，以及135 242 個藥物和細胞系反應(yīng)數(shù)據(jù)。其中藥物反應(yīng)指標(biāo)為IC50 值，即半抑制濃度，表示凋亡細胞與全部細胞數(shù)之比等于50%時所對應(yīng)的藥物濃度，反應(yīng)細胞對藥物的耐受程度。IC50 值越低，則說明細胞對藥物越敏感。此外，為了精確匹配相關(guān)細胞系的基因表達值，筆者從ArrayExpress網(wǎng)站獲得編號為E-MTAB-3610的數(shù)據(jù)，它包含GDSC 中所有細胞系的原始基因芯片［9］。接著，使用R 包Affy 讀取CEL 格式文件并提取為基因表達矩陣，同時利用該包中的Robust Multi-array Average(RMA)函數(shù)對原始數(shù)據(jù)進行標(biāo)準(zhǔn)化，以達到背景校正、分位數(shù)標(biāo)準(zhǔn)化、對探測強度進行l(wèi)og 轉(zhuǎn)換的目的。RMA標(biāo)準(zhǔn)化后，得到1 018個細胞系對應(yīng)的17 737 個基因表達值，其中去除基因名重復(fù)和空缺的數(shù)據(jù)后，還剩17 418個基因做后續(xù)分析。

1.2 細胞系樣本篩選和特征基因降維

對單個藥物處理來說，本文先對細胞系實驗樣本的IC50值進行歸一化處理：先計算均值，每個IC50值減去均值后再除以方差。然后對歸一化后的數(shù)據(jù)由小到大進行排序，前后20%占比的樣本分別視為有或無藥物應(yīng)答反應(yīng)，同時去除待驗證的中間樣本，將預(yù)測模型從回歸變?yōu)槎诸悊栴}。接著將每種藥物篩選后的數(shù)據(jù)集進行8/2 分，分別用來訓(xùn)練和驗證。此外，不同藥物實驗對應(yīng)的樣本數(shù)從48 到809不等，經(jīng)過以上步驟處理后有15 種藥物的樣本數(shù)小于20，樣本過少不利于建模和十折交叉驗證，因此這15種不在本次研究范圍。

接著，通過R語言匹配細胞系樣本的基因表達值作為待輸入特征，但選取的特征基因數(shù)過多，會提高計算復(fù)雜度、降低預(yù)測準(zhǔn)確性，因此對特征降維極其重要。本文將待訓(xùn)練樣本歸一化后的IC50值和單個基因的表達值進行皮爾遜相關(guān)性檢驗，通過R 中的Cor 函數(shù)選出相關(guān)系數(shù)絕對值大于0.3，P值小于0.05的基因。最后，對不同藥物來說，特征基因數(shù)從17 418降維到幾十至幾千不等。

1.3 基于集成學(xué)習(xí)的藥物敏感性預(yù)測建模

1.3.1 單模型建模本研究使用的單模型是4 種經(jīng)典的機器學(xué)習(xí)模型，分別為支持向量機（Support Vector Machine,SVM）、隨機森林（Random Forest,RF），極端隨機樹（Extremely Randomized Trees,ET）和基于Bagging 的決策樹（Decision Tree Based on Bagging,BT）。其中，本研究集成的主要基分類器為SVM。

SVM是基于結(jié)構(gòu)風(fēng)險最小化和VC維理論的監(jiān)督式學(xué)習(xí)方法，通過構(gòu)造分隔面將數(shù)據(jù)進行分離，在解決非線性、高維模式識別和大樣本分類中表現(xiàn)出獨特優(yōu)勢［10］。建模中，給定的腫瘤藥物預(yù)測的樣本集{(xi,yi),i=1,2,…,n}，xi表示預(yù)測因子基因特征向量（L維），yi表示某細胞系對某藥物是否有應(yīng)答。對于這種高維非線性預(yù)測問題，SVM通過將原始空間映射到高維特征空間來建立分類預(yù)測模型。分為以下幾步：（1）建立初始超平面：w·φ(xi)+b=0，其中，w為權(quán)值矢量，b為閾值，φ為目標(biāo)函數(shù)。（2）尋找最優(yōu)超平面，其可以歸納為求解如下優(yōu)化問題：s.t.(w·φ(xi)+b)-1+ξi≥0,ξi≥0，其中C為懲罰參數(shù)，ξi為非負松弛因子。（3）接著求解得到最優(yōu)分類決策函數(shù)：，其中，為最優(yōu)拉格朗日乘子，b*為分類閾值且為核函數(shù)，sgn(*)表示返回整型變量的函數(shù)，xj表示不同于xi的樣本點。

RF、ET、BT 都可看做Bagging 集成的變種，本文將其視為Boosting 集成的基分類器。RF 將bootstrap重抽樣方法和決策樹算法相結(jié)合，能對重要特征進行排序，是性能良好的用于分類或回歸的有監(jiān)督算法［11］。ET 與RF 相似，都由許多決策樹構(gòu)成，但它是完全隨機得到分叉值進而對回歸樹分叉，不同于RF在一個隨機子集內(nèi)得到最佳分叉屬性［12］。此外，該方法每棵回歸樹都使用全部訓(xùn)練樣本。BT是完全用Bagging分類器集成決策樹的方法。

[教學(xué)片段4]執(zhí)教老師在講不同國家贈送禮物的文化時，提到朋友結(jié)婚送禮的問題，讓學(xué)生給她建議。她設(shè)計了三類情況，中國人的婚禮、英國人的婚禮和中國女生與英國男生的婚禮，通過追問幫助學(xué)生認識贈送禮要隨當(dāng)?shù)氐牧?xí)俗，可以幫助學(xué)生理清關(guān)鍵信息的內(nèi)涵。提問如下：

1.3.2 Boosting 集成學(xué)習(xí)建模AdaBoost 是最常見的一種Boosting 算法，能有效應(yīng)用在兩類、多類單標(biāo)簽等分類問題中，具有較強自適應(yīng)能力、速度快和算法復(fù)雜度低等優(yōu)點［13］。其主要思想是給每個訓(xùn)練樣本賦予一個初始權(quán)值，再根據(jù)每次訓(xùn)練的樣本分類是否正確以及總體準(zhǔn)確率，來自動調(diào)整權(quán)值。接著將調(diào)整后的新數(shù)據(jù)集傳給下個分類器，最后對每次訓(xùn)練的基分類器進行加權(quán)組合得到最終決策分類器。

本研究將該方法用在腫瘤細胞系對藥物敏感性的分類中。設(shè)輸入的訓(xùn)練樣本為：{（x1,y1),(x2,y2），…，（xi,yi），…，（xn,yn）}，其中xi中表示輸入的帶有基因表達信息的細胞系樣本，yi∈{0,1}分別表示正樣本和負樣本，其中正樣本數(shù)為l，負樣本數(shù)m，n=l+m。具體步驟如下：

（1）初始化每個細胞系樣本的權(quán)重wi,i∈D（i），在每個弱分類器t中（t=1，…,T，T為總個數(shù)），其權(quán)重可以歸一化為一個概率分布：

（2）對每個特征訓(xùn)練一個弱分類器hj。由于本研究為二值分類，可通過遍歷得到閾值，該閾值使錯分的樣本最小且錯誤率最低。接著計算包含所有特征的弱分類器的加權(quán)錯誤率：

（3）選取錯誤率εt最小的弱分類器ht，并按照此最佳ht調(diào)整權(quán)重：

其中，εi=0 表示被正確地分類；εi=1 表示被錯誤地分類。

（4）加權(quán)得到最終的決策強分類器：

其中αt的計算公式為：經(jīng)過以上一系列的迭代和集成處理后，我們可以得到藥物敏感性預(yù)測的最終強分類器，來判定腫瘤細胞系能否對該藥物產(chǎn)生應(yīng)答反應(yīng)。AdaBoost 默認的基分類器為決策樹，本研究分別將SVM、RF、ET、BT 作為基分類器進行同質(zhì)集成?；诸惼鱾€數(shù)對預(yù)測精度會有影響，本文分別設(shè)置為3、10、20、50進行實驗，發(fā)現(xiàn)個數(shù)為10時達到理想效果，即使增加到50 準(zhǔn)確率沒有明顯提高，反而會加大內(nèi)存消耗和訓(xùn)練成本。因此，本文基分類器數(shù)量定為10。

1.3.3 十折交叉驗證交叉驗證（Cross Validation），也稱循環(huán)估計，是一種統(tǒng)計學(xué)上將樣本切割成較小子集的使用方法，常用于評估統(tǒng)計分析、機器學(xué)習(xí)算法對獨立于訓(xùn)練數(shù)據(jù)的數(shù)據(jù)集的泛化能力［14］。

十折交叉驗證（10-fold cross validation），將數(shù)據(jù)集分成10 份，輪流將其中9 份做訓(xùn)練，1 份做驗證，10次結(jié)果的均值作為對算法精度的估計。本研究建模過程中均使用了十折交叉驗證求均值的方法，其優(yōu)勢在于同時重復(fù)運用隨機產(chǎn)生的子樣本進行訓(xùn)練和驗證，以求藥物敏感性預(yù)測結(jié)果更精準(zhǔn)。

2 實驗驗證和結(jié)果分析

2.1 實驗評估指標(biāo)

為了評估模型好壞，對比4 種分類模型的預(yù)測精度和泛化能力。本研究采用以下常用評價指標(biāo)：AUC，Precision，Recall，F(xiàn)1值，Accuracy。在二分類中，“正例（positive）”和“負例（negative）”指分類標(biāo)簽，而“真（true）”和“假（false）”指模型預(yù)測結(jié)果是否與測試標(biāo)簽相對應(yīng)?？芍?，TP、TN 表示分類結(jié)果正確；FP、FN 表示分類結(jié)果錯誤。各項指標(biāo)的定義如下：

其中，TP 為真陽性，表示實驗細胞系樣本有藥物應(yīng)答，預(yù)測模型中也有。FP 為假陽性，表示實驗細胞系樣本無藥物應(yīng)答，但預(yù)測模型中有。FN為假陰性，表示實驗細胞系樣本有藥物應(yīng)答，但預(yù)測模型中沒有。ROC 曲線在坐標(biāo)系中越偏向左上角，即AUC 越大，則表示分類器的性能越好。

2.2 實驗結(jié)果與分析

2.2.1 藥物模型整體比較針對183 種藥物的建模結(jié)果，本文詳細比較了Boosting 集成方法和4 種單模型的各項評估指標(biāo)。由圖2 所有藥物的AUC 分布情況可知，AdaBoost+SVM 集成模型有108 種藥物集的AUC＞0.9，AUC 在＞0.95、0.9～0.95、0.85～0.9、0.8～0.85、0.75～0.8、＜0.75 等區(qū)間的藥物數(shù)占比依次為7.1%、51.9%、24.6%、11.4%、3.8%和1.1%，總體效果明顯優(yōu)于SVM、RF、ET、BT 和它們的集成模型。此外，其他3 種集成模型與各自的基學(xué)習(xí)器相比在大部分藥物的準(zhǔn)確率上也有所提升。

圖2 183種藥物在單模型和boosting集成模型的AUC分布情況Fig.2 AUC distribution of 183 drugs in the models based on a learner only and Boosting ensemble models

得到整體分析結(jié)果后，筆者選取在所有模型中都取得較高準(zhǔn)確率的藥物做進一步分析，分別為Venetoclax、OSI-027、Nutlin-3a(-)、AZD5991、Vorinostat、Oxaliplatin、Dasatinib、 Niraparib、 Paclitaxel、 Daporinad、Camptothecin、Cisplatin、MIRA-1。圖3是這13種藥物在集成模型和SVM中的AUC對比圖，可知，AdaBoost+SVM模型的預(yù)測效果最優(yōu)，且AUC都大于0.95。值得注意的是，對單個藥物來說，取得最優(yōu)效果的不一定都是AdaBoost+SVM，其他集成方法在個別藥物數(shù)據(jù)集中也可能成為最優(yōu)建模方法。

圖3 高準(zhǔn)確率藥物預(yù)測模型比較Fig.3 Comparison of high-accuracy drug prediction models

為了更詳細展示每種藥物在不同方法建模中的效果，筆者以藥物Nutlin-3a(-)為例進行分析。表1是該藥物在Boosting 集成模型和4 種單模型中的AUC、Accuracy、Precision、Recall以及F1值的具體數(shù)值和標(biāo)準(zhǔn)差。圖4為不同的模型進行十折交叉訓(xùn)練的ROC曲線，可以直觀發(fā)現(xiàn)該藥物在AdaBoost+SVM中效果最好，整體準(zhǔn)確率相比單模型高出約2.5%，相比其他集成模型高出約1.5%，且各項指標(biāo)的標(biāo)準(zhǔn)差相對最小。

表1 藥物Nutlin-3a(-)在所有模型中的結(jié)果比較（±s）Tab.1 Comparison of prediction results of Nutlin-3a(-)in all models(Mean±SD)

表1 藥物Nutlin-3a(-)在所有模型中的結(jié)果比較（±s）Tab.1 Comparison of prediction results of Nutlin-3a(-)in all models(Mean±SD)

模型AdaBoost+SVM AdaBoost+RF AdaBoost+ET AdaBoost+BT SVM RF ET BT AUC 0.977±0.017 0.966±0.023 0.963±0.028 0.971±0.022 0.952±0.037 0.955±0.04 0.956±0.031 0.962±0.017 Accuracy 0.925±0.043 0.883±0.061 0.879±0.067 0.922±0.046 0.867±0.072 0.87±0.077 0.889±0.069 0.912±0.034 Precision 0.945±0.063 0.933±0.04 0.917±0.068 0.927±0.074 0.951±0.055 0.913±0.074 0.935±0.057 0.92±0.069 Recall 0.908±0.08 0.823±0.115 0.836±0.123 0.911±0.058 0.77±0.136 0.816±0.125 0.835±0.124 0.908±0.046 F1值0.922±0.048 0.871±0.073 0.869±0.08 0.912±0.047 0.846±0.09 0.857±0.088 0.877±0.084 0.902±0.033

2.2.2 藥物特異性分析本文在所建模型基礎(chǔ)上以Nutlin-3a(-)為例進行了藥物特異性分析。從GDSC數(shù)據(jù)庫可知該抑制劑已驗證的藥物靶點為MDM2，作用通路為P53。如圖5所示，本文通過隨機森林模型選取該藥數(shù)據(jù)集的重要特征基因進行富集分析，從生物學(xué)和統(tǒng)計學(xué)兩方面驗證了該藥物的P53 靶向通路［15］。該生物功能分析為模型提供了可解釋性和生物理論支撐。

3 總結(jié)與展望

本文針對GDSC2 中的183 種抗癌藥物和809 個細胞系數(shù)據(jù)，經(jīng)過有效的樣本篩選和基因特征降維后，首次使用AdaBoost+SVM集成方法建模來預(yù)測不同細胞系對藥物的敏感性。此外，通過比較AUC、Accuracy、Precision、Recall 和F1值后，發(fā)現(xiàn)該模型相比4 種單模型和其他集成方法在整體藥物數(shù)據(jù)集中效果最優(yōu)。與現(xiàn)有研究的模型相比，不僅具有較高精確度，而且通過藥物特異性分析為模型提供了生物可解釋性。以上研究證實了將集成學(xué)習(xí)方法應(yīng)用在腫瘤藥物敏感性預(yù)測的可行性，能在臨床用藥的篩選中發(fā)揮作用。實際應(yīng)用中，可以先對病人的腫瘤組織進行基因測序，再使用模型去預(yù)測特定抗癌抑制劑對該病人是否有效，以達到提高用藥療效和減輕副作用的目的［16］。

圖4 Nutlin-3a(-)在AdaBoost+SVM 和其他模型中的ROC曲線Fig.4 ROC curves of Nutlin-3a(-)in AdaBoost+SVM and other models

圖5 Nutlin-3a(-)藥物集重要特征的KEGG富集結(jié)果Fig.5 KEGG enrichment result of important characteristics of Nutlin-3A(-)drug set

本研究的不足之處在于目前的模型只有基因表達數(shù)據(jù)作為輸入特征，雖然基于單一組學(xué)進行疾病研究已經(jīng)發(fā)現(xiàn)了諸多新的致病因子，但疾病的發(fā)生發(fā)展是一個復(fù)雜的網(wǎng)絡(luò)，基因拷貝數(shù)變異、體細胞突變、甲基化的異常等諸多因素都會影響生命體特征和藥物反應(yīng)［17］。因此，之后的研究會處理以上多組學(xué)數(shù)據(jù)作為訓(xùn)練特征，同時嘗試更多方法建模，以達到更好的腫瘤藥物預(yù)測效果。