基于生物信息学分析筛选重性抑郁障碍关键基因及构建诊断模型☆

刘丹 张敏 和申 彭代辉 肖泽萍

重性抑郁障碍（major depressive disorder，MDD）是常见的精神疾病，严重影响患者的生命健康[1]，但目前仍无单一机制能够全面解释该疾病的发生[2]。该疾病的诊断主要依靠症状学标准，尚无明确的用于诊断的生物标志物存在。因此，探索MDD 生物标志物至关重要。稳健排序聚合（robust rank aggregation，RRA）算法是一种基于次序统计整合排序列表的算法[3]。由于其对高通量数据的异常值及噪音具有鲁棒稳健性，且不依赖于整合不同数据集的基因表达水平，而是整合各个数据集的排序列表，尤其适用于对不同平台的测序数据进行meta 分析[3]。本研究拟首次使用RRA 算法对公共数据库中MDD 外周血转录组芯片数据进行整合分析，筛选关键基因并建立诊断模型，为探索MDD 的生物标志物提供参考。

1 资料与方法

1.1 数据来源从美国国家生物信息中心（NCBI）GEO数据库（https：//www.ncbi.nlm.nih.gov/geo/）中检索人类MDD 外周血转录组芯片数据集。纳入标准：①研究设计为病例对照研究，病例组患者诊断须为MDD；②患者除MDD外不共患其他疾病；③数据集注释后至少包含15000 个基因。筛选后共得到5 个MDD 外周血转录组芯片数据集，分别为GSE98793[4]、GSE76826[5]、GSE52790[6]、GSE38206[7]、GSE39653[8]。其中GSE98793 包含128 例MDD 患者及64 名健康对照样本；GSE76826 包含10 例MDD患者及12 名健康对照样本；GSE52790 包含10 例MDD 患者及12 名健康对照样本；GSE38206 包含9例MDD 患者及9 名健康对照样本；GSE39653 包含21例MDD患者及24名健康对照样本。

1.2RRA 差异表达基因（DEG）筛选RRA 算法通过对比每一个项目在列表中的实际位置及零假设随机排序的位置，对该项目进行显著性评分，用来重新排序并确定其显著性。当将该算法用于差异表达基因筛选时，其可通过对比每一个基因在每个数据集差异表达基因列表中的实际位置及零假设随机排序的位置，对该基因进行显著性评分，确定其差异表达的统计学意义。使用R4.2.0中limma数据分析包分别对经批次校正后的5 个数据集进行差异表达基因筛选，设置|log2 fold change（FC）|＞0.5和P＜0.05 为筛选阈值，得到每个数据集分别的差异表达基因列表。使用RRA 包对5 个数据集的差异表达基因结果进行整合分析，得出RRA 整合分析后的差异表达基因（P＜0.05）。

1.3 差异表达基因功能与信号通路富集分析使用在线工具Metascape 对差异表达基因进行基因本体（gene ontology，GO）功能信号通路富集分析，随后使用R软件进行可视化分析。检验水准α=0.05。

1.4 关键基因筛选及诊断模型构建Boruta 法筛选关键基因并采用多因素logistic 回归分析构建诊断模型。以样本量最大的GSE98793 为训练集，使用R 中Boruta 包从差异表达基因中进行关键基因筛选。Boruta 法是针对分类问题非常有效的特征选择方法，其依托于随机森林算法，能有效消除不重要的特征[9]。在筛选特征时，Boruta 法根据特征的重要性将特征分为“确认”、“待定”及“拒绝”三类，选取前两类基因作为关键基因，并使用rms 包采用多因素logistic 回归构建诊断模型。使用forestplot包绘制关键基因森林图，分析单个关键基因的表达水平与MDD 的关系。使用roc.test法比较单个关键基因及诊断模型在训练集中的诊断价值。

1.5 诊断模型诊断价值分析采用Bootstrap 重抽样法对GSE98793 重抽样1000 次作为内部验证。研究表明，相比于数据集拆分等其他方法，Bootstrap法在针对logistic 回归模型的内部验证方面具有更好的稳定性和更小的偏倚[10]。将剩余4 个数据集（GSE76826、GSE52790、GSE38206、GSE39653）作为外部验证集，绘制诊断模型区分MDD 患者和健康对照的受试者工作特征（receiver operating characteristic,ROC）曲线并计算曲线下面积（area under curve,AUC），对诊断模型的诊断价值进行外部验证。

本研究中数据分析处理及图像生成均采用R4.2.0。检验水准α=0.05。

2 结果

2.1 差异表达基因分析结果使用RRA 联合分析后共筛选出MDD 与健康对照差异表达基因31 个，其中上调基因20 个，下调基因11 个。差异表达基因热图见图1。

图1 MDD 与健康对照差异基因logFC 热图

2.2 差异表达基因的GO功能信号通路分析GO功能分析结果显示，差异表达基因主要富集的生物过程包括白细胞激活调控、损伤反应、刺激检测、细胞因子生成负调控，细胞组分主要包括突触前、特定颗粒腔及裂解空泡。可视化结果见图2。

图2 差异基因GO 富集分析

2.3 关键基因筛选及诊断模型构建使用GSE98793 作为模型训练集，使用Boruta 法对31 个差异表达基因进行筛选。Boruta 法共筛选出3 个“确认”基因（MMP8、TDRD9、FAM3B），4 个“待定”基因（LCN2、ARG1、NPTN、FANCF）。将上述7 个关键基因进行多因素logistic 回归分析并建立MDD 疾病诊断模型。关键基因森林图表明：FANCF（OR=0.277,95%CI：0.100～0.744）及FAM3B（OR=0.763,95%CI：0.602～0.967）为MDD 的保护因素，其表达水平升高可以降低MDD 的发病风险；TDRD9（OR=2.362,95%CI：1.434～4.071）为MDD 的危险因素，其表达水平升高可增大MDD 的发病风险（图3）。绘制该模型在训练集中的ROC 曲线（图4），AUC 为0.803（95%CI：0.740～0.867）。同时，各关键基因的ROC 曲线及AUC 值如图4。其中MMP8的AUC 值最大，为0.702（95%CI：0.625～0.779）。使用roc.test法比较模型及MMP8的AUC 值，结果显示模型的AUC 值高于MMP8的AUC 值，差异具有统计学意义（Z=2.793,P＜0.001）。

图3 关键基因森林图

2.4 诊断模型内部及外部验证内部验证结果显示AUC 为0.804（95%CI：0.757～0.851）。内部验证的ROC 曲线见图5，模型的校准曲线见图6。同时，该模型在外部验证集中表现出较好的诊断性能。在GSE76826 中的AUC 为0.781（图7A），在GSE38206中的AUC 为0.901（图7B），在GSE39653 中的AUC为0.722（图7C），在GSE52790 中的AUC 为0.725（图7D）。

图5 模型内部验证的ROC 曲线

图6 模型内部验证校正曲线图

图7 模型在各外部验证集中的ROC曲线

3 讨论

本研究使用RRA 的方法整合5 个MDD 外周血转录组数据集，得到31 个差异表达基因。将包含最大样本量共192 例样本的数据集GSE98793 作为训练集，从31个差异表达基因中通过Boruta法筛选出7 个基因作为关键基因建立诊断模型，该模型在训练集及为外部验证集中均表现出具有较好的预测能力。

既往基于转录组数据筛选MDD 生物标志物及构建诊断模型的研究存在一定局限，如HE 等[11]使用GSE98793 数据集建立了4 基因诊断模型，其AUC 为0.799，但其数据来源单一，仅使用了两个转录组数据集，未验证其模型在其他来源数据上的诊断价值。在转录组芯片数据集的meta分析方面，目前仅有1 项研究整合MDD 相关数据集构建了包含70个关键基因的诊断模型，该模型在外部验证集中AUC 为0.78[12]。而本研究构建的诊断模型仅包含7个基因，同时在外部验证集上表现较好，因此在可行性上明显优于上述模型。

对RRA 整合分析后得到的差异表达基因进行GO 功能分析，结果显示差异表达基因主要富集的生物过程包括白细胞激活调控、损伤反应、刺激检测及细胞因子生成负调控，这些通路均与免疫及免疫炎症相关。目前已有大量文献表明，免疫失调及免疫炎症反应与MDD 的病理生理机制及预后相关[13-14]。差异表达基因主要富集的GO 细胞组分包括突触前、特定颗粒腔及裂解空泡，其均与神经递质传递及神经元活动相关。

本研究诊断模型纳入MMP8、TDRD9、FAM3B、LCN2、ARG1、NPTN、FANCF等7 个基因。其中，NPTN基因编码两种亚型的跨膜糖蛋白，Np55 及Np65。Np65 是仅表达于大脑中的细胞黏附分子，主要表达于杏仁核、海马等与记忆及情绪调节相关的脑区中。在分子层面，其能与γ-氨基丁酸（γ-aminobutyric acid，GABA）受体结合，影响GABA能神经元的信息传递[15]。GABA 能系统参与包括抑郁症在内的多种神经精神疾病的发生[16-17]，作为新型抗抑郁药物的可能靶标之一，相关研究表明提升GABA 神经递质水平在抑郁症模型小鼠中具有抗抑郁效果[18]。同时，动物实验表明NPTN基因失活影响5-HT水平，导致小鼠皮质酮水平增高，增加小鼠抑郁样行为[19]。LCN2编码Lipocalin-2 蛋白，面对炎症刺激时诱导中枢神经系统产生趋化因子，参与先天性免疫反应，调节神经元炎症及神经元降解[20]。同时，LCN2还参与调节多种生物行为反应，如认知功能、疼痛敏感性、抑郁和焦虑[20-22]。研究表明，LCN2敲除小鼠出现焦虑及抑郁样行为，另有研究表明出现海马腹侧颗粒细胞和锥体细胞肥大，而这两种细胞的功能与情绪调节相关[22]。一项以老年MDD 患者为研究对象的研究发现老年MDD 患者的血浆Lipocalin-2 蛋白水平较健康对照显著上升[23]，另有研究表明血浆Lipocalin-2 蛋白水平与心衰患者的抑郁症状严重程度相关[24-25]。动物研究表明，LCN2参与调节炎症性肠病小鼠模型的抑郁症状[26]。两项全基因组关联分析研究表明MMP8与MDD 及围产期抑郁障碍的发生相关[27-28]。2021 年的一项研究表明在慢性轻度应激小鼠模型的海马中过表达白细胞介素4 可以上调ARG1表达，进而促进海马神经元生成并降低应激相关抑郁的风险[29]。MMP8、FAM3B及ARG1均参与免疫炎症反应及细胞因子调节，而如前文所述，免疫炎症反应与MDD 的发病机制密切相关[13,30]，因此，这3 个基因也可能通过调控免疫炎症反应介导抑郁的发生。目前尚未有研究探索FANCF或TDRD9与MDD 之间的关系，而由于本研究为基于生物信息学的研究，无法探究其作用机制，故仍需进一步的基础实验探索可能的作用机制。

综上所述，本研究通过对现有MDD 公共数据库数据进行整合分析，筛选出7 个MDD 关键基因，并构建出在内部及外部验证集中均表现出较好诊断能力的诊断模型。同时，7 个关键基因中2 个基因与MDD 的关系目前尚未见研究报告。因此，本研究为进一步探究MDD 的发病机制提供了新的可能方向，为进一步明确其潜在病理机制及基于生物标志物的诊断提供了依据。