单细胞ATAC分析标准流程（详细版）

背景：单细胞ATAC分析流程介绍。

1.为什么做单细胞ATAC分析？

单细胞转录组分析的是单细胞分辨率上不同细胞之间的转录组（mRNA）差异，由中心法则可以知道，生物发生的过程，从基因组、转录到蛋白质翻译的过程产生的任何差异都可能是细胞类型异质性的原因。ATAC的全称是Assay for Transposase Accessible Chromatin with high-throughput sequencing，从名字上可以知道，这是使用高通量研究染色质开放性的方法，实际上，研究染色质开放性的方法不止一种，而ATAC指的是使用Tn5转座酶研究染色质开放性。放在单细胞分辨率上，ATAC就成为了scATAC，sc即single-cell（单细胞）的意思。

2.scATAC的原理是什么？

很容易知道，scATAC研究对象是染色质，其属于基因的上游调控过程，对应于中心法则的第一步基因组位置。在基因的上游调控中，存在很多基因组的表观修饰，它们并没有改变染色质或者DNA的结构，却能对基因的表达产生一定的调控（增强或者减弱表达等）作用。表观修饰除了染色质可及性以外，还包括之前介绍过的DNA甲基化等。
那么染色质可及性的研究原理是什么呢？
从下图中可以看到：

image.png
这是一条染色体，其中部分基因组与组蛋白缠绕形成了核小体，部分基因组不与组蛋白缠绕形成了染色质开放性区域。
很容易理解，当基因组被组蛋白占据时候，它就不能被别的蛋白占据（即TF蛋白等），反过来，如果它处于开放性状态，那么它就能被其他的调控蛋白（TF）结合，从而影响近端甚至是远端基因的表达。
所以一定程度上讲，scATAC研究的关键就是染色质开放性影响的TF作用情况。
那么TF作用于基因组的条件是什么呢？这里总结2点：
（1）基因组上存在TF结合位点，即motif
（2）染色质处于开放性状态，可以被TF蛋白结合
整个调控过程应该是：
a.某个区域开放状态、某个区域存在转录因子X的motif-->b.转录因子X可以结合于这个区域、这个区域位于基因A的上游-->c.转录因子X可以调控基因A
另外注意，染色质开放区域中可能存在基因的启动子、增强子、沉默子和绝缘子等元件。

3.scATAC的实验原理是什么？

scATAC是基于Tn5转座酶并进行高通量测序。Tn5转座酶可以随机结合并切割开放性的DNA，并在切割位点插入其自身携带的二代测序接头。
在10x中导入样本必须是处理过的细胞核悬液，Tn5转座酶可以渗透入细胞核膜。实际上，对于所有的细胞核悬液，都有一些质控要求，例如核膜要完整、气泡尽量少，这很容易理解，如果核膜破损，很容易造成污染。
对于测序结果，有几个关键的名词：
（1）fragments：基于个人理解，这应该是Tn5转座酶切割到的所有序列。所有序列是什么意思呢？实际上，在Tn5切割的时候，得到的序列存在多种情况：仅仅是开放区域、包含一个核小体的开放区域、包含多个核小体的开放区域...
（2）peaks：基于个人理解，这里的peaks就是单纯的开放区域，没有包含核小体
经过10x测序和cellranger-arc分析后，得到的peaks矩阵代表了某个细胞某个peak的丰度。

4.scATAC的质控原理是什么？

对于任何一个组学来说，都会在具体的下游分析之前进行质量控制。scATAC也有一些质控的指标。
首先介绍一下黑名单，谈过恋爱的都知道，黑名单不是一个好的词语，在这里也是一样，ENCODE黑名单区域（blacklisted region）包含一些异常的、在二代测序中信号很高的区域，通常需要将处于这些区域的序列过滤；
其次，线粒体基因组因为没有染色质包装，更容易被Tn5酶切割，在某些情况下需要将线粒体序列过滤。
最后，就是配对错误和质量比较低的reads需要过滤。
那么，什么是质量比较低的reads呢？在scATAC的发展过程中，使用了以下的指标：
如下图所示：

image.png

它代表的是不同长度fragments的分布情况。前文介绍过，fragments代表所有的序列，其中，没有核小体的开放区域（nucleosome-free region, NFR）通常是<100bp，单个核小体区域_{200bp，双核小体}400bp，三个核小体区域~600bp。于是很容易从上图看出，主峰的长度大约是100多bp，对应了无核小体区域。
补充一点，一个核小体结合大约147bp（我看也有的文献说是146bp）的DNA。
此外，如下图所示：

image.png
它代表的是转录起始位点距离与信号强度的分布关系，这也很容易理解，NFR的片段应该富集在转录起始位点（TSS）附近，而结合核小体的区域在TSS位置应该缺失，在TSS两侧相对富集。
整个原理如下图所示：
image.png
补充TSS知识点：
TSS指的是转录起始位点：转录起始位点通常为一个基因的5`端转录的第一个碱基，它是与新生的RNA链第一个核苷酸相对应的DNA上的碱基，通常为一个嘌呤（A或者G）。TSS是基因转录的起点，通常也是启动子的一部分。
（1）TSS富集分数：所有基因TSS位点测序深度的平均值
参考基因组TSS位点文件：
image.png
由上图可以看到，每一行TSS都是一个位点的碱基。常规ATAC的TSS阈值如下图（不确定单细胞ATAC是否也一致）：
image.png
（2）另一种计算TSS富集分数
更常规的做法是提取TSS两侧区域的reads, 比如上下游各提取2kb的区域，划分等长窗口bin并统计每个窗口内的coverage, 最终会生成一个矩阵，行为基因，列为不同窗口，根据这个矩阵可以绘制TSS两侧reads分布图， 也可以计算TSS Enrichment score。
image.png

5.scATAC下游分析的原理是什么？

始终要记住，scATAC分析的重点就是开放区域的motif（转录因子结合位点）。

1. peak calling
   不论是bulkATAC还是scATAC，质控完成的第一步就是peak calling。
   peak calling就是识别开放性区域，通常使用的工具为MACS2分析工具（烈冰也是使用这个），当然，signac和ARCHR也是调用的这个工具。在常规的ATAC分析中，HMMRATAC（另一个软件）通过半监督隐马尔科夫模型把基因组分成高信号的活性染色质区域、中等信号强度的核小体区域和低信号的背景区域。
2. 差异peak分析
   其实在这一步之前，还包括数据的降维、聚类以及细胞类型注释等，请注意，scATAC特有的降维方式是LSI降维，这是由它的数据结构决定的。
   当细胞被分为不同的细胞类型后，与单细胞转录组类似，可以使用findMarkers()函数分析细胞类型之间的差异peaks。
   这一步的意义就是研究染色质可及性的差异对不同细胞类型的影响，就像单细胞转录组研究基因表达差异对不同细胞类型影响一样。
3. peak注释
   peak根据其最近的基因或者调控元件进行注释，之后可以用Gene Ontology（GO）、KEGG和Reactome等数据库做功能富集分析。
4. motif富集分析（重点）
   motif在这里代表的是TF对应的DNA结合位点，所以它应该是测序DNA序列上的一段，注意是测序数据本身的DNA序列。
   motif的富集需要使用TF数据库JASPAR2020，主要是寻找peaks中的motif序列。
   注意并不是所有的TF都有鉴定好的motif结合序列，motif的富集并不代表相应的TF已经结合，只是证明可以结合。
5. 转录因子的footprint（重点）
   得到了motif的富集，只能说明该位置可以结合TF，但实际上是否结合并不知道，因此需要做转录因子的footprint分析。footprint的原理是：一个TF结合在DNA上，阻止Tn5切割，在染色质开放区域留下一个相对缺失的位置。
   image.png
   如上图所示，当motif位置出现了一个低谷，那么就证明该motif序列已经被转录因子所结合。那么实际应用情况如何呢？如下图：
   image.png
   从图中可以看到，有两种细胞类型：CD8激活细胞和CD8祖细胞，它们在motif（EOMES和TBX21)位置出现了低谷，而CD8激活细胞在低谷两侧的Tn5富集分数更大，低谷更明显，说明在这个位点的TF结合到了motif上的可能性更大，这也是单细胞的footprint差异分析。
6. 染色质共可及性
   共可及性的原理是：利用单细胞染色质可及性数据预测基因组中更有可能在细胞核的物理上接近的区域，得到多个细胞类型的两个峰值之间的可及性关系。换句话说，当A峰在单细胞内可访问时，B峰通常也可访问，这产生了很强的相关性，但并不一定意味着这些峰值之间存在调控关系。
   image.png
   图中第四行Links即代表peaks之间的相关性。
7. 伪时间分析
   需要注意的是，单细胞转录组的伪时间分析可以刻画基因在时间过程中的表达变化，而scATAC能够刻画时间过程中的染色质开放性变化和motif活性变化。
   下图是motif轨迹：
   image.png
   下图是染色质可及性轨迹:
   image.png