全球范围内,颈椎病(cervical spondylosis,CS)患病率随人口老龄化及现代生活方式的改变而显著上升
[1],同时抑郁症(major depressive disorder, MDD)已成为导致残疾调整生命年的主要疾病之一
[2]。临床发现,CS患者的抑郁症患病率显著高于普通人群
[3]。另有研究表明,MDD患者的CS患病率高于普通人群
[4]。循证医学证据表明CS与MDD显著相关
[5]。利用遗传变异推断因果关系的研究也证实了二者的双向因果关联
[6]。提示二者的显著关联可能共享部分病理生理基础。然而,目前对该现象的理解仍停留于观察和推断层面,缺乏系统的分子机制研究。
现有研究虽尚未对CS与MDD共病现象提出研究假说,但在慢性疼痛与MDD关联的不同角度下提出了免疫假说
[7]、神经内分泌假说
[8]、神经假说
[9]等。尽管如此,CS与MDD之间的分子机制研究在系统基因组层面证据、关键信号通路和分子网络、可用于临床诊断的生物标志物等方面存在空白。鉴于此,现通过整合多种生物信息学方法,对CS和MDD的基因表达数据集进行分析,系统探索两种疾病的共同分子特征。除利用传统生物信息学分析外,通过机器学习算法筛选潜在的诊断生物标志物。从分子水平揭示CS与MDD共病的生物学基础,为理解这两种疾病的共同发病机制提供新的视角。
1 材料与方法
1.1 数据来源及预处理
从基因表达综合数据库(Gene Expression Omnibus, GEO,网址:
http://www.ncbi.nlm.nih.gov/geo)检索关键词“Cervical spondylosis”和“Depressive disorder”,获取两个独立的基因表达数据集用于整合分析。GSE98793数据集基于GPL570平台,包含MDD患者与健康对照者的样本。GSE153761数据集则基于GPL22120平台,包含CS相关疾病患者与健康对照者的样本。预处理过程则根据各平台特性分别进行背景校正、分数位标准化和表达量汇总,根据平台注释信息进行mRNA探针进行筛选,并进行基因注释与除重。随后匹配样本临床信息(分组、性别、年龄)与表达矩阵,确保后续分析准确性。
1.2 差异表达与交集分析
采用limma包对每个数据集分别进行差异表达分析。在线性模型中,纳入疾病分组、年龄和性别作为协变量以校正潜在影响,并利用经验贝叶斯方法进行推断统计。初始阶段优先采用统一常规阈值(
P<0.05且
log
2FC
>1(其中FC为差异倍数)进行差异表达基因 (differentially expressed genes, DEG)筛选,以保证分析标准一致性;实际分析发现,MDD (GSE98793)经该阈值筛选后无显著差异表达基因,无法支撑后续共病核心基因分析,CS(GSE153761)可稳定筛选差异基因。为在保证数据可靠性的前提下推进共病分析,MDD调整阈值为
P<0.05且
log
2FC
>0.25,同时补充敏感性分析验证该阈值可靠性:结果显示,筛选出的差异基因经GO(gene ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析可显著富集于MDD已知核心通路(神经炎症通路等),无随机富集现象,排除微弱变化基因过度筛选导致的偏差。
1.3 共表达网络分析及蛋白质互作网络分析
由于GSE153761未满足11例样本的最低要求
[10],因此仅对GSE98793表达数据进行加权基因共表达网络分析(weighted gene co-expression network analysis,WGCNA)。为识别在MDD中高度协同表达的基因模块,对GSE98793表达数据进行 WGCNA。筛选变异系数前50%的高变异基因构建网络。软阈值功率设定为5,以确保网络符合无尺度拓扑结构。采用动态剪枝法识别基因模块,最小模块大小为30。通过计算模块特征向量与疾病性状的相关性,鉴定与MDD最相关的关键模块,并提取其基因。
同时构建蛋白质-蛋白质相互作用网络(protein-protein interaction,PPI)。将双疾病交集基因、WGCNA关键模块基因等基因构建为候选基因集,通过biomaRt映射为蛋白编码基因,并利用检索相互作用基因/蛋白质的搜索工具STRING数据库(v12)进行映射。采用互作得分≥400或≥150的双阈值策略平衡网络的可靠性与完整性。并计算网络拓扑属性,并将度中心性排名前10%的节点定义为Hub基因。
1.4 多重异构网络和随机游走重启分析
在完成上述多层分析的基础上,构建整合多维度生物学信息的多重异构网络。该网络包含三个层次:①PPI互作层,边权重统一设置为1.0;②WGCNA共表达层,为同一模块内的基因对添加边,并赋予0.8的权重;③GO功能相似层,基于GO 生物学过程(biological process, BP)注释并计算基因对间的功能相似性,对得分超过0.6的基因对添加边,其权重即为相似性得分。在该网络的基础上应用随机游走重启(random walk restart,RWR)算法,以差异表达分析中的种子基因作为起点,模拟信息流动过程。算法重启概率设置为0.7,最终输出网络中每个节点的稳态概率(RWR得分),概率高的基因,被认为可能在连接两种疾病机制中存在可调控作用。
1.5 网络拓扑与社区结构分析
为系统分析网络拓扑属性,通过计算包括节点数、边数、密度、平均度、平均路径长度和聚类系数在内的全局参数以量化其连接模式与效率,并通过检验度分布的幂律特性来评估无尺度性。随后采用Walktrap算法对网络进行社区结构检测,识别内部连接密集的基因模块,并以模块度指标评估划分质量。社区结构分析为探索基因模块的功能聚类提供参考,后续将结合富集分析进一步解析。
1.6 机器学习筛选潜在关键基因
为进一步量化基因的共病风险,采用弹性网络回归模型进行机器学习建模。每个基因由多维特征向量表示,包括:|log2FC|、-lgP、WGCNA模块重要性、PPI度中心性、PPI中介中心性以及RWR得分。目标变量为这些特征归一化后的加权综合评分。通过5折交叉验证选择最优正则化参数λ(即惩罚系数,用于控制模型收缩强度),并基于模型系数评估各特征的重要性。根据预测得分对所有基因进行排序,排名靠前的基因被视为潜在的共病关键基因。
1.7 功能与富集分析
对机器学习模型筛选出的潜在关键基因进行GO、KEGG及Reactome通路富集分析,以揭示经过多组学数据整合与机器学习算法优先筛选出的高阶候选基因所参与的更精确的生物学过程和通路,从而为其在共病机制中的角色提供功能学解释。
1.8 统计分析和可视化
所有统计分析均在R语言(v4.5.1)环境完成。差异表达分析使用limma包,WGCNA使用WGCNA包,网络分析使用igraph和ggraph包,富集分析使用clusterProfiler包。可视化方面,火山图、气泡图、条形图等使用ggplot2包绘制;热图使用pheatmap包绘制;网络可视化使用igraph和ggraph包完成;和弦图使用circlize包生成。
2 结果
2.1 抑郁症与颈椎病共享差异表达基因的鉴定与分析
通过对GSE98793和GSE153761数据集开展DEG分析,分别在两个数据集中鉴定出330个和568个显著DEG。进一步的交集分析显示,二者共存在22个共享基因[
图1(a)], 其表达方向如
表1所示。该结果提示共享基因在两种疾病中的表达方向存在明显异质性:仅1个基因在两者中均上调,1个均下调,其余20个基因的表达方向在抑郁症与颈椎病中不一致。火山图[
图1(b)、
图1(c)]直观展示了各数据集中差异表达基因的分布情况,其中红色标记点表示达到显著性阈值的基因,为后续机制探索提供了清晰的可视化基础。
2.2 共表达网络与蛋白质互作网络揭示核心功能模块与枢纽基因
基于GSE98793数据集构建的WGCNA分析识别出7个核心基因模块。通过模块-性状相关性分析(
图2)可知,灰色模块与MDD相关性最为显著(相关系数=0.26,
P<0.05),因此被确定为关键模块,共包含1 510个基因。随后,将该模块映射至STRING数据库,构建蛋白质互作网络以解析其内部的基因互作关系。该网络共包含56个节点和59条相互作用边。PPI网络(
图3)识别出连接度最高的6个枢纽基因,分别为
PGLYRP1、
CEBPE、
CCDC80、
LAG3、
PKP2和
FAP,提示它们在MDD与CS共病的分子调控中可能处于核心地位。
2.3 多重异构网络与随机游走重启算法定位共病关键基因
通过整合PPI、WGCNA及功能相似性信息等多维度结果,构建包含10 529个节点的多重异构网络[
图4(a)]。并且在该网络上执行RWR分析,算法在迭代8次后迅速收敛。而基于稳态概率分布进行排序,显示
CEACAM8、
ADORA3、
OLFM4、
NAST1、
KLRG1和
TNFRSF17等基因的RWR得分较高[
图4(b)]。其中,
PGLYRP1和
DEFA4等高得分基因在网络可视化中处于重要枢纽位置[
图4(a)],且连接了不同功能模块。
2.4 机器学习模型整合多组学特征并量化基因共病风险
为筛选高置信度共病关键基因,通过整合差异表达倍数变化、WGCNA模块重要性、PPI网络中心性及RWR得分等维度特征构建弹性网络回归模型。基于该模型对候选基因进行风险排序,结果[
图5(a)]显示,
LTF、
OLFM4、
CEACAM8、
DEFA4及
LCN2等基因评分较高。特征重要性分析[
图5(b)]表明,RWR得分、差异表达倍数变化和PPI节点度是共病预测中贡献最大的三个特征,而DEG的
P值及WGCNA模块重要性的影响相对较小。且模型展现出良好预测性能,实际综合评分与预测风险评分之间呈高度线性拟合[
图5(c)],且均方误差(mean squared error,MSE)低至0.000 009。
2.5 网络拓扑与社区结构解析揭示共病网络的模块化特性
对构建的多重异构网络进行拓扑属性分析,该网络呈现出典型的稀疏特性,具有10 580个节点和398条边,网络密度为0.000 007。社区结构分析将网络划分为10 206个社区,模块度高达0.876,表明该共病分子网络具有高度的模块化结构。最大的连通组件包含51个基因。
2.6 功能富集分析阐明关键基因参与的生物学通路
对高风险基因进行多维度功能富集分析,结果显示,高风险基因显著富集于免疫与炎症相关的生物学过程及信号通路。KEGG分析[
图6(a)]表明,PI3K-Akt 信号通路、Ras 信号通路、钙信号通路、丝裂原活化蛋白激酶(mitogen-activated protein kinase,MAPK)信号通路、Rap1 信号通路以及中性粒细胞胞外陷阱形成等通路均表现出较高的基因富集数量与统计显著性。Reactome通路分析[
图6(b)]进一步表明,中性粒细胞脱颗粒、抗菌肽通路不仅包含较多基因,其富集强度显著。同时,细胞外基质降解、基质金属蛋白酶激活、补体级联调控等过程也被明显富集。
GO生物学过程分析[
图7(a)]显示,纳入基因主要参与抗菌体液免疫反应、对细菌 / 真菌的防御反应、趋化作用、细胞杀伤等免疫相关过程,在网络中形成密集的相互作用,核心免疫功能模块高度关联。Reactome通路互作网络分析[
图7(b)]显示,中性粒细胞脱颗粒、抗菌肽、细胞外基质组织、补体级联调控等通路之间存在紧密连接,通过多个关键基因形成协同调控网络。
2.7 多维度分析高频基因的综合验证与共病核心地位确认
为明确在多维度分析中持续发挥作用的核心基因,对MDD与CS差异表达基因筛选(
DEG_MDD、
DEG_CS)、WGCNA关键模块基因(
WGCNA_Key)、PPI枢纽基因(
PPI_Hub)、RWR中的Top传播基因(
RWR_Top)、机器学习的共病高风险基因(
ML_Risk)、多重异构网络中的Top基因(
Feature_Top)出现的潜在共病基因进行交叉统计,筛选前10个高频基因并可视化其分布(
图8)。结果显示,
MMP8、
CEACAM8和
OLFM4在各分析步骤中均稳定出现,且均被聚类于社区1中,结果呈现一致性。
RETN及
HP也在各重要基因中多次出现。
3 讨论
近90%的人到50岁时都将出现椎间盘退行性变,故CS正呈发病率高、病程长的趋势发展
[11]。除引起颈部疼痛、神经功能障碍等症状外,CS可引起焦虑、抑郁等并发症状
[12]。虽已有研究证实颈部疼痛存在抑郁、焦虑等精神病理风险因素
[13],但二者在分子层面的共病机制探索仍然有限。整合多种生物信息学分析方法,对MDD与CS的共病分子网络进行系统分析,鉴定多个在共病机制中持续发挥核心作用的候选基因。通过多重异构网络构建与RWR算法的联合分析,从10 529个节点的复杂网络中有效识别出了具有较强影响力的核心基因,机器学习模型也表明RWR得分、差异表达倍数变化和PPI节点度是本次研究贡献较大的特征。多维度分析策略在既往生物信息学分析方法的基础上进一步深入
[14-15],为揭示二者在分子水平上的交叉调控提供了重要依据。
MMP8是稳定表现的核心基因之一。该基因属基质金属蛋白酶家族,其主要功能是降解细胞外基质中的Ⅰ、Ⅱ和Ⅲ型胶原
[16]。椎间盘髓核组织以Ⅱ型胶原为主要成分,
MMP8可通过降解Ⅱ型胶原引起组织退行性变
[17]。研究表明,外泌体携带的
MMP8参与椎间盘退变的发生与进展
[18]。Ren等
[19]研究表明,
MMP8抑制剂通过调控肿瘤坏死因子-α/肿瘤坏死因子受体1/核因子κB(tumor necrosis factor-alpha/nuclear factor kappa-B, TNF-α/TNFR1/NF-κB)信号通路抑制NOD样受体热蛋白结构域相关蛋白3(NOD-like receptor family pyrin domain containing 3, NLRP3)炎症小体激活等多重机制,有效减轻慢性不可预见性轻度应激(chronic unpredictable mild stress, CUMS)诱导的神经炎症和抑郁状态。此外,
MMP8还可通过激活蛋白酶活化受体,进而调控KEGG富集中显著富集的MAPK、Ras和Rap1等信号通路,促进多种炎症因子的释放
[20]。已有研究表明,MDD患者血清中髓系来源的
MMP8表达显著升高,且
MMP8能够穿透血脑屏障进入中枢神经系统调节神经炎症反应
[21]。上述研究与Reactome通路分析中“中性粒细胞脱颗粒”“细胞外基质降解”等显著富集通路高度一致,提示
MMP8可能是连接CS外周组织损伤与MDD中枢神经炎症的关键分子桥梁。
CEACAM8在本次分析中表现同样稳定。该基因是在中性粒细胞表面广泛表达的一种蛋白,可介导中性粒细胞间的同型与异型黏附,在炎症部位促进免疫细胞的募集与活化
[22]。
CEACAM8在验证状态下过表达可过度激活中性粒细胞,促进促炎细胞因子的分泌,从而刺激局部血管扩张、通透性增加,形成“炎症级联反应”
[23]。该基因目前已被识别为多种炎症相关疾病的临床生物标志物,机器学习领域也已将其用于评估临床病理特征与疾病诊断之间的关系
[24]。多重异构网络显示,
CEACAM8在连接不同功能模块中处于重要的枢纽位置,提示其在MDD与CS共病中可能承担了外周免疫信号向中枢炎症传导的桥梁功能。
OLFM4功能与上述基因高度互补。该基因是NF-κB信号通路的靶基因,可通过结合含核苷酸结合寡聚化结构域蛋白1/2(nucleotide-binding oligomerization domain-containing protein 1/2, NOD1/NOD2)抑制其激活,形成负反馈调控路径,该机制与Ras/Rap1信号通路共同构成复杂的调节网络
[25]。
O’Sullivan等
[26]研究证实,
OLFM4是病理性炎症和终末器官损伤的标志物,也是介导这些病理过程的直接参与者。同时,OLFM4还可通过抑制中性粒细胞过度脱颗粒,减少
MMP8、抗菌肽、弹性蛋白酶等炎症介质的释放,避免炎症反应失控
[27]。该功能与Reactome通路富集分析中“中性粒细胞脱颗粒”“抗菌肽通路”的显著富集结果相吻合。然而,
OLFM4在MDD与CS中呈现相反的差异表达方向,有关神经与局部组织退行性变联系的具体生物机制值得深入研究。
同时,
HP和
RETN也多次出现在各分析步骤的高频基因列表中。研究指出,
HP是具有共病诊断价值的共享血液生物标志物,其表达与中性粒细胞浸润呈正相关
[28]。
RETN是巨噬细胞和炎症细胞分泌的炎症因子,可通过激活NF-κB等通路增加TNF-α、IL-6等促炎细胞因子的释放
[29],虽然现有研究主要聚焦于代谢性疾病领域,但抵抗素在慢性低度炎症和免疫调控中的关键作用提示其可能在MDD与CS共病的炎症微环境重塑中发挥一定作用。
功能富集分析进一步揭示了上述核心基因共同参与的关键生物学通路与免疫模块。中性粒细胞脱颗粒、抗菌肽反应、细胞外基质重塑以及补体级联调控等通路的显著富集,通过
MMP8、
OLFM4、
CEACAM8及
HP等基因形成了协同调控网络
[30]。研究指出,中性粒细胞不仅是终末效应细胞,还是一类具有高度适应性和组织重塑能力的动态免疫群体,能够参与无菌性炎症、组织修复乃至神经炎症调控
[31]。研究表明,中性粒细胞的脱颗粒程序具有刺激依赖性与时间分层特性,其异质性表达可精确调控局部炎症强度
[32]。为
MMP8、
OLFM4、
CEACAM8在中性粒细胞内激活所引起的级联反应提供了理论依据。基于上述分析结果,
MMP8、
CEACAM8和
OLFM4被确认为MDD与CS共病机制中最稳定的核心候选基因,它们可能通过调控中性粒细胞的功能极化,将外周组织退变引发的局部炎症信号转化为持续的低度全身性炎症,进而影响中枢神经系统的免疫微环境
[33],最终参与CS与抑郁症共病的发生发展。
PDF
(20033KB)
