单细胞分析(Signac)： PBMC scATAC-seq 聚类

引言

在本教学指南中，我们将探讨由10x Genomics公司提供的人类外周血单核细胞（PBMCs）的单细胞ATAC-seq数据集。

加载包

首先加载 Signac、Seurat 和我们将用于分析人类数据的其他一些包。

if (!requireNamespace("EnsDb.Hsapiens.v75", quietly = TRUE))
    BiocManager::install("EnsDb.Hsapiens.v75")

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)

library(ggplot2)
library(patchwork)

归一化和线性降维

在数据处理中，我们首先进行标准化处理，采用的是词频-逆文档频率（TF-IDF）方法。这个过程分为两步：首先在细胞层面进行标准化，以消除不同细胞测序深度带来的偏差；其次在峰值层面进行标准化，目的是让那些较为罕见的峰值获得更高的权重。

接下来是特征选择环节。由于单细胞ATAC-seq数据的动态范围较窄，这使得像单细胞RNA测序那样进行变量特征选择变得困难。因此，我们可以选择仅利用特征（峰值）的前n%来进行降维处理，或者通过FindTopFeatures()函数排除在少于n个细胞中检测到的特征。在本例中，我们将采用所有特征进行分析，尽管我们发现仅使用特征的一个子集（例如，将min.cutoff设置为‘q75’以选取所有峰值的前25%）也能获得非常相似的分析结果，并且计算速度更快。该函数还会自动将用于降维的特征标记为Seurat对象的VariableFeatures()。

然后，我们对选定特征的TD-IDF矩阵执行奇异值分解（SVD）。这一步骤将为我们提供数据的降维表示，对于熟悉单细胞RNA测序的用户来说，这个过程类似于主成分分析（PCA）的输出结果。

将TF-IDF标准化与SVD相结合的步骤被称为潜在语义索引（LSI），这种方法最初由Cusanovich等人在2015年引入，用于单细胞ATAC-seq数据的分析。

pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)

通常，潜在语义索引（LSI）的第一个主成分反映的是测序的深度（即技术层面的变异），而非生物学上的变化。如果确实如此，那么在后续分析中应该将这个成分剔除。我们可以通过DepthCor()函数来评估每个LSI主成分与测序深度之间的相关性。

DepthCor(pbmc)

我们观察到，第一个潜在语义索引（LSI）主成分与细胞的总计数数之间具有极强的相关性。因此，在后续分析中，我们将排除这个主成分，目的是避免基于细胞的总测序深度来进行细胞分组，而应依据细胞在特定细胞类型峰值上的可访问性模式来进行分组。

非线性降维和聚类

细胞已经嵌入到一个低维度的空间里，这使得我们可以采用分析单细胞RNA测序数据(scRNA-seq)时常用的方法来进行基于图的聚类分析，并通过非线性降维技术进行数据可视化。RunUMAP()、FindNeighbors()和FindClusters()这些函数都是Seurat软件包提供的。

pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
DimPlot(object = pbmc, label = TRUE) + NoLegend()

创建基因活动矩阵

通过UMAP的可视化技术，我们可以清晰地看到人类血液中存在许多不同的细胞群体。如果您对单细胞RNA测序(scRNA-seq)分析的外周血单个核细胞(PBMC)有所了解，您可能会在单细胞ATAC-seq数据中辨认出特定的骨髓细胞和淋巴细胞群体。然而，由于我们对非编码基因区域的功能作用了解较少，因此在单细胞ATAC-seq数据中对聚类进行注释和解读会更加困难。

我们可以通过评估与基因相关的染色质可及性来尝试量化基因组中每个基因的活性，并基于单细胞ATAC-seq数据开发出一种新的基因活性检测方法。这里，我们将采用一种简便的方法，即累加与基因本体及其启动子区域相交的片段（我们也推荐使用Cicero工具，它能够实现类似的目标，并且我们提供了一个指南，展示如何在Signac工作流中运行Cicero）。

为了构建一个基因活性矩阵，我们首先提取基因的坐标，并将其扩展到包括上游2kb的区域（因为启动子的可及性通常与基因表达相关）。接着，我们利用FeatureMatrix()函数计算每个细胞中对应这些区域的片段数量。这些步骤都可以通过GeneActivity()函数自动完成。

gene.activities <- GeneActivity(pbmc)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(pbmc$nCount_RNA)
)

现在我们可以通过展示标准标记基因的活性来帮助我们解读ATAC-seq的聚类结果。需要注意的是，这些活性数据的噪声会远大于单细胞RNA测序(scRNA-seq)的测量结果。这主要是因为它们来源于稀疏的染色质数据，并且基于一个假设——即基因体或启动子的可及性与基因表达之间存在普遍的对应关系，但这种假设并不总是成立。尽管如此，我们还是可以通过这些基因活性模式开始识别出单核细胞、B细胞、T细胞和自然杀伤细胞(NK细胞)的群体。然而，仅凭监督分析来进一步细分这些细胞类型是有难度的。

DefaultAssay(pbmc) <- 'RNA'

FeaturePlot(
  object = pbmc,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)