一直以来都不是太清楚scATAC-seq中fragments文件、peak calling，以及counts矩阵的关系，本文就来梳理一下scATAC-seq的上游数据处理中的一些基本概念。

---

fragment

• fragment是Tn5转座酶随机插入到开放染色质中剪切出的DNA片段，再经过测序和mapping到基因组上就得到了fragments文件。

Tn5转座酶剪切开放DNA示意图

• Tn5转座酶以同源二聚体的形式结合到DNA上，在两个Tn5分子间隔着9bp的DNA序列。根据这个情况，每个Tn5同源二聚体的结合事件会产生2个Insertions，中间隔着9bp。因此，真实的"开放"位置的中心在Tn5二聚体的正中间，而不是Tn5的插入位置。为了尽可能的还原真实情况，一般会对Tn5的Insertions进行了校正，即正链的插入结果往右移动4bp(+4 bp), 负链的插入结果往左偏移5bp(-5 bp)。

---

fragments文件

• Signac官方教程给出fragments文件的定义是：

fragments文件表示所有单个单元格中所有独特片段的完整列表。它是一个相当大的文件，处理速度较慢，并且存储在磁盘上（而不是内存中）。但是，保留此文件的优点是它包含与每个单个单元格关联的所有片段，而不是仅映射到峰值的片段。有关片段文件的更多信息可以在 10x Genomics网站或sinto 网站上找到。

• 以Signac官方给出的pbmc_10k的数据为例，来看看fragments文件和peak counts矩阵到底是什么关系。首先读入fragments文件

F <- read.table("/local/txm/txmdata/scATAC_seq/data/Signac/atac_v1_pbmc_10k_fragments.tsv.gz")
dim(F)
# 189803188         5   居然有将近两亿行！
head(F)
#     V1    V2    V3                 V4 V5
# 1 chr1 10066 10279 TTAGCTTAGGAGAACA-1  2
# 2 chr1 10072 10279 TTAGCTTAGGAGAACA-1  2
# 3 chr1 10079 10316 ATATTCCTCTTGTACT-1  2
# 4 chr1 10084 10340 CGTACAAGTTACCCAA-1  1
# 5 chr1 10085 10271 TGTGACAGTACAACGG-1  1
# 6 chr1 10085 10339 CATGCCTTCTCTGACC-1  1

• 可以看到fragments文件有5列，将近两亿行，第一列是染色体名称，第二列是单个fragment的起始位置，第三列是单个fragment的终止位置，第四列是检测到这个fragment的细胞barcode，第五列是这个fragment在这个细胞中检测到的的个数。这是一个类似于bed文件格式的文件。

---

peak_counts矩阵

• Signac官方教程给出peak_counts矩阵的定义是：

peak_counts矩阵类似于用于分析单细胞 RNA-seq 的基因表达counts矩阵。然而，矩阵的每一行代表基因组的一个区域（峰），而不是基因，预计代表开放染色质的区域。矩阵中的每个值代表映射在每个峰内的每个单个条形码（即一个单元格）的 Tn5 积分位点的数量。您可以在10X 网站上找到更多详细信息。

• 我们来读入counts矩阵，看看他和fragments文件的关系

peak_counts <- Read10X_h5(filename = '/local/txm/txmdata/scATAC_seq/data/Signac/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5')

• 首先看细胞barcode有什么关系

length(unique(F$V4))
# 527201   
ncol(peak_counts)
# 8728    
colnames(peak_counts)[5417] %in% F$V4
# TRUE   

• fragments中有527201（50w）个barcode，而peak_counts矩阵里只有8728个细胞， peak_counts中的所有barcode均在fragments里。说明 fragments中包含了peak_counts的所有信息。
• 接下来取出一个细胞，来看这个细胞的peak和fragment有什么关系：

peak_counts[1:11,1:4]
# 11 x 4 sparse Matrix of class "dgCMatrix"
#                    AAACGAAAGAGCGAAA-1 AAACGAAAGAGTTTGA-1 AAACGAAAGCGAGCTA-1 AAACGAAAGGCTTCGC-1
# chr1:565107-565550                  .                  .                  .                  .
# chr1:569174-569639                  .                  .                  .                  .
# chr1:713460-714823                  .                  .                  2                  8
# chr1:752422-753038                  .                  .                  .                  .
# chr1:762106-763359                  .                  .                  4                  2
# chr1:779589-780271                  .                  .                  .                  .
# chr1:793516-793741                  .                  .                  .                  1
# chr1:801120-801338                  .                  .                  .                  .
# chr1:804872-805761                  .                  .                  .                  4
# chr1:839520-841123                  .                  .                  2                  2
# chr1:841866-842572                  .                  .                  .                  2

• 取出第四列的AAACGAAAGGCTTCGC-1细胞，查看它在fragments文件中的fragments情况

AAACGAAAGGCTTCGC_fragment <- F[which(F$V4 == colnames(peak_counts)[4]),]
dim(AAACGAAAGGCTTCGC_fragment)
# 167297      5      这个细胞有167297个fragments
head(AAACGAAAGGCTTCGC_fragment,13)
#         V1     V2     V3                 V4 V5
# 1449  chr1 712868 713146 AAACGAAAGGCTTCGC-1  1
# 1893  chr1 713783 714045 AAACGAAAGGCTTCGC-1  2   属于chr1:713460-714823 \
# 2283  chr1 713944 713997 AAACGAAAGGCTTCGC-1  1   属于chr1:713460-714823 ｜ 共7个read，peak_counts矩阵中的值:8
# 3940  chr1 714144 714209 AAACGAAAGGCTTCGC-1  3   属于chr1:713460-714823 ｜
# 4603  chr1 714263 714732 AAACGAAAGGCTTCGC-1  1   属于chr1:713460-714823 /
# 6284  chr1 757378 757538 AAACGAAAGGCTTCGC-1  2
# 8887  chr1 762951 763224 AAACGAAAGGCTTCGC-1  1   属于chr1:762106-763359    共1个read，peak_counts矩阵中的值:2
# 9448  chr1 773225 773453 AAACGAAAGGCTTCGC-1  2
# 10524 chr1 793548 793847 AAACGAAAGGCTTCGC-1  3
# 10955 chr1 802879 803148 AAACGAAAGGCTTCGC-1  1
# 11403 chr1 805213 805358 AAACGAAAGGCTTCGC-1  1
# 12160 chr1 805561 805603 AAACGAAAGGCTTCGC-1  1
# 12931 chr1 823987 824159 AAACGAAAGGCTTCGC-1  2

• 总之，fragments文件是10x机器产生的原始数据，包含了所有细胞的所有的片段的信息。与fragments文件同名的.tbi文件是fragments的tabix索引，便于从任意基因组位置快速读取fragments文件。

---

• 从fragments文件得到peak_counts矩阵涉及到对细胞进行筛选（Cell Calling）以及对原始peak进行合并与筛选（Peak Calling）。

---

Peak Calling

• 目前存在五种Peak Calling方法：
  1. 直接使用数据库中的DNase-Seq和ChIP-Seq的Peak，例如：cisTopic；
  2. 使用整套数据集上的所有Read进行Peak Calling，例如：CellRanger和Cicero；
  3. 使用一个细胞系（Cell Line）上所有的Read进行Peak Calling，例如：chromVAR；
  4. 使用一个细胞类型（Cell Type）下的所有Read，例如：snapATAC；
  5. 两阶段方法，即先用其他注释（例如：Bin）得到Feature-Cell矩阵，并作聚类，然后把每一个类中所有的Read汇总起来做Peak Calling，例如：scATAC-pro和ArchR。

Cell Ranger v2.0 Peak Calling：原始移调事件用于生成具有 401bp 移动窗口总和的本地平滑信号轨道。在拟合和选择全局峰值阈值后，信号高于阈值（以绿色显示）的连续区域将作为候选峰值调用产生。 Cell Ranger v2.0 如何对候选区域中的单个假定峰值执行局部信噪比估计：信号的绿色部分显示了正在检查的推定峰值，峰值信号测量为绿色部分的中值。灰色部分被掩盖了，因为它们是其他假定的峰值，因此不用于估计局部背景。红色部分用于局部背景估计，峰值背景作为所有红色部分的中值。

---

Cell Calling

• 去除doublets和受污染的细胞等

---

参考

https://satijalab.org/signac/articles/pbmc_vignette.html
https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/algorithms/overview
https://www.jianshu.com/p/f7975da8e1e8