基于公共數(shù)據(jù)庫(kù)構(gòu)建甲狀腺癌差異基因Cox比例風(fēng)險(xiǎn)回歸模型對(duì)患者預(yù)后評(píng)估價(jià)值

姜雪麗，趙 辛，趙鐵錚

甲狀腺癌是一種好發(fā)于中老年女性的內(nèi)分泌系統(tǒng)常見(jiàn)惡性腫瘤，多數(shù)病理類型為甲狀腺乳頭狀癌[1-4]。據(jù)統(tǒng)計(jì)，2017年全球新發(fā)甲狀腺癌病例占所有惡性腫瘤新發(fā)病例3.4%，在女性惡性腫瘤新發(fā)病例中甲狀腺癌居第5位[5]。我國(guó)每年甲狀腺癌新發(fā)病例和死亡病例數(shù)占全世界15%左右[6]。近年甲狀腺癌患者病死率逐年增加[7-9]。目前，臨床上對(duì)甲狀腺癌患者預(yù)后判斷主要基于臨床病理特征，但由于惡性腫瘤有明顯異質(zhì)性，且患者個(gè)體差異較大，故此預(yù)測(cè)手段具有明顯局限性[10]。近年，腫瘤基因和蛋白表達(dá)預(yù)測(cè)患者預(yù)后正逐步受到重視。

本研究通過(guò)對(duì)GEO和TCGA數(shù)據(jù)庫(kù)中甲狀腺癌患者轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行聯(lián)合分析，篩選出甲狀腺癌中差異表達(dá)的基因;隨后通過(guò)單因素Cox回歸分析，篩選出與預(yù)后相關(guān)的差異基因;最后通過(guò)多因素Cox回歸分析成功構(gòu)建5個(gè)差異基因預(yù)測(cè)甲狀腺癌患者預(yù)后的模型，并結(jié)合患者臨床資料對(duì)模型進(jìn)行驗(yàn)證，用于判斷甲狀腺癌患者預(yù)后，現(xiàn)報(bào)告如下。

1 資料與方法

1.1資料收集和差異分析 在GEO數(shù)據(jù)庫(kù)中以“Thyroid cancer”為關(guān)鍵詞，選取人甲狀腺癌組織樣本轉(zhuǎn)錄組數(shù)據(jù)集GSE138198和GSE50901，使用R軟件及Sva功能包校正批次，利用Limma包分析差異基因，認(rèn)為符合|log2(Fold Change)|>0.585且P<0.05的基因?yàn)椴町惢?。在TCGA數(shù)據(jù)庫(kù)中下載人甲狀腺癌RNA-Seq數(shù)據(jù)，共含有58份正常甲狀腺樣本和509份甲狀腺癌樣本，利用R軟件及Limma包分析差異基因，認(rèn)為符合|log2(Fold Change)|>2且P<0.01的基因?yàn)椴町惢?。甲狀腺癌患者的臨床資料利用R軟件進(jìn)行整理，包括性別、年齡、臨床分期、風(fēng)險(xiǎn)值評(píng)分、生存時(shí)間和生存狀態(tài)等，去除臨床資料不完整患者，最終保留498例采用單因素和多因素Cox回歸分析進(jìn)行甲狀腺癌患者臨床病理特征與預(yù)后關(guān)系的研究。

1.2Cox比例風(fēng)險(xiǎn)回歸模型建立和驗(yàn)證 利用R軟件將TCGA數(shù)據(jù)庫(kù)中甲狀腺癌患者隨機(jī)分為兩組，即試驗(yàn)組和驗(yàn)證組，使用Survival包對(duì)在GEO和TCGA數(shù)據(jù)庫(kù)中均發(fā)生改變基因與患者預(yù)后的關(guān)系進(jìn)行單因素Cox回歸分析；使用Glmnet包對(duì)單因素Cox回歸分析結(jié)果有差異的基因(P<0.05)進(jìn)行LASSO回歸分析；使用Survival和Survminer包對(duì)經(jīng)LASSO回歸分析篩選后的基因進(jìn)行多因素Cox回歸分析，進(jìn)一步檢驗(yàn)具有風(fēng)險(xiǎn)預(yù)測(cè)能力基因并計(jì)算各基因回歸系數(shù)，從而構(gòu)建判斷患者預(yù)后的Cox比例風(fēng)險(xiǎn)回歸模型。風(fēng)險(xiǎn)值為樣本中各基因表達(dá)量與回歸系數(shù)乘積之和，計(jì)算試驗(yàn)組、驗(yàn)證組及整體組(TCGA數(shù)據(jù)整體)各樣本風(fēng)險(xiǎn)值，選取試驗(yàn)組風(fēng)險(xiǎn)值中位值為臨界值，將各組分別分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，繪制Kaplan-meier生存曲線，采用受試者工作特征(ROC)曲線評(píng)估模型預(yù)測(cè)甲狀腺癌患者生存率準(zhǔn)確性，借助R軟件及Ggplot2、Survminer、SurvivalROC、Pheatmap等功能包進(jìn)行數(shù)據(jù)可視化。

1.3GO和KEGG功能富集分析 將TCGA數(shù)據(jù)庫(kù)中甲狀腺癌基因表達(dá)矩陣按照風(fēng)險(xiǎn)值高低分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組兩組，使用R軟件及Limma包行差異基因篩選，符合|log2(Fold Change)|>0.585且P<0.05的基因?yàn)椴町惢颍肦軟件及Enrichplot、Org.Hs.eg.db、Cluster Profiler、Ggplot2等功能包進(jìn)行GO和KEGG功能富集分析及數(shù)據(jù)可視化。

1.4統(tǒng)計(jì)學(xué)方法 應(yīng)用R 3.6.1軟件及相關(guān)功能包對(duì)所有數(shù)據(jù)進(jìn)行處理分析，α=0.05為檢驗(yàn)水準(zhǔn)。

2 結(jié)果

2.1差異基因篩選 對(duì)GEO數(shù)據(jù)庫(kù)中甲狀腺癌數(shù)據(jù)集GSE138198和GSE50901進(jìn)行批次校正后發(fā)現(xiàn)，與正常甲狀腺組織相比，甲狀腺癌組織中有897個(gè)基因表達(dá)發(fā)生改變，其中506個(gè)基因上調(diào)，391個(gè)基因下調(diào)，見(jiàn)圖1a；結(jié)合TCGA數(shù)據(jù)庫(kù)分析發(fā)現(xiàn)，在此兩個(gè)數(shù)據(jù)庫(kù)甲狀腺癌組織中有241個(gè)基因均上調(diào)，207個(gè)基因均下調(diào)，見(jiàn)圖1b。

圖1 甲狀腺癌組織和正常甲狀腺組織差異基因篩選

2.2差異基因單因素Cox回歸分析 對(duì)試驗(yàn)組差異基因進(jìn)行單因素Cox回歸分析顯示，57個(gè)差異基因與甲狀腺癌患者預(yù)后相關(guān)(P<0.05)，其中27個(gè)差異基因?yàn)楦唢L(fēng)險(xiǎn)基因，30個(gè)差異基因?yàn)榈惋L(fēng)險(xiǎn)基因，見(jiàn)圖2。

圖2 甲狀腺癌患者差異基因單因素Cox回歸分析

2.3Cox比例風(fēng)險(xiǎn)回歸模型構(gòu)建 試驗(yàn)組單因素Cox回歸分析中與甲狀腺癌患者預(yù)后有關(guān)的差異基因經(jīng)LASSO回歸分析篩選出6個(gè)候選差異基因，對(duì)其進(jìn)行多因素Cox比例風(fēng)險(xiǎn)回歸模型構(gòu)建，最后發(fā)現(xiàn)PHLDA2、GPR137B、PORCN、MAPK4和TSPYL2共5個(gè)基因參與模型構(gòu)建，其中PHLDA2為低風(fēng)險(xiǎn)基因，GPR137B、PORCN、MAPK4和TSPYL2為高風(fēng)險(xiǎn)基因,見(jiàn)表1。風(fēng)險(xiǎn)值=(PHLDA2×-0.1028)+(GPR137B×0.0880)+(PORCN×0.1112)+(MAPK4×0.2403)+(TSPYL2×0.1465)。基于風(fēng)險(xiǎn)值將試驗(yàn)組分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組兩組，生存分析發(fā)現(xiàn)，高風(fēng)險(xiǎn)組生存時(shí)間和生存率低于低風(fēng)險(xiǎn)組，見(jiàn)圖3a，且此模型中試驗(yàn)組1年生存率的ROC曲線下面積為0.929，見(jiàn)圖3b。

表1 甲狀腺癌患者差異基因Cox比例風(fēng)險(xiǎn)回歸模型

圖3 甲狀腺癌患者差異基因Cox比例風(fēng)險(xiǎn)回歸模型構(gòu)建分析3a和3b均為試驗(yàn)組

2.4Cox比例風(fēng)險(xiǎn)回歸模型驗(yàn)證 利用驗(yàn)證組及整體組對(duì)Cox比例風(fēng)險(xiǎn)回歸模型進(jìn)行驗(yàn)證，生存分析發(fā)現(xiàn)，驗(yàn)證組和整體組高風(fēng)險(xiǎn)組生存時(shí)間和生存率低于低風(fēng)險(xiǎn)組，見(jiàn)圖4a和4b；此外，驗(yàn)證組和整體組1年生存率的ROC曲線下面積分別為0.680和0.773，見(jiàn)圖4c和4d。

圖4 甲狀腺癌患者差異基因Cox比例風(fēng)險(xiǎn)回歸模型驗(yàn)證4a和4c為驗(yàn)證組，4b和4d為整體組

2.5GO和KEGG功能富集分析 對(duì)TCGA數(shù)據(jù)庫(kù)高和低風(fēng)險(xiǎn)組進(jìn)行差異基因分析發(fā)現(xiàn)，152個(gè)基因在高和低風(fēng)險(xiǎn)組甲狀腺癌樣本中差異表達(dá)。差異基因GO功能富集分析發(fā)現(xiàn)，受體介導(dǎo)的內(nèi)吞作用和先天免疫應(yīng)答激活信號(hào)傳導(dǎo)等顯著富集，見(jiàn)圖5a；KEGG功能富集分析發(fā)現(xiàn)，mTOR信號(hào)通路、甲狀腺激素信號(hào)通路和HIF1信號(hào)通路等顯著富集，見(jiàn)圖5b。

圖5 甲狀腺癌樣本差異基因GO和KEGG功能富集分析5a為GO功能富集分析，5b為KEGG功能富集分析

2.6甲狀腺癌患者臨床病理特征與預(yù)后關(guān)系 采用單因素和多因素Cox回歸分析評(píng)估甲狀腺癌患者臨床病理特征與預(yù)后關(guān)系。單因素Cox回歸分析結(jié)果顯示，甲狀腺癌患者年齡、臨床分期、風(fēng)險(xiǎn)值評(píng)分與預(yù)后有關(guān)(P<0.01)；多因素Cox回歸分析結(jié)果顯示，甲狀腺癌患者年齡和風(fēng)險(xiǎn)值評(píng)分與預(yù)后相關(guān)(P<0.05或P<0.01)，見(jiàn)表2。

表2 甲狀腺癌患者臨床病理特征與預(yù)后關(guān)系

3 討論

目前，臨床上對(duì)惡性腫瘤患者預(yù)后判斷主要基于臨床病理特征[11]，具有一定主觀性，且患者預(yù)后個(gè)體差異較大，故構(gòu)建科學(xué)預(yù)后評(píng)估模型對(duì)惡性腫瘤治療效果評(píng)估具有重要意義[12]。Cox比例風(fēng)險(xiǎn)回歸模型是以生存時(shí)間和生存狀態(tài)為應(yīng)變量，可同時(shí)分析多種因素對(duì)生存情況影響的一種半?yún)?shù)回歸模型。近年來(lái)，越來(lái)越多的學(xué)者將其應(yīng)用于惡性腫瘤患者預(yù)后預(yù)測(cè)研究中[13]，通過(guò)對(duì)腫瘤樣本進(jìn)行測(cè)序分析，結(jié)合患者臨床病理特征，構(gòu)建基于microRNA、LncRNA、mRNA和蛋白等的預(yù)后預(yù)測(cè)模型[14]，在臨床診療中具有廣闊應(yīng)用前景，吸引了廣大科研工作者的注意。

本研究通過(guò)對(duì)公共數(shù)據(jù)庫(kù)中甲狀腺癌樣本的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行聯(lián)合分析，篩選出在甲狀腺癌組織和正常甲狀腺組織中差異表達(dá)的基因，對(duì)差異基因進(jìn)行單因素Cox回歸分析顯示，有57個(gè)差異基因與甲狀腺癌患者預(yù)后相關(guān)，其中27個(gè)差異基因?yàn)楦唢L(fēng)險(xiǎn)基因，30個(gè)差異基因?yàn)榈惋L(fēng)險(xiǎn)基因，隨后通過(guò)LASSO回歸分析和多因素Cox回歸分析構(gòu)建由PHLDA2、GPR137B、PORCN、MAPK4和TSPYL2共5個(gè)基因組成的Cox比例風(fēng)險(xiǎn)回歸模型。基于此模型，將試驗(yàn)組、驗(yàn)證組和整體組分別分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，生存分析發(fā)現(xiàn)各高風(fēng)險(xiǎn)組生存時(shí)間和生存率低于各低風(fēng)險(xiǎn)組，采用ROC曲線評(píng)估此模型預(yù)測(cè)甲狀腺癌患者生存率的準(zhǔn)確性，發(fā)現(xiàn)此模型具有較高準(zhǔn)確性。多因素Cox回歸分析結(jié)果顯示，甲狀腺癌患者年齡和風(fēng)險(xiǎn)值評(píng)分與預(yù)后相關(guān)。上述結(jié)果表明此模型具有一定臨床應(yīng)用價(jià)值，可用于甲狀腺癌患者預(yù)后判斷。此外，通過(guò)對(duì)高和低風(fēng)險(xiǎn)甲狀腺癌樣本進(jìn)行差異基因篩選和功能富集分析發(fā)現(xiàn)，受體介導(dǎo)的內(nèi)吞作用和先天免疫應(yīng)答激活信號(hào)傳導(dǎo)等顯著富集，mTOR信號(hào)通路、甲狀腺激素信號(hào)通路和HIF1信號(hào)通路等顯著富集。既往研究也表明，mTOR信號(hào)通路和HIF1信號(hào)通路等與甲狀腺癌細(xì)胞的增殖、遷移及侵襲等惡性行為密切相關(guān)[15]。

對(duì)本研究構(gòu)建模型的5個(gè)差異基因進(jìn)行文獻(xiàn)檢索發(fā)現(xiàn)，盡管其在甲狀腺癌中的作用報(bào)道較少，但在其他類型腫瘤中有較多研究。GPR137B在本研究模型中是一個(gè)高風(fēng)險(xiǎn)基因，在甲狀腺癌組織中高表達(dá)(數(shù)據(jù)未展示)，其在腫瘤中的作用以往未見(jiàn)報(bào)道，但敲低其同一家族基因GPR137已被證實(shí)可以抑制卵巢癌、胰腺癌和肝癌等細(xì)胞的增殖[16-17]；PORCN蛋白介導(dǎo)WNT蛋白棕櫚酰化，對(duì)于WNT蛋白的分泌及WNT信號(hào)通路的激活具有重要意義，而WNT信號(hào)通路激活與甲狀腺癌惡性進(jìn)展密切相關(guān)[18]；MAPK4過(guò)表達(dá)與肺腺癌、膀胱癌、低級(jí)別膠質(zhì)瘤和甲狀腺癌生存差異相關(guān)，過(guò)表達(dá)MAPK4通過(guò)激活A(yù)KT/mTOR信號(hào)通路促進(jìn)腫瘤惡性進(jìn)展[19]。本研究結(jié)果與上述研究結(jié)果基本相符，提示GPR137B、PORCN和MAPK4在甲狀腺癌中具有重要生理功能。然而，由于甲狀腺癌研究隊(duì)列的不足，本研究將TCGA數(shù)據(jù)庫(kù)中甲狀腺癌患者通過(guò)隨機(jī)分組的方式構(gòu)建模型并驗(yàn)證，故需要更多其他甲狀腺癌隊(duì)列研究來(lái)進(jìn)一步驗(yàn)證本研究所構(gòu)建的模型，且GPR137B、PORCN和MAPK4等基因?qū)谞钕侔┘?xì)胞增殖、遷移和侵襲等的影響需要后續(xù)細(xì)胞和動(dòng)物實(shí)驗(yàn)進(jìn)一步探討。

總之，本研究基于公共數(shù)據(jù)庫(kù)成功構(gòu)建5個(gè)差異基因的甲狀腺癌Cox比例風(fēng)險(xiǎn)回歸模型，具有較高的準(zhǔn)確性和可靠性，有助于臨床醫(yī)生判斷甲狀腺癌患者預(yù)后。