基于单细胞转录组测序分析骨关节炎软骨细胞分化的分子机制*

徐文飞 梅其杰 明春玉 段戡 袁长深 郭锦荣 胡琪 曾超 梁晓辉

骨关节炎是一种常见的慢性退行性疾病，常伴有关节肿胀和软骨损伤等变化[1]；因其具有高发病率和致残率，将会成为全球主要致残疾病之一[2-4]。然而目前骨关节炎的发病机制仍存在争议，阐明发病机制成为迫切需要解决的问题。而软骨细胞对于骨关节炎具有至关重要的作用，其中软骨细胞肥大或衰老担任着重要角色[5]，故软骨老化或过度使用将导致软骨细胞稳态丧失，进而加重骨关节炎病情[6]。因此，深入研究骨关节炎与软骨细胞之间关系，从分子水平阐述骨关节炎软骨细胞发病机制成为新热点和机遇。

单细胞转录组测序（single-cell RNA sequencing,scRNA-seq）是通过对组织中细胞进行独立分析，来揭示细胞亚群之间相互作用的关系，以及细胞内基因之间相互调控的过程，较传统的RNA 测序技术分析更加详细，对疾病研究具有重要作用[7-9]。当前单细胞转录组测序已广泛运用于多种疾病（包括骨关节炎），尚缺乏单细胞测序联合生物信息学（基因芯片）关于骨关节炎软骨细胞的分析，因此本研究在其他研究的基础上将两者相结合，共同探讨骨关节炎软骨细胞的作用机制，以期更好地治疗骨关节炎。

1 材料与方法

1.1 基因芯片数据来源

通过GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）下载骨关节炎软骨细胞基因芯片数据集GSE104782，其所采用平台为GPL20301，该芯片数据包括10 个骨关节炎样本，1 464个软骨单细胞测序信息。

1.2 质控与数据过滤

运用R 语言对GSE104782 进行质控与数据过滤，将重复基因取均值来进行数据质控，并转换为“Seurat”对象，进一步通过“Percentage Feature Set”函数计算线粒体基因百分比，进而将数据过滤，从而获得基因特征及测序深度，并对数据标准化，提取细胞间变异系数较大的基因。

1.3 主成分分析

主成分分析（principal component analysis, PCA）是一种最常用的线性降维方法，能确定可用维度并筛选相关基因。运用R语言对数据进行标准化预处理，进而使用PCA分析，可获得每个PCA成分相关基因、P值及均匀分布情况。

1.4 TSNE聚类分析与Marker基因

TSNE 聚类分析是可视化高维数据工具，可在二维或三维空间里展示聚类结果。采用R语言对单细胞测序所得数据分组，优化标准模块，进行可视化聚类分析；同时，通过寻找数据差异表达特征获得Marker基因，并对各个成分中所含Marker基因进行分析。

1.5 细胞注型及分化轨迹分析

基因芯片数据集（GSE104782）所包含细胞类型均为软骨细胞，采用R 语言中的SingleR、Monocle 包进行细胞注释，进一步验证该数据集细胞类型，并运用伪时间分析法来阐述软骨细胞分化轨迹，同时查阅相关文献深入分析骨关节炎软骨细胞的具体分化途径。

1.6 GO及KEGG分析

利用R 语言工具对软骨细胞的Marker 基因进行GO 富集及KEGG信号通路分析，获得软骨细胞相关功能及信号通路，其中显著性基因富集临界值设定为P＜0.05。

1.7 蛋白互助网络分析

选取细胞成分9 中所含Marker 基因进行蛋白互助网络分析，将Marker 基因导入STRING（Version 11.0）在线分析软件中（https://string-db.org/cgi/input.pl），获得蛋白质-蛋白质相互作用（protein-protein interaction,PPI），并将其结果导入Cytoscape 3.7.2 软件，运用插件cytoHubba 中Degree 算法获得前10 位关键靶基因并构建PPI 网络。

2 结果

2.1 质控与数据过滤结果

运用R 语言对基因芯片数据集GSE104782 进行分析，发现基因测序深度与基因特征具有显著联系，随着测序不断深入，基因特征也更加丰富（见图1）；该数据集获得基因17 380 种，筛选出COL1A1、IBSP、PRG4、TAC1、SPP1、HSPA6、COL10A1、CCL2等1 500 种重要基因（见图2）；同时分析每种样本基因特征及研究深度，发现基因特征主要2 500～5 000范围内，而研究深度在0 ～ 500 000之间（见图3）。

图1 测序基本特征与深度关联

图2 测序基因表达

图3 测序的基本特征与深度

2.2 PCA分析结果

通过PCA分析相关数据发现，软骨细胞由20个主成分组成，按P＜0.05 进行统计学分析，结果显示20 个主成分具有统计学差异，均可纳入研究（见图4）。以PC1 与PC2绘制主成分分析图形，发现软骨细胞主成分之间具有一定的关联（见图5）。同时，进一步分析PC1、PC2、PC3 及PC4 的相关基因，绘制主成分热图，展示样本中主成分所含基因分布及聚类情况（见图6、图7），结果显示PC1 中COL1A2基因、PC2中OGN基因、PC3中HLA-DRA基因及PC4中SOD2基因在样品组织中显著高表达。

图4 主成分分析结果

图5 PC1与PC2的关联

图6 前4位主成分基因分布曲线图

图7 前4位主成分基因分布热图

2.3 TSNE聚类分析结果

本研究进一步对单细胞测序所得的数据进行TSNE聚类分析，共获得9 类不同颜色的细胞聚类，分别代表软骨细胞中不同成分（见图8）；同时以校正P=0.05 为标准进行差异分析，获得1 886 个Marker 基因并以热图及气泡图的方式显示有些基因在不同成分中的表达情况（见图9、图10），同时围绕关键基因LECT1、TGFBI分析表达及分布情况，发现LECT1（logFC=2.33，P=0.045）、TGFBI（logFC=1.66，P=0.044）均在细胞成分9 种高表达（见图11）。

图8 软骨细胞成分分布

图10 Marker基因分布气泡图

图11 关键基因的表达与分布

2.4 细胞注释及细胞轨迹分析结果

运用R语言工具中SingleR包对获得的9种细胞成分进行注释，发现均为软骨细胞，与研究内容相符；通过伪时间分析法进一步阐述软骨细胞的分化轨迹（见图12），发现共有6种分化途径，其中显示Cluster0细胞聚类单一，为最早分化细胞，查阅相关文献发现Cluster0 可能为软骨祖细胞，经关键调控因子分化为增殖软骨细胞、肥大前软骨细胞和肥大软骨细胞等不同软骨细胞成分[10]。

图12 软骨细胞的分化轨迹

2.5 GO和KEGG通路富集分析

GO 是一组大量的基因注释术语，可用于识别高通量基因组或转录组数据的特征生物学属性。采用R 语言工具对Marker 基因进行GO 功能富集分析，以校正P＜0.05 进行筛选，主要涉及软骨发育、结缔组织发育等生物过程（biological process, BP）；涉及含胶原蛋白的细胞外基质、内质网内腔等细胞组分（cellular component, CC）；涉及肝素结合、糖胺聚糖结合等分子功能（molecular function,MF）（见图13）。KEGG 通路分析：主要与蛋白质消化吸收、ECM受体相互作用、人乳头瘤病毒感染、粘着斑、癌症中的蛋白多糖、补体与凝血级联反应、PI3K-Akt信号通路、矿物吸收等信号通路有关（见图14）。

图13 GO富集分析

图14 KEGG信号通路分析

2.6 差异基因的PPI网络分析

STRING 数据库是通过已知和预测的蛋白质相互作用数据所组成，本研究选取细胞成分9 种所含Marker 基因进行蛋白互助网络分析，运用STRING 数据库获得PPI 网络图，并且将其所得结果导入Cytoscape 3.7.2 软件，运用插件cytoHubba 中Degree 算法发现前5 位关键靶基因分别为COL1A1、COL2A1、VEGFA、COL1A2、DCN（见图15）。

图15 蛋白质互助网络分析

3 讨论

骨关节炎是一种由多种因素引起的慢性退行性疾病[11]，主要病理表现为软骨细胞凋亡和细胞外基质丢失，故维持软骨正常内环境稳定是治疗骨关节炎的关键所在[12]。软骨是一种无血管、淋巴和非神经支配的组织，而软骨细胞作为软骨唯一的细胞类型，具有内在修复潜力[13]；然而如何保护及修复软骨细胞仍然十分艰难，故本研究采用单细胞转录组进一步探究软骨细胞基因表达及分化途径，获得软骨保护及修复方法，为更好治疗骨关节炎提供新思路。

本研究通过单细胞测序及生物信息学分析发现17 380种软骨细胞基因，其中包含LECT1、TGFBI等1 886种关键基因。其中LECT1在骨关节炎软骨细胞中呈高表达状态，其作为软骨蛋白-1编码可抑制软骨细胞生长，进而影响骨关节炎病情[14]；而TGFBI不仅调节合成代谢和分解代谢标志物的表达，而且还降低炎症因子[15]。同时进一步研究发现软骨细胞共包含多种软骨细胞类型及分化途径，其中软骨祖细胞为最早分化细胞，可经关键调控因子分化为增殖软骨细胞、肥大前软骨细胞和肥大软骨细胞等不同成分，因此本研究推测干预软骨祖细胞可能为修复及保护软骨细胞提供新的方法。

采用GO、KEGG及蛋白互助网络分析软骨细胞具体机制，GO 富集分析发现软骨细胞基因涉及软骨发育、结缔组织发育等生物过程（BP）；涉及含胶原蛋白的细胞外基质、内质网内腔等细胞组分（CC）；涉及肝素结合、糖胺聚糖结合等分子功能（MF）。其中，软骨细胞及结缔组织发育可起着修复软骨作用，且关节软骨退变与细胞外基质降解密切相关[16]；分泌型肝素结合肽通过刺激细胞外基质合成、减少降解基质金属蛋白酶和诱导其抑制剂来刺激软骨细胞聚集和增殖，增强软骨形成[17]；糖胺聚糖链连接的蛋白聚糖存在于软骨细胞表面和细胞周围软骨基质中，可参与细胞-细胞和细胞-基质的相互作用[18-19]。因此，本研究认为从生物过程、细胞组分、分子功能角度保护软骨细胞可能会具有显著的疗效。

KEGG通路分析主要与蛋白质消化吸收、ECM受体相互作用、人乳头瘤病毒感染、粘着斑、癌症中的蛋白多糖、补体与凝血级联反应、PI3K-Akt信号通路、矿物吸收等信号通路有关；其中发现黄芪甲苷可通过上调ECM受体相互作用来促进软骨细胞的快速增殖，对关节软骨退变发挥了保护作用[20]；而补体与凝血级联及金黄色葡萄球菌感染主要参与天然免疫反应，可能提高软骨细胞的修复作用[21-22]；同时，在骨关节炎软骨细胞中PI3K-AKT 通路处于下调状态，通过激活PI3K/AKT 可减少软骨细胞损伤，进而促进骨关节炎软骨细胞再生[23]。

通过Cytoscape 3.7.2 软件获得5 个关键靶基因，包含COL1A1、COL2A1、VEGFA、COL1A2、DCN。在软骨受到损伤后，可促进胶原蛋白（COL1A1、COL2A1、COL1A2）下调，使得胞外囊泡-壳寡糖偶联物逆转软骨损伤而起到保护作用[24]；血管内皮生长因子（vascular endo‐thelial growth factor, VEGF）是血小板衍生生长因子家族的重要成员之一，主要分布于关节的成骨细胞、软骨细胞及滑膜内成纤维细胞等组织[25]，其主要通过促进血管内皮增殖、促进肿瘤淋巴管生长、上调抗凋亡蛋白Bcl-1表达促进骨关节炎的形成[26]；而核心蛋白聚糖（decorin, DCN）可延缓软骨蛋白聚糖碎片的丢失和软骨表面的纤维化，从而对改善软骨退变起到保护作用[27]。

综上所述，本研究运用单细胞测序联合生物信息学方法研究骨关节炎软骨细胞分化的具体机制，结果显示软骨细胞分为多种类型且具有多种分化途径，可通过调节关键基因（LECT1、TGFBI）来干预软骨细胞再分化，进而为保护软骨细胞起着重要作用。然而本研究仍存在不足之处，首先未对结果采用试验进行验证；同时纳入单细胞基因芯片样本较少，可能存在一定偏倚，需后续扩大样本进一步研究，因此期待后续能更深入研究软骨细胞相关机制，为治疗骨关节炎提供坚实基础。