基于GEO数据库构建神经母细胞瘤缺氧相关的预后模型

Construction of Hypoxia-Related Prognostic Signature in Neuroblastoma Based on GEO Database

摘要: 目的：基于基因表达综合(Gene Expression Omnibus, GEO)数据库分析缺氧相关基因与神经母细胞瘤(NB)临床特征及免疫环境的关系，构建缺氧相关预后模型。方法：从GEO数据库下载数据集GSE62564，从MSigDB数据库下载缺氧相关基因集。通过共识聚类算法将NB患者分为高低缺氧组，通过Kaplan-Meier曲线、卡方检验、“ESTIMATE”算法等比较两组间临床、免疫环境差别。通过LASSO回归建立缺氧相关预后模型，使用受试者工作特征(ROC)曲线下面积(AUC)评价模型的预测性能，并对模型进行独立预后分析。结果：根据缺氧相关基因将NB样本分为两组，两组患者的生存、年龄、肿瘤分期及MYCN状态均存在显著差异。两组样本的免疫浸润情况存在显著差异，B细胞、树突状细胞、NK细胞及T细胞等多种免疫细胞含量在两组间存在差异。通过LASSO回归得到了6个缺氧相关基因(CSRP2, DTNA, SAP30, NCAN, WSB1, NAGK)构成的预后模型。通过模型的风险评分中位数将患者分为高、低风险组，Kaplan-Meier曲线显示两组间预后存在显著差异。模型在1，3，5年时的AUC值分别为0.882，0.916及0.914，说明其具有良好的预测价值。多因素COX回归分析表明风险评分可作为NB的独立预后因素。结论：缺氧影响NB的免疫环境及患者预后，由缺氧相关基因构成的预后模型能较好评估NB患者预后，并为寻找新的治疗靶点提供帮助。

Abstract: Aims: Analyzing relationship between hypoxia-related genes and clinical features and immune en-vironment of neuroblastoma (NB) based on GEO (Gene Expression Omnibus) database to construct hypoxia-related prognostic signature. Methods: GSE62564 was downloaded from GEO database and hypoxia-related gene set was obtained from MSigDB database. Patients in GSE62564 were divided into high/low hypoxia subgroups via Consensus Clustering analysis. Differences in immune and clinical features between two groups were identified by Kaplan-Meier curve, Chi- square test and ESTIMATE algorithm. Hypoxia-related prognostic signature was constructed by Lasso-Cox regres-sion. AUC value of the signature was calculated for the evaluation of the prognostic model. Results: Clinical features including prognosis, age, tumor stage and MYCN status were significantly different in two subgroups. Immune environment and distribution of immune cell including B cells, T cells, NK cells and dendritic cells were also different between groups. A six-gene hypoxia-related prog-nostic signature (CSRP2, DTNA, SAP30, NCAN, WSB1, NAGK) was constructed using Lasso-Cox re-gression, and patients were divided into high-/low-risk group by median of risk score. Kaplan-Meier curve revealed the significant difference between prognosis of the two groups. The AUCs for the 1-, 3-, 5-year OS predictions for the signature were 0.882, 0.916 and 0.914 which revealed the prog-nostic value of the signature. Risk score was identified as an independent prognostic factor by the multivariate Cox regression analysis. Conclusion: Hypoxia has an effect on immune environment and prognosis of NB patients. A hypoxia-related prognostic signature was constructed, which will contribute to predicting prognosis of NB and finding new therapy target.

1. 引言

神经母细胞瘤(neuroblastoma, NB)是一种交感神经系统颅外实体瘤，源自神经嵴细胞。作为最常见的儿童实体肿瘤，90%的NB发生于5岁之前 [1] 。因NB死亡的人数约占儿童癌症相关死亡人数的10% [2] 。NB的临床表现具有很强的异质性。部分NB患者可自发缓解，而另一部分患者即使经过手术治疗，肿瘤仍会发生进展 [3] 。根据MYCN状态、年龄、分期及分级，可将NB患者分为高风险及低风险级别，例如Susan等在2009年的研究中描述的国际神经母细胞瘤风险组(INRG) [4] 。尽管手术、放化疗和免疫治疗等治疗手段在不断进步，高风险患者的五年生存率仍低于60% [5] [6] 。因此，寻找新的生物标志物和治疗靶点有重要的临床意义。缺氧是实体肿瘤的重要特征。肿瘤细胞快速增殖及肿瘤血管结构、功能异常都会导致肿瘤组织区域缺氧的发生 [7] 。在包括NB在内的多种肿瘤中，缺氧状态都被证实可以影响肿瘤的进展、侵袭和转移 [8] [9] [10] 。此外，肿瘤的缺氧状态与肿瘤免疫微环境也可发生相互作用 [11] 。本研究通过生物信息学分析缺氧相关基因与NB临床特征及免疫环境的关系，并构建缺氧相关预后模型，为NB的早期诊断与治疗提供新思路。

2. 资料与方法

2.1. 数据收集与处理

从基因表达综合(Gene Expression Omnibus, GEO)数据库下载数据集GSE62564。从分子特征数据库(Molecular Signatures Database, MSigDB)下载缺氧相关数据集，数据集内包含200个缺氧相关基因。

2.2. 缺氧相关基因的共识聚类分析

在R软件版本(4.3.0)中使用“ConsensuClusterPlus”包对GSE62564中的样本进行聚类分析。利用缺氧相关基因的表达量将样本分为两组。研究中使用了1000次迭代的共识聚类算法，在每次迭代中采样80%的数据。使用“survival”包对两组患者的预后情况进行分析并绘制Kaplan-Meier曲线。

2.3. 缺氧状态与肿瘤免疫环境的相关性分析

使用R软件中的“CIBERSORT”通过反卷积算法，估算肿瘤样本中22种免疫细胞的组成比例，并比较不同分组间免疫细胞比例差异。通过“Estimate”算法计算每个样本的免疫评分及间质评分，比较高、低缺氧分组间分数的差异。

2.4. 缺氧相关基因预后风险模型的构建与评价

下载GSE62564数据集中患者的生存信息，使用R中的“survival”包探究缺氧相关基因对NB患者预后的影响。P小于0.05的基因作为生存相关基因。使用“glmnet”包，基于LASSO回归算法建立缺氧相关预后模型，并通过预后模型为每个患者计算风险评分。接着，患者根据风险评分的中位数被分为高风险和低风险组。通过Kaplan-Meier曲线比较两组间预后情况。通过R中的“ROCR”包构建模型的受试者工作(receiver operating characteristic, ROC)曲线并计算曲线下面积(area under curve, AUC)值。最后通过单因素、多因素分析验证风险评分是否为独立预后因素。

2.5. 列线图模型的构建

使用R包“rms”，综合风险评分及年龄、临床分期、MYCN状态及肿瘤进展情况，构建列线图预后模型。使用“ROCR”包绘制该模型用于预测NB患者1，3，5年生存情况的ROC曲线。使用“rmda”包绘制临床决策曲线(DCA)。

3. 结果

3.1. 通过缺氧相关基因对NB进行分类

下载GSE62564数据集中498位患者的mRNA表达矩阵及临床信息。将从MSigDB数据库中下载的200个缺氧相关基因与GSE62564中的基因名重叠，共有184个缺氧相关基因。通过共识聚类算法，我们将498例NB样本根据缺氧相关基因表达量进行分组。当聚类系数k取2至6时，k = 2时取得较好结果，即组内相关性高，组间相关性低(图1(a))。据此，我们将498例样本分为聚类1 (cluster1)和聚类2 (cluster2)，分别包含342例和156例样本。缺氧相关基因在两组间的表达情况如图1(b)所示。Kaplan-Meier生存曲线显示cluster1中患者的预后显著好于cluster2 (P < 0.001)，表明缺氧相关基因与患者预后显著相关(图1(c))。

3.2. 不同缺氧状态下的免疫环境探究

肿瘤缺氧状态与肿瘤免疫微环境关系密切。为探究不同缺氧状态下肿瘤微环境的改变，我们使用R中的“ESTIMATE”算法评估不同样本免疫浸润的状况。如图2(a)所示，cluster1中样本的基质评分及免疫评分均高于cluster2，表明cluster1中样本具有更多的免疫细胞及基质细胞浸润。为进一步探究各样本中免疫细胞浸润的具体情况，我们使用“CIBERSORT”算法评估22种免疫细胞在样本中的含量。如图2(b)所示，B细胞、树突状细胞、NK细胞及T细胞等多种免疫细胞含量在两组间存在差异。

3.3. 探究缺氧状态与临床特征的关系

Cluster1及cluster2中患者的临床信息如表1所示，包括年龄、性别、临床分期及MYCN状态。使用卡方检验比较两组患者的临床特征。结果显示，两组患者的年龄、临床分期、MYCN状态存在显著差异(P < 0.05)。两组患者的性别无显著差异。

Figure 1. Dividing NB patients into high/low hypoxia groups by hypoxia-related genes

图1. 通过缺氧相关基因将NB患者分为高/低缺氧组

Figure 2. Differences in immune environment between cluster 1 and cluster 2

图2. 组1和组2间免疫环境的差别

Table 1. Relationship between hypoxia and clinical features of NB patients

表1. 缺氧与NB患者临床特征的相关性

3.4. 构建缺氧相关预后模型

为进一步完善对NB患者预后的评估，并筛选可能的治疗靶点，我们构建了缺氧相关预后模型。首先，通过单因素COX分析方法，结合498位患者的预后信息，筛选对NB患者生存具有显著影响的缺氧相关基因，共106个(P < 0.05)。将这106个生存相关基因纳入分析，通过LASSO回归算法构建预后模型(图3(a))。构建的预后模型由六个缺氧相关基因(CSRP2, DTNA, SAP30, NCAN, WSB1, NAGK)及相应系数组成。通过预后模型为每个患者计算风险评分，计算公式为：风险评分 = (CSRP2表达量 × 0.08370066) + (SAP30表达量 × 0.28005528) + (NCAN表达量 × 0.01927312) − (WSB1表达量 × 0.42357840) − (NAGK表达量 × 0.10539461) − (DTNA表达量 × 0.06373091)。

Figure 3. Construction of hypoxia-related prognostic signature by Lasso-Cox regression

图3. 通过LASSO回归构建缺氧相关预后模型

3.5. 对预测模型预测效果的评估

根据风险评分的中位数，将498位患者分为高风险组及低风险组。使用R中的“pheatmap”包绘制6个基因在各样本中的表达情况(图3(b))。通过“survival”包绘制Kaplan-Meier曲线(图4(a))，发现低风险组患者预后显著好于高风险组(P < 0.001)。箱线图(图4(b))说明生存患者的风险评分显著高于死亡患者(P < 0.001)。通过ROC曲线计算预测模型在预测患者1、3、5年时生存情况的AUC值，分别为0.882，0.916及0.914，说明该模型对患者预后具有较好的预测效果(图4(c))。

Figure 4. Prediction value of hypoxia-related prognostic signature

图4. 缺氧相关预后模型的预测价值

3.6. 临床指标及风险评分与预后的关系

收集NB患者的临床信息，包括性别、年龄、肿瘤分期、MYCN状态。使用多因素COX回归分析各临床指标及风险评分与预后的关系。如表2所示，年龄、分期、MYCN状态及风险评分与预后显著相关。多因素COX分析表明年龄、分期及风险评分可作为患者预后的独立预测指标。

3.7. 构建列线图模型

结合风险评分及与预后相关的年龄、肿瘤分期，构建用于预测NB患者预后的列线图模型(图5(a))。根据模型计算每位患者的评分，并通过ROC曲线计算模型在1、3、5年的AUC值，分别为0.895，0.921及0.919。比较列线图与肿瘤分期的AUC值，可见列线图的预测效率高于肿瘤分期(图5(b))。通过临床决策曲线比较列线图模型与肿瘤分期诊断的准确性(图5(c))。

Table 2. The multivariate Cox regression analysis of prognostic factors in NB patients

表2. 影响NB患者预后的多因素分析

(a) (b) (c)

Figure 5. Construction of prognostic nomogram

图5. 构建列线图模型

4. 讨论

神经母细胞瘤(NB)是一种在交感神经系统发育过程中，起源于神经嵴的实体肿瘤 [12] 。尽管NB的发病率并不高，但其仍是儿童最常见的颅外恶性实体肿瘤，占据儿童癌症死亡人数的10% [13] 。NB具有很强的异质性，其中高风险亚型对传统治疗反应差，患者预后较差 [14] 。肿瘤细胞增殖具有快速且不受控制的特点，限制了肿瘤区域氧气的使用。因此，血液供应不足及缺氧是几乎所有实体肿瘤的微环境特征 [15] 。虽然在肿瘤中有新生血管形成，但持续缺氧会引起肿瘤血管分布不规则，毛细血管之间距离增加，超过氧气的扩散能力 [16] [17] 。肿瘤细胞在缺氧环境中侵袭性及耐药性可能增强。缺氧诱导基因表达变化，对各种细胞的生理功能有重要影响，并最终影响患者预后 [18] 。例如在肝细胞癌中，缺氧可通过影响Treg细胞及树突状细胞，介导肝细胞癌免疫抑制 [19] 。在乳腺癌中，缺氧相关基因ADAM12的表达变化可介导肿瘤的侵袭与转移 [20] 。

作为明确的肿瘤特征，已有研究关注缺氧及缺氧相关基因在NB发生、发展中的作用。但缺氧状态对于NB免疫环境及患者预后的影响有待进一步探究。本研究通过共识聚类算法，依据MSigDB数据库中的缺氧相关基因集将GSE62564中的498位患者分为高、低缺氧分组，并比较两组间临床特征及免疫环境的区别。本研究中，不同缺氧状态的患者，其免疫微环境有显著差异。接着，使用LASSO回归算法建立了由CSRP2，DTNA，SAP30，NCAN，WSB1，NAGK六个缺氧相关基因组成的预后风险模型。该模型具有良好的预测性能，并可作为NB的独立预后因素。最后，将风险评分与临床信息结合，构建列线图模型，进一步提升预测效率。

本预后模型由6个缺氧相关基因组成。CSRP2属于富含半胱氨酸和甘氨酸蛋白(cysteine and glycine rich protein, CSRP)家族，编码LIM结构域蛋白，参与调节发育和细胞分化过程。CSRP2在多种肿瘤中起调控作用。在乳腺癌中，HIF-1可通过靶向CSRP2，促进乳腺癌细胞的侵袭 [21] 。而在结肠癌中，CSRP2可通过HIPPO信号通路抑制肿瘤细胞增殖，且CSRP2的低表达与结肠癌患者的不良预后相关 [22] 。DNTA主要编码肌营养不良蛋白相关蛋白复合物(dystrophin-associated protein complex, DAPC)中的支架蛋白，主要作用是维持心肌、骨骼肌的结构完整性 [23] 。Hu等 [24] 在研究中发现，DNTA可通过调控STAT3、TGFβ1和P53信号通路促进乙肝病毒相关肝癌的进展。在食管癌中，microRNA-301b-DNTA轴可影响肿瘤细胞的生长、侵袭和迁移 [25] 。SAP30编码的蛋白是组蛋白脱乙酰酶复合物的组成成分，已被报道可通过UHRF1-SAP30-MXD4轴影响急性髓性白血病的发展 [26] 。Tang等 [27] 的研究建立了肝细胞癌缺氧相关预后模型，其中SAP30可作为肝细胞癌不良预后的标志物。这与本研究中SAP30的趋势相同。神经蛋白聚糖(neurocan, NCAN)主要参与调节细胞的粘附和迁移。在NB中，NCAN可促进肿瘤细胞的增殖与迁移，并与NB患者的不良预后显著相关 [28] 。此外，NCAN也可在肝等非神经组织中表达，并可能是酒精性肝病引起肝癌的风险因素 [29] 。此前多项关于NB预后模型的研究也将NCAN纳入其中，说明其对NB预后的影响有较大的意义，值得进一步研究 [30] [31] 。NAGK (N-acetylglucosamine kinase)属于乙酰己糖胺激酶家族。在NB患者中，血浆MYCN与NAGK的比值可用于检测MYCN的扩增情况，并可用于判断肿瘤负荷及预后情况 [32] [33] 。此外，NAGK可通过干预细胞内糖代谢影响胰腺癌肿瘤生长及患者预后 [34] [35] 。WSB1属于E3泛素连接酶，参与了多种细胞过程，如甲状腺激素调节、免疫调节和调节细胞对缺氧的反应等 [36] 。据研究，WSB1可通过调控ATM激酶降解来促进肿瘤进展 [37] 。WSB1可通过介导pVHL降解促进肿瘤转移。这在Poujade等关于乳腺癌的研究中得到证实 [38] [39] 。据Keren Shichrur等 [40] 的研究，沉默WSB1可抑制NB细胞的生长。但Chen等 [41] 的研究中指出，WSB1的高表达与较好的临床结局相关，这与本研究相符。因此，WSB1在NB中的作用有待进一步研究。

综上所述，本研究基于GEO数据库，分析了缺氧状态与NB免疫环境及临床特征的关系，并构建了由6个缺氧相关基因组成的预后风险模型，以期更好地判断NB患者的预后，并为NB的治疗提供新的思路。

基金项目

山东省自然科学基金资助项目(No. ZR2020MH213)。

NOTES

^*通讯作者。