单细胞笔记11-scATAC-seq上游数据处理

作者: 江湾青年 | 来源:发表于2021-11-03 19:42 被阅读0次

单细胞笔记11-scATAC-seq上游数据处理
单细胞-目录索引
单细胞转录组及空间转录组学习资源
单细胞转录组及空间转录组学习资源
单细胞scRNA-seq学习笔记2-数据预处理
从10X单细胞bam文件中提取某一个细胞类型的bam
单细胞上游分析
单细胞pepline的上游是怎么构建的？
NCBI公共数据库下载10x单细胞数据进行上游分析的全流程
scRNAseq上游分析

一直以来都不是太清楚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。