帕金森病诊断基因的生物信息学分析

陈忠，连世忠

1 山西医科大学第一临床医学院，太原 030000；2 山西医科大学第一医院神经外科

帕金森病（PD）是与多巴胺能神经元丧失相关的常见的神经退行性疾病之一，其发病率随着年龄的增长而稳步上升［1］。PD 患者运动迟缓并伴有至少一种以上表现（肌肉僵硬、静止性震颤或姿势不稳）［2］。其主要特征是黑质致密部多巴胺能神经元早期变性死亡，细胞内广泛存在α-突触核蛋白聚集［1，3］。迄今为止，PD 的病理机制尚未明确，目前研究 较 多 的 是SNCA、LRRK2、VPS35、EIF4G1、DNAJC13、CHCHD2基因突变，其治疗仍以左旋多巴为首的药物对症治疗［2，4］。2021 年10 月，本研究通过生物信息学分析，对不同数据集中的PD基因信息进行整合，从而找到PD 的关键基因并建立诊断模型，预测靶向调控关键基因的微小RNA（miRNA）和转录因子（TF），为探索PD 的遗传学病因和发病机制提供参考。

1 资料与方法

1.1 资料 从美国国家生物信息中心（NCBI）基因表达综合数据库（GEO，https：//www. ncbi. nlm. nih.gov/geo/）下载3个PD数据集（GSE20146、GSE20153、GSE20141）作为内部训练集，3 个PD 数据集平台注释文件为Affymetrix 人基因表达阵列（GPL570）。GSE20146 数据集包括10 例PD 患者样本和10 例正常样本；GSE20153包括8例PD患者样本、8例正常样本；GSE20141包括10例PD患者样本、8例正常样本。此外从NCBI 基因表达综合数据库另外下载2 个PD数据集（GSE20291、GSE20292）作为外部验证集，2个PD 数据集平台注释文件是Affymetrix 人基因表达阵列（GPL96）。GSE20291数据集包括15例PD 患者的样本、20例正常样本；GSE20292包括11例PD患者的样本、18例正常样本。

1.2 差异表达基因（DEG）筛选 用R4.0.2 软件（https：//www.R-project.org）将3 个PD 数据集合并作为内部训练集，进行预处理（背景校正、归一化、log2 转化）。当多个探针对应1 个共同基因时，取其平均值作为其表达值。用“sva”包消除3 个数据集之间的批次效应。用“limma”包筛选P＜0.05 的DEG，以|log2 Fold change（FC）|＞0.5 作为选择DEG的截止点（FC为倍性变化）。

1.3 加权基因共表达网络分析（WGCNA）及模块鉴定 为了探讨基因间的相互作用，应用系统生物学方法WGCNA 构建基因共表达网络。将样本中有25%以上变异的基因整合的数据集导入WGCNA；为保证网络构建结果的可靠性，剔除离群样本；用“pick-Soft-Threshold”函数，由共同表达式相似度得到的软阈值功率β 计算邻接度；将邻接关系转化为拓扑重叠矩阵（TOM），并计算相应的不相似度（1-TOM）；通过分层聚类和动态树切割函数对模块进行检测；为将表达谱相似的基因分类到基因模块中，对基因树形图进行平均连锁层次聚类，采用“TOMbased”差异测量方法，最小基因组为50；对与临床属性相关的模块，计算模块隶属度（MM，特定基因与模块特征基因之间的相关性）和基因显著性（GS，特定基因与临床变量之间的相关性）。最后，对特征基因网络进行可视化，进一步分析模块中的基因信息［5］。从内部训练集提取的DEG 与WGCNA 重要模块中的基因取交集得到交集基因。

1.4 关键基因筛选及诊断模型构建 逐步回归法选择关键基因并采用多因素逻辑回归分析构建诊断模型。对所有样本中交集基因的表达量取中位数，基因表达量高于中位数则定义为高表达，反之则为低表达。用SPSS23.0 统计软件对数据进行统计分析，用逐步回归筛选出关键基因；用“rms”包采用多因素逻辑回归分析用关键基因构建诊断模型并绘制列线图。

1.5 诊断模型的校准 使用“rms”包绘制校准曲线以评估PD 诊断模型列线图，用“rms”包测量C 指数量化列线图的辨别性能。用R 软件对2 个PD 数据集（GSE20291、GSE20292）进行合并，作为外部验证集，预处理方法同1.2。用外部验证集验证诊断模型，计算C 指数，绘制受试者工作特征（ROC）曲线并计算曲线下面积。

1.6 miRNA-TF-mRNA 调控网络的构建 利用miRTarBase、Starbase 和Targetscan 数据库预测关键基因的靶向miRNA。为了提高预测的准确性，只保留3 个数据库预测的miRNA。使用rich 数据库（http：//amp. pharm. mssm. edu/Enrichr/）预测关键基因的靶向TF。选取P＜0.05 的结果作为截断值。在获得miRNA-TF-mRNA 调控关系后，使用Cytoscape3.7.2 软件对miRNA-TF-mRNA 调控网络进行可视化。

1.7 统计学方法 采用SPSS23.0 统计软件和R4.0.1 软件进行数据处理，用R4.0.1 软件和Cytoscape3.7.2软件进行图像生成处理。P＜0.05为差异有统计学意义。

2 结果

2.1 DEG 筛选结果 DEG 共405个，其中表达上调基因191个，表达下调基因214个。表达上调基因前5 位为KDM6B、GRINA、IFI35、TRIM36、NME4，表达下 调 基 因 前5 位 为BASP1、NCKAP1、CNTNAP2、ZNF536、SYT1。

2.2 WGCNA 分析和模块鉴定结果 所有基因按照方差由大到小进行排序，选取了方差前25%（3 100 个）的基因进行分析。用“flashClust”工具包进行聚类分析，将阈值设置为100，检测并删除3 个离群样本，保留了51 个样本。用“WGCNA”包中的“pickSoftThreshold”函数筛选出1~20 的功率参数。选择β=12 的幂作为软阈值，以保证网络的无标度。将阈值设置为0.20，合并集群树中的类似模块，共获得8 个模块，其中基因具有相似的共表达性状。tan 模块的模块特征基因（ME）与PD 的正相关性最高（r=0.34，P＜0.05），red 模块的ME 与PD 的负相关性最高（r=-0.31，P＜0.05）。因此，包含100 个基因的tan 模块和包含180 个基因的red 模块被鉴定为PD的重要模块。在tan模块中获得基因的MM与GS呈正相关（cor=0.28，P＜0.05）；在red 模块中获得基因的MM 与GS 呈正相关（cor=0.24，P＜0.05）。从内部训练集筛选的DEG 分别与tan 模块和red 模块中的基因取交集得到19个交集基因。

2.3 关键基因及诊断模型 用逐步回归法从19 个交集基因中最终筛选出5 个关键基因，包括PNMA3、AEBP1、PABPC1、GCA、GSTM2。将5 个关键基因进行多因素逻辑回归分析，见表1，建立包含上述5 个关键基因的诊断模型，并将其表示为列线图。该诊断模型的校正曲线显示一致性良好，计算C 指数为0.890，并进一步绘制了ROC 曲线，曲线下面积为0.890。

表1 关键基因的多因素逻辑回归分析

2.4 miRNA-TF-mRNA 调控网络分析结果 得到24 个miRNA-TF-mRNA 调控关系。该调控网络由Cytoscape3.7.2软件构建，包括8个miRNA、9个TF、5 个mRNA，见 图1。8 个miRNA 分 别 为hsa-miR-532-3p、hsa-miR-423-3p、hsa-miR-485-5p、hsa-miR-149-5p、hsa-miR-520g-3p、hsa-miR-429、hsa-miR-200c-3p、hsa-miR-145-5p；9 个TF 分 别 为SNAI2、ATF4、PLAU、SP3、JUN、NFIA、LTF、TCF3、SNAI1；5个mRNA 分别为PNMA3、AEBP1、PABPC1、GCA 和GSTM2。

图1 miRNA-TF-mRNA调控网络

2.5 诊断模型验证结果 外部验证集验证诊断模型，计算C 指数为0.752，并绘制了ROC 曲线，曲线下面积为0.752。

3 讨论

PD 是与多巴胺能神经元丧失相关的最常见神经退行性疾病之一，其发生机制尚未完全明确。在本研究中，我们整合了3 个PD 数据集作为内部训练集，使用两种不同的方法（WGCNA 分析和差异表达分析）得到交集基因。通过差异表达分析，我们在PD 样本和正常对照组之间共鉴定出405 个DEG，其中表达上调基因191 个，表达下调基因214 个。前5位表达上调基因为KDM6B、GRINA、IFI35、TRIM36、NME4，而前5 位表达下调基因为BASP1、NCKAP1、CNTNAP2、ZNF536、SYT1。有文献报道，CNTNAP2可作为参与特发性和携带LRRK2 基因G2019S 突变的PD 患者的潜在候选基因［6］。也有研究发现，在PD 小鼠模型中miR-34-5p 可以靶向调控SYT1 从而发挥神经保护作用［7］。本研究WGCNA 分析发现，tan模块的ME与PD的正相关性最高，red模块的ME与PD 的负相关性最高，因此包含100 个基因的tan模块和包含180 个基因的red 模块被鉴定为PD 的重要模块。从内部训练集提取的DEG与WGCNA重要模块中的基因取交集得到了19个交集基因。

通过对54 个样本中19 个交集基因的表达量取中位数，则基因表达量高于中位数则定义为高表达，反之则为低表达。用逐步回归选择关键基因，筛选出5 个关键基因，包括PNMA3、AEBP1、PABPC1、GCA 和GSTM2，并采用多因素逻辑回归分析构建诊断模型，该诊断模型的校正曲线显示了良好的一致性，计算C 指数为0.890，并进一步绘制了ROC 曲线，曲线下面积为0.890，说明了模型的良好分辨力。外部数据集也再一次验证了其良好的诊断能力，计算C 指数为0.752，曲线下面积为0.752。然而暂无文献证实这5 个基因对PD 的诊断有重要意义。但有研究表明，AEBP1 为胶质瘤的潜在致癌驱动基因，对其治疗和干预具有潜在影响［8］。也有文献报道，PABPC1 过表达可抑制胶质母细胞瘤的恶性进展［9］，而且GCA 对于检测铜绿假单胞菌和其他假单胞菌属物种是否存在有重要意义［10］。研究还表明，GSTM2 是胰腺癌化疗优化和预后生物标志物的潜在靶标［11］。这5 个基因对于PD 的诊断价值需要进一步验证。

通过对关键基因进行miRNA-TF-mRNA 调控网络分析得到24 个miRNA-TF-mRNA 调控关系，包括8 个miRNA 和9 个TF。研究发现，抑制miR-429 可以减弱小胶质细胞的炎症反应和创伤性脑损伤介导的脑损伤［12］。LTF 在脑出血后血肿解毒中也有重要作用［13］。NFIA 则被发现是一种胶质的生成开关，能够从多能干细胞中快速衍生出功能性星形胶质细胞［14］。也有研究表明，MALAT1与miR-200c-3p结合可以上调SIRT1 表达，从而诱导脑微血管内皮细胞自噬并保护脑微血管内皮细胞免受氧和葡萄糖侵害［15］。抑制miR-145-5p 表达可减小大鼠急性脑缺血中的梗死体积［16］。并且ATF4 的表达可加深原发性脑肿瘤的恶性程度，促进肿瘤细胞增殖和肿瘤血管生成［17］。然而以上miRNA 和TF 如何对关键基因进行调控，则需进一步验证。

综上所述，通过综合生物信息学分析得到5 个关 键 基 因，即PNMA3、AEBP1、PABPC1、GCA 和GSTM2，进而建立了诊断PD 的模型，并且证实该模型具有良好的诊断能力。通过对关键基因进行miRNA-TF-mRNA 调控网络分析得到24 个miRNATF-mRNA 调控关系，包括8 个miRNA 和9 个TF。然而关键基因对于PD 的诊断及miRNA 和TF 如何对关键基因进行调控，进而影响PD 的进展，则需进一步研究。以上关键基因、miRNA 和TF 也许会成为PD 的生物标志物，同时也为今后探明PD 的发病机制提供了更多方向。