何東生 曾先捷

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

1 資料與方法

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

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

1.2 差異基因的篩選

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

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

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

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

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

2 結(jié)果

2.1 差異基因的篩選

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

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

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

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

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

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

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

圖2 ROC曲線判斷患者預(yù)后的準確性

3 討論

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

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

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