神经胶质瘤(Glioma)是最常见和最致命的原发性脑肿瘤,约占中枢神经系统恶性肿瘤的80%,具有明显的临床异质性,其中胶质母细胞瘤最具侵袭性和致命性。肿瘤的快速生长和异质性是其侵袭性进展的重要原因,表现为神经功能障碍和认知能力下降等症状。传统治疗包括术后放疗和替莫唑胺联合治疗,然后以替莫唑胺辅助治疗
[1]。然而,肿瘤的侵袭性及其在脑组织中的深层位置使完全切除变得非常困难,即使手术成功,残留的肿瘤细胞仍然可能导致复发。此外,Glioma经常对传统放疗和化疗的药物产生耐药。血脑屏障的存在进一步阻碍了治疗剂对肿瘤组织的输送
[2]。因此,需要迫切寻求新的诊疗方式。铁死亡是一种新型细胞程序性死亡,它依赖于铁,细胞内铁、ROS积聚、谷胱甘肽耗竭以及脂质过氧化的累积都是铁死亡的关键因素
[3]。近期研究表明,在小鼠胶质瘤模型中,铁基纳米材料可以产生ROS诱导脂质过氧化引起铁死亡
[4‑6]。这表明胶质瘤与铁死亡存在密切联系。
脂质过氧化和氧化脂质也可诱导MAP1LC3周转和自噬体形成。越来越多的证据表明,过量的自噬和溶酶体活性可以通过铁积累和脂质过氧化来促进铁死亡
[7,8]。铁死亡依赖性诱导自噬也被认为是抗肿瘤的有效策略。自噬是通过溶酶体去除和降解细胞内成分的过程,如未使用的蛋白质和受损的细胞器
[9]。自噬货物被隔离在吞噬载体中,吞噬载体形成自噬体,然后与溶酶体融合形成自噬‑溶酶体,它调节了细胞稳态。这一过程在哺乳动物模型中得到了广泛的研究。尽管诱导自噬的铁死亡是一种很有潜力的治疗方式,但它们在胶质瘤中的作用机制尚未完全阐明。因此,本研究目的是确定与Glioma的发展相关的自噬和铁死亡基因。目前,尚缺乏关于Glioma的自噬‑铁死亡基因机制的生物信息学方面的研究。本研究使用生物信息学分析了Glioma的公众数据集,在基因表达综合数据库(Gene Expression Omnibus,GEO)中下载数据集,然后将分析其与自噬‑铁死亡(autophagy ‑ferroptosis,AF)相关的基因。旨在找到具有差异表达的特征基因,为胶质瘤疾病中的自噬和铁死亡提供新的思路。
1 资料和方法
1.1 数据收集
本研究中从4个关于Glioma患者及其正常组织样本的高通量测序数据集:GSE137902、GSE21354、GSE31262和GSE90598从GEO基因表达数据库(
https://www.ncbi.nlm.nih.gov/geo/)中选取数据集
[10,11],这些数据集中的基因表达谱包含48个Glioma患者和22个健康对照组。
1.2 Glioma的差异表达基因的分析
收集的数据集通过使用R软件4.4.1“limma”包进行分位数归一化和log2转换
[10]。在使用了R软件中“sva”包合并这些数据集后,使用“Combat”算法来消除批处理的影响
[11]。主成分分析法(principal components analysis,PCA)可以消除数据集之间的批次效应。接下来,使用R软件的“limma”包对Glioma患者和正常组织之间的差异表达基因(differentially expressed genes,DEGs)进行标准化和差异分析。DEGs的统计学显著性标准是“调整后
P<0.05和|log2差异倍数(FC)|>0.585差异倍数(foldchange,FC)”。此外,使用“ggplot2”包、“pheatmap2.0”包分别构建DEGs的火山图和热图
[10‑12]。
1.3 免疫浸润分析
使用“CIBERSORT”算法(
https://ciber sortx.Stanf Ord.edu/)评估每个样品的特异性免疫细胞富集的相对比例,排除所有样品中具有零值的免疫细胞类型
[11]。然后,通过使用R软件中的“vioplot”包用箱线图显示了22个免疫细胞亚型之间的差异,免疫细胞类型的浸润水平在Glioma和正常组织中差异
P<0.05,有统计学意义。通过使用“corrplot”包展示了 22 种免疫细胞亚型之间的相关性
[13]。
1.4 Glioma和自噬‑铁死亡相关差异基因(DEAFGs)的筛选
先从Genecards数据库(
https://www.genecards.org/cgibin/carddisp.pl)中铁死亡基因和自噬基因进行筛选,以“相关系数(Relevance score)>1”作为筛选的阈值,收集到907个铁死亡基因和3 926个自噬基因
[14]。然后,使用 “VennDiagram”软件包将这些基因组合成AF相关基因。然后,将这些基因与Glioma的DEGs进行交集,以确定与Glioma相关DEAFGs,最后将结果进行可视化
[15]。
1.5 KEGG/GO富集分析
在京都基因和基因组百科全书(KEGG)和基因本体(GO)的功能富集分析中,使用R软件的“clusterProfiler”包对Glioma相关的DEAFGs进行了功能富集分析。为了研究差异表达基因的功能作用,并注释其生物学过程和机制通路,统计学显著性为
P<0.05
[16]。
1.6 机器学习方法鉴定特征基因
使用“glmnet”包的最小绝对收缩和选择算子(LASSO)和“randomForest”包的随机森林分析(RF)算法,寻找与Glioma相关的自噬‑铁死亡特征基因
[17]。通过对基因的重要性进行排序和交集,找到了与Glioma相关的自噬‑铁死亡特征基因。差异性分析检验用于确定筛选的特征基因的表达差异。此外,使用工作特征曲线(ROC)来评估这些基因的性能。首先,使用“pROC”包确定临界值,然后计算出曲线下面积(AUC),评估特征基因的临床诊断意义及其样本准确性(AUC大于0.7表示具有较高的准确性),最后对上述结果进行可视化
[18]。
1.7 免疫组织化学(IHC)分析
使用小鼠胶质瘤GL261细胞系建立C57BL/6小鼠原位胶质瘤模型。使用IHC分析检测正常脑组织和小鼠Glioma组织中特征基因蛋白的表达水平。先将提取的组织进行固定、包埋、切片,然后将它们放入柠檬酸缓冲液中进行抗原修复。然后在4 ℃下与一抗稀释液一起孵育过夜。然后,它们与二抗在37 ℃孵育30 min。之后用苏木精进行核复染,最后封片后观察。
1.8 miRNAs‑genes调控网络构建
根据两种机器学习方法获取的6个特征基因,利用3个在线数据库TargetScan(
http://www.targetscan.org/)、miRDB(
http://www.mirdb.org/)和miRanda(
http://www.microrna.org/)去预测与这些特征基因相关的miRNA靶基因,并构建miRNAs‑genes调控网络,同时使用这3个数据库预测时,相关性的可靠性更高
[19]。并使用 Cytoscape软件进行可视化。
2 结果
2.1 Glioma相关基因的分析结果
将4个Glioma数据集(GSE137902、GSE21354、 GSE31262和 GSE90598)合并,排除了批次效应。结果显示,批次矫正之前,不同4个数据集是明显分散的(
图1A),但在使用PCA降维去除批次效应后,4个数据集中的数据分布相当均匀且在一个范围内的(
图1B),说明消除了批次效应对后续分析的影响。接下来,通过比较正常对照组和Glioma组,设置
P<0.05,|log2FC|≥0.585,得到1 376个差异基因,其中763个下调基因和613个上调基因(
图1C、D)。
2.2 免疫侵润分析
进一步了解Glioma的免疫细胞浸润的相关性,使用CIBERSORT算法评估了免疫细胞浸润的比例,排除了预测免疫细胞浸润矩阵中所有样品中具有零值的细胞类型。结果显示了22种类型免疫细胞的浸润图谱(
图2A)。差异分析显示,Glioma组与对照组相比表现出更高水平的CD4记忆激活的T细胞、CD4记忆终止的T细胞、单核细胞和嗜酸性粒细胞的浸润水平(
P<0.05)(
图2B)。相反,在CD8
+T细胞和NK细胞激活中表现出较低的水平(
P<0.05)(
图2B)。此外,22个免疫细胞之间的关系(
图2C)中,CD4幼稚T细胞和嗜酸性粒细胞都与激活的NK细胞(
r=0.77)、单核细胞(
r=0.77)呈显著正相关,但也都与调节性T 细胞(Tregs)(
r=-0.77)呈负相关等。
2.3 Glioma和DEAFGs的分析
将来自Genecards数据库的907个铁死亡基因和3 926个自噬基因合并成AF‑related genes数据集。然后将这些数据集与来自GEO的DEGs进行交集,以确定关于Glioma的AF‑related genes的DEAFGs。使用Venn图(
图3A)识别出12个与关于胶质瘤的DEAFGs。接下来,对交集到的DEAFGs进行差异分析。分别发现了10个上调基因(
CD44、TNFAIP3、YAP1、NEDD4、VIM、MYC、BIRC5、TP53、AIM2、HIF1A)和2个下调基因(
NEDD4L、FBXW7),这些基因在
图3B的热图中显示。
2.4 DEAFGs功能富集途径及分析
DEAFGs进行GO功能和KEGG信号通路富集分析。生物学过程(BP)、 细胞组分(CC)和分子功能(MF)是GO富集分析的3个部分(
图4A~C)。GO分析发现这些基因显著富集在BP类的凋亡信号通路、线粒体和巨自噬细胞的调节、有丝分裂的调节和对辐射的反应。富集的主要MF包括泛素蛋白结合、DNA结合转录因子结合及DNA结合转录阻遏物活性。主要富集CC包括RNA聚合酶Ⅱ转录调节复合物、生殖细胞核。
KEGG是一个数据库资源,有助于进行信号通路分析。通过对特定基因进行KEGG信号通路富集分析,发现这些基因在Epstein-Barr病毒感染、癌症中的中枢碳代谢、癌症蛋白聚糖、癌症微小RNA、结直肠癌、泛素介导的蛋白水解、海马信号通路、甲状腺癌、乙型肝炎、膀胱癌、Kaposi sarcoma相关疱疹病毒感染等通路显著富集(
图4D、E)。
2.5 特征基因的筛选
筛选Glioma相关自噬‑铁死亡的特征基因,在LASSO分析中,共选择了6个特征性DEAFGs(
图5A、B)。另外,随机森林分析发现了12个具有相对重要性>1的特征DEAFGs(
图5C、D),这些特征基因见
表1。然后使用Venn图将两种算法交集,最终确定了6个特征基因(
图5E)。这些基因包括
HIF1A、TNFAIP3、YAP1、NEDD4L、BIRC5、AIM2。
2.6 特征基因鉴定分析
将选择的6个特征基因进行差异性分析,条件是
P < 0.05,|log2FC|≥1,这些特征基因中有5个上调基因(
HIF1A、TNFAIP3、YAP1、BIRC5、AIM2)和1个下调基因(
NEDD4L),见
图6A。此外,进行了ROC曲线分析(
图6B),以研究这6个特征基因的潜在准确性价值,ROC曲线的水平和垂直坐标分别表示灵敏度和特异性,结果显示这六个特征基因
HIF1A、TNFAIP3、YAP1、NEDD4L、BIRC5、AIM2的AUC分别为0.867、0.822、0.809、0.829、0.747和0.723,AUC均大于0.7。
2.7 IHC分析验证
为了验证上述结果的可靠性,本研究建立了GL261细胞系的小鼠原位胶质瘤模型,将随机选取的2个特征基因(
YAP1和
BIRC5)进行 IHC 分析,验证Glioma与正常脑组织之间的特征基因蛋白表达之间的相关性。IHC染色定量评分结果(
图7)显示,Glioma细胞中
YAP1和
BIRC5表达水平显著增高,与正常脑组织相比呈显著差异(
P<0.05)。
2.8 miRNA和特征基因的相互作用
使用网络分析,进一步了解特征基因和miRNAs之间的相互作用。根据TargetScan、miRDB和miRanda 3个在线数据库,对6个特征基因(
HIF1A、TNFAIP3、YAP1、BIRC5、AIM2和NEDD4L)的预测,然后将预测的miRNA进行交集,并构建miRNA‑genes相互调控网络(
图8A、B)。结果显示,这个调控网络中存在共同的枢纽基因,共包含49个节点,其中6个特征基因与43个miRNA(如
hsa‑miR‑27b‑3p、hsa‑miR‑141‑3p、hsa‑miR‑513a‑5p等)存在相互作用。
3 讨论
近些年,铁死亡被视为肿瘤学中一个有希望的前沿领域,它不仅可以抑制肿瘤生长,还可以增强免疫治疗反应并克服癌症治疗的耐药性,有研究表明自噬引起的铁死亡在Glioma发展中发挥着重要作用
[20,21],但其机制尚不清楚。通过对GEO数据库数据集的综合分析、生物学功能预测、免疫细胞浸润分析、特征基因筛选以及miRNA‑genes调控网络的构建,本研究旨在探讨自噬和铁死亡在Glioma中的潜在作用机制,并为未来抗肿瘤策略开发潜在靶点提供新的见解。
本研究从GEO数据库中提取了Glioma数据集,对它们相关的DEGs进行分析,筛选出了613个上调和763个下调的差异表达基因。通过分析胶质瘤中免疫细胞浸润的动态变化,发现在Glioma中,CD4记忆终止的T细胞、CD4记忆激活的T细胞、单核细胞和嗜酸性粒细胞的比例都较高。可能是因为胶质母细胞瘤是所有癌症组织学中T细胞亚群最丰富的肿瘤之一,CD4
+ T细胞通过增强抗肿瘤反应能激活中间免疫亚群,然后刺激表达相容性复合物(MHCII)的亚型(如单核细胞和小胶质细胞)进而通过分泌刺激性细胞因子(如IFN‑γ)来启动免疫反应,招募嗜酸性粒细胞
[22]。此外,发现Glioma患者的 CD8
+T细胞和激活的NK细胞比例下调,可能是由于程序性细胞死亡导致CD8
+ T细胞和NK细胞丢失,或由于细胞毒性淋巴细胞能通过杀死 T 细胞、NK 细胞和抗原呈递细胞来降低免疫激活
[23]。表明免疫细胞在Glioma中发挥了重要免疫功能。
本研究收集并差异性分析胶质瘤的DEGs与自噬‑铁死亡相关基因。发现12个与Glioma发展、进展和预后有关的DEAFGs,即
CD44、TNFAIP3、YAP1、NEDD4L、VIM、MYC、BIRC5、TP53、AIM2、NEDD4、HIF1A和
FBXW7。其中,CD44是Glioma细胞膜上高表达的一种多功能跨膜糖蛋白,是透明质酸(HA)的受体,当CD44和HA结合时,CD44的构象会发生变化,这会促进细胞信号通路来调控肿瘤细胞的增殖和迁移
[24]。波形蛋白在Glioma中的表达受生存素(BIRC5)的调节,BIRC5是高级别胶质瘤不良预后因素
[25]。肿瘤抑癌基因
TP53的突变会与癌蛋白相互作用,诱导星形细胞瘤发生
[26]。此外,NEDD4 是 E3 泛素蛋白连接酶,负责泛素‑蛋白酶体系统中的底物特异性,它的失调会促进肿瘤的进展
[27]。FBXW7作为F‑box 蛋白家族的一员,是一种重要的肿瘤抑制因子。由于p53突变,FBXW7表达会降低,这会导致致癌蛋白c‑MYC的积累,并促进神经胶质瘤的发展
[28,29]。本研究对这些基因进行了功能富集分析,以预测它们潜在的生物学功能和通路,发现这些基因大部分参与传染性病毒感染通路、癌症相关的通路、泛素介导的蛋白水解等信号通路。有研究报道在Glioma中检测到几种病毒
[30],包括人巨细胞病毒、人瘤病毒、EB 病毒、卡波西肉瘤相关疱疹病毒、乙型和丙型肝炎病毒,原因是传染性病原体与Glioma的形成有关联。此外,这些基因的表达在多种癌症通路中都出现。例如,Yes相关蛋白1在食道癌、肝癌以及膀胱癌中的表达促进了疾病进展
[31‑33]。
为了进一步探究AF与Glioma之间的特征基因,使用两种机器学习—Lasso回归和RF分析发现了6个与Glioma相关的特征基因。使用差异性分析揭示了5个上调(
HIF1A、TNFAIP3、YAP1、BIRC5、AIM2)和1个下调(
NEDD4L)特征基因。ROC曲线表明这些特征基因识别Glioma的准确性较高。其中缺氧诱导因子(HIF1A)是转录因子,可以直接或间接激活转录因子引起的肿瘤生物学及其微环境的改变,从而增强肿瘤的侵袭能力
[34]。肿瘤坏死因子α诱导蛋白3(TNFAIP3),也称为A20,是一种泛素编辑酶,在Glioma组织中高表达,它具有泛素连接酶和去泛素酶活性,可以通过RIP1中去除K63 泛素链并将K48泛素链添加到RIP1中进行降解从而抑制神经胶质瘤细胞的凋亡
[35]。YAP1是Hippo 信号通路的关键效应子,Hippo/YAP 通路是一个激酶级联反应,磷酸化两种转录共激活因子YAP/TAZ。YAP通过抑制GSK3β然后激活 β‑catenin 促进Glioma的发展,还通过上调HMGB1 介导自噬促进胶质瘤进展
[36]。BIRC5是细胞凋亡蛋白抑制剂(IAP)家族中新发现的成员,在恶性癌症中高表达,但在正常组织中低表达,IAP 蛋白家族至少存在一种凋亡重复结构域(BIR),其中的经典成员 XIAP可以和BIRC5相互作用,然后直接结合并抑制半胱天冬酶活性,这有助于阻止胶质瘤细胞凋亡
[37]。黑色素瘤缺乏因子2(AIM2)是炎症小体的调控蛋白之一。它介导炎细胞因子IL‑1β和IL‑18的成熟和释放,IL‑1β 和 IL‑18 还有助于形成瘤内旁分泌环,并支持血管生成和肿瘤增殖
[38]。神经前体细胞表达发育性下调的4样(NEDD4L)是 NEDD4 家族的成员,它的过表达会介导肿瘤癌基因鞘氨醇激酶 2的泛素化,从而抑制恶性神经胶质瘤的发展
[39]。此外,通过IHC 实验分析证实了这些特征基因的可靠性。因此,这些特征基因很可能成为Glioma的诊断、治疗及预后的靶点。
miRNA作为内源性非编码RNA小分子,可以靶向mRNA抑制翻译或触发mRNA降解
[40]。miRNA与脑部肿瘤(包括高、低级别神经胶质瘤)的发病机制密切相关
[41]。构建的miRNA和特征基因调控网络显示
miR‑27b‑3p、miR‑141‑3p、miR‑485–5p、miR‑513a‑5p和
miR‑497‑5p等miRNA 参与特征基因的调控。其中miR‑27b‑3p、miR‑141‑3p通过调节YAP1控制胶质瘤细胞的增殖,迁移和凋亡发挥抗癌功能
[42,43]。
NEDD4L被验证是miR‑513a‑5p的直接靶基因,由于miR‑513a‑5p 的过表达, NEDD4L 介导的胶质瘤细胞对替莫唑胺的敏感受到很大影响
[44]。此外,在常氧下过表达miR‑485–5p会抑制胶质瘤细胞的活力,而缺氧时HIF1A会与启动子区域结合,从而阻止miR‑485‑5p 转录
[34]。这些研究结果支持了本研究生物信息学预测。因此,挖掘与这些特征基因相关的miRNA可以为Glioma的诊断和治疗提供新的靶点,但仍需要继续探索其具体机制。
尽管本研究方法新颖,但它仍然存在一些限制。本研究结果基于数据库的Glioma数据集进行分析,由于样本量有限,未进行具体分型。总之,这项研究通过生物信息学分析发现了与Glioma相关的AF的潜在生物标志物和靶点,为Glioma发病机制和临床治疗的研究提供了新的思路。
作者贡献度说明:
杨小红:数据筛选、数据分析、文章绘图及撰写论文;陈旭永:数据分析、参与文献查找;应卓鹏、方贾栩:给予相关修改意见;陈建强、战跃福、关莹:指导文章结构设计、论文修改及内容校审。
所有作者声明不存在利益冲突关系。
国家自然科学基金资助项目(82260343)
海南省重点研发计划项目(ZDYF2023SHFZ142)
海南省卫生健康科技创新面上项目(WSJK2024MS136)
海南省自然科学基金面上项目(821RC692)