楊仁義, 彭書旺, 王永恒, 董宇軒, 段姍姍
(1.湖南中醫藥大學第一附屬醫院胃腸甲狀腺血管外科,湖南 長沙 410007;2.湖南中醫藥大學研究生院,湖南 長沙 410208)
國家癌癥中心發布的《2022 年全國癌癥報告》[1]顯示:甲狀腺癌(thyroid cancer,TC)新發病例為20.3 萬人,其中女性新發病例(15.3 萬人)明顯高于男性(5.0 萬人),是發病率最高的內分泌惡性腫瘤,也是發病率上升最快且位列第7 位的常見惡性腫瘤。95%的TC 為分化型TC,外科手術是目前首選的治療方式,同時放射性碘治療(radioactive iodine,RAI)和內分泌輔助治療也是提高患者生存預后和生活質量的有效手段[2-3],但治療后遠期復發率仍較高,針對免疫檢查點或分子靶向治療藥物開發大多數僅停留在療效觀察階段[4-5],尚未取得分子機制層面的突破。因此,研究TC 潛在的分子機制,尋找預后生存的生物標志物和治療靶點,構建有價值的預后評估模型是亟待解決的問題。鐵死亡是一種主要依賴于鐵介導的氧化損傷和隨后的細胞膜損傷的調節性細胞死亡方式[6],可通過外源性或轉運蛋白依賴,以及內源性或酶調節2 條主要途徑啟動,同時受鐵積累增加、自由基產生、脂肪酸供應和脂質過氧化共同調控[7],在一定程度上抑制TC細胞的增殖、侵襲及轉移。研 究[8-10]顯 示:鐵 死 亡 相 關 基 因(ferroptosisrelated genes,FRGs)是TC 預后分子及潛在治療靶點,能精準預測TC 患者生存預后,有效指導TC 患者的個體化治療。WANG 等[11]鑒定出醛酮還原酶家族1 成員 C3 (Aldo-Keto reductase family 1 member C3,AKR1C3)、BH3 相互作用域死亡激動劑(BH3 interacting domain death agonist,BID)、F-Box和 WD-40結構域蛋白7(F-Box and WD repeat domain containing 7,FBXW7)、谷胱甘肽過氧化物酶4 (glutathione peroxidase 4,GPX4)和絲裂原活化蛋白激酶激酶5(mitogen-activated protein kinase kinase kinase 5,MAP3K5)具有預后價值,同時敲低AKR1C3 可增強TC 細胞的增殖、侵襲和遷移能力,為TC 預后預測提供了新思路,但未對其相關分子機制進行深入探討。本研究采用生物信息學方法,提取GeneCards 和FerrDb 數據庫中FRGs,篩選TC 差異預后鐵死亡相關基因(prognosis and ferroptosis related genes,PFRGs),構建風險評分及預后模型,評估TC患者生存預后,進行共表達及富集分析,研究TC的FRGs調控網絡及作用機制。
1.1 TC 差 異 預 后 基 因(TC-differentiallyprognostic gene,TC-DPGs)獲取數據來源:從TCGA 數據庫(https://portal.gdc.cancer.gov/)下載TCGA-THCA 隊列中TC 轉錄組及臨床數據。差異分析:采用R 軟件DESeq2 包[12],對TCGATHCA 轉錄組數據進行差異分析,以|log2FC|>1且校正P<0.05 為篩選條件,篩選TC 差異表達基因 (thyroid cancer differentially expressed genes,TC-DEGs)。預 后 分 析:采 用R 軟件survminer 包和survival 包[13],對TCGA-THCA 臨 床 數 據 進 行預后分析,以TC 患者總生存(overall survival,OS)時間為預后參數,以Cox 回歸P<0.05 為篩選條件進行單因素Cox 回歸分析,篩選TC 預后基因(TC prognosis gene,TC-PGs)。通過韋恩圖在線取交集工具(http://bioinformatics.psb.ugent.be/webtools/Venn/)取“TC-DEGs”與“TCPGs”兩者的交集,獲取TC-DPGs。
1.2 TC 差異PFRGs 篩選數 據 來 源:通 過GeneCards 數 據 庫 (https://www.genecards.org/)[14],以“Ferroptosis”為 檢 索 詞,同 時 從FerrDb 數 據 庫 (http://www.zhounan.org/ferrdb)[15]下載FRGs,整合2 個數據庫基因,最終獲取FRGs。通過韋恩圖在線取交集工具,取“TC-DPGs”與“FRGs”交 集,獲 取TC 差 異PFRGs。
1.3 TC 差異PFRGs 表達分析數據來源:從TCGA 數據庫下載TCGA-THCA 隊列中TC 轉錄組數據,從GTEx 數據庫(https://gtexportal.org/home/)下載正常甲狀腺組織表達數據。將TCGA 數據庫中癌旁組織樣本與GTEx 數據庫中正常甲狀腺組織樣本作為對照組,將TCGA 數據庫中TC 組織樣本作為腫瘤組,進行非配對樣本秩和檢驗;對包含腫瘤組織與癌旁組織的患者樣本進行配對樣本秩和檢驗,比較腫瘤組與對照組TC 差異PFRGs [CD44、膜 聯 蛋 白A1 (Annexin A1,ANXA1)和 核 受 體 亞 家 族4A 類 成 員1 (nuclear receptor sub family 4 group A mumber 1,NR4A1)]表達的差異。為進一步分析TC 差異PFRGs(CD44、ANXA1 和NR4A1)蛋白在組織中表達情況,通過人類蛋白圖譜(The Human Protein Atlas,HPA)數 據 庫[16](https://www.proteinatlas.org/),檢索CD44、ANXA1 和NR4A1蛋白的免疫組織化學切片圖像,比較腫瘤組織與正常組織蛋白表達差異。
1.4 TC 差異PFRGs 生存預后分析采用R 軟件timeROC 包,繪制時間依賴性受試者工作特征(time-receiver operating characteristic,time-ROC)曲線,分別評估TC 差異PFRGs CD44、ANXA1和NR4A1 對TC 患者的1、3 和5 年生存預測效能。當ROC 曲線下面積(area under curve,AUC)>0.5 時,AUC 越 接 近1,對TC 患 者 的1、3 和5 年生存預測效能越好,分子的表達促進死亡發生;當AUC<0.5時,AUC越接近0,對TC患者的1、3 和5年生存預測效能越好,分子的表達促進生存發生。采用R 軟件survival 包和survminer 包,繪制Kaplan-Meier 曲線,根據基因表達的中位數分為高表達和低表達2 組,分別評估TC 差異PFRGs CD44、ANXA1 和NR4A1 對TC 患者生存狀況的影響。

1.6 TC 差異PFRGs Nomogram 圖構建及驗證采用R 軟件rms 包,使用多因素Cox 回歸分析方法建立Nomogram 預后模型,以評估和預測TC患者OS 的概率[17]。綜合分析“1.5”中單因素Cox 回歸分析結果,納入多因素Cox 回歸分析中P<0.05 的臨床特征作為因子,構建預測TC 患者1、3 和5 年OS 概率的Nomogram 預后模型。采用R 軟件rms 包和survival 包,對Nomogram 預后模型進行Calibration 分析驗證,以評估Nomogram 預后模型的預測能力。
1.7 TC 差異PFRGs 與蛋白共表達相關性分析
采 用R 軟 件limma 包[18],對TC 差 異PFRGs CD44、ANXA1 和NR4A1 進 行Spearman 相 關 性 分析,以相關系數R>0.7 且P<0.05 為篩選條件,篩選TC 差異PFRGs 的共表達基因,繪制共表達熱圖。整合CD44、ANXA1 和NR4A1 的共表達基因,通過STRING(https://cn.string-db.org/)[19]數據庫進行蛋白互作分析,采用Cytoscape 3.9.0 繪制蛋白-蛋白互作(protein-protein interaction,PPI)網絡圖,利用cytoHubba 插件[20]進行度中心性(degree)、接近中心性(closeness)和中介中心性(betweenness)拓撲分析,探尋TC 鐵死亡核心網絡。
1.8 TC 差 異 PFRGs 基 因 本 體 論(Gene Ontology,GO)和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes,KEGG)富集分析采 用R 軟 件clusterProfiler[21]、org.Hs.eg.db 包,對TC 差 異PFRGs CD44、ANXA1 和NR4A1 的共表達基因進行GO 和KEGG富集分析,以校正P<0.05 且q-value 的值<0.2 為篩選條件,聚焦TC 鐵死亡相關的生物進程(biological process,BP)、細 胞 組 分 (cell component,CC)、分子功能(molecular function,MF)和KEGG 通路,并構建“KEGG-靶點”網絡圖。
1.9 統計學分析采用R 軟件(3.6.3)中的Bioconductor 軟件包進行數據處理和統計學分析。TCGA 數據庫中轉錄組數據不符合正態分布,采用Wilcoxon檢驗評估腫瘤組和對照組之間PFRGs表達的差異,采用單因素和多因素Cox 回歸分析評估PFRGs 表達和臨床數據等與TC 患者OS 之間的相關性。以P<0.05 為差異有統計學意義。
2.1 TC 差異預后基因獲取結果以|log2FC|>1 且校正P<0.05 為篩選條件,鑒 定 出TC-DEGs 共6 773 個,其中3 317個在TC中表達上調,3 456個在TC 中表達下調,繪制火山圖(圖1A)及差異排序圖(圖1B)。以P.cox<0.05 為篩選條件,鑒定出TC-PGs 共1 444 個,其 中HR>1 的 危 險 基 因755 個,HR<1 的保護基因有689 個。取“TCDEGs”與“TC-PGs”兩者的交集,韋恩圖可視化,獲取TC-DPGs 共343 個(圖1 C)。

圖1 TC 差異預后基因篩選結果Fig.1 Screening results of TC differential prognosis genes
2.2 TC 差異PFRGs 篩選結果從GeneCards 數據庫中獲得FRGs 442 個,從FerrDb 數據庫中獲得FRGs 384 個,2 個 數 據 庫 存 在167 個 相 同FRGs,保留唯一值后,共獲取659 個FRGs。取“TCDPGs”與“FRGs”兩者的交集,韋恩圖可視化,獲 取PFRGs 共3 個(圖2),分 別 為CD44、ANXA1 和NR4A1,差異分析及預后分析結果見表1。TC 中CD44 和ANXA1 差異表達上調(|log2FC|>1)是保護性因素(HR<1),NR4A1 差異表達 下 調 (|log2FC| <-1)是 危 險 性 因 素(HR>1)。

表1 TC 差異PFRGs 差異和預后分析Tab.1 Differences and prognostic analysis on TC differential PFRGs

圖2 TC 差異PFRGs 韋恩圖Fig.2 Venn diagram of TC differential PFRGs
2.3 TC 差異PFRGs 表達情況配對或非配對樣本t檢驗結果(圖3A 和3B)顯示:與正常組織比較,TC 差 異PFRGs CD44 和ANXA1 在TC 中 表達上調,NR4A1 在TC 中表達下調。同時HPA 數據庫免疫組織化學圖片結果(圖4)顯示:TC 組織樣本中CD44 和ANXA1 蛋白過表達,NR4A1 蛋白低表達,此結果與TCGA 數據庫結果相互驗證。

圖3 TC 差異PFRGs 的表達Fig.3 Expressions of TC differential PFRGs

圖4 免疫組織化學法檢測TC差異PFRGs蛋白表達(×400)Fig.4 Protein expressions of TC differential PFRGs genes detected by immunohistochemistry(×400)
2.4 TC 差異PFRGs 生存預后情況time-ROC 生存預測結果(圖5A~5C)顯示:TC 患者1、3 和5 年生存預測中,ANXA1 均表現出相對最佳的預測效能,其中5 年的預測效價(AUC=0.200、0.221 和0.198)最佳;CD44 預測5 年生存預后具有最好效能,基因表達在TC 患者早期(1 年)促進死亡發生,在中后期(3 和5 年)促進生存發生;ANXA1 表 達 在TC 患 者1、3 和5 年 中 均 促 進 生 存發生;NR4A1 表達在TC 患者1、3 和5 年中均促進死亡發生。
Kaplan-Meier 曲線生存狀況結果(圖5D~5F)顯示:TC 差異PFRGs CD44、ANXA1 和NR4A1的表達與生存預后相關(P=0.048、0.005 和0.036),其中CD44 和ANXA1 的表達是保護性因素,NR4A1 的表達是危險因素,即低表達CD44 和ANXA1 與高表達NR4A1 具有較差的生存結局。

圖5 TC 差 異PFRGs 的time-ROC 和Kaplan-Meier 曲 線Fig.5 Time-ROC and Kaplan-Meier curves of TC differential PFRGs
2.5 TC 差異PFRGs 獨立預后分析結果對TC差異PFRGsCD44、ANXA1 和NR4A1 的表達進行多因素Cox 回歸分析,計算出3 個基因比例風險回歸模型的Coef 相關系數,即ANXA1(coef)=-0.361,NR4A1(coef)=0.316,CD44(coef)=0.008。將Coef 相關系數與3 個基因的表達值納入風險評分計算公式,計算TC 患者的風險評分。根據患者風險評分中位數(-0.049),將TCGA-THCA 甲狀腺患者分為高風險組255 例和低風險組255 例。
高和低風險組生存差異比較(圖6),Log-rank檢 驗 結 果 顯 示:HR=6.87,95%CI:2.58~18.30,P=0.003;Cox 回歸分析結果顯示:HR=6.88,95%CI:1.56~30.28,P=0.011,高風險組和低風險組患者生存預后比較差異有統計學意義,同時高風險組較低風險組具有更差的生存預后。time-ROC 顯示出風險評分能良好地預測TC患者1、3 和5 年的生存預后,1、3 和5 年的AUC值分別為0.761、0.767 和0.722,其中預測3 年的生存預后效價最優。

圖6 低風險組和高風險組TC 患者Kaplan-Meier 曲線(A)及time-ROC 曲線(B)Fig.6 Kaplan-Meier curve(A) and time-ROC curve(B) of TC patients in low risk group and high risk group
為進一步評估風險評分與其他臨床特征的預后關系,采用單因素和多因素Cox 回歸分析,以單因素Cox 回歸分析的P<0.1 作為納入多因素Cox 回歸分析的條件,對臨床特征進行獨立預后分析(圖7):風險評分是TC 患者生存預后的獨立因素(HR=8.882,95%CI:1.561~50.547,P=0.014),可獨立于其他臨床特征評估TC 患者生存預后。年齡(HR=1.078,95%CI:1.016~1.145,P=0.013)與TC 患者生存預后相關。

圖7 TC 患者臨床特征單因素和多因素Cox 回歸分析森林圖Fig.7 Forest plot of univariate and multivariate Cox regression analysis on clinical characteristics of TC patients
2.6 TC 差異PFRGs Nomogram 圖構建及驗證結果根據多因素Cox 回歸分析結果,以年齡和風險評分作為因子構建Nomogram 圖,預測TC 患者1、3 和5 年的OS 的概率。綜合Nomogram 圖中年齡及風險評分的總體得分,預測TC 患者的生存概率(圖 8A),該 Nomogram 圖 的 C-index=0.938(0.923~0.952),具 有 較 好 的 預 測 能 力。Calibration 分析結果顯示:Nomogram 圖預測TC 患者1、3 和5 年生存概率與實際觀察結果吻合良好,見圖8 B~D。

圖8 TC 差異PFRGs 的Nomogram 圖 和Calibration 圖Fig.8 Nomogram and Calibration charts of TC differential PFRGs
2.7 TC 差異PFRGs 與蛋白共表達相關性分析結果Spearman 相關性分析基因共表達結果顯示:有9 個編碼基因與CD44 表達呈正相關關系(P<0.001),有95 個編碼基因與ANXA1 表達呈正相關關系(P<0.001),有17 個編碼基因與NR4A1 表達呈正相關關系(P<0.001)。分別繪制CD44、ANXA1 和NR4A1 相關系數前5 位的共表達熱圖(圖9):與CD44 最相關的前5 位基因依次是SFTA3 (cor=0.773)、SSBP2 (cor=0.770)、RNF146 (cor=0.755)、IKZF2 (cor=0.732)和ESYT3(cor=0.718);與ANXA1 最相關的前5 位基因依次是DUSP6(cor=0.815)、SDC4(cor=0.812)、TMEM43 (cor=0.803)、BID (cor=0.793)和RAD23B (cor=0.783);與NR4A1 最相關的前5 位基因依次是FOSB (cor=0.889)、NR4A3 (cor=0.845)、NR4A2 (cor=0.840)、CSRNP1(cor=0.837)和ZFP36(cor=0.826)。

圖9 TC 差異PFRGs 共表達熱圖Fig.9 Heat maps of co-expression of TC differential PFRGs
PPI 結果顯示(圖10A):PPI 網絡中共有71 個節點,132 條邊,即71 種蛋白質之間相互作用,共有132 組對應關系。拓撲分析結果顯示:JUN、FOS 和ATF3 等 關 鍵 蛋 白 與CD44、ANXA1 和NR4A1 調控TC 鐵死亡有關(圖10B~D 和表2)。

表2 TC 差異PFRGs 關鍵蛋白拓撲分析Tab.2 Topological analysis on TC differential PFRGs key proteins

圖10 TC 差異PFRGs PPI 網絡圖及拓撲分析圖Fig.10 PPI network and topology analysis diagrams of TC differential PFRGs
2.8 TC 差異PFRGs GO 和KEGG 富集分析結果GO 和KEGG 富集分析結果顯示:在滿足AdjustedP<0.05 且q<0.2 條件下,TC 鐵死亡主要富集于58 個BP、4 個MF 和5 個KEGG,主要與細胞周期的調控有關。聚焦于細胞周期BP、MF和KEGG(圖11A)結果顯示:TC 鐵死亡主要富集于MAPK 活性失活、內肽酶活性的正向調節、凋亡通路的正向調節、ERK1 和ERK2 級聯等生物進程;富集于MAPK 活性、p-MAPK 活性、生長因子結合和轉化生長因子β 結合分子功能;主要富集于MAPK信號通路、Th17細胞分化和細胞凋亡等信號通路。
構建KEGG-靶點網絡圖(圖11 B):TC 鐵死亡主要與MAPK 信號通路中雙特異性磷酸酶(dual specificity phosphatase,DUSP1)1、DUSP5、DUSP6、Fos 原 癌 基 因 (Fos proto-oncogene,FOS)、白細胞介素1 受體輔助蛋白(interleukin 1 receptor accessory protein,IL1RAP)、Jun 原癌基因(Jun proto-oncogene,JUN)、MET 原癌基因(MET proto-oncogene,MET)、Ras 蛋白特異性鳥嘌呤核苷酸釋放因子 1 (Ras protein specific guanine nucleotide releasing factor 1,RASGRF1)、轉化生長因子α(transforming growth factor alpha,TGFA)、轉 化 生 長 因 子β 受 體1 (transforming growth factor beta receptor 1,TGFBR1)和TNF受體超家族成員 1A (TNF receptor superfamily member 1A,TNFRSF1A)分子靶點的調控有關。

圖11 TC 鐵死亡GO 和KEGG 氣泡圖(A)及KEGG-靶點網絡圖(B)Fig.11 TC ferroptosis GO and KEGG bubble diagram(A)and KEGG-target network diagram(B)
TC 是好發于中老年人,以高發病率和低死亡率為特點的內分泌惡性腫瘤[22]。目前,TC 根治術后腫瘤復發和不可切除患者缺乏有效全身治療手段是TC 患者生存預后不良的主要原因[23]。為進一步提升TC 患者生存獲益,突破臨床治療瓶頸,在腫瘤分子生物學研究的基礎上,本研究結果顯示:分子靶向藥物可為碘難治性TC、不可切除性TC 和手術切除后復發TC 患者帶來生存獲益。隨著腫瘤分子生物學研究及分子靶向藥物研發的不斷深入,索拉非尼、樂伐替尼和侖伐替尼相繼問世,在新輔助治療、靶向聯合免疫治療及姑息治療方面極大地提高了TC 患者的生存獲益[24]。鐵死亡是由鐵代謝、脂質代謝和氨基酸代謝共同介導的細胞程序性死亡模式,貫穿TC 發生發展的始終。研究[25]顯示:針對腫瘤細胞鐵死亡的分子靶向藥物(索拉非尼)可明顯延長肝癌、腎癌和食管癌等多種惡性腫瘤患者的生存期。本研究篩選出TC 差異PFRGs,分析其與TC 患者生存預后的關系,并構建預后風險評估系統及預測模型。
本研究通過TCGA、FerrDb 和GeneCards 數據庫,差異分析鑒定出3 317 個在TC 中上調基因和3 456 個在TC 中下調基因,單因素Cox 回歸分析鑒定出343 個差異基因與TC 患者的生存預后相關,其中包括3 個FRGs(CD44、ANXA1 和NR4A1)。結合基因表達譜及臨床病理驗證了CD44、ANXA1在TC 組織中明顯上調,NR4A1 在TC 組織中明顯下調,與差異分析結果相互佐證。CD44 是一種細胞黏附糖蛋白,在腫瘤進展和轉移中發揮作用,被廣泛用于TC 和其他惡性腫瘤的腫瘤干細胞標志物[26-27],LI 等[28]采用雙熒光素酶報告基因及免疫共沉淀實驗證實miR-205-5p 能靶向GGCT,介導TC 細胞中CD44 的表達下調,抑制TC 細胞增殖、遷移和侵襲。ANXA1 是一種磷脂結合的膜定位蛋白,在TC 組織和細胞(SW579 細胞)中表達上調,可通過外泌體介導TC 和甲狀腺濾泡上皮細胞(Nthy-ori3-1 細胞)之間的細胞通訊,促進TC 細胞及甲狀腺濾泡上皮細胞的增殖、侵襲、上皮間質轉化和腫瘤生長[29]。NR4A1 是核受體亞家族4A1組成員轉錄因子,入核調控下游靶基因轉錄,參與腫瘤細胞增殖、細胞凋亡和鐵死亡等多個進程。JIANG 等[30]發 現:NR4A1 直 接 與 LEF1 的 啟 動子區域結合,并導致與組蛋白乙酰化和 DNA 去甲基化的串擾,從而在轉錄上上調LEF1 表達,促進PTC 中下游生長相關基因的表達和腫瘤細胞增殖。
本研究采用Kaplan-Meier 曲線和time-ROC 曲線分析結果顯示:CD44、ANXA1 和NR4A1 的表達與生存預后相關(P=0.048、0.005 和0.036),對TC 患者均具有良好的1、3 和5 年生存預測作用;采用單因素和多因素Cox 回歸分析構建3 基因風險評分系統,計算TC 患者風險評分,結果顯示:風險評分是獨立于TC 其他臨床特征的預后因子(HR=8.882,95%CI:1.561~50.547,P=0.014),風險評分越高,TC 患者1、3 和5 年生存預后越差(AUC=0.761、0.767 和0.722);聯合TC 患者臨床特征構建預測TC 患者1、3 和5 年的生存概率的Nomogram 圖(C-index=0.938),對TC 患者的生存獲益具有較好的預測作用。
本研究中Spearman 相關性分析結果顯示:CD44 與SFTA3、SSBP2、RNF146、IKZF2 和ESYT3 最 相 關;ANXA1 與 DUSP6、SDC4、TMEM43、BID 和RAD23B 最 相 關;NR4A1 與FOSB、NR4A3、NR4A2、CSRNP1 和ZFP36 最相關。聯合3 個TC PFRGs 進行PPI 及拓撲分析結果顯示:TC 鐵死亡與JUN、FOS 和ATF3 等關鍵蛋白的相互作用有關。為研究TC 鐵死亡介導的下游作用機制,本研究進行GO 和KEGG 富集分析,發現TC 細胞鐵死亡主要與其共表達基因(DUSP1、DUSP5、DUSP6、FOS、IL1RAP、JUN、MET、RASGRF1、TGFA、TGFBR1 和TNFRSF1A)介導MAPK 信號通路,影響MAPK活性和p-MAPK 活性等分子功能,調控MAPK 失活有關。MAPK 信號通路與活性氧(reactive oxygen species,ROS)誘導的細胞鐵死亡有關,細胞內過量的ROS 以脂質過氧化的形式介導鐵死亡[31],參與多種腫瘤細胞的增殖、遷移和侵襲過程[32-33]。
綜上所述,本研究篩選得到的差異PFRGs(CD44、ANXA1 和NR4A1)與TC 發生發展和預后相關;成功構建了3 個基因聯合的風險評分系統及預測模型,結合TC 患者臨床特征能有效預測TC 患者1、3 和5 年生存期,TC 鐵死亡的作用機制可能與介導共表達基因、激活MAPK 信號通路和調控相關分子功能有關。