文献解读|Cell Discov（33.5）：临床相关骨质疏松亚型的多模式分子决定因素

由于全球人口迅速老龄化，骨质疏松症和相关的骨折风险已成为一个广泛传播的公共卫生问题。然而，骨质疏松症具有很大的异质性，现有的标准诊断措施不足以准确识别所有有骨质疏松性骨折风险的患者并指导治疗。

研究团队利用来自中国社区骨质疏松症队列 (CCCO) 的 366 个样本，构建了迄今为止最大的骨质疏松症队列的第一个多组学图谱（图1a-b）。这些生物大数据包括用于基于遗传和环境因素全面表征骨质疏松风险的基因组学、表观基因组学和宏基因组学数据（图1a-b）。他们首先通过高效的深度潜在空间融合(DLSF)分析，基于多组学景观和相应的多模式分子特征(M3S)，识别了中国人群中的临床干预相关亚型(CIS)，并且CIS在临床过程中表现出不同的敏感性。

图1. 基于数据密集型分析生成的多组学图谱研究骨质疏松症及其亚型的设计和工作流程。

(a) 以数据密集型分析为原则的研究设计和工作流程。 (b) 骨质疏松症分子景观的整体视图。

研究团队利用来自中国社区骨质疏松症队列 (CCCO) 的 366 个样本，构建了迄今为止最大的骨质疏松症队列的第一个多组学图谱（图1a-b）。这些生物大数据包括用于基于遗传和环境因素全面表征骨质疏松风险的基因组学、表观基因组学和宏基因组学数据（图1a-b）。他们首先通过高效的深度潜在空间融合(DLSF)分析，基于多组学景观和相应的多模式分子特征(M3S)，识别了中国人群中的临床干预相关亚型(CIS)，并且CIS在临床过程中表现出不同的敏感性。

图2. 骨质疏松症的多组学分子景观表征。

(a) PLSDA 空间中的样本组分布。(b) 不同组学数据级别的区分特征的百分比。(c) 不同组学数据级别的表型相关特征的百分比。 (d) 不同样本组在相应组学数据水平上的特征共表达网络。(e) 特定组学数据水平上不同样本组的特征及其表型关联的共表达网络模块。(f) 通过 RV 指数估计的不同组学数据集之间的共有信息。 (g) 基于联合自动编码器和自我表达技术的深层潜在空间融合模型，用于骨质疏松症多组学图谱的多组学无监督学习。

与其他传统的多组学分析方法相比，DLSF 使研究者能够识别与更多临床表型相关的样本聚类，包括相似网络融合 (SNF) 和用于数据集成和疾病亚型分析的扰动聚类 (PINS) （图3a），DLSF 是分析骨质疏松症多组学分子景观的有效方法（图3b）。

根据分子表型确定的这种亚型，骨质疏松个体在统一的潜在数据空间中可以清楚地分为两组（图3c），并且相同的样本在不同的个体观察数据空间中表现出相似的组，特别是在微生物群和甲基化模式中。这两个样本聚类/亚型与许多临床表型有很强的相关性（图3d）。例如，他们显示了两个亚型之间的许多表型分布差异。值得注意的是，与亚型 II 相比，亚型 I 内的个体往往具有较低的总髋骨矿物质密度 (TH BMD) 和其他相关指标，例如较低的股骨颈骨矿物质密度 (FN BMD)。同时，亚型 II 内的个体具有较低的腰椎骨矿物质密度 (L1–L4) (L1–L4 BMD) 和其他指数（包括较低的 L1–L4 BMD Z分数）。

两种亚型的骨转换标记物(BTM) 也表现出差异性变化（图3e-f）。与正常组和骨质减少组相比，I亚型I型前胶原N端前肽(PlNP)和骨钙素(OST)有显著变化；同时，亚型II与正常组相比PlNP和OST有显著变化。值得注意的是，I亚型的碱性磷酸酶(ALP)存在一定差异，两种亚型的含有交联C端肽(β-CTX)的I型胶原的β-CrossLaps的趋势发生了变化，分布也有所不同。这些结果表明，亚型 I 和 II 均表现出 BTM 功能障碍。特别是，通过对 CCCO 队列中同一个体的随访骨矿物质密度(BMD) 结果来评估这些骨质疏松症亚型对治疗的反应。与亚型II相比，亚型I的许多BMD临床指标都有所改善（图3g）。根据病史调查结果，观察到接受钙补充剂的个体百分比在骨质疏松症和其他患者之间存在显著差异，而不是在骨质疏松症亚型之间存在显著差异（图3h）。这表明骨质疏松症个体可以通过补充必要的钙来改善，不同的骨质疏松症亚型在改善过程中具有不同的敏感性，例如CIS1（亚型I）是敏感组，而CIS2（亚型II）往往是低敏感或不敏感组。

这些骨质疏松症亚型通过DLSF有效地鉴定出来，DLSF是一种发现符合生物学中心法则的分子亚型的上游方法（图3i）。因此，下一步的下游程序是通过反向生物学解释来表征亚型（即CIS）（图3i），主要包括差异关联矩阵分析(DAM)表征表型相关的分子特征和多组学定量性状位点(xQTL)表征表型相关的基因型，并在骨质疏松症诊断、亚型预后、骨折预测等方面实现临床应用（图3i）。

图3. 通过深层潜在空间融合模型识别的骨质疏松症的多模态亚型。

(a) 基于表型关联数量的 DLSF 和其他模型的效率评估。(b) DLSF 和其他模型的稳健性评估。 (c) 不同组学数据水平显示的综合骨质疏松症亚型及其判别空间。(d) 与骨质疏松症相关的几种表型和两种已确定的亚型（I 型和 II 型）。(e-f) 分别表示正常、骨质减少和骨质疏松亚型 I 和亚型 II 之间的 BTM 差异。 (g) 亚型 I 和亚型 II 的两年随访评估。 (h) 骨质疏松病史调查，表明接受特定医疗干预（例如钙补充剂）的个体百分比，以及骨质疏松症与其他疾病之间或两种骨质疏松症亚型之间的分布差异。 (i) 骨质疏松症及其亚型的分析流程。

为了了解已识别的骨质疏松症亚型的分子特征，研究了每种模式的 M3S 分子特征。首先，提取了 M3S 的微生物群和代谢物特征，因为微生物群会发生相互作用并对宿主代谢产生功能影响。微生物群在亚型之间表现出比代谢物更大的区分（图4a），这表明其在分析中占主导地位。微生物群特征与不同亚型的许多代谢物特征具有关联（图4b）。例如，小韦荣球菌和3-羟基丁酸之间的相关性在CIS1中为正相关，而在CIS2中为负相关（图4b）。

他们还通过 DAM 发现Klebsiella quasipneumoniae和苏糖酸之间的关联在 CIS1 中呈阴性，但在 CIS2 中呈阳性（图4b）。肺炎克雷伯菌与Klebsiella quasipneumoniae具有相似的代谢关联，并且可以诱导炎症免疫细胞的产生，这些细胞可以在遗传易感宿主的环境中迁移到远处的炎症组织，并参与非肺炎的发生和进展。此外，3-羟基丁酸浓度在1型糖尿病和糖尿病昏迷患者中升高，苏氨酸参与钙吸收，表明它们在骨质疏松症和亚型中的潜在区别作用和功能相关性。

这些特征在骨质疏松症及其亚型中的共同和特殊作用也得到了它们与骨质疏松症相关功能的关联的验证（图4c）。一般来说，微生物群特征比代谢物特征显示出更多的功能关联，表明微生物和宿主引起的常见和特定的代谢功能障碍。他们还发现这些特征在 CIS1 中比在 CIS2 中具有更多的功能关联。代谢特征分析揭示了额外的病原体信息（图4d）。通过富集分析，这些代谢物主要属于β羟基酸，参与脂肪酸生物合成和丙酸代谢。与正常组和骨质减少组相比，Klebsiella quasipneumoniae和小韦荣球菌在骨质疏松症CIS1中有特异性表达，3-羟基丁酸在骨质疏松症CIS1中有上调的趋势，而苏氨酸在骨质疏松症CIS2中有上调的趋势（图4e）。

图4. 与骨质疏松症及其亚型相关的微生物群和代谢物特征。

(a) 两种骨质疏松症亚型期间微生物群和代谢物特征的表达谱。 (b) 两种骨质疏松亚型期间微生物群和代谢物特征之间的关联矩阵。(c) 微生物群/代谢物特征以及注释的功能和途径之间的亚型特异性关联矩阵。 (d) 代谢物特征的功能富集。 (e) 几种具有代表性的微生物群和代谢物特征，表明骨质疏松症及其亚型。

在多组学图谱的甲基化模式中，表观基因组关联研究(EWAS)分析证实了表观遗传信息对不同临床表型的贡献（图5a）。与年龄或身体质量指数(BMI)相比，甲基化谱的不同主成分(PC)往往对基因型-表型关联有很多影响，并且不同的PC表现出不同的偏好。例如，与所有表型的平均值相比，CpG PC2和CpG PC3与FN BMD具有更大的相关性，而CpG PC7与L1-L4 BMD具有更大的相关性（图5b）。这些发现表明甲基化谱能够提供常见和特定的信息和特征来定义骨质疏松症及其亚型。

他们筛选出了一组与骨质疏松症及其亚型相关的CpG位点（图5a）。这些甲基化特征与上述代谢物和微生物群特征具有相应的关联（图5c-d），表明它们与骨质疏松症期间代谢改变的潜在相关性。例如，cg18008345在CIS1中与唾液链球菌和副血链球菌呈正相关，但在CIS2中相关性降低，而cg09057517在CIS1中与Gemmiger formicilis呈负相关，在CIS2中与副流感嗜血杆菌呈负相关（图5c）。此外，cg18243574在CIS1中显示出与微生物群特征的许多关联，包括肺炎克雷伯菌，而cg18249968在CIS2中与齿双歧杆菌具有显著的负相关性（图5c）。这些结果表明，宿主表观遗传信息可以解释在外部因素（如微生物群）中观察到的许多常见和特定特征。类似地，这些甲基化特征显示出与代谢物特征的亚型特异性相关性，包括己二酸、5Z-十二碳烯酸和9E-十四碳烯酸，表明表观遗传和代谢组水平上的互补信息和表征（图5d）。值得注意的是，cg18249968 在 CIS1 中对大多数代谢特征呈阴性，但在 CIS2 中对大多数代谢特征呈阳性。

这些跨模式关联表明甲基化特征也必须与临床表型具有共同和特定的关联（图5e）。他们发现了一些共同的关联，包括cg18008345与不同亚型的L1-L4 BMD指数显著负相关，cg23088318与其正相关，而cg21813369与两个亚型的许多指数呈负相关（图5e） 。同时，他们还观察到亚型特异性关联，包括CIS2中的cg09057517与吸烟史呈正相关，CIS1中的cg18243574与糖尿病家族史呈负相关，与中药(TCM)用药史呈正相关。

值得注意的是，在包括EWASdb在内的公共EWAS研究中可以检测到这些甲基化特征，并且许多特征位于DNA序列上的甲基化调节基因，这些基因与骨质疏松症具有功能关系（图5f）。例如，cg18008345是POLR2E中的位点。该基因在破骨细胞发生过程中表达上调，其表达受到抑制会抑制破骨细胞发生，而其过表达在体外表现出完全相反的作用，其基因产物与CREB1相互作用调节破骨细胞骨吸收。cg21813369位于CSN1S2B位点，在2型糖尿病大鼠骨骼中，CSN1S2蛋白通过骨形态测量蛋白信号传导显示出对骨质疏松的保护作用。

他们还过滤掉了许多与骨质疏松症相关的表型，以连接目标位点和基因（图5g）。例如，以下基因中的 CpG 变异与许多骨质疏松症相关表型有关：CSN1S2B和GNB5与肥胖、精神分裂症和骨肉瘤相关，PARK2与已知骨质疏松影响因素的克罗恩病相关，RUNX1与髓母细胞瘤相关，包括最常见的神经外转移部位包括骨骼和骨髓。鉴于这些功能相关性和结果，他们证实了M3S甲基化特征的生物学和生物医学意义，其表达在已识别的骨质疏松症亚型中表现出显著的下调变化（图5h）。与正常组和骨质减少组相比，cg21813369 (CSN1S2B) 在骨质疏松 CIS1 中的表达较低。相反，cg18243574 (CLEC3A) 在骨质疏松症 CIS2 中表达较低。这些特征代表了新的骨质疏松症亚型在表观遗传水平上的分子特征。

图5. 与骨质疏松症及其亚型相关的甲基化特征。

(a) 甲基化特征的CpG。(b) EWAS 结果摘要。 (c) 甲基化和微生物特征之间的亚型特异性关联矩阵。(d) 甲基化和代谢特征之间的亚型特异性关联矩阵。(e) 甲基化特征和临床表型之间的亚型特异性关联矩阵。(f) 有关 EWAS DB 中甲基化特征的公开数据。(g) 公开报道的与甲基化特征相关的表型。(h) 几个代表性的甲基化特征，表明骨质疏松症及其亚型。

他们总结了正常组、骨质减少组和骨质疏松组中SNP/SNV分布的差异（图6a），发现两种骨质疏松亚型的整体SNP发生显著变化，表明骨质疏松症及其亚型中存在基因型变化。位于snpGene INPP4B 的多个 SNP在 CIS1 中往往具有一致的变化，这对于ATP6V0A4和GLP1R来说是相似的。位于另一个snpGene PTPRD的几个SNP在CIS2中具有一致的变化，其中PDGFA和SCARA3具有类似的变化。INPP4B 和ATP6V0A4 是已知的破骨细胞分化调节剂，并且具有人类骨质疏松症的预后位点，而已知SCARA3可以调节脂肪细胞和成骨细胞分化之间的转换，因此代表了骨质流失和骨质疏松症的潜在治疗靶点。这些结果表明 SNP 对分子表型的不同因果关系。大多数SNP基因显示出密集的功能蛋白相互作用或关联（图6b），并且在钙信号通路和Rap1信号通路中富集（图6c）。钙信号参与破骨细胞的分化和功能，其功能障碍可导致破骨细胞骨吸收异常，从而导致骨质疏松。Rap1 在破骨细胞的再吸收功能中发挥着关键作用，其在成熟破骨细胞中的选择性抑制可延缓病理性骨质流失。值得注意的是，在比较毒物基因组学数据库中报道了许多整体snp基因（或作为CTD基因），它们在许多骨质疏松相关的推断网络中具有致病作用，包括雌二醇和白藜芦醇（图6c）。这些发现与之前的研究结果相吻合，即总雌二醇缺乏的老年男性患骨质疏松症，而骨质疏松症患者更有可能缺乏总雌二醇。

另一方面，他们研究了两种骨质疏松亚型对应的xQTL模型，以推断基因型与分子表型（如M3S）之间的定量关联。从QTL SNP变异型的分布来看，再次观察到了这两种亚型的共同和特定的分子特征（图6d）。对于具有定量关联的 M3S 的基因型特征，变异的主要分布是基因间和内含子，这对于不同类型的 QTL（包括甲基 QTL、代谢物 QTL 和微生物群 QTL）是一致的。相比之下，还观察到了SNP变异类型的具体特征，与CIS2相比，CIS1在微生物群QTL中具有更多的内含子变异，在代谢物QTL中具有更少的内含子变异。事实上，这种涉及xQTL的局部snp基因与上述整体snp基因有很多重叠（图6b）。相比之下，CIS1在甲基QTL和全局QTL之间具有最相似的snp基因，而CIS2在微生物群QTL和整体QTL之间具有最相似的snp基因（图6e）。这些因素再次支证明了骨质疏松症多组学分析和多模态模型的必要性和重要性。

CIS1与参与雌激素和AMPK信号通路的基因有相对较多的xQTL关联，是骨质疏松症个体中的一组敏感应答者，注解为“阳”型，代表临床干预期间骨密度有较大改善的预期。相反，CIS2与大量钙信号通路相关基因存在基因型-分子表型关联，在骨质疏松症个体中属于非敏感应答者，并注释为“阴”型，表明BMD改善预期较低（图6f）。

与 CIS1 相比，CIS2 在这些基因上往往具有更高的 SNV 等位基因频率，以及一些新的相关基因，如钙敏感受体(CASR)（图6g），表明候选基因CIS2 和 CIS1 之间补钙敏感性差异的遗传决定因素。此外，他们还发现了这种敏感性差异的许多候选表观遗传决定因素，并且钙相关基因的许多CpG在CIS2中表现出高甲基化（图6h）。例如，与 CIS1 相比，ESR1、VDR、CASR 和成骨细胞相关基因 BMPR1B 在 CIS2 中具有高甲基化。相比之下，破骨细胞相关基因 NFATC1 和 HDAC4可能在 CIS1 中过度甲基化（图6h）。

图6. 骨质疏松症亚型多模态分子特征M3S的基因组特征。

(a) 骨质疏松症亚型的整体 SNV 模式。(b) 整体snpGenes之间的蛋白质关联网络。(c) 整体snpGenes的功能富集。(d) 不同骨质疏松症亚型对应的 xQTL SNP 的 SNV 类型偏好。(e) xQTL 模型中不同 QTL 之间共有和特定的局部 snpGene。 (f) “阴-阳”骨质疏松亚型，其特征在于亚型特异性基因型-分子表型关联。 (g) 代表性钙相关基因的等位基因频率。(h) CIS 之间的甲基化差异，以少数钙和骨质疏松相关基因为例。

他们对数百个正常和骨质疏松样本构建了三种逻辑回归(LR)风险模型。第一个是基线 I，它基于年龄、体重指数以及饮酒和吸烟史。第二个是基线 II，即基线 I + BTM（即 β-CTX、OST 和 PINP）。第三个是本项研究的新模型，它是Baseline II + M3S。根据分析数据的决策曲线，该新模型在诊断骨质疏松症方面优于其他两个基线模型（图7a），表明M3S引入的临床辨别能力得到了改善。由于女性和男性之间的骨质疏松症风险可能有所不同，因此，性别也可以包含在相关风险模型中。

接下来，他们在一个独立（金山）队列中验证了 M3S、BTM 及其关联，该队列包含数百个正常和骨质疏松样本，其中 87 个样本与多组学数据相匹配。他们发现M3S在新队列中确实具有一致的关联，特别是对于微生物群和代谢物模式（图7b），这再次表明骨质疏松症的分子异质性。因此，在校正组学数据中存在的潜在批次效应后，基于这两种模式的风险模型再次分析，并且它仍然取得了比两个基线模型更高的性能（图7c）。

另一方面，另一种基于M3S的亚型评估风险模型可以有效识别独立金山队列中潜在的亚型个体（图7d）。M3S 模型预测的个体亚型存在许多表型关联（例如，通过相关表型的数量进行评估）（图7d）。

BMD 与骨折结果相关，BMD 越低，骨折的风险就越高。当追踪2023年收集的队列骨折事件时，发现骨质疏松组的骨折比例实际上高于正常组（图7e），支持BMD与骨折之间的因果关系。同时，在骨质疏松亚型中，CIS1 个体的骨折率低于 CIS2 个体。临床上普遍认为，骨折的人再次骨折的可能性较高。与其他情况相比，没有 PF（发生过骨折）的 CIS1 个体的骨折率最低；相反，无论是否考虑PF，CIS2个体的骨折率都大于10%（图7f-g）。

通过对指标值进行排序，可以评估用于预测骨折风险的生物标志物。以FN BMD为例（图7h），不同个体的值可以按升序排列（例如灰点），并且可以标记骨折事件（例如红点）。显然，在整个人群（即ALL组，包括正常、骨质减少和骨质疏松个体）的统计中，骨密度较低的个体确实有发生骨折事件的可能性。相比之下，在OP人群（即由骨质疏松症患者组成的OP组）中，骨折事件在BMD值之间的分布是分散的，这意味着用骨密度来预测骨折风险的效果可能不佳，这可能会成为临床上的一个新问题。（图7h-i）。

M3S实际上可以提供有效的骨质疏松症特异性骨折风险生物标志物（图7j）。例如，cg23088318 (AUC > 0.65) 和 cg155008935 (AUC ~ 0.7) 等甲基化位点是骨质疏松症特异性的。同时，无论骨密度如何，齿双歧杆菌和克雷伯菌都将成为 ALL 和 OP 人群的潜在生物标志物。鉴于 M3S 中的甲基化位点对骨折风险显示出显著的预测能力，他们筛选了 M3S 中一小部分 CpG 的整个组合，并发现了一些候选复合指数，通过成员 CpG 丰度的简单总和进行量化。例如，cg15320980 + cg18249968 + cg15508935 + cg23088318的组合可以在发现队列(CCCO)中实现0.73的AUC，甚至在独立金山队列中实现0.7的验证AUC（图7k）。

图7. M3S 的临床应用，包括骨质疏松症风险模型和骨折风险评估。

(a) 不同风险模型（骨质疏松症与其他）的决策曲线。 (b) CCCO 和金山队列之间表型关联的独立评估。 (c) 与两种基线方法相比，基于 M3S 的风险模型（骨质疏松症与其他）的独立性能评估。 (d) 与其他组件模型和基线方法相比，基于 M3S 的风险模型（阳 vs 阴子类型）的独立效应评估。 (e) CCCO 队列 4 年随访时骨折结果的统计数据。(f) 不同骨质疏松症亚型的骨折率。 (g) 通过既往骨折 (PF) 调整不同骨质疏松症亚型的骨折率。(h) 通过AUC排序分析进行骨折风险评估的病例。 (i) 骨折风险传统临床指标的 AUC。(j) 分别是 ALL 和 OP 群体中每个 M3S 特征的原始 AUC 和调整后的 AUC。(k) 分别根据发现队列(CCCO)和独立验证队列（金山）的原始和调整 AUC 预测骨折风险评估的综合指数。

本项研究构建了迄今为止最大的骨质疏松症队列的多组学图谱（来自 366 名参与者在三个时间点的纵向数据），并实施了一个可解释的数据密集型分析模型，通过DLSF，本项研究在中国个体中鉴定出两种具有相应分子表型的骨质疏松人群亚型，即CIS，其中骨密度有利于2年随访样本中对钙补充剂的反应。许多与这些分子表型相关的 snpGenes 揭示了骨质疏松症潜在的多种候选生物学机制，其中骨质疏松症及其亚型的 xQTL 偏好表明对不同生物学域的全基因效应。最后，根据4年的随访数据，发现这两种亚型与既往骨折具有不同的相关性和不同的骨折风险。因此，在临床应用中，M3S 可以帮助研究者进一步开发改进质疏松症诊断和治疗策略，并确定新的骨折预测综合指数，该指数在独立队列（166 名受试者）中得到了验证。

 + + + + +