摘要

对于所有类型的高通量单细胞数据，细胞类型分配（cell type assignment）是一个主要挑战。在许多情况下，这种分配需要反复手动使用外部和补充数据源。为了提高跨合作、平台的统一分配细胞类型的能力，我们开发了Cellar，这是一种软件工具，为分配和数据集比较过程中涉及的所有不同步骤提供交互式支持。我们将讨论Cellar实现的不同方法，如何将这些方法用于不同的数据类型，如何组合互补的数据类型，以及如何分析和可视化数据。我们通过使用Cellar注释来自多组学单细胞测序和空间蛋白质组学研究的几个HuBMAP数据集，展示了Cellar的优势。Cellar是开源的，包括几个带注释的HuBMAP数据集。

引言

包括人类生物分子图谱计划（HuBMAP，Human BioMolecular Atlas Program）在内的许多大型团队都专注于在单细胞水平上分析组织、器官和整个人体。这些团队使用几种不同的技术来研究单细胞的分子组成，包括单细胞RNA测序、单细胞ATAC测序、单细胞空间转录组学和单细胞空间蛋白质组学。除了这些大型团队，个别实验室还使用部分或全部这些模式生成数据。

在过去几年中，已经开发了许多方法来分配单细胞数据中的细胞类型。在大多数情况下，来自同一团队的不同小组，甚至在处理多种类型的单细胞数据时，同一小组，都依赖于不同的工具集。这使得整合和比较这些小组的数据变得困难，因为研究人员经常使用不同的分配技术、markers，甚至细胞类型的命名约定。（目前对单细胞数据的管理缺乏一个统一）

为了实现跨越不同单细胞组学平台和模式的大规模合作、集成和比较，我们开发了Cellar，这是一个交互式和图形化的细胞类型分配的web服务。Cellar实现了一套全面的方法，包括现有的方法和新的方法，涵盖了细胞类型分配过程中涉及的所有步骤。这些方法包括降维（dimensionality reduction）和表示（representation）、聚类（clustering）、基于参考的比对（reference-based alignment）、差异表达基因的识别（identification of differentially expressed genes）、功能集和标记集的交叉（intersection with functional and marker sets）、管理会话和输出结果的工具（tools for managing sessions and exporting results），以及同时分析和比较两个数据集的双重模式（dual mode for analyzing and comparing two datasets simultaneously）。由于细胞类型分配通常需要用户以领域知识的形式输入，Cellar采用半自动解决方案，允许用户根据自身需要去干预和修改每个处理步骤。为了实现这种交互式分析，Cellar提供了在空间单细胞图像（spatial single-cell images）中进行半监督聚类（semi-supervised clustering）和表达簇投影（projection of expression clusters）的方法。图1概述了Cellar的工作流程。过去一年，HuBMAP的成员对Cellar进行了测试，并将其用于注释来自不同器官、平台和模式的多个单细胞数据集。
fig2

图1：Cellar的工作流程。a–c为预处理（可选）：Cellar可以根据表达的基因过滤细胞，根据很少表达的基因的数量过滤基因。接下来，输入被normalization。d和e：降维和可视化。作为Cellar的一部分，实现了几种降维方法。然后，通过运行另一种（可能是相同的）降维方法，将降维后的数据可视化。f–i为聚类：Cell支持几种无监督和半监督聚类方法。它还实现了有监督的 label transfer 方法。j–l为细胞类型分配：Cellar支持使用多个功能注释数据库（functional annotation databases）来分配细胞类型。

结果

Analysis of scRNA-seq data

我们使用Cellar分析了11个HuBMAP的seq数据集（10x genomics），平均有来自5个不同组织（肾脏、心脏、脾脏、胸腺、淋巴结）的7500个细胞，所有这些都可以在Cellar中获得。Cellar首先通过去除不可靠的细胞和低计数基因来进行质量控制。根据用户标准应用额外的normalization和缩放scaling。然后，Cellar对数据的低维表示进行聚类，并为可视化进一步降低维度。我们通过分析含有5273个细胞的脾脏数据集（Cellar ID:HBMP3spleen-CC2）来证明这一基本pipeline。我们使用PCA，然后使用UMAP进行降维，使用Leiden算法进行聚类，得到总共16个簇（补充图1a）。
fig3

补充图1a

对于每一个簇，Cellar都识别出了top差异基因（top differential genes）。使用前500个差异基因，functional enrichment analysis 功能富集分析（GO、KEGG、MSigDB）将簇0确定为B细胞（例如，“B-Cell Activation”（ $q$ value=0）和“B-Cell Receptor Signaling Pathway”（ $q$ value=0）分别是GO和KEGG的前两类）。通过可视化两种已知B细胞 markers CD79A和TNFRSF13C的同时表达，进一步支持了这一任务。

功能富集分析就是把一个基因列表中，具有相似功能的基因放到一起，并和生物学表型关联起来；

GO与KEGG
为了解决将基因按照功能进行分类的问题，科学家们开发了很多基因功能注释数据库。其中比较有名的就是Gene Ontology（基因本体论，GO）和Kyoto Encyclopedia of Genes and Genomes（京都基因与基因组百科全书，KEGG）。

其中，GO是基因本体论联合会建立的一个数据库，旨在建立一个适用于各种物种的、对基因和蛋白功能进行限定和描述的，并能随着研究不断深入而更新的语义词汇标准。GO注释分为三大类，分别是：分子生物学功能（Molecular Function，MF）、生物学过程（Biological Process，BP）和细胞学组分（Cellular Components，CC），通过这三个功能大类，对一个基因的功能进行多方面的限定和描述。

而KEGG，大多数听说过KEGG的人都会把它当做一个基因通路（Pathway）的数据库，其实功能远不止于此。KEGG是一个整合了基因组、化学和系统功能信息的综合数据库。KEGG下属4个大类和17和子数据库，而其中有一个数据库叫做KEGG Pathway，专门存储不同物种中基因通路的信息，也是用的最多的一个，所以，久而久之，KEGG就被大家当做是一个通路数据库了。

除了无监督聚类外，Cellar还实现了基于参考数据集的有监督分配方法。这些可以直接利用双模式和在Cellar中实现的其他方法。例如，这种形式的分配可以与Cell的半监督聚类选项结合使用，以纠正标签转移过程中的噪声。为了说明这种用法，我们应用Scanpy的 Ingest function（可在Cellar中使用）来集成两个专家注释的脾脏数据集（Cellar ID:HBMP2-Splege-2和HBMP3-Splege-CC3）。我们使用HBMP3-CC3作为GT，并将标签从它转移到HBMP2-2。然后，我们将标签转移的结果与HBMP2-2的GT注释进行比较，并观察到调整后的 rand score（ARI）为0.39。相比之下，在HBMP2-2上运行Leiden会导致更低的ARI分数0.27。然后，我们使用Leiden的半监督自适应算法来改进标签转移的结果，其中选择噪声最小的簇作为约束，并且在算法迭代期间不允许改变。这导致了更好的ARI得分0.66，证明了标签转移和半监督聚类的好处。这些结果如补充图2所示。
fig4

补充图2：a：HBMP2-2脾脏数据集的GT注释。b：使用默认参数（ARI:0.27）运行vanilla Leiden后的簇分配。c：使用Scanpy Ingest（ARI:0.39）从HBMP3-CC3传播标签。在运行半监督Leiden之后，在从c（ARI:0.66）传输的标签上使用默认参数进行d里的簇分配。

Analysis of scATAC-seq data

虽然scRNA-Seq是目前使用最广泛的数据模式，但其他几种分子数据类型也在单细胞水平上进行分析。为了说明Cellar对此类数据的使用，我们使用它来注释scATACseq。Cellar可以通过两种不同的方式处理scATAC-seq数据：cell-by-gene和cell-by-cistopic。前者基于与所有基因附近区域相关的开放染色质可及性，而后者则依赖于顺式Topic，它使用潜在的Dirichlet分配来模拟顺式调节。生成的逐基因细胞或逐细胞顺题矩阵用于下游分析，如可视化和聚类。我们使用Cellar对scATAC-seq数据集进行注释，该数据集使用逐基因表示法分析外周血单个核细胞（Cellar ID:PBMC 10k Cell-By-Gene）。结果如补充图3所示。对簇0和簇4的DE分析确定了自然杀伤（NK）细胞的KLRD1 marker。
fig5

补充图3：b是KLRD1标记基因的表达水平。该基因是自然杀伤（NK）细胞的已知标记，在簇0和簇4中高度表达。

Analysis of spatial transcriptomics data (CODEX)

除了测序分析，最近的成像分析还可以提供单细胞水平上基因或蛋白质表达的信息。Cellar可以通过提供表达式簇和空间组织的并排视图来分析此类数据。为了说明这一点，我们分析CO-Detection by indEXing (CODEX)空间蛋白质组学数据。我们使用了一个包含46840个细胞的淋巴结数据集（Cellar ID:19-003 lymph node淋巴结R2）。聚类结果如图2所示，以及带有投影聚类注释的这些细胞的对应图块。鉴于该数据集中分析的蛋白质数量较少（19），并非所有的聚类都可以被分配到唯一的类型，尽管有几个聚类是基于Cellar中的DE基因分析分配的。Cellar将聚类和空间图像中的细胞颜色匹配起来，从而更容易识别特定的组织及其与分析细胞类型的关系。图2中的空间图显示B细胞紧密地聚集在一起，并被T细胞和淋巴中的其他细胞类型包围。B细胞簇也包含增殖细胞（proliferating cells）。
fig6

图2：a：lymph node CODEX数据集的UMAP可视化表示，包含46840个细胞，通过leiden聚类。b：在 spatial CODEX image上的分配投影，可在Cellar中并排显示。簇分配是从a复制的。考虑到只测量了几十个蛋白质水平，并不是所有的聚类都能被分配到唯一的细胞类型，尽管有几个聚类是根据Cellar中的差异基因分析分配的。淋巴中的B细胞簇被T细胞和其他类型的细胞包围。B细胞簇也包含增殖细胞（proliferating cells）。

Joint analysis of multiple modalities

最后，我们使用Cellar联合分析来自两种不同模式的数据。为此，我们使用了SNARE-seq肾脏数据集，该数据集对31758个细胞的转录组和染色质可及性进行了分析（Cellar ID: kidney SNARE ATAC/RNA 20201005）。在这里，我们首先在染色质形态上运行cisTopic，并通过在推断的cis调节主题上运行Leiden来确定簇分配（图3a）。我们使用这些标签来可视化图3b中的表达数据。这可以通过Cellar的双模式（dual mode）轻松实现，该模式允许基于细胞ID的标签从一种模式转移到另一种模式。Cellar鉴定了差异基因，我们用这些基因来绘制细胞类型图。例如，根据已知标记（SLC5A12， $p$ value=0）和GO term analysis（“Apical Plasma Membrane”， $p$ value=1e-4）分配簇1，这表明存在近端小管细胞（Proximal Tubule Cells）。
fig7

图3：a：肾脏SNAREseq数据集31758个细胞染色质形态的UMAP图。首先，我们通过运行cistopic获得一个由cistopic矩阵生成的单元，然后通过Leiden聚类定义聚类。b：表达矩阵（expression matrix）的相应UMAP图，以及从（图a）Cellar的双模式复制的集群分配，允许基于细胞ID的标签从一种模式转移到另一种模式。