SEACells：元细胞分析

元细胞是从单细胞测序数据中衍生的细胞分组，代表高度精细的不同细胞状态。在这里，作者介绍了单细胞细胞状态聚集 (SEACells)，这是一种用于识别元细胞的算法，它克服了单细胞数据的稀疏性，同时保留了传统细胞聚类所掩盖的异质性。SEACells 在识别 RNA 和ATAC模态分析中的全面、紧凑且分离良好的元细胞表现上优于现有算法，作者展示了如何使用 SEACells 来改善gene-peak关联、计算 ATAC 基因分数并推断分化过程中关键regulators的活动。元细胞级分析可扩展到大型数据集。作者利用元细胞揭示造血分化过程中染色质景观的表达动态，并唯一地识别与 2019 年冠状病毒病 (COVID-19) 患者群体中的疾病发作和严重程度相关的 CD4 T 细胞分化和活化状态。

来自：SEACells infers transcriptional and epigenomic cellular states from single-cell genomics data

方法概述

SEACells 力图将单个细胞聚合成代表不同细胞状态的元细胞，这种方式与数据模态无关。使用计数矩阵作为输入，它为每个元细胞提供每个细胞的权重、每个元细胞的硬分配以及每个元细胞的聚合计数作为输出。SEACells可以捕获数据中的全部细胞状态，包括较罕见的状态。基于几个关键假设建立了 SEACells：(1) 单细胞分析数据可以用低维流形（表型流形）来近似；(2) 观察到的细胞间变异性大部分是由于采样不完整造成的；(3) 大多数细胞可以分配到一组有限的细胞状态，每个状态都以不同的活性基因调控程序组合为特征。

SEACells 利用基于图的流形学习算法，该算法已被证明能够忠实而稳健地捕捉单细胞基因组数据中的细胞状态图。该算法首先构建一个表示表型流形的最近邻图。然后，它应用原型分析迭代细化元细胞。最后，它将计数聚合到一组输出元细胞中。流形构造针对每种数据模态定制，之后算法可以以与数据类型无关的方式进行。作者使用来自早期人类造血的 CD34+ 细胞来演示方法（图 1）。使用最小-最大采样进行初始化，它识别出一组均匀分布在表型流形中的代表性细胞状态（图 1e），并且特别擅长处理密度差异，从而确保捕获稀有状态。这些采样细胞是waypoints（每种细胞类型有多个），它们在邻居图中定义清晰的结构；然而，细胞状态本身仍然有些分散（图 1f）。

为了细化元细胞，作者采用了核原型分析-kernel archetypal analysis（图 1g）。原型分析-archetypal analysis是一种稳健的矩阵分解技术，已应用于数据矩阵，以识别细胞表型空间边界处的极端细胞状态。相反，作者将原型分析应用于细胞间相似性核矩阵。该核将细胞投影到更高维空间中，其中两个细胞只有当它们共享邻居并且与共享邻居的距离相似时才相似。这种转换所施加的更严格的相似性条件将高度相似的细胞投射到微小的簇中，使得核空间中的原型成为每个独特细胞状态的良好代表。因此，核原型分析将细胞划分为高度相似的细胞的紧密簇（图 1g），沿细胞-细胞相似性矩阵的对角线赋予紧密的块，这些块代表不同的细胞状态（图 1h）。
fig1

图1：a，6,800 个 CD34+ 分类的造血干细胞和祖细胞 (HSPC) 的 scRNA-seq UMAP。细胞按簇着色。b，每个簇的轮廓图突出显示密度并表明每个簇内存在多种细胞状态。插图：基因-基因协方差矩阵显示每个状态包含多个不同的基因表达程序。c，左：UMAP 表示 MEP（巨核细胞-红细胞祖细胞）簇。右：当 MEP 簇根据发育进程分为三个大小相等的箱体时（顶部，G1 到 G3），它反映了 GATA2（已知的 MEP 谱系驱动因素）的推断表达（底部）。d，覆盖图显示所有 MEP（顶部）、单个 MEP 细胞（底部）和 c 中的三个箱体中的 GATA2 可访问性。右：相应细胞中的 GATA2 表达。突出显示的峰值展示了可访问性动态如何跟踪表达动态。有关动态的信息在簇级别被掩盖，而单个细胞中的峰值识别太嘈杂。e，UMAP 与 a 中的一样，按细胞类型着色。用于元细胞识别的 SEACells 算法由waypoints（大红色圆圈）初始化，waypoints是采样以均匀覆盖表型景观的细胞子集。f，使用自适应高斯核计算的细胞间亲和力矩阵。g，核原型分析示意图。核矩阵分解为原型矩阵 B 和嵌入矩阵 A。根据矩阵 A 上的逐列最大值来识别元细胞成员资格。插图：核原型分析将细胞划分为高度相似的细胞簇，使其非常适合识别稳健的细胞状态。h，左：来自 f 的细胞-细胞亲和力矩阵，但按元细胞分配排序（到这里，体现出用于发现罕见簇的潜力）。右：覆盖来自 e 的 UMAP 的轮廓图，突出显示元细胞的分布；细胞和轮廓通过元细胞分配来着色。

实验

SEACells元细胞代表准确和稳健的细胞状态

作者首先在 10x Genomics 的外周血单核细胞 (PBMC) 公共多组学 (同时进行 scRNA-seq 和 ATAC-seq) 数据集上评估了 SEACells 的性能，这是一个经过充分研究的系统，具有不同的细胞群体。可以发现 SEACells 元细胞全面，在细胞类型中分布良好 (图 2a、b)。

元细胞有助于克服稀疏性，数据稀疏性在 scATAC-seq 中非常严重。作者发现每个 SEACells 元细胞都比单个细胞提供了更完整的分子表征——例如，通过揭示主要细胞类型已知marker基因的可及性。元细胞（但不是大多数单个细胞）的可及性和表达可以准确区分淋巴亚群（图 2c ）。因此，元细胞的粒度足以区分细胞类型内的状态；并且可以使用经典免疫marker进行查询。
fig2

图2：a、(i) 使用 RNA 数据得出的 10x Genomics 多组学数据集中人类 PBMCs 的 UMAP，突出显示细胞类型和 SEACells 元细胞。(ii) RNA 模态下每种细胞类型的元细胞分布。(iii) 细胞类型纯度分布（每种元细胞中最具代表性的细胞类型的频率）。高纯度代表更准确的元细胞。b、与 a 中一样，使用来自多组学数据集的 ATAC 数据得出人类 PBMCs 的 UMAP、元细胞和细胞类型纯度分布。c、CD4 和 CD8A 的元细胞可及性 (i) 和表达 (iii) 准确区分 CD8（绿色）和 CD4（橙色）T 细胞。TYROBP 和 CD8A 的元细胞可及性 (ii) 和表达 (iv) 区分 NK（棕色）和 CD8（绿色）T 细胞。插图：相应的单细胞可及性太稀疏，无法实现相同的区分。

元细胞的意义：获得一个新的细胞数据，这个数据不稀疏，并且保留了原有的异质性

SEACell促进的调控推断

可以通过识别 ATAC-seq 读取计数峰值内的假定转录因子 (TF) 结合基序来推断基因调控，这些基序代表开放或可访问的染色质区域。scATAC-seq 提供了许多观测结果（细胞），能以精细分辨率推断更复杂的基因调控模型，但数据稀疏性严重限制了其实用性。SEACells 元细胞能克服稀疏性，从而实现各种基因调控推断任务。

典型的 SEACells 元细胞包含 120 万个读数，与单个细胞中的 25,000 个读数相比有显著改善，但仍远少于典型bulk样本中的 5000 万个读数。为了提高 ATAC 峰值调用中的信噪比，作者利用了 ATAC-seq 片段长度分布（一种表示方式），其中第一和第二模式分别代表无核小体 (NFR) 片段（可能富含 TF 结合事件）和核小体。

调控推断的下一个任务是将每个基因与调控它的元素关联起来。跨细胞的可及性和表达之间的相关性已被用于预测scRNA-seq 和 scATAC-seq 中调控每个基因的峰集，但数据稀疏阻碍了单细胞水平上的稳健性。使用来自 CD34+ 骨髓 ATAC 数据中的 SEACells 元细胞，作者计算了核心造血基因集中每个基因 ±100 kb 范围内每个 NFR 峰的基因表达与可及性之间的相关性。使用 ATAC 元细胞的最相关峰的可及性跟踪基因表达，与单细胞数据相比有了显着改善（图 3a）。例如，关键的红系谱系调节剂TAL1的峰值可及性与元细胞中表达之间的相关性为0.82，而在单细胞水平上的相关性为0.03（图3a）。

对于关键的红细胞因子 GATA2，单细胞数据仅恢复了使用元细胞检测到的 11 个关联中的2个（图 3b）。为了系统地探索预测的peak-gene关联的准确性，作者通过汇总所有显着相关峰的可访问性并将其与基因表达进行比较来计算基因得分。SEACells 基因得分的相关性明显优于使用所有相关峰的汇总得出的得分。因此，SEACells 元细胞清楚地识别了与表达显着相关并可能调节相应基因的顺式调控元件。
fig3ab

图3a：ATAC 元细胞（顶部）或单细胞（底部）基因表达与 TAL1（红细胞）、MPO（髓系）和 IRF（树突状细胞）marker基因中最相关峰的可及性之间的 Spearman 相关性，基于 CD34 多组学数据计算。每个元细胞和单细胞根据细胞类型着色。注意：不同模态也可以计算Spearman 相关性（用于衡量两者的单调性）。
图3b：使用 NFR（顶部）或所有 ATAC（底部）片段绘制 HSC、MEP 和红细胞 (Ery) 中的红细胞因子 GATA2 的可及性景观。将染色质可及性分析限制在 NFR 片段可提高峰值分辨率和调节元件与基因的关联。弧线由峰值-基因 Spearman 相关性着色（右侧颜色值介于 0 和 1 之间），使用 SEACells ATAC 元细胞确定。

准确推断peak-gene关联可以促进将ATAC数据转换为gene score。

方法

SEACells

SEACells 是一种从单细胞数据定义元细胞的算法。SEACells 算法假设生物系统由明确定义且有限的细胞状态集组成。观察到的单细胞数据被认为是这些细胞状态的稀疏且嘈杂的测量值（当前最先进的单细胞测量技术只能捕获 <10% 的转录本或 <5% 的开放染色质区域）。尽管噪声程度很高，但由于定义细胞状态的基因表达模式和调控机制，从相同状态中采样的细胞被认为具有密切相关的表型。SEACells 旨在将密切相关的细胞聚集到代表它们的元细胞中，从而克服单细胞数据的稀疏性。此外，由于稀疏性，scATAC-seq 数据的实用性特别有限。SEACells 元细胞还提供了一种可扩展的表示，可以有效处理大规模单细胞数据。尽管聚类被广泛用于克服稀疏性，但聚类掩盖了数据中存在的大量异质性。SEACells 元细胞实现了保留异质性同时克服单细胞数据稀疏性。

SEACells 的输入包括 (1) 原始计数矩阵（例如，RNA 的转录本计数、ATAC 的peak或bins计数）；(2) 使用模态适当的预处理（例如 RNA 的主成分分析 PCA）得出的数据的低维表示；以及 (3) 要识别的元细胞数量。作为下游分析的输出，SEACells 生成代表元细胞的细胞分组、聚合的元细胞计数矩阵和代表高度相关细胞组的软分配。该算法可在 https://github.com/dpeerlab/SEACells 免费获取。

SEACells大致包含5个步骤：

使用在低维嵌入空间中计算的细胞之间的欧几里得距离构建KNN图，以表示表型流形（phenotypic manifold）
使用最近邻图导出细胞间相似性的亲和矩阵。使用自适应高斯核将graph中的距离转换为相似性，以解释表型流形中的细胞密度差异。亲和矩阵或核矩阵编码细胞之间的非线性关系
使用核矩阵作为原型分析的输入(图1g)。虽然原型分析通常应用于数据矩阵，但作者将其应用于核矩阵，该核矩阵将细胞划分为高度相似的细胞簇，并能够表征整个表型流形，使其非常适合识别稳健的细胞状态。原型分析将数据分解为原型矩阵，原型矩阵包括代表表型流形上细胞状态的细胞线性组合，以及将单个细胞重构为原型线性组合的隶属度矩阵(图1g)。该方法对数据进行分区，使细胞-细胞相似性矩阵沿对角线具有紧密的块结构;每个分区是一组最能代表细胞状态并定义一个metacell的cells。元细胞的数量被指定为原型分析的输入
将通过原型分析确定的分组标记为SEACells元细胞，并相应地汇总单个细胞原始计数，以导出metacell-by-feature矩阵
归一化计数矩阵，可用于所有下游分析任务，如聚类、可视化、数据整合、轨迹推断和基于ATAC-seq的调控推理

Toolkit for ATAC

目前已经开发出大量强大的工具来解释来自bulk ATAC-seq 数据的开放染色质数据。然而，由于稀疏性，它们不能直接应用于单细胞数据。SEACells 元细胞是紧密相关细胞的聚集体，因此在忠实保留数据的异质性和结构的同时，稀疏性大大降低。在这里，作者描述了一个适用于 scATAC-seq 数据的强大工具包，该工具包改编自bulk数据分析工具。

Peak calling

使用ArchR执行。ArchR首先对单细胞数据进行聚类，并使用MACS2 peak caller分别识别每个簇的峰值。然后将每个峰的大小调整为500个碱基，峰顶位于中心，并合并跨不同簇的重叠峰。合并的峰值再次调整为500个碱基。

ATAC-seq数据提供了跨越TF结合区域和非抑制区域核小体的开放染色质区域的概况。ATAC-seq数据的片段大小分布包含反映该信息多样性的特征模式。因为第一种模式代表NFR，作者修改了ArchR管道，只使用NFR片段(片段长度< 147)来识别峰值，而不是使用所有片段的默认值。这种变化导致对调控元件的识别更加敏感。

Peak-gene关联和gene score

尽管NFR片段的使用提高了被称为峰的灵敏度，但并非所有被识别的峰都代表调节基因表达的TF结合事件。已有研究提出利用整合ATAC和RNA数据的峰可达性与基因表达的相关性来识别可能调控基因表达的峰。SEACells元细胞为计算这些关联提供了理想的分辨率，当使用稀疏的单细胞数据计算时，这些关联是不可靠的。作者使用ATAC模态鉴定的元细胞来构建峰基因关联。

作者采用Ma等人的程序来确定显著的峰-基因关联（Chromatin potential identified by shared single-cell profiling of RNA and chromatin）。对于每个基因，使用归一化的元细胞表达和归一化的ATAC可及性，计算基因上游100kb和下游100kb内的每个峰值的Pearson相关性。为了评估峰-基因相关性的重要性，采样了100个峰的经验背景，这些峰与GC含量和所考虑的峰的可及性相匹配。根据气相色谱含量和样品经验背景可及性，将各峰分别分成100个bins。任何名义 $P < 1 × 10^{−1}$ 的峰都被认为是显著的峰-基因关联。使用NFR片段鉴定的峰用于本分析。与基因相关的所有峰的总可达性用于确定元细胞基因评分。

推断TF活性

为了利用峰基因关联，作者提供了一种简单的基因调控网络(GRN)方法来推断TF活性，用于识别与不同细胞类型相关的关键TF。