骨质疏松是一类常见的代谢性疾病，在绝经后妇女及老年人中常见。骨矿物质的流失导致骨密度量下降及骨微结构失衡和退化，骨骼的脆性和骨折风险大大增加^[1]。人骨髓源性间充质干细胞（human bone marrow derived mesenchymal stem cells，hBMSCs）作为一类多能干细胞，其成骨和成脂分化的方向选择对骨质疏松症进展有较大影响^[2]。

近年来，随着测序技术的发展，多种基因及其产物可以通过基因表达谱的分析和可视化与疾病相关联^[3]。在不同疾病或同一疾病不同亚型中，基因表达都有较大差异。m⁶A甲基化相关蛋白在多种疾病中存在一定的调控能力，在hBMSCs中也发挥着重要的作用^[4]。

本研究利用生物信息学技术，分析老年骨质疏松症患者hBMSCs的基因芯片，获得差异表达基因（differentially expressed genes，DEGs）谱，并构建m⁶A相关蛋白差异表达谱，为老年骨质疏松症的诊治和预防提供新的思路与方向。

1 材料与方法 1.1 数据收集与筛选

从NCBI基因表达综合（gene expression omnibus，GEO）数据库（http://www.ncbi.nlm.nih.gov/geo/）下载基因表达芯片GSE35956。GSE35956采用样本为从老年骨质疏松患者和老年非骨质疏松患者骨髓提取的hBMSCs，平台为GPL570，Affymetrix Human Genome U133 Plus 2.0 Array。

1.2 DEGs数据分析

使用R软件的Limma软件包（版本：3.40.2）分析数据集中mRNA的差异表达情况^[5]。数据选取adjusted P < 0.05且|logFC| > 2（logFC被定义为阈值mRNA差异表达的筛选）为DEGs。DEGs热图采用R软件包pheatmap绘制。

1.3 基因功能富集分析

为进一步筛选并确认DEGs的潜在功能，通过基因本体论（gene ontology，GO）对数据进行了功能富集分析。GO是一种广泛使用的工具，通过分子功能（molecular function，MF）、生物途径（biological pathway，BP）和细胞成分（cellular constituent，CC）来注释具有潜在功能的基因^[6]。京都基因和基因组数据库（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析是一种用于注释基因功能以及相关信号通路的分析方法。本研究使用R语言软件中的ClusterProfiler程序包分析DEGs的GO功能富集以及KEGG通路，气泡图使用R软件包ggord进行绘制。

1.4 m⁶A相关蛋白表达谱分析

候选的20种m⁶A相关蛋白基因来源于多种系及多种疾病的m⁶A相关分子及机制研究。m⁶A相关蛋白基因差异分析采用wilcox检验，箱线图通过R软件包ggplot2进行实现；热图通过R软件包pheatmap进行绘制。

1.5 蛋白质-蛋白质相互作用（protein-protein interaction，PPI）网络构建、集簇分析和关键（Hub）基因获取

STRING（Search Tool for the Retrieval of Interacting Genes）在线分析工具可以实现PPI网络可视化分析^[7]。在“Creative Commons BY 4.0”中，数据可以开放获取。导入已筛选DEGs进入STRING（版本11.0，https://string-db.org/），设置置信度 > 0.4，互作最大值=0，进行数据分析。获得所有结果由Cytoscape（版本3.5.1）进行可视化分析。MCODE（Molecular Complex Detection，版本1.4.2）为Cytoscape内置分析软件，可以从PPI中筛选出连接最为紧密的集簇。设置参数degree cutoff=2，node score cutoff=0.2，k-core =2，max.depth=100，取score值> 4，得到最为紧密的集簇。使用Cytoscape内置cytoHubba（版本0.1）软件进行Hub基因分析及可视化，筛选条件hubba nodes ranked取degree=10，得到10个Hub基因。

1.6 统计学分析

采用R语言进行数据分析。差异分析采用wilcoxon检验，P < 0.05为差异有统计学意义。

2 结果 2.1 DEGs的筛选和鉴定

本研究选用GEO数据集GSE35956，经Limma分析和条件筛选，共获得292个DEGs，包括93个低表达基因和199个高表达基因。通过火山图和热图进行可视化，得到明显的差异表达，DEGs个数较多，此处分别展示差异改变最大的50个上调基因和50个下调基因，见图 1。

2.2 基因功能及通路富集分析

通过对DEGs的GO及KEGG分析，针对292个基因进行了注释，分别得到了上调和下调DEGs的GO富集结果及KEGG分析结果。在GO分析中，上调的DEGs主要富集在成骨过程、结缔组织发展、细胞趋化、胞内运输的调节、细胞外结构组织等；下调的DEGs主要富集在tRNA相关代谢过程、突触的组织等。上调和下调的DEGs均可富集在成骨过程中。KEGG通路中，高表达基因富集在PI3K-Akt信号通路以及黏着斑相关通路中，低表达基因富集在Wnt信号通路、Hippo信号通路以及一些免疫调控通路中。见图 2。

2.3 m⁶A相关蛋白表达谱分析

20个候选的m⁶A相关蛋白中，通过差异分析得到ALKBH5、FTO、IGF2BP2、RBMX、VIRMA、YTHDC2、YTHDF2、ZC3H13等8个基因（图 3A），并获得基因差异分析的箱线图（图 3B）。

2.4 PPI网络及集簇分析和hub基因筛选

将STRING分析的结果导入Cytoscape软件进行进一步分析与绘制，得到PPI网络。通过MCODE（Molecular Complex Detection，版本1.4.2），设置参数为degree cutoff=2，node score cutoff=0.2，k-core =2，max.depth=100，score > 4，得到2个集簇模块，模块1包含13个节点（node），37条连线（edge），score=6.167（图 4A）；模块2包含5个节点（node），10条连线（edge），score=5（图 4B）。cytoHubba（版本0.1）进行Hub基因筛选，获得了10个Hub基因（图 4C），颜色越深，代表ranked score排名越靠前。10个Hub基因分别为CCR5、TAC1、GNGT1、RAC2、VAMP2、PXN、CSF1、TRIP10、RAB11B、ADRA2A。

3 讨论

近年来，随着研究的不断深入，越来越多的证据表明老年骨质疏松的进展与成骨分化和破骨作用失衡有关。这一过程中，hBMSCs发挥着重要的作用。本研究通过生物信息学技术，研究了老年骨质疏松患者和非骨质疏松患者的hBMSCs DEGs谱，通过表达谱的GO富集分析及KEGG通路分析，发现上调及下调的差异基因均可以富集在骨形成相关通路及细胞生长分化相关通路中；而骨小梁形成过程富集在下调的差异基因中。骨小梁的形成是骨皮质在松质骨内的延伸部分，对于成骨细胞及血管形成具有重要作用^[8]。骨小梁的减少会极大地影响成骨细胞分化和骨再生，导致骨的脆性增加^[9]。GO富集的结果显示在老年骨质疏松中，hBMSCs的差异基因对骨生成具有重要的调控作用。KEGG通路结果显示，PI3K-Akt信号通路和Wnt通路等被高度富集。PI3K-Akt信号通路可以通过调控下游不同的分子，促进间充质干细胞的成骨分化^[10-12]，但TAN等^[13]研究发现，TAZ通过抑制PI3K-Akt信号通路促进骨质疏松大鼠模型间充质干细胞成骨分化。这提示PI3K-Akt信号通路的复杂调控模式具有双向性。此外，Wnt通路及Hippo信号通路已被证明可以促进成骨分化^[14-15]，其相关基因在本研究中低表达，表明成骨分化过程被明显抑制。

m⁶A作为真核生物RNA中最丰富的甲基化修饰，在间充质干细胞成骨分化的表观遗传上具有重要的作用^{[6, 16]}。GUO等^[17]研究发现，FTO作为去甲基化酶，不仅可以对脂肪调控，还可以在骨质疏松中发挥调控作用。本研究通过对m⁶A甲基化蛋白的筛选和差异分析，获得了ALKBH5、FTO、IGF2BP2、RBMX、VIRMA、YTHDC2、YTHDF2、ZC3H13等8个基因的表达差异，包括了m⁶A甲基化蛋白的各种类型。这为进一步研究老年骨质疏松m⁶A甲基化修饰提供了思路。

通过PPI网络的构建、集簇分析和hub基因筛选得到最密集集簇和Hub基因（CCR5、TAC1、GNGT1、RAC2、VAMP2、PXN、CSF1、TRIP10、RAB11B、ADRA2A）。这些Hub基因中，CCR5、TAC1、VAMP2、CSF1等均和hBMSCs的免疫及炎症反应相关^[18]。而骨质疏松的发生发展过程中，炎症反应及其趋化因子对间充质干细胞的迁移和归巢起到趋化作用^[19]。这也提示了骨质疏松和免疫系统的相关性。Hub基因的获得为研究间充质干细胞在老年骨质疏松中蛋白互相作用调控以及信号通路提供了条件。

综上所述，本研究构建了老年骨质疏松hBMSCs DEGs谱，并分析了具有潜在靶向的信号通路及细胞生物学过程，还得到了m⁶A相关mRNA表达的差异性特征。同时，分析获得的相关差异基因在老年骨质疏松的hBMSCs成骨分化中具有较高的组织特异性和功能特异性。