基于生信分析构建骨质疏松氧化应激相关miRNA-mRNA调控网络*

莫 亮, 王章正, 马 超, 何 伟, 周 驰, 刘予豪△

骨质疏松症(osteoporosis,OP)是以骨量降低和骨组织微结构破坏为特征的代谢性骨病[1]。随着对氧化应激在骨代谢疾病中作用的深入研究,目前认为,由活性氧物质(reactive oxygen species,ROS)介导的氧化应激不但能造成成骨细胞活性下降,诱导成骨细胞和骨细胞凋亡从而减少骨形成,也能通过激活破骨细胞前体的分化和提高破骨细胞活性,加强骨吸收[2],导致OP的发生。微小核糖核酸(micro-ribonucleic acid,microRNA,miRNA)是小的非编码RNA,能够通过抑制信使核糖核酸(messenger RNA,mRNA)翻译或促进mRNA降解来调节转录后基因表达,许多正常的生理过程和疾病发生有赖于miRNA的调控[3]。研究表明miRNA可以调节骨重塑并在OP中发挥重要的作用[4-5]。而miRNA的转录、生物合成、转运和功能与ROS高度相关,ROS可以诱导或抑制miRNA的表达,并通过调控靶基因来调节下游生物学功能[6-7]。部分miRNA能通过影响骨细胞氧化应激水平,促进OP的发展[8-9]。目前,针对氧化应激导致OP的相关基因和(或)通路的研究以体外或体内实验为主,国内外尚未见从生物信息学角度对氧化应激和OP的关系展开综合分析的研究报道。本文通过对公共数据库中OP和氧化应激相关基因表达数据进行整合分析,构建OP氧化应激相关miRNA-mRNA调控网络,分析氧化应激导致OP的潜在机制,以期为抗氧化应激治疗OP提供相关靶点和通路。

1 资料和方法

1.1 数据来源

从GEO公共数据库(https://www.ncbi.nlm.nih.gov/geo/)下载基因表达芯片GSE74209、GSE35958和GSE35956。GSE74209包含OP和非OP患者的新鲜股骨颈骨小梁miRNA测序数据,平台为GPL20999。GSE35958以及GSE35956将来自于OP患者股骨头骨髓中分离的骨髓间充质干细胞和来自于非OP患者的骨髓间充质干细胞进行比较,平台为GPL570。此外,在AmiGO数据库(http://amigo.geneontology.org/)通过关键词“oxidative stress”获得17个氧化应激相关miRNA及446个氧化应激相关mRNA。AmiGO可用于搜索和浏览Gene Ontology数据库(http://geneontology.org/),允许用户查询相关基因功能及其注释数据[10-11]。基因表达芯片具体信息见表1。

表1 基因表达芯片信息Table 1 Information of GEO data

1.2 差异基因分析

使用GEO2R在线分析工具对芯片GSE74209和GSE35958的数据进行差异分析,获得差异miRNA(differentially expressed miRNAs,DEMs)和差异基因(differentially expressed genes,DEGs)。将GSE74209得到的数据取adj.Pvalue<0.05且|logFC|≥1为DEMs;GSE35958得到的数据取adj.Pvalue<0.05且|logFC|≥1为DEGs,并将得到的数据运用R语言进行可视化。GSE35956作为验证集下载其相应的表达矩阵。

1.3 氧化应激相关DEGs和DEMs的筛选

将GSE35958筛选得到的DEGs和AmiGO数据库得到的氧化应激相关基因取交集获得OP氧化应激相关DEGs。同样将GSE74209筛选得到的DEMs和AmiGO数据库得到的氧化应激相关miRNA取交集获得OP氧化应激相关DEMs。通过R语言进行韦恩图可视化。

1.4 功能富集分析

为了探讨氧化应激相关DEGs的潜在作用,对其进行GO(Gene Ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes)功能富集分析,GO富集分析可以通过注释基因或其产物、识别基因芯片数据而分析得到生物学特征结果[12],本研究重点关注GO分析结果中的生物学过程。KEGG分析是常用的基因通路分析方法,可以将基因注释和富集在信号通路上进行研究[13]。本研究应用Metascape在线网站(http://metascape.org/)进行富集分析,Metascape是一个整合了GO、KEGG、Uniprot等多个数据库的工具,可用于GO分析和KEGG信号通路分析,并且具有时效性和准确性[14]。筛选条件为最小重叠(min overlap)=3和最小富集(min enrichment)=1.5,P<0.01被认为差异具有统计学意义。通过R语言绘制柱状图展示前15条通路。

1.5 蛋白互作网络构建和分析

为了阐明氧化应激和OP的关系,进行蛋白质之间功能相互作用分析。STRING(http://strin g-db.org)是一款在线分析工具,可以用来构建蛋白互作网络(protein-protein interaction network,PPI网络)和进行在线分析[15]。将筛选得到的氧化应激相关DEGs导入STRING数据库,设置置信度(confidence score)≥0.4,得到相应数据,并使用Cytoscape软件构建PPI网络。Cytoscape(版本3.7.1)是一个用于可视化分子相互作用网络的开源生物信息学软件平台[16]。

1.6 骨质疏松中关键氧化应激基因的筛选和验证

通过使用Cytoscape中cytoHubba插件来识别关键基因。使用了7种常用算法(MCC,MNC,Degree,Closeness,Radiality,Stress,EPC)来评估和选择关键基因。随后通过GeneMANIA(http://www.genemania.org/)构建了这些关键基因的共表达网络,GeneMANIA是一个分析基因内部关联的在线网站[17]。随后通过数据集GSE35956验证关键基因的表达。

1.7 构建miRNA-mRNA互作用网络

miRTarbase(http://mirtarbase.mbc.nctu.edu.tw/php/index.php)是一个经过实验验证的miRNA-mRNA相互作用数据库,其中包含有实验数据(报告分析、蛋白质印迹、微阵列或pSILAC)支持的4076个miRNA和23045个靶基因相互作用数据[18]。将氧化应激相关miRNA导入miRTarbase数据库获得相应靶基因,并和氧化应激相关DEGs相映射构建miRNA-mRNA调控网络,通过Cytoscape软件进行可视化。

2 结果

2.1 OP氧化应激相关DEGs的筛选及PPI网络构建

对数据集GSE35958进行GEO2R差异表达分析,经过条件筛选得到2527个DEGs,通过火山图进行可视化(图1A)。将其与AmiGO数据库中收集到的443个氧化应激相关基因取交集,最终得到86个OP氧化应激相关DEGs,包含69个上调基因和17个下调基因(图1B和表2)。通过STRING和Cytoscape软件分析与绘制,得出PPI网络(图1C)。

A:GSE35958差异表达火山图,红色点代表上调DEGs,蓝色点代表下调DEGs;B:氧化应激相关基因和差异表达基因韦恩图;C:氧化应激相关DEGs蛋白互作网络图图1 骨质疏松氧化应激相关DEGs筛选Fig.1 The selection of oxidative stress-related DEGs in osteoporosis

A:pFastBacTM Dual +[β2m+MR1/IgG1 Fc]重组质粒结构;B:MR1二聚体结构模拟图;C:将MR1二聚体装载到与兔抗人IgG Fc抗体结合的dynalbeads上,形成MR1 aAPC图1 MR1/IgG1 Fc表达载体、二聚体及MR1 aAPC结构模拟图Fig.1 Structure paradigms of pFastBacTM Dual +[β2m+MR1/IgG1 Fc]plasmid,dimeric MR1/IgG1 Fc fusion protein and MR1 aAPC

表2 骨质疏松氧化应激相关DEGsTable 2 Oxidative stress-related DEGs in osteoporosis

2.2 GO功能与KEGG信号通路富集分析

根据GO和KEGG分析,以P<0.01为筛选条件并按P值从小到大排列,将GO和KEGG分析结果中排名前15的富集通路用柱状图表示(图2)。GO分析中,DEGs主要富集在氧化应激和活性氧相关生物过程。KEGG分析中,DEGs主要富集在活性氧通路、癌症通路、破骨细胞分化通路和细胞凋亡等通路中。

A:KEGG富集分析;B:GO富集分析图2 氧化应激相关DEGs功能富集分析Fig.2 Enrichment analysis of oxidative stress-related DEGs

M:DNA marker;1:β2m(片段长度为360 bp);2:MR1/IgG1 Fc(片段长度为1902 bp);3:Bacmid;4:Bacmid+[β2m+MR1/IgG1 Fc](片段长度为4822 bp)图2 重组质粒的菌液PCR扩增产物DNA电泳鉴定Fig.2 DNA bands of PCR products from recombinant plasmid

2.3 关键基因的筛选及验证

CytoHubba是Cytoscape软件中用于识别关键基因的插件,提供了12种分析算法来计算PPI网络图中的关键基因。本文利用其中7种算法筛选出12个重叠的关键基因(图3A),其中FOXO3、FYN、AKT1、TNF、MAPK1、MAPK3、TP53和EGFR的表达在GSE35958中较对照组上调,TXN、JUN、CTNNB1和APP的表达在GSE35958中较对照组下调。基于GeneMANIA数据库,对这些关键基因进行了共表达网络和相关功能分析,结果表明关键基因及其共表达基因主要与氧化应激和活性氧功能相关(图3B)。

A:韦恩图表示通过7种算法筛选出12个重叠的关键基因;B:关键基因及其共表达基因通过GeneMANIA进行分析图3 关键基因的筛选Fig.3 The selection of hub genes

M:marker;A:10 μg /mL羊抗人IgG Fc抗体包被ELISA检测孔,HRP标记羊抗人IgG Fc抗体(1∶15000稀释)作为检测抗体进行ELISA;B:蛋白质电泳胶及考马斯亮蓝染色,箭头指示为MR1/IgG1 Fc二聚体;C:Western blot检测用HRP标记羊抗人IgG Fc抗体(1∶20000稀释)的MR1二聚体图3 ELISA及Western blot检测Sf9细胞MR1/IgG1 Fc融合蛋白Fig.3 MR1/IgG1 Fc fusion protein expression in Sf9 cells detected by ELISA and Western blotting

为了验证关键基因的表达水平,选择另外一个包含OP和非OP的骨髓间充质干细胞基因表达数据集(GSE35956)作为验证,结果表明EGFR、FOXO3、FYN、AKT1在OP组显著上调,APP和TXN在OP组中显著下调,同GSE35958数据集中的表达趋势相同。而CTNNB1和MAPK1则和数据集GSE35958中的表达趋势相反(图4)。

A:流式细胞术结果代表图;B:MR1二聚体包被率统计结果;Blank:未经处理的dynalbeads;Control:PE标记同型对照抗体作为对照抗体排除非特异性染色;Fc:APC标记羊抗兔IgG检测兔抗人IgG Fc结合在dynalbeads效率;MR1:PE标记的MR1抗体检测MR1/IgG1 Fc二聚体的包被效率图4 流式细胞术检测MR1/IgG1 Fc二聚体Fig.4 MR1/IgG1 Fc dimer detected by flow cytometry

OP:骨质疏松组;Control:对照组;两组数据之间的比较使用均值t检验,*P< 0.05,**P< 0.01,***P< 0.001图4 数据集GSE35956中关键基因的表达差异Fig.4 Expression difference of hub genes in GSE35956 database

2.4 差异miRNA及氧化应激相关miRNA的筛选和鉴定

经GEO2R初步分析和条件筛选,在芯片GSE74209共获得142个DEMs,包括43个下调miRNA和99个上调miRNA。和AmiGO数据库获得的17个氧化应激相关miRNA取交集后获得5个共同miRNA,包括hsa-miR-675-5p、hsa-miR-155-5p、hsa-miR-210-3p、hsa-miR-92a-3p和hsa-miR-17-3p,都属于表达上调的miRNA。利用R语言将得到的数据通过火山图和热图进行可视化(图5)。

A:流式细胞检测代表结果;B:流式细胞检测统计结果;Control:未经处理的PBMCs;DH5α supernatant/MR1 aAPCs:加载DH5α培养上清的MR1 aAPCs与PBMCs共培养组;DH5α water-soluble/MR1 aAPCs:加载DH5α裂解后水溶性代谢物的MR1 aAPCs与PBMCs共培养组;DH5α lipid-soluble/MR1 aAPCs:加载DH5α裂解后脂溶性代谢物的MR1 aAPCs与PBMCs共培养组;*P<0.05图5 加载DH5α代谢物的MR1 aAPCs促进MAIT细胞活化Fig.5 MR1 aAPCs loaded with DH5α metabolites promotes MAIT cell activation

A:氧化应激相关miRNA和差异表达miRNA韦恩图;B:GSE74209中DEMs火山图;C:GSE74209中5个氧化应激相关miRNA表达热图图5 骨质疏松氧化应激相关miRNA筛选Fig.5 The selection of oxidative stress-related miRNAs in osteoporosis

2.5 miRNA-mRNA调控网络

通过miRTarbase数据库获得5个氧化应激相关miRNA预测靶基因,将其和氧化应激相关DEGs进行比对和匹配,使用Cytoscape软件构建miRNA-mRNA调控网络,并通过logFC值判断其表达是上调还是下调,得到可视化网络(图6)。has-miR-210-3p通过映射无相应靶基因,因此不在调控网络中。has-miR-92a-3p调控下游12个DEGs,has-miR-155-5p调控下游13个DEGs,has-miR-17-3p调控下游11个DEGs,has-miR-675-5p调控下游2个DEGs。其中has-miR-155-5p不但调控最多的下游DEGs,而且调控最多的关键基因(JUN、FOXO3、EGFR、TXN、TNF)。考虑到FOXO3、EGFR、TXN已在数据集GSE35956中得到验证,因此推测,has-miR-155-5p-FOXO3EGFRTXN可能是氧化应激诱导OP的重要途径。

红色代表表达上调,蓝色代表表达下调,V形代表DEMs,椭圆形代表DEGs,三角形代表关键基因图6 miRNA-mRNA网络调控图Fig.6 miRNA-mRNA regulatory network

A:流式细胞术检测MAIT细胞增殖情况代表图;B:流式细胞术检测结果统计;Control:染色CFSE的PBMCs;DH5α supernatant/MR1 aAPCs:加载DH5α培养上清的MR1 aAPCs与PBMCs共培养组;DH5α water-soluble/MR1 aAPCs:加载DH5α裂解后水溶性代谢物的MR1 aAPCs与PBMCs共培养组;DH5α lipid-soluble/MR1 aAPCs:加载DH5α裂解后脂溶性代谢物的MR1 aAPCs与PBMCs共培养组;*P<0.05图6 加载DH5α代谢物的MR1 aAPC促进MAIT细胞增殖Fig.6 MR1 aAPC loaded with DH5α metabolites promotes proliferation of MAIT cells

3 讨论

高通量分析技术和生物信息学方法正被应用于生物医学研究的所有领域,对数据库中存储的大量生物数据进行分析,有利于加深对疾病的理解,助力疾病治疗和诊断[19]。现有研究尝试从抗氧化应激方向上探索治疗OP的新型药物[20],然而氧化应激导致OP的具体分子机制目前尚不明确。本研究利用生物信息学整合生物数据库信息构建OP氧化应激相关miRNA-mRNA调控网络,以期为抗氧化治疗OP提供有希望的靶点和通路。

首先,通过对数据集GSE35958进行分析和筛选得到了86个氧化应激相关DEGs。GO分析证实了这些DEGs主要富集在氧化应激相关的生物过程。KEGG分析发现DEGs富集在活性氧(ROS)、癌症、破骨细胞分化和细胞凋亡等通路上。氧化应激是人体在代谢过程中ROS的产生和清除失衡所致,当机体由于衰老、疾病等原因导致ROS大量积聚时,氧化应激随之产生。破骨细胞介导的骨吸收大于成骨细胞介导的骨形成是发生OP的直接原因,ROS可促进破骨细胞的分化和增殖[21]。抗氧化剂预处理后的破骨细胞在减少ROS产生的同时,可显著抑制破骨细胞前体的分化和成熟[22]。ROS同样可调节成骨细胞凋亡信号导致成骨细胞的凋亡[23],而抗氧化剂可逆转氧化应激诱导的成骨细胞凋亡并增加成骨细胞活性[24]。因此,这些参与了氧化应激的相关DEGs可能通过介导破骨细胞分化或成骨细胞凋亡等生物过程促进OP的发展,而抗氧化治疗可以减缓这个过程。

为了进一步挖掘氧化应激导致OP的核心靶点,在构建PPI互作网络后,通过7种算法在PPI网络中筛选出了12个关键基因。其中上调基因包括FOXO3、FYN、AKT1、TNF、MAPK1、MAPK3、TP53和EGFR,下调基因包括TXN、JUN、CCTNB1和APP。通过进一步验证,发现EGFR、FOXO3、FYN、AKT1、APP和TXN在数据集GSE35956中呈现同样的表达趋势。最新研究发现EGFR信号通路是机体调节骨代谢的重要途径,在成骨细胞分化过程中,EGFR的激活降低了成骨细胞转录因子Runx2和Osterix的表达[25]。同时EGFR信号通路刺激破骨细胞前体细胞的增殖,EGFR基因缺陷小鼠表现出破骨细胞募集延迟、增殖和生成减少[26]。FOXO家族蛋白在氧化应激状态下通过上调抗氧化酶的表达来降低ROS产生[27]。在骨髓间充质干细胞成骨分化过程中,FOXO3基因通过诱导线粒体自噬调节细胞内活性氧水平,从而影响其成骨分化[28]。FYN在破骨细胞前体细胞中表达,通过介导FcRγ/DAP12的磷酸化调控RANKL刺激的破骨细胞形成[29]。AKT1是骨细胞中主要存在的AKT亚型,在骨髓间充质干细胞中敲除AKT1基因,会导致成骨转录因子Runx2的表达增加,刺激成骨细胞分化[30]。而AKT的下降同时会激活GSK3β,减少下游NFATc1入核,从而抑制破骨细胞分化[31]。这表明AKT1既是成骨细胞分化的负调节因子,也是成骨细胞偶联破骨细胞生成的强阳性介质。APP在各类细胞中普遍表达,研究发现APP基因缺陷小鼠骨小梁和皮质骨骨量减少,进一步研究发现APP基因缺陷小鼠成骨细胞线粒体受损,并伴有ROS和细胞凋亡的增加,提示APP与ROS诱导的OP相关[32]。TXN在细胞抗氧化防御中扮演重要角色,研究发现TXN表达水平的升高,与骨髓间充质干细胞的成骨分化相关,提示TXN编码蛋白可能正向促进成骨分化[33]。综上所述,现有研究表明这些氧化应激相关的关键基因与骨代谢存在直接或间接联系,可能在OP中发挥重要作用,然而笔者发现对于它们在人体中功能的研究相对较少,需要进一步的研究验证它们对人OP的潜在作用。

最后,通过对数据集GSE74209的分析,以及进一步的筛选得到了5个在OP中表达上调且与氧化应激相关的miRNA。通过miRTarbase数据库获得的miRNA预测靶点和氧化应激相关DEGs相匹配成功构建miRNA-mRNA调控网络。结果发现大部分下游DEGs被has-miR-17-3p、has-miR-92a-3p和has-miR-155-5p所调控,表明它们可能在氧化应激诱导OP过程中发挥了重要作用。已有体外实验发现miR-17-3p可通过靶向调节SRY-box转录因子6(Sox6)从而抑制骨髓间充质干细胞的成骨分化[34]。miR-92a-3p被证实可通过调节ROS水平参与炎性疾病的发生和发展[35],在骨折小鼠中靶向抑制miR-92a的表达可增加血管形成并加快骨折愈合[36]。对于miR-155-5p,在miRNA-mRNA调控网络中受其调控的下游DEGs最多,且包括多数关键基因,因此其可能与氧化应激诱导OP的关系更加密切。研究发现miR-155-5p可通过调控SOCS1增强TNF-α抑制的成骨分化[37],同时影响破骨细胞分化[38]。此外,has-miR-155-5p在绝经后女性血清中显著上调,被认为可能作为绝经后女性OP的诊断标志物[39]。同时,体内和体外实验已证实miR-155可抑制破骨细胞的活性,因此具有抗OP的治疗潜力[40]。可见,上述miRNA有望成为未来抗氧化应激预防OP的潜在靶点。

综上所述,本研究通过构建氧化应激诱导OP的miRNA-mRNA的互作网络,发现了has-miR-155-5p-FOXO3EGFRTXN等具有OP治疗潜力的靶向通路,为抗氧化应激治疗OP提供了有希望的靶点和通路。