韓 睿 郭成安

在精準(zhǔn)腫瘤學(xué)研究領(lǐng)域，研究人員目前正針對如何能為每位患者量身定制最佳的治療方案進行深入探索[1].實施精準(zhǔn)腫瘤學(xué)醫(yī)療技術(shù)的一個重要先決條件是能根據(jù)患者腫瘤的特定分子特征正確識別與所用藥物反應(yīng)相關(guān)的生物標(biāo)記物.然而，藥物反應(yīng)與腫瘤分子的生物特征之間的關(guān)系十分復(fù)雜，目前尚未提出有效模型進行準(zhǔn)確描述.當(dāng)前的機器學(xué)習(xí)技術(shù)、不斷積累的醫(yī)療大數(shù)據(jù)、以及將二者結(jié)合而開展的研究工作，為解決這一問題提供一個有效途徑.近年來，隨著個體基因測序技術(shù)的日益完善，針對患者個人的腫瘤衍生分子數(shù)據(jù)，對用于治療腫瘤的藥物活性進行無偏見、精確地建模與預(yù)測，已成為實現(xiàn)腫瘤患者個體化、精準(zhǔn)治療的關(guān)鍵.因此，該項研究不僅可促進精準(zhǔn)腫瘤學(xué)的發(fā)展，還能顯著降低臨床實施的成本，使更多的患者受益.

目前，一些國際醫(yī)學(xué)研究機構(gòu)已公開發(fā)布臨床藥物基因組學(xué)數(shù)據(jù)集，為科研人員從事抗腫瘤藥物反應(yīng)預(yù)測研究提供重要的基礎(chǔ)數(shù)據(jù).2004年，英國威康信托基金會桑格研究所發(fā)布COSMIC數(shù)據(jù)庫[2].該數(shù)據(jù)庫匯集人類腫瘤細(xì)胞系的基因突變信息，其中一個專門組成部分CCLP(COSMIC Cell Lines Project)是對最常使用的1 000多種腫瘤細(xì)胞系的深入分析數(shù)據(jù).2009年，英國桑格研究院推出腫瘤藥物敏感性基因組學(xué)GDSC(Genomics of Drug Sensitivity in Cancer)數(shù)據(jù)庫[3]，提供1 000多種腫瘤細(xì)胞系的染色體拷貝數(shù)變異、基因表達等基因組學(xué)數(shù)據(jù)及265種抗腫瘤藥物作用在腫瘤細(xì)胞系的半抑制濃度(Half Maximal Inhibitory Concentration),即IC50值.

抗腫瘤藥物反應(yīng)預(yù)測研究大多基于高通量測序技術(shù)收集的不同基因組水平信息[4].相比上萬維的基因組學(xué)特征，腫瘤細(xì)胞系的數(shù)量通常只有一千個，而在全部的基因組學(xué)特征中，可能只有一部分與某一抗腫瘤藥物反應(yīng)相關(guān).因此，如何從原始的基因組學(xué)數(shù)據(jù)中提取合適的特征子集是預(yù)測抗腫瘤藥物反應(yīng)的關(guān)鍵.Xu等[5]提出AutoBorutaRF,采用自編碼器網(wǎng)絡(luò)篩選一個適當(dāng)維度的特征，再利用Boruta算法進一步確定隨機森林的輸入特征.Vougas等[6]提出DLNNs(Deep Learning Neural Networks),基于關(guān)聯(lián)規(guī)則挖掘算法,從基因表達、拷貝數(shù)變異和基因突變?nèi)N基因組學(xué)數(shù)據(jù)中提取單尺度的基因組學(xué)特征，并從中隨機挑選200維特征作為預(yù)測抗腫瘤藥物反應(yīng)的每個子分類器的輸入特征.Sharifi-Noghabi等[7]采用3個編碼子網(wǎng)絡(luò)，分別學(xué)習(xí)基因表達、基因突變、拷貝數(shù)變異三種基因組學(xué)數(shù)據(jù)的特征表示，并將上述單組學(xué)特征連接成一個多組學(xué)特征,作為預(yù)測子網(wǎng)絡(luò)的輸入特征.

用于抗腫瘤藥物反應(yīng)預(yù)測的機器學(xué)習(xí)方法主要包括：邏輯回歸、彈性網(wǎng)絡(luò)、隨機森林、KNN(K-Nearest Neighbor)，SVM(Support Vector Machine)及淺層神經(jīng)網(wǎng)絡(luò)等.近年來，深度學(xué)習(xí)成為藥物反應(yīng)預(yù)測領(lǐng)域中一種值得研究的方法[8]，可分為兩類.1)基于深度神經(jīng)網(wǎng)絡(luò)(Deep Neural Networks, DNN)的方法.Vougas等[6]將多個子分類器的輸出加權(quán)平均后,獲得最終預(yù)測抗腫瘤藥物反應(yīng)的輸出.每個子分類器均采用包含3個隱層的DNN結(jié)構(gòu).Sharifi-Noghabi等[7]提出MOLI(Multi-omics Late Integra-tion)，采用一個集成的多組學(xué)特征作為輸入，預(yù)測靶向藥物家族的藥物反應(yīng).其中的編碼子網(wǎng)絡(luò)和預(yù)測子網(wǎng)絡(luò)均采用單隱層的DNN結(jié)構(gòu).Li等[9]提出DeepDSC(A Deep Learning Method to Predict Drug Sensitivity of Cancer Cell Lines)，將基因表達特征和藥物的化學(xué)特征一起作為輸入，用于預(yù)測細(xì)胞系的藥物敏感性，預(yù)測網(wǎng)絡(luò)采用包含4個隱層的DNN結(jié)構(gòu).Zhu等[10]提出遷移學(xué)習(xí)框架，建立一般藥物反應(yīng)的預(yù)測模型，采用7個隱層的DNN結(jié)構(gòu).2)基于深度卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks, CNN)的方法.Chang等[4]提出CDRscan(Cancer Drug Response Profile Scan),實現(xiàn)對抗腫瘤藥物的IC50值的預(yù)測，其中2個輸入CNN分支分別對細(xì)胞系的基因突變和藥物的分子結(jié)構(gòu)信息進行建模.在文獻[4]的基礎(chǔ)上，Bai等[11]進一步針對聯(lián)合采用基因測序信息和藥物分子結(jié)構(gòu)信息的方法進行研究，提出采用雙輸入和雙輸出的深度CNN，實現(xiàn)對所用藥物的敏感性及IC50值的預(yù)測，取得較優(yōu)效果.

雖然現(xiàn)有方法在解決抗腫瘤藥物反應(yīng)預(yù)測的關(guān)鍵問題上取得一定進展，但仍存在如下問題.1)如何從原始的腫瘤細(xì)胞系提取合適的基因組學(xué)特征.Vougas等[6]基于關(guān)聯(lián)規(guī)則挖掘，提取單尺度的基因組學(xué)特征，從中隨機挑選多組特征子集分別作為待選特征，雖然取得一定效果，因其隨機性較大，難以覆蓋全部基因組學(xué)特征，所以無法保證具有穩(wěn)定的性能.同時由于該方法需要采用多組待選特征分別進行訓(xùn)練，需要的預(yù)測子網(wǎng)絡(luò)數(shù)目較多，造成訓(xùn)練耗時過長.另外，采用單一規(guī)則提取基因組學(xué)特征，可能會排除某些對預(yù)測模型有用的輸入特征.2)在已有用于抗腫瘤藥物反應(yīng)研究的醫(yī)學(xué)數(shù)據(jù)樣本中，正樣本數(shù)據(jù)相對偏少，對于訓(xùn)練過深的網(wǎng)絡(luò)模型容易產(chǎn)生過擬合問題，泛化能力嚴(yán)重下降.

針對上述的特征提取和預(yù)測模型設(shè)計等問題，本文提出抗腫瘤藥物反應(yīng)預(yù)測方法，運用機器學(xué)習(xí)技術(shù)，對患者腫瘤基因測序數(shù)據(jù)進行處理、特征提取及建模，預(yù)測各種不同的抗腫瘤藥物的療效反應(yīng).首先，提出基于多尺度關(guān)聯(lián)規(guī)則的數(shù)據(jù)挖掘方法，對基因組學(xué)數(shù)據(jù)進行不同尺度的特征挑選.進而通過累積窗函數(shù)對挑選后的基因組學(xué)數(shù)據(jù)進行局部累積，進一步進行數(shù)據(jù)壓縮，提取具有較強整體表達性的基因特征信息.然后，以多層全連接神經(jīng)網(wǎng)絡(luò)為模型、以提取的多尺度累積基因特征為輸入樣本，進行訓(xùn)練和建模.最后，分別采用特征融合和決策融合，實現(xiàn)某一腫瘤基因測序數(shù)據(jù)對于各種不同抗腫瘤藥物反應(yīng)結(jié)果的預(yù)測.在COSMIC、GDSC數(shù)據(jù)庫上的仿真實驗表明，本文方法在敏感性、特異性、準(zhǔn)確率、特性曲線面積值等關(guān)鍵性能指標(biāo)上均取得較優(yōu)值.

1 抗腫瘤藥物反應(yīng)預(yù)測方法

1.1 總體設(shè)計方案

針對抗腫瘤藥物反應(yīng)預(yù)測研究中的關(guān)鍵問題，本文提出抗腫瘤藥物反應(yīng)預(yù)測方法，運用機器學(xué)習(xí)技術(shù)，對患者腫瘤基因測序數(shù)據(jù)進行處理、特征提取及建模，實現(xiàn)對各種不同的抗腫瘤藥物的療效反應(yīng)進行預(yù)測的方法.本文方法主要由數(shù)據(jù)預(yù)處理、多尺度局部累積特征提取、多層全連接神經(jīng)網(wǎng)絡(luò)和融合與決策四個模塊組成，具體結(jié)構(gòu)框架如圖1所示.

圖1 抗腫瘤藥物反應(yīng)預(yù)測方法框圖Fig.1 Flow chart for prediction of antitumor drug response

在數(shù)據(jù)預(yù)處理模塊中，本文對COSMIC、GDSC數(shù)據(jù)庫提供的基因組學(xué)數(shù)據(jù)和藥物反應(yīng)數(shù)據(jù)進行篩選和數(shù)值化處理.在多尺度局部累積特征提取模塊中，首先采用多尺度關(guān)聯(lián)規(guī)則挖掘方法對人體腫瘤細(xì)胞系的基因表達、基因突變和拷貝數(shù)變異數(shù)據(jù)進行特征挑選，保留其中與藥物反應(yīng)存在關(guān)聯(lián)的有用信息,去除無關(guān)信息.然后，在多尺度特征的基礎(chǔ)上采用局部累積處理，選擇最優(yōu)的局部累積特征觀察窗提取每種尺度的局部累積特征.在多層全連接神經(jīng)網(wǎng)絡(luò)模塊中，以3個3層全連接神經(jīng)網(wǎng)絡(luò)為模型，以提取的多尺度累積基因組學(xué)特征為輸入樣本，分別進行訓(xùn)練和建模.在融合與決策模塊中，分別采用決策融合和特征融合，設(shè)計具體的融合和決策算法，實現(xiàn)某一腫瘤基因測序數(shù)據(jù)對于各種不同抗腫瘤藥物反應(yīng)結(jié)果的預(yù)測.

1.2 數(shù)據(jù)預(yù)處理

本文采用CCLP、GDSC數(shù)據(jù)庫作為源數(shù)據(jù).CCLP數(shù)據(jù)庫提供1 001株人類腫瘤細(xì)胞系的基因表達、拷貝數(shù)變異和基因突變這3種常見的基因組學(xué)(特征個數(shù)大于60 000)的深入分析數(shù)據(jù).GDSC數(shù)據(jù)庫提供這1 001株腫瘤細(xì)胞系對265種抗腫瘤藥物的藥物反應(yīng)測量值(如ln(IC50)值).為了不影響方法的預(yù)測性能，本文選取其中的196種藥物，去除藥物反應(yīng)數(shù)據(jù)缺失比例大于20%的69種.

CCLP數(shù)據(jù)庫提供的基因組學(xué)數(shù)據(jù)包括：16 445維基因表達譜、24 858維拷貝數(shù)變異和19 384維基因突變.文獻[6]采用的預(yù)處理方法是將上述數(shù)據(jù)量化處理為符號表達的基因狀態(tài)信息.為了與DLNNs進行公平對比，本文采用與其相同的基因狀態(tài)信息.同時依照文獻[12]對上述基因狀態(tài)的符號表達進行數(shù)值化處理，以便后續(xù)的網(wǎng)絡(luò)訓(xùn)練.具體為，在基因表達狀態(tài)中，采用1表示基因的高表達，-1表示低表達，0表示介于二者之間.在拷貝數(shù)變異狀態(tài)中，采用1表示超出重復(fù)片段正常值的基因，-1表示有缺失片段的基因，0表示正常.在基因突變狀態(tài)中，采用1表示攜帶體細(xì)胞除沉默突變之外的所有單點突變，其余突變類型用0表示.

在GDSC數(shù)據(jù)庫中，文獻[6]采用零-均值規(guī)范化的ln(IC50)值劃分預(yù)測網(wǎng)絡(luò)的輸出類別.其中，閾值-1用于劃分藥物的敏感性(Sensitive)或不敏感性(Non-sensitive)，閾值1用于劃分藥物的耐藥性(Resistant)或非耐藥性(Non-Resistant).本文采用的理想標(biāo)簽是與對藥物的敏感性、耐藥性進行預(yù)測有關(guān).如需對抗腫瘤藥物的敏感性進行分類預(yù)測，根據(jù)文獻[6]，在訓(xùn)練過程中：對數(shù)據(jù)庫上標(biāo)有敏感的細(xì)胞系，將標(biāo)簽值置為1，作為正樣本；將其余未標(biāo)有敏感的細(xì)胞系的標(biāo)簽值置為0，作為負(fù)樣本.如需對抗腫瘤藥物的耐藥性進行分類預(yù)測：對標(biāo)有耐藥的細(xì)胞系,將標(biāo)簽值置為1，作為正樣本；將其余未標(biāo)有耐藥的細(xì)胞系的標(biāo)簽值置為0，作為負(fù)樣本.

1.3 多尺度累積特征提取

CCLP數(shù)據(jù)庫提供超過6萬維的基因組學(xué)數(shù)據(jù)，如何從中提取合適的特征集并有效壓縮特征維度是預(yù)測抗腫瘤藥物反應(yīng)的關(guān)鍵問題之一.文獻[6]和文獻[13]通過設(shè)定關(guān)聯(lián)規(guī)則中一組確定的支持度閾值S0和可信度閾值C0，采用Apriori算法[14]提取與藥物反應(yīng)密切關(guān)聯(lián)的基因組學(xué)特征.這種在單一S0和C0下提取的基因組學(xué)特征稱為單尺度特征.單尺度特征提取方法采用單一的特征挑選規(guī)則，可能會遺漏原始基因組學(xué)數(shù)據(jù)中的某些有用信息.針對這一問題，本文提出多尺度累積特征提取方法，盡可能地從原始基因組學(xué)數(shù)據(jù)中提取更多有用的基因組學(xué)特征，進一步提高基因特征的整體表達性.

1.3.1 多尺度關(guān)聯(lián)規(guī)則挖掘

區(qū)別于單一的特征挑選規(guī)則(設(shè)定關(guān)聯(lián)規(guī)則中單一的S0和C0)，通過分析文獻[6]和文獻[13]，本文發(fā)現(xiàn)設(shè)定多組S0和C0，可提取不同尺度的基因組學(xué)特征.這相當(dāng)于既細(xì)致又全面地分析原始的基因組學(xué)數(shù)據(jù).采用多尺度關(guān)聯(lián)規(guī)則挖掘方法提取多尺度基因組學(xué)特征的一般步驟可描述為：首先設(shè)定V種不同的強關(guān)聯(lián)規(guī)則，每種強關(guān)聯(lián)規(guī)則均包含一組確定的支持度閾值Sv和可信度閾值Cv，v=1,2,…,V.在確定的Sv和Cv下，采用Apriori算法從訓(xùn)練集T中挖掘NFv個形如“基因?某一種藥物反應(yīng)”的強關(guān)聯(lián)規(guī)則(滿足支持度大于等于Sv，可信度大于等于Cv的條件).提取NFv個強關(guān)聯(lián)規(guī)則中的基因，組成與某一抗腫瘤藥物反應(yīng)密切關(guān)聯(lián)的一維基因組學(xué)特征Xv.當(dāng)v≤V時，重復(fù)上述操作，直到獲得一組多尺度基因組學(xué)特征X=[X1,X2,…,Xv].文獻[6]設(shè)定S0=C0=0.004 48，對應(yīng)本文采用的訓(xùn)練樣本數(shù)669，其值為3/669.因此，本文選擇1/669、2/669和3/669這3種尺度的支持度閾值和可信度閾值，用于提取多尺度基因組學(xué)特征，即V=3.而0/669這一尺度相當(dāng)于對原始基因組學(xué)數(shù)據(jù)未做任何特征提取，本文并未采用.

對于尺度為v的基因組學(xué)特征數(shù)據(jù)集Xv={XEv,XMv,XCv}，由基因表達特征XEv、基因突變特征XMv和拷貝數(shù)變異特征XCv組成.不同尺度下得到XEv、XMv和XCv的特征維度不同，大小與特征提取時采用的S0和C0呈反比例關(guān)系.例如，在較大的S0和C0下，特征挑選規(guī)則要求“基因Gi高表達”和“藥物Dj的敏感性”同時發(fā)生的概率較高.在此特征挑選規(guī)則下，原始基因組學(xué)數(shù)據(jù)中大量與藥物反應(yīng)關(guān)聯(lián)程度較低的基因被刪除，只保留重要的基因組學(xué)特征，此時Xv的特征維度較低.而在較小的S0和C0下，提取的基因組學(xué)特征會保留大量的原始基因組學(xué)數(shù)據(jù)，此時Xv的特征維度較高.而以高維的基因組學(xué)特征為輸入，會造成后續(xù)預(yù)測子網(wǎng)絡(luò)的處理難度增加.針對這一問題，本文采用局部累積特征提取方法，可進一步降低Xv的特征維度.

1.3.2 局部累積特征提取

為了進一步降低基因組學(xué)特征的維度，緩解后續(xù)預(yù)測子網(wǎng)絡(luò)的處理難度，文獻[6]采用從基因組學(xué)特征Xv中有放回地隨機挑選一小部分特征的方法.此方法具有很強的隨機性，很難保證多次隨機挑選的基因組學(xué)特征會全部覆蓋Xv.此外，由于人類基因組中通常只有小部分基因發(fā)生突變，基因組學(xué)數(shù)據(jù)在很大程度上會呈現(xiàn)稀疏性.針對上述問題，本文提出局部累積特征提取方法，對每種尺度的基因組學(xué)特征進行以局部累積特征觀察窗為單元的累積處理.設(shè)基因序列集為

{xk(n)|k=0,1,…,K-1,n=0,1,…,N-1}，

其中,xk(n)為第k個基因序列，K為全部基因序列個數(shù)，N為基因序列長度.設(shè)Fk(i)為xk(n)的第i段累積特征，則

(1)

其中,wM(m)為寬度為M的累積窗函數(shù)，LM為累積特征序列的長度.wM(m)的取值可根據(jù)同一種類型的基因組學(xué)數(shù)據(jù)內(nèi)部的重要性確定.由于本文采用的基因組學(xué)數(shù)據(jù)是經(jīng)量化后的基因狀態(tài)信息，該量化操作會損失基因組學(xué)數(shù)據(jù)內(nèi)部的重要性信息，因此，本文同等對待在同種類型的基因組學(xué)數(shù)據(jù)中發(fā)生突變的基因，將wM(m)取為如下矩形窗：

對于提取基因突變序列的局部累積特征(記為FMk(i))，將式(1)中的xk(n)采用相應(yīng)的基因突變序列進行計算即可.在提取局部累積特征的過程中，若最后一段包含的基因組學(xué)數(shù)據(jù)長度不足M，對缺少的數(shù)據(jù)進行補0操作.此外，考慮到基因表達狀態(tài)包含1和-1兩種非0情況，為了防止累積運算中發(fā)生正負(fù)抵消的情況，本文將基因表達序列(仍用xk(n)表示)分解成2個序列：

然后分別計算局部累積特征：

(2)

拷貝數(shù)變異序列與基因表達序列相同,都是取值為1或-1的序列，局部累積特征提取過程與基因表達序列一致，仍按照式(2)求取局部累積特征.

從上述累積特征提取過程可看出，若在特征觀察窗范圍內(nèi)的一段基因組學(xué)特征中不存在發(fā)生突變的基因，則這一段的累積特征輸出為0.相反，若存在發(fā)生突變的基因個數(shù)越多，則這一段的累積特征輸出會越大.這一結(jié)果與腫瘤的發(fā)生是由多個相關(guān)基因突變共同作用引起的思想吻合.局部累積特征提取的實質(zhì)是在多尺度特征的基礎(chǔ)上采用局部累積處理，突出基因區(qū)域(或基因段)的整體影響，而不是僅考慮個別基因位的影響，從而可進一步提高特征的整體表達性.該方法不僅能反映腫瘤細(xì)胞系的特性，還在保證基因組學(xué)特征全覆蓋的同時降低特征維度和緩解稀疏性.

1.4 基于神經(jīng)網(wǎng)絡(luò)的抗腫瘤藥物反應(yīng)預(yù)測

采用機器學(xué)習(xí)方法有效捕獲生物標(biāo)記物與抗腫瘤藥物反應(yīng)之間復(fù)雜的非線性關(guān)系是預(yù)測抗腫瘤藥物反應(yīng)的另一個關(guān)鍵問題，現(xiàn)有研究方法已取得一定成果，但也存在網(wǎng)絡(luò)結(jié)構(gòu)較復(fù)雜的情況[6,9-10].為了簡化網(wǎng)絡(luò)結(jié)構(gòu),滿足多尺度局部累積特征輸入的需求，進一步提升藥物反應(yīng)的預(yù)測性能，本文設(shè)計一套抗腫瘤藥物反應(yīng)預(yù)測方法.該系統(tǒng)為每個單尺度的局部累積特征構(gòu)建一個基于3層全連接神經(jīng)網(wǎng)絡(luò)的預(yù)測子網(wǎng)絡(luò)，并采用決策融合和特征融合，分別融合3個子網(wǎng)絡(luò)的預(yù)測輸出或特征輸出,實現(xiàn)對抗腫瘤藥物的敏感性或耐藥性進行分類預(yù)測.

1.4.1 預(yù)測子網(wǎng)絡(luò)的設(shè)計

針對腫瘤細(xì)胞系某一尺度的局部累積特征Fk，本文基于3層全連接神經(jīng)網(wǎng)絡(luò)設(shè)計1個預(yù)測子網(wǎng)絡(luò),用于預(yù)測某一種抗腫瘤藥物對Fk的藥物反應(yīng).多層全連接神經(jīng)網(wǎng)絡(luò)具有結(jié)構(gòu)靈活的特點，滿足為每種藥物設(shè)計適宜的網(wǎng)絡(luò)結(jié)構(gòu)的目的要求.本文設(shè)計的每個預(yù)測子網(wǎng)絡(luò)均采用N0-N1-N2-N3的3層全連接神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，N1>N2,其中,N0為輸入特征維度，N1、N2為隱層節(jié)點個數(shù)，N3為輸出節(jié)點個數(shù).預(yù)測子網(wǎng)絡(luò)中隱層的激活函數(shù)采用正切函數(shù)，輸出層的激活函數(shù)采用softmax.為了防止預(yù)測子網(wǎng)絡(luò)過擬合，每個隱層均使用失活層和L2正則化.

如果預(yù)測子網(wǎng)絡(luò)是針對藥物的敏感反應(yīng)進行預(yù)測，那么網(wǎng)絡(luò)的兩類輸出節(jié)點分別表示藥物作用在某一細(xì)胞系上表現(xiàn)為敏感或不敏感的可能性.對抗腫瘤藥物反應(yīng)的理想標(biāo)簽1和0進行One-Hot編碼，得到向量[1,0]T或[0,1]T作為預(yù)測子網(wǎng)絡(luò)的理想輸出.若藥物作用在某一細(xì)胞系的藥物反應(yīng)越敏感，預(yù)測其敏感性的輸出越接近于1；否則，越接近于0.兩類輸出節(jié)點之和為1，節(jié)點中取輸出最大對應(yīng)的類別作為預(yù)測子網(wǎng)絡(luò)最終的預(yù)測結(jié)果.

1.4.2 損失函數(shù)的設(shè)計

分類問題中經(jīng)常使用的損失函數(shù)是交叉熵?fù)p失函數(shù)，顯著特性是數(shù)量多的簡單樣本損失值會淹沒數(shù)量少的困難樣本損失值，使網(wǎng)絡(luò)模型的優(yōu)化方向并不是預(yù)期所想.例如，在目前已有的用于抗腫瘤藥物反應(yīng)研究的醫(yī)學(xué)數(shù)據(jù)樣本中，腫瘤細(xì)胞系的用藥反應(yīng)表現(xiàn)為敏感或耐藥的正樣本個數(shù)通常是缺少的，這樣的數(shù)據(jù)分布會造成網(wǎng)絡(luò)模型的預(yù)測性能尤其是正樣本的預(yù)測性能下降.為了解決這一問題，本文引入由Lin等[15]提出的改進的交叉熵?fù)p失函數(shù)Focal Loss，核心思想是為正負(fù)樣本分配不同的損失權(quán)重，并在優(yōu)化過程中逐漸降低簡單樣本的影響.Focal Loss優(yōu)點是不需額外增加困難樣本，便可有效訓(xùn)練所有樣本.Focal Loss的表達式為

FL(d,p)=-α(1-p)γdlog2p

-(1-α)pγ(1-d)log2(1-p),

(3)

其中,d為訓(xùn)練數(shù)據(jù)中的某一細(xì)胞系(基因序列)對某種抗腫瘤藥物表現(xiàn)為敏感(或不敏感)的理想標(biāo)簽(d=1表示敏感，d=0表示不敏感)，p為所用方法對d的預(yù)測結(jié)果(p∈[0, 1]).

由式(3)可看出，F(xiàn)ocal Loss通過最小化交叉熵?fù)p失優(yōu)化模型的訓(xùn)練過程.在交叉熵?fù)p失函數(shù)中：增加一個加權(quán)因子α，平衡正負(fù)樣本的比例；增加一個調(diào)節(jié)因子(1-p)r，調(diào)節(jié)簡單樣本損失降低的速率.α和1-α分別作為正、負(fù)樣本的損失權(quán)重.隨著γ的增加，調(diào)節(jié)因子能減小簡單樣本對損失的貢獻，擴大樣本接收低損失的范圍，增加糾正錯誤分類樣本的重要性.

1.4.3 融合與決策

不同尺度的基因組學(xué)特征包含不同的特征信息，并對最終的預(yù)測結(jié)果起到不同的作用.本文針對如何合理有效地融合不同尺度這一問題進行研究.按照特征融合和網(wǎng)絡(luò)預(yù)測的先后順序，分為早融合和晚融合.早融合是先融合多尺度特征，再在融合后的特征上訓(xùn)練分類器.晚融合是對以不同尺度的特征為輸入的預(yù)測子網(wǎng)絡(luò)的輸出結(jié)果進行融合.因此，本文提出如下兩種方法，提高單尺度特征的抗腫瘤藥物反應(yīng)預(yù)測性能，具體描述如下.

方法1為決策融合，采用加權(quán)融合方法(Weighted Fusion, WF)，下文中均采用Proposed-WF表示該方法.加權(quán)融合方法是對以不同單尺度特征為輸入的預(yù)測子網(wǎng)絡(luò)的輸出pv進行加權(quán)后，得到多尺度特征對某種抗腫瘤藥物反應(yīng)的預(yù)測結(jié)果:

最簡單的wv取值為1/V，即單尺度預(yù)測子網(wǎng)絡(luò)的輸出pv對最終藥物反應(yīng)預(yù)測結(jié)果p的貢獻是平等的.此外，wv取值還可根據(jù)每個預(yù)測子網(wǎng)絡(luò)采用訓(xùn)練集進行預(yù)測性能評估時得到的歸一化AUC(Area under Curve)值或ACC(Accuracy)值等性能指標(biāo)確定.此時，單尺度預(yù)測子網(wǎng)絡(luò)的輸出pv對最終藥物反應(yīng)預(yù)測結(jié)果p的貢獻與預(yù)測子網(wǎng)絡(luò)的綜合性能指標(biāo)呈正比關(guān)系.

2 實驗及結(jié)果分析

2.1 實驗環(huán)境

本文的抗腫瘤藥物反應(yīng)預(yù)測方法在CPU+GPU平臺上實現(xiàn).實驗平臺配置為：Intel?CoreTMi7-7700@3.6 GHz，內(nèi)存32.00 GB，NVIDIA GeForce GTX 1080Ti.操作系統(tǒng)為Windows 7，開發(fā)環(huán)境為Rstudio 1.1.453+Pycharm Edu 3.5.1+Keras框架.Apriori算法的實現(xiàn)使用R 3.4.0中的Arules 1.6.3程序包，預(yù)測模塊的實現(xiàn)使用Python3.5中的Keras 2.2.4庫，后端為TensorFlow1.4.0.

本文采用的評價指標(biāo)為：敏感性(Sensitivity, SENS)、特異性(Specificity, SPEC)、約登指數(shù)(YD’s Index)，全體樣本準(zhǔn)確率(ACC)和ROC(Receiver Operating Characteristic Curve)曲線下的面積AUC值.

為了定量計算上述性能評價指標(biāo)，根據(jù)各評價指標(biāo)在抗腫瘤藥物反應(yīng)預(yù)測問題上的生物學(xué)含義，在實驗中首先分別計算出真陽性(True Positive, TP)、假陽性(False Positive, FP)、真陰性(True Negative, TN)、假陰性(False Negative, FN)4個實驗測試統(tǒng)計值.其中，TP定義為：在對某種抗腫瘤藥物反應(yīng)進行預(yù)測的實驗中，預(yù)測為敏感的所有基因序列中其真實標(biāo)簽也為敏感的基因序列的個數(shù).FP定義為：預(yù)測為敏感的所有基因序列中其真實標(biāo)簽為不敏感的基因序列的個數(shù).TN定義為預(yù)測為不敏感的所有基因序列中其真實標(biāo)簽也為不敏感的基因序列的個數(shù).FN定義為：預(yù)測為不敏感的所有基因序列中其真實標(biāo)簽為敏感的基因序列的個數(shù).在此基礎(chǔ)上分別計算SENS、SPEC、YD’Index、ACC、AUC值：

AUC值指標(biāo)是計算ROC特性曲線與橫坐標(biāo)軸圍成的面積，是與分類閾值無關(guān),對樣本類別不平衡問題的魯棒性的評價指標(biāo).

2.2 藥物反應(yīng)預(yù)測性能分析

本文設(shè)計一套綜合性的評估方案，用于分析不同特征提取方法或不同網(wǎng)絡(luò)結(jié)構(gòu)對抗腫瘤藥物反應(yīng)預(yù)測性能的影響.方案的目的是合理評估本文方法在預(yù)測腫瘤細(xì)胞系對藥物反應(yīng)表現(xiàn)為敏感或耐藥方面的性能.

為了公平對比，本文與針對GDSC數(shù)據(jù)庫中上百種藥物進行抗腫瘤藥物反應(yīng)預(yù)測的方法進行對比分析，具體包括：基線模型DLNNs,AutoBorutaRF[5]和MOLI[7].為了與DLNNs和MOLI進行公平對比，本文采用相同的數(shù)據(jù)劃分及預(yù)處理方式，以基因表達、拷貝數(shù)變異和基因突變3種基因組學(xué)數(shù)據(jù)的狀態(tài)信息作為細(xì)胞系的特征，網(wǎng)絡(luò)輸出為分類腫瘤細(xì)胞系對抗腫瘤藥物的敏感性或耐藥性.DLNNs和MOLI的網(wǎng)絡(luò)模型是依照文獻[6]和文獻[7]中提供的程序代碼編程實現(xiàn).AutoBorutaRF是以基因表達譜、拷貝數(shù)變異的實值和基因突變的狀態(tài)信息3種基因組學(xué)數(shù)據(jù)作為細(xì)胞系的特征，模型輸出為分類腫瘤細(xì)胞系對抗腫瘤藥物的敏感性.

為了避免本文方法過學(xué)習(xí)或欠學(xué)習(xí)，實驗時本文分別采用三折交叉驗證和五折交叉驗證.根據(jù)腫瘤細(xì)胞系的組織來源進行訓(xùn)練集和測試集的劃分，這是為了保證本文方法可學(xué)習(xí)到所有組織上的腫瘤細(xì)胞特征.具體以三折交叉驗證為例，在每次實驗過程中，將來自同一組織的腫瘤細(xì)胞系隨機劃分3份，其中2份用于訓(xùn)練，1份用于測試.該過程重復(fù)3次，每次留不同的部分作為測試集，最終的預(yù)測性能指標(biāo)是3次測試的平均結(jié)果.若某種組織只包含2個細(xì)胞系，訓(xùn)練集和測試集均分.若某種組織只包含1個細(xì)胞系，分配到訓(xùn)練集中.經(jīng)上述劃分，訓(xùn)練集共包含669個細(xì)胞系樣本，測試集共包含332個細(xì)胞系樣本，比例接近2∶1.

各方法的預(yù)測性能對比如表1所示，表中黑體數(shù)字表示最優(yōu)值.AutoBorutaRF針對GDSC數(shù)據(jù)庫上98種藥物，經(jīng)十折交叉驗證后得到平均預(yù)測性能.DLNNs、MOLI及本文方法是針對GDSC數(shù)據(jù)庫上196種藥物，經(jīng)三折交叉驗證后得到平均預(yù)測性能.

表1 各方法的預(yù)測性能對比Table 1 Prediction performance comparison of different methods

由表1可知，本文方法在全部性能指標(biāo)上均有所提升.由于在GDSC數(shù)據(jù)庫上存在正負(fù)樣本不平衡問題，模型學(xué)習(xí)到正樣本的特征不足，造成敏感性指標(biāo)值較低.DLNNs的敏感性指標(biāo)值明顯低于特異性指標(biāo)值也說明這一點.為了克服正負(fù)樣本不平衡問題，AutoBorutaRF采用如下策略：將負(fù)樣本分成T份，每份負(fù)樣本Ti的樣本數(shù)近似等于正樣本R.每份Ti和R訓(xùn)練一個RF(Random Forest)分類器.根據(jù)每個RF投票獲得最終的抗腫瘤藥物反應(yīng)預(yù)測結(jié)果.由表1可看出，AutoBorutaRF的敏感性指標(biāo)值高于DLNNs.Proposed-WF和Proposed-FF采用Focal Loss損失函數(shù)解決這一問題，其中，F(xiàn)ocal Loss的參數(shù)設(shè)置為γ=2，α=0.2.Proposed-WF和Proposed-FF得到的敏感性指標(biāo)值均超過對比方法，比次優(yōu)的AutoBorutaRF分別提高0.101和0.097.Proposed-WF和Proposed-FF得到的ACC指標(biāo)值比次優(yōu)的MOLI分別提高0.056和0.059，AUC指標(biāo)值比次優(yōu)的MOLI分別提高0.138和0.136.

在GDSC數(shù)據(jù)庫的196種抗腫瘤藥物中，采用Proposed-WF和Proposed-FF進行預(yù)測，獲得AUC值最高的藥物是達拉菲尼(Dabrafenib)，分別為0.970和0.964.達拉菲尼用于治療轉(zhuǎn)移性黑色素瘤，是BRAF基因突變的抑制劑.本文采用的多尺度累積特征提取方法在挑選基因組學(xué)特征的過程中，證實達拉菲尼可能會對BRAF基因產(chǎn)生重要的影響.因此，BRAF基因也作為預(yù)測達拉菲尼藥物反應(yīng)的輸入特征.

2.3 預(yù)測性能的影響參數(shù)

2.3.1 不同基因組學(xué)特征對預(yù)測性能的影響

為了分析不同特征提取方法對抗腫瘤藥物反應(yīng)預(yù)測性能的影響，本節(jié)在藥物反應(yīng)預(yù)測模型基本相似的前提下，開展隨機挑選特征、單尺度特征、單尺度局部累積特征、多尺度特征和多尺度局部累積特征的分解實驗.實驗中隨機挑選特征是指文獻[6]的特征提取方法.3組單尺度基因組學(xué)特征X1、X2、X3是通過分別設(shè)置3條關(guān)聯(lián)規(guī)則中的最小支持度閾值和最小可信度閾值為1/669、2/669和3/669得到的.單尺度局部累積特征是指采用累積窗函數(shù)從單尺度基因組學(xué)特征Xv中提取的局部累積特征Fv.本文采用7種不同窗寬的累積窗函數(shù),分別對輸入特征X1、X2和X3提取局部累積特征，即M=1,5,10,20,30,40,50.實驗發(fā)現(xiàn)：相比其它M，單尺度特征X3在M=10下提取的局部累積特征F3的藥物反應(yīng)預(yù)測性能相對較優(yōu)，而X2和X1在M=30下提取的F2和F1的預(yù)測性能相對較優(yōu).因此，在多尺度局部累積特征融合方法中，僅采用這些預(yù)測性能相對較優(yōu)的單尺度局部累積特征進行特征融合.多尺度特征MX是指采用Proposed-FF分別將單尺度基因組學(xué)特征X1、X2、X3的N2層特征輸出后級聯(lián)得到的特征.多尺度局部累積特征MF是指采用Proposed-FF分別將單尺度局部累積特征F1、F2、F3的N2層特征輸出后級聯(lián)得到的特征.經(jīng)五折交叉驗證后，本文統(tǒng)計上述不同的特征提取方法對抗腫瘤藥物反應(yīng)預(yù)測的性能，結(jié)果如表2所示，表中黑體數(shù)字表示最優(yōu)值.

表2 不同基因組學(xué)特征的預(yù)測性能對比Table 2 Prediction performance comparison of different genomic features

由表2可知，采用以多尺度特征MX為輸入，對抗腫瘤藥物的敏感性或耐藥性進行分類預(yù)測的性能優(yōu)于采用單尺度特征X1、X2或X3.同樣，采用多尺度局部累積特征MF的預(yù)測性能也優(yōu)于采用單尺度局部累積特征F1，F(xiàn)2或F3.例如，采用多尺度局部累積特征MF進行預(yù)測，得到最優(yōu)性能的ACC值為0.835，AUC值為0.888，相比相對較優(yōu)的采用單尺度局部累積特征F1，ACC值和AUC值分別提高0.005和0.005.這說明相比單尺度特征，多尺度特征提取方法能提取更多對藥物反應(yīng)預(yù)測性能有幫助的基因組學(xué)特征，進一步提高每種單尺度特征的抗腫瘤藥物反應(yīng)預(yù)測性能.此外，采用局部累積特征對藥物反應(yīng)進行分類預(yù)測的各項性能指標(biāo)值均優(yōu)于不采用局部累積特征的.例如，相比不采用局部累積的多尺度特征MX，采用多尺度局部累積特征MF進行預(yù)測得到的ACC值和AUC值，分別提高0.065和0.122.這說明考慮多個相關(guān)基因的共同作用，在特征觀察窗內(nèi)提取基因組學(xué)局部區(qū)域的累積特征，可進一步提高基因組學(xué)特征的整體表達性.

2.3.2 不同網(wǎng)絡(luò)結(jié)構(gòu)對預(yù)測性能的影響

網(wǎng)絡(luò)4采用N0-N1-N2-N3-N4-N5-N6的6層結(jié)構(gòu),N1～N6節(jié)點數(shù)大小先增大后減小，取值分別為

經(jīng)五折交叉驗證后，本文統(tǒng)計網(wǎng)絡(luò)1～網(wǎng)絡(luò)4的抗腫瘤藥物反應(yīng)預(yù)測性能，結(jié)果如表3所示，表中黑體數(shù)字表示最優(yōu)值.

表3 不同網(wǎng)絡(luò)結(jié)構(gòu)的預(yù)測性能對比Table 3 Prediction performance comparison of different network structures

由表3可看出，雖然隨著神經(jīng)網(wǎng)絡(luò)的層數(shù)增加，藥物反應(yīng)預(yù)測性能在特異性指標(biāo)上有一定程度的提升，但敏感性指標(biāo)卻下降明顯.網(wǎng)絡(luò)4的特異性指標(biāo)為0.896，相比網(wǎng)絡(luò)1，提升0.054；敏感性指標(biāo)為0.650，相比網(wǎng)絡(luò)1，下降0.115.此外，由于醫(yī)學(xué)數(shù)據(jù)的正負(fù)樣本不平衡問題，當(dāng)樣本量占比較大的負(fù)樣本的正確率提升時，全體樣本的準(zhǔn)確率也隨之提升.采用網(wǎng)絡(luò)4的ACC值為0.864，相比網(wǎng)絡(luò)1，提升0.029.但此時另一關(guān)鍵的綜合性能指標(biāo)AUC值卻隨著神經(jīng)網(wǎng)絡(luò)的層數(shù)增加而有不同程度的下降.例如，本文采用的3層的神經(jīng)網(wǎng)絡(luò)1獲得的AUC值為0.888，相比采用6層的神經(jīng)網(wǎng)絡(luò)4，提升0.013.從網(wǎng)絡(luò)模型的參數(shù)量角度分析，網(wǎng)絡(luò)4的參數(shù)量為

網(wǎng)絡(luò)1的參數(shù)量為

例如，當(dāng)N0=100時，網(wǎng)絡(luò)4的參數(shù)量為3 122，網(wǎng)絡(luò)1的參數(shù)量為872.這說明采用更多層數(shù)的神經(jīng)網(wǎng)絡(luò)需要更多的計算代價，卻未得到更好的預(yù)測性能.在綜合權(quán)衡考慮下，本文決定采用網(wǎng)絡(luò)1的3層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu).

2.4 實用性價值評估

有研究指出：AUC值在0.5～0.7時，模型的預(yù)測準(zhǔn)確性較低.AUC值在0.7～0.9時，模型具有較好的預(yù)測準(zhǔn)確性.AUC值大于0.9時，模型具有較高的預(yù)測準(zhǔn)確性[16].AutoBorutaRF也采用AUC≥0.7作為合格的模型預(yù)測性能評判標(biāo)準(zhǔn).本文統(tǒng)計當(dāng)前研究方法中AUC≥0.7的比例，與Proposed-WF和Proposed-FF進行對比.AutoBorutaRF中AUC≥0.7的比例是參考該文給出的實驗結(jié)果，統(tǒng)計GDSC數(shù)據(jù)庫上98種藥物獲得的，而DLNNs和MOLI是統(tǒng)計GDSC數(shù)據(jù)庫上196種藥物獲得的.經(jīng)三折交叉驗證后，本文統(tǒng)計上述不同方法獲得預(yù)測性能AUC≥0.7的比例，具體如下:AutoBorutaRF為60.2%，MOLI為70.4%，DLNNs為72.0%，Proposed-WF為100%，Proposed-FF為100%.由仿真結(jié)果可看出，本文方法所有抗腫瘤藥物反應(yīng)預(yù)測模型的AUC值均大于0.7，這表明基于多尺度局部累積特征和多層全連接神經(jīng)網(wǎng)絡(luò)的抗腫瘤藥物反應(yīng)預(yù)測模型具有較優(yōu)的預(yù)測準(zhǔn)確性，均是合格的預(yù)測模型.以Proposed-FF為例，具體AUC值比例分布如下：AUC值在0.7～0.8的比例為0.3%，AUC值在0.8～0.9的比例為76.4%，AUC值在0.9～1的比例為23.3%.

抗腫瘤藥物反應(yīng)預(yù)測的目的是要面向臨床應(yīng)用，針對不同患者,通過該模型提供個性化治療的藥物推薦.在預(yù)測準(zhǔn)確性方面，全球疾病誤診率高達30%，即ACC≥0.7.由于腫瘤存在異質(zhì)性，醫(yī)生的評估準(zhǔn)確率在實際中更低.經(jīng)三折交叉驗證后，本文方法與AutoBorutaRF、DLNNs、MOLI進行對比，各方法獲得預(yù)測性能ACC≥0.7的比例如下：DLNNs為73.5%，AutoBorutaRF為84.8%，MOLI為86.5%，Proposed-WF為100%，Proposed-FF為100%.Auto-BorutaRF中ACC≥0.7的比例是參考該文給出的實驗結(jié)果，統(tǒng)計GDSC數(shù)據(jù)庫上98種藥物獲得的，而DLNNs和MOLI是統(tǒng)計GDSC數(shù)據(jù)庫上196種藥物獲得的.

本文方法所有抗腫瘤藥物反應(yīng)預(yù)測模型的ACC值均大于0.7，具體ACC值比例分布如下：ACC值在0.7～0.8的比例為17.99%，ACC值在0.8～0.9的比例為80.16%，ACC值在0.9～1的比例為1.85%.因此，本文方法有利于確定藥物與腫瘤細(xì)胞系基因組學(xué)特征之間的關(guān)系，有望在個性化醫(yī)療領(lǐng)域發(fā)揮重要作用.

3 結(jié) 束 語

本文結(jié)合現(xiàn)有公開發(fā)布的人類腫瘤細(xì)胞系的基因組學(xué)數(shù)據(jù)，提出基于多尺度局部累積特征和神經(jīng)網(wǎng)絡(luò)的抗腫瘤藥物反應(yīng)預(yù)測方法.采用多尺度關(guān)聯(lián)規(guī)則挖掘方法提取基因組學(xué)數(shù)據(jù)的多尺度特征，避免由于采用單一尺度特征提取時可能會遺漏某些有用信息而存在的問題.提出多尺度局部累積特征提取方法，突出基因區(qū)域(或基因段)的整體影響，而不是僅考慮個別基因位的影響，進一步提高特征的整體表達性.提出一套由3個3層子網(wǎng)絡(luò)，并且分別通過采用決策融合或特征融合方法進行融合的預(yù)測模型，結(jié)構(gòu)簡單、計算量較小、性能穩(wěn)定、泛化性較好.本文工作是在COSMIC、GDSC數(shù)據(jù)庫提供的1 001株細(xì)胞系的基因組學(xué)數(shù)據(jù)和196種抗腫瘤藥物反應(yīng)數(shù)據(jù)的基礎(chǔ)上進行的.大量實驗表明，本文方法在敏感性、特異性、準(zhǔn)確率及AUC值等關(guān)鍵性能指標(biāo)上均較優(yōu).

隨著基因組學(xué)研究的發(fā)展，還將有更多的基因組數(shù)據(jù)出現(xiàn)，這為進一步采用深度卷積神經(jīng)網(wǎng)絡(luò)模型進行抗腫瘤藥物反應(yīng)預(yù)測提供研究基礎(chǔ)，也是作者下一步的主要研究方向之一.此外，文獻[17]的研究表明沉默突變也有可能會導(dǎo)致人類癌癥，因此如何利用沉默突變信息進行抗腫瘤藥物反應(yīng)的預(yù)測也是一個值得研究的方向.