FeRGs生物信息學分析及肺腺癌預后模型構(gòu)建

楊彩珍,何杰,柳廣南

(1 廣西醫(yī)科大學,廣西 南寧 530021； 2 廣西醫(yī)科大學第二附屬醫(yī)院)

肺腺癌(LUAD)約占肺癌病例的40%，其起病隱匿，多數(shù)病人檢出時往往已處于晚期，總生存期低于5年[1-2]。LUAD病因復雜且具有異質(zhì)性[3]，其預后差且難以預測，因此需要一個預測其預后的模型。鐵死亡是新近發(fā)現(xiàn)的一種細胞程序性死亡過程，機制可能主要與脂質(zhì)及脂質(zhì)過氧化物蓄積等有關[4]，誘發(fā)癌細胞鐵死亡成為一種新的治療方法[5]。有研究表明，鐵死亡在LUAD的發(fā)生發(fā)展中起重要作用，鐵死亡相關基因(FeRGs)參與了LUAD的發(fā)生[6]。然而，F(xiàn)eRGs是否與LUAD的預后相關尚不清楚。本文旨在用生物信息學方法研究FeRGs對LUAD的預后價值并建立預后模型，為LUAD病人的預后評估提供新方法。

1 材料與方法

1.1 數(shù)據(jù)來源

從公共數(shù)據(jù)庫癌癥基因組圖譜(https://portal.gdc.cancer.gov)下載585例病人的RNA-seq表達矩陣及其臨床數(shù)據(jù)，數(shù)據(jù)截止時間為2020年9月。從相關文獻中收集到60個FeRGs[7-10]。

1.2 預后模型的構(gòu)建

應用R語言中的“l(fā)imma”包[11]篩選腫瘤組織與癌旁組織的差異表達FeRGs，篩選條件為偽發(fā)現(xiàn)率(FDR)<0.05。應用單因素Cox回歸分析篩選出與預后相關的FeRGs。取差異表達FeRGs和預后相關的FeRGs的交集，得到預后相關的差異表達FeRGs。為了最大限度地降低過度擬合的風險，采用“glmnet”包[11]中的LASSO回歸分析構(gòu)建預后模型。模型的懲罰參數(shù)(λ)根據(jù)十折交叉驗證來確定。根據(jù)每個基因的表達水平及其相應的回歸系數(shù)計算病人的風險評分，以評分的中位數(shù)為界，將病人分為高風險組和低風險組。采用“stats”和“Rtsne”包[12]進行主成分分析(PCA)和t-分布隨機鄰近插入(t-SNE)分析，以評估模型的病人分類能力。用單因素和多因素Cox回歸驗證風險評分的獨立性。

1.3 功能富集分析

用“clusterProfiler”包[11]對高低風險組的差異表達基因進行KEGG通路富集分析和基因本體(GO)分析，以了解這些差異表達基因的作用位點和信號通路。用“GSVA”和“GSEABase”包[11]進行單樣本基因集富集分析(ssGSEA)，分析高低風險組免疫細胞浸潤情況、免疫相關通路活性的差異。

1.4 統(tǒng)計學分析

應用R軟件(版本4.0.2)進行統(tǒng)計學分析。不同組織基因表達水平的比較采用非參數(shù)秩和檢驗。用Mann-Whitney檢驗比較組間免疫細胞以及免疫通路的ssGSEA結(jié)果。用對數(shù)秩檢驗和生存曲線分析比較兩個風險組的生存率。以P<0.05為差異有統(tǒng)計學意義。

2 結(jié) 果

2.1 預后相關的差異表達基因

本文60個FeRGs中，有45個(75%)基因存在差異表達，高、低表達分別為11、34個。單因素Cox回歸分析顯示，60個FeRGs中有15個與預后相關(圖1)。將45個差異表達的FeRGs與15個預后相關的FeRGs取交集后得到14個預后相關的差異表達FeRGs，見表1和圖2。其相關性網(wǎng)絡以及蛋白質(zhì)相互作用網(wǎng)絡結(jié)果見圖3。

2.2 構(gòu)建的預后模型

應用LASSO回歸分析從上述14個基因中篩選出11個(ALOX15、ATP5MC3、CISD1、DPP4、FANCD2、GLS2、GSS、PHKG2、ACSL3、PEBP1、PGD)建立預后模型，基因篩選的過程見圖4。風險評分的計算公式如下:風險評分=(-0.054)×ALOX15表達量+0.050×ATP5MC3表達量+0.204×CISD1表達量+(-0.021)×DPP4表達量+0.070×FANCD2表達量+(-0.364)×GLS2表達量+0.157×GSS表達量+0.179×ACSL3表達量+(-0.188)×PHKG2表達量+(-0.346)×PEBP1表達量+0.038×PGD表達量。模型的回歸系數(shù)及與預后相關程度見表2。以風險評分的中位數(shù)為界，將病人分為高風險組、低風險組(圖5A)，風險評分越大，病人的預后越差(圖5B)。PCA分析和t-SNE分析顯示，高、低風險組的病人分布在兩個方向，說明本文建立模型能很好地區(qū)分兩個分組的病人(圖5C、D)。Kaplan-Meier曲線分析顯示，兩組生存曲線存在差異，低風險組生存率較高(χ2=25.734，P<0.05)(圖5E)。ROC曲線分析顯示，1、2、3年曲線下面積(AUC)分別為0.688、0.697、0.691(圖5F)。

圖1 15個與LUAD預后相關的FeRGs森林圖

表1 14個預后相關的差異表達FeRGs的差異分析

A：14個預后相關的差異表達FeRGs維恩圖；B：14個預后相關的差異表達FeRGs森林圖；C:14個預后相關的差異表達FeRGs熱圖，N為癌旁組織，T為腫瘤組織。

A：基因相關性網(wǎng)絡圖，連線代表基因之間存在共表達的關系，紅色線代表正相關，藍色線代表負相關；B：蛋白質(zhì)相互作用網(wǎng)絡圖，玻璃球代表基因，連線的顏色代表不同蛋白存在相互作用的證據(jù)來源和類型。淡藍色線代表證據(jù)來自現(xiàn)有的數(shù)據(jù)庫，粉紅色線代表證據(jù)來自實驗驗證，青色線代表證據(jù)來自文本挖掘，黑色線代表蛋白之間存在共表達關系。

A：根據(jù)十折交叉驗證結(jié)果選擇最優(yōu)λ值，在λ值最小時篩選出11個基因；B：LASSO回歸系數(shù)分布圖，λ值取最小值時所對應的基因數(shù)為11個。

表2 11個模型基因的回歸系數(shù)及與預后相關的程度

2.3 模型的獨立預后分析

單因素Cox分析顯示，風險評分與LUAD的預后有明顯的聯(lián)系(HR=3.202,95%CI=2.169～4.726,P<0.001)；多因素Cox分析顯示，風險評分為LUAD總生存率的獨立預測因子(HR=2.886,95%CI=1.943～4.285,P<0.001)。見圖6。

2.4 功能富集分析

GO分析顯示，高、低風險組的差異表達基因參與了染色體分離、免疫反應和多肽酶活性調(diào)節(jié)等生物過程。KEGG結(jié)果顯示，這些差異表達基因主要富集在在鐵死亡、IL-17信號通路和脂肪酸代謝等通路(圖7)。ssGSEA結(jié)果顯示，高風險組與低風險組參與抗原呈遞過程的細胞組分明顯不同，巨噬細胞和自然殺傷細胞在高風險組的浸潤程度高于低風險組(z=-3.671、5.074，P<0.01)，而肥大細胞的浸潤程度則在低風險組更高(z=-5.318，P<0.01)(圖8A)。免疫功能分析顯示，高低風險組之間的炎癥促進、副炎癥作用也明顯不同(z=-2.736、-3.213，P<0.05)(圖8B)。

A：兩個風險組風險評分分布曲線；B：風險評分和生存狀態(tài)散點圖；C：兩個風險組PCA分析；D：兩個風險組t-SNE分析；E：高低風險組生存曲線對比圖；F：3年總生存率ROC曲線及AUC。

A：單因素Cox分析初步篩選LUAD預后的預測指標；B：多因素Cox分析進一步篩選預測指標，風險評分可作為LUAD的獨立預測因子。

A：GO分析；B：KEGG分析。

A：免疫細胞浸潤分析結(jié)果；B：免疫相關通路分析結(jié)果。***P<0.001；**P<0.01;*P<0.05;ns表示P>0.05。

3 討 論

相關統(tǒng)計顯示，2015年國內(nèi)肺癌發(fā)病率、死亡率均排首位[13]。LUAD在肺癌中占比最大，同時具有異質(zhì)性和預后差的特點。鐵死亡為新的研究熱點，但是關于FeRGs在LUAD預后中作用的研究較少。本文研究建立一個基于FeRGs的新型預后模型以評估LUAD的預后。

本文研究通過閱讀文獻獲得了60個FeRGs，用生物信息學的方法分析顯示，有45個(75%)基因在LUAD組織和癌旁組織中存在差異表達，說明鐵死亡與LUAD密切相關，具有進一步研究的價值。應用LASSO回歸分析在14個預后相關的差異表達FeRGs中篩選出11個基因用于建立模型，包括ALOX15、ATP5MC3、CISD1、DPP4、FANCD2、GLS2、GSS、PHKG2、ACSL3、PEBP1和PGD。這11個基因分別為鐵代謝相關的基因(CISD1、FANCD2、PHKG2等)、與脂質(zhì)代謝相關的基因(PEBP1、ACSL3、ALOX15、DPP4等)、與能量代謝相關的基因(PGD、GLS2)、參與谷胱甘肽合成的基因(GSS)[14]，以及目前需要進一步進行研究的基因(ATP5MC3)[8]。通過查閱文獻，發(fā)現(xiàn)這些基因直接或者間接參與了鐵死亡過程的調(diào)節(jié)。CISD1的表達水平降低導致線粒體中鐵的積累，可誘導肝癌細胞發(fā)生鐵死亡[15]；FANCD2參與DNA損傷修復，對鐵死亡具有負調(diào)節(jié)作用[16]；PHKG2過表達可以引起鐵的積累，從而引起鐵死亡的發(fā)生[17]；PEBP1可與15-脂氧合酶結(jié)合形成復合物導致過氧氫多不飽和磷脂酰乙醇胺的生成，而過氧氫多不飽和磷脂酰乙醇胺的堆積可以引起細胞發(fā)生鐵死亡[18]；ACSL3可將外源單不飽和脂肪酸轉(zhuǎn)化為脂肪?；o酶A，促使細胞對鐵死亡產(chǎn)生抵抗作用[19]；ALOX15的表達對鐵死亡具有正調(diào)節(jié)作用；降低DPP4的活性可降低細胞鐵死亡的發(fā)生[20]；PGD通過調(diào)控PI3K/Akt通路誘導非小細胞肺癌細胞發(fā)生鐵死亡；GLS2的高表達可以促進胃癌細胞發(fā)生鐵死亡；GSS表達降低可以間接引起細胞發(fā)生鐵死亡[21]。本文建立模型的基因中有4個對鐵死亡有負調(diào)節(jié)作用(CISD1、FANCD2、ACSL3、GSS)，6個基因(ALOX15、PHKG2、PEBP1、PGD、GLS2、DPP4)可以促進鐵死亡的發(fā)生。本文進一步分析發(fā)現(xiàn)，基因多數(shù)存在相互作用，例如PGD與GSS、GLS2、PEBP1等存在相互作用。同時，本文建立的模型具有良好的風險病人分類功能，可區(qū)分高風險組和低風險組的病人，而且預測性能較好；Cox回歸分析顯示，風險評分與預后相關，且具預測獨立性；通過ssGSEA分析發(fā)現(xiàn)，高、低風險病人的預后存在差異。其原因可能為：①高、低風險組脂肪酸代謝、鐵死亡、免疫反應和IL-17信號通路等生物過程和通路存在差異；②高、低風險組的免疫細胞浸潤類型和程度存在差異，免疫通路活性如炎癥通路的活性、副炎癥通路的活性也存在差異。

綜上所述，本文構(gòu)建了一個由11個FeRGs組成的新型LUAD預后模型，該模型為LUAD預后的評估提供了參考。