甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)模型的構(gòu)建

何東生 曾先捷

甲狀腺癌(thyroid cancer,TC) 作為全球發(fā)病率增速最快的腫瘤，是內(nèi)分泌系統(tǒng)中最常見(jiàn)的惡性腫瘤[1]。雖然甲狀腺癌患者5年生存率相對(duì)其他惡性腫瘤患者較高，但是在過(guò)去的十年里，甲狀腺癌患者的死亡率仍持續(xù)增加[2-3]。2017年，全球甲狀腺癌新發(fā)病例數(shù)占所有腫瘤新發(fā)病例的3.4%，在女性腫瘤中發(fā)病率上升到第五位[4-5]。甲狀腺癌已經(jīng)成為1個(gè)不容忽視的全球性問(wèn)題[6]。若能進(jìn)一步明確甲狀腺癌的獨(dú)立預(yù)后因子，對(duì)判斷甲狀腺癌患者的預(yù)后、尋找新的治療靶點(diǎn)和提高患者總體生存率將會(huì)起到至關(guān)重要的作用。目前，預(yù)測(cè)甲狀腺癌的生存預(yù)后主要通過(guò)患者的臨床-病理特征進(jìn)行判斷[7]。然而，隨著對(duì)甲狀腺癌分子機(jī)制研究的深入，通過(guò)基因預(yù)測(cè)甲狀腺患者預(yù)后正逐步受到重視[8]。本研究對(duì)TCGA數(shù)據(jù)庫(kù)中甲狀腺癌患者的mRNAs測(cè)序數(shù)據(jù)及相應(yīng)的生存數(shù)據(jù)進(jìn)行分析。在對(duì)與甲狀腺癌生存預(yù)后顯著相關(guān)的mRNAs進(jìn)行cox回歸分析之后，篩選出6個(gè)潛在的獨(dú)立預(yù)后因子并構(gòu)建預(yù)測(cè)甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)的多因素cox回歸模型，為甲狀腺癌患者的治療及改善預(yù)后提供參考依據(jù)。

1 資料與方法

1.1 數(shù)據(jù)資料收集

下載TCGA數(shù)據(jù)庫(kù)(https://portal.gdc.cancer.gov/)中人類(lèi)甲狀腺癌的RNA-Seq數(shù)據(jù)和507例患者臨床信息。RNA-Seq數(shù)據(jù)中包括58例正常甲狀腺樣本和510例甲狀腺癌樣本的測(cè)序數(shù)據(jù)。通過(guò)R軟件(https://www.r-project.org/)對(duì)所下載數(shù)據(jù)進(jìn)行整理。

1.2 差異基因的篩選

利用edgeR函數(shù)包[9]對(duì)甲狀腺正常樣本和癌樣本進(jìn)行分析，篩選出差異表達(dá)的mRNA。通過(guò)除異均值化M值法(trimmed mean of M values,TNM)對(duì)RNAs的表達(dá)量進(jìn)行標(biāo)準(zhǔn)化處理后[10]，采用錯(cuò)誤發(fā)現(xiàn)率(false discocery rate,FDR)對(duì)RNAs的顯著性進(jìn)行校正[11]。差異表達(dá)的RNAs的cut-off值設(shè)定為：①FDR≤0.01;②∣log2fold change (FC)∣≥2。

1.3 單因素和多因素Cox回歸分析

排除5例缺失mRNA測(cè)序數(shù)據(jù)的甲狀腺癌患者資料，將502例甲狀腺癌患者的生存數(shù)據(jù)與相應(yīng)的差異表達(dá)的mRNA測(cè)序數(shù)據(jù)相合并。利用survival函數(shù)包對(duì)差異表達(dá)的mRNA和甲狀腺癌患者的生存預(yù)后之間的關(guān)系進(jìn)行單因素Cox回歸分析[12]。與甲狀腺癌患者生存預(yù)后顯著相關(guān)的mRNAs(P<0.001)被納入多因素cox回歸分析以進(jìn)一步檢驗(yàn)這些基因判斷甲狀腺癌患者預(yù)后的價(jià)值。通過(guò)以上分析，構(gòu)建判斷甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)的cox回歸模型：風(fēng)險(xiǎn)值 (Risk Score)=expmRNA1×βmRNA1+expmRNA2×βmRNA2+……+expmRNAn×βmRNAn(exp:表達(dá)水平；β：多因素cox回歸分析的回歸系數(shù))[13]。

1.4 甲狀腺癌患者的高低風(fēng)險(xiǎn)劃分以及ROC 曲線(xiàn)的繪制

根據(jù)預(yù)測(cè)甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)的cox回歸模型，計(jì)算502個(gè)甲狀腺癌患者的風(fēng)險(xiǎn)值。以中位風(fēng)險(xiǎn)值為界，將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組并進(jìn)行Kaplan-Meier生存曲線(xiàn)的繪制和log-rank分析。使用survivalROC函數(shù)包繪制ROC (receiver operating characteristic)曲線(xiàn)[14]，判斷通過(guò)cox回歸模型預(yù)測(cè)甲狀腺癌患者5年生存率的準(zhǔn)確性。

2 結(jié)果

2.1 差異基因的篩選

對(duì)TCGA數(shù)據(jù)庫(kù)中58個(gè)正常甲狀腺癌樣本和510個(gè)甲狀腺癌樣本的測(cè)序數(shù)據(jù)進(jìn)行分析后，發(fā)現(xiàn)17651個(gè)基因中的1097mRNAs差異表達(dá)(233個(gè)下調(diào),864上調(diào))。

2.2 單因素cox和多因素cox回歸分析

通過(guò)survival函數(shù)包進(jìn)行單因素cox回歸分析，對(duì)甲狀腺癌中1097個(gè)差異表達(dá)的mRNAs與患者生存預(yù)后之間的關(guān)系進(jìn)行探索。如表1所示，11個(gè)差異表達(dá)的mRNA與甲狀腺癌患者的生存預(yù)后顯著相關(guān)(P<0.001)。通過(guò)多因素cox回歸分析對(duì)11個(gè)基因預(yù)測(cè)患者預(yù)后的能力進(jìn)行進(jìn)一步檢驗(yàn)后，我們發(fā)現(xiàn)RIPPLY3、PCOLCE2、FAM111B、ZSCAN4、SALL3、DLK1六個(gè)基因具有充當(dāng)獨(dú)立預(yù)后因子的功能。判斷甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)值的多因素cox回歸模型為：風(fēng)險(xiǎn)值=(RIPPLY3×-0.6194) + (PCOLCE2×0.5299) + (FAM111B×-0.3838) + (ZSCAN4×-0.4329) + (SALL3×0.3508) + (DLK1×0.2216)。RIPPLY3、FAM111B、ZSCAN4在甲狀腺癌患者低風(fēng)險(xiǎn)組的表達(dá)量高于高風(fēng)險(xiǎn)組，PCOLCE2、SALL3、DLK1在甲狀腺癌低風(fēng)險(xiǎn)組的表達(dá)量低于高風(fēng)險(xiǎn)組。

表1 與甲狀腺癌患者生存預(yù)后顯著相關(guān)的差異表達(dá)基因

2.3 甲狀腺癌患者的高低風(fēng)險(xiǎn)劃分以及ROC曲線(xiàn)的繪制

根據(jù)判斷甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)值的多因素cox回歸模型，分別計(jì)算502個(gè)甲狀腺癌患者的預(yù)后風(fēng)險(xiǎn)值，以中位數(shù)為界，將患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。Kaplan-Meier生存分析結(jié)果顯示甲狀腺癌高風(fēng)險(xiǎn)患者與低風(fēng)險(xiǎn)患者的生存率之間存在顯著差異，高風(fēng)險(xiǎn)組甲狀腺癌患者的生存預(yù)后顯著劣于低風(fēng)險(xiǎn)組甲狀腺癌患者，見(jiàn)圖1。高風(fēng)險(xiǎn)組甲狀腺癌患者的五年生存率為84.7%，低風(fēng)險(xiǎn)組甲狀腺癌患者的生存率為100%。ROC曲線(xiàn)下面積為0.955(圖2)，確認(rèn)了基于多因素cox回歸模型預(yù)測(cè)甲狀腺癌患者預(yù)后風(fēng)險(xiǎn)的準(zhǔn)確性。

圖1 甲狀腺癌患者高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組的Kaplan-Meier生存分析

圖2 ROC曲線(xiàn)判斷患者預(yù)后的準(zhǔn)確性

3 討論

甲狀腺癌作為內(nèi)分泌系統(tǒng)最常見(jiàn)的惡性腫瘤，在形態(tài)學(xué)特征和預(yù)后方面存在著巨大的異質(zhì)性[6]。由于甲狀腺癌復(fù)發(fā)率的增高，導(dǎo)致甲狀腺癌患者不可治愈率以及死亡率有所增加[15]。目前，預(yù)測(cè)腫瘤患者生存的預(yù)后主要依靠患者臨床特征，如年齡，性別，腫瘤位置、大小、分期、轉(zhuǎn)移等[7-8]。然而，隨著對(duì)甲狀腺癌分子層面研究的深入，這一觀(guān)點(diǎn)正不斷受到質(zhì)疑[8]。本研究通過(guò)對(duì)甲狀腺癌中異常表達(dá)的基因進(jìn)行探索，篩選出與患者總體生存預(yù)后顯著相關(guān)的核心基因，從基因表達(dá)方面對(duì)甲狀腺癌患者的預(yù)后風(fēng)險(xiǎn)進(jìn)行預(yù)測(cè)。

近幾年，微陣列和高通量測(cè)序技術(shù)的發(fā)展在揭示癌癥發(fā)生發(fā)展過(guò)程中基因的異常表達(dá)和識(shí)別與癌癥診斷和預(yù)后相關(guān)的生物標(biāo)記物方面起到至關(guān)重要的作用[16]。本研究對(duì)TCGA數(shù)據(jù)庫(kù)中甲狀腺癌患者的基因表達(dá)譜及相應(yīng)的臨床信息進(jìn)行綜合性分析。通過(guò)對(duì)58例正常甲狀腺樣本和510例甲狀腺癌樣本的測(cè)序數(shù)據(jù)進(jìn)行分析，我們篩選出了233個(gè)下調(diào)和864上調(diào)的mRNAs。在對(duì)1097個(gè)差異表達(dá)的mRNAs進(jìn)行單因素cox回歸分析篩選出與甲狀腺癌患者預(yù)后顯著相關(guān)的mRNAs之后，將其納入多因素cox回歸分析并篩選出6個(gè)甲狀腺癌的獨(dú)立預(yù)后因子：RIPPLY3、PCOLCE2、FAM111B、ZSCAN4、SALL3和DLK1。基于這6個(gè)基因構(gòu)建具有預(yù)測(cè)患者預(yù)后風(fēng)險(xiǎn)功能的多因素cox回歸模型。RIPPLY3、FAM111B、ZSCAN4基因的多因素cox回歸系數(shù)均小于0，提示這三個(gè)基因高表達(dá)患者的預(yù)后會(huì)優(yōu)于低表達(dá)的患者。PCOLCE2、SALL3、DLK1基因的多因素cox回歸系數(shù)均大于0，提示這3個(gè)基因低表達(dá)患者的預(yù)后會(huì)優(yōu)于高表達(dá)的患者。在對(duì)502例甲狀腺癌患者進(jìn)行風(fēng)險(xiǎn)值的評(píng)分后，將患者劃分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組并進(jìn)行Kaplan-Meier生存分析，結(jié)果證明高風(fēng)險(xiǎn)組甲狀腺癌患者的總體生存率顯著劣于低風(fēng)險(xiǎn)組患者。ROC曲線(xiàn)驗(yàn)證了通過(guò)6個(gè)甲狀腺癌獨(dú)立預(yù)后因子判斷甲狀腺癌患者風(fēng)險(xiǎn)值的可靠性。通過(guò)預(yù)測(cè)甲狀腺癌患者的預(yù)后風(fēng)險(xiǎn)值，可以為甲狀腺癌患者進(jìn)行個(gè)性化治療提供參考，從而降低患者死亡率，改善患者預(yù)后。本研究的不足之處在于所有數(shù)據(jù)均來(lái)自TCGA數(shù)據(jù)庫(kù)提供的甲狀腺癌mRNA測(cè)序數(shù)據(jù)，可能無(wú)法完全代表其蛋白水平表達(dá)的情況。在后續(xù)的實(shí)驗(yàn)中，應(yīng)該結(jié)合免疫組化、蛋白印記以及細(xì)胞功能等試驗(yàn)進(jìn)一步討論和驗(yàn)證篩選出來(lái)的6個(gè)甲狀腺癌獨(dú)立預(yù)后因子在甲狀腺癌中的表達(dá)情況以及在預(yù)測(cè)患者生存預(yù)后方面所起的作用。雖然6個(gè)基因在甲狀腺癌的蛋白表達(dá)水平以及功能有待進(jìn)一步研究，但在其他癌種中已有相關(guān)報(bào)道。ZSCAN4通過(guò)調(diào)節(jié)端粒的長(zhǎng)度維持基因組穩(wěn)定，在腫瘤細(xì)胞以及胚胎干細(xì)胞中，ZSCAN4在延伸端粒長(zhǎng)度方面起到重要調(diào)節(jié)作用[17-18]；PCOLCE2在卵巢癌中起到了促癌因子的作用，不利于患者臨床預(yù)后；FAM111B的低表達(dá)被認(rèn)為與子宮頸癌患者的遠(yuǎn)處轉(zhuǎn)移有關(guān)[19-20]。此外，DLK1被認(rèn)為可作為胚胎癌的治療靶點(diǎn)[21]。

總之，本研究通過(guò)對(duì)甲狀腺癌患者的基因表達(dá)譜及其臨床信息進(jìn)行綜合分析，發(fā)現(xiàn)了6個(gè)最有可能具有預(yù)測(cè)患者生存預(yù)后功能的獨(dú)立因素并構(gòu)建預(yù)后風(fēng)險(xiǎn)預(yù)測(cè)模型，對(duì)甲狀腺癌患者的診斷和治療和改善預(yù)后提供了參考依據(jù)。