教程链接：https://satijalab.org/signac/articles/pbmc_vignette.html
发布时间：March 08, 2022

signac的github链接：https://github.com/timoast/signac
signac的官网：https://satijalab.org/signac/index.html

在本教程中，我们将分析 10x Genomics 提供的人类外周血单个核细胞 (PBMC) 的单细胞 ATAC-seq 数据集。教程中用到的所有文件均可通过 10x Genomics 网站获得：

• Raw data (atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5)
• Metadata (atac_v1_pbmc_10k_singlecell.csv)
• fragments文件 (atac_v1_pbmc_10k_fragments.tsv.gz)
• fragments文件索引 (atac_v1_pbmc_10k_fragments.tsv.gz.tbi)

一、下载数据

在shell端运行下面脚本下载所有必须的文件：

wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi

二、加载分析所用的R包

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

library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Hsapiens.v75)
library(ggplot2)
library(patchwork)
set.seed(1234)

• Peak/Cell matrix. 该文件类似于用于分析单细胞RNA-seq的基因表达计数矩阵。然而，矩阵的每一行不是基因，而是代表一个基因组区域（“peak region”，一个峰），该区域被预测为一个开放的染色质区域。矩阵中的每个值代表每个峰内映射的每个barcode（即一个细胞）的 Tn5 整合位点的数量。你可以在 10X 网站上找到更多详细信息。

• Fragment file. 该文件储存了所有单个细胞中唯一片段（unique fragments）的完整列表。该文件非常大，处理速度较慢，并且存储在磁盘上（而不是内存中）。但是，保留此文件的优点是它包含与每个细胞相关的所有片段，而不是仅映射到峰的片段。有关片段文件的更多信息，请访问 10x Genomics 网站或 sinto 网站。

首先，我们使用由cellranger-atac生成的Peak/Cell矩阵和细胞的metadata创建一个Seurat对象，并将磁盘上的fragment文件的路径存储在Seurat对象中：

# file1-读取Peak/Cell矩阵
counts <- Read10X_h5(filename = "./vignette_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
# 查看peaks和细胞的数目
dim(counts)
# [1] 89796  8728
counts[1:5,1:5]
# 5 x 5 sparse Matrix of class "dgCMatrix"
# AAACGAAAGAGCGAAA-1 AAACGAAAGAGTTTGA-1 AAACGAAAGCGAGCTA-1
# chr1:565107-565550                  .                  .                  .
# chr1:569174-569639                  .                  .                  .
# chr1:713460-714823                  .                  .                  2
# chr1:752422-753038                  .                  .                  .
# chr1:762106-763359                  .                  .                  4
# AAACGAAAGGCTTCGC-1 AAACGAAAGTGCTGAG-1
# chr1:565107-565550                  .                  .
# chr1:569174-569639                  .                  .
# chr1:713460-714823                  8                  .
# chr1:752422-753038                  .                  .
# chr1:762106-763359                  2                  .

# file2-读取细胞注释信息
metadata <- read.csv(
  file = "./vignette_data/atac_v1_pbmc_10k_singlecell.csv",
  header = TRUE,
  row.names = 1
)
# 查看metadata信息
names(metadata)
# [1] "total"                            "duplicate"                       
# [3] "chimeric"                         "unmapped"                        
# [5] "lowmapq"                          "mitochondrial"                   
# [7] "passed_filters"                   "cell_id"                         
# [9] "is__cell_barcode"                 "TSS_fragments"                   
# [11] "DNase_sensitive_region_fragments" "enhancer_region_fragments"       
# [13] "promoter_region_fragments"        "on_target_fragments"             
# [15] "blacklist_region_fragments"       "peak_region_fragments"  
head(metadata)[,1:5]
# total duplicate chimeric unmapped lowmapq
# NO_BARCODE         8335327   3922818    71559   746476  634014
# AAACGAAAGAAACGCC-1       3         1        0        1       0
# AAACGAAAGAAAGCAG-1      14         1        0        4       2
# AAACGAAAGAAAGGGT-1       7         1        0        1       0
# AAACGAAAGAAATACC-1    9880      5380       79      106    1120
# AAACGAAAGAAATCTG-1       2         0        0        0       0

# 创建ChromatinAssay
chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  genome = 'hg19',
  fragments = './vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz',
  min.cells = 10,
  min.features = 200
)

# 构建seurat对象
pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

pbmc
# An object of class Seurat 
# 87561 features across 8728 samples within 1 assay 
# Active assay: peaks (87561 features, 0 variable features)

ATAC-seq数据使用自定义的assay— ChromatinAssay 进行存储。这使得一些专门的功能可以用于分析单细胞的基因组assay，例如 scATAC-seq。通过打印该assay，我们可以看到染色质检测中包含的一些附加信息，包括motif基序信息、基因注释和基因组信息。

pbmc[['peaks']]
# ChromatinAssay data with 87561 features for 8728 cells
# Variable features: 0 
# Genome: hg19 
# Annotation present: FALSE 
# Motifs present: FALSE 
# Fragment files: 1 

另外，我们可以调用 Seurat 对象上的 granges()函数，并将 ChromatinAssay 设置为活动分析（或在 ChromatinAssay 上），以查看与对象中每个特征相关的基因组范围。有关 ChromatinAssay 类的更多信息，请参阅对象交互教程。

granges(pbmc)
## GRanges object with 87561 ranges and 0 metadata columns:
##           seqnames            ranges strand
##              <Rle>         <IRanges>  <Rle>
##       [1]     chr1     565107-565550      *
##       [2]     chr1     569174-569639      *
##       [3]     chr1     713460-714823      *
##       [4]     chr1     752422-753038      *
##       [5]     chr1     762106-763359      *
##       ...      ...               ...    ...
##   [87557]     chrY 58993392-58993760      *
##   [87558]     chrY 58994571-58994823      *
##   [87559]     chrY 58996352-58997331      *
##   [87560]     chrY 59001782-59002175      *
##   [87561]     chrY 59017143-59017246      *
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

我们还可以为pbmc对象的人类基因组信息添加基因注释。这使得下游函数直接从对象中提取基因注释信息。

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
# change to UCSC style since the data was mapped to hg19
seqlevelsStyle(annotations) <- 'UCSC'
# add the gene information to the object
Annotation(pbmc) <- annotations

四、计算QC指标

下面，我们为前面的scATAC-seq对象计算一些 QC 指标。我们建议使用以下指标来评估数据质量。跟scRNA-seq一样，这些参数的预期值范围将根据样本的生物系统、细胞活力和其他因素而有所不同。

• 核小体条带模式 (Nucleosome banding pattern)：DNA 片段大小（DNA fragment sizes）的直方图（由配对末端测序read确定），根据包裹在单个核小体周围的DNA长度显示很强的核小体条带模式。我们计算每个单个细胞，并量化单核小体与无核小体片段的大致比率（存储为 nucleosome_signal）；
• 转录起始位点富集分数(Transcriptional start site (TSS) enrichment score) ：ENCODE 项目根据以 TSS 为中心的片段与 TSS 侧翼区域中的片段的比率定义了 ATAC-seq 靶向分数（参见 https://www.encodeproject.org/data-standards/terms/）。较差的 ATAC-seq 实验通常具有较低的 TSS 富集分数。我们可以使用 TSSEnrichment() 函数为每个细胞计算此指标，并将结果存储在列名 TSS.enrichment下的元数据中；
• peaks中的片段总数 (Total number of fragments in peaks)：衡量细胞测序深度/复杂性的指标。由于测序深度低，质控可能需要排除reads数很少的细胞。reads数极高的细胞可能代表了双细胞、细胞核团块；
• peaks中片段的比例 (Fraction of fragments in peaks)：表示落在ATAC-seq peaks内的总片段比例。比例较低的细胞（即<15-20％）通常表示在质控中应删除的低质量细胞或技术误差。请注意，此值可能对所使用的peaks值敏感。
• 比对到基因组"blacklist"区域中的reads比率 (Ratio reads in genomic blacklist regions)：ENCODE 项目提供了基因组上的"blacklist"区域列表，该区域通常与人工信号相关的读取。reads比对到这些区域的比例较高的细胞（与reads比对到peaks区域比例相比）通常代表了技术误差，应将其删除。 Signac包中包含针对人类（hg19 和 GRCh38）、小鼠 (mm10)、果蝇 (dm3) 和秀丽隐杆线虫 (ce10)物种 在ENCODE “blacklist”区域。

请注意，最后三个指标可以从CellRanger的输出中获得（存储在对象元数据中），但也可以使用Signac为非10x数据集计算该指标。

pbmc <- NucleosomeSignal(object = pbmc)

# compute TSS enrichment score per cell
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

# add blacklist ratio and fraction of reads in peaks
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments

我们可以根据分数对细胞进行分组，并绘制所有TSS位点的可及性信号，从而检查TSS富集分数。
在TSSEnrichment()函数中，设置 fast=TRUE选项将仅计算TSS富集分数，而不按每个细胞的 Tn5插入频率，存储每个细胞的位置矩阵，可以节省内存。但是，设置 fast=TRUE 将不允许使用 TSSPlot() 函数对不同细胞组的 TSS 富集信号进行下游绘图，如下所示：

pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()

image.png

我们还可以观察所有细胞的fragment长度规律性变化，并按核小体信号强度高低进行细胞分组。可以看到，单核小体/无核小体比率的异常值（基于上图）的细胞具有不同的核小体带型。剩下的细胞表现出ATAC-seq实验成功的典型模式。

pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')

image.png

用小提琴图呈现以上指标在每个细胞的分布情况；

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

image.png

然后，我们在质控中，去除QC指标异常的细胞。

# 根据不同QC指标进行数据过滤
pbmc <- subset(
  x = pbmc,
  subset = peak_region_fragments > 3000 &
    peak_region_fragments < 20000 &
    pct_reads_in_peaks > 15 &
    blacklist_ratio < 0.05 &
    nucleosome_signal < 4 &
    TSS.enrichment > 2
)
pbmc
## An object of class Seurat 
## 87561 features across 7060 samples within 1 assay 
## Active assay: peaks (87561 features, 0 variable features)

五、归一化和线性降维

• 归一化 (Normalization)：Signac是基于词频-逆文档频次 (term frequency-inverse document frequency，TF-IDF) 方法进行的归一化。它是将进行两步归一化操作，既跨细胞归一化以校正细胞测序深度的差异，又跨峰（peaks）归一化以为更罕见的峰提供更高的值。

TF-IDF是一种用于信息检索与文本挖掘的常用加权技术。TF-IDF是一种统计方法，用以评估一字词对于一个文件集或一个语料库中的其中一份文件的重要程度。字词的重要性随着它在文件中出现的次数成正比增加，但同时会随着它在语料库中出现的频率成反比下降。TF-IDF加权的各种形式常被搜索引擎应用，作为文件与用户查询之间相关程度的度量或评级。
TF-IDF的主要思想是：如果某个单词在一篇文章中出现的频率TF高，并且在其他文章中很少出现，则认为此词或者短语具有很好的类别区分能力，适合用来分类。

• 特征选择 (Feature selection)：scATAC-seq数据的低动态范围使得筛选高可变特征选择具有挑战性，就像我们对 scRNA-seq 所做的那样。相反，我们可以仅使用前 n% 的特征（峰值）进行降维，或者使用FindTopFeatures()函数删除少于n个细胞中的特征。在这里，我们将使用所有特征，尽管我们注意到，在仅使用特征子集时，我们也得到非常相似的结果（尝试将 min.cutoff 设置为“q75”以使用前 25% 的所有峰值），运行时间更快。此函数自动将用作降维的特征集，即Seurat对象的VariableFeatures()。

• 降维 (Dimension reduction)：接下来，我们使用上面选择的特征（峰值）在 TD-IDF 矩阵上运行奇异值分解（SVD）。这将返回该对象的低维数据表示结果（对于更熟悉 scRNA-seq的用户，可以将其视为类似于PCA的输出）。

# 使用RunTFIDF函数进行数据归一化
pbmc <- RunTFIDF(pbmc)
# 使用FindTopFeatures函数选择可变的features
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
# 使用RunSVD函数进行线性降维
pbmc <- RunSVD(pbmc)

线性降维后，第一个LSI成分通常捕获测序深度（技术差异），而不是生物学差异。在这种情况下，应从下游分析中删除该成分。我们可以使用DepthCor函数评估每个LSI成分与测序深度之间的相关性：

DepthCor(pbmc)

image.png

由上图，可以看到第一个LSI组分与细胞的总计数counts 之间存在非常强的相关性，因此我们将去除此组分执行下游步骤。

六、非线性降维和聚类

数据线性降维后，我们可以使用分析scRNA-seq数据的常规方法来对细胞进行基于图的聚类，并进行非线性降维（如UMAP等）可视化。其中，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()

image.png

七、创建基因表达活性矩阵（gene activity matrix）

从上图UMAP降维可视化的结果中可以看出，在人的外周血中存在多个细胞类群。如果我们还对PBMC的scRNA-seq分析结果熟悉的话，甚至还可以在scATAC-seq数据中发现某些髓样和淋巴样细胞群体的存在。但是，在scATAC-seq数据中对不同的细胞cluster进行注释和解释更具有挑战性，因为我们对非编码基因组区域功能的了解远少于对蛋白质编码区域（基因）的了解。

但是，我们可以尝试通过评估与每个基因相关的染色质可及性来量化基因组中每个基因的表达活性，并创建一个基于scATAC-seq数据的新基因活性测定方法。在这里，我们将使用一种简单的方法来对与基因体和启动子区域相交的reads数进行求和（我们也建议使用Cicero工具，它可以实现类似的目标）。

为了创建基因表达活性矩阵（gene activity matrix），我们从EnsembleDB中提取了人类基因组的基因坐标，并将其扩展到TSS上游2kb区域（因为启动子区域的可及性通常与基因表达相关）。然后，我们使用FeatureMatrix函数计算映射到这些区域的每个细胞的片段数。这将获取任何一组基因组坐标，计算与每个细胞中的这些坐标相交的reads次数，并返回一个稀疏矩阵。这些步骤由 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) 
)

现在，我们可以对一些典型的marker基因的“活性”进行可视化展示，以帮助我们解释scATAC-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 
)

image.png

八、和scRNA-seq数据整合

为了帮助解释scATAC-seq数据，我们可以基于来自同一生物系统（人PBMC）的scRNA-seq实验数据对细胞进行分类。我们利用此处描述的跨模态集成和标签转移方法，，将scATAC-seq数据和scRNA-seq数据进行整合。我们基于scATAC-seq数据的基因活性矩阵和scRNA-seq数据集的基因表达矩阵识别它们之间共享的相关模式，以鉴定两种模式下匹配的生物学状态。该过程根据每个scRNA-seq定义的簇标签返回每个细胞的分类评分。

在这里，我们利用跨模式整合和标签传输的方法，将scATAC-seq数据和scRNA-seq数据进行整合。我们基于scATAC-seq数据的基因活性矩阵和scRNA-seq数据集的基因表达矩阵识别它们之间共享的相关模式，以鉴定两种模式下匹配的生物学状态。该过程根据每个scRNA-seq定义的簇标签返回每个细胞的分类评分。

在这里，我们加载了人类 PBMC 的预处理 scRNA-seq 数据集，同样由 10x Genomics 提供。可以从 10x 网站下载此实验的原始数据，并在 GitHub 上查看用于构建此对象的代码。或者，您可以在此处下载预处理的Seurat对象。

# Load the pre-processed scRNA-seq data for PBMCs 
pbmc_rna <- readRDS("./vignette_data/pbmc_10k_v3.rds") 
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

image.png

你可以看到基于scRNA的分类与UMAP可视化完全一致。我们现在可以轻松地注释我们的 scATAC-seq 生成的cluster（或者我们可以使用RNA本身的分类）。
我们注意到cluster14 被标记成 CD4 记忆 T 细胞，但它是一个非常小的细胞群，具有较低的 QC 指标。由于该群可能代表低质量的细胞，我们将其从下游分析中删除。

pbmc <- subset(pbmc, idents = 14, invert = TRUE)
pbmc <- RenameIdents(
  object = pbmc,
  '0' = 'CD14 Mono',
  '1' = 'CD4 Memory',
  '2' = 'CD8 Effector',
  '3' = 'CD4 Naive',
  '4' = 'CD14 Mono',
  '5' = 'DN T',
  '6' = 'CD8 Naive',
  '7' = 'NK CD56Dim',
  '8' = 'pre-B',
  '9' = 'CD16 Mono',
  '10' = 'pro-B',
  '11' = 'DC',
  '12' = 'NK CD56bright',
  '13' = 'pDC'
)

九、查找cluster间差异可及性peaks区域

为了鉴定出不同细胞簇之间的差异可及性区域，我们可以进行差异可及性（differential accessibility，DA）检验分析，并使用不同的函数进行可视化展示。正如Ntranos等人所建议的，我们使用逻辑回归模型（LR）进行DA分析。并添加片段总数作为潜在变量，以减轻差异测序深度对结果的影响。对于稀疏数据（例如 scATAC-seq），我们发现通常需要将 FindMarkers() 中的min.pctthreshold设置比默认值设置的更低（0.1是为 scRNA-seq 数据设计的）。

下面，我们将重点比较Naive CD4细胞和CD14 monocytes细胞。我们还可以在小提琴图、特征图、点图、热图或 Seurat 中的任何可视化工具上可视化这些标记峰（ marker peaks）。

# change back to working with peaks instead of gene activities
DefaultAssay(pbmc) <- 'peaks'

da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14 Mono",
  min.pct = 0.05,
  test.use = 'LR',
  latent.vars = 'peak_region_fragments'
)

head(da_peaks)
##                                  p_val avg_log2FC pct.1 pct.2     p_val_adj
## chr14-99721608-99741934  1.622701e-307  0.7746012 0.879 0.023 1.420853e-302
## chr14-99695477-99720910  3.713144e-241  0.7611656 0.825 0.023 3.251266e-236
## chr17-80084198-80086094  1.789358e-231  0.9220662 0.695 0.006 1.566780e-226
## chr7-142501666-142511108 1.284659e-230  0.6751046 0.767 0.032 1.124860e-225
## chr2-113581628-113594911 2.820671e-212 -0.7572235 0.055 0.693 2.469808e-207
## chr6-44025105-44028184   3.989239e-206 -0.7447979 0.063 0.654 3.493017e-201

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

plot1 | plot2

image.png

fc <- FoldChange(pbmc, ident.1 = "CD4 Naive", ident.2 = "CD14 Mono")
head(fc)
##                     avg_log2FC pct.1 pct.2
## chr1-565107-565550  0.03917220 0.012 0.004
## chr1-569174-569639  0.09343254 0.032 0.008
## chr1-713460-714823  0.16852080 0.421 0.208
## chr1-752422-753038 -0.29489502 0.010 0.094
## chr1-762106-763359  0.11296934 0.290 0.162
## chr1-779589-780271  0.22080483 0.051 0.002

峰坐标很难单独解释。我们可以使用 ClosestFeature() 函数提供EnsDb中基因的注释信息，找到与每个peak最接近的基因。如果我们查看基因列表，会发现naive T细胞中的开放峰与BCL11B和GATA3（T细胞分化的关键调控因子）等基因接近，而单核细胞中开放的峰与CEBPB（单核细胞分化的关键调控因子）基因接近。我们可以通过对 ClosestFeature() 返回的基因集进行基因GO富集分析来进一步跟进这个结果，有许多R包可以做到这一点（例如，参见 GOstats 包）。

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

closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)
closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono)
head(closest_genes_cd4naive)
##                           tx_id   gene_name         gene_id   gene_biotype type
## ENST00000443726 ENST00000443726      BCL11B ENSG00000127152 protein_coding  cds
## ENST00000357195 ENST00000357195      BCL11B ENSG00000127152 protein_coding  cds
## ENST00000583593 ENST00000583593      CCDC57 ENSG00000176155 protein_coding  cds
## ENSE00002456092 ENST00000463701       PRSS1 ENSG00000204983 protein_coding exon
## ENST00000455990 ENST00000455990       HOOK1 ENSG00000134709 protein_coding  cds
## ENST00000557595 ENST00000557595 AE000662.92 ENSG00000259003 protein_coding  cds
##                           closest_region             query_region distance
## ENST00000443726  chr14-99737498-99737555  chr14-99721608-99741934        0
## ENST00000357195  chr14-99697682-99697894  chr14-99695477-99720910        0
## ENST00000583593  chr17-80085568-80085694  chr17-80084198-80086094        0
## ENSE00002456092 chr7-142460719-142460923 chr7-142501666-142511108    40742
## ENST00000455990   chr1-60280790-60280852   chr1-60279767-60281364        0
## ENST00000557595  chr14-23027629-23027664  chr14-23012859-23028219        0

head(closest_genes_cd14mono)
##                           tx_id     gene_name         gene_id   gene_biotype 
## ENST00000432018 ENST00000432018          IL1B ENSG00000125538 protein_coding 
## ENSE00001638912 ENST00000455005 RP5-1120P11.3 ENSG00000231881        lincRNA 
## ENST00000445003 ENST00000445003 RP11-290F20.3 ENSG00000224397        lincRNA 
## ENST00000568649 ENST00000568649         PPCDC ENSG00000138621 protein_coding 
## ENST00000484822 ENST00000484822          RXRA ENSG00000186350 protein_coding 
## ENST00000560869 ENST00000560869       TNFAIP2 ENSG00000185215 protein_coding 
##                 type            closest_region              query_region 
## ENST00000432018  cds  chr2-113593760-113593806  chr2-113581628-113594911 
## ENSE00001638912 exon    chr6-44041650-44042535    chr6-44025105-44028184 
## ENST00000445003  gap   chr20-48884201-48894027   chr20-48889794-48893313 
## ENST00000568649  cds   chr15-75335782-75335877   chr15-75334903-75336779 
## ENST00000484822  gap  chr9-137211331-137293477  chr9-137263243-137268534 
## ENST00000560869  utr chr14-103589798-103590288 chr14-103587259-103591113 
##                 distance 
## ENST00000432018        0 
## ENSE00001638912    13465 
## ENST00000445003        0 
## ENST00000568649        0 
## ENST00000484822        0 
## ENST00000560869        0

十、基因组区域绘图

我们可以使用 CoveragePlot() 函数绘制按cluster、细胞类型或存储在对象中的任何其他元数据分组的细胞绘制跨基因组区域的Tn5整合频率。这些代表pseudo-bulk可及性轨迹，将来自组内所有细胞的信号已被平均在一起以可视化区域中的DNA可访问性（感谢 Andrew Hill 在他出色的博客文章中为该功能提供灵感）。除了这些可及性轨迹外，我们还可以可视化其他重要信息，包括基因注释、峰坐标和基因组链接（如果它们存储在对象中）。

# set plotting order 
levels(pbmc) <- c("CD4 Naive","CD4 Memory","CD8 Naive","CD8 Effector","DN T","NK CD56bright","NK CD56Dim","pre-B",'pro-B',"pDC","DC","CD14 Mono",'CD16 Mono') 
CoveragePlot( 
  object = pbmc, 
  region = rownames(da_peaks)[1], 
  extend.upstream = 40000, 
  extend.downstream = 20000 
)

image.png

参考资料：

使用Signac包进行单细胞ATAC-seq数据分析（一）：Analyzing PBMC scATAC-seq