《Genetics Research》：Single-Cell Transcriptomics and Integrated Bioinformatic Analysis Reveal Critical Biomarkers and Immune Infiltration Characteristics in Osteoarthritis
编辑推荐：
本文通过单细胞RNA测序（scRNA-seq）和整合生物信息学分析，系统揭示了骨关节炎（OA）软骨组织的细胞异质性、差异表达基因（DEGs）及免疫浸润特征。研究鉴定出八个软骨细胞亚群，并发现TNF、TGF-β和PI3K–Akt等信号通路显著富集。通过筛选公共数据集（如GSE114007），确定了2247个DEGs，并与ChEMBL数据库中的药物靶点基因取交集，最终锁定26个关键OA相关基因（如NR4A2、BMP1、AVPR1A）。ROC曲线验证了NR4A2、BMP1和AVPR1A在多个队列中具有卓越的诊断价值（AUC > 0.70）。分子对接显示bexarotene与上述靶点具有良好结合亲和力，提示其作为潜在治疗候选物的前景。此外，免疫浸润分析揭示了关键基因与多种免疫细胞（如T细胞亚群、巨噬细胞）的显著相关性，强调了OA炎症微环境的重要性。功能富集分析进一步表明，ECM重构、细胞增殖及HIF-1等通路在OA进展中发挥核心作用。本研究为OA的精准诊断和治疗靶点开发提供了新的分子视角。
单细胞RNA测序揭示膝关节软骨组织中的独特细胞亚群
为了解析骨关节炎（OA）与非OA（正常）个体膝关节软骨的细胞异质性，研究对软骨样本（GSE220243）进行了单细胞RNA测序（scRNA-seq）。经过严格的质量控制（QC），排除了表达基因数少于50或基因在少于3个细胞中表达的细胞。使用Seurat包中的PercentageFeatureSet函数计算线粒体和核糖体基因含量。共有11个样本通过QC，获得高质量转录组数据集。通过FindVariableFeatures函数识别高变基因，并对数据进行归一化、缩放和主成分分析（PCA）以降维。
随后使用FindNeighbors和FindClusters函数（维度=15，分辨率=0.2）进行聚类分析，揭示了七个主要簇和一个额外的次要细胞群体。通过SingleR辅助进行细胞类型注释，将主要簇分类为软骨细胞。研究进一步根据转录特征细分了这些软骨细胞簇，产生了八个亚群（Chondrocytes_1至Chondrocytes_8）。值得注意的是，正常组和OA组对每个亚群均有贡献，但比例不同。通过FindAllMarkers函数进行的标记基因分析证实了每个软骨细胞亚群具有独特的转录谱，并通过代表性标记基因的亚群特异性表达模式进行说明。
接下来，使用ClusterProfiler包对每个亚群进行富集分析。KEGG通路过度表征分析显示，七个软骨细胞亚群在包括“TNF信号通路”、“TGF-β信号通路”和“PI3K-Akt信号通路”在内的通路中表现出显著的功能富集。相比之下，一个富含成纤维细胞的簇（"Fibroblasts_4"）在统计阈值下未显示显著富集。这些结果强调了软骨细胞亚群之间的功能多样性，突出了OA发病机制中可能涉及的特定信号级联。
总的来说，单细胞转录组分析揭示了膝关节软骨内八个转录不同的亚群，正常和OA组织在亚群组成和基因表达程序上存在显著差异。这项工作为OA受累软骨的细胞景观提供了关键见解，并揭示了疾病进展的潜在通路。
骨关节炎关键差异表达基因的鉴定
为了阐明OA的转录改变，研究在使用limma包中的normalizeBetweenArrays函数对GSE114007数据集进行归一化后重新进行了分析。共检测到2247个DEGs，包括1311个上调基因和936个下调基因。通过火山图进行视觉检查，确认了OA软骨中几个基因的表达发生显著变化，例如COL1A1、POSTN和IL11在OA软骨中上调最为显著，而HILPDA、ADM和ATF3则是OA软骨中下调最显著的基因。热图强调了OA和对照样本之间不同的转录特征，OA组织与正常组织分开聚类，表明病变软骨中存在深刻的分子差异。
考虑到DEGs的治疗意义，研究接下来从ChEMBL数据库中收集了205个推定的OA治疗靶点基因，并将其与我们的DEG集取交集，产生了26个被指定为关键OA相关药物靶点的基因，包括PPARG、NCOR2、NR4A2、PRKAG2、RARB、RARA、BMP1、GDF15、RPS6KA4、SSTR5、BCL6、AVPR1A、MYOC、THRA、NR4A1、CRHR1、NDUFB3、NDUFB1、NDUFB2、NDUFA12、KCNH2、SGK1、GNRHR、HUNK、PSMA4和NR4A3。随后，使用并行scRNA-seq数据集（特别是GSE220243）中识别出的标记基因，以归一化的GSE114007数据为参考，进行了基因集变异分析（GSVA）。对Chondrocytes_1亚群GSVA得分的箱线图分析表明，OA组织相对于正常组织显著降低（p = 0.0023），表明该簇的转录特征在OA中发生了深刻改变。
最后，评估了这26个关键DEGs与所有九个软骨细胞亚群（Chondrocytes_1-8和一个额外的次要簇）的GSVA得分之间的相关性。Pearson相关分析揭示了特定亚群GSVA得分与26个候选药物靶点基因之间存在显著的正相关和负相关（p < 0.05 或 p < 0.01），如相关热图所示。这些发现突出了OA发病机制中推定的分子驱动因素和潜在治疗靶点，值得进一步的功能验证。
关键差异表达基因的验证
为了评估我们整合分析中确定的26个关键OA相关基因的诊断效用，我们在三个独立的验证队列（GSE114007、GSE82107和GSE98918）中检查了它们的表达谱和相应的临床分类（OA与正常）。在26个候选基因中，NR4A2、BMP1和AVPR1A在所有三个数据集中 consistently 表现出强大的区分能力（AUC > 0.70），强调了它们作为OA诊断可靠生物标志物的潜力。其余基因表现不一，有些仅在特定队列中达到诊断显著性。综上所述，这些结果突出了NR4A2、BMP1和AVPR1A作为在OA背景下值得进一步功能和转化研究的有力候选基因。
功能富集分析
为了研究我们发现的更广泛的功能意义，我们首先使用194个精选通路对GSE114007数据集进行了单样本基因集富集分析（ssGSEA）。然后将ssGSEA富集分数与我们26个候选基因的表达进行相关性分析，并将得到的相关系数可视化在密度图和热图中。值得注意的是，NR4A2、BMP1和AVPR1A——它们在多个队列中表现出一致的诊断能力——也与多个OA相关通路表现出强相关性，加强了它们在疾病发病机制中的潜在作用。
同时，我们使用ClusterProfiler包对GSE114007中的2246个DEGs进行了Gene Ontology（GO）和Kyoto Encyclopedia of Genes and Genomes（KEGG）通路富集分析（p.adjust < 0.05）。结果产生了796个显著富集的生物过程（BP）条目，33个细胞组分（CC）条目，24个分子功能（MF）条目和25个KEGG通路。每个类别中前10个最显著条目的点图说明了OA病理生物学的主要功能主题，包括骨矿化、软骨发育以及关键信号通路，如Wnt、PI3K–Akt和HIF-1。这些分析共同表明，改变的细胞外基质（ECM）稳态和多个关键通路的失调可能促进OA进展。
免疫浸润特征分析
为了探索OA的免疫学背景，我们首先使用GSVA包对GSE114007数据集中的26个免疫细胞类型基因集进行了ssGSEA。然后将得到的免疫浸润分数与我们26个关键OA相关基因的表达水平通过Pearson相关性进行分析。值得注意的是，多个基因——包括NR4A2、BMP1和AVPR1A——与某些T细胞亚型、先天免疫区室和其他免疫细胞群体表现出强相关性，提示它们在OA发病机制中可能具有调节作用。
接下来，我们利用ESTIMATE包计算了GSE114007数据集中每个样本的StromalScore、ImmuneScore、ESTIMATEScore和TumorPurity值。我们再次将这些分数与26个关键基因进行相关性分析，并在第二个热图中可视化结果。与免疫细胞类型数据一致，几个基因与基质和免疫微环境分数显示出显著相关性，强调了它们对OA中软骨退变和局部免疫浸润之间复杂相互作用潜在影响。这些发现共同指出了OA病理学的免疫调节维度，并突出了在关节炎症和组织重塑背景下可能作为治疗或诊断生物标志物的候选基因。
NR4A2、AVPR1A和BMP1与Bexarotene和甲氯芬那酸的分子对接
为了进一步探索我们顶级候选基因——NR4A2、AVPR1A和BMP1——的治疗潜力，我们根据它们与这些基因产物报告的相互作用，从ChEMBL中选择了两种小分子化合物：bexarotene（DB00307）和meclofenamic acid（DB00939）。NR4A2、AVPR1A和BMP1的蛋白质结构从蛋白质数据库（PDB）获得，两种配体的分子文件（MOL格式）从PubChem检索。使用AutoDock（版本1.5.7）进行分子对接。对接复合物表现出不同的结合能和相互作用残基。
值得注意的是，bexarotene对所有三个靶点都显示出有利的结合构象和中等结合能，突出了通过小分子抑制剂或调节剂靶向NR4A2、AVPR1A和BMP1的可行性。甲氯芬那酸也形成了稳定的复合物，尽管相对于bexarotene结合能稍低。这些发现强调了通过小分子抑制剂或调节剂靶向NR4A2、AVPR1A和BMP1的可行性。未来的体外和体内研究需要验证这些对接预测并评估这些化合物在OA中的治疗效果。
Hub基因、富集通路和治疗药物的相互作用网络
为了说明OA中关键通路和推定药物靶点之间的功能和药理学关系，我们构建了一个相互作用网络，包含26个hub基因、前20个最显著富集的GO条目（BP、MF和CC）、前20个富集的KEGG通路以及两种候选小分子化合物（bexarotene和甲氯芬那酸）。如图所示，每个红色倒三角形代表26个与OA发病机制相关的hub基因之一；绿色圆圈表示BP条目，淡蓝色圆圈表示MF条目，黄色圆圈表示CC条目，紫色圆圈对应KEGG通路。通过ChEMBL和计算机对接识别的两种药物以蓝绿色三角形表示。几个hub基因，包括NR4A2、BMP1和AVPR1A，显示出与骨代谢调节剂、糖皮质激素激素和各种功能通路的多重联系，反映了OA中炎症、ECM重塑和软骨细胞稳态之间复杂的相互作用。这个整合网络强调了靶向这些关键基因和通路的潜在治疗策略，提供了驱动OA进展的分子结构的全面视图。
Hub基因在mRNA表达水平上的验证
为了验证关键OA相关基因的表达水平，对OA患者和健康对照的软骨样本进行了定量实时PCR（qRT-PCR）。结果显示，与正常组织相比，OA软骨组织中NR4A2、BMP1和AVPR1A的mRNA表达水平显著上调（p < 0.05）。这些发现表明，这些基因可能通过促进软骨微环境中的炎症、代谢重编程和内皮功能障碍来促进OA进展。qRT-PCR结果与基于GSE114007和GSE220243数据集的生物信息学预测一致，并在临床标本的mRNA水平上得到进一步验证。这些结果支持了我们多组学方法的可靠性，并突出了NR4A2、BMP1和AVPR1A在OA中的诊断和治疗潜力。
讨论
在这项研究中，我们应用scRNA-seq和整合生物信息学来描绘OA的细胞和分子基础。我们的结果揭示了膝关节软骨内的显著异质性，识别了若干DEGs，并通过通路富集、免疫浸润分析和分子对接揭示了潜在的治疗靶点。这些发现不仅证实了软骨退变的既定机制，而且为疾病发病机制和潜在干预策略提供了新的见解。
我们的scRNA-seq分析揭示了膝关节软骨中八个不同的软骨细胞亚群，正常和OA组织在它们的转录谱和相对丰度上存在显著差异。这些发现与先前的报告一致，表明软骨细胞不是同质群体，而是包含在软骨稳态和退变中具有不同功能的功能离散亚型。
值得注意的是，我们观察到在大多数软骨细胞簇中，涉及炎症和组织重塑的通路，如TNF、TGF-β和PI3K–Akt信号通路，显著富集。这些通路长期以来被认为与OA进展有关，其中促炎细胞因子如TNF-α驱动分解代谢过程，而失调的TGF-β和PI3K-Akt信号通路导致异常的软骨细胞增殖和ECM转换。我们的数据通过揭示特定亚群的软骨细胞可能优先参与这些信号级联，增加了粒度性，强调了在OA治疗中亚群水平靶向的重要性。
通过重新分析GSE114007批量转录组数据集，我们识别了2247个DEGs，其中许多反映了OA的标志性过程——即ECM重塑、软骨稳态和炎症信号。值得注意的是，几个胶原和基质相关基因（例如COL1A1和POSTN）是OA软骨中上调最显著的基因之一，这与表征该疾病的既定ECM失调一致。相反，如HILPDA、ADM和ATF3等基因显著下调，提示细胞应激反应和代谢调节可能存在缺陷。
通过将这些DEGs与ChEMBL中的205个推定OA药物靶点取交集，我们将焦点缩小到26个关键基因，包括NR4A2、BMP1和AVPR1A。ROC曲线分析证实了这三个基因在多个队列中具有强大的诊断效用，突出了它们作为临床相关生物标志物的前景。尽管BMP1先前已与软骨和骨稳态相关联，但我们的发现通过证明其在OA中的持续上调和与炎症信号的潜在关联扩展了其重要性。类似地，NR4A2（一种核受体）和AVPR1A（一种加压素受体）在炎症和组织重塑中都有文献记载的作用，强调了它们作为OA驱动因素的合理性。
对已识别DEGs的GSEA进一步支持了OA涉及ECM调节、细胞增殖以及关键信号通路（如Wnt、PI3K–Akt和HIF-1）的多因素扰动。这些通路共同控制软骨细胞存活、代谢和基质合成，它们的失调与进行性软骨破坏一致。同时，我们的免疫浸润分析突出了软骨细胞与各种免疫细胞类型（包括T细胞亚群和先天免疫细胞）之间潜在的相互作用。观察到的NR4A2、BMP1和AVPR1A表达与免疫细胞浸润分数之间的相关性尤其引人注目。这一发现支持了OA不仅是一种机械磨损疾病，而且也是一种由低度炎症驱动的疾病的新兴观点。先前的研究已经记录了OA滑膜和软骨中的免疫细胞浸润，但我们的结果表明某些软骨细胞亚群可能主动调节或响应这些免疫元素。
为了解决我们关键候选基因的转化潜力，我们进行了与两种小分子化合物——bexarotene和甲氯芬那酸——的分子对接，预测它们靶向NR4A2、AVPR1A和BMP1。两种配体都与这些蛋白质形成了稳定的复合物，尽管结合能和相互作用残基不同。尽管这些发现是初步的，但它们为未来的实验验证奠定了基础。Bexarotene显示出有利的结合构象，表明它可以调节OA病理生理学中的关键分子相互作用。
尽管有这些有希望的见解，必须承认某些局限性。首先，我们的scRNA-seq数据集来自相对较少的患者，可能无法完全捕捉OA的个体间变异性。需要更大、更多样化的队列来验证我们发现的普遍性。其次，我们对接实验的计算机性质需要严格的生化和临床前测试来确认实际的结合亲和力和治疗效果。第三，尽管进行了免疫浸润分析，但需要进一步的工作来描绘关节微环境中软骨细胞基因表达变化与免疫细胞招募或激活之间的因果关系。
最后，尽管我们的分析 pinpointed 关键通路和潜在治疗靶点，但OA的多因素病因意味着单靶点方法可能不足以阻止疾病进展。同时解决机械应力、炎症和ECM稳态的联合疗法可能对稳健的临床获益是必要的。
结论
总之，这项整合研究以单细胞分辨率提供了OA软骨的全面描绘，识别了关键分子驱动因素，并为治疗开发提出了新的方向。通过结合scRNA-seq、差异表达分析、免疫谱分析和分子对接，我们强调了OA发病机制中软骨细胞亚群、炎症和ECM失调之间复杂的相互作用。我们的发现突出了NR4A2、BMP1和AVPR1A作为高优先级生物标志物和药物靶点，为未来开发能够改善或预防OA软骨退变的精准干预措施的努力提供了宝贵的基础。
