基于生信分析对阿尔茨海默病炎症相关关键基因的筛选及鉴定
2023-09-28司蕊豪刘羽茜朱仲康王旭侯文晓赵丹玉
司蕊豪,刘羽茜,朱仲康,王旭,侯文晓,赵丹玉
(辽宁中医药大学中西医结合学院1.生物化学教研室;2.生理学教研室;3.组织与胚胎学教研室,沈阳 110847)
阿尔茨海默病(Alzheimer′s disease,AD)是一种常见的慢性神经退行性疾病。据世界卫生组织报道,截至2020 年9 月,全球约有5 000 万AD 患者[1]。由于人口老龄化的加剧,预计到2050 年AD 病例数将增至1.52 亿[2]。AD 的主要临床表现为记忆力减退和认知功能障碍,主要病理特征为β 淀粉样蛋白(amyloid β-protein,Aβ)沉积、tau 蛋白过度磷酸化、神经纤维缠结等[3]。大量研究表明,在AD 发病过程中,各种炎症及炎症相关细胞因子在不同的大脑区域中异常表达,并且与疾病进展显著相关[4-5]。因此,研究炎症相关基因(inflammation-related genes,IRGs)的表达对于诊断与防治AD 尤为重要。本文旨在通过生信分析鉴定AD 过程中与IRGs 相关的特征基因,分析特征基因所富集的生物学途径,筛选并鉴定关键基因,为AD 的临床诊断和治疗提供潜在的靶点。
1 对象与方法
1.1 数据的获取及差异表达IRGs 的筛选 通过GEO 数据库(https://www.ncbi.nlm.nih.gov/geo/)搜索“Alzheimer”并且下载GSE144459[6-7](GPL21163 Agilent-074809SurePrint G3 Mouse GE v2 8x60K Microarray 数据更新至Aug 25,2020)和GSE135999[8](GPL21810 Agilent -074809 SurePrint G3 Mouse GE v2 8x60K Microarray[Feature Number version]数据更新至Aug 28,2021)中的AD 模型小鼠海马和正常小鼠海马转录组数据。使用limma[9]包分别分析两个数据集中的差异表达基因(differentially expre ssed genes,DEGs)。筛选标准为:P<0.05。并 使 用ggplot2[10]对DEGs 进行可视化。然后在NCBI 数据库中的Gene 子数据库中搜索关键词“Inflammatory”及“Musmusculus”获得IRGs。运用在线韦恩图工具(http://bioinformatics.psb.ugent.be/webtools/Venn/) 对上述获得的差异表达基因(DEG1、DEG2)与NCBI 数据库获得IRGs 取交集并获得差异表达IRGs。
1.2 差异表达IRGs 的功能富集分析 基于GO 及KEGG 数据库,使用clusterProfiler[11]包对差异表达IRGs 进行功能富集分析,并对富集结果进行可视化。GO 富集分析包含生物过程、分子功能、细胞组成。本研究中筛选条件为P<0.05,count≥1。选取排名前15 的信号通路探索差异表达IRGs 的生物学功能。
1.3 PPI 网络构建及候选基因筛选
1.3.1 差异表达IRGs 的PPI 网络 为探究差异表达IRGs 之间是否存在互作关系,利用STRING 数据库(https://string-db.org)对106 个差异表达IRGs 进行PPI 网络构建。置信度为0.4(Confidence=0.4),去除离散的蛋白。用Cytoscape 软件对蛋白网络图进行美化,得到PPI 网络图。
1.3.2 PPI 网络中候选基因筛选-cytoHubba 利用cytoHubba 插件来筛选PPI 中的关键节点。根据nodes在网络中的属性进行排名。以MCC、DMNC、MNC、DEGREE、EPC 5 种算法分别进行排序,然后对5 种算法的TOP30 基因取交集,获得候选基因。
1.4 特征基因的筛选及鉴定
1.4.1 LASSO 回归分析筛选特征基因 本研究以GSE144459(Healthy=16,AD=16)作为训练集,以GSE135999(Healthy=6,AD=7)作为验证集,来进行LASSO 回归分析。设置参数为:family ="binomial",type.measure="class";LASSO 分 析 工 具 为R 语 言glmnet[12]包,以lambda 最小值取0.000 时的模型为最佳诊断模型。
1.4.2 鉴定特征基因并绘制受试者工作特征(ROC)曲线 将筛选得到的特征基因输入文件进行主成份分析(principle component analysis,PCA),以确定对疾病及正常小鼠的区分能力。使用pROC[13]包绘制每个特征基因的ROC 曲线,判断特征基因在AD 中的预测诊断能力。ROC 曲线下面积(AUC)的值一般在0.5~1.0 之间。在0.5~0.7 时预测诊断准确性较低,在0.7~0.9 时具有一定的预测诊断意义,在0.9 以上越接近于1.0,说明其预测诊断效果越好。
1.5 LASSO、SVM 、GSVA 分类器构建 分别以GSE144459 及GSE135999 作为训练集和验证集,采用LASSO、SVM、GSVA 3 种不同的算法分别计算8个特征基因,使用R 语言pROC 包绘制训练集及验证集的ROC 曲线,并计算AUC,一般AUC 值越大,说明预测的结果较为准确,对应的分类器具有较高的分类意义。
1.6 特征基因功能网络 为进一步探索特征基因(adj.P<0.05)所靶向的通路及其发挥的功能,本研究使用clusterProfiler 包对其进行基于GO 和KEGG 的富集分析。并绘制特征基因与通路的网络图并且对富集结果进行可视化。
1.7 特征基因的表达验证 用AD 相关的数据集GSE122063(AD=56,Control=44)验证8 个特征基因(Fcgr2b、Cd14、Fcgr4、Fcgr1、Lyz2、Fcgr3、Ly86 及Themis2)。采用秩和检验和R 语言ggplot2 包对特征基因的表达进行分析和可视化。
2 结果
2.1 差异表达 IRGs 的筛选 用limma 包筛选GSE135999 数据集AD/WT 小鼠海马组织中共存在2 101 个差异表达基因(DEG1),其中920 个上调基因,1 181 个下调基因;GSE144459 数据集AD/WT小鼠海马组织中共存在9 239 个差异表达的基因(DEG2),其中3 883 个上调基因,5 256 个下调基因。如图1A、1B 所示。对上述获得的差异表达基因(DEG1、DEG2)与NCBI 数据库获得IRGs 取交集,以获得差异表达的IRGs。结果共获得90 个上调表达IRGs,16 个下调表达IRGs ,如图1C、D 所示。
图1 差异表达的IRGs 的筛选Fig 1 Screening of differentially expressed IRGs
2.2 功能富集分析结果 为探索差异表达IRGs 的生物学功能,本研究对106 个差异表达IRGs 进行KEGG 功能富集分析后,使用气泡图对TOP15 通路进行展示(图2A)。得到差异表达IRGs 参与了金黄色葡萄球菌感染、中性粒细胞外陷阱形成、破骨细胞分化、吞噬体、Toll 样受体、FcγR 介导吞噬作用、趋化因子等信号通路。此外对106 个差异表达IRGs 进行GO 分析,发现在BP 中这些基因参与细胞对生物刺激的反应、单核细胞增殖、淋巴细胞增殖、白细胞迁移、炎症反应的调节等方面。在CC 中与溶酶体、吞噬囊泡、膜微区、膜筏、炎症体复合体等显著相关。在MF 中与G 蛋白耦联受体结合、免疫球蛋白结合、Toll 样受体结合、CCR 趋化因子受体结合、细胞因子受体结合等生物学功能相关(图2B)。
图2 差异表达IRGs 功能富集分析Fig 2 Functional enrichment analysis of differentially expressed IRGs
2.3 PPI 网络构建及候选基因筛选
2.3.1 差异表达IRGs 的PPI 网络 利用STRING数据库(https://string-db.org)对106 个差异表达IRGs进行PPI 网络构建,得到98 个蛋白的互作网关系,包括98 个节点,815 条边。然后用Cytoscape 插件对蛋白网络图进行美化,如图3。图中绿色表示下调表达,红色表示上调表达,颜色越深,上调/下调的越明显。
图3 差异表达IRGs的PPI网络Fig 3 PPI networks for differential expression of IRGs
2.3.2 PPI 网络中候选基因筛选-cytoHubba 基于cytoHubba 插件中的算法,对5 种算法的TOP30 基因取交集,获得17 个候选基因。分别是Tyrobp、Fcgr3、Fcgr2b、Fcgr1、C1qa、Cd14、Cd86、Ctss、Irf8、Cd68、Fcgr4、Ccl4、Lyz2、Ly86、Clec7a、Ccl3、Themis2。然后使用UpSetR[14]包进行可视化,得到UpSet 图,如图4。
图4 MCC、DMNC、MNC、DEGREE、EPCTOP30 基因交集Upset 图Fig 4 MCC,DMNC,MNC,DEGREE,EPC TOP30 gene intersection Upset diagram
2.4 特征基因的筛选及鉴定
2.4.1 特征基因的筛选 用LASSO 回归对训练集GSE144459(Healthy=16,AD=16)和验证集GSE135999(Healthy=6,AD=7)进行分析,得到基因系数的图形和交叉验证的误差图。如图5,以lambda 最小值取0.000 时的模型为最佳诊断模型,可见LASSO 回归分析筛选得到由8 个基因构成的诊断模型,分别为Fcgr2b、Cd14、Fcgr4、Fcgr1、Lyz2、Fcgr3、Ly86 及Themis2。
图5 LASSO 回归分析筛选特征基因Fig 5 Selection of characteristic genes by LASSO regression analysis
2.4.2 鉴定特征基因并绘制ROC 曲线 对8 个特征基因进行PCA 分析得出无论是在GSE135999还是GSE144459 中,特征基因均可较好的区分AD组和正常组(图6)。此外,通过绘制特征基因的ROC 曲线(图7) 可知每个基因的AUC 值均在0.7以上,由此说明特征基因对AD 具有较强的预测能力。
图6 特征基因对样本的PCAFig 6 PCA of characteristic gene pair samples
图7 特征基因的鉴定Fig 7 Identification of characteristic genes
2.5 LASSO、SVM、GSVA 分类器构建 由LASSO模型(图8A)中可以看出,训练集及验证集的AUC均大于等于0.8;SVM 模型中(图8B)训练集及验证集的AUC 值均为1;GSVA 模型中(图8C)中训练集及验证集的AUC 值均在0.9 以上。由此说明3 种模型均为较好的分类器,对AD 具有较高的预测能力,因此这8 个基因可作为AD 的特征基因。
图8 LASSO、SVM、GSVA分类器的ROC 曲线Fig 8 ROC curve of LASSO,SVM,GSVA classifiers
2.6 特征基因功能网络 在生物过程方面,共获得151 个Terms,特征基因与抗原刺激的急性炎症反应、B 细胞介导免疫的调节、免疫球蛋白介导免疫应答的调节等功能相关(图9A)。在分子功能方面,共获得9 个Terms,主要与酰胺结合、β-淀粉样蛋白结合、水解酶活性、免疫受体活性、溶菌酶活性等功能相关(图9B)。在细胞组成方面,共获得7 个Terms,与高尔基体顺面内胞浆网、质膜外侧固有成分、膜微区、膜筏等功能相关(图9C)。此外特征基因参与了急性髓系白血病、FcγR-介导吞噬作用、利什曼病、中性粒细胞外陷阱形成、系统性红斑狼疮、肺结核等疾病相关通路(图9D)。
图9 特征基因功能网络Fig 9 Characteristic gene function network
2.7 特征基因的表达验证 如图10 所示,8 个特征基因中的5 个基因(Fcgr2b、Cd14、Fcgr1、Fcgr3、Ly86)在AD 脑组织中表达上调且有统计学意义。因此可作为AD 发展过程中与炎症相关的关键基因。
图10 5 个特征基因在相关的GSE122063 数据集的脑组织样本中的表达Fig 10 Expression of five characteristic genes in brain tissue samples from the associated GSE122063 dataset
3 讨论
AD 是最常见且发病率最高的痴呆类型,严重影响着人们的生活质量,已经被世界卫生组织列为全球公共卫生重点之一,同时也是造成老年人死亡的第三主要病因[2]。目前尚未找到有效治愈的药物或方法。临床上被诊断具有轻度认知障碍的AD 患者其病情往往已经进展到了较晚的阶段,此时发生了不可逆转的脑损伤,不能通过保护神经或清除Aβ沉积来治愈,这也是目前AD 治疗效果不理想的重要原因[15]。因此,早期诊断以及寻找有效的治疗靶点对AD 的防治具有重要的意义。大量研究表明,神经炎症在AD 发生、发展中起着重要作用。本研究通过生物信息学分析方法鉴定出Fcgr2b、Cd14、Fcgr1、Fcgr3、Ly86 为AD 神经炎症过程中的关键基因,为AD 的机制研究及诊断治疗提供新的思路。
神经炎症在AD 发生中涉及Aβ 的生成、神经元丢失、氧化应激等多方面病理过程。作为脑内固有免疫细胞的小胶质细胞在神经炎症的发生、发展过程中起着至关重要的作用。AD 发病早期,在炎症因子和Aβ 等病理产物的作用下,小胶质细胞发挥趋化和吞噬作用对病理物质进行清除,并减轻炎症。然而随着疾病的逐渐加重,小胶质细胞被持续刺激导致过度激活,其趋化和吞噬能力严重受损,对病理产物的清除能力下降,使得Aβ 大量积聚,此时促炎因子白细胞介素-1β、白细胞介素-6、肿瘤坏死因子-α 等产生持续增多,抑炎因子分泌减少,神经炎症不断加重,产生神经毒性,进一步加重神经损伤[16]。AD 发病初期神经炎症的原因目前尚未完全阐明。在本研究中,KEGG 富集分析结果提示AD 神经炎症的发生与金黄色葡萄球菌感染、中性粒细胞外陷阱形成、细胞因子与细胞因子受体相互作用、FcγR 介导吞噬作用、Toll 样受体等信号通路相关。GO 富集分析结果提示AD 的发生与白细胞迁移、炎症的调节、白细胞增殖、对细菌源性分子的反应、对外部刺激反应的积极调节、细胞因子产生的正调节等生物学功能相关。因此研究神经炎症机制对AD 的诊断及防治具有至关重要的意义。
研究表明,Fc 片段结合受体(fcreceptors,FcRs)能够与IgG 的Fc 部分相结合。根据其不同的亲和力,FcRs 可分为激活型受体和抑制型受体。其中Fcgr1、Fcgr3 为激活型受体,Fcgr2b 是唯一的抑制型受体[17]。激活型受体在与免疫复合物(immune complex,IC)交联后,Src 家族激酶(src family kinase,SFK)使免疫受体酪氨酸激活基序(immunoreceptor tyrosine-based activation motif,ITAM)上的酪氨酸残基磷酸化,募集脾酪氨酸激酶(spleen tyrosine kinase,SYK),介导信号的激活从而调控下游分子表达[18]。Fcgr1、Fcgr3 广泛存在于人和小鼠的小胶质细胞中。AD 状态下,两者在小胶质细胞中表达量增加。小胶质细胞在激活型受体的作用下,对Aβ 的吞噬能力增强[19。Fcgr2b 虽然不含ITAM 基序,但其含有免疫受体酪氨酸抑制基序(Immunoreceptor tyrosinebased inhibition motif,ITIM),能够通过其携带的ITIM 基序来拮抗激活型受体所引起的促进吞噬、释放炎性介质等信号,转导抑制型信号[20]。有研究提出Fcgr2b 可能通过调控吞噬作用来参与心肌缺血再灌注损伤的病理过程[21]。此外,Fcgr2b 是寡聚体Aβ的受体,寡聚体Aβ 会导致在原代培养的小胶质细胞中Fcgr2b 聚集增多[22]。Fcgr2b 过表达则会介导Aβ对神经细胞的损伤作用,导致小胶质细胞凋亡[19]。总之,Fcgr1、Fcgr3、Fcgr2b 在神经炎症AD 中的具体机制仍需要进一步研究。本研究通过对GSE135999 和GSE144459 数据集差异表达的IRGs 分析表明,在AD 小鼠海马中,Fcgr1、Fcgr3、Fcgr2b 均为上调基因,在外部数据集GSE122063 中人脑组织中得以验证,并且有较高的诊断效能,因此可能成为AD 发病过程中神经炎症相关的关键基因。
Cd14 作为白细胞分化抗原的一种,主要分布在单核细胞、巨噬细胞等细胞表面[23],其生物学活性主要表现在识别并结合脂多糖复合物。Cd14 作为脂多糖受体,可以参与并调控脂多糖性细胞反应[24]。研究报道,多种损伤相关分子模式(damage-associated molecular patterns,DAMPs)可以作用于Cd14,促进炎性细胞分泌促炎因子[22]。Cd14 还作用于RNA 剪接、转运,参与蛋白酶体蛋白质的分解过程。Cd14 也被证实与AD 的发病相关,其作用机制可能是小胶质细胞表面激活的Cd14 促进细胞摄取Aβ,同时调节免疫应答,改变脑内的炎症,影响AD 的发生、发展[25]。本研究结果表明Cd14 在AD 小鼠的海马以及人脑组织中均表现为上调基因,与之前学者研究一致[25]。
Ly86 基因所编码的蛋白质为MD-1,它与RP105(Toll 样受体家族蛋白)相结合形成免疫复合物RP105/MD-1,在B 细胞、巨噬细胞和树突细胞等免疫细胞中表达[26]。有研究猜测Ly86 可能参与免疫功能异常相关疾病的调节过程[27]。已证实Ly86 基因的DNA 甲基化在肥胖、胰岛素抵抗和炎症中发挥重要作用[28]。有研究发现,与野生型相比,MD-1 基因敲除可有效缓解高脂饮食诱导的小鼠脂肪组织炎症及胰岛素抵抗。感染等内质网应激触发因素可诱导血清和尿液中MD-1 的产生[29]。本研究结果发现Ly86 在AD 小鼠的海马和人脑组织中均为上调基因,其机制还有待研究。
综上所述,本研究基于生物信息学分析筛选出与AD 炎症相关的8 个特征基因,其中5 个关键基因(Fcgr2b、Cd14、Fcgr1、Fcgr3、Ly86)可能成为AD潜在的诊断和治疗靶点,其作用机制有待进一步探究。本文为探索AD 神经炎症发生、发展的分子机制、诊断和治疗提供了依据。