ATAC-seq 数据分析实战

一、 ATAC-seq原理和基础知识

1. ATAC-seq原理

真核生物的DNA并不是裸漏的，而是组蛋白和染色体/染色质结合。DNA一圈一圈的缠绕在8个组蛋白上，形成核小体。一个个核小体构成串珠式的结构，然后进一步折叠、聚合，并在其他架构蛋白的协助下，形成染色体。经过一系列操作就将超长的DNA链，折叠成很小的结构，塞进小小的细胞核内。

基因的转录，需要将DNA的高级结构打开，但是不需要DNA链全部解开，只需要打开一部分，也就是基因表达的区域解开即可。这一过程，主要由染色体组蛋白的修饰（尤其是乙酰化）来实现的。这部分打开的染色质，就叫做开放染色质（染色体和染色质是同一种物质的两种形态，染色质是伸展状态，染色体是高度螺旋的状态）。而染色质一旦打开，就允许一些调控蛋白（比如转录因子）跑过来与之结合。而染色质的这种特性，就叫做染色质的可及性，所以说染色质的可及性反应的是调控因子与开放染色质结合的状态，与转录调控密切相关。
ATAC-seq是如何找开放染色质区域的呢？
使用了转座酶Tn5：DNA转座是一种由DNA转座酶介导，把DNA序列从染色体的一个区域插入到另一个区域的现象，类似”粘贴复制“。这个过程也是需要插入位点的染色质是开放的。
既然转座酶Tn5容易结合在开放染色质上，只要人为的将NGS接头连接到转座酶，携带这些接头的转座酶进入细胞核后，切开染色质开放区域，使染色质断裂并将这些接头插入到开放的染色质区域中，这样裂解细胞、破碎DNA后，利用已知序列的测序标签进行NGS测序，就知道哪些区域是开放区域了。

2. Tn5转座子

1. 转座概念

可移动的DNA片段即可移动因子在基因组上自由转移称为转座，DNA与所插入的基因位点可以是非同源的。转座是产生基因多样性的重要机制，可移动因子可产生插入、缺失、倒置以及染色体融合突变。
转座需要通过转座酶来催化。原核生物的转座分为两种方式，复制转座和保守转座：

• 复制转座的供体DNA完整，把通过复制的DNA片段插入基因位点上
• 保守转座则是从供体DNA上分离一段DNA，以转座酶为中介，连接到目标DNA上而实现的

2. 参与分子

• 转座子（Transposon） ： 可移动DNA片段
• 转座酶（Transposase / TNP）：催化转座的蛋白质；野生型Tn5转座酶是一种活性极低的蛋白质
• 目标DNA（Target DNA）： 可以与转座子在同一个DNA分子上，甚至转座子内；或在另一个DNA分子上

1. 转座子

（1） 简化的转座子结构

包含合成Tnp的DNA序列，两个19bp长的末端以及任意DNA序列。

• 末端是两个19bp长的片段，将Tnp和任意DNA序列包含在其中。
• 常见的末端有三种：外末端（outside end / OE），内末端（inside end / IE）和镶嵌性尾端（mosaic end / ME）。组合方式有两个反向的OE，或者两个反向的IE，或两个反向的ME，又或是两组反向的 OE和IE组合

（2） Tn5转座子的结构

Tn5转座子由两个反向的插入片段 IS50 以及两组 OE 和 IE 构成

• IS50 包括三个抗生素抗性基因。 IS50R 负责编码 Tnp 和转座抑制物（Inh），而 IS50L 负责编码两个低活性蛋白

IS（insertion sequence）： 插入序列，很小（< 2.5 kb）DNA片段，可以在不同的基因位点跳跃，或自我复制。通常存在于细菌与古细菌基因中，但也存在于真核生物的转座元素中。编码的基因一般只与移动有关。

2. 转座酶

• Tn5 Tnp是一种转座酶，可以将DNA片段从一个位置移动到另一个位置，来自大肠杆菌，全长477个氨基酸。
• Tn5 Tnp可以与特异性DNA识别和结合，特异性DAN是指Tn5或IS50的末端反向重复序列。
• Tn5 Tnp的主要功能区有三个，N末端、催化结构域和C末端：
  1. N末端是特异性结合DNA结构域，可以识别和结合Tn5或IS50的末端反向重复序列
  2. 催化结构域是转座反应的核心，可以切割和连接DNA，并形成双聚体
  3. C末端是合成复合体的必需部分，可以促进Tn5 Tnp之间的相互作用，并影响转座效率

3. 转座过程

Tn5转座对目标DNA的特异性要求不高，可以插入到任何双链DNA上。但是，Tn5也有一些偏好性，比如倾向于插入到AT富的区域，或者靠近某些特定的序列。Tn5转座酶（Tnp）的突变也可以改变其对目标DNA的结合特异性和亲和力

3. fastq文件格式

fastaq格式是一种基于文本的存储生物序列和对应碱基质量的文件格式。
下面为一个illumina平台测序的真实数据，其中包含了一条reads的信息

@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,

FASTQ格式存储的序列信息，一条reads信息可以分为四行：

1. 第一行主要存储序列测序时的坐标等信息
   指测序仪上的物理位置，用来标识不同的读段，通常有：
   1. 流动池（flowcell）编号
   2. 流动池通道（lane）编号
   3. 照片编号
   4. 照片中X坐标
   5. 照片中Y坐标
   @ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
   其中@是开始的标记符号
   ST-E00126:128:HJFLHCCXX是测序仪唯一的设备名称
   2是lane的编号
   1101是照片的编号
   7405是在照片中的X坐标
   1133是在照片中的Y坐标
2. 第二行是测序得到的序列信息，一般用ATCGN来表示，其中N表示荧光信号干扰无法判断具体碱基
3. 第三行是以+号开始，用来存储一些附加信息，一般是空的
4. 第四行是质量信息，与第二行的碱基序列是一一对应的，每一个符号对应的ASCII值可以理解为对应位置的碱基的质量值，越大说明测序的质量越好，不同版本对应的而不同
   在测序仪进行测序时，会自动根据荧光信号的强弱给出一个参考的测序错误概率P，为了节省存储这个信息的空间，所以进行了一下操作：
   1. 将P取 log10 再乘以 -10，得到Q，如果P=1%，则Q=20
   2. 再将这个Q加上33或64称为Phred，再将Phred对应到ASCII字符上。如Q=20，则Phred为53，对应符号为 ‘5’

4. shell获取文件每行长度

cat test.sh | awk '{print length($0)}'

5. shell 字符串的匹配

批量操作通常 ： ls | grep 筛选 | while read id;do 操作 ;done
在while内部的操作中，通常会用到文件名，则可以使用字符串匹配，例如文件名为 kaikobase.1.ht2

文件名： id=kaikobase.1.ht2
		用. 来分割字符串，如果是用其他字符分割，则下面例子中.应该相应改变
1. 左删除_最小匹配
${id#*.}
结果为： 1.ht2
2. 左删除_最大匹配（贪婪匹配）
${id##*.}
结果为： ht2
3. 右删除_最小匹配
${id%.*}
结果为： kaikobase.1
4. 右删除_最大匹配（贪婪匹配）
${id%%.*}
结果为：kaikobase

总结：1. *表示要删除的字符
	 2. 左删除是 #号, 右删除是 % 号 （记忆：在键盘位置上，#在$左边，所以是左删除；%在$右边，所以是右删除）
	 3. 删除的最小和最大匹配是指从左边或右边数第一个或最后一个符号。
	 4. 例如 ${id##*_} 表示左删除，最大匹配。即从左边第一个字符开始至字符串中最后一个_符号，都会被删除

二、数据比对和过滤

1. 比对

1. hisat2

HISAT2是一款用于将二代测序数据（DNA和RNA）比对到基因组数据的快速和敏感的比对软件。相比于Bowtie/TopHat2等软件，HISAT2具有更高的敏感性和更快的运算速度，这得益于其优化了索引建立的策略，并采用了新的比对策略。