張志勇 唐愛華 李雙蕾 許淑華 高 美

(1 廣西中醫(yī)藥大學(xué)研究生學(xué)院，南寧市 530000，電子郵箱：1916633986@qq.com；2 廣西中醫(yī)藥大學(xué)第一附屬醫(yī)院內(nèi)分泌科，南寧市 530000)

甲狀腺癌是一種較常見的惡性腫瘤，近年來該病在全球范圍內(nèi)的發(fā)病率急劇上升[1]。甲狀腺癌主要有甲狀腺乳頭狀癌、甲狀腺濾泡癌、甲狀腺髓樣癌、甲狀腺間變性癌4種類型，其中甲狀腺乳頭狀癌最為常見，占甲狀腺癌的90%以上[2]。甲狀腺癌的預(yù)后與其組織學(xué)類型密切相關(guān)，如甲狀腺乳頭狀癌的10年總生存率約為98%，而甲狀腺間變性癌的中位生存時(shí)間只有3～5個(gè)月[3]。對于甲狀腺癌，手術(shù)切除是最重要的治療方法，而不同的輔助治療對某些病理類型的甲狀腺癌是有效的，例如放射性碘治療分化型甲狀腺癌，以及化療治療甲狀腺間變性癌。雖然大多數(shù)甲狀腺癌患者預(yù)后良好，但約10%分化良好的甲狀腺癌患者對放射性碘治療無反應(yīng)，而腫瘤分化程度差或未分化的患者更容易出現(xiàn)復(fù)發(fā)和死亡[4-5]。因此，尋找治療甲狀腺癌的新方法非常必要。糖酵解是腫瘤細(xì)胞獲取能量的重要途徑之一，阻斷糖酵解途徑有望成為治療腫瘤的潛在方法。故本研究通過生物信息學(xué)技術(shù)預(yù)測甲狀腺癌預(yù)后相關(guān)糖酵解高風(fēng)險(xiǎn)基因及其與患者預(yù)后的關(guān)系。

1 材料與方法

1.1 轉(zhuǎn)錄組數(shù)據(jù)及臨床數(shù)據(jù)的下載 從腫瘤基因組圖譜(The Cancer Genome Atlas,TCGA)數(shù)據(jù)庫(https://portal.gdc.cancer.gov)下載甲狀腺癌相關(guān)的轉(zhuǎn)錄組數(shù)據(jù)，即甲狀腺癌相關(guān)基因及基因在臨床樣本(包括癌旁組織樣本和腫瘤樣本)中的表達(dá)量；對所獲數(shù)據(jù)進(jìn)行ID轉(zhuǎn)換，即將樣本名轉(zhuǎn)換為基因名，獲取每個(gè)樣本中所包含的基因。同時(shí)下載患者的臨床數(shù)據(jù)，包括患者的生存時(shí)間、生存狀態(tài)(生存或死亡)、年齡、性別、臨床分期等基本信息，剔除樣本中信息不全的患者，得到最終的臨床樣本，獲取其基因表達(dá)數(shù)據(jù)集文件及基因表型數(shù)據(jù)文件。

1.2 基因集富集分析 在基因集富集分析(Gene Set Enrichment Analysis，GSEA)數(shù)據(jù)庫(http://www.gsea-msigdb.org/gsea/index.jsp)中搜集所有與糖酵解相關(guān)的基因集，將獲取的基因表達(dá)數(shù)據(jù)集文件、基因表型數(shù)據(jù)文件及基因集文件導(dǎo)入GSEA軟件，進(jìn)行GSEA，篩選出P<0.05的基因集[6]。

1.3 甲狀腺癌預(yù)后相關(guān)糖酵解基因的獲取 應(yīng)用Perl程序語言提取癌旁組織樣本和腫瘤樣本中糖酵解基因集(即1.2中篩選出的P<0.05的基因集)的表達(dá)量，采用Wilcoxon檢驗(yàn)驗(yàn)證癌旁組織樣本和腫瘤樣本上述基因的表達(dá)量是否具有差異，篩選出P<0.05的差異性表達(dá)基因，輸出基因表達(dá)量文件。將差異性表達(dá)基因的表達(dá)量與生存數(shù)據(jù)(生存狀態(tài)、生存時(shí)間)合并，采用R軟件對基因表達(dá)數(shù)據(jù)和生存數(shù)據(jù)進(jìn)行相關(guān)性分析，得到與甲狀腺癌預(yù)后相關(guān)的糖酵解基因。

1.4 網(wǎng)絡(luò)模型的構(gòu)建 根據(jù)甲狀腺癌預(yù)后相關(guān)的糖酵解基因，采用R軟件構(gòu)建風(fēng)險(xiǎn)函數(shù)相關(guān)的Cox回歸模型，得到糖酵解預(yù)后模型，由R軟件自動計(jì)算出每個(gè)臨床樣本的風(fēng)險(xiǎn)值及所有臨床樣本風(fēng)險(xiǎn)值的中位值，當(dāng)樣本風(fēng)險(xiǎn)值高于中位值即判斷為高風(fēng)險(xiǎn)樣本，以此將患者分為高、低風(fēng)險(xiǎn)組，采用R軟件繪制出兩組患者的生存曲線圖。并通過Cox回歸模型分析結(jié)果獲得各基因的風(fēng)險(xiǎn)比，當(dāng)風(fēng)險(xiǎn)比>1時(shí)認(rèn)為對應(yīng)的基因即為高風(fēng)險(xiǎn)基因。

1.5 受試者工作特征曲線圖、熱圖、時(shí)間-生存散點(diǎn)圖的繪制及預(yù)后分析 采用R軟件繪制出糖酵解預(yù)后模型預(yù)測甲狀腺癌患者預(yù)后的受試者工作特征 (receiver operating characteristic，ROC)曲線圖，當(dāng)曲線下面積>0.7時(shí)，表示糖酵解預(yù)后模型預(yù)測的準(zhǔn)確性高。根據(jù)糖酵解預(yù)后模型預(yù)測得到與甲狀腺癌預(yù)后相關(guān)的基因，繪制熱圖及所有樣本的時(shí)間-生存散點(diǎn)圖，并進(jìn)行預(yù)后分析。

2 結(jié) 果

2.1 臨床樣本的基本信息 共獲得504例甲狀腺癌臨床樣本，剔除3例信息不全的臨床樣本，共獲得501例臨床樣本及56 753個(gè)基因數(shù)據(jù)?；颊叩纳鏁r(shí)間為0～5 423 d，平均生存時(shí)間為1 031 d，中位生存時(shí)間為723 d；存活患者487例，死亡患者14例；發(fā)病年齡為15～89歲，平均發(fā)病年齡為45歲，中位發(fā)病年齡為46歲；女性患者366例，男性患者135例；臨床分期Ⅰ期285例，Ⅱ期52例，Ⅲ期111例，Ⅳ期53例。

2.2 GSEA結(jié)果 通過GSEA數(shù)據(jù)庫共搜索到5個(gè)與糖酵解相關(guān)的基因集，基因集名稱分別為BIOCARTA_GLYCOLYSIS_PATHWAY、HALLMARK_GLYCOLYSIS、KEGG_GLYCOLYSIS_GLUCONEOGENESIS、REACTOME_GLYCOLYSIS、REACTOME_REGULATION_OF_GLYCOLYSIS_BY_FRUCTOSE_2_6_BISPHOSPHATE_METABOLISM。對這些基因集進(jìn)行GSEA，結(jié)果顯示只有HALLMARK_GLYCOLYSIS基因集的P值和校正P值小于0.05(見表1，僅列出該基因集的GSEA結(jié)果)，且富集分析圖顯示該基因在腫瘤樣本中表達(dá)活躍(見圖1)。

表1 GSEA結(jié)果

圖1 HALLMARK_GLYCOLYSIS基因集的富集分析圖

2.3 與甲狀腺癌預(yù)后相關(guān)的糖酵解基因 獲得147個(gè)差異性表達(dá)基因，將其表達(dá)量文件與生存數(shù)據(jù)合并后進(jìn)行相關(guān)分析，共獲得與甲狀腺癌預(yù)后相關(guān)的糖酵解基因48個(gè)。

2.4 網(wǎng)絡(luò)模型構(gòu)建結(jié)果 將與甲狀腺癌預(yù)后相關(guān)的48個(gè)糖酵解基因納入Cox回歸模型，構(gòu)建糖酵解預(yù)后模型，得到與甲狀腺癌預(yù)后相關(guān)的基因共15個(gè)，其中高風(fēng)險(xiǎn)基因共8個(gè)，分別為轉(zhuǎn)化生長因子β誘導(dǎo)基因(transforming growth factor beta-induced gene,TGFBI)、溶質(zhì)載體家族16成員3(solute carrier family 16 member 3,SLC16A3)、血小板型磷酸果糖激酶(phosphofructokinase platelet, PFKP)、絲氨酸/蘇氨酸蛋白磷酸酶2A催化亞基β(serine/threonine-protein phosphatase 2A catalytic subunit beta，PPP2CB)、丙酮酸激酶M(pyruvate kinase M, PKM)、DNA損傷誘導(dǎo)轉(zhuǎn)錄子4(DNA damage inducible transcript 4,DDIT4)、碳水化合物磺基轉(zhuǎn)移酶6(carbohydrate sulfotransferase 6，CHST6)、Quiescin巰基氧化酶1(Quiescin sulfhydryl oxidase 1，QSOX1)，見表2。根據(jù)臨床樣本的風(fēng)險(xiǎn)值及其中位值，將患者劃分為高風(fēng)險(xiǎn)組、低風(fēng)險(xiǎn)組，兩組患者的生存曲線見圖2。由圖可知：在疾病的早期，隨著時(shí)間的延長患者的生存率降低，但在疾病后期，隨著時(shí)間的延長患者的生存率逐漸趨于穩(wěn)定，且高風(fēng)險(xiǎn)組的總體生存率低于低風(fēng)險(xiǎn)組(P=0.007)。

表2 甲狀腺癌預(yù)后相關(guān)的高風(fēng)險(xiǎn)基因

圖2 生存曲線圖

2.5 ROC曲線分析結(jié)果 ROC曲線分析結(jié)果顯示糖酵解預(yù)后模型預(yù)測甲狀腺癌患者預(yù)后的曲線下面積為0.783，提示根據(jù)糖酵解預(yù)后模型預(yù)測甲狀腺癌患者預(yù)后的準(zhǔn)確率較高。見圖3。

圖3 ROC曲線

2.6 熱圖及時(shí)間-生存散點(diǎn)圖 根據(jù)糖酵解預(yù)后模型得到的15個(gè)基因繪制出熱圖(見圖4)和所有樣本的時(shí)間-生存散點(diǎn)圖(見圖5)。從圖4可以看出，表達(dá)量最高的4個(gè)基因分別是PKM、PPP2CB、PFKP、DDIT4，從圖5可以看出，隨著根據(jù)Cox回歸模型計(jì)算出的風(fēng)險(xiǎn)值的增加，死亡的患者數(shù)增多。

圖4 熱圖

圖5 時(shí)間-生存散點(diǎn)圖

3 討 論

既往研究顯示，腫瘤細(xì)胞比正常細(xì)胞需要更多的葡萄糖，即使在氧氣充足的條件下，腫瘤細(xì)胞也是通過糖酵解進(jìn)行葡萄糖代謝，而不是通過線粒體氧化磷酸化產(chǎn)生腺嘌呤核苷三磷酸，這就是“瓦伯格效應(yīng)”[7]。然而，腫瘤細(xì)胞特殊的代謝模式不是被動的改變，而是遺傳表達(dá)的積極改變，其改變了腫瘤發(fā)生和侵襲性的代謝模式。有研究表明，甲狀腺癌的發(fā)生與葡萄糖攝取增加有關(guān)[8]。腫瘤細(xì)胞膜上的單羧酸轉(zhuǎn)運(yùn)蛋白(monocarboxylate transporter，MCT)與甲狀腺癌的預(yù)后有關(guān)[9]，提示腫瘤代謝的關(guān)鍵分子可能是影響甲狀腺癌預(yù)后的因素。

本研究通過生物信息學(xué)分析預(yù)測出與甲狀腺癌預(yù)后相關(guān)的糖酵解高風(fēng)險(xiǎn)基因共8個(gè)，從ROC曲線可以看出，通過Cox回歸構(gòu)建的糖酵解預(yù)后模型預(yù)測甲狀腺癌患者預(yù)后的準(zhǔn)確性較高。研究顯示，缺氧條件下源自腫瘤的外泌體包含信息性miRNA，其參與癌細(xì)胞與癌旁細(xì)胞的相互作用，從而促進(jìn)腫瘤微環(huán)境的組織重塑[10]。與正常甲狀腺濾泡細(xì)胞相比，從缺氧的甲狀腺乳頭狀癌、BCPAP細(xì)胞及KTC-1細(xì)胞中分離出的外泌體增強(qiáng)了人臍靜脈內(nèi)皮細(xì)胞的血管生成；其次，在缺氧條件下，甲狀腺乳頭狀癌BCPAP細(xì)胞的外泌體中miR-21-5p表達(dá)上調(diào)，而miR-21-5p可直接靶向并抑制TGFBI基因的表達(dá)，從而增加內(nèi)皮細(xì)胞血管的形成，從敲除miR-21-5p的低氧BCPAP細(xì)胞中分離出的外泌體促進(jìn)血管生成的作用減弱，提示缺氧的甲狀腺乳頭狀癌細(xì)胞可通過外泌體miR-21-5p/TGFBI途徑來促進(jìn)血管生成[11]。國外學(xué)者發(fā)現(xiàn)，與濾泡腺瘤相比，TGFBI在濾泡癌中過表達(dá)[12]。PKM基因廣泛分布于多種組織中，有M1和M2兩種亞型。在腫瘤細(xì)胞中PKM2被PKM1替代以逆轉(zhuǎn)warburg效應(yīng)，從而減少乳酸的產(chǎn)生和增加氧氣的消耗[13]。有學(xué)者發(fā)現(xiàn)，人甲狀腺癌BCPAP細(xì)胞和TPC1細(xì)胞中PKM2 mRNA的表達(dá)水平高于非腫瘤細(xì)胞[14]。Bikas等[15]也發(fā)現(xiàn)PKM2在甲狀腺癌FTC-133細(xì)胞和BCPAP細(xì)胞中過表達(dá)，且呈糖酵解依賴性。毒理學(xué)研究表明，DDIT4基因在甲狀腺癌中非?；钴S[16]。SLC16基因家族有14個(gè)成員，其中SLC16A1、SLC16A3、SLC16A7和SLC16A8這4個(gè)基因分別編碼MCT1、MCT4、MCT2和MCT3，參與各種組織中L-乳酸、丙酮酸、短鏈脂肪酸和單羧酸鹽藥物的轉(zhuǎn)運(yùn)[17]，使用MCT1阻斷劑能特異性靶向甲狀腺癌ATC細(xì)胞，阻斷癌細(xì)胞對丙酮酸的攝取[18-19]。但關(guān)于SLC16A3、PFKP、PPP2CB、CHPF2、QSOX1基因與甲狀腺癌的相關(guān)性，目前尚未見文獻(xiàn)報(bào)告，需進(jìn)一步研究探討。

腫瘤的糖酵解代謝能力對患者的預(yù)后有重要意義。糖酵解旺盛的腫瘤生長快、侵襲性強(qiáng)、復(fù)發(fā)迅速，患者的預(yù)后往往更差[20]。本研究的生存曲線顯示，高風(fēng)險(xiǎn)組的生存率較低風(fēng)險(xiǎn)組下降，表明攜帶高風(fēng)險(xiǎn)基因的甲狀腺癌患者的預(yù)后較差。國外有研究顯示，MCT1可作為甲狀腺癌預(yù)后的獨(dú)立預(yù)測指標(biāo)，抑制MCT1表達(dá)可使腫瘤中的惡性低氧細(xì)胞“間接饑餓”[19]。然而影響腫瘤預(yù)后的因素錯綜復(fù)雜，糖酵解代謝只是其中的一個(gè)環(huán)節(jié)，故通過糖酵解基因來判斷患者的預(yù)后并不一定完全正確，抑制糖酵解高風(fēng)險(xiǎn)基因的表達(dá)，可能有助于改善甲狀腺癌患者的預(yù)后，但仍需大量實(shí)驗(yàn)研究驗(yàn)證。

綜上所述，甲狀腺癌糖酵解途徑與TGFBI、PKM、DDIT4、SLC16等高風(fēng)險(xiǎn)基因密切相關(guān)，抑制高風(fēng)險(xiǎn)基因的表達(dá)，可能有助于改善甲狀腺癌患者的預(yù)后。但本研究僅在數(shù)據(jù)庫的基礎(chǔ)上挖掘甲狀腺癌糖酵解途徑的基因，有待今后進(jìn)一步臨床實(shí)驗(yàn)研究驗(yàn)證結(jié)果的可靠性。