骨质疏松诊断标志物的筛选及其与免疫浸润的关系

计 超, 陈一霄, 刘 凯, 苏啸尘, 滕梦豪, 郭 燕, 姬文晨, 李 萌

(1. 西安交通大学第一附属医院 骨科, 陕西 西安, 710061;2. 西安交通大学第二附属医院 生物诊断治疗国家地方联合工程研究中心, 陕西 西安, 710004;3. 西安交通大学 生命科学与技术学院生物医学信息与基因组中心, 陕西 西安, 710049)

骨质疏松(OP)是最常见的代谢性骨病之一[1-2], 临床特征包括骨密度和骨质量下降、骨微结构破坏和骨脆性增加,从而导致骨折易感性增加[3-4], 因此OP的早期诊断对于改善治疗和预后至关重要。目前, OP的常规诊断多依靠临床表现和影像学检查监测骨密度[5]。然而,这种测量方法的敏感性较差,可能会忽视潜在骨折的早期症状[6]。近年来,为OP的临床管理开发的生物标志物已显示出检测高骨折风险个体的敏感性和可靠性[7]。这些生物标志物对于预测OP风险、确定潜在治疗靶点和探索机制很有价值。另一方面,免疫系统在骨科疾病中的作用得到证实,能够促使“骨免疫学”的产生和发展[8]。OP是最常见的炎症性骨质流失病症之一。研究[9]表明,免疫细胞浸润在OP的发生、发展中起重要作用,如Th17细胞是T淋巴细胞破骨细胞的一个亚群,血液及外周组织中的Th17细胞可作为OP的重要标志物。CIBERSORT是一种软件工具,用于分析表征由多种细胞类型组成的免疫细胞亚群[10]。因此,从免疫系统的角度,评估免疫细胞的浸润情况,对于阐明OP的分子机制以及开发新的治疗靶点至关重要。本研究首先筛选得到OP患者的差异表达基因(DEGs), 然后使用机器学习筛选其潜在的诊断标记,并应用CIBERSORT分析OP患者和正常受试者的22个免疫细胞亚群差异,此外还分析了筛选得到的诊断标记与免疫浸润细胞的关系,以更好地了解OP发生、发展的分子及免疫机制。

1 材料与方法

1.1 数据获取

从基因表达综合数据库(GEO)( https://www.ncbi.nlm.nih.gov/geo/)中下载OP的表达谱数据集GSE56815、GSE13850、GSE35959和GSE7158。

1.2 数据预处理和差异基因筛选

R语言(版本3.6.2)和R软件包用于数据分析[11]。整合GSE56815和GSE13850基因数据集表达矩阵,使用“sva 3.34.0”包消除批次间差异。使用二维主成分分析(PCA)聚类来证明样本之间的校准效果。通过“limma 3.42.0”包过滤DEGs, 使用“ggplot2 3.3.0”包绘制DEGs火山图,展示DEGs的差异表达情况。设置P<0.05和|log2FC|>1为统计学显著性的标准。

1.3 富集分析

利用基因集富集分析(GSEA)平台对DEGs进行基因本体论(GO)和京都基因和基因组百科全书(KEGG)富集分析[12]。FDR<0.25和P<0.05是显著富集的标准。“cluster Profiler 3.14.3”包用于分析和绘制结果。

1.4 建模及筛选诊断标记物

最小绝对值收敛和选择算法(LASSO)逻辑回归和支持向量机递归特征消除(SVM-RFE)用于OP诊断标记的特征选择和筛选[13-14]。对GSE35959和GSE7158数据集的表达矩阵进行质量控制后,使用“sva 3.34.0”包去除批次间的差异。然后将数据集合并成一个独立的数据集,并基于独立的数据集验证所获得的诊断标记的组合诊断效率。“glmnet 4.0.2”包用于LASSO逻辑回归分析。通过“e1071 1.7.6”包运行SVM-RFE以进一步确定这些生物标志物在OP中的诊断价值。判断筛选结果的显著性阈值为P<0.05。最后,将LASSO逻辑回归和SVM-RFE算法识别出的基因结合起来进行后续分析。

1.5 免疫细胞浸润情况评估

CIBERSORT工具用于分析基因表达矩阵数据; 筛选P<0.05的样本,得到免疫细胞浸润矩阵。然后使用“ggplot2 3.3.0”包对免疫细胞浸润矩阵数据进行PCA聚类分析,绘制出二维PCA聚类图。“corrplot 0.88”包用于绘制相关热图,以可视化22种浸润免疫细胞之间的相关性。“ggplot2 3.3.0”包绘制小提琴图以直观反映免疫细胞浸润差异。

1.6 诊断标志物与免疫细胞浸润的相关性分析

“ggstatsplot 0.2.0”包用于分析诊断标志物和免疫浸润细胞之间的Pearson相关分析。“ggplot2 3.3.0”包用于可视化结果。

2 结 果

2.1 DEGs的筛选

首先,对GSE56815和GSE13850数据集的基因表达矩阵进行归一化和批量校正。数据集归一化前后的二维PCA聚类图如图1A和图1B所示。结果表明,归一化后的2组样本聚类更加明显。数据预处理后,使用R语言软件从基因表达矩阵中共提取331个DEGs, 见图1C。ZER1、RPL10、FOLR1、DVL1和EFCAB1等基因的表达水平显著上调;MAP3K3、IL1R2、RPL3L、SLC41A3、EGFR和WNT1等基因的表达水平显著下调(图1D)。

A: GSE56815和GSE13850数据集样本校正前PCA聚类图; B: GSE56815和GSE13850数据集样本校正后PCA聚类图(红色代表GSE56815数据集,蓝色代表GSE13850数据集); C: DEGs火山图(蓝色代表上调差异基因,绿色代表无显著差异基因,红色代表下调差异基因); D: 上调及下调的前50个DEGs的热图(红色代表上调的差异基因,绿色代表下调的差异基因)。

2.2 功能富集分析

GO富集分析中的细胞成分分析发现,大多数DEGs富集于“核仁”“内吞囊泡”“锚定点”和“内质网”相关通路(图2A)。此外,对分子功能通路的分析表明,已筛选得到的DEGs可能通过调节“RNA结合”“酶结合”和“翻译调节活性”相关通路来发挥其分子功能(图2B)。GO的生物过程通路的富集结果表明, DEGs主要与“含有复杂亚基组织的蛋白质”“细胞凋亡过程”“大分子分解代谢过程”和“生物合成过程的正调控”的调控有关(图2C)。KEGG通路富集分析的结果见图2D所示。富含DEGs的通路主要涉及子宫内膜癌、前列腺癌等癌症相关通路、胞吞相关通路和Wnt信号通路。以上结果表明, DEGs主要集中在激素调节相关通路(如“内吞小泡”“内吞参与通路”“子宫内膜癌等癌症相关通路”)和免疫反应过程(如“内质网”“翻译调节剂的活性”和“生物合成的正向调节”),提示其在OP的发病机制中起着重要作用。

A: OP的DEGs GO富集分析的前9个细胞成分分析; B: OP的DEGs GO富集分析前9个分子功能分析;C: OP的DEGs GO富集分析前9个生物学过程分析; D: KEGG富集分析OP的 DEGs的前9条富集通路。

2.3 诊断标志物的筛选和验证

LASSO逻辑回归和SVM-RFE算法用于从DEGs中筛选基因作为OP的诊断标志物。2种算法得到的基因标志物重叠,最终筛选出5个诊断相关基因为关键基因:SKAP2、SLC30A3、TDRD12、RPL10和CSPP1。为了进一步探索这些关键基因的诊断效能,将GSE35959和GSE7158数据集合并用作验证集(图3A、3B)。当SKAP2、SLC30A3、TDRD12、RPL10、CSPP1作为变量输入时,验证集的诊断效率达到了较高水平[曲线下面积(AUC)为0.786](图3C), 提示由SKAP2、SLC30A3、TDRD12、RPL10联合CSPP1构成的预测模型具有较高的诊断价值。

A: 合并验证集的GSE35959和GSE7158数据集前的PCA聚类图; B: 合并验证集的GSE35959和GSE7158数据集后的PCA聚类图; C: 5个诊断标志物拟合后验证诊断效果的ROC曲线(深蓝色曲线的AUC用于衡量模型性能, AUC越接近1, 表示模型性能越好)。

2.4 免疫细胞浸润分析

22种免疫细胞的相关热图显示,活化的肥大细胞与浆细胞呈正相关,静息肥大细胞也与嗜酸性粒细胞呈正相关。此外,活化的CD4记忆T细胞与调节性T细胞呈负相关,巨噬细胞与记忆B细胞呈负相关(图4A)。免疫细胞浸润差异分析显示, OP患者和正常对照者的M2巨噬细胞和静息CD4记忆T细胞差异有统计学意义(P<0.05), 见图4B。

A: 22种免疫细胞的相关热图(彩色方块大小表示相关性的强度,红色表示正相关,蓝色表示负相关,颜色越深表示相关性越强); B: 22种免疫细胞占比的小提琴图。

2.5 关键基因与浸润免疫细胞的相关性分析

进一步研究已筛选的OP诊断分子标记与免疫细胞浸润之间的相关性(图5A)。结果发现,SKAP2与静息CD4记忆T细胞呈正相关(r=0.238,P=0.009), 见图5B;CSPP1与幼稚B淋巴细胞呈正相关(r=0.198,P=0.030)(图5C);RPL10与嗜酸性粒细胞呈负相关(r=-0.179,P=0.050)(图5D)。

A: 5个关键基因与浸润免疫细胞之间相关性的热图(红色代表正相关,蓝色代表负相关); B: SKAP2与免疫浸润细胞之间的相关性;C: CSPP1与免疫浸润细胞之间的相关性; D: RPL10与免疫浸润细胞之间的相关性(圆点大小代表基因与免疫细胞相关性的强弱。圆点越大则相关性越强; 圆点越小则相关性越弱。圆点颜色深浅代表P值的大小。颜色越蓝则P值越小; 颜色越绿则P值越大)。

3 讨 论

OP是一种慢性骨科疾病,其特征是骨密度和骨质量下降,骨微结构破坏,导致骨脆性增加[15]。筛选OP生物标志物有助于改进骨折风险的评估方法并改善其早期诊断效果。既往研究[8]表明,免疫细胞浸润在OP的发生、发展过程中具有重要作用。因此,特异性诊断标志物的筛选和免疫细胞浸润的分析对于改善OP患者的预后具有深远意义。此外, CIRBESORT工具为分析疾病免疫细胞浸润模式提供了便利。本研究确定了OP的早期诊断标志物,并进一步探讨了免疫细胞浸润在OP中的作用,通过差异表达分析共鉴定出331个DEGs。有研究[16]使用这些OP细胞群进行生物信息学分析, GO富集分析显示, DEGs主要与内质网活性、翻译调节因子和生物合成过程的正调控有关,提示免疫反应在OP的发病机制中具有重要意义。研究[17]表明,雌激素通过抑制破骨细胞的活性在OP的发生过程中发挥重要作用, Wnt信号通路可以控制成骨细胞的形成[18]。

本研究结合LASSO逻辑回归和SVM-RFE算法确定SKAP2、SLC30A3、TDRD12、RPL10和CSPP1作为OP的诊断标志物。研究[19]表明,SKAP2基因编码src激酶相关磷蛋白2参与不同生理过程,包括整合素信号转导、细胞迁移和癌症进展。研究[20]发现,SKAP2在类风湿性关节炎和1型糖尿病中的表达降低。本研究结果表明, OP中SKAP2表达降低,可能与调节内分泌激素有关。SLC30A3的编码产物是锌转运蛋白3, 主要参与突触小泡中锌的积累。本研究发现, OP患者SLC30A3的表达显著上调,这与雌激素调节对OP患者的作用一致。此外,研究[21]发现,RPL10在脊椎动物软骨内骨发育过程中呈差异表达。本研究发现,RPL10在OP患者中表达上调,推测RPL10通过调节内分泌激素影响其骨骼发育。此外,对TDRD12基因的研究尚未得到统一结论。研究[22]发现,TDRD12参与JAK-STAT通路,在脂质代谢调控中发挥作用,而脂质代谢紊乱与骨吸收和骨形成失衡密切相关,可能是绝经后OP发生的基础。LI Q H等[23]发现,CSPP1衍生的circ-CSPP1促进卵巢癌细胞的增殖和迁移,可以推测OP患者中下调的CSPP1也可能受到内分泌激素的影响。因此,SKAP2、SLC30A3、TDRD12、RPL10和CSPP1可能参与OP的发生、发展,可作为其诊断标志物。

为了进一步探讨免疫浸润细胞在OP中的作用,采用CIBERSORT对OP的免疫浸润进行综合评估。结果表明,与对照组相比,OP中M2巨噬细胞浸润增加,静息CD4记忆T细胞浸润减少。研究[24]发现,活化的T淋巴细胞是RANKL和肿瘤坏死因子-α(TNF-α)的主要来源,可导致各种病理和炎症条件下的骨破坏。由此提示, CD4记忆T细胞和巨噬细胞在OP的发病机制中起重要作用,应成为进一步研究重点。此外,关键基因与免疫细胞相关性分析显示,SKAP2与静息CD4记忆T细胞呈正相关,与M2巨噬细胞呈负相关;CSPP1与幼稚B淋巴细胞呈正相关,与巨噬细胞呈负相关;RPL10与静息NK细胞呈正相关,与嗜酸性粒细胞呈负相关;SLC30A3与CD8+T细胞呈正相关,与活化的CD4记忆T细胞和活化的NK细胞呈负相关;TDRD12与M0巨噬细胞呈正相关,与中性粒细胞呈负相关。上述结果表明,本研究筛选的分子标记基因可能在增加CD4+T细胞和M2巨噬细胞或减少CD8+T细胞方面发挥作用,参与OP的发生进展,以上结果还需要在后续研究中进一步阐明。

综上所述,本研究使用LASSO逻辑回归和SVM-RFE确定了SKAP2、SLC30A3、TDRD12、RPL10和CSPP1作为OP的诊断标志物。此外, CIBERSORT被用于分析OP相关的免疫细胞浸润,结果表明,静息CD4+T细胞和M2巨噬细胞等免疫浸润细胞可能参与OP的发生和发展。此外,本研究还发现,SKAP2与静息CD4记忆T细胞和巨噬细胞相关,SLC30A3与激活的CD4记忆T细胞和激活的NK细胞相关,TDRD12、CSPP1与巨噬细胞相关,RPL10与静息NK细胞相关。这些免疫细胞可能在OP的发生、发展中起着关键作用,对这些免疫细胞的进一步探索可能会证实其是OP免疫治疗的靶点。本研究筛选的OP相关诊断分子标志物有助于为OP的早期诊疗提供新的思路和方向,为药物靶点提供新的可能性。