scATAC-seq 数据分析: Signac（官方流程复现 + 步骤详解）

一、安装并读取所需的R包

install.packages("Signac")
install.packages("BiocManager")
if (!requireNamespace("EnsDb.Hsapiens.v75", quietly = TRUE))
    BiocManager::install("EnsDb.Hsapiens.v75") #如果下载速度过慢导致失败直接去BioManger官网下载包再安装
install.packages("Seurat")
install.packages("spatstat.utils")
install.packages("spatstat.geom", dependencies = TRUE)
install.packages("hdf5r")
BiocManager::install("biovizBase")


library(EnsDb.Hsapiens.v75) #这个包提供了人类基因组版本Ensembl v75的注释信息
library(patchwork)
library(Signac)
library(Seurat)
library(ggplot2)
library(hdf5r)
library(spatstat.geom)
library(biovizBase)

二、去官网获取分析所需的数据

        atac_v1_pbmc_10k_singlecell.csv

        atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5

        atac_v1_pbmc_10k_fragments.tsv.gz

        atac_v1_pbmc_10k_fragments.tsv.gz.tbi

        官网下载地址:   atac_v1_pbmc_10k -Datasets -Single Cell ATAC -Official 10x Genomics Supporthttps://support.10xgenomics.com/single-cell-atac/datasets/1.0.1/atac_v1_pbmc_10k?

        根据 Signac 开发人员的建议，包含片段文件的文件夹也应包含片段索引文件 （.tbi）。例如，对于人类 10X PBMC 小插图数据 atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz 和 atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi 应位于同一文件夹中。

三、数据预处理，构建Siganc对象

        Siganc所需输入文件介绍

          ① Peak/Cell matrix ：

                行名为Peak区域（指染色质上的高信号区域，对应于基因组中的开放区域），列名为Barcode（细胞），矩阵内的数字代表在此Peak中Tn5酶的结合数。

             ② Fragment file ：

                指的是染色质测序（ChIP-seq、ATAC-seq 等）中记录测序片段位置和信息的文件。这些片段文件是染色质测序数据的核心部分，用于描述基因组上的测序片段的位置和性质。每行代表一个测序片段，包含染色体名称、片段的起始位置和结束位置等信息。

         ③ csv描述数据属性文件 ：

            这是一个包含染色质测序数据的元数据表格，其中记录了每个细胞的各种统计信息。这些信息来自于染色质测序数据的处理和分析，用于描述每个细胞的质量和属性。

• barcode：细胞条码，用于标识每个细胞
• total：总的片段数目
• duplicate：重复的片段数目
• chimeric：融合的片段数目
• unmapped：未映射的片段数目
• lowmapq：映射质量低的片段数目
• mitochondrial：位于线粒体的片段数目
• passed_filters：通过筛选的片段数目
• cell_id：细胞的唯一标识符
• is__cell_barcode：是否是细胞条码
• TSS_fragments：与转录起始位点相关的片段数目
• DNase_sensitive_region_fragments：与DNase敏感区域相关的片段数目
• enhancer_region_fragments：与增强子区域相关的片段数目
• promoter_region_fragments：与启动子区域相关的片段数目
• on_target_fragments：与目标区域相关的片段数目
• blacklist_region_fragments：位于黑名单区域的片段数目
• peak_region_fragments：位于峰值区域的片段数目

预处理代码：

# 读取Peak/Cell 矩阵
counts <- Read10X_h5(filename = "../vignette_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")


# 读取描述csv文件到metadate
metadata <- read.csv(
  file = "../vignette_data/atac_v1_pbmc_10k_singlecell.csv",
  header = TRUE,
  row.names = 1
)


# 使用Seurat 的 CreateChromatinAssay 函数创建一个染色质测序（Chromatin Assay）对象
chrom_assay <- CreateChromatinAssay(
# counts 参数指定计数矩阵
  counts = counts,
  sep = c(":", "-"),
# genome 指定基因组版本，常见的基因组版本包括人类基因组的 hg19（GRCh37） 和 hg38（GRCh38）版本
  genome = 'hg19',
  fragments = 'atac_v1_pbmc_10k_fragments.tsv.gz',
# 指定了最少的细胞数目。片段必须在至少这么多的细胞中被检测到，才会被包含在创建的染色质测序对象中
  min.cells = 10,
# 指定了最少的特征（片段）数目。片段必须在至少这么多的细胞中被检测到，才会被包含在创建的染色质测序对象中。
  min.features = 200
)


# 使用 Seurat 的 CreateSeuratObject 函数创建一个 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)


# pbmc[['peaks']] 打印了关于染色质测序数据的基本信息
pbmc[['peaks']]
# ChromatinAssay data with 87561 features for 8728 cells
# Variable features: 0 
# Genome: hg19 
# Annotation present: FALSE 
# Motifs present: FALSE 
# Fragment files: 1 

#granges(pbmc) 是一个 GRanges 对象，这个 GRanges 对象存储了染色质测序数据的区域信息，即片段的位置和相关属性
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对象添加人类基因组的基因注释

        下列代码的目的是将基因注释信息添加到 Seurat 对象中，以便在后续的分析中可以将染色质测序数据与基因注释进行比较和重叠分析

# 从 EnsDb 数据库中提取基因注释信息
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)

# 将注释数据的染色体编号改为 UCSC 格式（添加 'chr' 前缀）
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))

# 设置基因组版本为 "hg19"，以符合数据映射到 hg19 的情况
genome(annotations) <- "hg19"

# 将基因注释信息添加到 Seurat 对象中
Annotation(pbmc) <- annotations

        运行后 Seurat 对象 pbmc 具有了与基因注释信息相关的属性，可以在后续的分析中使用这些注释信息来深入探索染色质测序数据与基因之间的关系。

五、QC（质量控制）

        推荐使用以下指标来评估数据质量：

1. Nucleosome banding pattern（核小体条带分布模式）: 核小体条带分布模式是指DNA片段大小的直方图（通过配对末端测序读数确定）呈现出明显的核小体带模式，对应于包裹在单个核小体周围的DNA长度。对每个单细胞计算这个分布模式，并量化单核小体与无核小体片段的大致比率（存储为 nucleosome_signal）。核小体信号可以显示核小体在细胞中的分布情况（分布模式）。一个好的scATAC-seq实验通常会呈现出明显的核小体带模式，其中峰值对应于DNA包裹在单个核小体周围的长度。如果核小体信号没有明显的带状模式，可能表明数据存在噪声、污染或其他问题。

2. Transcriptional start site (TSS) enrichment score（转录起始位点富集分数）：ENCODE项目基于片段在TSS中心和TSS相邻区域之间的比率定义了一个ATAC-seq定位分数。这个分数反映了细胞中TSS附近的ATAC-seq片段相对于TSS相邻区域的富集程度。TSS是基因的转录起始位点，通常会有更多的开放染色质和ATAC-seq片段。因此，高富集分数意味着细胞中的ATAC-seq片段更多地集中在基因的TSS区域，这可能表示这些细胞在转录调控和基因表达方面具有更高的活性和功能。反之，低TSS富集分数可能意味着数据质量较差，或者细胞处于不活跃的状态，染色质在TSS附近较为封闭，ATAC-seq片段分布较均匀，可能不具有明确的TSS信号。通过计算TSS富集分数，我们可以从ATAC-seq数据中获得有关细胞的转录调控状态和基因表达活性的信息，进而用于细胞的质量控制和生物学解释。

3. Total number of fragments in peaks（峰值区域中的片段总数）：这是衡量细胞测序深度和复杂性的指标。读数非常少的细胞可能需要被排除在外，因为其测序深度较低。而具有极多片段数的细胞可能代表着双体细胞、细胞核团聚或其他伪迹。测序深度是指从一个细胞中获得的测序读数的总数，它反映了对细胞基因组的覆盖程度。测序深度较低的细胞可能会导致数据质量不足，很难进行准确的分析和解释。因此，具有非常少读数的细胞可能需要在后续分析中排除，以确保分析的可靠性。双体细胞是指在一个测序读数中捕获了两个细胞的信号，因此会有异常高的片段数。类似地，细胞核团聚或其他技术伪迹也可能导致片段数异常升高。这些细胞可能会影响数据的解释和分析，因此需要进行进一步的检查或排除。

4. Fraction of fragments in peaks（片段在峰值区域的比例）：表示所有片段中落在ATAC-seq峰值区域内的比例。具有低值的细胞（即<15-20%）通常代表低质量细胞或应该被排除的技术伪迹。需要注意的是，这个值对使用的峰值集合非常敏感。这个指标反映了测序数据中的片段有多少落在ATAC-seq的峰值区域内，也就是开放染色质区域。峰值区域通常与基因的调控元件和功能相关，因此高质量的细胞通常会有较高的片段比例落在这些区域内。低片段比例可能暗示细胞质量较差，或者可能存在一些技术伪迹。例如，数据处理过程中可能出现片段损失，或者细胞中的ATAC-seq片段分布不均匀。这可能会影响后续分析的结果可靠性。因此，在分析过程中通常会选择性去除低片段比例的细胞，以确保结果的可靠性和一致性。

5. Ratio reads in genomic blacklist regions（基因组黑名单区域中读数的比例）：ENCODE项目提供了一个黑名单区域列表，其中包含与伪迹信号经常相关的读数。具有大量读数映射到这些区域（相对于映射到峰值区域的读数）的细胞通常代表技术伪迹，应该被移除。Signac软件包中包含了人类（hg19和GRCh38）、小鼠（mm10）、果蝇（dm3）和秀丽隐杆线虫（ce10）的ENCODE黑名单区域。

        以下是代码示例，用于计算细胞的核小体信号分数、TSS富集分数，并添加黑名单比例和片段在峰值区域的比例（便于后续筛除低质量细胞）。然后使用函数DensityScatter()对对象元数据中存储的变量之间的关系进行可视化。

# 计算细胞的核小体信号（NS）分数
pbmc <- NucleosomeSignal(object = pbmc)

# 计算细胞的TSS富集分数
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

# 添加黑名单比例和片段在峰值区域的比例
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments

# 可视化变量之间的关系
DensityScatter(pbmc, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

元数据中存储的变量间关系的可视化图

     ①  转录起始位点富集分数作图

# 这段代码通过条件判断，将TSS富集分数大于3的细胞标记为'High'组，其他细胞标记为'Low'组
# 并调用TSSPlot()函数，绘制不同组细胞在TSS位点上的可及性信号图

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

      ②  核小体条带分布模式作图

# 计算每个细胞的核小体信号强度（NS），即计算位于 147bp 到 294bp 之间（单核小体）的片段与小于 147bp（无核小体）的片段的比率
# 将核小体信号分数大于4的细胞标记为'NS > 4'组，其他细胞标记为'NS < 4'组
# 核小体信号分数是之前计算的，用于反映细胞的核小体条带分布情况
# 官方根据多次实验，认为 NS > 4 的细胞质量不佳

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

      ③ 对第五步开头提到用于评估质量的5个指标进行可视化

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 > 3000 &             # 确保每个细胞的峰值片段数在3000以上
    nCount_peaks < 30000 &                   # 确保每个细胞的峰值片段数在30000以下
    pct_reads_in_peaks > 15 &                # 确保在峰值内的片段占比大于15%
    blacklist_ratio < 0.05 &                 # 确保位于黑名单区域的片段占比小于5%
    nucleosome_signal < 4 &                  # 确保核小体信号得分小于4
    TSS.enrichment > 3                       # 确保转录起始位点富集得分大于3
)

六、标准化和降维

     ① 标准化方式

        标准化是scATAC-seq数据分析中的重要步骤，旨在确保数据的可比性和适用性，以便进行后续分析。在Signac中，采用了一种特定的标准化方法，称为词频-逆文档频率（TF-IDF）标准化。

Signac中的TF-IDF标准化包括两个步骤：

• 细胞间标准化：此步骤旨在校正细胞测序深度的差异。确保具有较高测序深度的细胞不会在下游分析中占主导地位。每个细胞的ATAC-seq峰值会根据其整体测序深度进行标准化

• 峰值间标准化：此步骤涉及对数据在峰值间进行标准化。更罕见的峰值（即在较少细胞中出现的峰值）会获得较高的值。这有助于突显特定细胞群或条件下独特的峰值

        总之，TF-IDF标准化允许消除因测序深度和峰值频率差异而产生的偏差。这确保数据得到适当缩放，可以用于后续分析，如降维、聚类和差异可及性分析

     ② 特征选取方式

        与单细胞RNA测序（scRNA-seq）数据相比，单细胞ATAC测序（scATAC-seq）数据的动态范围较低，这使得进行变量特征选择变得更加困难。通常情况下，会对scRNA-seq数据进行变量特征选择，以便在降维和聚类等分析中仅使用重要的基因。然而，在处理scATAC-seq数据时，可以采取不同的策略。

        其中一种策略是选择仅使用前n％的特征（峰值）进行降维分析，或者使用FindTopFeatures()函数移除在不到n个细胞中出现的特征。这是因为scATAC-seq数据的动态范围较小，很难确定哪些特征具有明显的差异和变化。因此，通过选择一部分最显著的特征，可以减少降维的计算负担。

        在实际操作中，有时会使用所有的特征，但也可以只使用部分特征，例如仅使用前25％的峰值。这样做可以加快运行速度，而且在使用降维结果进行后续分析时，效果可能会相似。被用于降维的特征会自动由函数设置为Seurat对象的VariableFeatures()，以便在后续分析中使用。

        在本流程中我们选取所有峰值作为降维和分析的变量特征，不进行特征筛选。这可能会在某些情况下保留数据中所有的信息，但在其他情况下可能会增加计算的复杂性。

    ③ 降维方式

        降维的目的是将高维的数据表示转换为低维的表示，以便在保留主要信息的同时减少数据的复杂性。这一步中使用了一种叫做奇异值分解（SVD）的数学方法，应用在TF-IDF矩阵上，TF-IDF矩阵是由①中进行的TF-IDF标准化所得。

        SVD会将原始的TF-IDF矩阵分解成三个矩阵的乘积：U、S和V^T。其中，U矩阵代表数据点在降维空间中的表示，S是一个对角矩阵，包含了每个维度上的奇异值，V^T矩阵代表特征（峰值）在降维空间中的表示。

        这个过程类似于PCA（主成分分析），它在scRNA-seq领域经常使用。通过SVD，我们可以将高维的TF-IDF矩阵投影到低维的空间，从而得到一个新的、更紧凑的数据表示。这对于后续的可视化、聚类和分析等任务非常有帮助。

    ④ 代码实现

# 对pbmc对象运行TF-IDF标准化
pbmc <- RunTFIDF(pbmc)

# 在归一化后的数据上选择顶部特征（峰值），此处选取了所有特征
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')

# 在选择的特征上运行奇异值分解（SVD）进行降维
pbmc <- RunSVD(pbmc)

        在上述3条代码中均用到了LSI算法，LSI代表Latent Semantic Indexing（潜在语义索引），是一种用于文本挖掘和信息检索的技术。它的主要目标是通过降低维度和发现隐藏在数据中的语义关系来改善文本数据的表示和分析。在scATAC-seq数据中，LSI被用来进行降维，类似于主成分分析（PCA），以便将高维度的数据转换为较低维度的表示，同时保留数据的主要变异模式。

        然而，在某些情况下，第1个LSI组件可能会捕获与测序深度相关的技术变异，而不是真实的生物学变异。因此，对于后续的分析，通常需要评估LSI组件与测序深度之间的相关性，以确定是否应该将这个组件从分析中排除。

        我们可以使用DepthCor()函数来评估每个LSI组件与测序深度之间的相关性

DepthCor(pbmc)

       通过上图分析发现，第1个LSI组件与细胞的总计数（测序深度）之间存在强烈的相关性，而其他组件的分布均在可接受范围内。这意味着第一个LSI组件可能主要捕获了技术性的变异，而不是生物学上的变异。所以在后续的分析中，会忽略第1个LSI组件，以避免将细胞基于测序深度而分组在一起。这将有助于更好地捕获细胞的生物学变异和特征。

七、非线性降维和聚类

        之前使用的SVD降维（线性降维）可能无法捕捉到数据中的非线性结构和复杂模式。这就是为什么在许多单细胞分析中，包括scATAC-seq分析，研究人员会将SVD降维与非线性降维方法（如UMAP）结合使用的原因之一。UMAP是一种非线性降维方法，能够更好地捕捉细胞之间的非线性关系和相似性。它在保留数据的局部和全局结构方面表现出色，通常能够提供更有意义的降维可视化。

# 运行UMAP算法进行非线性降维，reduction 参数指定了要降维的数据类型（这里是'lsi'，为SVD降维后的低维表示）
pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)

# 寻找细胞之间的邻居关系
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)

# 使用图聚类算法进行细胞聚类，algorithm参数指定了用于执行聚类的算法。algorithm = 3表示使用Leiden聚类算法
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)

# 生成二维降维可视化图
DimPlot(object = pbmc, label = TRUE) + NoLegend()

八、创建基因表达活性矩阵

        scRNA-seq 分析中，通常可以通过对基因的表达水平来了解细胞类型和状态。而在 scATAC-seq 中，我们无法直接测量基因的表达水平，因为该技术主要测量了染色质区域的可及性，而不是基因的RNA表达。然而，我们仍然可以尝试通过染色质的可及性来估计基因的活性，从而为每个细胞计算一个基因活性分数。

       首先通过提取基因坐标，将这些坐标扩展到包括上游 2 kb 的区域（扩展基因坐标，特别是在基因的上游区域，可以更全面地捕获基因的启动子区域。这是因为在启动子区域上游可能存在一些调控元件，如增强子，它们对基因的调控也很重要。通过将基因坐标向上游扩展，我们可以涵盖这些可能对基因活性产生影响的区域，从而更准确地估计基因的活性水平）。然后使用FeatureMatrix()函数计算映射到每个基因区域的片段数。这个过程可以得到一个基因活性矩阵，矩阵的行代表细胞，列代表基因，每个单元格的值表示对应基因区域的片段数。

# 利用scATAC-seq数据，计算基因活性局矩阵
gene.activities <- GeneActivity(pbmc)

# 将基因活性矩阵作为一个新的测序数据（assay）添加到Seurat对象中
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)

# 对基因活性数据进行对数归一化处理
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',    # 使用'LogNormalize'方法来进行标准化
  scale.factor = median(pbmc$nCount_RNA)    # scale.factor 参数用于控制归一化的缩放因子，通常会选择数据集的中位数
)

# 查看基因活性矩阵
pbmc[['RNA']]

        经过以上步骤，可以对一些典型的 maker 基因进行可视化，以帮助解释 scATAC-seq 的聚类结果。需要注意的是，与 scRNA-seq 相比，该基因活性会产生更多的噪声。这是因为它们由开放染色质数据测量得到，同时也假设基因体/启动子的可及性与基因表达之间存在一般对应关系，而这种对应关系可能并不总是成立。尽管如此，我们还是可以基于此数据区分单核细胞、B细胞、T细胞和NK细胞等细胞群体。但仅通过监督性分析很难进一步细分这些细胞的类型。

        下段代码中，首先将 Seurat 对象默认的测序数据类型设置为 'RNA '，以便后续的分析。然后，使用 FeaturePlot 函数绘制了几个典型 marker 基因的活性，这些基因通常与特定细胞类型相关联。例如，'MS4A1' 是B细胞标记基因，'CD3D' 是 T 细胞标记基因，'LEF1' 是 T 和 B 细胞标记基因，'NKG7' 是 NK 细胞标记基因，'TREM1' 是单核细胞标记基因，'LYZ' 是单核细胞和单核吞噬细胞标记基因。

# 将默认测序数据类型设置为'RNA'
DefaultAssay(pbmc) <- 'RNA'

# 使用FeaturePlot函数可视化经典标记基因的活性
# 设置max.cutoff为'q95'，表示在可视化时仅显示95%分位点内的数据点
FeaturePlot(
  object = pbmc,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)

九、整合scRNA-seq数据进行分析

        获取整合所需的 rds 文件 ： https://signac-objects.s3.amazonaws.com/pbmc_10k_v3.rdshttps://signac-objects.s3.amazonaws.com/pbmc_10k_v3.rds

        很难直接从 scATAC-seq 数据中分辨细胞类型。为了帮助解释这些数据，可以借助与其同一生物系统（人类PBMC）相关的 scRNA-seq 实验信息来进行分类，从而更好区分细胞类型。 

        为了将 scRNA-seq 数据中的细胞分类应用到 scATAC-seq 数据中，我们使用了交叉模态整合和标签转移的方法。这意味着我们要寻找两种不同类型的数据之间的相关模式，以便在一个数据集中找到细胞类型的标签，然后将这些标签传递到另一个数据集中。

        为了执行这一任务，我们加载了一个预处理过的人类PBMCs的scRNA-seq数据集。这个数据集来自10x Genomics，并且已经进行了预处理，使其适合进行交叉模态整合和标签转移。

# 加载预处理的人类PBMC scRNA-seq数据
pbmc_rna <- readRDS("pbmc_10k_v3.rds")

# 寻找scRNA-seq和scATAC-seq数据之间的转移锚点
transfer.anchors <- FindTransferAnchors(
  reference = pbmc_rna,           # 参考的scRNA-seq数据
  query = pbmc,                   # 待传递细胞类型标签的scATAC-seq数据
  reduction = 'cca'               # 降维方法选择CCA
)

# 使用转移锚点从scRNA-seq传递细胞类型标签到scATAC-seq数据
predicted.labels <- TransferData(
  anchorset = transfer.anchors,     # 转移锚点
  refdata = pbmc_rna$celltype,      # scRNA-seq中的细胞类型标签作为参考
  weight.reduction = pbmc[['lsi']], # 使用scATAC-seq的LSI进行权重降维
  dims = 2:30                       # 计算转移时使用的维度范围
)

# 将预测的细胞类型标签添加为元数据到scATAC-seq对象中
pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)

# 创建scRNA-seq数据的维度降低图
plot1 <- DimPlot(
  object = pbmc_rna,
  group.by = 'celltype',  # 按照细胞类型进行分组
  label = TRUE,  # 显示标签
  repel = TRUE  # 防止标签重叠
) + NoLegend() + ggtitle('scRNA-seq')  # 禁用图例，添加标题

# 创建scATAC-seq数据的维度降低图
plot2 <- DimPlot(
  object = pbmc,
  group.by = 'predicted.id',  # 按照预测的细胞群体标签进行分组
  label = TRUE,  # 显示标签
  repel = TRUE  # 防止标签重叠
) + NoLegend() + ggtitle('scATAC-seq')  # 禁用图例，添加标题

# 将两个维度降低的图叠加在一起进行比较
plot1 + plot2

         可以看出，基于 scRNA-seq 的聚类结果与使用 scATAC-seq 数据计算的 UMAP 聚类结果一致。请注意，在 scATAC-seq 数据集中没有预测到任何血小板细胞，这是符合预期的，因为血小板不具有细胞核，无法通过 scATAC-seq 检测到。

        下段代码的目的是将 scATAC-seq 的聚类结果与基于标签传递的预测结果相结合，以获得更准确的细胞类型标识。通过将每个聚类中的细胞按照预测标签重新命名，我们可以更好地将 scATAC-seq 数据中的细胞与 scRNA-seq 数据中的预测标签对应起来。这样做有助于更准确地识别细胞类型，从而为后续的功能注释和解释提供更有价值的信息。

# 用预测的标签重新分配每个聚类的名称
for (i in levels(pbmc)) {
  # 获取当前聚类中的细胞
  cells_to_reid <- WhichCells(pbmc, idents = i)
  # 找到当前聚类中预测标签数量最多的标签
  newid <- names(which.max(table(pbmc$predicted.id[cells_to_reid])))
  # 将当前聚类中的细胞重新赋予新的预测标签
  Idents(pbmc, cells = cells_to_reid) <- newid
}

十、寻找不同细胞类型之间的差异可及性峰值（Peaks）

        寻找不同细胞类型之间的差异可及性峰值的意义在于，它可以帮助我们识别在不同细胞类型之间存在的基因座的可及性差异。这些差异可能反映了不同细胞类型之间的功能差异，如基因的调控状态、转录因子结合位点的变化等。通过分析这些差异可及性峰值，我们可以更好地理解不同细胞类型之间的生物学差异和功能特征，从而揭示基因座在细胞类型特定性和功能中的作用。

     ① 寻找差异可及性峰值       

         要在细胞聚类之间找到差异可及性区域，可以执行差异分析（Differential Accessibility，DA）测试。使用逻辑回归进行DA测试，这是根据Ntranos等人在2018年提出的方法，适用于 scRNA-seq 数据，并且将总的片段数作为潜在变量加入，以减轻差异测序深度对结果的影响。在这里，我们将重点比较 Naive CD4 细胞和 CD14 单核细胞之间的差异，使用这些方法其实可以比较任何细胞组。我们还可以将这些标记峰值可视化在小提琴图、特征图、点图、热图或 Seurat 中的任何可视化工具上。

# 切换回使用峰值数据而不是基因活性数据
DefaultAssay(pbmc) <- 'peaks'

# 寻找在"CD4 Naive"和"CD14+ Monocytes"细胞类型之间的差异可及性峰值
# 运行时间较长，可能需要 1 h
da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14+ Monocytes",
  test.use = 'LR',  # 使用逻辑回归作为差异分析方法
  latent.vars = 'nCount_peaks'  # 引入总的峰值数量作为潜在变量
)

# 显示差异可及性峰值的前几行
head(da_peaks)

##                                  p_val avg_log2FC pct.1 pct.2     p_val_adj
## chr14-99721608-99741934  3.914345e-277   5.644662 0.868 0.020 3.427440e-272
## chr14-99695477-99720910  1.987922e-218   5.100138 0.797 0.021 1.740645e-213
## chr17-80084198-80086094  1.033863e-216   7.093091 0.668 0.004 9.052607e-212
## chr7-142501666-142511108 2.064224e-215   4.954426 0.754 0.026 1.807455e-210
## chr2-113581628-113594911 2.932572e-184  -5.055830 0.035 0.666 2.567789e-179
## chr6-44025105-44028184   9.481165e-179  -4.414731 0.046 0.622 8.301803e-174

        上方输出是通过差异可及性分析找到的前几个差异峰值的结果，这些峰值在"CD4 Naive"和"CD14+ Monocytes"细胞类型之间有显著差异。每行结果表示一个峰值，以下是各列的解释：

• chr14-99721608-99741934：峰值的名称，通常由染色体名和坐标构成
• p_val：差异分析的P值，表示在这两个细胞类型之间差异可及性的显著性
• avg_log2FC：平均对数折叠变化，表示两个细胞类型之间的平均差异
• pct.1：在"CD4 Naive"细胞类型中，该峰值在所有峰值中的百分比
• pct.2：在"CD14+ Monocytes"细胞类型中，该峰值在所有峰值中的百分比
• p_val_adj：经过FDR校正的P值，用于控制多重比较的误差

        这些结果提供了在这两种细胞类型之间存在显著差异的峰值的信息，包括差异的显著性水平、平均差异大小以及这些峰值在不同细胞类型中的相对丰度。这些差异峰值可能涉及到不同细胞类型之间的特定功能差异。

        下方代码进行可视化：

# 使用 VlnPlot 函数绘制小提琴图，显示不同细胞群体中特定峰值的信号强度分布
plot1 <- VlnPlot(
  object = pbmc,
  features = rownames(da_peaks)[1],      # 特定的峰值名称
  pt.size = 0.1,
  idents = c("CD4 Naive","CD14+ Monocytes")  # 要绘制的细胞群体标签
)

# 使用 FeaturePlot 函数绘制特征在不同细胞上的分布情况（类似散点图）
plot2 <- FeaturePlot(
  object = pbmc,
  features = rownames(da_peaks)[1],      # 特定的峰值名称
  pt.size = 0.1
)

# 将两个图垂直排列，以便进行比较和分析
plot1 | plot2

• 第一张小提琴图显示了特定峰值的信号强度分布情况。图中的每个小提琴表示一个细胞群体（"CD4 Naive" 和 "CD14+ Monocytes"），并展示了该细胞群体中每个单个细胞的峰值信号强度分布。
• 第二张特征图（FeaturePlot）显示了特定峰值在每个细胞上的信号情况，类似于散点图。图中的每个点代表一个细胞，其水平位置对应于该细胞在特定峰值上的信号强度。通过这张图，你可以观察到不同细胞之间特定峰值的变化程度，以及细胞群体之间的差异。

        另一种在两组细胞之间寻找差异可及性区域的方法是查看两组细胞之间的可及性折叠变化（Fold Change）。与运行更复杂的差异可及性分析方法相比，这种方速度更快，但不能考虑细胞之间的潜在变量，例如总测序深度的差异，并且不执行任何统计检验。但这仍然是一种快速探索数据的有用方法，可以使用 Seurat 中的 FoldChange()函数进行实现。

# 速度较快
fc <- FoldChange(pbmc, ident.1 = "CD4 Naive", ident.2 = "CD14+ Monocytes")
# 按折叠变化排序
fc <- fc[order(fc$avg_log2FC, decreasing = TRUE), ]

head(fc)

##                          avg_log2FC pct.1 pct.2
## chr6-28416849-28417227     11.45411 0.067     0
## chr7-110665002-110665493   11.41080 0.073     0
## chr1-65010850-65011196     11.25947 0.054     0
## chr8-19317420-19317942     11.22146 0.061     0
## chr1-172836553-172836955   11.20312 0.058     0
## chr2-191380525-191380926   11.15811 0.056     0

        当然，也可以使用以下代码寻找  一群 / 每群  细胞的差异可及性峰值

# 示范用，不用运行
# 这段代码使用差异可及性分析来寻找在"CD4 Naive"细胞类型中表达显著高于其他细胞类型的峰值
da_peaks_single <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  test.use = 'LR',  
  latent.vars = 'nCount_peaks'  
)

# 示范用，不用运行
# 这段代码在整个 scATAC-seq 数据集中，进行差异可及性分析
da_peaks_multiple <- FindMarkers(
  object = pbmc,
  test.use = 'LR',  
  latent.vars = 'nCount_peaks'  
)

     ② 根据峰值寻找关联基因

        峰值的坐标本身可能很难单独解释。我们可以使用 ClosestFeature()函数找到每个峰值最近的基因。如果探索基因列表，你会发现在 Naive T 细胞中打开的峰值靠近 BCL11B 和 GATA3 等基因（T细胞分化的关键调节因子），而在单核细胞中打开的峰值靠近 CEBPB 等基因（单核细胞分化的关键调节因子）。我们可以通过对 ClosestFeature()返回的基因集进行基因本体富集分析来进一步研究这些结果，有许多R包可以实现这一功能（例如，GOstats包）。

# 从差异可及性峰值分析结果中选择平均折叠变化大于3的峰值，这代表在“CD4 Naive”细胞群体中的可及性增加的峰值区域
open_cd4naive <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])

# 从差异可及性峰值分析结果中选择平均折叠变化小于-3的峰值，这代表在“CD14+ Monocytes”细胞群体中的可及性增加的峰值区域
open_cd14mono <- rownames(da_peaks[da_peaks$avg_log2FC < -3, ])

# 使用ClosestFeature函数，找到与在“CD4 Naive”细胞群体中可及性增加的峰值最接近的基因
closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)

# 使用ClosestFeature函数，找到与在“CD14+ Monocytes”细胞群体中可及性增加的峰值最接近的基因
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
## ENST00000546420 ENST00000546420    CCDC64 ENSG00000135127 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
## ENST00000546420 chr12-120427684-120428101 chr12-120426014-120428613        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
## ENST00000409245 ENST00000409245         TTC7A ENSG00000068724 protein_coding
## ENST00000484822 ENST00000484822          RXRA ENSG00000186350 protein_coding
##                 type           closest_region             query_region distance
## ENST00000432018  cds chr2-113593760-113593806 chr2-113581628-113594911        0
## ENSE00001638912 exon   chr6-44041650-44042535   chr6-44025105-44028184    13465
## ENST00000445003  gap  chr20-48884201-48894027  chr20-48889794-48893313        0
## ENST00000568649  cds  chr15-75335782-75335877  chr15-75334903-75336779        0
## ENST00000409245  cds   chr2-47300841-47301062   chr2-47297968-47301173        0
## ENST00000484822  gap chr9-137211331-137293477 chr9-137263243-137268534        0
Plotting genomic regions

• tx_id: 转录本的唯一标识符
• gene_name: 基因名称
• gene_id: 基因的唯一标识符
• gene_biotype: 基因的生物学类型
• type: 关联的类型（如CDS、exon、gap等）
• closest_region: 最近的区域
• query_region: 查询的区域（差异峰值区域）
• distance: 基因与查询区域之间的距离

        利用 ClosestFeature() 进行峰值特征距离分析具有依赖距离计算、不考虑功能关系、基因注释可能不完整版、只考虑距离最近的基因等缺点。

        当分析差异可及性峰值时，除了使用ClosestFeature()函数找到最近的基因外，还有其他常用且更精准的策略可以帮助解释和分析这些峰值的生物学意义。一些常用策略包括：

• 基因本体富集分析（GO Enrichment Analysis）：对差异峰值附近的基因进行基因本体富集分析，以确定这些峰值可能参与的生物学过程、细胞组分和分子功能。常用的R包包括GOstats、clusterProfiler等。

• 通路富集分析（Pathway Enrichment Analysis）：类似于基因本体富集分析，但是关注的是差异峰值附近基因的通路富集，从而揭示这些峰值可能影响的信号通路和生物学功能。常用的R包包括ReactomePA、KEGGprofile等。

十一、基因组ATAC信号可视化

  CoveragePlot()函数的作用是在基因组的不同区域上绘制 Tn5 整合频率，通过将细胞按照聚类、细胞类型或其他储存在对象中的元数据进行分组。这些图像呈现了伪批量的可及性轨迹，其中来自同一组细胞的信号已经被平均在一起，用来可视化特定区域内的DNA可及性。除了DNA可及性轨迹外，如果数据存在，CoveragePlot()函数还可以可视化基因注释、峰值坐标和基因组链接等重要信息。

        延伸基因组区域的作用是为了更全面地展示所选区域的Tn5整合频率分布。通过将基因组区域向上游和下游延伸，我们可以捕捉到更多可能与该区域相关的信号，进而更好地了解这个区域的DNA可及性特征

# 设置绘图顺序，重新排列细胞类型
levels(pbmc) <- c("CD4 Naive", "CD4 Memory", "CD8 Naive", "CD8 effector", "Double negative T cell", "NK dim", "pre-B cell", "B cell progenitor", "pDC", "Dendritic cell", "CD14+ Monocytes", "CD16+ Monocytes")

# 使用CoveragePlot()函数绘制特定基因组区域的可及性轨迹
CoveragePlot(
  object = pbmc,                             # 要绘制的Seurat对象
  region = rownames(da_peaks)[1],            # 要绘制的基因组区域，选择差异峰值中的第一个峰值，此代码中为BCL11B基因，BCL11B在细胞分化和发育中扮演着重要的角色，尤其是在淋巴细胞分化和免疫调控中
  extend.upstream = 40000,                   # 向上游延伸的基因组区域长度
  extend.downstream = 20000                  # 向下游延伸的基因组区域长度
)

• X 轴坐标: 表示基因组的坐标位置，从上游向下游逐渐增加

• Y 轴坐标: 表示 ATAC-seq 信号强度，即在给定基因组位置的开放性水平

• Peaks 线条: 在图中，会出现表示 ATAC-seq 峰值的突出线条。这些线条的高度表示在该基因组位置的可及性信号强度

• Coverage 曲线: 连接峰值的高度形成的曲线表示在基因组区域内 ATAC-seq 信号的覆盖度分布情况

• 延伸区域: 由于设置了 extend.upstream 和 extend.downstream 参数，图中还包括了延伸的上游和下游基因组区域。这有助于观察峰值及其周围区域的可及性情况

        上图可以直观地看到不同细胞群体在特定基因组区域内的 ATAC-seq 信号强度和所选取基因（BCL11B）在染色体（14号染色体）上的位置，两者相互匹配，从而帮助了解不同细胞类型或群体之间的染色质开放性差异。这可以用来解释细胞类型之间的差异、调查调控元件的功能以及研究基因表达的调控机制。

参考文献：

Analyzing PBMC scATAC-seq • Signac (stuartlab.org)https://stuartlab.org/signac/articles/pbmc_vignette.html