基于生物信息學預測甲狀腺癌預后相關高風險糖酵解基因及其與患者預后的關系▲

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

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

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

1 材料與方法

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

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

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

1.4 網絡模型的構建 根據甲狀腺癌預后相關的糖酵解基因，采用R軟件構建風險函數相關的Cox回歸模型，得到糖酵解預后模型，由R軟件自動計算出每個臨床樣本的風險值及所有臨床樣本風險值的中位值，當樣本風險值高于中位值即判斷為高風險樣本，以此將患者分為高、低風險組，采用R軟件繪制出兩組患者的生存曲線圖。并通過Cox回歸模型分析結果獲得各基因的風險比，當風險比>1時認為對應的基因即為高風險基因。

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

2 結 果

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

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

表1 GSEA結果

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

2.3 與甲狀腺癌預后相關的糖酵解基因 獲得147個差異性表達基因，將其表達量文件與生存數據合并后進行相關分析，共獲得與甲狀腺癌預后相關的糖酵解基因48個。

2.4 網絡模型構建結果 將與甲狀腺癌預后相關的48個糖酵解基因納入Cox回歸模型，構建糖酵解預后模型，得到與甲狀腺癌預后相關的基因共15個，其中高風險基因共8個，分別為轉化生長因子β誘導基因(transforming growth factor beta-induced gene,TGFBI)、溶質載體家族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損傷誘導轉錄子4(DNA damage inducible transcript 4,DDIT4)、碳水化合物磺基轉移酶6(carbohydrate sulfotransferase 6，CHST6)、Quiescin巰基氧化酶1(Quiescin sulfhydryl oxidase 1，QSOX1)，見表2。根據臨床樣本的風險值及其中位值，將患者劃分為高風險組、低風險組，兩組患者的生存曲線見圖2。由圖可知：在疾病的早期，隨著時間的延長患者的生存率降低，但在疾病后期，隨著時間的延長患者的生存率逐漸趨于穩定，且高風險組的總體生存率低于低風險組(P=0.007)。

表2 甲狀腺癌預后相關的高風險基因

圖2 生存曲線圖

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

圖3 ROC曲線

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

圖4 熱圖

圖5 時間-生存散點圖

3 討 論

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

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

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

綜上所述，甲狀腺癌糖酵解途徑與TGFBI、PKM、DDIT4、SLC16等高風險基因密切相關，抑制高風險基因的表達，可能有助于改善甲狀腺癌患者的預后。但本研究僅在數據庫的基礎上挖掘甲狀腺癌糖酵解途徑的基因，有待今后進一步臨床實驗研究驗證結果的可靠性。