王英倫，何孜灝，黎思穎，鄭子豪，董博文，鄧展程，陳秋銳，沈晗，2

（1.廣東藥科大學(xué)生命科學(xué)與生物制藥學(xué)院，廣東 廣州 510006；2.廣東省生物技術(shù)候選藥物研究重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510006）

黑色素瘤是來(lái)源于黑色素細(xì)胞的一類高度惡性腫瘤，多發(fā)于皮膚，早期發(fā)現(xiàn)并手術(shù)切除原發(fā)性皮膚黑色素瘤具有較好療效，然而其發(fā)病具有隱匿性強(qiáng)、易轉(zhuǎn)移等特點(diǎn)，多數(shù)患者在晚期發(fā)生轉(zhuǎn)移后才被診斷，預(yù)后較差[1]。同時(shí)，原發(fā)黑色素瘤與轉(zhuǎn)移黑色素瘤在治療手段上也存在一定差異，因此，及時(shí)準(zhǔn)確地判斷黑色素瘤的轉(zhuǎn)移特征對(duì)臨床治療具有重要的指導(dǎo)意義。近年來(lái)，隨著二代高通量測(cè)序技術(shù)的發(fā)展和應(yīng)用，已積累了大量腫瘤相關(guān)的測(cè)序數(shù)據(jù)，其中以TCGA（The Cancer Genome Atlas）數(shù)據(jù)庫(kù)和GEO（Gene Expression Ommibus）數(shù)據(jù)庫(kù)為代表，這些數(shù)據(jù)信息為腫瘤相關(guān)研究提供了極大便利。然而，面對(duì)海量的生物信息數(shù)據(jù)，傳統(tǒng)研究方法已無(wú)法適應(yīng)要求，AI（artificial intelligence，人工智能）技術(shù)的發(fā)展讓科學(xué)家可以借用電腦實(shí)現(xiàn)數(shù)據(jù)的高效快速分析，機(jī)器學(xué)習(xí)是AI 中最具智能特征、最前沿的研究領(lǐng)域之一[2]。本文選取TCGA 數(shù)據(jù)庫(kù)及GEO數(shù)據(jù)庫(kù)中皮膚黑色素瘤的轉(zhuǎn)錄組數(shù)據(jù)及臨床特征信息，運(yùn)用機(jī)器學(xué)習(xí)XGBoost（eXtreme Gradient Boosting）算法，構(gòu)建皮膚惡性黑色素瘤轉(zhuǎn)移特征預(yù)測(cè)模型，旨在為幫助判斷皮膚黑色素瘤患者的臨床分型及分期提供有益的借鑒。

1 方法

1.1 數(shù)據(jù)下載與整理

從TCGA 數(shù)據(jù)庫(kù)網(wǎng)站（https://portal.gdc.cancer.gov/）下載470 例皮膚惡性黑色素瘤樣本數(shù)據(jù)，包含原始RNA-Seq Count 數(shù)據(jù)、RNA 表達(dá)FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）數(shù)據(jù)和臨床信息數(shù)據(jù)，其中103 例為原發(fā)性腫瘤（TP）性狀，367 例為轉(zhuǎn)移性腫瘤（TM）性狀，將FPKM 數(shù)據(jù)整理轉(zhuǎn)換成TPM（Transcripts Per Kilobase of exon model per Million mapped reads）數(shù)據(jù)。從GEO 數(shù)據(jù)庫(kù)網(wǎng)站（https://www.ncbi.nlm.nih.gov/gds/）下 載GSE190113 數(shù) 據(jù)集，其包含103 例皮膚惡性黑色素瘤樣本的RNA 表達(dá)TPM 數(shù)據(jù)和臨床信息數(shù)據(jù)，所有數(shù)據(jù)分析前進(jìn)行標(biāo)準(zhǔn)化處理。

1.2 差異表達(dá)基因分析

分別利用R 語(yǔ)言（version 4.0.5）中的edgeR 包[3]與DESeq2 包[4]，以TCGA 原始RNA-Seq Count 數(shù)據(jù)制成的基因表達(dá)矩陣和對(duì)應(yīng)的臨床信息數(shù)據(jù)作為輸入資料，進(jìn)行兩次差異分析，完畢后，將兩者所得結(jié)果以P<0.05 且|log2FoldChange|>1 為篩選閾值，獲取表達(dá)差異顯著基因。并將兩種算法獲得的差異表達(dá)基因取交集。為驗(yàn)證結(jié)果，使用交集差異表達(dá)基因的TPM 數(shù)據(jù)，使用R 語(yǔ)言中的FactoMineR包[5]進(jìn)行了主成分分析（principal component analysis,PCA）。

1.3 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（WGCNA分析）

WGCNA 分析旨在找出協(xié)同表達(dá)的基因模塊，探索基因網(wǎng)絡(luò)與轉(zhuǎn)移特征表型之間的關(guān)聯(lián)，從而篩選出相關(guān)程度高的基因模塊以及網(wǎng)絡(luò)中的核心基因。利用R 語(yǔ)言中的WGCNA 包[6]，為加快運(yùn)算速度同時(shí)減少噪音，我們計(jì)算了基因表達(dá)量在樣本內(nèi)的方差，選取方差排名前25%的基因制成表達(dá)矩陣進(jìn)行層析聚類，剔除明顯離群樣本。導(dǎo)入篩選過(guò)后的基因表達(dá)矩陣與轉(zhuǎn)移相關(guān)性狀（TM 與TP），選取恰當(dāng)?shù)摩轮底鳛檐涢撝担瑯?gòu)建拓?fù)渲丿B矩陣，實(shí)現(xiàn)無(wú)標(biāo)度網(wǎng)絡(luò)。隨后進(jìn)行動(dòng)態(tài)樹剪切，將相似模塊融合，完成WGCNA分析。

1.4 生物學(xué)功能和通路富集分析

利用R 語(yǔ)言中的clusterProfiler 包[7]，對(duì)轉(zhuǎn)移特征相關(guān)模塊中的基因進(jìn)行GO 富集分析（生物學(xué)進(jìn)程BP、細(xì)胞成分CC、分子功能MF），并利用GSEA軟件[8]（version 4.1.0）對(duì)全部基因進(jìn)行GSEA 富集分析驗(yàn)證GO分析結(jié)果。

1.5 ROC分析

利用R 語(yǔ)言中的pROC 包[9]，以TP 和TM 為特征，對(duì)全部基因進(jìn)行批量化ROC 分析，繪制每個(gè)基因的ROC 曲線并計(jì)算該曲線下的面積AUC 值，并篩選AUC>0.7的基因。

1.6 機(jī)器學(xué)習(xí)XGBoost分類模型的構(gòu)建

對(duì)差異分析中的差異表達(dá)基因、WGCNA 分析中轉(zhuǎn)移特征相關(guān)模塊內(nèi)的基因和ROC 分析中AUC值大于0.8 的基因進(jìn)行交集分析，篩選交集作為訓(xùn)練模型中使用的特征。利用R語(yǔ)言中的mlr包[10]，選擇XGBoost（eXtreme Gradient Boosting，極限梯度boosting）算法并選用上述特征，構(gòu)建出一個(gè)可區(qū)分轉(zhuǎn)移性或非轉(zhuǎn)移性惡性黑色素瘤的機(jī)器學(xué)習(xí)模型，具體方法如下：將470 個(gè)樣本隨機(jī)抽樣分為兩部分?jǐn)?shù)據(jù)集，4/5 用作訓(xùn)練集，1/5 用作內(nèi)部測(cè)試集（所有的分割數(shù)據(jù)集步驟中均確保兩種類型的樣本的抽取比例與原始比例相等）；輸入訓(xùn)練集，使用XG‐Boost 算法進(jìn)行五折交叉驗(yàn)證、隨機(jī)搜索1000 次的超參數(shù)調(diào)優(yōu)獲取最佳超參數(shù)組合,以MMCE（mean misclassification error，平均誤分類誤差）作為評(píng)估指標(biāo)；將最優(yōu)超參數(shù)組合應(yīng)用于XGBoost 算法，并用其來(lái)學(xué)習(xí)訓(xùn)練集數(shù)據(jù)，獲得機(jī)器學(xué)習(xí)分類模型；為驗(yàn)證包括超參數(shù)調(diào)優(yōu)步驟在內(nèi)的整個(gè)模型構(gòu)建過(guò)程，評(píng)估該算法在訓(xùn)練集上的性能，采用嵌套交叉驗(yàn)證，其外循環(huán)為五折交叉驗(yàn)證，內(nèi)循環(huán)為超參數(shù)調(diào)優(yōu)過(guò)程。至此，獲得以惡性黑色素瘤中基因表達(dá)量為特征判別非轉(zhuǎn)移性或轉(zhuǎn)移性的模型。最后，使用該模型預(yù)測(cè)內(nèi)部測(cè)試集中以測(cè)試其泛化能力，并進(jìn)一步選取GEO 數(shù)據(jù)庫(kù)中的GSE190113 數(shù)據(jù)集作為外部測(cè)試集，驗(yàn)證該模型預(yù)測(cè)能力。

2 結(jié)果

2.1 TP組與TM組差異表達(dá)基因篩選

將TCGA 數(shù)據(jù)庫(kù)來(lái)源的皮膚惡性黑色素瘤樣本按轉(zhuǎn)移特征分為原發(fā)腫瘤組（TP 組）和轉(zhuǎn)移腫瘤組（TM 組），使用R 語(yǔ)言的edgeR 包與DESeq2 包進(jìn)行差異分析，分別獲得2041 個(gè)與3846 個(gè)差異表達(dá)基因，將兩者結(jié)果取交集后獲得1804 個(gè)差異表達(dá)基因。如圖1所示為edgeR包分析火山圖，如圖2所示為DESeq2 包分析火山圖，以TP 組作為對(duì)照組，藍(lán)色點(diǎn)表示下調(diào)基因，紅色點(diǎn)表示上調(diào)基因，灰色點(diǎn)則表示表達(dá)差異不顯著的基因。交集后差異表達(dá)基因的PCA 分析結(jié)果顯示TP 組和TM 組存在差異表達(dá)現(xiàn)象（圖2）。

圖1 edgeR包分析差異基因火山圖Figure1 Differential gene volcano plots analyzed with the edgeR package

圖2 DESeq2 R包分析差異基因火山圖Figure 2 Differential gene volcano plots analyzed with the DESeq2 R package

2.2 WGCNA分析

WGCNA 分析結(jié)果顯示，發(fā)現(xiàn)35 個(gè)轉(zhuǎn)移性狀相關(guān)模塊（圖4），其中black 模塊與TM 性狀的負(fù)相關(guān)性最高，r=-0.46，同時(shí)black 模塊與TP 性狀的正相關(guān)性最高，r=0.46，該模塊包含529個(gè)基因，在TM 組中呈下調(diào)表達(dá)。

圖3 交集差異表達(dá)基因的PCA分析Figure3 PCA analysis of intersection differentially expressed genes

圖4 模塊與轉(zhuǎn)移性狀的WGCNA分析Figure4 WGCNA analysis of modular and metastatic traits

2.3 GO與GSEA富集分析

GO 富集分析結(jié)果提示（圖5），在生物學(xué)進(jìn)程（BP）方面，black 模塊中基因主要參與了表皮發(fā)育、分化等進(jìn)程，在細(xì)胞成分（CC）方面，模塊基因與細(xì)胞-細(xì)胞連接、角質(zhì)化包膜、中間絲細(xì)胞骨架等組分密切相關(guān)，在分子功能（MF）方面，模塊基因主要參與了細(xì)胞內(nèi)肽酶、絲氨酸水解酶等酶活性調(diào)節(jié)及細(xì)胞骨架構(gòu)建等功能。

圖5 black模塊基因的GO富集分析Figure 5 GO enrichment analysis of black module genes

為進(jìn)一步分析轉(zhuǎn)移特征基因集的分布情況，對(duì)TM 組和TP 組進(jìn)行了GSEA 富集分析，圖6 結(jié)果顯示了在FDRq-val<0.25 條件下，標(biāo)準(zhǔn)化富集得分（NES）最高的10 個(gè)功能基因集和NES 最低的10 個(gè)功能基因集，圖7 顯示了與TM 組負(fù)相關(guān)的代表性GSEA 富集詳細(xì)信息圖，主要涉及細(xì)胞骨架中間絲的構(gòu)成、黏附以及細(xì)胞連接通道活性等生物學(xué)過(guò)程，提示細(xì)胞骨架結(jié)構(gòu)的減弱與腫瘤轉(zhuǎn)移特征的發(fā)生有密切關(guān)系。

圖6 與TM組相關(guān)的GSEA富集分析Figure 6 GSEA enrichment analysis associated with the TM group

圖7 代表性的標(biāo)準(zhǔn)化GSEA富集評(píng)分（NES）與FDR值Figure 7 Representative normalized GSEA enrichment score(NES)vs.FDR values

2.4 ROC分析

基于單個(gè)基因的RNA 表達(dá)TPM 值，設(shè)定TP 與TM為判別特征，進(jìn)行批量化ROC分析，繪制針對(duì)每個(gè)基因的ROC 曲線，計(jì)算其下面積AUC 值，篩選AUC>0.7的基因，共獲得806個(gè)基因，利用此類基因RNA 表達(dá)TPM 值對(duì)判斷TP 及TM 特征的準(zhǔn)確性在70%以上，圖8展示了代表性基因的ROC曲線。

圖8 代表性基因的ROC曲線Figure 8 ROC curves of representative genes

2.5 機(jī)器學(xué)習(xí)XGBoost分類模型的建立及性能評(píng)估

基于差異表達(dá)基因、black模塊基因和AUC>7.0基因，三者取交集后共獲得143 個(gè)基因（圖9），其中AL049555.1基因可能是污染基因，將其排除，剩余142 個(gè)基因作為機(jī)器學(xué)習(xí)的特征基因，具體基因名如下：

圖9 三部分分析結(jié)果的交集基因維恩圖Figure 9 Venn diagram of intersection genes of three-part analysis results

使用XGBoost 算法構(gòu)建機(jī)器學(xué)習(xí)模型，為驗(yàn)證包括超參數(shù)調(diào)優(yōu)步驟在內(nèi)的整個(gè)模型構(gòu)建過(guò)程，評(píng)估該算法在訓(xùn)練集上的性能進(jìn)行的嵌套交叉驗(yàn)證結(jié)果顯示MMCE 值約為0.12，說(shuō)明該算法在訓(xùn)練集上表現(xiàn)良好。

使用構(gòu)建出的分類模型對(duì)內(nèi)部測(cè)試集進(jìn)行預(yù)測(cè)，內(nèi)部測(cè)試集共有95 例樣本，預(yù)測(cè)結(jié)果繪制混淆矩陣圖[10]，如圖10 所示，內(nèi)部測(cè)試集中82 例被正確分類，13例被錯(cuò)誤分類，MMCE值約為0.14，ROC曲線分析，AUC=0.88。

圖10 內(nèi)部測(cè)試集混淆矩陣及ROC曲線Figure 10 Confusion matrix and ROC curve of internal test set

使用構(gòu)建出的分類模型對(duì)外部測(cè)試集GSE190113 進(jìn)行預(yù)測(cè)，外部測(cè)試集共有105 例樣本，如圖11 所示，外部測(cè)試集中86 例被正確分類，19 例被錯(cuò)誤分類，MMCE 值約為0.18，ROC 曲線分析，AUC=0.85。

圖11 外部測(cè)試集混淆矩陣及ROC曲線Figure 11 Confusion matrix and ROC curve of external test set

3 討論

惡性黑色素瘤占皮膚惡性腫瘤第3位，近年來(lái)，其發(fā)生率和死亡率逐年升高。由于皮膚惡性黑色素瘤易發(fā)生隱匿性轉(zhuǎn)移，往往轉(zhuǎn)移后病程進(jìn)入中晚期，預(yù)后較差。因此，及早診斷及早局部手術(shù)切除是爭(zhēng)取治愈的最佳方法。目前臨床上，如何有效判斷黑色素瘤的原發(fā)灶和轉(zhuǎn)移灶特征，尚缺乏行之有效的手段。

近年來(lái)，隨著生物信息技術(shù)的飛速發(fā)展，目前已累積了較為可觀的腫瘤相關(guān)生物大數(shù)據(jù)，利用AI技術(shù)進(jìn)行大數(shù)據(jù)分析利用，在腫瘤輔助診斷[11]、預(yù)后模型構(gòu)建[12]及腫瘤標(biāo)志物篩查[13]等方面已取得了長(zhǎng)足進(jìn)步。

本研究篩選TCGA 數(shù)據(jù)庫(kù)中，470 例皮膚惡性黑色素瘤樣本，根據(jù)其臨床特征分為TP 組和TM組，應(yīng)用R語(yǔ)言的edgeR包與DESeq2包共篩選出兩組間的差異表達(dá)基因1804 個(gè)，利用WGCNA 分析篩選出與TM 及TP 性狀相關(guān)性最高的black 模塊，該模塊包含529 個(gè)轉(zhuǎn)移特征相關(guān)基因，并利用ROC曲線分析獲得806 個(gè)AUC>0.7 的轉(zhuǎn)移特征判別基因，3 種篩選方法取交集，最后確定了142 個(gè)轉(zhuǎn)移特征基因，接著利用XGBoost 算法搭建機(jī)器學(xué)習(xí)模型，以142 個(gè)基因的RNA 表達(dá)TPM 值為輸入對(duì)象，對(duì)皮膚惡性黑色素瘤的轉(zhuǎn)移特征進(jìn)行預(yù)測(cè)，通過(guò)內(nèi)部驗(yàn)證集和外部驗(yàn)證集評(píng)估模型性能，結(jié)果證實(shí)本模型能有效判斷黑色素瘤的轉(zhuǎn)移特征，內(nèi)部測(cè)試集的預(yù)測(cè)準(zhǔn)確率為88%，外部驗(yàn)證集的預(yù)測(cè)準(zhǔn)確率為85%。

本機(jī)器學(xué)習(xí)模型尚存在一些不足之處，首先，由于不同機(jī)構(gòu)采用的檢測(cè)方法存在差異，導(dǎo)致轉(zhuǎn)錄組數(shù)據(jù)可能存在批次效應(yīng)，會(huì)對(duì)模型預(yù)測(cè)結(jié)果產(chǎn)生影響。其次，本模型采用的數(shù)據(jù)來(lái)自TCGA 數(shù)據(jù)庫(kù)，樣本為歐美白人人種，對(duì)于我國(guó)為代表的東亞人種是否適用尚需檢驗(yàn)。最后，由于存在數(shù)據(jù)量和數(shù)據(jù)特征的局限性，在實(shí)際應(yīng)用中可能會(huì)存在偏置性問(wèn)題。因此，數(shù)據(jù)的標(biāo)準(zhǔn)化、采用國(guó)內(nèi)患者樣本數(shù)據(jù)及擴(kuò)大數(shù)據(jù)量等措施將是本研究需要進(jìn)一步改進(jìn)的方向。