衰老是一种自然的生物学过程，其标志是所有活生物体的生理功能逐渐退化，并与死亡的脆弱性增加有关［6］。同时衰老与多种疾病的风险增加有关，如心血管疾病、癌症、神经退行性疾病和中风［7-9］。尽管衰老相关基因在许多癌症的发展中起重要作用，但关于衰老基因在胶质瘤患者预后的应用中却少见报道。
因此本研究通过中国脑胶质瘤基因组图谱（Chinese glioma genome atlas，CGGA）建立了衰老相关多基因预后模型，并经验证确定该模型对胶质瘤患者的预后具有优秀的预测能力。
1 实验部分
1.1 实验材料
从CGGA数据库中下载TCGA RNA seq与mRNAseq_325两个数据集。其中TCGA RNA seq数据集包含702个胶质瘤患者的转录组信息与临床信息作为训练集。mRNAseq_325包含325个胶质瘤患者的转录组信息与临床信息作为验证集。得到的表达谱数据均进行log（FPKM+1）标准化，并剔除半数表达值小于1的样本与无生存时间的样本。从GSEA（gene set enrichment analysis，https：//www.gsea-msigdb.org）网站输入“aging”查找衰老相关基因集“GOBP_AGING”。该基因集包含313个基因。图像绘制采用用网页工具Sangerbox（http：//sangerbox.com/）。
1.2 方 法
1.2.1 预后相关基因筛选 使用R软件包“survival”，整合生存时间、生存状态和衰老相关基因表达数据，利用单因素回归方法评估每个基因的预后显著性。
1.2.2 风险预测模型的构建与验证 使用R软件包“glmnet”，利用拉索（least absolute shrinkage and selection operator，LASSO）回归的方法进行回归分析，进一步缩小目标范围。在得到相应目标基因后，采用多因素回归的方法，评估这些基因的预后显著性，根据相关系数建立风险预测模型，并计算各个样本的风险评分。以风险评分的中位数作为截值，把患者分为高风险组和低风险组，分析比较两组之间的生存状态，并通过受试者工作特征曲线（receiver operating characteristic curve，ROC）分析评价风险模型的准确性。将风险预测模型导入mRNAseq_325数据集计算各个样本的风险评分。以风险评分的中位数作为截值，把患者分为高风险组和低风险组，分析比较两组之间的生存状态，并通过ROC分析评价风险模型的准确性。
1.2.3 筛选影响胶质瘤患者的临床亚型 整合患者的年龄、临床、肿瘤分型、风险得分，突变情况。绘制临床亚型与风险因素关系的箱线图，从而筛选出影响胶质瘤患者生存的临床亚型。
1.2.4 富集分析 从GSEA （http：//software.broad- institute.org/gsea/index.jsp）网站获得GSEA软件（version 3.0），根据风险评分将样本分成高低两组，并从Molecular Signatures Database（http：//www.gsea-msigdb.org/gsea/downloads.jsp）下载c2.cp.kegg.v7.4.symbols.gmt子集合，用以评估相关途径和分子机制，基于基因表达谱和表型分组，设定最小基因集为5，最大基因集为5 000，1 000次重抽样（p<0.05）。
2 结果与分析
2.1 衰老相关多基因预后模型的建立
使用单因素分析得到214个与预后相关的基因，如图1所示。LASSO回归降维后发现，当λ=0.05时，模型具有较好的收缩效果。此时，共有21个基因可用于风险模型的构建，如图2所示。将入选的21个基因进行多因素回归分析，得到4个构建预后模型的基因（Casp9、Bmpr1a、Foxm1、Igfbp5）。其中，Casp9、Bmpr1a是保护因子，Foxm1、Igfbp5为风险因子，如图3所示。
基于基因表达量的预后风险模型公式为：Riskscore= Casp9×-3.577+Bmpr1a×-3.573+ Foxm1×
1.582+ Igfbp5×1.734。
通过公式，计算胶质瘤患者的风险系数，并按照中位数将患者分为高、低风险组，ROC显示：1、3、5 a的ROC下面积分别为0.91，0.95和0.95如图4（a）。反应出该模型具有优秀的准确度（ROC值>0.80）。预后热图显示：随着风险分数的增加，患者的生存率明显下降，正如预期的那样，Casp9，Bmpr1a是保护性因素，随着风险分数的增加表达呈现下调趋势，Foxm1，Igfbp5是危险因素，随着风险分数的增加表达呈现上调趋势，如图4（b）所示。生存概率（kaplan-meier，KM）生存曲线显示：低风险组患者的预后显著优于高风险组（p<0.001），如图4（c）所示。
<G:\武汉工程大学\2023\第3期\周宇能-2.tif>[-14 -12 -10 -8 -6 -4 -2
lg2（λ）][800
600
400
200
0][偏似然偏差][λ=0.05 λ=0.12]
图2 LASSO回归模型筛选预后基因
Fig. 2 LASSO regression model for screening
prognostic genes
2.2 衰老相关多基因预后模型的验证
通过风险模型公式计算mRNAseq_325中各样本的风险得分，按照中位数将患者分为高、低风险组。结果显示，ROC分析显示，1、3、5 a的ROC下面积分别为0.73，0.80和0.80，表明该风险模型在不同数据集中也具有准确性，图5（a）所示。Casp9，Bmpr1a随着风险分数的增加表达呈现下调趋势，Foxm1，Igfbp5随着风险分数的增加表达呈现上调趋势。与之前所得结果一致，图5（b）所示。低风险组患者的预后显著优于高风险组（p<0.001），图5（c）。
2.3 临床亚型风险评分的差异
将各临床亚型分组后发现，病理分级越高，风险得分也越高（p<0.001），如图6（a）所示。1p/19q染色体缺失亚型比未缺失亚型的风险得分低（p<0.001），如图6（b）所示。IDH突变状态显示，野生型较突变型具有更高的风险得分（p<0.001）如图6（c）所示。在肿瘤组织分类中，多形胶质母细胞瘤在所有胶质瘤组织分类中具有最高的风险评分（p<0.001），如图6（d）所示。性别与风险评分之间并无显著关系（p>0.05），如图6（e）所示。年长的患者风险得分显著偏高（p<0.001），如图6（f）所示。
2.4 GSEA富集结果
用GSEA比较高危组和低危组之间的京都基因和基因组百科全书基因集（Kyoto encyclopedia of gene and genomes，KEGG）基因集富集差异。结果表明，JAK_STAT_信号路径、P53_信号通路、泛素介导的蛋白水解富集在高风险组（p<0.05）。长期电位富集在低风险组。如图7所示。
<G:\武汉工程大学\2023\第3期\周宇能-7.tif>[0 5 000 10 000 15 000 20 000
数据集排序][1.2
0.8
0
-0.8
-1.2][1.0
0.5
0
-0.5
-1.0][JAK_STAT_信号路径（ES=-0.580，NP=0.006）
p53_信号路径（ES=-0.529，NP=0.039）
泛素介导的蛋白水解（ES=-0.320，NP=0.040）
长期电位（ES=0.625，NP=0.015）
][富集得分][排名表指标]
图7 GSEA信号通路
Fig. 7 GSEA signaling pathway
3 结 论
在近几十年来，尽管在胶质瘤的诊断、化疗、放疗和手术切除治疗等方面取得了进展，但患者的总生存期并未得到显著改善，平均生存期仅为12~15月，并且大多数患者在治疗后依旧会复发，因此在现阶段胶质瘤依旧被认定为不可治疗的疾病［10-12］。但胶质瘤难以治愈的原因是多元的，包括血脑屏障对药物治疗的阻碍，胶质瘤异质性导致的免疫逃逸等［10］。由于当前关于胶质瘤的高度异质性导致其不可治愈已成为当前研究者的共识。因此研究者对胶质瘤的异质性进行了不同层面的表征，以希望能够全面了解胶质瘤的特性。
研究表明，胶质瘤可以通过多种通路对自身存活、治疗反应以及患者预后产生影响。如胶质瘤细胞受到高度缺氧环境的刺激后，可以通过缺氧通路对低氧分压环境做出反应，增加血管内皮生长因子（vascular endothelial growth factor，VEGF）等相关基因的表达，促进血管生成，为肿瘤细胞提供所需的氧气和营养物质，从而保证肿瘤细胞的存活，最终对患者的生存状况产生了重要影响［13］。各种代谢通路的改变则使得胶质瘤细胞可以利用多种物质产生能量，以适应各种不利的生存环境［14］。而在关于人类衰老过程对胶质瘤的影响中，临床研究表明，与其他癌症一样，老年人占胶质瘤患者的多数。同时，老年患病群体表现出更强的肿瘤侵袭性，更大的治疗耐药性，以及更差的治疗有效性［15-16］。基于分子水平的研究表明，衰老通路上游靶点Nlrp3基因可以控制细胞衰老进程，这一通路可以通过外界辐射激活。细胞衰老的进程可以通过释放衰老相关的分泌表型促进肿瘤的进展和复发。当Nlrp3基因受到抑制时，胶质瘤细胞增长速度降低，且对放射治疗更加敏感［17］。因此，衰老相关基因对于胶质瘤的治疗与预后具有深远影响。
本文在通过生物信息学全面探讨了衰老相关基因与胶质瘤患者预后之间的关系后建立了包含Casp9、Bmpr1a、Foxm1、Igfbp5这些衰老相关基因的预测模型。其中Casp9是半胱氨酸天冬氨酸蛋白酶9的编码基因，是内源性凋亡的起始caspase［18］。许多研究证明Casp9在p53依赖性凋亡中起决定性作用，并能够通过蛋白质水解作用切割并激活半胱天冬酶，如半胱氨酸蛋白酶-3（Casp3）和半胱氨酸蛋白酶-7（Casp7）等，诱导癌细胞的凋亡［19- 20］。当其发生突变或激活失败时，往往会导致癌细胞具有更高的死亡阈值［21］。因此Casp9属于保护因子，与所预测的结果一致。骨形态发生蛋白受体1A型（Bmpr1a），可以与细胞表面的两种丝氨酸/苏氨酸激酶受体（I型和II型）结合并形成异聚复合物，这两种受体都是信号转导所必需的［22］。因此，Bmpr1a与多种肿瘤的生长、分化或凋亡相关。但Bmpr1a是否是癌基因在不同的研究中表现出不同的结论［23］。有研究使用MiR-656直接靶向Bmpr1a，在BMP-2存在或不存在的情况下，Bmpr1a反过来抑制胶质瘤细胞增殖、神经球形成、迁移和侵袭［24］。因此，关于Bmpr1a在胶质瘤中的具体作用与机制尚需进一步探究。Foxm1基因编码的蛋白质是参与细胞增殖的转录激活剂，能够调节细胞周期基因的表达，如Cdc25B和有丝分裂细胞周期蛋白［25］。同时，也是人类实体癌中最常见的过度表达蛋白之一［26］。许多研究表明，Foxm1的过度表达与肿瘤细胞的耐药性、癌症分级以及患者的不良预后息息相关［27-28］。在胶质瘤的相关研究中发现，Foxm1能够通过参与DNA修复的复制因子C亚单位5（Rfc5）的反式激活而导致GBM的化疗耐药性［29］。当Foxm1的表达受到抑制时，与DNA修复相关的基因（包括Mre11和Rad51）的表达水平显著降低，从而使GBM细胞对放疗更加敏感［30］。因此，正如分析结果，Foxm1是影响胶质瘤患者生存的高风险因素。许多研究表明，Igfbp5与多种癌症的发展有关。但Igfbp5在不同癌症类型或相同组织类型癌症中的表达水平和功能存在差异。一方面，有研究表明，Igfbp5可以通过促进Pknox2介导的p53上调抑制胃癌细胞的生长与迁移［31］。另一方面，其他研究表明Igfbp5可以作为癌基因促进肿瘤转移［32］。在近期关于胶质瘤的研究中发现，Igfbp5能够通过调节EMT促进GBM细胞侵袭，但同时通过抑制AKT的磷酸化抑制GBM细胞增殖［33］。因此，Igfbp5对GBM的影响具有复杂性，更为细致的作用机制有待后续研究。
综上所述，本研究基于CGGA数据库鉴定了4个衰老相关基因作为反映胶质瘤预后的潜在标志物，并通过不同数据集对所构建的模型进行了验证。结果表明，该模型能够准确地预测胶质瘤患者的预后，且准确性高于绝大多数已发表的模型，如自噬相关基因预后模型［34］，缺氧相关基因预后模型［35］等。但仍需要后续临床研究验证模型的诊断能力。