高低峰值骨量人群生物标志物差异及其在骨质疏松中的诊疗价值

林适 袁嘉尧 林贤灿 杨彬彬 林燕平 黄佳纯 连晓航 万雷 黄宏兴*

1.广州中医药大学第三临床医学院，广东 广州 510000 2.广州中医药大学第三附属医院，广东 广州 510240

骨质疏松症(osteoporosis，OP)是一种与年龄相关，以骨强度下降为特征的全身性骨骼疾病[1]。OP发病机制主要是成骨和破骨进程间的平衡被打破，骨吸收增加而骨形成减少，从而导致骨量降低、骨组织微结构恶化，骨脆性增加，以致骨折发生率升高[1-2]。OP是老年人骨折最常见的原因，OP相关性骨折好发于脊柱、髋部及腕部，构成巨大的医疗、公共卫生和经济负担[2-3]。由于人口老龄化的加剧和过去二三十年来我国生活方式的巨大变化，OP患者数量预计在不久的将来会激烈增长，我国公共卫生将面临严峻挑战[4]。高通量技术的最新进展使许多生物化学骨生物标志物的识别成为可能，有助于OP患者的诊断和管理[5]。目前，确定更多有用的生物标志物或靶标，以改善OP的诊断和管理仍是当务之急[6]。

峰值骨量(peak bone mass，PBM)是人体整个生长发育期最大的骨积累量，是个体未来骨质疏松症和骨折风险的重要预测指标[7]。一般来说，男性和女性的骨量在前20年显著增加，并在青春期晚期或青年期达到稳定[8]。成人的骨量等于成年早期的峰值骨量减去成年后的骨质流失量，提高PBM对于预防或者延缓OP发病的重要性不言而喻。流行病学研究结果显示，PBM值每增加10%，可使50岁后骨折风险降低50%[9]。低PBM值与更高的OP发病率和骨折风险相关，而因为高PBM值能够为成年人和老年人提供更大量的骨质储备及更高的骨强度，能明显减少或延迟OP的发生[8]。因此，如何在儿童及青少年时期实现PBM值的提高，为中老年时期提供更高的骨量储备，实现高骨密度和骨强度积累更有利于OP及OP相关骨折的预防是目前的一个热门研究话题。

目前大多数研究都着眼于骨质疏松人群和健康人群生物标志物的差异，鲜有研究分析高低峰值骨量人群生物标志物差异。鉴于峰值骨量在骨质疏松发病及防治中的重要地位，本研究基于生物信息手段分析高低峰值骨量人群间生物标志物的差异，并研究其在骨质疏松诊断治疗中的价值。

1 资料和方法

1.1 高低峰值骨量基因芯片下载

通过基因组学存储数据库GEO(Gene Expression Omnibus，https://www.ncbi.nlm.gov/gds)检索“osteoporosis”“peak bone mass”等关键词，设置样本来源为人，得到符合条件的数据集GSE7158，平台文件为：GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array。该数据集招募了878名20～45岁的健康中国女性，平均27.3岁，使达到并维持峰值骨量，根据PBM表型分布进行受试者筛选，最后分别选中12名极低和14名极高PBM受试者进行测序分析。

1.2 数据预处理及差异表达分析

通过Perl脚本(www.perl.org/)处理下载得到的数据集数据，并获得基因表达矩阵文件。利用R软件中的limma工具包进行差异表达基因(DEGs)筛选。DEGs筛选标准为：∣logFC∣>2并且P<0.01。

1.3 GO功能和KEGG信号通路分析

为了明确筛选得到的差异基因的生物学作用，利用R语言中的clusterProfiler工具包进行差异基因GO(gene ontology)功能富集分析。GO是基因功能国际标准分类体系，它旨在建立适用于各种物种的，对基因和蛋白质功能进行限定和描述的，并能随着研究不断深入而更新的语言词汇标准[10]。GO分为分子功能(molecular function，MF)、生物过程(biological process，BP)和细胞组成(cellular component，CC)3个部分[10]。采用KOBAS 3.0在线工具(http://kobas.cbi.pku.edu.cn/)进行KEGG通路富集分析[11]。

1.4 蛋白互作(PPI)网络分析及关键基因筛选

通过String数据库(https://www.string-db.org/)进行差异基因PPI网络分析[12]，设置去除无连接的节点，得到PPI网络可视化图谱及PPI网络注释文件。将注释文件导入Cytoscape(版本3.7.1)软件中，并利用cytoHubba插件10种算法计算网络种节点得分[13]。利用R语言中的UpSetR工具包对每个算法分值最高的前25个基因进行交叉验证，最终得到各个算法均定义为关键基因的列表。然后再进行关键基因PPI网络可视化展示。

1.5 关键基因主成分分析及组间表达比较

为进一步明确筛选得到的关键基因是否能够代表高低峰值骨量人群基因表达的差异，本研究进行了关键基因在高低峰值骨量两组人群间的表达主成分分析(principal component analysis，PCA)。PCA是多元统计中的一种统计方法，是最常用的线性降维方法，具有降低数据维度同时保留较多的原始数据点的特性[14]。本研究中PCA的运用，可以有效降低基因表达特征集维数的同时，也可以很好地解释各个基因表达特征集与不同样本的相关性。

1.6 关键基因和骨质疏松的相关性

为验证关键基因和骨质疏松之间的相关性，从GEO数据库中以“osteoporosis”为关键词，设置样本来源为人，检索并下载得到基因表达数据集GSE35959，平台文件和高低骨量数据集为同一平台(GPL570)，通过该数据集验证关键基因在骨质疏松患者和健康者之间的表达差异。

2 结果

2.1 高低峰值骨量人群差异表达基因筛选

通过R语言中的limma工具包进行差异表达基因筛选，设置筛选条件为：∣logFC∣>2并且P<0.01。通过分析，一共得到182个差异基因，其中73个下调差异基因，109个上调差异基因，差异表达热图见图1A及火山图见图1B。

2.2 差异基因GO和KEGG富集分析结果

为发掘差异基因的生物学作用，采用R语言中的clusterProfiler工具包进行了GO功能注释，对各GO功能显著性靠前的5条功能注释进行展示(图2)，KOBAS3.0在线工具进行了KEGG通路富集分析，并进行相似通路聚类分析(图3)。差异基因GO功能注释结果显示，生物过程主要涉及DNA的复制、氨基酸的生物合成、脱氧核糖核酸酶活性的调控和蛋白质泛素化的负调控等；细胞组分主要定位在基底外侧质膜、神经元细胞体膜、细胞体膜、核复制叉和MCM复合体等；分子功能主要涉及蛋白磷酸酶结合、病毒受体活性、外源性蛋白结合、错配修复复合物结合和磷酸酶结合等。KEGG富集分析结果显示富集通路分成8个聚类，值得注意的是，C1聚类中的破骨细胞分化通路(osteoclast differentiation)、C3聚类中的PI3K-AKT信号通路(PI3K-Akt signaling pathway)、C7聚类中的糖尿病并发症中的AGE-RAGE信号通路(AGE-RAGE signaling pathway in diabetic complications)和其他聚类中的铁死亡信号通路(Ferroptosis)。

图3 KEGG通路富集分析Fig.3 KEGG pathways enrichment of DEGs

2.3 差异基因PPI网络建立及分析

将所有差异基因名称导入STRING在线数据库中，按照默认设置并选择去除无连接节点，构建差异基因PPI网络(图4)。将PPI注释文件导入Cytoscape软件中，利用cytoHubba插件计算PPI网络中10种算法下各节点分值，利用R语言中的UpSetR工具包进行交叉验证得到11个关键基因(图5A)，然后提取11个关键基因相互作用网络(图5B)：节点颜色越深代表其在网络中的degree值越高。其中CDK16为下调差异基因，其他为上调差异基因。

图4 差异基因PPI网络分析Fig.4 PPI network analysis of DEGs

图5 关键基因筛选Fig.5 Identification of key genes

2.4 关键基因主成分分析

为进一步明确筛选得到的关键基因是否能够代表高低峰值骨量人群基因表达的差异，本研究进行了关键基因在高低峰值骨量人群间表达水平的PCA分析。PCA分析显示组间差异显著，说明关键基因可以将高低峰值骨量人群区分开(图6)。

图6 关键基因在高低峰值骨量人群间表达水平的PCA分析Fig.6 Results of principal component analysis

2.5 关键基因和骨质疏松症的相关性

为验证关键基因和骨质疏松症的相关性，通过下载和高低峰值骨量人群基因表达芯片为同一平台文件(GPL570)的数据集GSE35959，进行关键基因表达水平验证，探究关键基因在骨质疏松患者和健康者之间的表达差异。GSE35959包含了5例骨质疏松患者和14名健康者的基因表达数据。通过验证发现CCR5、CDK16、RBBP7和SRSF7在骨质疏松患者和健康者中表达水平差异显著(图7)。另外，主成分分析也显示这4个基因在骨质疏松患者和健康者组间差异显著，组内差异较小，能显著区分OP人群和健康人群(图8)，具有潜在的诊疗价值。

图7 验证关键基因与OP患者中的表达水平Fig.7 Validation of expression level of key genes in OP patients

图8 4个验证基因主成分分析Fig.8 Principal component analysis of 4 validated genes

3 讨论

随着高通量测序技术及生物信息学的发展，基因芯片分析已经成为识别疾病和健康状态间差异表达基因的有效方法，是识别疾病相关生物标志物、作用途径以及药物靶标的强大技术手段[15-16]。多项研究利用基因芯片分析确定了骨质疏松发病机制中的关键基因，有助于更准确地识别骨质疏松发病中基因生物标志物。Xia等[17]通过基因芯片分析发现VPS35、PHF20、NFKB2和HIRA在骨形成和骨代谢调控中发挥重要作用，可以作为潜在的诊断及治疗靶点。Qian等[18]通过WGCNA等生物信息学方法对基因芯片进行分析，确定PPWD1可以作为绝经后骨质疏松潜在诊断生物标志物。Zhang等[19]研究通过mRNA和lncRNA芯片分析，揭示了JUN和ACSL5可能是骨质疏松潜在生物标志物和临床治疗靶点，并构建了骨质疏松ceRNA调控网络。

骨质疏松症的发病往往呈隐匿性，早期阶段并不出现相关症状，以致不容易被发现及早期治疗。对于高低峰值骨量人群生物标志物的差异分析，并发掘其在骨质疏松症中的潜在诊疗价值具有重要意义。本研究充分发挥生物信息学优势，对GEO数据库中高低峰值骨量人群的基因表达芯片进行了分析，最终筛选得到182个差异基因，其中73个下调差异基因，109个上调差异基因。GO和KEGG富集分析显示，这些差异基因可以通过多个生物进程及信号通路发挥作用。值得注意的是，KEGG通路分析结果中C1聚类的破骨细胞分化通路、C3聚类的PI3K-AKT信号通路、C7聚类的糖尿病并发症中的AGE-RAGE信号通路和其他聚类的铁死亡信号通路和骨质疏松症发病机制密切相关。林院等[20]研究表明丹酚酸B能够通过激活PI3/AKT信号通路，上调MC3T3细胞ALP、OPN、COLⅠ和OCN表达，刺激MC3T3细胞骨形成。Li等[21]证实miR-363-3p可以通过激活PI3K/AKT信号通路，促进破骨细胞分化，抑制成骨细胞分化，促进骨质疏松症的发病。AGE-RAGE间相互作用可以影响细胞功能、运动和新陈代谢，与骨重塑过程密切相关[22]。Yi等[23]研究发现莫诺苷可以通过抑制AGE-RAGE信号通路减轻高糖介导的骨髓间充质干细胞功能障碍，并可能成为糖尿病骨质疏松症的潜在治疗手段。铁死亡是最近的研究热点[24]，越来越多的研究关注铁死亡途径在骨质疏松症防治中的作用机制。Yang等[25]的研究中通过体内和体外实验证实，内皮细胞分泌的外泌体可以通过抑制铁蛋白吞噬依赖性的铁死亡，逆转糖皮质激素引起的骨质疏松症。Ni等[26]全面研究了破骨细胞铁死亡的体内外作用，发现靶向调控HIF-1α和铁蛋白，诱导破骨细胞铁死亡可能是治疗骨质疏松症的一种替代疗法。

本研究通过对差异基因进行PPI网络分析，确定了11个在高低峰值骨量人群中差异表达的关键基因，PCA分析显示关键基因在高低峰值骨量两组组间差异显著，能显著区分高低峰值骨量人群。另外，通过进一步验证发现CCR5、CDK16、RBBP7和SRSF7在骨质疏松患者和健康者中表达水平差异显著；同时，PCA也显示这4个基因在骨质疏松患者和健康者组间表达差异显著，组内差异较小，具有潜在诊疗价值。CCR5可以通过调节破骨细胞功能，抑制骨吸收。Lee等[27]通过体外实验证实敲除CCR5表达，可抑制破骨细胞活性及骨吸收过程，抵抗RANKL诱导的骨质疏松症。CDK16是细胞周期依赖激酶家族中的一员，参与转录及细胞周期的调节[28]。RBBP7是组蛋白修饰和染色质重塑复合体成员之一，具有调节细胞周期和细胞增殖分化功能[29]。SRSF7是富含丝氨酸/精氨酸(SR)蛋白家族中的一员，调节各种mRNA前体的组成剪接和选择性剪接，在细胞凋亡中发挥重要作用[30]。目前为止，未有直接证据表明CDK16、RBBP7和SRSF7参与调控骨质疏松症的发病过程，是未来的潜在研究方向。

综上所述，本研究利用生物信息学方法鉴定了高低峰值骨量人群生物标志物的差异，并通过进一步验证发现差异标志物中CCR5、CDK16、RBBP7和SRSF7可能与骨质疏松症的发病密切相关。这4个差异表达基因有可能成为早期筛选骨质疏松症高风险人群的生物标志物，对骨质疏松症的防治发挥重要作用。与此同时，本研究存在一些局限性，例如：GEO芯片样本量不足，缺乏实验证据支持。但是，本研究利用生物信息学的优势，发掘了潜在的生物标志物，为后续进一步实验研究及临床治疗提供了有效依据。