基于GEO数据库芯片的老年骨骼肌减少症关键基因筛选与生物信息学分析

张涛 安云 尹露 陈霄 闫慧新 马书杰 严隽陶

(1上海中医药大学附属岳阳中西医结合医院，上海 200437；2上海市第二康复医院)

Screening and bioinformatics analysis of hub genes of sarcopenia based on GEO database chip

ZHANG Tao,AN Yun,YIN Lu,etal.

Yueyang Integrated Traditional Chinese and Western Medicine Hospital,Shanghai University of Traditional Chinese Medicine,Shanghai 200437,China

【Abstract】ObjectiveTo analyze the differential expression of hub genes in sarcopenia by using bioinformatics.MethodsThe microarray data of lateral femoral muscle samples in patients with sarcopenia was downloaded from the data platform of GEO and the differentially expressed genes(DEGs) were screened by using bioinformatics.The DEGs were subjected to gene ontology (GO) enrichment analysis and Kyoto encyclopedia of genes and genomes (KEGG) signaling pathway analysis with David online database．The protein-protein interaction (PPI) network was constructed and analyzed using the String11.0 online database. Cytoscape v3.6.1 and its plugin MCODE were used to analyze the central node protein of PPI network and hub genes in the process of sarcopenia were screened out.ResultsThrough analysis of the GSE1428 data set,1 001 DEGs were screened,of which 19 were up-regulated and 982 were down-regulated;these genes were mainly involved in the biological processes such as extrinsic apoptotic signaling pathway in absence of ligand,negative regulation of intrinsic apoptotic signaling pathway,signal transduction,peptidyl-tyrosine phosphorylation,positive regulation of epithelial cell migration,activation of protein kinase B activity.DEGs were found to be mainly concentrated in the pathways in cancer,signaling pathways regulating pluripotency of stem cells,HIF-1 signaling pathways,Rap1 signaling pathways by KEGG pathway analysis.After analyzing the PPI network,SRC,VEGFA,JUN,UBC,CASP3,CYCS,IGF1,KRAS,F2,IL-4 were identified as hub genes.ConclusionsThe selected hub genes may become an important marker for diagnosing sarcopenia or a potential therapeutic target,and provide a new theoretical basis for the pathogenesis of the sarcopenia.

【Keywords】 Sarcopenia;Bioinformatics;Differentially expressed genes;Enrichment analysis;Protein-protein interaction

骨骼肌减少症(sarcopenia)又称为肌少症、骨骼肌衰减征，是一种常见于中老年人的慢性疾病，可增加流行人群跌倒、骨折及住院的风险〔1,2〕，严重影响老年人正常的生活活动能力〔3〕，且住院、护理等费用也明显高于一般人〔4〕。该病的发病与年龄密切相关，健康人骨骼肌总质量从40～80岁可以下降大约40%，60岁后肌肉功能每年则会下降3%〔5〕。神经功能下降、生长激素分泌的改变、炎症通路激活、胰岛素抵抗、慢性疾病、脂肪浸润和营养不良等引起的肌肉蛋白合成和蛋白水解不平衡是该病的主要病因，相关的分子通路研究则围绕细胞增殖、细胞凋亡及线粒体衰退等生物学过程进行〔6,7〕。事实上，骨骼肌维持、发育的分子机制涉及多种信号通路、基因和蛋白质之间的相互作用(PPI)，维持肌肉蛋白合成和蛋白水解之间的平衡〔8〕。骨骼肌减少症发病和维持发展机制复杂，现有的研究大多围绕在个体基因研究上，但全局数据的探索上仍有很多未知。本研究通过生物信息学分析老、幼骨骼肌之间的差异基因、相关通路及关键蛋白的PPI，以确定骨骼肌减少症关键基因、通路和中枢蛋白，进一步探究骨骼肌减少症的发生发展机制。

1 对象与方法

1.1GEO芯片数据来源 以“sarcopenia”作为搜索词，在美国国立生物中心(NCBI)的GEO数据库(https://www.ncbi.nlm.nih.gov/geo)搜索已公布的人类骨骼肌减少症的基因芯片数据集。本研究使用GSE1428数据集〔9〕，该系列基于GPL96平台(HG-U133A)Affymetrix Human Genome U133A Array，包括10例年轻受试者(19～25岁)和12例年龄较大受试者(70～80岁)的股外侧肌整体基因表达谱。

1.2数据处理及筛选差异表达基因 将GSE1428数据集中的数据分组后使用GEO自带的GEO 2R分析工具进行分析，利用R语言limma〔10〕工具包输出同时满足|log2 fold change(log2FC)|>1且P<0.05的差异表达基因(DEGs)，然后根据相应平台文件所对应的R语言hgu133a.db包将探针名转化为基因名；然后使用heatmap.2工具包绘制热图，将每个DEGs的差异表达情况进行直观展示。

1.3差异基因基因本体(GO)功能富集和京都基因与基因组百科全书(KEGG)通路富集分析 将获得的DEGs上传到 David 在线数据库(https://david.ncifcrf.gov)进行 GO 富集分析及 KEGG 信号通路分析。GO 功能分析从分子功能(MF)、生物学过程(BP)及细胞组分(CC)3个部分进行富集分析〔11〕。通过 KEGG数据库对差异基因进行通路富集分析，了解疾病状态下显著改变的代谢通路。将P<0.05设定为差异显著临界值，分别选取GO富集分析及KEGG信号通路分析中基因富集数最多的10个功能进行分析。

1.4差异基因编码蛋白互作(PPI)网络分析及关键基因表达验证 采用STRING在线数据库(https://string-db.org)，将参与GO功能富集和KEGG通路富集的差异基因所编码的蛋白进行PPI分析。并利用Cytoscape软件对PPI网络进行拓扑学分析，以度值前5位的差异基因作为PPI网络中的关键基因。最后使用Cytoscape软件中MCODE插件对PPI网络进行模块分析，以MCODE分数>5分作为显著性模块的筛选标准，筛选出模块关键基因〔12〕。

2 结 果

2.1差异基因的筛选 通过R语言3.6.2软件将原始芯片归一化处理，利用limma包对骨骼肌减少症和正常肌肉组织归一化的数据进行差异表达分析，最终得到1 001个差异基因。其中上调19个，下调982个，见图1。使用R语言中的gplots包绘制差异基因热图，见图2。

2.2差异基因GO分析和KEGG通路分析 DEGs参与的生物学进程主要包括缺乏配体的外在凋亡信号通路、内在凋亡信号通路的负调控、信号转导、肽基酪氨酸磷酸化、蛋白激酶B活性的激活等；细胞组分分析显示DEGs主要参与原生质膜的组成、细胞表面、HFE-转铁蛋白受体复合物、微管相关复合物、线粒体小核糖体亚基、细胞溶质等的组成；参与的分子功能主要包括蛋白质结合、蛋白质均二聚活性、转运活性、过氧化物酶体增殖物激活受体结合、肌动蛋白结合、生长因子活性、蛋白磷酸酶抑制剂活性等。

进行KEGG通路分析后发现，DEGs主要集中在癌症通路、调节干细胞多能性的信号通路、HIF-1信号通路、碳代谢、Rap1信号通路等，见图3，图4。

图1 骨骼肌减少症差异表达基因火山图

绿色为下调，红色为上调图2 骨骼肌减少症最大的前50 个基因热图

圆圈的大小表示参与基因的数量；颜色表示参与该功能的显著性图3 骨骼肌减少症差异基因GO功能分析和显著参与的KEGG信号通路

图4 骨骼肌减少症相关KEGG通路分析及相应基因

2.3差异基因PPI网络 将上述筛选出的DEGs上传到String v11.0在线数据库，转换差异基因编码PPI网络，共包含有711个节点和2 956条边。通过Cytoscape v3.6.1软件构建PPI网络；进行拓扑学分析，筛选度值(Degree)>20的蛋白分子进行互作网络分析，共包含64个节点和385条边，见图5。度值前10的hub基因分别为原癌基因，非受体酪氨酸激酶(SRC)、血管内皮生长因子(VEGF)、原癌基因，AP-1转录因子亚基(JUN)、泛素(UB)C、半胱氨酸天冬氨酸蛋白酶(CASP)3、细胞色素c体细胞(CYCS)、胰岛素样生长因子(IGF)1、GTPase，原癌基因,GTPase(KRAS)、凝血因子Ⅱ，凝血酶(F2)、白细胞介素(IL)-4，其相互作用关系见图6。使用Cytoscape中MCODE插件对子模块进行分析，可发现一些在骨骼肌减少症局部调控中起重要作用的关键基因，见图7。

图5 蛋白质相互作用网络分析

圆圈的大小表示该基因的度值；红色为上调基因，绿色为下调基因图6 骨骼肌减少症hub基因差异PPI 分析

颜色越深代表度值越大图7 骨骼肌减少症差异基因蛋白网络的子网络分析(包含3个模块)

3 讨 论

骨骼肌减少症是一种以肌肉质量和肌肉功能丧失为特征的年龄相关性疾病，会使老年人跌倒、功能受限、残疾和死亡等不良后果的风险增加〔13〕。正确理解该病的发病机制有助于制定新的诊断和治疗策略。本研究通过对微阵列数据进行生物信息学分析，识别出老年患者和年轻人骨骼肌之间的差异基因，并通过进行GO富集分析和KEGG途径富集分析揭示DEGs的功能，探索老年人骨骼肌减少症的相关基因、蛋白和通路。

通过KEGG通路分析发现影响骨骼肌减少症发病的相关基因富集最多的信号通路为癌症通路，该通路是过氧化物酶体增殖剂激活受体(PPARs)、丝裂原活化蛋白激酶(MAPK)、HIF-1、PI3K-Akt-mTOR、转化生长因子(TGF)-β、JAK-STAT等骨骼肌减少症的发生相关信号通路的共同交汇通路，参与骨骼肌细胞的增殖、应激、炎症、分化、转化、凋亡等多种生物过程。调节干细胞多能性的信号通路也参与了其中某些过程。此外，研究结果显示DEGs还在HIF-1信号通路上富集，该通路活性的调控机制受到活性氧和氮氧化物等多种因素的影响，在骨骼肌细胞凋亡和自噬过程中同样发挥了重要作用。通过激活调节葡萄糖转运和糖酵解的基因，HIF-1还可以提高肌细胞对缺氧的适应性，并增强VEGF及其受体的表达，以应对氧供应减少〔14〕。研究表明，在暴露于严重的低压低氧前10 min抑制HIF-1转录因子，可以减轻骨骼肌损伤〔15〕。Rap1信号通路在多种细胞类型的细胞黏附和整联蛋白功能中起重要作用，在细胞间和细胞-细胞外基质相互作用的控制中发挥主导作用〔16〕。调控MAPK或PI3K-AKT通路活性促进细胞的增殖〔17,18〕，调节VEGF诱导的血管生成〔19〕，并参与肌动蛋白细胞骨架的生成〔20〕。本研究显示Rap1信号通路相关基因在骨骼肌减少症样本中出现下调，体现了其对于骨骼肌稳态的积极作用。

通过构建PPI网络，研究确定了骨骼肌减少症中关键基因。研究显示，SRC激酶在骨骼肌减少症样本中下调。作为第1个被识别并与细胞增殖调节相关的酪氨酸激酶，通过参与激活MAPK、PI3K等信号通路促进细胞增殖〔21〕。同时，SRC激酶可以通过抑制糖原合成酶(GSK)-3来发挥其抗凋亡作用，防止细胞蛋白质的过度分解〔22〕。此外，SRC还参与细胞黏附生长的信号传导，促进肌动蛋白细胞骨架的调节，维持其正常的生理功能〔23〕。KRAS的活化能够激活一系列下游通路，尤其是RAF/MEK/ERK和PI3K/AKT通路，也可促进细胞分裂增殖和抑制细胞凋亡〔24,25〕，KRAS表达的降低可能是肌肉减少症的诱发因素。JUN作为JAK/STAT信号通路成员，负调控骨骼肌细胞的增殖，在骨骼肌减少症中表达上调。有报道证实了JAK/STAT信号通路在老年人骨骼肌肌卫星细胞中的激活〔26〕。但也有研究表明，JUNB能够抑制蛋白质的降解，增加JunB表达可显著提高肌纤维的横截面积〔27〕。

尽管肌细胞的凋亡会引起骨骼肌的氧化损伤和肌纤维的萎缩坏死,但这一生物过程在降解受损炎性细胞、限制组织损伤、恢复组织完整性等方面同样发挥积极作用。细胞正常的凋亡过程在骨骼肌的正常发育、自稳平衡方面起非常重要的作用〔28〕，本研究也发现了一些促进细胞凋亡的关键基因在骨骼肌减少症肌肉样本中表达下调。CYCS从线粒体内释放进入到细胞质，在ATP作用下激活凋亡蛋白酶活化因子(Apaf)-1，Apaf-1活化后依次激活Caspase-9和Caspase-3，是肌细胞发生凋亡的早期诱因〔29〕，而CYCS基因在本研究中显示表达下调。UBC可以通过酶促反应相互连接，进而介导靶蛋白降解，泛素化蛋白降解途径可介导细胞凋亡〔30〕。目前没有文献证明泛素蛋白酶体系统对于骨骼肌减少症的促进作用，但UBC可以识别、标记并降解错误折叠蛋白，减轻细胞的损害，在应对细胞应激反应时可发挥重要作用〔31〕。IL-4可以通过上调Bcl-2蛋白家族调控骨骼肌细胞的自噬和凋亡，在骨骼肌内稳态发挥重要作用〔32〕。而IL-4表达的降低可能会使Bcl-2控制细胞程序性死亡的作用减弱,降低骨骼肌线粒体更新速率及抗氧化水平〔33,34〕。本研究结果提示，维持细胞正常凋亡的相关基因、蛋白表达减少在骨骼肌减少症的发病中有一定促进作用。

作为SRC和KRAS的上游调控因子，VEGFα具有促进内皮细胞增生、迁移及增加血管通透性的作用，可维持血管正常的形态和完整性〔35〕。研究发现，VEGF表达降低是老年毛细血管稀疏和血管生成减弱的潜在原因，而微循环是否正常也是影响骨骼肌发育的重要因素，阻断VEGF会破坏由超载或剪应力引起的骨骼肌血管生成〔36〕。同时，VEGF还可能通过VEGF受体Flk-1和Flt-1作用于肌卫星细胞，刺激肌卫星细胞迁移并抑制肌卫星细胞凋亡，维持骨骼肌的稳态，促进骨骼肌再生〔37〕。因此，在骨骼肌减少症的研究和治疗中，应重视微循环对肌肉可塑性、修复和再生的作用。结合之前通路分析，衰老相关VEGF表达减少导致血管生成受损的主要原因可能是老年人HIF-1、Rap1活性缺陷〔38〕。

肌肉减少症的发展与多种合成代谢激素水平下降有关，而30岁以后睾丸素的生成速度将会以每年1%的速度减少〔39〕,而这种由衰老而出现生长激素下降也会使IGF-1的表达降低，最终影响骨骼肌蛋白的合成。IGF-1通过参与IGF1/PI3K/Akt/PKB/mTOR信号通路增强肌肉蛋白质合成〔40〕，并负性调控转录因子FoxO家族以减少蛋白质的分解〔41〕。实验显示肌肉特异性胰岛素和IGF-1受体双敲除小鼠会出现FoxO诱导的骨骼肌大量减少〔42〕。除此之外，本研究发现影响骨骼肌减少症生成的基因还参与了肌动蛋白细胞骨架通路的调节。细胞骨架蛋白作为骨骼肌细胞重要的结构，可以连接和锚定肌细胞内的结构成分〔43〕。肌动蛋白作为细胞骨架网中的一员，在骨骼肌减少症的发生发展中有十分重要的作用。关键基因F2通过与细胞表面F2受体和整合素结合激活肌动蛋白细胞骨架通路〔44〕，并与其受体结合增加了肌动蛋白应力纤维的形成，调节了肌动蛋白细胞骨架组织〔45〕。

本研究通过模块分析还发现一些新的与骨骼肌减少症有关的基因，如中心体蛋白(CEP)290、细胞分裂周期(CDC)6、G蛋白耦联受体(GPR)17、G蛋白亚基(GNB)5、溶血磷脂酸受体(LPAR)2、甲酰基肽受体(FPR)2等，可能作为骨骼肌减少症潜在诊断和治疗的靶点。这些基因在骨骼肌减少症中研究较少，仅有少量研究证实CDC6主要参与DNA的复制，对于维持细胞周期的正常运转和细胞增殖具有重要作用〔46〕。有学者证实了CDC6对于肌卫星细胞增殖有促进作用〔47〕。

综上，影响骨骼肌减少症发病的某些基因涉及了癌症通路，可能与癌症的发生有密切联系。此外，在诊断和治疗骨骼肌减少症时，应重视细胞的正常增殖与凋亡、生长激素水平及局部微循环对骨骼肌的修复和再生作用。本研究的局限性在于所选数据纳入样本量较少，且未进行不同性别的差异基因分析。所筛选出的关键基因是否能成为骨骼肌减少症潜在的靶点，还需要进一步进行分子机制的研究证实。