姜雪麗，趙 辛，趙鐵錚

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

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

1 資料與方法

1.1資料收集和差異分析 在GEO數(shù)據(jù)庫中以“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ù)庫中下載人甲狀腺癌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)值評分、生存時(shí)間和生存狀態(tài)等，去除臨床資料不完整患者，最終保留498例采用單因素和多因素Cox回歸分析進(jìn)行甲狀腺癌患者臨床病理特征與預(yù)后關(guān)系的研究。

1.2Cox比例風(fēng)險(xiǎn)回歸模型建立和驗(yàn)證 利用R軟件將TCGA數(shù)據(jù)庫中甲狀腺癌患者隨機(jī)分為兩組，即試驗(yàn)組和驗(yàn)證組，使用Survival包對在GEO和TCGA數(shù)據(jù)庫中均發(fā)生改變基因與患者預(yù)后的關(guān)系進(jìn)行單因素Cox回歸分析；使用Glmnet包對單因素Cox回歸分析結(jié)果有差異的基因(P<0.05)進(jìn)行LASSO回歸分析；使用Survival和Survminer包對經(jīng)LASSO回歸分析篩選后的基因進(jìn)行多因素Cox回歸分析，進(jìn)一步檢驗(yàn)具有風(fēng)險(xiǎn)預(yù)測能力基因并計(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)曲線評估模型預(yù)測甲狀腺癌患者生存率準(zhǔn)確性，借助R軟件及Ggplot2、Survminer、SurvivalROC、Pheatmap等功能包進(jìn)行數(shù)據(jù)可視化。

1.3GO和KEGG功能富集分析 將TCGA數(shù)據(jù)庫中甲狀腺癌基因表達(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)功能包對所有數(shù)據(jù)進(jìn)行處理分析，α=0.05為檢驗(yàn)水準(zhǔn)。

2 結(jié)果

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

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

2.2差異基因單因素Cox回歸分析 對試驗(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)基因，見圖2。

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

2.3Cox比例風(fēng)險(xiǎn)回歸模型構(gòu)建 試驗(yàn)組單因素Cox回歸分析中與甲狀腺癌患者預(yù)后有關(guān)的差異基因經(jīng)LASSO回歸分析篩選出6個(gè)候選差異基因，對其進(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)基因,見表1。風(fēng)險(xiǎn)值=(PHLDA2×-0.1028)+(GPR137B×0.0880)+(PORCN×0.1112)+(MAPK4×0.2403)+(TSPYL2×0.1465)?；陲L(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)組，見圖3a，且此模型中試驗(yàn)組1年生存率的ROC曲線下面積為0.929，見圖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)證組及整體組對Cox比例風(fēng)險(xiǎn)回歸模型進(jìn)行驗(yàn)證，生存分析發(fā)現(xiàn)，驗(yàn)證組和整體組高風(fēng)險(xiǎn)組生存時(shí)間和生存率低于低風(fēng)險(xiǎn)組，見圖4a和4b；此外，驗(yàn)證組和整體組1年生存率的ROC曲線下面積分別為0.680和0.773，見圖4c和4d。

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

2.5GO和KEGG功能富集分析 對TCGA數(shù)據(jù)庫高和低風(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)答激活信號傳導(dǎo)等顯著富集，見圖5a；KEGG功能富集分析發(fā)現(xiàn)，mTOR信號通路、甲狀腺激素信號通路和HIF1信號通路等顯著富集，見圖5b。

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

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

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

3 討論

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

本研究通過對公共數(shù)據(jù)庫中甲狀腺癌樣本的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行聯(lián)合分析，篩選出在甲狀腺癌組織和正常甲狀腺組織中差異表達(dá)的基因，對差異基因進(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)基因，隨后通過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曲線評估此模型預(yù)測甲狀腺癌患者生存率的準(zhǔn)確性，發(fā)現(xiàn)此模型具有較高準(zhǔn)確性。多因素Cox回歸分析結(jié)果顯示，甲狀腺癌患者年齡和風(fēng)險(xiǎn)值評分與預(yù)后相關(guān)。上述結(jié)果表明此模型具有一定臨床應(yīng)用價(jià)值，可用于甲狀腺癌患者預(yù)后判斷。此外，通過對高和低風(fēng)險(xiǎn)甲狀腺癌樣本進(jìn)行差異基因篩選和功能富集分析發(fā)現(xiàn)，受體介導(dǎo)的內(nèi)吞作用和先天免疫應(yīng)答激活信號傳導(dǎo)等顯著富集，mTOR信號通路、甲狀腺激素信號通路和HIF1信號通路等顯著富集。既往研究也表明，mTOR信號通路和HIF1信號通路等與甲狀腺癌細(xì)胞的增殖、遷移及侵襲等惡性行為密切相關(guān)[15]。

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

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