基于生物信息學(xué)分析和實驗驗證篩選影響肝癌患者預(yù)后的鐵死亡調(diào)控基因

姜宇朗，王 錚，孫明瑜

（1.上海中醫(yī)藥大學(xué)附屬曙光醫(yī)院，上海中醫(yī)藥大學(xué)肝病研究所肝腎疾病病證教育部重點(diǎn)實驗室，上海 201203；2.上海中醫(yī)藥大學(xué)，上海 201203）

原發(fā)性肝癌是目前我國發(fā)病率第四和死亡率第二的惡性腫瘤，嚴(yán)重威脅著我國人民的生命和健康安全［1］。由于肝癌起病并無特異性的臨床癥狀，大部分患者在確診時往往已經(jīng)處于中晚期，喪失了最佳的手術(shù)治療機(jī)會［2］，導(dǎo)致五年生存率較差，因此早診斷早治療是肝癌防治的一大目標(biāo)，對于肝癌診斷預(yù)后具有良好敏感性和特異性的預(yù)后生物學(xué)標(biāo)志物仍亟待開發(fā)［3］。

鐵死亡是近年來發(fā)現(xiàn)的一種依賴于鐵離子的調(diào)控且有別于凋亡、焦亡、程序性壞死、自噬的新型調(diào)控性細(xì)胞死亡［4］，其主要的細(xì)胞形態(tài)學(xué)變化特點(diǎn)是線粒體的改變、包括線粒體的皺縮、膜密度的增高、嵴的減少等，生化學(xué)方面改變主要包括脂質(zhì)過氧化物的蓄積，活性氧（reactive oxygen species，ROS）的增多，最終對磷脂膜造成損傷并使細(xì)胞死亡［4，5］。鐵死亡作為一種代謝性細(xì)胞死亡方式受到多種代謝途徑的調(diào)控，其中包括氧化應(yīng)激、鐵離子代謝、線粒體能量代謝和胱氨酸代謝、多不飽和脂肪酸的合成以及葡萄糖的代謝密切相關(guān)［5，6］。研究表明鐵死亡在肝癌的發(fā)生、發(fā)展中都扮演著重要的角色，誘導(dǎo)肝癌細(xì)胞鐵死亡有望發(fā)展為肝癌治療領(lǐng)域的一顆新星［7，8］。

研究發(fā)現(xiàn)，鐵死亡相關(guān)基因在癌癥的預(yù)后情況和預(yù)防診斷中也存在著重大的價值，有研究發(fā)現(xiàn)17個鐵死亡相關(guān)的長鏈非編碼RNA 可以用來預(yù)測結(jié)腸癌的轉(zhuǎn)歸［9］，8 個在胃腺癌患者中差異表達(dá)的預(yù)后FRGs 可以預(yù)測胃腺癌患者的生存預(yù)后狀態(tài)［11］。但是診斷肝癌預(yù)后轉(zhuǎn)歸的差異鐵死亡基因仍然未被有效發(fā)掘，本文旨在通過生物信息學(xué)和實驗驗證的方式，找出影響肝癌患者預(yù)后的差異FRGs，從而為臨床肝癌預(yù)后新的生物標(biāo)志物的開發(fā)提供線索，并為肝癌鐵死亡療法發(fā)掘新的治療效果判斷的標(biāo)志物。本研究的流程和具體方法見圖1。

圖1 生物信息學(xué)分析的工作流程示意圖Fig 1 Schematic diagram of the workflow of bioinformatics analysis

1 材料和方法

1.1 基因表達(dá)量的下載和鐵死亡相關(guān)基因的表達(dá)量的整理

肝癌患者的基因測序結(jié)果從癌癥基因組圖譜計劃（The Cancer Genome Atlas， TCGA）（https：//portal.gdc.cancer.gov/）數(shù)據(jù)庫下載獲得，結(jié)果共獲得422 個數(shù)據(jù)集，其中包括374 例肝癌患者和50 例正常樣本。其中鐵死亡相關(guān)基因從Ferrdb 數(shù)據(jù)庫（http：//www.zhounan.org/ferrdb/）獲取得到。將FRGs 和TCGA 數(shù)據(jù)庫基因表達(dá)相關(guān)聯(lián)，得到FRGs 在每個樣本中的表達(dá)數(shù)據(jù)。

使用limma 包將肝癌組織和正常組織的基因表達(dá)量進(jìn)行非參數(shù)檢驗，設(shè)定差異條件為錯誤發(fā)現(xiàn)率（false discovery rate， FDR）＜0.05 并 且| log2Fold-Change|＞1，獲得在肝癌樣本和正常組織中差異表達(dá)的相關(guān)FRGs。

通過綜合分析患者的生存時間和生存狀況，并結(jié)合FRGs 的差異表達(dá)，使用survival 包進(jìn)行單因素COX 回歸分析，從而確定與肝癌患者預(yù)后有關(guān)的FRGs（P＜0.05）。

1.2 影響肝癌患者預(yù)后的鐵死亡相關(guān)差異表達(dá)基因的確定和PPI 網(wǎng)絡(luò)的構(gòu)建

首先使用venn 包來篩選出差異表達(dá)的FRGs和預(yù)后生存相關(guān)的FRGs 的共同部分，這樣可以獲得影響肝癌患者預(yù)后的差異FRGs。然后，使用pheatmap 包繪制出這些差異表達(dá)的預(yù)后FRGs 的熱圖。最后，通過survival 包繪制相應(yīng)的森林圖。將上述基因?qū)氲鞍谆プ骶W(wǎng)絡(luò)（protein protein interation network， PPI）工具（http；//string-db.org/）可視化蛋白之間的互作關(guān)系。為進(jìn)一步體現(xiàn)蛋白質(zhì)之間的表達(dá)關(guān)聯(lián)性，計算出核心蛋白的表達(dá)相關(guān)性并使用igraph 包繪制蛋白質(zhì)互作網(wǎng)絡(luò)圖。

1.3 風(fēng)險模型的構(gòu)建及其驗證

利用glmnet 包中的LASSO 回歸比例風(fēng)險模型，構(gòu)建模型并獲取基因風(fēng)險系數(shù)，計算風(fēng)險評分。根據(jù)風(fēng)險值中位數(shù)，將樣本分為高低風(fēng)險組。通過survival 和survminer 包，比較兩組生存和預(yù)后差異。使用timeROC 包繪制時間依賴ROC 曲線，評估模型預(yù)測準(zhǔn)確性，利用pheatmap 包繪制風(fēng)險曲線。

1.4 PCA 和tSNE 降 維 分 析

使用Rstne 包進(jìn)行主成分（principal component analysis， PCA）和 T 分 布 隨 機(jī) 鄰 嵌 入 分 析（t-distributed stochastic neighbor embedding， tSNE）觀察患者根據(jù)高低風(fēng)險模型的分組情況的分布。

1.5 獨(dú)立預(yù)后因素分析

將患者的腫瘤分期、腫瘤等級、年齡、性別、風(fēng)險評分使用survival 包進(jìn)行單因素和多因素COX 回歸分析，獨(dú)立預(yù)測因素分析，并繪制單因素和多因素對生存時間和生存狀態(tài)的影響的森林圖。

1.6 風(fēng)險差異基因分析

使用limma 包設(shè)定錯誤條件發(fā)現(xiàn)率為FDR 值小于0.05 且| log2FoldChange|＞1，將風(fēng)險值分別屬于高低組的人群進(jìn)行基因差異表達(dá)分析，得到各個樣本在TCGA 數(shù)據(jù)庫中高低風(fēng)險組中的基因的平均表達(dá)量。

1.7 GO 和KEGG 富集分析和免疫相關(guān)分析

使用enrich-plot 包將風(fēng)險差異基因進(jìn)行GO 和KEGG 富集分析，借助ggplot 包使用氣泡圖將結(jié)果可視化。使用GESAbase 和GSVA 包進(jìn)行免疫細(xì)胞和功能的免疫浸潤相關(guān)分析，并用ggplot 包繪制高低風(fēng)險組患者的箱線圖。

1.8 差異的預(yù)后鐵死亡基因的免疫組化驗證

在人類蛋白質(zhì)圖譜數(shù)據(jù)庫（https：//www.proteinatlas.org/）下載TOP20 基因在正常肝組織和肝癌組織的高分辨率免疫組化結(jié)果圖。

1.9 差異的預(yù)后鐵死亡基因的細(xì)胞層面驗證

裂解正常肝細(xì)胞L02 和肝癌細(xì)胞HepG2，提取細(xì)胞總RNA，測定其濃度和純度后，將RNA 逆轉(zhuǎn)錄為cDNA，以cDNA 為模板進(jìn)行qPCR 擴(kuò)增。以正常組作為對照，采用RQ 值來顯示目的基因mRNA表達(dá)水平。比較正常肝細(xì)胞和肝癌細(xì)胞的目的基因的表達(dá)量。

1.10 差異的預(yù)后鐵死亡基因組織層面驗證

選取來自曙光醫(yī)院的六名臨床確診為肝細(xì)胞癌的患者，手術(shù)切除樣本組織后稱取肝癌和癌旁組織各50 mg，勻漿機(jī)勻漿， 提取組織總RNA， 進(jìn)行后續(xù)PCR 檢測。比較正常肝組織和癌旁組織目的基因的表達(dá)量。

1.11 統(tǒng)計學(xué)方法

數(shù)據(jù)用GraphPad Prism 9 軟件進(jìn)行統(tǒng)計，兩組間計量資料符合正態(tài)分布且方差齊時，選用用t檢驗，多組間比較用One-way-ANOVA 單因素方差分析，進(jìn)一步多組間比較采用Tukey 檢驗。若數(shù)據(jù)不符合正態(tài)分布則采用非參數(shù)秩和檢驗。所有數(shù)據(jù)以（±s）表示，P＜0.05 為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 差異的肝癌預(yù)后FRGs 的篩選結(jié)果

從TCGA 數(shù)據(jù)庫共獲得374 例肝癌樣本和50例正常樣本的RNA-seq 數(shù)據(jù)和對應(yīng)的臨床數(shù)據(jù)集。從ferrdb 數(shù)據(jù)庫上共獲得FRGs249 個，通過分析得到84 個差異表達(dá)的FRGs 和90 個預(yù)后相關(guān)的FRGs，對兩者取交集可以得到42 個和肝癌預(yù)后相關(guān)的差異表達(dá)的FRGs（圖2A）。依據(jù)上述42 個基因表達(dá)數(shù)據(jù)繪制熱圖（圖2C）。同時篩選出P值排名靠前的20 個基因，其風(fēng)險比率（hazard ratio，HR）及其95%置信區(qū)間如表1 所示，同時繪制森林圖（圖2B），確定肝癌的風(fēng)險基因，尤其以ZFP69B、EIFS21、ATG3、ATG7等鐵死亡相關(guān)基因的致病風(fēng)險效應(yīng)最為明顯。

表1 單因素COX 回歸分析得出的與肝癌預(yù)后相關(guān)的FGRs（P值排名前20）Tab 1 FGRs related to liver cancer prognosis obtained from univariate COX regression analysis（top 20 P value）

圖2 篩選出的影響肝癌預(yù)后的FRGsFig 2 Screened FRGs affecting prognosis in liver cancer

2.2 蛋白質(zhì)相互作用網(wǎng)絡(luò)的構(gòu)建結(jié)果

將42 個交集與肝癌預(yù)后相關(guān)的差異FRGs 導(dǎo)入STRING 數(shù)據(jù)庫，剔除游離靶點(diǎn)，得到基于String數(shù)據(jù)庫的蛋白質(zhì)互作網(wǎng)絡(luò)（圖3A），其中相關(guān)性最高的核心靶蛋白主要有MAPK、NRAS、SRC、HRAS和CDKN2A。蛋白相關(guān)性作用網(wǎng)絡(luò)（圖3B）結(jié)果顯示，ALB 與其他蛋白之間呈現(xiàn)出負(fù)向調(diào)控的關(guān)系，而其余基因則基本表現(xiàn)出正相關(guān)。

圖3 蛋白互作網(wǎng)絡(luò)的結(jié)果示意圖Fig 3 Schematic representation of the results of the protein interaction network

2.3 風(fēng)險預(yù)后模型的建立以及檢驗

LASSO 回歸分析得到肝癌預(yù)后的風(fēng)險模型，最 終 獲 得11 個 基 因（G6PD、HRAS、SLC1A5、MT3、SRXN1、SLC7A11、ZFP69B、SLC2A1、STMN1、RRM2、SQSTM1）。根據(jù)風(fēng)險值的中位數(shù)，將患者分為高、低風(fēng)險組，結(jié)果顯示高風(fēng)險組患者的生存率和生存時間顯著低于高風(fēng)險組患者（P＜0.01）（圖4A）?；颊呷晟媛实腞OC 曲線面積均大于0.6，說明本研究建立的風(fēng)險模型可信（圖4B）。同時根據(jù)患者的風(fēng)險值及其分組繪制患者人數(shù)和風(fēng)險值的關(guān)系以及存活時間和患者人數(shù)的關(guān)系（圖4C，D），可以看出，隨著風(fēng)險值的增加病人的死亡數(shù)明顯增加。為了進(jìn)一步研究高低風(fēng)險組人群的分布情況，使用PCA（圖4E）和t-SNE（圖4F）可視化降維分析。最后再進(jìn)行單因素（圖4G）和多因素（圖4H）COX 回歸分析來檢測能單獨(dú)影響肝癌患者預(yù)后的因素，研究發(fā)現(xiàn)，肝癌分期和風(fēng)險評分可以作為肝癌患者預(yù)后的獨(dú)立預(yù)測因子（P＜0.01）。

圖4 影響肝癌患者預(yù)后的模型的構(gòu)建及效能檢驗和肝癌患者獨(dú)立預(yù)后因素的分析Fig 4 Construction and efficacy test of the model affecting the prognosis of liver cancer patients with independent prognostic factors

2.4 風(fēng)險差異基因的GO、KEGG 富集分析結(jié)果

使用limma 包獲得差異FRGs 在高低風(fēng)險組的患者的平均表達(dá)量，并做差異分析，最終得到風(fēng)險差異FRGs，GO 分析表明風(fēng)險預(yù)后FRGs 主要和體液免疫應(yīng)答、細(xì)胞吞噬作用、淋巴細(xì)胞介導(dǎo)免疫等生物學(xué)過程有關(guān)，KEGG 富集分析顯示風(fēng)險差異基因參與包括細(xì)胞周期、PI3K/AKT 通路、脂質(zhì)氧化、癌癥中的蛋白多糖的表達(dá)等（圖5A，B），聚類分析結(jié)果顯示，這些肝癌的風(fēng)險FRGs 和免疫的關(guān)系極為密切，這些基因可能在調(diào)控腫瘤的免疫微環(huán)境中起作用。

圖5 差異預(yù)后FRGs 的Go 和KEGG 富集分析結(jié)果氣泡圖Fig 5 Bubble plots of Go and KEGG enrichment analysis results for differential prognostic FRGs

2.5 免疫浸潤分析

進(jìn)一步為了探究差異預(yù)后FRGs 和免疫功能的關(guān)系，通過單樣本基因集富集分析（ssGSEA）確定樣本高風(fēng)險組和低風(fēng)險組的免疫細(xì)胞浸潤評分和免疫相關(guān)功能評分（圖6）。相關(guān)評分顯示，高低風(fēng)險組患者的aCDs、Macrophages、NK 細(xì)胞表達(dá)差異具有統(tǒng)計學(xué)意義，在免疫功能方面主要和抗原呈遞細(xì)胞的激活、免疫檢查點(diǎn)、MHC 分子、副炎癥、2 型干擾素的應(yīng)答有關(guān)。

圖6 高低風(fēng)險組患者的免疫細(xì)胞和免疫功評分Fig 6 Immune cell and immunocompetence scores of patients in high and low risk groups

2.6 正常組織和肝癌組織的差異預(yù)后FRGs 免疫組化結(jié)果

從人類蛋白質(zhì)圖譜數(shù)據(jù)庫中下載得到了16 個基因在肝組織和肝癌組織中的免疫組化結(jié)果（圖7）（部分基因在數(shù)據(jù)庫中并無蛋白質(zhì)表達(dá)數(shù)據(jù)）。分析目的基因在兩者之間的相對表達(dá)量發(fā)現(xiàn)G6PD、PRDX1、TXNRD1、ATG3 在肝癌組織中表達(dá)明顯上 調(diào)，而SLC1A5、NRAS、EIF2S1、SLC38A1、ATG7在肝癌組織中表達(dá)明顯下調(diào)，而其余基因表達(dá)量變化并無明顯差異。

圖7 TOP20 基因在正常肝組織和肝癌組織中的相對表達(dá)的免疫組化結(jié)果圖Fig 7 Immunohistochemical results of TOP20 gene expression in normal and hepatocellular carcinoma tissues.

2.7 肝細(xì)胞與肝癌細(xì)胞的差異預(yù)后FRGs 的基因表達(dá)結(jié)果

為進(jìn)一步驗證上述生物信息分析結(jié)果的可靠性，在肝細(xì)胞和肝癌細(xì)胞中進(jìn)行PCR 擴(kuò)增反應(yīng)后，和肝細(xì)胞比較肝癌細(xì)胞中部分差異FRGs 發(fā)生顯著改變（圖8），肝癌細(xì)胞中G6PD、STMN1、SLC7A1 1、TXNRD1表 達(dá) 水 平 顯 著 增 加（P＜0.05），而SLC1A5、ATG3水平顯著下調(diào)（P＜0.05）。

圖8 基因在肝細(xì)胞和肝癌細(xì)胞差異預(yù)后FRGs 的相對基因表達(dá)量Fig 8 Relative gene expression of TOP20 gene in differentially prognostic FRGs of hepatocytes and hepatocellular carcinoma cells

2.8 肝癌和癌旁組織的差異預(yù)后FRGs 的相對表達(dá)結(jié)量

和癌旁組織比較，肝癌組織的G6PD、SLC7A1 1、TXNRD1水平顯著增加（P＜0.05），而SLC2A1、SLC1A5、ATG3表達(dá)水平顯著降低（P＜0.05）。這與細(xì)胞實驗的結(jié)果基本一致（圖9）。

圖9 TOP20 基因在肝癌和癌旁組織差異預(yù)后FRGs 的相對表達(dá)量Fig 9 Relative expression of TOP20 gene in differentially prognostic FRGs of hepatocellular carcinoma and Paracarcinomatous tissues

3 討論

基于生物信息學(xué)的分析結(jié)果發(fā)現(xiàn)G6PD、ZFP69B、SLC2A1、SLC1A5、HILPDA、MAFG、STMN1、SRXN1、SLC7A11、NRAS在肝癌預(yù)后的風(fēng)險比當(dāng)中占比較高。意味著這些基因有進(jìn)一步開發(fā)為肝癌預(yù)后的標(biāo)志物的價值。根據(jù)多個差異表達(dá)的FRGs 基因構(gòu)建的風(fēng)險預(yù)后模型可以較好的將患者區(qū)分為高風(fēng)險和低風(fēng)險人群，為肝癌的診斷預(yù)防提供參考，這些基因也可能是預(yù)測靶向鐵死亡療法治療肝癌的有效性預(yù)測的潛在生物標(biāo)志物［12］。進(jìn)一步的結(jié)合免疫組化和細(xì)胞與組織的PCR 實驗結(jié)果發(fā)現(xiàn)其中G6PD、SLC1A5、TXNRD1的變化趨勢和生物信息學(xué)的推測結(jié)果一致。這提示在這些差異預(yù)后的FRGs 當(dāng)中，這些靶標(biāo)可能位于整個肝癌鐵死亡調(diào)控網(wǎng)絡(luò)的核心位置。

其中G6PD是葡萄糖-6-磷酸脫氫酶的編碼基因，其主要功能是產(chǎn)生 NADPH，而NADPH 是抵御氧化劑和還原性生物合成反應(yīng)中的關(guān)鍵電子供體［13］。G6PD 是葡萄糖代謝的關(guān)鍵酶，在肝癌患者中G6PD 往往處于高水平，這也意味著肝癌中葡萄糖高代謝狀態(tài)，而葡萄糖高代謝又會誘導(dǎo)鐵死亡的發(fā)生，G6PD也可以作為肝癌患者預(yù)后的獨(dú)立風(fēng)險因素來對患者進(jìn)行定期監(jiān)測［14］。因此G6PD有成為肝癌鐵死亡預(yù)后基因的潛力，但仍然需要進(jìn)一步基礎(chǔ)實驗和臨床實驗佐證。

SLC1A5作為人溶質(zhì)載體家族的一員，是一種鈉離子依賴性的中性氨基酸的轉(zhuǎn)運(yùn)體，可以轉(zhuǎn)運(yùn)包括和鐵死亡密切相關(guān)的谷氨酰胺、甘氨酸、谷氨酸、異亮氨酸等［15］。肝癌細(xì)胞氨基酸高代謝狀態(tài)是腫瘤高侵襲力的體現(xiàn)之一，代謝重編程是腫瘤細(xì)胞的特征性印記之一［16，17］。靶向SLC1A5可以阻止谷氨酰胺進(jìn)入細(xì)胞，從而降低內(nèi)腫瘤細(xì)胞內(nèi)谷氨酰胺介導(dǎo)的高新陳代謝狀態(tài)，在一定程度上“餓死”腫瘤以達(dá)到癌癥治療的目的［18，19］。

在經(jīng)典的鐵死亡防御體系中，胱氨酸通過細(xì)胞內(nèi)唯一的胱氨酸/谷氨酸反向轉(zhuǎn)運(yùn)體進(jìn)入細(xì)胞，隨后在TXNRD1 酶的作用下迅速被還原為半胱氨酸，其作為谷胱甘肽的底物之一促進(jìn)隨后的谷胱甘肽合成［20，21］。最后GPX4利用谷胱甘肽作為底物來還原胞內(nèi)過多的脂質(zhì)過氧化物，從而抑制鐵死亡的發(fā)生。因此TXNRD1可以間接反應(yīng)胞內(nèi)GSH 水平，也能預(yù)測患者對于鐵死亡的敏感性，從而評估鐵死亡藥物的臨床療效［22］。

篩選出的三個基因都和腫瘤的代謝密切相關(guān)，尤其是氨基酸的吸收、分布、轉(zhuǎn)化方面，而氨基酸代謝作為鐵死亡三大組成部分之一和鐵死亡的關(guān)系不言而喻。諸多氨基酸的代謝都參與了鐵死亡的調(diào)控，而其中可能可以進(jìn)一步發(fā)掘出誘導(dǎo)鐵死亡治療疾病的新靶點(diǎn)。但是本文并未進(jìn)一步研究這三者之間內(nèi)在的調(diào)控關(guān)系，它們之間是否存在某些內(nèi)在或者外在的關(guān)聯(lián)性需要進(jìn)一步探究。

總而言之，本研究篩選出的和肝癌相關(guān)FRGs可能用來預(yù)測常規(guī)病理因素之外的肝癌預(yù)后，與鐵死亡相關(guān)的這些特異的生物標(biāo)志物也有進(jìn)一步發(fā)展成為肝癌的診斷性指標(biāo)的可能，同時針對這些靶點(diǎn)的鐵死亡途徑也可能為肝癌的治療提供新的方向。鑒于鐵死亡和免疫的關(guān)系密切，未來癌癥治療的研究重點(diǎn)應(yīng)放在免疫阻斷劑和靶向治療結(jié)合誘導(dǎo)細(xì)胞鐵死亡［23，24］。

作者貢獻(xiàn)度說明：

姜宇朗設(shè)計文章的思路與框架、撰寫論文；王錚參與生信分析；孫明瑜指導(dǎo)論文撰寫并審校。

所有作者聲明不存在利益沖突關(guān)系。