帕金森病小脑枢纽基因及诊断模型

包成政，王功俊，罗雪莲，王洁，蓝学群，王喻，陈金梅，李雪斌,,4

帕金森病（Parkinson’s disease，PD）是第二大退行性病变，临床主要表现为静止性震颤、运动迟缓、肌强直、姿势步态异常等症状[1]。随着疾病进展，患者可因治疗效果差、并发症多等导致生活质量严重下降，甚至死亡。深入研究PD 的发病机制，寻找更多早期诊断标志物及潜在干预靶点使PD得以早诊断、早治疗是现阶段的热门研究方向。目前PD的研究主要集中在基底节区。有学者提出小脑有可能参与PD震颤的发生或发展过程[2]，也有学者认为小脑-丘脑-皮质回路有可能参与PD 震颤[3]。小脑的功能异常有可能是PD冻结步态的基础[4]，使用神经调节技术可改善PD 患者运动症状[5]，刺激小脑可以降低PD 患者左旋多巴治疗诱发的运动障碍的严重程度，小脑可能参与PD的疾病进程[6]。小脑在PD 的发病机制中具有重要地位，但目前相关研究较少，同时探索PD的诊断方式也是现行研究的热点。

本研究基于生物信息学方法，从GEO数据库获取PD小脑总RNA芯片GSE20314、GSE28894 的数据，通过对该数据进行limma差异分析和加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA），并通过VennDiagram 寻找枢纽（Hub）基因。再使用逻辑回归分析构建诊断模型，并通过2 个血液样本数据集进行验证，计算受试者工作特征曲线（receiver operator characteristic curve，ROC）的曲线下面积（area under the curve，AUC）评估其诊断性能。探索小脑在PD中的作用及构建具有应用价值的诊断模型。

1 资料与方法

1.1 数据采集

数据集来自NCBI 基因表达综合公共数据库（GEO）。GSE20314 和GSE28894 是在小脑中表达的RNA数据集，GPL96注释的GSE20314包括4个对照组人群和4 个PD 患者的小脑RNA 注释样本；GPL6104注释的GSE28894 包括14 个对照组人群和15 个PD 患者的小脑RNA注释样本。在分析前将2样本去批次及归一化，见图1。GSE18838和GSE6613是外周血RNA表达的2个数据集，其中GPL5175注释的GSE18838包括11个正常人样本和17个PD患者样本，GPL96注释的GSE6613包括22个正常人样本和50个PD患者样本。

图1 GSE20314、GSE28894箱线图和Density图

1.2 方法

1.2.1 差异表达分析 Limma是一种基于广义线性模型的差异表达筛选方法[7]。在分析前将GSE20314 和GSE28894 利用“inSilicoMerging”包进行合并，再使用经典贝叶斯方法进行去除批次效应[8]。我们使用R 软件limma包（version 3.40.6）进行差异分析。P＜0.05和差异倍数1.1 倍的基因被视为差异基因（DEGs）。DEGs的筛选及热图和火山图绘制使用Sangerbox在线分析工具（http://vip.sangerbox.com/home.html）。

1.2.2 加权基因共表达网络分析 基因之间共表达及基因与表型之间的关系探索使用在线分析工具Sangerbox（http://vip.sangerbox.com/home.html）。首先分别计算每个基因的MAD（median absolute deviation），剔除MAD 最小的前50%的基因。利用R 软件包WGCNA的goodSamplesGenes方法去除离群的基因和样本，进一步构建scale-free co-expression network。首先对所有配对基因分别采用Pearson相关矩阵和平均连锁法，然后用幂函数A_mn=|C_mn|^β（C_mn=基因m和基因n之间的Pearson相关关系，A_mn=基因m和基因n 之间的邻接关系）构造加权邻接矩阵（β是一个能强调基因之间强相关性和惩罚弱相关性的软阈值参数）。选择适当的幂后，将邻接关系转化为一个拓扑重叠矩阵（TOM），该矩阵可以度量一个基因的网络连通性，定义为该基因与所有其他基因的邻接关系之和作为网络基因比，并计算相应的不相似度（1-TOM）。为了将具有相似表达谱的基因划分为基因模块，根据基于tom的不相似度度量进行平均连锁层次聚类。为进一步计算模块特征基因的不同相似度，选择一条切线作为模块树状图，合并一些模块，并对每个模块赋予不同的颜色描绘。模块隶属度（MM）是基因表达值与模块ME 的相关系数，表示该基因与该模块的相关性。基因显着性（GS）是基因表达值与表型之间的相关系数，代表基因与表型之间的关联。

1.2.3 枢纽（Hub）基因鉴定 使用“VennDiagram”包获取DEGs 和WGCNA 获得的关键模块基因的交集，进而获得与PD较为密切相关的Hub基因，Hub基因在对照组人群和PD患者小脑中的表达使用箱线图表示。

1.2.4 逻辑回归模型 逻辑回归是一种广义线性回归分析模型，可用于疾病的自动诊断。本研究采用2 个响应变量的逻辑回归，响应变量为1 时表示PD 样本，响应变量为0时表示正常样本。逐步回归分析最初用于消除对响应变量不显著的因素，仅保留显著的因素以简化模型。逐步回归从模型中迭代地添加或删除变量，直到保留最显著因素。之后使用逻辑回归来拟合这些显著因素与响应变量之间的关系。最后使用受试者工作特征曲线（ROC）和ROC 曲线下面积（AUC）评估模型的诊断效果。

1.3 统计分析

所有分析在Sangerbox、SPSS 23.0、Graphpad Prism 9.4.0中进行，根据是否符合正态分布选择t检验或秩和检验；P＜0.05为差异有统计学意义。

2 结果

2.1 Limma差异基因

DEGs以P＜0.05、差异倍数1.1倍的基因作为筛选条件，共得到116 个DEGs，使用火山图及热图（显示t值前10位及最后10位基因）表示，见图2。

图2 Limma差异分析火山图及热图

2.2 WGCNA结果

从2个数据集合集选取基因并用于构建加权基因共表达网络。当软阈值功率设置为16时，尺度独立性达到0.89，平均连接值为0.29。当最小模块大小设置为30时，合并距离＜0.25的模块，最终获得13个共表达模块；值得注意的是grey模块被认为无法被分配给任何模块的基因集合。然后对每个模块与临床特征进行相关性分析。purple 模块与PD 的正相关性最高（r=0.35，P=0.03）。计算与基因的表达相关性获得GS，同时计算模块特征向量与基因的表达相关性获得MM，基于截止标准（|MM|＞0.6及|GS|＞0.1），purple模块中具有高连通性的188个重要基因被提取出来用于进一步分析，见图3。

图3 WGCNA结果

2.3 Hub基因

使用VennDiagram 获得到7 个Hub 基因，使用Veen图表示，见图4；使用箱线图表示各基因在各组的表达，见图5。

图4 Veen图

图5 GSE20314和GSE28894小脑组织中的Hub基因

2.4 诊断模型的构建和验证

使用逻辑回归算法构建基因预测模型，逐步回归分析提示SBNO1、VPS13C 对PD 的诊断影响最显著。在小脑样本数据集GSE20314、GSE28894 中，SBNO1+VPS13C ROC曲线AUC值为0.9386，见图6A。在血液样本数据集GSE18838、GSE6613中进一步验证该模型提示具有良好的诊断性能，AUC 值分别为0.8663、0.6564，见图6B、C。

图6 SBNO1+VPS13C在各组相应的ROC曲线和AUC值

3 讨论

PD 主要为中脑黑质多巴胺能神经元缺失和纹状体多巴胺神经递质减少，导致出现运动迟缓、肌强直等症状，治疗上常以补充多巴胺为主要治疗手段[9]。现有研究表明PD 震颤产生与小脑有较为明显的关系，PD患者腹侧被盖区小脑连接性高于黑质小脑连接性，在正常人群中这种现象则会相反[10]。Kawabata等[11]发现小脑和基底神经节间的功能连接不仅与运动相关，与认知表现也有联系。小脑的功能或形态调节不仅与运动症状相关，与一些非运动症状也有不可忽视的联系，在PD 中小脑发挥着代偿的作用[12]。早期PD 患者中，小脑-丘脑-皮质回路的募集与运动特征的进展呈正相关，同样也参与PD的补偿机制[13]。甚至有学者提出小脑可能参与PD 的猝死机制[14]。Rusholt 等[15]的体视学研究，揭示了PD的白质改变情况。基于现有小脑参与PD 发病机制的学说，相关药物的开发也逐步开展，如姜黄素膳食补充剂可能对PD 的小脑相关症状具有保护作用[16]。越来越多证据表明小脑在PD中的重要性。

人类VPS13基因的突变是造成神经发育和神经退行性疾病的原因，VPS13可能在细胞器之间的脂质交换中起作用。作为VPS13之一的液泡蛋白13C（VPS13C）与PD 有着较为密切的联系，而VPS13A 与舞蹈症棘皮细胞增多症存在联系[17]。通过全外显子组测序已成功识别导致家族性PD 的基因，发现VPS13C 位列其中[18]。已有研究提示VPS13C 变异与PD 相关[19]。PD 和路易体痴呆均为路易体病，VPS13C 与两者都有联系。VPS13C中罕见的错义突变与路易体病相关，隐性复合杂合错义突变可能对VPS13C的表达和功能产生不同的影响。罕见的隐性复合杂合错义突变的组合可降低VPS13C 表达并导致路易体病风险增加[20]。由此可见VPS13C与PD间具有密切关系。在PD细胞模型中，可见VPS13C 部分定位于线粒体的外膜。VPS13C 沉默与较低的线粒体膜电位、线粒体片段化、呼吸频率增加、PINK1/Parkin 依赖的线粒体恶化以及响应线粒体损伤的PARK2转录上调有关。VPS13C功能丧失会增加线粒体对压力的脆弱性，从而激活PINK1/Parkin 依赖性线粒体质量控制途径[21]，或许由此参与PD 的发病过程。不同PD 患者群之间VPS13C存在着些许区别。VPS13C 基因的单核苷酸多态性rs2414739 与等位基因模型的PD 易感性存在显著联系，但这不包括中国人群[22,23]。但在另外一项研究中筛选中国早发性PD患者队列中的VPS13C突变，结果发现PD 患者中仍有携带VPS13C 的复合杂合突变[24]。VPS13C也与东亚人群有着显著联系[25]。在伊朗人中，VPS13C rs2414739 在PD 与正常人群中也存在显著差异[26]。VPS13C在不同人群中是否有着区别，需进一步扩大样本去验证。综上研究可见VPS13C 与PD 间密切的关系，这与我们的研究结果相符。我们的研究结果进一步发现VPS13C 与SBNO1 基因共同构建基于小脑的PD 诊断模型具有良好的诊断性能。在外周血液中，其诊断性能较脑组织有所下降，这或许与脑组织样本往往比血液样本更能代替PD 的病理改变有关。综合本研究结果及现有研究报道，我们推测该模型对PD患者的诊断或许存在一定程度上的指导意义。

鉴于PD 与小脑间存在着相应联系，我们认为SBNO1、VPS13C等Hub基因有可能参与PD的发生、发展过程，未来或可作为PD的潜在干预靶点，且SBNO1+VPS13C 基因构建的诊断模型具有较好的诊断性能。但其中的相关机制及应用价值仍需进一步探讨及研究。