基于数据挖掘的股外侧肌基因表达谱的增龄性变化研究

杨若愚，李晓凡，李馨逸，朱利扬，王碧云

随着我国人口老龄化的不断加剧，老年人的健康状况越来越受到重视。有相关研究表明，老年人肌肉衰减综合症的发病率呈上升趋势，严重影响到了老年人群的生活质量[1]。 肌肉衰减综合症是一种以肌肉质量、力量、功能丧失为特征的老年综合征，可能会导致老年人肌肉萎缩、 肌力下降、 行动能力受限、跌倒、骨折、丧失自主生活能力甚至死亡，严重影响患者的身心健康，并给家庭社会带来沉重负担[2]。

近年来的研究数据显示，肌肉衰减综合症作为一种增龄性疾病，与多种发病因素相关，这些因素主要包括：个体运动量消减、神经肌肉功能衰退、蛋白质的摄取及合成降低、个体激素水平波动、细胞凋亡及微环境变化、骨骼肌线粒体功能异常、活性氧水平增高、 慢性炎症反应、 骨骼肌自噬性程序性细胞死亡、钙的稳定状态失衡、自由基氧化受损、缺乏热量与蛋白质的摄取、 骨骼肌修护功能损坏及基因与人种等[3-4]。 肌肉衰减综合症的高发人群为老年人，多伴有骨质疏松症[5]，且与跌倒、活动障碍及代谢紊乱密切相关[6]，并易造成慢性心肺疾病和吞咽功能障碍，是老年人生理功能逐渐减退的表现之一和重要原因。

肌肉衰减综合症的发病机制目前尚不完全明确，但其与年龄、激素水平、炎症及促细胞分解因子、神经肌肉调节失衡、营养等因素有关[7]。 随着高通量测序技术的不断发展，高通量测序技术在人类疾病诊断、文库筛选、转录调控、植物反转录研究等方面得到了越来越普遍的应用，并取得了巨大的成功，使得通过高通量测序技术研究肌肉衰减综合症成为了一种可能。 通过筛选与肌肉衰减综合症可能相关的差异表达基因，为之后针对肌肉衰减综合症进行有效的靶基因治疗，提供肌肉衰减综合症可靠的诊断和预后生物学标志物的研究及诊断模型的建立打下基础。本研究利用生物信息学的方法，从美国国立生物 中 心 （National Center for Biotechnology Information，NCBI）的基因芯片公共数据库（GEO）下载并整理股外侧肌增龄性的转录组高通量数据，筛选差异表达基因，寻找可能与骨骼肌增龄性因素相关的分子标记作为肌肉衰减综合症相关的诊断、治疗、预后的生物学标志物。

1 研究对象和方法

1.1 研究对象

GSE1428 数据集包含的10 名年轻 （19～25 岁）和12 名老年（70～80 岁）男性受试者，均为高加索人种。这些受试者健康状况良好，体检及临床实验室检查正常。在采集样本前的6 个月内，受试者没有服用过处方药，也没有进行过阻抗或耐力训练[8]。

1.2 研究方法

1.2.1 数据获取及标准化处理

利用GEO 下载与肌肉衰减综合症相关的基因表达数据集GSE1428，用于筛选差异表达基因，该组数据采用Affymetrix Human Genome U133A Array（GPL96）平台进行检测。利用R 语言平台对GSE1428芯片数据集进行数据评估和标准化处理。

1.2.2 差异基因的筛选

利用R 语言平台读取探针平台文件GPL96，安装并加载limma 包进行差异基因的分析，安装并加载BiomaRt 包对基因ID 进行转换和注释。 分析年轻男性股外侧肌样本和老年男性股外侧肌样本获得差异基因，包括表达上调基因和表达下调基因，显著差异基因的筛选条件设定为adj.P＜0.05 且|log FC|＞0.6（adj.P 为调整后的P，FC 为差异倍数，指老年人基因表达与青年人基因表达的比值），并利用R 语言软件包绘制相关图表。

1.2.3 差异基因的富集分析

利用R 语言平台和Cytoscape 软件进行GO 和KEGG pathway 富集分析并绘制相关图表。 其中，GO富集分析包括生物过程 （Biological Process，BP）、分子 功 能 （Molecularfunction，MF） 和 细 胞 组 分（Cell Components，CC），adj.P＜0.05 为具有统计学意义。

1.2.4 蛋白互作网络制作及核心基因的确定

利用STRING 在线工具和Cytoscape 软件制作蛋白互作网络，并筛选相关核心基因。

2 研究结果

2.1 股外侧肌基因表达差异分析

经过对GSE1428 芯片数据进行基因表达差异分析及差异基因的筛选，共获得差异表达基因2 028个，其中有2 020 个表达下调基因和8 个表达上调基因，根据adj.P 升序排序，将其中排序前25 个差异表达的基因列入表1。

表1 差异表达基因的筛选Table1 Screening of Differentially Expressed Genes

在筛选出差异基因后，对数据进行可视化分析。利用R 语言平台制作差异表达基因的热图（图1）和火山图（图2）。 如图1 所示，热图中的每个小方格表示每个基因，其颜色表示该基因表达量大小，颜色越深，其表达量就越大（红色为上调，绿色为下调）。 每行代表不同样本中每个基因的表达量情况，每列表示每个样品中所有基因的表达量情况。 可从图2 直观地看到上调基因和下调基因差异倍数与adj.P 的分布情况。

2.2 基因表达差异的富集分析

图1 差异基因热图Figure1 Heat Map of Differential Genes

2.2.1 GO 富集分析

在进行差异基因的筛选之后，利用R 语言平台和Cytoscape 以及Cytoscape 的插件BinGO 进行GO富集分析，GO 分析具体过程包括生物过程、分子功能和细胞组分。 图3 为GO 功能富集分析的结果。

图2 差异基因火山图Figure2 Volcano Map of Differential Genes

如图3 所示，差异表达基因主要参与信号转导（signal transduction）、RNA 聚 合 酶Ⅱ启 动 子 转 录 的正调控 （positive regulation of transcription from RNA polymerase II promoter）、炎 症 反 应（inflammatory response）等生物过程。 差异表达基因主要参与的分子功能方面有转录激活子活性、RNA 聚合酶II 核心启动子近端区序列特异性结合（transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding）、蛋白质异源二聚化活性（protein heterodimerization activity）、受 体 结 合（receptor binding）等。 而在细胞组分方面，差异基因主要与细胞质膜（plasma membrane）、质膜的组成成分（integral component of plasma membrane）、 细胞外基质（extracellular exosome）等方面有关。

图3 GO 功能富集分析气泡图Figure3 Bubble Diagram of GO Functional Enrichment Analysis

对其KEGG 通路富集分析，得到了KEGG 通路相关数据信息。 图4 为KEGG 通路富集分析的气泡图。其中富集差异基因较多的通路有代谢信号转导通路（metabolic pathways）、 神经活性配体受体相互作用通路（neuroactive ligand-receptor interaction）、抗菌素生物合成通路（biosynthesis antibotics）和细胞因子-细胞因子受体相互作用通路（cytokine-cytokine receptor interaction）等。

图4 KEGG 通路富集分析的气泡图Figure4 Bubble Diagram of KEGG Pathway Enrichment Analysis

2.2.2 蛋白互作网络及核心基因的筛选

基于差异基因的筛选结果和KEGG 通路相关数据信息，使用STRING 在线数据库将差异基因导入，将最低互作分值（minimum required interaction score）设置成最高可信（highest confidence 0.9），删除与其他蛋白质无相互作用的节点后，完成肌肉衰减综合症的蛋白互作网络，分析差异基因所编码的蛋白质之间的相互作用，得到蛋白互作网络图（图5）。

图5 差异基因的蛋白互作网络图Figure5 Map of Protein Interaction Networks of Differential Genes

将图5 所显示的蛋白互作网路的数据导入Cytoscape 软件进行分析，计算节点的度（degree），利用蛋白互作网络核心基因互相作用的具体算法得出基因在互作网络中的权重值，数值越大说明该基因在网络中的作用越大，并列出degree 排名前10 位的中心节点蛋白，然后与差异基因的KEGG 富集分析所得到的显著性基因取交集，得到核心基因（hub genes）。图6 为排名前10 的核心基因的相互作用关系。

图6 排名前10 位的核心基因互作网络关系图Figure6 Top 10 Core Gene Interaction Networks

表2 显示了图6 中所涉及的核心基因的具体信息及其对应的degree。 这些基因对于肌肉衰减综合症都是下调基因，说明这10 个基因表达水平的下调可能与肌肉衰减综合症的发生和发展存在着一定的关系。

表2 核心基因的具体信息Table2 Basic Information of Hub Gene

3 分析与讨论

肌肉衰减综合症作为一种增龄性的退行性疾病，其发生与发展有着循序渐进的过程，很多基因在肌肉衰减综合症的发生和进展过程中起到至关重要的作用[9]。 肌肉衰减综合症作为一种在老年人群中常见的疾病，对老年人的日常生活活动能力造成了严重的危害，但却又往往被忽视，所以探究与肌肉衰减综合症有关的差异基因则具有重要的临床意义与研究价值[10-11]。虽然有许多关于肌肉衰减综合症发生机制的研究，但其发生机制仍未完全明确，随着新技术特别是高通量基因测序技术的不断发展，通过从基因层面对肌肉衰减综合症进行进一步的讨论和研究，可对肌肉衰减综合症有更深层次的了解。本研究通过对基因芯片数据再挖掘分析，寻找与肌肉衰减综合症有关的差异表达基因，为之后研究肌肉衰减综合症的机制、诊断、治疗和预后提供参考。

本研究首先通过对芯片数据集GSE1428 进行分析，最终发现有2 028 个差异表达基因，其中有2 020个下调基因和8 个上调基因，随后对这些差异基因进行GO 功能富集分析，发现差异基因主要富集于信号转导、RNA 聚合酶Ⅱ启动子转录的正调控、 炎症反应等生物过程。 在分子功能方面主要富集于转录激活子活性，RNA 聚合酶II 核心启动子近端区序列特异性结合、蛋白质异源二聚化活性、受体结合等。 而在细胞组分方面，差异基因富集于细胞质膜、质膜的组成成分、细胞外基质等。 结果提示这些差异基因可能通过不同的信号通路、生物学过程产生了关键的作用，从而影响到肌肉衰减综合症的发生和发展。

在筛选所得到的差异基因KEGG 富集通路结果中，神经活性配体-受体相互作用这一通路具有最高的显著性，而钙信号通路同时具有较高显著性和较多的基因富集。 肌肉衰减综合症患者患病期间肌肉含量的减少和肌肉生理功能的衰退可能与多种神经递质及其受体相关，而神经活性配体受体相互作用信号通路是质膜上所有与细胞内外信号通路相关的受体和配体的集合[12]。

肌肉衰减综合症患者的部分临床表现如肌肉肌力的减弱和引发肌肉收缩的刺激阈值增大都或许与钙信号通路相关。Ca2+通道的化学本质是载体蛋白，在达到刺激阈值之后，Ca2+可与此载体蛋白相结合从而被转运，进而引发后续生理功能的实现，这种Ca2+和载体蛋白的结合过程类似酶与底物的结合过程，具有一定的专一性。生物体中刺激以神经冲动的形式通过运动神经传递至效应器即可引起肌肉收缩，即冲动由神经肌肉接头传达于肌纤维细胞膜，使肌纤维细胞膜产生一个可传导的动作电位，即兴奋-收缩偶联，从而引发横桥运动，引起肌肉的收缩，收缩后的肌肉必须发生舒张后才具有下一次收缩的条件[13]。 在发生兴奋-收缩偶联过程中，由刺激产生的动作电位沿横管系统进入三联管，此时横管膜去极化使终池大量释放Ca2+，Ca2+通过钙离子通道的转运[14]，才可以产生一系列的肌肉功能运动。除了肌细胞本身，衰老、缺氧、缺血所造成的损伤会引起机体过度释放兴奋性神经递质，从而使突触后膜处于一种持续性去极化状态, 造成大量钙离子内流，引发细胞内一系列钙离子的生化反应，导致神经细胞凋亡[15]。

本研究得到的可能与增龄性骨骼肌肌力和功能衰减相关的排名前10 的核心基因，可为后续寻找肌肉衰减综合症相关分子标记物提供线索，但因为是通过对组学数据进行的生物信息学分析，并没有进行相关的功能验证，所以如果直接进行应用还存在一定问题，如果能在模型动物上进行相关的基因功能验证则较为稳妥。

4 小结

骨骼肌增龄性变化存在大量的基因表达的变化，其中绝大多数为基因表达的下调；差异基因可能通过多条重要的分子通路影响老年人群骨骼肌功能的下降和衰退，调控肌肉衰减综合症的发生与发展；筛选出的核心基因可为后续肌肉衰减综合症分子标记的相关研究工作提供新的线索和思路。