论文解读 | 用公共数据库做单细胞RNA数据挖掘：单一表型下不同

单细胞测序揭示了骨关节炎发展的潜在发病机制

一、文章基本信息

（1）题目：单细胞测序揭示了骨关节炎发展的潜在发病机制（2020）
（2）期刊名：Gene
（3）2019年影响因子：2.984
（4）作者单位：温州医科大学附属医院风湿免疫科
（5）摘要：骨关节炎（OA）是一种慢性退行性改变，发病率很高，导致生活质量下降和社会经济负担增加。 这项研究旨在探讨与OA相关的潜在关键基因和通路，这些基因和途径可用作早期治疗的潜在生物标记。 从公共数据库（GSE104782和GSE109449）下载了OA中1464个软骨细胞和192个成纤维细胞的单细胞基因表达谱，用于后续分析。 使用Seurat和SingleR软件对OA软骨细胞和OA成纤维细胞中的细胞亚群进行了鉴定，共发现8种软骨细胞亚群和3种成纤维细胞亚群。此外，鉴定了44个成纤维样软骨细胞和成纤维细胞之间的共同标记基因，并且通过功能富集分析进一步将粘着斑通路确定为OA的重要潜在机制。 此外，组织和细胞水平的逆转录实时定量PCR（RT-qPCR）实验证实，两个关键标记基因（COL6A3和ACTG1）可能参与了OA的发展。 综上所述，我们推断OA中的软骨细胞可能通过上皮黏附途径上调COL6A3和ACTG1的表达以完成成纤维细胞转化。这些发现有望进一步了解OA纤维化过程的发展，并为早期OA的治疗提供有希望的目标。
（6）研究目的：探索OA软骨细胞纤维化的过程、关键基因和通路。
（7）主要发现：

1. 鉴定了OA中8种软骨细胞亚群和3种成纤维细胞亚群
2. 鉴定了44个成纤维样软骨细胞和成纤维细胞之间的共同标记基因
3. Marker基因的功能富集分析指出了粘着斑通路是OA的重要潜在机制
4. 逆转录实时定量PCR（RT-qPCR）验证了关键标记基因（COL6A3和ACTG1）可能参与了OA的发展
5. 总结，OA中的软骨细胞可能上调COL6A3和ACTG1的表达，通过粘着斑途径完成成纤维细胞的转化

二、单细胞数据库来源：GEO（GSE104782和GSE109449）

（1）GSE104782

1. 源数据发表的文章：单细胞RNA-seq分析揭示了人类骨关节炎（OA）的进展
2. 源数据发表期刊：Ann Rheum Dis 2018 （IF2019: 16.192)
3. 源数据情况：10例行膝关节置换术的OA患者的1464个软骨细胞
4. 源文章的主要发现：
   1）在OA软骨中鉴定出7个软骨细胞亚群，包括三个功能不同的新表型。
   2）发现在增殖性软骨细胞、增生前软骨细胞和肥大软骨细胞（HTCs）之间存在一种潜在的转换，并在HTCs中定义了一个新的亚群。
   3）我们揭示了软骨祖细胞（CPCs）的新标记物，并通过计算分析证明了CPCs和纤维软骨细胞之间的关系。
   4）得出了与临床结果相关的预测指标，阐明了不同细胞类型在OA早期诊断和治疗中的作用。

（2）GSE109449

1. 源数据发表的文章：类风湿关节炎（RA）中功能不同的疾病相关成纤维细胞亚群
2. 源数据发表期刊：Nat Commun 2018 （IF2019: 12.121)
3. 源数据情况：上样细胞使用蛋白质表面标记物对滑膜成纤维细胞进行门控，数据包括来自2名骨关节炎患者和2名类风湿关节炎患者的384个细胞。Li et al.的研究采用了其中2名OA患者的192个滑膜成纤维细胞细胞。
4. 源文章主要发现：
   1）使用靶向亚群的bulk转录组学和单细胞转录组学，揭示人类滑膜组织成纤维细胞亚群之间的功能和转录差异。
   2）确定七个具有不同表面蛋白表型的成纤维细胞亚群，通过整合转录组数据将其归纳为三个亚群。
   3）相对于OA患者，RA患者的成纤维细胞亚群的特征是蛋白质Podoplanin，THY1膜糖蛋白和cadherin-11的高表达，但缺乏CD34。相对于OA，RA患者滑膜的成纤维细胞亚群数目也成倍增长。
   4）单细胞转录组和免疫组化共同显示，成纤维细胞位于发炎的滑膜中的血管周围区域，分泌促炎细胞因子，具有增生作用，并具有浸润细胞的体外表型特征。

三、分析思路和结果

1 研究基础

1）软骨细胞是健康软骨的主要成分，产生并维持软骨基质，主要由胶原蛋白和蛋白聚糖组成。 成纤维细胞是结缔组织最常见的细胞，可合成胶原蛋白和细胞外基质。 在OA开始时，软骨细胞的合成代谢能力大大减弱，从而损害软骨修复。 在晚期阶段，软骨细胞分化成软骨细胞和成纤维细胞之间的纤维软骨细胞，产生异常成分，例如纤连蛋白碎片。 此外，根据体外研究的最新发现（Deroyer等人，2019），软骨细胞被报告为调节过程增殖和转分化为“软骨肌成纤维细胞”。 随着OA的发展，处于同一病态组织中的处于不同病理阶段的细胞将会出现，在疾病的发展中起着独特的作用。 因此，通过scRNA-seq分析鉴定不同病理阶段细胞簇的研究工作将有助于更好地了解OA的病因和进展。

图1 设计的分析和实验流程

2 方法与结果

（1）数据获取和质控：

1. 方法步骤
   1)对OA软骨细胞和成纤维细胞分别进行质控，分别取每百万转录本（TPM）值和log2（TPM值+ 1）值用于后续分析；
   2)去除表达基因在200—7500之外的细胞，表达量在10以下的细胞和线粒体基因表达超过5%的也被去除。
   3分别对两个数据集进行PCA降维。
2. 结果
   OA软骨细胞剩余1464个细胞，192个成纤维细胞.

(2) 分群和鉴定各软骨细胞亚群的标记基因：

图2 分群和鉴定各软骨细胞亚群的标记基因

1. 方法步骤
   1）使用聚类和T-SNE分析鉴定出细胞亚群；
   2）鉴定了每个亚群的标记基因，用每个亚群的前10个Marker基因绘制了差异基因表达热图；
   3）对每个亚群进行了top2 Marker基因的t-SNE结果映射可视化；
   4）进行了拟时序分析；
2. 结果
   1）通过t-SNE分析，在1,464个细胞中发现了8个不同的簇，通过Marker基因标记发现其中7种为软骨细胞，1种为纤维样软骨细胞（图2A）；
   2）热图显示8种亚群可以被Marker基因很好的区分开（图2B）；
   3）每个亚群的top2 Marker基因和亚群的分布基本一致(APOD, SDHA, AKR1C2, PRG4, CCL2, S100A2, FOS, HES1, COL10A1, IBSP, KRT17, KRT16, COL1A1, IFI27, SLC5A12, and S100A9)（图2C）；
   4）拟时序分析表明软骨细胞功能分化出现了4个分支，成纤维样软骨细胞仅分布在分支3后面的轨迹末端，这暗示了从软骨细胞向成纤维细胞去分化的潜在过程。 （图2D）

（3）纤维细胞样软骨细胞亚群中标记基因的功能富集分析和蛋白质相互作用分析

图3 基于GO的差异基因富集分析气泡图

1. 方法步骤（GO）
   1）（基因本体GO分析）为了进一步说明与成纤维样软骨细胞亚群的标记基因相关的生物学过程（BP），使用 g：Profiler在线工具进行功能富集分析并可视化不同簇之间的相互作用。
2. 结果
   1）标记基因在生物学过程（BP）中显着富集，包括细胞外基质组织（GO：0010033），细胞外结构组织（GO：0001501），对有机物的响应（GO：0071310），骨骼系统发育 （GO：0070887），细胞对有机物的反应（GO：0009888）等，它们在细胞外基质中起着至关重要的作用。 ECM是由胶原蛋白，酶和糖蛋白组成的细胞外大分子，ECM的合成代谢/分解代谢平衡破坏被认为是OA发病机理的主要事件。
   2）对于细胞成分（CC），发现这些标记基因主要富含细胞外基质（GO：0031012），含胶原的细胞外基质（GO：0062023），细胞外区域部分（GO：0044421），细胞外空间（GO：0005615）和细胞外区域（GO：0005576）。此外，标记基因显着丰富了细胞外基质结构成分（GO：0005201），结构分子活性（GO：0005198），赋予抗张强度的细胞外基质结构成分（GO：0030020），生长因子结合（ GO：0019838）和在分子功能组中赋予抗压性的细胞外基质结构成分（GO：0030021） 。细胞外结构组织和细胞外基质组织也与软骨细胞的生长以及胶原蛋白的重塑密切相关（Bella and Hulmes，2017）。
   3）对于分子功能（MF），Marker基因显着丰富了细胞外基质结构成分（GO：0005201），结构分子活性（GO：0005198），赋予抗张强度的细胞外基质结构成分（GO：0030020），生长因子结合（GO：0019838）和具有抗压性的细胞外基质结构成分（GO：0030021）。GO分析表明细胞成分主要位于细胞外空间，细胞外区域部分和细胞外区域，与ECM的位置一致。

图4 KEGG

1. 方法步骤（KEGG）
   1）（KEGG分析）

2. 结果
   1）KEGG通路分析表明，在OA的成纤维样软骨细胞亚群中，7种通路表现出明显的富集，包括局灶黏附，ECM受体相互作用，蛋白质消化和吸收，TGF-β信号通路， 癌症中的核糖体和蛋白多糖通路（图4）。功能富集分析（包括BP和KEGG通路分析）表明，标记基因与粘着斑通路，ECM-受体相互作用和TGF-β信号通路密切相关，表明这些标记基因可能参与了 OA中软骨细胞转化为成纤维细胞和胶原纤维化的必要性。
   2）据报道，粘着斑复合物是细胞ECM粘附的先决条件，通过细胞粘附和迁移，粘着斑激酶（FAK）的活化可能与OA有关，而关于粘着斑途径与粘连的关系尚缺乏具体的研究。 OA（Shahrara等，2007； Prasadam等，2013）。

   
   图5 PPI网络
   图6 八种软骨细胞亚群中7个候选基因的表达情况

1.方法步骤（PPI）
1）纤维样软骨细胞亚群中，粘着斑途径的text-mining（STRING）和PPI网络分析（图5）
差异基因在不同细胞亚群中的表达情况（图6）

2. 结果
   1）通过text-mining（STRING），确定了粘着斑途径是重要的途径。
   2）粘着斑途径的PPI网络表明在成纤维样软骨细胞簇中使用相关的标记基因紧密相互作用，其中COL1A1，COL1A2，COMP，CHAD，COL6A1和COL6A3表现出高度相关的核心位置。
   3）如图6，值得注意的是，在成纤维样软骨细胞簇中，一些未被报道的标记基因（RHOA、FN1、COL6A1、CHAD、CAV1、ACTG1和COL6A3）被鉴定为候选基因，并且它们也富集在粘着斑途径中高表达，特别是RHOA、COL6A1和COL6A3。此外，FN1、CHAD和CAV1在每个簇之间的表达水平基本上是无差别的，表明簇特异性较低。然而，其他4个基因的表达水平存在显著差异，COL6A3在成纤维细胞样软骨细胞群中高表达。

（4）候选基因的基因集富集分析（GSEA）

图7 Marker基因GSEA和RT-qPCR结果. （A）双峰分组候选基因得分阈值的双正态分布模型。红线代表低分组候选基因得分分布，绿线代表高分组候选基因得分分布。（B）GSEA显示了与炎症相关的基因集，这些基因集富含具有标记基因的样品，这些标记基因包括IgA产生的免疫网络，自然杀伤细胞介导的细胞毒性，细胞粘附分子（CAM），细胞因子-细胞因子受体相互作用，T细胞受体信号传导途径和造血细胞谱系；

1. 方法步骤
   1）GSEA基因集富集分析：为了进一步研究候选基因在OA机制中的潜在作用，根据候选基因水平将软骨细胞分为高亚组和低亚组。 我们将候选基因的表达转化为候选基因分数，并将2.6e-04作为根据双正态分布模型分组的阈值(图7a)。
2. 结果
   1）GSEA结果表明，候选基因的低表达并没有在任何途径有效富集（p>0.05）。而高表达的候选基因与炎症途径有关，例如产生IgA的肠道免疫网络（标准化富集分数（NES）=2.091，FDR=9.96E−04）、自然杀伤细胞介导的细胞毒性（NES=2.056，FDR=1.36E−03）、细胞粘附分子（CAMs）（NES=2.115，FDR=2.05E−03），细胞因子-细胞因子受体相互作用（NES=2.002，FDR=2.85E−03）、T细胞受体的信号转导(NES = 1.974, FDR = 4.06E−03)和造血细胞谱系(NES=2.133，FDR=4.10E−03)。（图7b）

（5）OA中每种成纤维细胞亚群的标记基因（第二个数据集）

图8a-d 共同Marker基因鉴定
图8e KEGG分析

1. 方法步骤
   1）t-SNE分析
   2）拟时序分析
   3）维恩图展示1、2、3号成纤维细胞亚群分别与纤维样软骨细胞亚群共同Marker基因鉴定
   4）共同Marker基因的KEGG通路分析
   鉴定纤维样软骨细胞亚群中未被报道的候选基因与共同Marker基因的共同关键基因
2. 结果
   1）192个成纤维细胞分为3个不同亚群
   2）拟时序分析表明2号成纤维细胞亚群与其他两种亚群在发育上存在较大差别。
   3）共同Marker基因鉴定表明，2号成纤维细胞亚群和纤维样软骨细胞亚群之间成功鉴定出44种常见标记基因，而在其他亚群中几乎没有共同的标记基因。
   4）KEGG通路分析了2号成纤维细胞亚群和纤维样软骨细胞亚群共同的44种常见标记基因，明显富集到了胶原沉积相关的途径，尤其是在粘着斑粘附和ECM受体相互作用途径。
   5）44个常见标记基因与7个候选基因之间重合的结果确认了四个未被报道的关键基因（COL6A1，COL6A3，ACTG1和RHOA）。

（6）通过RT-qPCR验证OA中的关键基因

图9 与OA组和对照组相比，COL6A3（a）（e）和ACTG1（c）（f）的相对表达差异显着，而滑膜组织和软骨细胞中COL6A1（b）和RHOA（d）则无显着差异

1. 方法步骤
   1）使用RT-qPCR对比OA患者和对照组滑膜组织中的COL6A1，COL6A3，ACTG1和RHOA表达，选出显著差别的基因A等
   使用RT-qPCR对比OA患者和对照组软骨组织中A等基因的表达，选出显著差别的基因B等。
2. 结果
   1）OA滑膜组织中ACTG1和COL6A3的表达水平均显着高于对照组滑膜组织，但RHOA、COL6A1表达差异不显著。
   OA软骨组织中ACTG1和COL6A3的表达水平均显着高于对照组软骨组织，表明ACTG1和COL6A3参与了软骨组织的纤维化过程。

写在最后的

此文章是一篇典型的单一表型下不同组织/细胞大类间的单细胞测序分析论文，其分析思路（故事思路）基本如下：
（1）首先，在公开数据库中检索关键词寻找拟研究细胞大类（本文是OA软骨细胞和OA成纤维细胞）的单一表型下同一组织/细胞大类的单细胞表达谱类文章。从而获得两种同一表型下两种或多种组织/细胞大类的基因-barcodes表达矩阵。
（2）其次，两种数据集分别进行质控、降维、聚类、二维嵌入可视化、注释和差异基因分析、拟时序分析的基础分析，目的是根据先验知识标注细胞亚群，找到两种组织/细胞大类中有明显分化联系的细胞亚群，即有很多共同DEGs的，且分别在两个拟时序分析图像上独立（这说明要么处于分化末要么位于分化初，可能有分化联系）的两个亚群（话说能否根据Marker基因把不同数据集的几个亚群投射到一个拟时序分析图谱上，应该可以）。
（3）接下来，对发育起点的细胞亚群（本文中是成纤维样软骨细胞亚群）进行差异表达基因的功能富集分析（GO和KEGG），找到富集的生物学过程、细胞成分和分子功能（GO），以及富集的通路（KEGG），并根据先验知识进行解释，解释的目的是增强研究目的与找到的拟研究Marker基因之间关联的说服力，若之前已有人研究过该课题，那么可以尝试寻找没有人解释过的通路。在本文中，作者们重点研究了黏着斑通路。
（4）在确定了黏着斑通路之后，就可以获得该通路相关的DEGs，利用蛋白质—蛋白质相互作用网络（PPI）寻找DEGs相互作用情况和核心基因。粘着斑途径的PPI网络表明，与该通路相关的成纤维样软骨细胞簇Marker基因呈现紧密相互作用，其中COL1A1，COL1A2，COMP，CHAD，COL6A1和COL6A3表现出高度相关的核心位置。（这步其实对讲故事关系不大，就是做个图，关键在下面）
（5）在确定了黏着斑通路之后，就可以获得该通路相关的DEGs，寻找其中未被报道过与软骨组织纤维化有关的基因作为下一步分析的候选基因，在本文中RHOA、FN1、COL6A1、CHAD、CAV1、ACTG1和COL6A3被鉴定为候选基因。
（6）分析候选基因在软骨组织各个细胞亚群中的表达情况，将表达差异明显的几个基因作为重点分析对象（表达上调或下调都是重点，重点去查他们的文献，因为用他们讲故事会比较漂亮），当然其他几个基因也得查。
（7）接下来做GSEA，分别研究高低表达的候选基因在各个通路上的富集情况，这步也是增强结果的可解释性。
（8）接下来联合分析两个数据集/组织，（3）—（7）研究了分化起点的亚群的Marker基因，确定了候选基因。下一步进行另一个数据集中承接成纤维样细胞亚群的关联细胞亚群。在本文中，是根据拟时序分析和差异表达基因相似性确定的，即2号成纤维细胞亚群。
（9）找到44个共同的Marker基因后，在做一遍KEGG，又重新确认了这些基因富集到的通路和功能的确和纤维化有关。
（10）（3）—（7）中我们找到了7个候选基因，在这里研究者又找了44个共同基因和7个候选基因之间的重合基因，一共4个，作为最后的候选基因做RT-qPCR验证。（不过我认为就算不重合也能解释的通，因为成纤维细胞的差异基因不是与软骨细胞亚群比较得来的，3个不重合的基因只是在2号成纤维细胞与其他两种成纤维细胞比较时不显著差异，不代表2号成纤维细胞中这3个基因与其他软骨组织细胞亚群没有差异）
（11）采集OA和健康人的临床样本，用RT-qPCR验证4个候选基因是否在OA和健康人中差异表达。结论是只有2个基因差异表达，那么这两个基因就参与了软骨细胞的成纤维化。
（不过，这里还有个问题，4个候选基因是12例OA患者的测序数据经过了各种统计学分析得出的科学结论，这个结论的偏差似乎可以用机器学习中的训练集和真实世界测试集的分布偏差来解释。然而，RT-qPCR中也只有40例（其中20例OA），这里得到的2个基因无效，是否需要进行进一步的证据合并来得出最终结论呢？比如meta-analysis？