基于XGBoost對肺鱗癌和肺腺癌的分類預(yù)測

冷 菲 李 巍

(國家兒童醫(yī)學(xué)中心 首都醫(yī)科大學(xué)附屬北京兒童醫(yī)院 遺傳與出生缺陷防治中心 北京市兒科研究所 出生缺陷遺傳學(xué)研究北京市重點(diǎn)實(shí)驗(yàn)室 兒科重大疾病研究教育部重點(diǎn)實(shí)驗(yàn)室，北京100045)

肺癌是最常見的惡性腫瘤之一，2012年全球新發(fā)病例約182萬例，病死率約為159萬[1]。盡管近年來已經(jīng)取得了治療方法的進(jìn)步，包括微創(chuàng)手術(shù)方法，化學(xué)療法和靶向治療，但肺癌患者的5年生存率大多數(shù)情況下僅為10%～20%[2]。大多數(shù)肺癌由兩種主要病理亞型構(gòu)成：肺鱗狀細(xì)胞癌和肺腺癌，二者在臨床表現(xiàn)上非常相似，但是在發(fā)病機(jī)制以及治療和預(yù)后方面都有明顯不同[3]。

肺鱗癌和肺腺癌來自不同的細(xì)胞，不僅在生物模式方面，而且在分子特征方面，以及最重要的治療策略方面均存在差異。例如，激活表皮生長因子受體的突變和間變性淋巴瘤激酶(anaplastic lymphoma kinase，ALK)融合蛋白的突變通常發(fā)生在肺腺癌中，而不是肺鱗癌中，使針對這些基因的藥物對肺鱗癌無效[4]。因此，利用這兩種不同肺癌亞型的分子特征和機(jī)制的差異預(yù)測亞型分類，可以有針對性地對不同亞型疾病進(jìn)行干預(yù)，同時，通過深入研究兩種亞型的主要差異特征，將有助于更深入地理解和鑒定新的肺癌治療分子靶向策略。

本研究應(yīng)用機(jī)器學(xué)習(xí)分類器，極限梯度增強(qiáng)算法(extreme gradient boosting，XGBoost)，通過RNA表達(dá)譜區(qū)分肺鱗癌和肺腺癌患者。選擇這種方法是因?yàn)樗哂酗@著的優(yōu)點(diǎn)，包括：可以處理缺失值，需要數(shù)據(jù)縮放，提示梯度增強(qiáng)算法中的有效變異，在各種競賽中取得突出的成績，并已成功應(yīng)用于其他研究領(lǐng)域。通過使用XGBoost區(qū)分肺鱗癌和肺腺癌，并且識別能夠區(qū)分兩者的主要RNA分子。目前，沒有研究使用XGBoost根據(jù)RNA特征客觀的對兩種亞型進(jìn)行分類。

本研究利用癌癥基因組圖譜(The Cancer Genome Atlas，TCGA)數(shù)據(jù)庫發(fā)布的高通量數(shù)據(jù)[5]，構(gòu)建了包含474例樣本的肺鱗癌和491例樣本的肺腺癌數(shù)據(jù)庫，包含每個樣本的RNA表達(dá)數(shù)據(jù)。接下來，鑒定了肺鱗癌和肺腺癌的差異表達(dá)RNA。之后，基于XGBoost模型開發(fā)肺鱗癌和肺腺癌分類預(yù)測模型，并對預(yù)測中的重要特征進(jìn)行分析。最后基于預(yù)測中特征重要性，分析了顯著影響肺鱗癌和肺腺癌分類預(yù)測的RNA分子。

1 資料與方法

1.1 數(shù)據(jù)來源

從TCGA數(shù)據(jù)門戶檢索肺鱗癌樣本和肺腺癌樣本的RNA序列數(shù)據(jù)。 TCGA數(shù)據(jù)集(https://portal.gdc.cancer.gov/)由超過2 PB的基因組數(shù)據(jù)組成，可公開獲得，這種基因組信息有助于癌癥研究界改進(jìn)預(yù)防，診斷和治療癌癥。本研究符合TCGA提供的出版指南。由于數(shù)據(jù)來自TCGA數(shù)據(jù)庫，因此道德委員會無須進(jìn)一步批準(zhǔn)。

1.2 數(shù)據(jù)集構(gòu)建

從TCGA數(shù)據(jù)門戶下載肺鱗癌和肺腺癌的RNA-seq數(shù)據(jù)，數(shù)據(jù)來源于IlluminaHiSeq_RNASeq測序平臺；所有數(shù)據(jù)都可以免費(fèi)下載。分別去除兩個數(shù)據(jù)集中的正常樣品，共得到474例肺鱗癌樣品和491例肺腺癌腫瘤樣本。去除僅在肺鱗癌樣本出現(xiàn)或僅在肺腺癌樣本中測得的mRNA，得到最終mRNA數(shù)據(jù)集。在研究中，主要使用Python和R語言編寫的程序代碼來分析和處理RNA數(shù)據(jù)。

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

使用Ensembl數(shù)據(jù)庫(http://www.ensembl.org/index.html，版本95)鑒定了mRNA[6]，排除了未包括在數(shù)據(jù)庫中的mRNA。之后，刪除平均讀數(shù)≤1的mRNA，以過濾掉所有未表達(dá)的mRNA。

1.4 特征處理

由于mRNA數(shù)據(jù)種類太多，作為特征維度太大，需進(jìn)行特征處理，即從特征中選擇出特征子集。為了便于后期對分類影響的評價(jià)，將不對前期數(shù)據(jù)進(jìn)行特征變換處理。本文基于基因表達(dá)量在亞型之間的差異度對特征進(jìn)行篩選。

使用EdgeR對兩種疾病亞型的mRNA數(shù)據(jù)進(jìn)行表達(dá)差異分析，采用M值修剪均值(trimmed mean of M-values，TMM)算法進(jìn)行標(biāo)準(zhǔn)化，并獲得差異表達(dá)的mRNA。所有P值使用錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)來校正多重測試的差異有統(tǒng)計(jì)學(xué)意義。對于倍數(shù)變化(log2絕對值)>2且FDR為P<0.001的表達(dá)差異被認(rèn)為是顯著的。使用R平臺中的ggplots和熱圖包生成火山圖。刪除差異無統(tǒng)計(jì)學(xué)意義表達(dá)的mRNA，將差異有統(tǒng)計(jì)學(xué)意義表達(dá)的mRNA作為模型的特征集合。

1.5 模型算法

基于XGBoost算法建立模型。不同于傳統(tǒng)集成決策樹算法，XGBoost在損失函數(shù)里加入了正則項(xiàng)，在控制模型的復(fù)雜度的同時，可以防止模型過擬合。即得到目標(biāo)函數(shù)為：

1.6 統(tǒng)計(jì)分析與結(jié)果評價(jià)

采用Python 3.7.1進(jìn)行統(tǒng)計(jì)學(xué)分析。為了更加準(zhǔn)確和全面的評估模型，采用準(zhǔn)確率、曲線下面積(area under curve，AUC)和洛倫茲曲線(kolmogorov-smirnov，KS)作為預(yù)測結(jié)果的評價(jià)標(biāo)準(zhǔn)。

1.7 模型訓(xùn)練和驗(yàn)證

根據(jù)以上算法，建立分類模型。將處理后的數(shù)據(jù)進(jìn)行多次訓(xùn)練，對不同參數(shù)進(jìn)行多次調(diào)整，最終得到最佳結(jié)果。

2 結(jié)果

2.1 特征處理

從TCGA數(shù)據(jù)庫中下載肺鱗癌和肺腺癌臨床數(shù)據(jù)和mRNA表達(dá)數(shù)據(jù)，剔除正常組織樣本和信息缺失樣本，得到474例肺腺癌樣本和491例肺鱗癌樣本。分別下載474例肺腺癌腫瘤組織和491例肺鱗癌腫瘤組織中mRNA的表達(dá)量數(shù)據(jù)?；诒狙芯刻卣魈幚順?biāo)準(zhǔn)，得到1 099個mRNA在肺鱗癌和肺腺癌中差異有統(tǒng)計(jì)學(xué)意義的表達(dá)。圖1通過火山圖顯示了-log(FDR)和logFC兩個維度上所有差異表達(dá)的mRNA的分布。將所有mRNA表達(dá)水平標(biāo)準(zhǔn)化為樣品平均值。剔除表達(dá)差異無統(tǒng)計(jì)學(xué)意義的基因，選擇差異表達(dá)的基因作為模型的特征子集，建立疾病亞型預(yù)測模型。

圖1 基因差異表達(dá)分析結(jié)果Fig.1 Gene differential expression analysis results

2.2 預(yù)測結(jié)果及模型性能比較

本次實(shí)驗(yàn)在Linux環(huán)境下進(jìn)行，選擇 Python 語言進(jìn)行實(shí)現(xiàn)，將1 099個差異表達(dá)的mRNA作為特征，基于XGBoost進(jìn)行建模，預(yù)測結(jié)果與數(shù)據(jù)庫中真實(shí)數(shù)據(jù)進(jìn)行比較，最終得出評判結(jié)果。計(jì)算分類預(yù)測的準(zhǔn)確率為96.55%，AUC值為99.04%，Kappa值為0.92，肺鱗癌陽性預(yù)測值為0.97，肺腺癌陽性預(yù)測值為0.96(表1)。

表1 預(yù)測結(jié)果統(tǒng)計(jì)Tab.1 Statistics of the prediction results

為了比較模型性能，除了XGBoost模型之外，本研究同時還采用了邏輯回歸算法和支持向量機(jī)算法進(jìn)行建模和預(yù)測。對3種模型分別進(jìn)行參數(shù)優(yōu)化，利用相同方法劃分訓(xùn)練集和測試集，分別對肺鱗癌和肺腺癌的分類進(jìn)行預(yù)測，3種模型最終準(zhǔn)確率分別為96.55%，89.97%和89.34%(表2)，XGBoost模型預(yù)測準(zhǔn)確率明顯高于邏輯回歸模型(χ2=37.66，P<0.001)和支持向量機(jī)模型(χ2=49.8，P<0.001)，前者比后兩者的準(zhǔn)確率分別提高近6.58和7.21個百分點(diǎn)。 3種模型的預(yù)測AUC結(jié)果見圖2。XGBoost模型的受試者操作特征(receiver operating characteristic，ROC)曲線，分別高于邏輯回歸模型和支持向量機(jī)模型近4個百分點(diǎn)。

表2 3種預(yù)測模型預(yù)測結(jié)果比較Tab.2 Comparison of prediction results of three prediction models

圖2 不同模型預(yù)測曲線下面積結(jié)果比較Fig.2 Area under curve comparison for different prediction models

為了評估模型將正、負(fù)樣本區(qū)分開的程度，計(jì)算KS值(圖3)XGBoost模型的KS值明顯高于logistic regression模型和SVM模型，基于mRNA的肺鱗癌和肺腺癌的分類，XGBoost模型的預(yù)測性能優(yōu)于另外2種模型。

圖3 不同模型預(yù)測洛倫茲曲線結(jié)果比較Fig.3 Kolmogorov-Smirnov comparison of different prediction models

2.3 肺癌亞型的分子標(biāo)記因子

根據(jù)特征重要性，對1 099個mRNA進(jìn)行打分，權(quán)重越大代表特征對預(yù)測結(jié)果貢獻(xiàn)越大。統(tǒng)計(jì)所有特征的重要性分?jǐn)?shù)的分布情況(圖4)，發(fā)現(xiàn)932個mRNA特征的重要性分?jǐn)?shù)為0，即大多數(shù)mRNA對疾病亞型的分類沒有影響。

圖4 基因特征重要性打分的分布情況Fig.4 Distribution of feature scores

為了進(jìn)一步篩選特征，根據(jù)特征的權(quán)重對特征從大到小排序，依次累加特征，形成逐漸增大的特征子集。計(jì)算不同大小特征子集下預(yù)測結(jié)果的準(zhǔn)確率，即這些特征對預(yù)測結(jié)果的累計(jì)貢獻(xiàn)率，結(jié)果如圖5。僅使用貢獻(xiàn)率最大的mRNA進(jìn)行預(yù)測，準(zhǔn)確率已經(jīng)達(dá)到91.5%。隨著特征數(shù)目的增加，準(zhǔn)確率逐漸升高，當(dāng)特征數(shù)為11個時，準(zhǔn)確率達(dá)到95.9%。之后，準(zhǔn)確率的數(shù)值趨于穩(wěn)定，不再根據(jù)特征數(shù)的增加而顯著增長，甚至偶爾存在輕微下降的情況。

圖5 不同特征子集的預(yù)測準(zhǔn)確率Fig.5 Prediction accuracy of different feature subsets

因此，選擇F-score排名前11的mRNA作為區(qū)分肺鱗癌和肺腺癌的分子標(biāo)記因子，這11個mRNA的基因名稱和重要性評估結(jié)果見表3。

表3 排名前11的基因名稱及重要性得分Tab.3 Top 11 gene names and their importance scores

3 討論

肺鱗癌和肺腺癌是最常見的兩種非小細(xì)胞肺癌，二者臨床表現(xiàn)相似，但致病機(jī)制和治療預(yù)后不同，尤其存活率差異很大，而目前尚無有效的早期診斷方法。本研究基于轉(zhuǎn)錄組數(shù)據(jù)深入研究兩種肺癌亞型的分子特征，首先通過差異分析找到1 099個在肺鱗癌和肺腺癌中表達(dá)差異有統(tǒng)計(jì)學(xué)意義的基因，將它們作為模型特征，之后結(jié)合XGBoost建立模型，對癌癥亞型進(jìn)行預(yù)測。準(zhǔn)確率達(dá)到96%以上?；谧罱K建立的模型，對特征重要性進(jìn)行評估。此外，還比較了XGBoost模型和其他模型的預(yù)測結(jié)果。選擇了分類問題中比較經(jīng)典的邏輯回歸算法和支持向量機(jī)算法建立模型。

邏輯回歸使用Sigmoid函數(shù)，將線性模型的結(jié)果壓縮到0～1，使其擁有概率意義。支持向量機(jī)則是在特征空間中尋找使正類負(fù)類間隔最大的超平面的線性分類器。根據(jù)結(jié)果可知，這兩種模型的準(zhǔn)確率均較XGBoost偏低。由于XGBoost模型用到了損失函數(shù)的二階泰勒展開，因此與損失函數(shù)更接近，收斂更快。并且在損失函數(shù)中加入了正則項(xiàng)，可以有效控制模型復(fù)雜度，防止模型過擬合。根據(jù)特征重要性評估，找到11個基因作為區(qū)分肺鱗癌和肺腺癌的分子標(biāo)記因子，其中貢獻(xiàn)前3名的基因依次是MACC1、KRT5和SPRR2E。MACC1蛋白定位于人類的7號常染色體上(7P21.1)，具有廣泛的生物學(xué)功能，特別是在調(diào)控惡性腫瘤的侵襲和轉(zhuǎn)移等方面具有不可代替的重要功能[7]。MACC1 在多種惡性腫瘤如結(jié)腸癌、卵巢癌、肺癌、肝癌、胃癌等組織中表達(dá)異常增高，與腫瘤臨床分期、有無遠(yuǎn)處轉(zhuǎn)移密切相關(guān)，有作為腫瘤轉(zhuǎn)移和預(yù)后判斷的潛在獨(dú)立指標(biāo)[8]。Keratin 5，也稱為KRT5，K5或CK5，是由KRT5基因在人體內(nèi)編碼的蛋白質(zhì)。這種蛋白質(zhì)涉及多種疾病，包括大皰性表皮松解癥和乳腺癌和肺癌[9]。SPRR2E是染色體1q21上的人表皮分化復(fù)合物的一部分，在一項(xiàng)關(guān)于表皮鱗狀細(xì)胞癌的研究[10]中，與正常相比，惡性角質(zhì)形成細(xì)胞系中SPRR2E呈現(xiàn)低表達(dá)，提示終末分化缺陷，而這正是致癌轉(zhuǎn)化的特征。以上研究[9-10]表明SPRR2E主要功能與皮膚角質(zhì)形成細(xì)胞的分化密切相關(guān)，但與肺癌發(fā)生的作用尚不明確。有研究[11]顯示支氣管上皮細(xì)胞經(jīng)氡染毒后，會引起SPRR2E表達(dá)量的顯著變化。

綜上所述，本文建立了預(yù)測肺鱗癌和肺腺癌的模型，具有較高的準(zhǔn)確性和良好的穩(wěn)定性，為其早期診斷和治療提供理論依據(jù)。同時找到一些特征變量作為區(qū)分肺鱗癌和肺腺癌的分子標(biāo)記，后續(xù)研究中將會進(jìn)一步探索這些特征變量對分類的影響機(jī)制。