基于公共數據庫構建甲狀腺癌差異基因Cox比例風險回歸模型對患者預后評估價值

姜雪麗，趙 辛，趙鐵錚

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

本研究通過對GEO和TCGA數據庫中甲狀腺癌患者轉錄組數據進行聯合分析，篩選出甲狀腺癌中差異表達的基因;隨后通過單因素Cox回歸分析，篩選出與預后相關的差異基因;最后通過多因素Cox回歸分析成功構建5個差異基因預測甲狀腺癌患者預后的模型，并結合患者臨床資料對模型進行驗證，用于判斷甲狀腺癌患者預后，現報告如下。

1 資料與方法

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

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

1.3GO和KEGG功能富集分析 將TCGA數據庫中甲狀腺癌基因表達矩陣按照風險值高低分為高風險組和低風險組兩組，使用R軟件及Limma包行差異基因篩選，符合|log2(Fold Change)|>0.585且P<0.05的基因為差異基因，利用R軟件及Enrichplot、Org.Hs.eg.db、Cluster Profiler、Ggplot2等功能包進行GO和KEGG功能富集分析及數據可視化。

1.4統計學方法 應用R 3.6.1軟件及相關功能包對所有數據進行處理分析，α=0.05為檢驗水準。

2 結果

2.1差異基因篩選 對GEO數據庫中甲狀腺癌數據集GSE138198和GSE50901進行批次校正后發現，與正常甲狀腺組織相比，甲狀腺癌組織中有897個基因表達發生改變，其中506個基因上調，391個基因下調，見圖1a；結合TCGA數據庫分析發現，在此兩個數據庫甲狀腺癌組織中有241個基因均上調，207個基因均下調，見圖1b。

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

2.2差異基因單因素Cox回歸分析 對試驗組差異基因進行單因素Cox回歸分析顯示，57個差異基因與甲狀腺癌患者預后相關(P<0.05)，其中27個差異基因為高風險基因，30個差異基因為低風險基因，見圖2。

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

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

表1 甲狀腺癌患者差異基因Cox比例風險回歸模型

圖3 甲狀腺癌患者差異基因Cox比例風險回歸模型構建分析3a和3b均為試驗組

2.4Cox比例風險回歸模型驗證 利用驗證組及整體組對Cox比例風險回歸模型進行驗證，生存分析發現，驗證組和整體組高風險組生存時間和生存率低于低風險組，見圖4a和4b；此外，驗證組和整體組1年生存率的ROC曲線下面積分別為0.680和0.773，見圖4c和4d。

圖4 甲狀腺癌患者差異基因Cox比例風險回歸模型驗證4a和4c為驗證組，4b和4d為整體組

2.5GO和KEGG功能富集分析 對TCGA數據庫高和低風險組進行差異基因分析發現，152個基因在高和低風險組甲狀腺癌樣本中差異表達。差異基因GO功能富集分析發現，受體介導的內吞作用和先天免疫應答激活信號傳導等顯著富集，見圖5a；KEGG功能富集分析發現，mTOR信號通路、甲狀腺激素信號通路和HIF1信號通路等顯著富集，見圖5b。

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

2.6甲狀腺癌患者臨床病理特征與預后關系 采用單因素和多因素Cox回歸分析評估甲狀腺癌患者臨床病理特征與預后關系。單因素Cox回歸分析結果顯示，甲狀腺癌患者年齡、臨床分期、風險值評分與預后有關(P<0.01)；多因素Cox回歸分析結果顯示，甲狀腺癌患者年齡和風險值評分與預后相關(P<0.05或P<0.01)，見表2。

表2 甲狀腺癌患者臨床病理特征與預后關系

3 討論

目前，臨床上對惡性腫瘤患者預后判斷主要基于臨床病理特征[11]，具有一定主觀性，且患者預后個體差異較大，故構建科學預后評估模型對惡性腫瘤治療效果評估具有重要意義[12]。Cox比例風險回歸模型是以生存時間和生存狀態為應變量，可同時分析多種因素對生存情況影響的一種半參數回歸模型。近年來，越來越多的學者將其應用于惡性腫瘤患者預后預測研究中[13]，通過對腫瘤樣本進行測序分析，結合患者臨床病理特征，構建基于microRNA、LncRNA、mRNA和蛋白等的預后預測模型[14]，在臨床診療中具有廣闊應用前景，吸引了廣大科研工作者的注意。

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

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

總之，本研究基于公共數據庫成功構建5個差異基因的甲狀腺癌Cox比例風險回歸模型，具有較高的準確性和可靠性，有助于臨床醫生判斷甲狀腺癌患者預后。