一直以來都不是太清楚scATAC-seq中fragments文件延赌、peak calling，以及counts矩陣的關(guān)系叉橱，本文就來梳理一下scATAC-seq的上游數(shù)據(jù)處理中的一些基本概念挫以。

fragment

• fragment是Tn5轉(zhuǎn)座酶隨機(jī)插入到開放染色質(zhì)中剪切出的DNA片段，再經(jīng)過測序和mapping到基因組上就得到了fragments文件窃祝。

Tn5轉(zhuǎn)座酶剪切開放DNA示意圖

• Tn5轉(zhuǎn)座酶以同源二聚體的形式結(jié)合到DNA上掐松，在兩個Tn5分子間隔著9bp的DNA序列。根據(jù)這個情況粪小，每個Tn5同源二聚體的結(jié)合事件會產(chǎn)生2個Insertions大磺，中間隔著9bp。因此探膊，真實的"開放"位置的中心在Tn5二聚體的正中間杠愧，而不是Tn5的插入位置。為了盡可能的還原真實情況逞壁，一般會對Tn5的Insertions進(jìn)行了校正流济，即正鏈的插入結(jié)果往右移動4bp(+4 bp), 負(fù)鏈的插入結(jié)果往左偏移5bp(-5 bp)锐锣。

fragments文件

• Signac官方教程給出fragments文件的定義是：

fragments文件表示所有單個單元格中所有獨特片段的完整列表。它是一個相當(dāng)大的文件袭灯，處理速度較慢刺下，并且存儲在磁盤上（而不是內(nèi)存中）。但是稽荧，保留此文件的優(yōu)點是它包含與每個單個單元格關(guān)聯(lián)的所有片段，而不是僅映射到峰值的片段工腋。有關(guān)片段文件的更多信息可以在 10x Genomics網(wǎng)站或sinto 網(wǎng)站上找到姨丈。

• 以Signac官方給出的pbmc_10k的數(shù)據(jù)為例，來看看fragments文件和peak counts矩陣到底是什么關(guān)系擅腰。首先讀入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的細(xì)胞barcode沐绒，第五列是這個fragment在這個細(xì)胞中檢測到的的個數(shù)。這是一個類似于bed文件格式的文件旺坠。

peak_counts矩陣

• Signac官方教程給出peak_counts矩陣的定義是：

peak_counts矩陣類似于用于分析單細(xì)胞 RNA-seq 的基因表達(dá)counts矩陣乔遮。然而，矩陣的每一行代表基因組的一個區(qū)域（峰）取刃，而不是基因蹋肮，預(yù)計代表開放染色質(zhì)的區(qū)域。矩陣中的每個值代表映射在每個峰內(nèi)的每個單個條形碼（即一個單元格）的 Tn5 積分位點的數(shù)量璧疗。您可以在10X 網(wǎng)站上找到更多詳細(xì)信息坯辩。

• 我們來讀入counts矩陣，看看他和fragments文件的關(guān)系

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

• 首先看細(xì)胞barcode有什么關(guān)系

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

• fragments中有527201（50w）個barcode崩侠，而peak_counts矩陣?yán)镏挥?728個細(xì)胞漆魔， peak_counts中的所有barcode均在fragments里。說明 fragments中包含了peak_counts的所有信息啦膜。
• 接下來取出一個細(xì)胞有送，來看這個細(xì)胞的peak和fragment有什么關(guān)系：

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細(xì)胞，查看它在fragments文件中的fragments情況

AAACGAAAGGCTTCGC_fragment <- F[which(F$V4 == colnames(peak_counts)[4]),]
dim(AAACGAAAGGCTTCGC_fragment)
# 167297      5      這個細(xì)胞有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機(jī)器產(chǎn)生的原始數(shù)據(jù)八拱，包含了所有細(xì)胞的所有的片段的信息阵赠。與fragments文件同名的.tbi文件是fragments的tabix索引涯塔，便于從任意基因組位置快速讀取fragments文件。

• 從fragments文件得到peak_counts矩陣涉及到對細(xì)胞進(jìn)行篩選（Cell Calling）以及對原始peak進(jìn)行合并與篩選（Peak Calling）清蚀。

Peak Calling

• 目前存在五種Peak Calling方法：
  1. 直接使用數(shù)據(jù)庫中的DNase-Seq和ChIP-Seq的Peak匕荸，例如：cisTopic；
  2. 使用整套數(shù)據(jù)集上的所有Read進(jìn)行Peak Calling枷邪，例如：CellRanger和Cicero榛搔；
  3. 使用一個細(xì)胞系（Cell Line）上所有的Read進(jìn)行Peak Calling，例如：chromVAR东揣；
  4. 使用一個細(xì)胞類型（Cell Type）下的所有Read践惑，例如：snapATAC；
  5. 兩階段方法嘶卧，即先用其他注釋（例如：Bin）得到Feature-Cell矩陣尔觉，并作聚類，然后把每一個類中所有的Read匯總起來做Peak Calling芥吟，例如：scATAC-pro和ArchR侦铜。

Cell Ranger v2.0 Peak Calling：原始移調(diào)事件用于生成具有 401bp 移動窗口總和的本地平滑信號軌道。在擬合和選擇全局峰值閾值后钟鸵，信號高于閾值（以綠色顯示）的連續(xù)區(qū)域?qū)⒆鳛楹蜻x峰值調(diào)用產(chǎn)生钉稍。

Cell Ranger v2.0 如何對候選區(qū)域中的單個假定峰值執(zhí)行局部信噪比估計：信號的綠色部分顯示了正在檢查的推定峰值，峰值信號測量為綠色部分的中值携添〖廾ぃ灰色部分被掩蓋了，因為它們是其他假定的峰值烈掠，因此不用于估計局部背景羞秤。紅色部分用于局部背景估計，峰值背景作為所有紅色部分的中值左敌。

Cell Calling

• 去除doublets和受污染的細(xì)胞等

參考

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