基于Signac的scATAC-seq流程学习

在进行染色质数据的预处理时，Signac会使用来自CellRanger生成的两个相关的数据。峰/细胞矩阵（Peak/Cell matrix）：这一矩阵类似于在单细胞RNA测序分析中使用的基因表达计数矩阵。不过，不同于以“基因”为行，这里的每一行代表的是基因组中的一个区域（即一个峰，peak），该区域被预测为开放染色质区域。矩阵中的每个数值表示：在对应的单细胞（由唯一条形码标识）中，落在该峰区域内的Tn5转座整合位点的数量。片段文件（Fragment file）：该文件包含了所有单细胞中所有唯一片段（fragments）的完整列表。它的体积通常比峰/细胞矩阵大得多，读取和处理速度较慢，并且是存储在磁盘上（而不是内存中）进行操作的。保留该文件的优点在于：它不仅包含了那些映射到峰区域的片段，还包含了每个单细胞关联的全部片段信息。

原始数据下载

wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5
wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv
wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz
wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz.tbi

分析流程：

1.导入

rm(list=ls())
options(timeout = 6000000000)  
library(Signac)
library(Seurat)
library(GenomicRanges)
library(ggplot2)
library(patchwork)

2.数据预处理

# 使用cellranger-atac生成的峰/细胞矩阵和细胞元数据创建一个Seurat对象，并将磁盘上的片段文件路径存储在Seurat对象
counts <- Read10X_h5(filename = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv",
  header = TRUE,
  row.names = 1
)

# 如果没有h5文件也有解决办法：https:///signac/articles/pbmc_vignette
chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz",
  min.cells = 10,
  min.features = 200
)

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

pbmc[['peaks']]

granges(pbmc)

# 移除与染色体支架（例如(KI270713.1)）或其他序列对应的特征，而不是(22+2)标准染色体。
peaks.keep <- seqnames(granges(pbmc)) %in% standardChromosomes(granges(pbmc))
pbmc <- pbmc[as.vector(peaks.keep), ]

3.向人类基因组pbmc对象添加基因注释

# 每个基因组组装都会发布多个补丁，在处理映射数据时，建议使用与映射操作相同的组装补丁注释；
library(AnnotationHub)
ah <- AnnotationHub()

# 在“注释库”网站上搜索“Ensembl 98 EnsDb”（人类基因组注释数据库）
query(ah, "EnsDb.Hsapiens.v98")
#AnnotationHub with 1 record
# snapshotDate(): 2025-04-08
# names(): AH75011
# $dataprovider: Ensembl
# $species: Homo sapiens
# $rdataclass: EnsDb
# $rdatadateadded: 2019-05-02
# $title: Ensembl 98 EnsDb for Homo sapiens
# $description: Gene and protein annotations for Homo sapiens based on Ensembl version 98.
# $taxonomyid: 9606
# $genome: GRCh38
# $sourcetype: ensembl
# $sourceurl: http://www.
# $sourcesize: NA
# $tags: c("98", "AHEnsDbs", "Annotation", "EnsDb", "Ensembl", "Gene", "Protein",
#   "Transcript") 

# 使用object[["AH75011"]]来提取记录
ensdb_v98 <- ah[["AH75011"]]

# 从EnsDb中提取基因注释信息
annotations <- GetGRangesFromEnsDb(ensdb = ensdb_v98)
# 由于数据已映射至hg38版本，故改为UCSC格式。
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"
# 将基因信息添加到该对象中
Annotation(pbmc) <- annotations

4.计算质控指标

计算QC指标，像scRNA-seq一样需要进行质控
核小体带型：DNA片段大小的直方图（通过成对末端测序读取确定）应表现出强烈的核小体带型，对应于环绕在单个核小体上的DNA长度。每个细胞进行一次计算，并量化单核小体片段与无核小体片段的大致比率（存储为 nucleosome_signal）。
转录起始位点（TSS）富集评分：ENCODE项目根据以TSS为中心的片段与TSS邻近区域片段的比例定义了ATAC-seq靶向评分（参见 https://www./data-standards/terms/）。质量较差的ATAC-seq实验通常会具有低的TSS富集评分。可以使用TSSEnrichment() 函数为每个细胞计算该指标，结果存储在元数据中，列名为 TSS.enrichment。
峰中的片段总数：这是细胞测序深度/复杂度的一个指标。具有非常少的读取的细胞可能需要因测序深度过低而被排除。具有极高读取数的细胞可能代表双细胞、细胞核聚集或其他伪影。
峰中片段的比例：表示所有片段中落在 ATAC-seq 峰内的片段比例。具有低值（即<15-20%）的细胞通常代表低质量细胞或应被移除的技术伪影。请注意，这个值可能对所使用的峰集很敏感。
在基因组黑名单区域中的读取比例：ENCODE项目提供了一份黑名单区域的清单，这些区域表示通常与伪信号相关的读取。高比例的读取映射到这些区域（相对于映射到峰值的读取）的细胞通常代表技术伪影，应被剔除。ENCODE为人类（hg19 和 hg38）、小鼠（mm9 和 mm10）、果蝇（dm3 和 dm6）以及线虫（ce10 和 ce11）提供了黑名单区域，这些区域已包含在 Signac包中。FractionCountsInRegion函数可以用来计算每个细胞在给定区域集内的所有计数的比例。我们可以使用这个函数和黑名单区域来找到每个细胞内的黑名单计数比例

# 计算每个细胞的核小体信号评分。
pbmc <- NucleosomeSignal(object = pbmc)

# 计算每个细胞的转录起始位点（TSS）富集评分。
pbmc <- TSSEnrichment(object = pbmc)

# 添加峰值中的读取比例。
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100

# 添加黑名单比例。
pbmc$blacklist_ratio <- FractionCountsInRegion(
  object = pbmc, 
  assay = 'peaks',
  regions = blacklist_hg38_unified
)

5.DensityScatter可视化

DensityScatter(pbmc, x = 'nCount_peaks', y = 'TSS.enrichment', 
               log_x = TRUE, quantiles = TRUE)
# 还可以查看所有细胞的片段长度周期性，并按高或低核小体信号强度对细胞进行分组。可以看到，在单核小体/无核小体比率上表现为极端值的细胞（基于上述图表）具有不同的核小体条带模式。其余细胞表现出典型成功的ATAC-seq实验的模式。
pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')

还可以查看所有细胞的片段长度周期性，并按高或低核小体信号强度对细胞进行分组。可以看到，在单核小体/无核小体比率上表现为极端值的细胞具有不同的核小体条带模式。其余细胞表现出典型成功的ATAC-seq实验的模式。

6.小提琴展示QC指标分布情况并过滤

VlnPlot(
  object = pbmc,
  features = c('nCount_peaks', 'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal', 'pct_reads_in_peaks'),
  pt.size = 0.1,
  ncol = 5
)

pbmc <- subset(
  x = pbmc,
  subset = nCount_peaks > 9000 &
    nCount_peaks < 100000 &
    pct_reads_in_peaks > 40 &
    blacklist_ratio < 0.01 &
    nucleosome_signal < 4 &
    TSS.enrichment > 4
)
pbmc

7.归一化

归一化：Signac执行词频-逆文档频率（TF-IDF）归一化。这是一种两步归一化程序，既能在细胞间进行归一化以校正细胞测序深度的差异，也能在峰之间进行归一化，以赋予更罕见的峰更高的值。
特征选择：scATAC-seq，数据的动态范围较低，这使得进行变量特征选择变得具有挑战性，就像对scRNA-seq所做的那样。相反，可以选择仅使用前n%的特征（峰）进行降维，或者使用 FindTopFeatures()函数移除在少于n个细胞中出现的特征。在这里，将使用所有特征，尽管使用特征子集（尝试将 min.cutoff 设置为“q75”以使用所有峰的25%）时看到了非常相似的结果，并且运行时间更快。用于降维的特征将自动通过此函数设置为Seurat对象中的 VariableFeatures() 。
降维：接下来对TD-IDF，矩阵运行奇异值分解（SVD），使用上述选择的特征（峰）。这将返回对象的降维表示（对于更熟悉 scRNA-seq 的用户，可以将其视为 PCA 输出的类似物）。

pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)
# 第一个LSI组件通常捕捉的是测序深度（技术变异），而不是生物变异。如果是这种情况，则应从下游分析中移除该组件。可以使用 `DepthCor()` 函数评估每个 LSI 组件与测序深度之间的相关性：
DepthCor(pbmc)

这里看到，第一个 LSI 成分与细胞的计数的总数之间存在非常强的相关性。因此将不使用这个成分进行后续步骤，因为不想根据测序深度将细胞分组，而是想根据它们在细胞类型特异性峰上的可及性模式进行分组。

8.非线性降维和聚类

# 现在细胞已嵌入低维空间，可以使用常用于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()

9.创建基因活性矩阵

UMAP可视化揭示了人类血液中存在多个细胞群。如果熟悉PBMC的scRNA-seq分析，可能会在scATAC-seq数据中识别出某些髓系和淋巴系细胞群。然而，在scATAC-seq数据中注释和解释簇更具挑战性，因为关于非编码基因组区域的功能作用知之甚少，而关于蛋白质编码区域（基因）的功能作用则了解得更多。
可以通过评估与基因相关的染色质可及性来量化基因组中每个基因的活动，并从scATAC-seq数据中创建一个新的基因活性检测方法。在这里，将采用一种简单的方法，即对与基因体和启动子区域相交的片段进行求和（还推荐探索Cicero工具，该工具可以实现类似的目标，在此提供了一个 vignette，展示了如何在 Signac 工作流程中运行Cicero）。
要创建基因活性矩阵，提取基因坐标并将其扩展以包含上游2kb区域（因为启动子可及性与基因表达通常相关）。然后，使用 FeatureMatrix()函数统计每个细胞中映射到这些区域的片段数量。这些步骤由 GeneActivity() 函数自动完成：

gene.activities <- GeneActivity(pbmc)
# 将基因活性矩阵添加到 Seurat 对象中作为一个新测定，并进行归一化。
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(pbmc$nCount_RNA)
)

# 现在可以可视化典型标记基因的活性，以帮助解释ATAC-seq簇。请注意，这些活性将比 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
)

10.与scRNA-seq数据整合

# Load the pre-processed scRNA-seq data for PBMCs
pbmc_rna <- readRDS("pbmc_10k_v3.rds") # 来源于官网
pbmc_rna <- UpdateSeuratObject(pbmc_rna)
# 整合映射
transfer.anchors <- FindTransferAnchors(
  reference = pbmc_rna,
  query = pbmc,
  reduction = 'cca'
)

predicted.labels <- TransferData(
  anchorset = transfer.anchors,
  refdata = pbmc_rna$celltype,
  weight.reduction = pbmc[['lsi']],
  dims = 2:30
)

pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)
plot1 <- DimPlot(
  object = pbmc_rna,
  group.by = 'celltype',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')

plot2 <- DimPlot(
  object = pbmc,
  group.by = 'predicted.id',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')

plot1 + plot2

可以看到，基于scRNA的分类与使用scATAC-seq数据计算的UMAP可视化结果是一致的。然而，在scATAC-seq 数据集中，有一小部分细胞被预测为血小板。这是令人差异的，因为血小板是无核细胞，应该无法通过 scATAC-seq检测到。被预测为血小板的这些细胞有可能实际上是血小板前体——巨核细胞，巨核细胞主要位于骨髓中，但在健康患者的外周血中极少发现，例如这些PBMC来自的个体。考虑到正常骨髓中巨核细胞本身就非常稀有（< 0.1%），这种情况似乎不太可能。
如果要确定这群血小板细胞是什么就需要进一步去查看血小板细胞，同时如果数量较少的话可以选择去除

11.去除不恰当的细胞

predicted_id_counts <- table(pbmc$predicted.id)

# 识别出预测的id值，这些值对应的细胞数量超过20个
major_predicted_ids <- names(predicted_id_counts[predicted_id_counts > 20])
pbmc <- pbmc[, pbmc$predicted.id %in% major_predicted_ids]

# 将细胞身份更改为每个细胞的预测标签
Idents(pbmc) <- pbmc$predicted.id

12.查找细胞类型之间差异可接近的峰值

为了找到细胞簇之间的差异可接近区域，可以进行差异可接近性（DA）测试。一种简单的方法是进行Wilcoxon秩和检验，而presto包中实现了一种极快的Wilcoxon检验，可以在Seurat对象上运行。
改回使用峰值而不是基因活性。

DefaultAssay(pbmc) <- 'peaks'

# wilcox是test.use的默认选项。
da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14+ Monocytes",
  test.use = 'wilcox',
  min.pct = 0.1
)

head(da_peaks)

plot1 <- VlnPlot(
  object = pbmc,
  features = rownames(da_peaks)[1],
  pt.size = 0.1,
  idents = c("CD4 Naive","CD14+ Monocytes")
)
plot2 <- FeaturePlot(
  object = pbmc,
  features = rownames(da_peaks)[1],
  pt.size = 0.1
)

plot1 | plot2

13.使用ClosestFeature() 函数找到每个峰最近基因

open_cd4naive <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_cd14mono <- rownames(da_peaks[da_peaks$avg_log2FC < -3, ])

closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)
closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono)

14. 绘制基因组区域

# 为了绘图目的，将相关的细胞类型放在一起会很方便。我们可以通过运行 SortIdents() 函数，根据注释细胞类型之间的相似性自动排序绘图顺序
pbmc <- SortIdents(pbmc)

# 可视化 CD4 未成熟细胞和 CD14 单核细胞中开放的 DA 峰，这些峰位于这些细胞类型的一些关键标记基因附近，分别是 CD4 和 LYZ。在这里我们将用灰色突出显示 DA 峰区域。
# 找到与目标基因重叠的DA峰值
regions_highlight <- subsetByOverlaps(StringToGRanges(open_cd4naive), LookupGeneCoords(pbmc, "CD4"))

CoveragePlot(
  object = pbmc,
  region = "CD4",
  region.highlight = regions_highlight,
  extend.upstream = 1000,
  extend.downstream = 1000
)

regions_highlight <- subsetByOverlaps(StringToGRanges(open_cd14mono), LookupGeneCoords(pbmc, "LYZ"))

CoveragePlot(
  object = pbmc,
  region = "LYZ",
  region.highlight = regions_highlight,
  extend.upstream = 1000,
  extend.downstream = 5000
)

参考资料：

Signac:https:///signac/articles/pbmc_vignette

微精选