楊彩珍,何杰,柳廣南

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

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

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

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

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

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

1.3 功能富集分析

用“clusterProfiler”包[11]對(duì)高低風(fēng)險(xiǎn)組的差異表達(dá)基因進(jìn)行KEGG通路富集分析和基因本體(GO)分析，以了解這些差異表達(dá)基因的作用位點(diǎn)和信號(hào)通路。用“GSVA”和“GSEABase”包[11]進(jìn)行單樣本基因集富集分析(ssGSEA)，分析高低風(fēng)險(xiǎn)組免疫細(xì)胞浸潤(rùn)情況、免疫相關(guān)通路活性的差異。

1.4 統(tǒng)計(jì)學(xué)分析

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

2 結(jié) 果

2.1 預(yù)后相關(guān)的差異表達(dá)基因

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

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

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

圖1 15個(gè)與LUAD預(yù)后相關(guān)的FeRGs森林圖

表1 14個(gè)預(yù)后相關(guān)的差異表達(dá)FeRGs的差異分析

A：14個(gè)預(yù)后相關(guān)的差異表達(dá)FeRGs維恩圖；B：14個(gè)預(yù)后相關(guān)的差異表達(dá)FeRGs森林圖；C:14個(gè)預(yù)后相關(guān)的差異表達(dá)FeRGs熱圖，N為癌旁組織，T為腫瘤組織。

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

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

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

2.3 模型的獨(dú)立預(yù)后分析

單因素Cox分析顯示，風(fēng)險(xiǎn)評(píng)分與LUAD的預(yù)后有明顯的聯(lián)系(HR=3.202,95%CI=2.169～4.726,P<0.001)；多因素Cox分析顯示，風(fēng)險(xiǎn)評(píng)分為L(zhǎng)UAD總生存率的獨(dú)立預(yù)測(cè)因子(HR=2.886,95%CI=1.943～4.285,P<0.001)。見(jiàn)圖6。

2.4 功能富集分析

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

A：兩個(gè)風(fēng)險(xiǎn)組風(fēng)險(xiǎn)評(píng)分分布曲線；B：風(fēng)險(xiǎn)評(píng)分和生存狀態(tài)散點(diǎn)圖；C：兩個(gè)風(fēng)險(xiǎn)組PCA分析；D：兩個(gè)風(fēng)險(xiǎn)組t-SNE分析；E：高低風(fēng)險(xiǎn)組生存曲線對(duì)比圖；F：3年總生存率ROC曲線及AUC。

A：?jiǎn)我蛩谻ox分析初步篩選LUAD預(yù)后的預(yù)測(cè)指標(biāo)；B：多因素Cox分析進(jìn)一步篩選預(yù)測(cè)指標(biāo)，風(fēng)險(xiǎn)評(píng)分可作為L(zhǎng)UAD的獨(dú)立預(yù)測(cè)因子。

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

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

3 討 論

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

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

綜上所述，本文構(gòu)建了一個(gè)由11個(gè)FeRGs組成的新型LUAD預(yù)后模型，該模型為L(zhǎng)UAD預(yù)后的評(píng)估提供了參考。