基于XGBoost對肺鱗癌和肺腺癌的分類預測

冷 菲 李 巍

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

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

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

本研究應用機器學習分類器，極限梯度增強算法(extreme gradient boosting，XGBoost)，通過RNA表達譜區分肺鱗癌和肺腺癌患者。選擇這種方法是因為它具有顯著的優點，包括：可以處理缺失值，需要數據縮放，提示梯度增強算法中的有效變異，在各種競賽中取得突出的成績，并已成功應用于其他研究領域。通過使用XGBoost區分肺鱗癌和肺腺癌，并且識別能夠區分兩者的主要RNA分子。目前，沒有研究使用XGBoost根據RNA特征客觀的對兩種亞型進行分類。

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

1 資料與方法

1.1 數據來源

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

1.2 數據集構建

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

1.3 數據預處理

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

1.4 特征處理

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

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

1.5 模型算法

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

1.6 統計分析與結果評價

采用Python 3.7.1進行統計學分析。為了更加準確和全面的評估模型，采用準確率、曲線下面積(area under curve，AUC)和洛倫茲曲線(kolmogorov-smirnov，KS)作為預測結果的評價標準。

1.7 模型訓練和驗證

根據以上算法，建立分類模型。將處理后的數據進行多次訓練，對不同參數進行多次調整，最終得到最佳結果。

2 結果

2.1 特征處理

從TCGA數據庫中下載肺鱗癌和肺腺癌臨床數據和mRNA表達數據，剔除正常組織樣本和信息缺失樣本，得到474例肺腺癌樣本和491例肺鱗癌樣本。分別下載474例肺腺癌腫瘤組織和491例肺鱗癌腫瘤組織中mRNA的表達量數據。基于本研究特征處理標準，得到1 099個mRNA在肺鱗癌和肺腺癌中差異有統計學意義的表達。圖1通過火山圖顯示了-log(FDR)和logFC兩個維度上所有差異表達的mRNA的分布。將所有mRNA表達水平標準化為樣品平均值。剔除表達差異無統計學意義的基因，選擇差異表達的基因作為模型的特征子集，建立疾病亞型預測模型。

圖1 基因差異表達分析結果Fig.1 Gene differential expression analysis results

2.2 預測結果及模型性能比較

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

表1 預測結果統計Tab.1 Statistics of the prediction results

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

表2 3種預測模型預測結果比較Tab.2 Comparison of prediction results of three prediction models

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

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

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

2.3 肺癌亞型的分子標記因子

根據特征重要性，對1 099個mRNA進行打分，權重越大代表特征對預測結果貢獻越大。統計所有特征的重要性分數的分布情況(圖4)，發現932個mRNA特征的重要性分數為0，即大多數mRNA對疾病亞型的分類沒有影響。

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

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

圖5 不同特征子集的預測準確率Fig.5 Prediction accuracy of different feature subsets

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

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

3 討論

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

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

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