基底膜相关基因风险模型在胃癌中的分析应用

李鸿昊，许鸣,，韩美霞，邓媛媛，黎丽旋,

【提要】目的 分析基底膜(BM)相关基因在胃癌中的表达意义。方法 从公共数据库下载BM相关基因、胃癌转录组数据与临床信息；基于TCGA筛选BM差异基因后进行COX和Lasso回归分析构建风险模型，分析风险评分与生存预后以及临床特征之间的关系， GSE26253队列验证其预测生存预后的准确性。功能富集、GSVA分析风险基因潜在的作用通路。计算免疫细胞浸润与风险评分之间的关系，基于风险基因模型发掘胃癌相关的免疫检查点。结果 基底膜相关基因高风险评分是胃癌患者生存预后的独立危险因素；风险基因与免疫功能密切相关(P<0.05); 风险基因相关的免疫检查点及相关通路存在关联(P<0.05)。结论 基底膜相关风险基因模型在胃癌中具有重要的临床应用价值。

据最新调查显示，胃癌的发生率和死亡率均位列于所有肿瘤类型的前列[1-2]。近些年来，随着高通量测序技术不断进步，越来越多的基因在胃癌中的表达作用被逐渐挖掘[3]。多数胃癌患者在被诊断时已处于中晚期阶段，预后较差。因此，积极探索与胃癌相关的潜在预测及预后基因具有重要的价值[4-5]。基底膜 (BM) 是最古老的动物细胞外基质 (ECM)，其周围围绕大多数组织和上皮细胞，形成片状结构。其中包含层粘连蛋白和胶原分子，对于细胞的粘附和转移起到至关重要的作用[6]。研究表明BM相关基因与肿瘤细胞脱落、转移存在紧密关联因此，理解BM相关基因在肿瘤的发展进程中具有重要意义[7]。

本次研究基于公共数据库的发掘，首次从BASE数据库中下载BM相关基因[8]，研究其在胃癌中的表达意义。首先筛选在胃癌中差异倍数较大的BM相关基因，构建风险模型。基于风险模型对患者进行风险评分，从而评价风险模型在胃癌中的临床应用。同时采用基因功能富集分析其在生物学功能方面的应用，以及探讨风险模型与胃癌免疫系统之间的关系，发掘潜在的免疫检查点，旨在临床方面为胃癌的治疗寻求新的思路。

1 材料与方法

1.1 使用数据库网站

TCGA(https://www.cancer.gov/)；BASE (https://bmbase.manchester.ac.uk)；GEO(https://www.ncbi.nlm.nih.gov/geo)；MSigDB(https://www.gsea-msigdb.org/gsea/msigdb)；GEPIA(http://gepia.cancer-pku.cn/)；KM-plot(http://kmplot.com/analysis/)。

1.2 数据下载

从TCGA官网STAD-counts平台下载胃癌转录组数据及临床信息。从GEO数据库下载GSE26253芯片数据。从BASE数据库中下载BM相关基因。

1.3 筛选基底膜相关差异基因

基于TCGA队列基因表达数据，利用R中limma包进行数据整理，筛选BM相关差异基因。设置条件为|LogFC|>2，矫正后的P值FDR<0.05。使用R对差异基因进行单因素COX回归预测，对具有预后价值的差异基因进行LASSO-Cox回归分析构建风险预后模型。

1.4 建立风险预后模型

基于上述差异基因，采用 R 软件 glmnet包进行 LASSO-Cox回归分析，构建风险模型。风险模型通过以下公式计算：风险模型=(系数1 * mRNA1基因表达量) + (系数2 * 基因表达量 mRNA2) + (系数n * mRNAn基因表达量)。 根据风险模型表达公式，对TCGA数据库中每个胃癌病人进行风险评分。根据风险评分的中位数将患者分为高风险组和低风险组，绘制Kaplan-Meier生存曲线。此外，分析TCGA队列患者临床特征与风险评分之间的相关性。对GSE26253芯片数据进行整理，得到433个胃癌术后患者的基因表达数据和术后无瘤生存时间，使用风险评估模型对每个胃癌患者进行风险评分，根据风险评分的中位数将患者分为高风险组和低风险组，绘制Kaplan-Meier生存曲线。

1.5 受试者操作特征(ROC)与列线图

利用R中survivalROC包绘制受试者工作特征(ROC)曲线判读风险模型预测患者第1、3、5年生存率的敏感性和特异性，评估风险模型预测的准确性。PCA 分析探讨高低风险组病人的分布范围。使用临床特征与风险模型构建列线图，预测胃癌患者第1-，3-，5-年的生存率，C指数与校准曲线评估列线图的准确性。

1.6 基因功能富集分析与GSVA

通过R中ClusterProfiler包进行GO分析，包括分子功能(MF)，生物学途径(BP)，细胞组分(CC)。并使用同样的方法计算KEGG富集通路分析。筛选条件为LogFC>1，FDR<0.05。GSVA进一步探索风险基因潜在的作用通路。

1.7 免疫浸润与免疫检查点

肿瘤的进展与免疫浸润密切相关，我们对于免疫浸润在高低风险组的病人进行了进一步分析。CIBERSORT计算基于风险模型评分分组的高低风险组病人的免疫细胞浸润水平间的差异，分析免疫功能与风险评分之间的关系。另从文献中整理已发表的免疫检查点，用于分析风险模型与免疫检查点的关联。

1.8 统计学分析

所有的统计学分析基于R软件包(R 4.2.1.)计算，P<0. 05表示差异具有统计学意义。

2 结果

2.1 从BASE数据库下载224个BM相关基因

从TCGA中下载胃癌的转录组数据集临床信息(临床信息见表格1)，得到375例胃癌组织及32例正常组织样本。从GEO中下载了从GEO数据库下载GSE26253芯片数据，得到433个胃癌患者组织样本的基因表达数据及胃癌患者术后随访时间。

表1 TCGA胃癌转录组的临床样本特征

2.2 筛选BM相关差异基因

如表2所示，在TCGA队列中，通过与正常组织做比较，在224个BM相关差异基因中，共筛选出22个差异基因(|LogFC|>2， FDR<0.05)，其中COL4A6、UNC5D、AMTN在胃癌组织中表达下调，19个基因在胃癌组织中表达上调。

表2 TCGA-STAD队列中差异表达4倍及以上的BM相关差异基因

图1 胃癌患者的BM相关差异基因的筛选 A：22个BM相关基因构成的基因热图；B：22个BM相关基因构成的火山图

2.3 构建基底膜相关基因风险模型

使用R survival对上述22个差异基因进行单因素COX回归预测(图2A)，对具有预后价值的差异基因进行LASSO-Cox回归分析构建了BMS-GE风险模型(图2B和2C)：风险模型=(VCAN*0.13435+ADAMTS18*0.06122)。采用风险模型计算公式，在TCGA-STAD队列中，对每个样本进行风险评分，以风险评分中位数为临界值，将队列患者分为高风险组(n=185)和低风险组(n=186)，Kaplan-Meier曲线显示高风险组生存时间明显缩短(图3A)；使用风险评估模型对GSE26253队列中样本进行风险评分，以风险评分中位数为临界值，结果显示高风险组患者无瘤生存期明显低于风险组患者(图3B)。分析风险模型与临床特征的相关性，结果显示：G3组风险评分明显大于G2组；N3组风险评分大于N1组风险评分；StageⅡ/StageⅢ组风险评分均大于StageI组风险评分；T2/T3/T4组风险评分均大于T1组风险评分，T4组评分大于T3组风险评分(P<0.05)。

图2 COX-LASSO回归分析 A：COX回归森林图 22个基底膜相关差异基因进行分析得到4个具有预后价值的基因；B：Lasso回归选择最合适的λ值，4个危险基因的Lasso系数路径图；C：Lasso回归交叉验证

图3 风险模型与胃癌生存预后及临床特征之间的关系 A.TCGA队列中高风险组的总生存率明显短于低风险组患者;B.GSE26253队列中，高风险组的术后无瘤生存率明显短于低风险组患者。C：风险模型与TCGA队列临床特征之间的关系

2.4 ROC 曲线和PCA分析

绘制ROC曲线结果显示：1年时的曲线下面积(AUC)为0.606，3年时AUC为0.608，5 年时AUC为0.773(图 4A)。PCA分析能够准确区分高低风险组的分布范围(图4B)。结果提示BM相关基因风险模型在预测胃癌患者的生存预后方面具有一定价值。

图4 ROC曲线与主成分分析 A：风险模型评价ROC曲线，AUC曲线下面积均大于60%；B：胃癌患者的生存时间及生存状态；C：PCA分析能够准确区分高低风险组患者的分布

2.5 独立预后分析

COX回归分析风险评分是否能够成为独立预后危险因素。单因素结果显示年龄、肿瘤分期、风险评分为生存预后相关因素(图5A)；多因素结果显示年龄、 肿瘤分期、风险评分为独立预后危险因素(图5B)。

图5 对胃癌临床特征及风险评分急性COX回归分析 A：单因素COX回归森林图;B：多因素COX回归森林图

2.6 构建列线图

基于TCGA队列，纳入多个预后因素(年龄、性别、分级、分期以及T分期、M期、N分期、风险评分)，采用风险模型计算公式对患者进行风险评分，预测胃癌患者的生存率。从图中可以看出第1年、5年的预测能力与预测曲线几乎重合，第3年的预测能力与预测曲线略有偏差(图6A)。校准曲线评估列线图的准确性，C指数为0.681，提示列线图具有较准确的预测能力(图6B)。

图6 列线图评估胃癌患者预后生存率 A:列线图预测第1-，3-，5-的生存率；B:校准曲线预测第1-，3-，5-的生存率

2.7 功能富集分析

GO分析显示来自TCGA队列风险组间的差异基因主要富集于含胶原细胞外基质组织和细胞外部组织结构(表3和图7)。KEGG 途径分析显示DEG主要富集在Neuroactive ligand-receptor interaction(脑组织神经活动配体-受体相互作用通路信号通路)及PI3K-Akt signaling pathway(丝氨酸/苏氨酸激酶信号通路)(表4和图7)。

表3 GO富集分析

表4 KEGG通路富集分析

图7 GO富集分析与KEGG富集分析预测风险基因富集的相关通路 A：GO富集分析圈图；B：GO富集分析柱状图； C：GO富集分析气泡图(BP代表生物学过程；CC代表细胞成分；MF代表分子功能 D：Kegg富集分析圈图；E：Kegg富集分析柱状图；F：Kegg富集分析气泡图

2.8 GSVA

GSVA进一步阐明基于风险组的特征的分子机制。在TCGA队列中，GSVA结果显示TGFβ信号通路、粘着斑激酶信号通路、间隙连接通路、糖胺聚糖生物合成硫酸软骨素通路主要富集在高风险组。药物代谢-其他酶通路、视黄醇代谢通路、抗坏血酸和醛酸代谢通路、戊糖和葡萄糖醛酸相互转化通路主要富集在低风险组。

图8 GSVA/Kegg分析结果

表5 GSVA分析结果前十条通路集合

2.9 风险模型与免疫细胞浸润之间的关系

根据CIBERSORT算法，计算高低风险组评分与免疫浸润之间的关系(图9)。CIBERSORT结果显示M2巨噬细胞、单核细胞在高风险组中呈现高表达(P<0.05)。此外我们进一步探索了高低风险组与免疫相关功能之间的关系，趋化因子受体(CCR)在高风险中高表达，副炎症相关反应(Parainflammation)、T细胞抑制功能(T_cell_co-inhibition)、I型干扰素相关功能在高风险组中高表达。考虑到免疫检查点的重要性，我们进一步探索了高低风险组与免疫检查点之间的关系,我们发现在高风险组中，多个免疫检查点呈现高表达趋势(ADORA2A、BTLA、CD27、CD28、CD40、CD40LG、CD44、CD48、CD8、CD86、CD400、CD200、CD200R1、CD244、CD274、CD276、CTLA4、HAVCR2、ICOS、LAIR1、NRP1、PDCD1LG2、TIGIT)，其中包括TNF超家族/TNF受体超家族的多个成员(TNFRSF4、TNFRSF8、TNFRSF9、TNFSF4、TNFSF14、TNFSF18 TNFRSF5/CD40，TNFRSF7/CD27)。

图9 风险基因模型与免疫浸润之间的关系 A：免疫细胞浸润在高低风险组间的差异；B：免疫功能在高低风险组间的差异

2.10 ADAMST18在胃癌中具有一定研究价值

GEPIA数据库与Km-plotter数据库均提示在胃癌中ADAMST18高表达预示较差的生存预后。

图11 A:Km-plotter数据库;B:GEPIA数据库

3 讨论

肿瘤的转移是恶性肿瘤的基本生物学特征，是临床上绝大多数肿瘤患者生存质量较差的原因[9]。肿瘤细胞脱落侵入基质是肿瘤转移的重要环节，肿瘤细胞从原发灶脱离涉及多种复杂基质，多项研究表明基底膜在肿瘤细胞脱落过程中充当十分重要的介质[10-11]。胃癌是最常见的恶性肿瘤类型之一，预后较差。因此我们对基底膜相关基因在胃癌中的表达意义做出研究。

基于TCGA数据库的研究，我们首次采用COX回归及Lasso回归分析方法成功构建了包含VCAN、ADAMTS18两种基底膜相关基因的风险模型。使用GSE26253队列中患者进行风险评分验证，得出高风险评分预示较差预后生存。构建列线图及ROC曲线发现风险模型可以成为胃癌的独立预后风险因素且可以预测胃癌患者的生存期，此外风险评分与胃癌的分期、分级等临床特征之间相关，具有临床应用价值。

研究表明VCAN在消化系统肿瘤中具有重要的表达意义，其与TGF-β信号通路、EMT通路等相互作用，参与肿瘤细胞的血管生成、免疫抑制等多种生物学功能从而促进肿瘤的迁移与转移，在胃癌中的研究也相对较多，同时有研究表明VCAN可作为胃癌的潜在治疗靶点[12-14]。据报道，ADAMTS是一类具有锌离子依赖性的分泌型基质金属蛋白酶，其家族成员在不同类型肿瘤的形成过程包括肿瘤细胞的侵袭、转移发挥作用以及与部分重要的癌症相关的信号通路之间存在紧密关联，ADAMTS18在食管癌、肝癌、肺癌、肾癌等均具有表达意义，但在胃癌中的研究报道相对较少[15-16]。我们使用KM-plotter数据库与GEPIA数据库初步探索了单基ADAMTS18在胃癌方面的生存预后价值，发现ADAMTS18高表达提示较差的胃癌生存预后，由此可见ADAMTS18在胃癌中具有潜在的研究价值，值得进一步探讨。

进一步研究风险基因模型在生物功能学方面的作用时，我们发现风险基因主要富集于含胶原细胞外基质组织和细胞外部组织结构。GSVA显示TGF/β通路、粘着斑激酶信号通路、缝隙连接通路参与胃癌发生发展，以上通路均与肿瘤细胞的脱落、侵袭与迁移紧密关联[17-18]。KEGG 途径分析显示风险基因主要富集在脑组织神经活动配体-受体相互作用通路及PI3K-Akt信号通路。EMT过程作为肿瘤微环境的重要组成成分，是肿瘤加速进程的重要中间过程。研究表明，PI3K-Akt信号通路参与多种类型肿瘤的发展，其是EMT过程中的重要介质，肿瘤上皮细胞经历EMT过程后，细胞间的粘附能力下降、侵袭能力增强，促进肿瘤细胞随血液循环运动，加速肿瘤向全身转移[19-20]。肿瘤细胞、微环境与宿主在动态演变中相互联系、相互促进，加快肿瘤的发展。在这一过程中，细胞外基质充当重要角色，肿瘤的运动侵袭性、调节肿瘤微环境的能力以及在侵袭组织中定植生长均与基质相互关联，降解细胞外基质是肿瘤侵袭和转移的必经阶段[21]。按分布部位来说，基底膜属细胞外基质的一部分。通过以上的分析，因此我们理论上认为基底膜相关风险基因或许通过影响细胞外基质的行为功能而影响肿瘤微环境的变化促进胃癌的发生发展。此外，我们发现了脑组织神经活动配体-受体相互作用通路在胃癌中有所涉及，而脑组织神经活动配体-受体相互作用通路多在神经系统疾病包括脑神经胶质瘤中报道，在胃癌中少有报道[22-23]。早有研究报道胃部出现异常神经可导致胃癌的形成[24]，由此可见神经系统与胃癌的发生具有一定关联，脑组织神经活动配体-受体相互作用通路在神经系统疾病发挥重要作用，那么该通路是否在胃癌的发展进程中具有重要作用值得日后进一步研究。

免疫微环境作为肿瘤微环境有机整体的主要组成成分，是肿瘤领域研究的热点，同样在我们的研究中发现，M2巨噬细胞、单核细胞在高风险组中呈现高表达趋势，已知巨噬细胞、单核细胞的预示胃癌患者的不良预后，这可能是高风险组预后差的重要影响因素；在免疫功能方面，高风险组中副炎症相关反应(Parainflammation)、T细胞抑制功能(T_cell_co-inhibition)、II型干扰素相关功能高表达[25-26]。以上说明风险模型在胃癌免疫细胞浸润中发挥十分重要的作用。肿瘤坏死因子受体超家族 (tumor necrosis factor receptor superfamily, TNFRSF)是一类可以通过半胱氨酸富集区域 (cysteine-rich domains, CRDs)与肿瘤坏死因子超家族 (tumor necrosis factors superfamily, TNFSF)结合的细胞因子受体超家族，TNF超家族/TNF受体超家族相互作用在多种类型肿瘤中的发生进程中发挥重要角色，其与免疫系统及免疫疾病的形成中也发挥重要角色[27]。当我们探索风险模型与免疫检查点之间的关系时，我们发现在高风险组中，TNF超家族/TNF受体超家族的多个成员均呈现高表达。此外，在高表达风险组中，多种CD分化抗原呈现高表达趋势，我们知道CD分化抗原多分布在一些免疫细胞上，在肿瘤进程中通过调节免疫系统成分而改变肿瘤微环境，促进肿瘤的发生发展。如 CD27作为肿瘤坏死因子受体超家族一员，表达在T 细胞、髓质胸腺细胞、B 细胞亚群和自然杀伤细胞，共刺激 B 细胞和 T 细胞激活，从而发挥重要的免疫功能[28]。因此我们推测风险模型可与TNF超家族/TNF受体或CD分化抗原存在某种关联，共同促进肿瘤的发生发展。

综上所述，基于数据库建立的基底膜相关基因风险模型对胃癌的生存预后以及其与临床特征之间的关系具有重要的临床应用价值。在生物功能学方面，通过以上结果分析，我们推测，基底膜风险基因或许可通过与多条通路相互作用调节肿瘤微环境的变化而促进胃癌的发生发展。最后我们发现基底膜风险基因与胃癌的多个免疫检查点相互作用，尤其与CD分化抗原与TNF超家族/TNF受体超家族紧密关联，而CD分化抗原与TNF家族均是免疫系统的重要组成成分，基底膜风险基因是否是免疫系统的重要组成成分有待实验进一步探索。脑组织神经活动配体-受体相互作用通路与ADAMTS18在胃癌领域具有潜在的研究价值。但本文的不足之处在于本研究基于公共数据库进行分析，生物功能学分析多是基于已有文献报道进行分析，缺乏一定的实验支撑及机制探索，我们拟在未来会在实验基础上对机制进行进一步的探讨。