ATAC基本概念
ATAC(Assay for Transposase Accessible Chromatin with hight-throughput sequencing),指的是转座酶可及的染色质区域的高通量测序。 转座酶是由转座子编码的一种酶，在NGS中用于文库构建，最常用的是Tn5转座酶，其随机性好，稳定性高，插入位点易测。通过转座酶， 只需一步反应就可以实现DNA的片段化，末端修复，接头连接，大大加快了文库构建的过程。 通过Tn5转录酶富集开放染色质区域的DNA序列，然后进行高通量测序。通过将测序得到的序列mapping回基因组，可以识别开放的区域在基因组上的具体位置，即peak区域。
ATAC peak shift(偏移)
peak calling之前必须的操作。 ATAC使用Tn5转座酶来完成文库的构建工作，Tn5转座酶在连接adapter序列时，会存在9bp的gap。为了填补gap区域，有一个5分钟的延伸过程，gap补齐之后在进行PCR扩增。在下机数据中，序列是经过了gap补齐的，不是最初始的断裂点了。利用补齐之后的下机序列，得到的peak位置和真实的切割位置之间存在了偏倚，为了还原真实的染色质开放区域，在进行peak calling之前，需要对reads比对的参考基因组位置进行偏移。
# 原理：对于正链上的reads需要向右偏移4bp, 比对的起始位置加4，对于负链上的reads, 则向左偏移5bp, 比对的起始位置减5bp。
# 操作：对于原始的bam文件，首先利用bedtools转换成bed文件，只保留reads比对上参考基因组的位置，然后再进行比对位置的偏移，
def tn5_shift_ta(ta, out_dir):
prefix = os.path.join(out_dir,
os.path.basename(strip_ext_ta(ta)))
shifted_ta = '{}.tn5.tagAlign.gz'.format(prefix)
cmd = 'zcat -f {} | '
cmd += 'awk \'BEGIN '
cmd += ' '
cmd += 'else if ($6 == "-") print $0}}\' | '
cmd += 'gzip -nc > {}'
cmd = cmd.format(
ta,
shifted_ta)
run_shell_cmd(cmd) return shifted_ta
插入片段长度分布图质控
ATAC文库，其插入片段存在典型的分布规律，每200bp会存在一个峰，这个周期性波动反应的是核小体的个数。在ATAC_seq的数据分析中，会对插入片段长度分布进行可视化，观察其是否符合这样的周期性规律，一定程度可以反映文库构建的质量。
针对双端测序，直接提取bam文件的第9列进行绘图
samtools view input.bam | \
awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)} {print $1"\t"abs($9)}' | \
sort | uniq | cut -f2 > fragment.length.txt
# R语言绘图
data <- read.table("fragment.length.txt", header = F)
# 设置插入片段长度的阈值，过滤掉太长的片段
length_cutoff <- 1200
fragment <- data$V1[data$V1 <= length_cutoff]
# 利用直方图统计频数分布，设置柱子个数
breaks_num <- 500
res <- hist(fragment, breaks = breaks_num, plot = FALSE)
# 添加坐标原点
plot(x = c(0, res$breaks),
y = c(0, 0, res$counts) / 10^2,
type = "l", col = "red",
xlab = "Fragment length(bp)",
ylab = expression(Normalized ~ read ~ density ~ 10^2),
main = "Sample Fragment sizes")
ATAC文库质控(ENCODE)
ATAC测序饱和度分析
ATACseqQC包针对ATAC文库进行QC。 ATACseqQC的饱和度分析思路如下，对effectice fragment进行随机抽样，比如10%，20%到100%， 对于每个梯度分别进行peak calling, 统计其peak个数，然后以effectice fragment数目为横坐标，peak个数为纵坐标绘制散点图，并进行拟合，如果拟合曲线的末端斜率上升平缓，表明继续加大测序数据量，新识别到的peak个数也不会新增太多，就可以认为当前的测序数据量是比较饱和的，不比再加测数据量了。如果曲线末端仍处于一个斜率急剧上升的阶段，说明测序数据量不够，还需加大测序数据量。
## 第一步：对effectice fragment进行抽样
## 所谓effectice fragment，即过滤掉PCR重复，去除了线粒体之后的有效序列。在实际分析中，质控后的clean reads首先比对到参考基因组上，生成原始的bam文件，也称为raw bam。然后对原始的bam文件进行过滤，去除MAPQ较低的alignment, 去除PCR重复，去除比对到线粒体的序列，经过重重筛选，得到了一个filter之后的bam文件，这个filter bam文件中保存的就是effectice fragment。
## 提取sam文件的header部分作为一个单独的文件；
## 使用linux内置的shuf命令对sam文件的正文部分随机抽取一定比例的行数，追加到header文件中，生成一个新的sam文件
bash downsample.sh GL1.bam\ # 利用文章附带的脚本
## 利用排序好的bam文件可以直接进行peak calling, 然后统计每个抽样百分比对应的effectice fragment数目和peak个数，就可以绘制测序饱和度图了。在ATACseqQC中，提供了saturationPlot函数，可以直接读取一系列peak文件，然后绘制测序饱和度图
PCR重复序列判定
- samtools
输入文件为经过samtools fixmate命令处理、按照比对上染色体坐标位置排序之后的bam文件
# 第一步，按照read name排序bam文件 samtools sort -n -o namesort.bam input.bam # 第二步，运行fixmate命令 samtools fixmate -m namesort.bam fixmate.bam # 第三步，按照coordinate排序bam文件 samtools sort -o positionsort.bam fixmate.bam #第四步，运行markdup命令 samtools markdup positionsort.bam markdup.bam - picard MarkDuplicates
输入文件为按照比对位置排序后的bam文件
# 第一步，按照coordinate排序bam文件 samtools sort -o positionsort.bam input.bam # 第二步，运行MarkDuplicate命令 java -jar picard.jar MarkDuplicate \ I=positionsort.bam \ O=markdup.bam \ M=markdup.metrc.csv - sambamba
# 第一步，按照coordinate排序bam文件 sambamba sort -o positionsort.bam input.bam # 第二步，运行markdup命令 sambamba markdup positionsort.bam markdup.bam
参考链接
生信修炼手册(1) 生信修炼手册(2) Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position 生信修炼手册(3)
