亚洲人阿尔兹海默症miRNA-mRNA网络的生物信息学分析

杨泽若，张 燚，胡 柳，温 轶

(浙江养生堂天然药物研究所， 杭州 310024)

阿尔兹海默症(Alzheimer’s disease，AD)是一种慢性神经退行性疾病，在临床上主要表现为记忆丧失与日常行动功能障碍。据阿尔兹海默协会(Alzheimer’s Association)报道[1]，截至2018年， 阿尔兹海默症是造成美国65岁以上人群死亡的第五大原因。根据中国疾病预防控制中心发布的数据显示[2]，作为一种高致残率、影响独立生活能力的疾病，AD所导致的死亡人数在中国所有疾病中排名第五。AD的病因众说纷纭，但确切的发病机制尚未阐明，迄今也无特效治疗或逆转疾病进展的药物。

随着高通量测序技术和全基因组关联分析研究的发展，多个与AD患病相关的基因被发现，其中已证实有三个基因——淀粉样前体蛋白(APP)、早老素1(PSEN1)和早老素2(PSEN2)是家族性阿尔兹海默症(Familial Alzheimer disease，FAD)的致病基因。目前针对AD的研究主要集中在高加索人群，而对亚洲人群以及中国人群的研究相对较少。Jia等[3]招募了404个家系的1 330例AD或轻度认知障碍(Mild cognitive impairment，MCI)患者并检测PSENs / APP突变，鉴定了11个新的突变位点，新的PSENs / APP突变表明中国人与其他种族之间AD发病机制可能存在异质性。Zheng等[4]对汉族FAD患者进行了全基因组测序，鉴定出了汉族人特有的基因外显子错义突变rs3792646，表明Complement C7是汉族人患AD的新风险基因。Han等人指出[5]，CLU rs11136000多态性与白种人AD发病显著相关，而它与亚洲人群AD发病并没有显著性关联。由此我们猜测，不同人种之间的AD遗传/发病机制存在差异，因此有必要对研究相对匮乏的亚洲人群进行研究，以挖掘亚洲人群特有的AD致病基因、生物标志物以及风险因素。

microRNA(miRNA)是一类保守的非编码RNA小分子，它们主要通过抑制靶mRNA翻译或促使其降解实现对靶基因的转录后调控[6]。miRNA在神经发育、组织分化和突触形成过程中起着重要作用[7]。此外，越来越多的证据表明血浆miRNA含量可能作为潜在的疾病标志物，比如血浆miRNA-206水平升高可预测认知能力下降和痴呆程度[8]。目前，有关miRNA 在亚洲人群AD发生过程中的作用研究甚少，而且关于亚洲人AD的表达谱研究往往只针对mRNA 或者miRNA表达谱，并没有将二者进行关联分析，结合miRNA 及mRNA表达谱的亚洲人AD系统分析尚未见报道。

因此，收集现有公开的亚洲人AD与正常人mRNA/miRNA表达谱数据，通过生物信息学手段，筛选核心差异基因与差异miRNAs，构建miRNA-mRNA的差异表达网络。从转录组水平系统的对亚洲人群可能的AD发生机制进行分析，为亚洲人AD的诊疗提供新的潜在靶标。

1 材料与方法

分析流程主要包含以下环节：数据收集与下载，数据集预处理与单个数据集差异分析，可合并脑区的荟萃分析(Meta分析)和合并分析，核心差异基因选取，潜在靶基因预测，核心miRNA-mRNA网络构建，差异基因/靶基因富集分析，靶基因相关蛋白的互作分析(见图1)。

图1 分析流程图Fig.1 Analysis flow chart

1.1 基因表达数据收集

亚洲人AD患者与正常人的表达谱数据来自基因表达综合数据库(Gene Expression Omnibus, GEO)[9](https://www.ncbi.nlm.nih.gov/geo/)。通过限定关键词“Alzheimer AND HAN”或者“Alzheimer AND ASIAN”， 选取了包含额叶 (Frontal Cortex, FC) 、颞叶 (Temporal Cortex，TC)、海马体 (Hippocampus，HP) 与内嗅皮质 (Entorhinal Cortex，EC) 等AD发病相关脑区的数据。截止2020年6月，一共收集4个表达谱数据集，分别是GSE131617[10]、GSE36980[11]、GSE139384[12]和GSE120584[13]。为确保数据质量，剔除了小于50岁的样本以及在主成分分析中明显离群的样本。最终一共获得了212个AD患者脑组织mRNA数据(76个FC，71个TC，58个EC以及7个HP)与92个正常人(Control，CT)脑组织mRNA数据(34个FC，35个TC，13个EC以及10个HP)。此外，miRNA血清数据集中包含1 021个AD以及288个CT。4个数据集的具体信息(见表1)。

表1 从GEO收集的亚洲人AD患者与正常对照的表达谱数据集Table 1 Expression data sets for AD patients and normal controls in Asian populations collected from GEO

1.2 数据集下载、预处理以及单个数据集差异分析

首先用R包GEOquery[14]下载数据集，利用探针注释的文件对表达矩阵进行探针ID转换。接着利用R包limma[15]的线性模型，对单个表达谱数据进行归一化处理、性别校正与差异分析：其中miRNA数据集进行Benjamini-Hochberg(BH)法[16]校正，P校正值 < 0.05为显著差异；mRNA数据集取P值 < 0.05 为显著差异。由于GSE131617数据集根据布拉克分期(Braak Staging)[17]将AD样本分为多组，将该数据集中布拉克分期为0的样本划分为正常人，将其余布拉克分期样本统一合并为AD患者。

1.3 FC/TC脑区多个mRNA数据集整合分析

由于多个数据集中均包含FC和TC两个脑区的数据，同时采用以下两种算法进行整合分析：1)合并同脑区多个数据集，去除批次效应后进行差异分析；2)meta分析中的鲁棒排名聚集算法(Robust Rank Aggregation，RRA)[18]。

1.3.1 FC/TC脑区mRNA数据集的合并分析

将不同研究的芯片数据合并进行归一化，可以移除不同研究的批次效应，同时保留真实的生物学差异[19]。通过R包sva中的ComBat函数[20]，将三个mRNA数据集的表达矩阵进行合并，并且去除批次效应，然后用limma的线性模型对合并之后的表达矩阵进行性别校正，最后对合并之后的表达矩阵进行差异分析。

1.3.2 FC/TC脑区mRNA数据集的RRA分析

对单个数据集中FC或TC脑区得到的差异基因(P< 0.05)，按照绝对值表达差异倍数(|logFC|)从大到小进行排列，并进行RRA运算，以RRAscore< 0.05 作为筛选阈值，获得三个mRNA数据集数据中排名相对靠前的差异基因。

1.4 核心差异基因选取

对于存在于多个数据集中的脑区(FC/TC)，通过1.3中两种方法的交集来获得核心差异基因。对于仅存在于单个数据集的脑区(EC/HP)，通过进一步降低P值与增加绝对值表达差异倍数(|logFC|)的筛选阈值来获得核心差异基因。

1.5 差异miRNA潜在靶基因预测

利用miRWalk在线工具[21](http://zmf.umm.uniheidelberg.de/apps/zmf/mirwalk/)进行差异miRNA的靶基因预测。选择靶基因3端UTR作为miRNA的作用区域，并选择PicTar5[22]、 miRWalk2.0、 TargetScan[22]、DIANA[23]、miRDB[24]以及miRanda[25]等6个靶基因数据库，选取至少被其中2个靶基因数据库收录的靶基因作为潜在靶基因。

1.6 靶基因选取、miRNA-mRNA网络构建、靶基因蛋白互作分析

将1.5得到的潜在靶基因与1.4得到的各脑区核心差异基因进行交集，得到四个脑区的靶基因。将靶基因与对应的差异miRNA以及差异表达倍数等导入Cytoscape软件[26]，构建四个脑区的miRNA-mRNA差异表达网络。用geneMANIA[27](https://genemania.org/)数据库，对靶基因进行蛋白互作分析。

1.7 基因本体论(GO)以及京都基因组百科全书数据库(KEGG)通路富集分析

通过R包clusterprofiler[28]对四个脑区的差异基因与靶基因进行GO以及KEGG通路富集分析。

2 结 果

2.1 AD患者与正常对照组miRNA/mRNA差异表达结果

2.1.1 miRNA差异表达结果

以P校正值< 0.05和|logFC| > 0.5作为差异miRNA的选择标准，共筛选出5个差异miRNAs， 其中4个(hsa-mi-22、hsa-miR-24、hsa-miR-125b-1-3p、hsa-mi-125-3p)在AD组下调，1个(hsa-miR-208a-5p)在AD组上调。

2.1.2 mRNA差异表达结果(差异基因与核心差异基因)

差异基因与核心差异基因的具体筛选要求及差异个数统计数据(见表2)。对于存在于多个数据集中的脑区(FC/TC)，采用方法1.31的合并分析进行差异基因筛选，以上差异基因再与1.32中RRA分析结果取得交集获得核心差异基因；对于仅存在于单个数据集的脑区(EC/HP)，通过常规筛选条件获得差异基因，并通过提高筛选条件来获得核心差异基因。

表2 差异基因与核心差异基因的筛选条件与个数统计Table 2 Thresholds for DEGs and KDEGs selections and summary statistics

如表2所示，额叶、颞叶、海马体与内嗅皮质中分别存在346、662、461与234个差异基因(见图2)，研究发现NPTX2、SLC14A1与GJA1在四个脑区中均为差异基因，且上述基因在AD患者(与正常对照相比)的四个脑区中呈现出相同的上调或下调表达趋势。额叶、颞叶、海马体与内嗅皮质中分别存在15、29、63与20个核心差异基因(见表2)。SLC14A1是唯一一个在四个脑区中都被选入核心差异基因的基因。此外，NPAS4同时出现在FC、TC与EC脑区；GMPR、ITGB4、APLNR与NPTX2同时出现在FC、TC与HP脑区。

2.2 miRNA潜在靶基因预测和miRNA-mRNA网络构建

利用miRWalk等6个数据库，对5个差异miRNAs进行潜在靶基因预测，并与前文所述四个脑区的核心差异基因进行交集，获得靶基因。共有4个差异miRNAs的潜在靶基因与核心差异基因存在交集，分别是hsa-mi-22、hsa-miR-24、hsa-miR-125b-1-3p、hsa-mi-125-3p。4个差异miRNAs与四个脑区靶基因的调控关系见miRNA-mRNA网络(见图3)。

2.3 四个脑区靶基因统计及GeneMANIA的蛋白互作分析

额叶、颞叶、海马体与内嗅皮质中分别存在6、10、25与7个靶基因(见表3)。其中，SLC14A1同时出现在四个脑区的miRNA-mRNA网络中，NPTX2出现在FC、TC与HP脑区的网络中，AQP1与ANTXR1出现在FC与TC脑区的网络中，TXNIP出现在TC与HP脑区的网络中。此外，DDIT3、FOSB为FC脑区特有靶基因；HLA-DMA、TPT1、BRMS1、DLG4、ITPKA为TC脑区特有靶基因；CD44、EMP1、PRRX1、HCN1、KCNH5、HTR3B、CBLN4、RAB15、GFRA2、GABRA3、HPCA、PTPN3、IL12RB2、CACNG3、CTXN3、LRRC2、NCALD、HTR2A、GABRA1、RGS4、PCSK1、MET为HP脑区特有靶基因；FOS、SOD2、PDYN、NAMPT、PLAC8、MT1M为EC脑区特有靶基因。靶基因的GeneMANIA结果表明，四个脑区中的靶基因间均存在复杂的相互作用，包括共表达相关性、共享蛋白结构域、基因相互作用、物理相互作用与共同通路等等。值得注意的是，HP靶基因数目最多、差异倍数更大、具有更复杂的蛋白相互作用关系、且多在HP中特异表达；EC靶基因与其他脑区靶基因交集最少，SLC14A1为唯一交集。

表3 靶基因统计以及geneMANIA蛋白互作分析结果Table 3 Target genes and geneMANIA results of the PPI network analyses

2.4 四个脑区中差异基因/靶基因GO与KEGG通路富集结果

四个脑区差异基因的GO与KEGG富集结果(见图4)。结果表明FC、TC与HP脑区的GO富集结果比较接近，主要集中在神经递质传导、化学突触传递等神经信号传递相关通路，而EC脑区的GO富集结果更偏向于免疫相关，如中性粒细胞脱颗粒等。KEGG分析中， FC与TC脑区都富集到了神经性退行性疾病通路，如亨廷顿氏病、渐冻症、帕金森症等，HP富集到了神经递质传导相关通路，如神经活性配体受体相互作用信号通路、钙离子信号通路、GABA能突触信号、谷氨酸能突触信号，而EC主要富集在炎症、免疫相关通路。

四个脑区靶基因的GO与KEGG富集结果(见图5a,5b)。从GO分析结果来看，四个脑区各自的靶基因可被同时富集到转运相关的生物学过程中；FC、TC与HP的靶基因被同时富集到神经递质受体活性调控的生物学过程；HP的靶基因可被特异性地富集到离子跨膜运输调控过程中；EC的靶基因可被特异性地富集到低温应激、衰老、神经元死亡、活性氧应激等4个生物学过程中。从KEGG分析结果来看，四个脑区的靶基因未能富集到同一通路上，但FC、TC与EC脑区的靶基因能被同时富集到可卡因成瘾通路。

3 讨 论

本研究收集了GEO数据库中亚洲人AD患者与正常人mRNA/miRNA数据，通过荟萃分析与合并分析，从转录组水平系统地对亚洲人群可能的AD发病机制进行了分析，首次构建了亚洲人四个脑区AD的miRNA-mRNA差异表达网络。

目前有多项生信研究利用已有AD表达谱数据对AD机制进行了系统性研究。Xu等[29]运用合并分析方法，整合了20个表达谱数据，构建了AD患者四个脑区的mRNA转录谱，发现了YAP1等AD上游调控因子。Shi等[30]通过分析1个RNA-seq数据集，鉴定了AD患者上游调控的lncRNA，构建了lncRNA-mRNA共表达网络。Moradifard等[31]分析了1个miRNA以及6个mRNA芯片数据，通过荟萃分析发现了AD中重要的mRNA/miRNA相互作用关系。以上工作往往只采用一种整合分析方法，如荟萃分析或合并分析，或者只进行了单组学分析。本研究对多数据集脑区的数据同时应用了两种整合分析方法，可以更加准确地得到多数据集脑区中的差异基因。此外，上述生信研究的表达谱数据均为高加索人群，而不同人种之间的AD遗传/发病机制可能存在差异，因此本研究旨在收集亚洲人AD表达谱数据，构建亚洲人AD的miRNA-mRNA差异表达网络。

发现5个差异miRNAs，已有文献报道miR-22与miR-24在AD病人中低表达，与本研究中的发现一致：Jovicic等[32]指出miR-22在阿尔兹海默症患者中低表达；Hu等[33]通过荟萃分析发现miR-24在AD患者的脑脊液中低表达。尽管Ylmaz等[34]指出hsa-mi-125-3p和hsa-miR-125b-5p在AD患者中(相比于正常人)并没有显著差异，但是其研究样本较少(AD 172个，CT 109个)，且亚洲人群占比较少。因此，这两个差异miRNAs及其靶基因在亚洲人群AD发生发展中的作用值得进一步的研究。此外，本研究发现 NTPX2、GJA1与SLC14A1在四个脑区中均为差异基因。Xiao等[35]表明NPTX2在 AD患者皮质中大量减少，且与患者认知能力和海马体积密切相关。Kajiwara等[36]指出星形胶质细胞特异性表达的GJA1在AD病人大脑上调，并发现GJA1与AD淀粉样蛋白、tau病理以及认知功能密切相关。Kerstin等[37]指出，SLC14A1的表达在AD病人和小鼠模型中均显著上调。上述基因在本研究的表达趋势与文献报道一致，进一步证明了本研究分析策略和结果的可靠性。

SLC14A1(尿素转运蛋白1)是唯一在四个脑区中共有的靶基因。通过蛋白互作分析，发现SLC14A1和AQP1(水通道蛋白1)存在蛋白相互作用，同时AQP1在AD患者中表达上调，且四个脑区靶基因的GO都富集到液体转运等生物学过程。有研究表明SLC14A1在亨廷顿舞蹈病中显著增高，从而破坏大脑中的尿素平衡，导致纹状体神经元大量死亡[38]。此外，Hansmannel等[39]通过AD患者额叶皮层的转录组测序，指出脑中尿素平衡与AD发病相关。因此猜测SLC14A1与AQP1在AD患者中的表达上调与相互作用，可能会引起脑中尿素循环与水循环紊乱，从而造成神经损伤并诱导AD发病。此外，通过艾伦大脑图谱(Allen Brain Map)[40]中的单细胞转录组数据，发现AQP1、SLC14A1都在星形胶质细胞(Astro L1-6 FGFR3)中特异表达。 Misawa等[41]也指出，星形胶质细胞中AQP1上调可能与AD患者脑中Aβ沉积相关。因此本研究进一步推测，AD患者星形胶质细胞SLC14A1、AQP1的差异表达，将破坏星形胶质细胞的水平衡与尿素平衡，从而引起细胞肿胀破裂，最终影响神经元功能。

通过整合四个脑区靶基因和miRNA-mRNA的差异网络构建，对单脑区的特异性靶基因进行了分析，特别是HP脑区，其靶基因数目最多、差异倍数更大、多数在海马中特异表达、并且主要是配体门控性离子通道相关基因，包括GABRA1(A型γ-氨基丁酸受体亚基Alpha1)，GABRA3(A型γ-氨基丁酸受体亚基Beta3)，HTR3B(5-羟色胺受体3B)与HTR2A(5-羟色胺受体2A)等。Agenor等[42]表明GABA受体 α1、β3亚型在AD患者的海马中表达下调。Garcia-Alloza等[43]发现胆碱能与5-羟色胺能系统之间的不平衡，可能与AD的认知障碍相关。此外，Pedro 等[44]证明Cbln4参与形成并维持GABA能连接，敲除Cbln4后神经元GABA能连接大大减少，且Cbln4在AD小鼠的海马区中表达显著下调。上述基因在本研究中的HP脑区中均表达下调，因此，本研究推断亚洲人群AD患者的海马体GABA受体亚型与5-羟色胺能受体亚型表达下调，可能破坏兴奋性能与抑制性能神经元之间信号平衡，最终造成海马体中神经元的功能障碍。

此外，EC脑区靶基因与其他脑区靶基因交集只有SLC14A1，且其GO/KEGG富集通路与其他三个脑区有较大区别。Flynn等[45]表明SOD2(线粒体超氧化物歧化酶2)可以将超氧化物还原为过氧化氢，保护神经元细胞免受谷氨酸诱导的氧化应激和细胞毒性。Xing等[46]指出，NAMPT(烟酰胺磷酸核糖基转移酶)减少导致NAD +减少，从而导致NAD+ / NADH的比率降低与线粒体功能障碍。在本研究中，SOD2、NAMPT等基因在EC中下调，因此本研究猜测亚洲人群AD患者EC中SOD2与NAMPT的下调可能引起线粒体功能障碍，继而造成突触受损、神经细胞凋亡，最终促使AD的病情发展。

本研究主要基于GEO数据库中的微阵列芯片数据，因此具有一些局限性。首先，微阵列芯片数据受制于其技术的局限性，例如技术动态范围较窄、灵敏度较低等。其次，本研究使用的芯片数据来自不同条件下的多个研究，样本量不均并且许多患者临床细节未知(如用药情况等)，增加了本研究不确定性。尽管本研究采用两种整合分析方法，也无法完全消除各项研究的异质性。此外，本研究受限于亚洲人数据集的样本量，某些脑区(如海马体)样本量较少，其结果可能不具备普适性。值得注意的是，上述芯片数据的基因表达量为单脑区所有细胞的平均值，并没有达到单细胞精确度，因此我们无法直接证实SLC14A1是否在AD患者的星形胶质细胞特异差异表达，并影响AD患者星形胶质细胞的功能。最后，本研究没有详细探究FC与TC脑区特异性变化的靶基因，可能忽略了这些靶基因在上述脑区的特异性作用。本研究虽然对SLC14A1在亚洲人AD中的潜在作用进行了推测，但仍然需要构建相应的细胞和动物模型进一步研究其确切的作用机制。此外，本研究希望通过单细胞转录组技术，核实本研究的推测，并进一步研究星形胶质细胞在亚洲人AD中的作用机制。

4 结 论

收集现有公开的亚洲人AD与正常人mRNA/miRNA表达谱数据，通过生物信息学手段，从转录组水平系统地对亚洲人群可能的AD发生机制进行了分析，最终筛选出差异miRNAs与核心差异基因，首次构建了亚洲人四个脑区的关键miRNA-mRNA差异表达网络，发现了SLC14A1基因在四个脑区均差异表达。本研究推测SLC14A1和AQP1蛋白相互作用，影响了AD患者星形胶质细胞的体液转运功能，最终影响神经元细胞，从而诱导AD的发病。此外，研究发现亚洲人AD海马体与内嗅皮质中的特异表达靶基因，分别与神经元间信号传递与线粒体功能障碍相关。以上研究为亚洲人AD的诊疗提供新的潜在靶标，也为组学分析提供了新的思路。