ATAC-Seq 分析流程

本文分享我的 ATAC-Seq 分析流程，为了控制篇幅，关键的 Peak calling 步骤会详细介绍，而 Motif 分析等一些步骤则拆分到以后文章，在未来几天推送。虽然如此，篇幅还是挺长的，建议一次读完。

简单来说，基因的转录需要 DNA 从缠绕在组蛋白的紧密状态释放，成为松散的开放区域，这样转录复合体才能结合上去。ATAC-Seq 技术能高通量分析开放的染色质区域，揭示基因的转录状态。

只有开放的染色质区域 Tn5 才能够携带测序接头结合并将接头转座上去，PCR 扩增后就可以进行测序。

基本的实验要求

一、实验最好有生物学重复，无论什么技术都会有一定的随机误差。理论上还需要一个实验的技术对照，类似于 Chip-Seq 实验的 Input 对照。即用普通的破碎法(如超声)打断 DNA 然后同样加接头测序，作为与转座法的对照，确定基因组哪些区域是更难进行测序的。一般不做这个对照实验。
二、依据不同的物种和估计的染色质开放程度改变，一般人种要求 50 M 比对上 Reads，TF footprinting 分析要 200 M 以上 Reads. 合格的实验唯一比对率应该在 80% 以上。

For inferring differences in open chromatin within human samples we generally use >50M mapped reads and for transcription factor foot-printing we use >200M mapped reads (Neph et al., 2012).

本教程主要的分析步骤是质控、比对、Peak calling、差异分析、注释、Motif finding和可视化。下面一一给出步骤示范。

fastq 质控

先用 fastqc 检查质量，fastp 进行过滤和去接头。有人认为除了去接头，不额外进行其他截短之类操作，防止影响后续去重。也有人认为根据染色质开放区域实际片段大多在 150 bp 以下，如果是太长的测序，应该进行截短，如统一截短到 75bp.

典型的 ATAC-Seq 片段分布，小片段占比高

在我这次数据我进行了测试，统一截短到 75bp 在 fastqc 对单端检测重复率下降 1% 左右，用 fastp 对双端的检测重复率不变。同时使用 Bowtie2 的比对率也没变化或只有轻微改变(< 1%)。 也有可能因为过短的插入片段因为导致 Reads 3' 端太多 G 碱基且低质量，早被 fastp 过滤了，所以无法体现出影响。总的来说，自己的数据 2 方案都试试看看情况决定。

比对和过滤

用 BWA, bowtie2 等常规软件比对。这里 -k 设置一条 reads 至多报告的匹配次数，默认是多个匹配报告最佳；-X 设置插入片段最大长度，默认 500. 同时用 samtools 排序输出。
哈佛教程推荐用 -k 10 是因为他们的 Genrich 支持利用多比对信息。如果用其他软件做 Peak calling 还是默认不用 -k 操作。使用 Genrich 要求用 -n 进行排序。

  bowtie2 --very-sensitive -k 10 -X 2000 --threads 4 -x ${GRCh38} \ 
  -1 ${CleanDir}/${Sample}_R1.fq.gz -2 ${CleanDir}/${Sample}_R2.fq.gz \ 
  -S ${BamDir}/${Sample}.sam

  samtools view -b -h ${BamDir}/${Sample}.sam | samtool sort -n -O BAM \ 
  -o ${BamDir}/${Sample}_sorted.bam

用 picard 对 PCR 重复进行标记。后续用 Genrich 的话不需要标记，Genrich 自带移除 PCR 重复参数。

gatk MarkDuplicates -I ${BamDir}/${Sample}_sorted.bam -M ${BamDir}/${Sample}_dup_matrics.txt -O ${BamDir}/${Sample}_MD.bam

顺便可以看看插入片段分布情况。插入片段分布应该是随着长度下降的，但是在 <100,200,400,600bp 处有峰，分别对应着 NER(nucleosome-free regions) 和 单、双、三核小体。因为一个核小体缠绕的 DNA 长度约为 200bp.

 gatk CollectInsertSizeMetrics -H ${BamDir}/${Sample}_InsertSize.pdf -I ${BamDir}/${Sample}_MD.bam -O ${BamDir}/${Sample}_InsertSize.txt

测序插入片段形成典型的计数分布原因 一个样本的实际片段分布案例

线粒体 DNA 与核基因组不同，他没有缠绕组蛋白形成核小体，因此都是开放的。一般来说会导致测序结果里有不少线粒体信号，需要移除。同时移除基因组 Blacklist 区域。
移除后过滤低质量的比对， Bowtie2 设置比对质量大于 30. 要注意不同比对软件对 MAPQ 值设置不同。同样，如果使用 Genrich 进行 Peak calling 那么这些步骤不需要进行，他有相应参数实现。

  bedtools intersect -nonamecheck -v -a ${BamDir}/${Sample}_MD.bam -b ${Blacklist} > ${BamDir}/${Sample}_RMBL.bam
  samtools view -h -F 1804 -q 30 -O SAM ${BamDir}/${Sample}_RMBL.bam | grep -v -e "^chrM" | samtools sort -O BAM -o ${BamDir}/${Sample}_Filtered.bam 
  samtools index ${BamDir}/${Sample}_Filtered.bam 

移动 Reads 位置

如下图所示 Tn5 酶结合和切割的位点是开放染色质区域，所以我们分析出来的峰中心应该在这个切割位点，也就是 reads 的 5' 端。

Tn5 结合的地方就是测序读长开始的地方

如下面示意图，Tn5 建库时会插入 9bp 的序列。为了正确表示峰区域位置，我们需要将比对到正链（+）的 reads 移动 +4 bp, 比对负链（-）的 reads 移动 -5 bp.

Tn5 建库示意图

这个移动可以用 deeptools 的 alignmentSieve 命令完成。

alignmentSieve --numberOfProcessors 8 --ATACshift -b ${BamDir}/${Sample}_Filtered.bam -o ${BamDir}/${Sample}_Shift.bam

同样，Genrich 会完成这个移动。如果不做 TF footprinting 分析，这个影响不大，可以不进行这步操作。

FRiP 计算

Fraction of reads in peaks (FRiP) 指落入峰区域的 reads 占比，应该要高于 0.3 较好。这个指标也不用硬性要求，所以计算也不需要那么严谨，简单化处理。

samtools view -c ${BamDir}/KO1_ATAC_ShiftSort.bam
bedtools intersect -u -a ${BamDir}/KO1_ATAC_ShiftSort.bam -b ${PeakDir}/KO1_ATAC_peaks.narrowPeak | samtools view -c

得到数值分别是 24634649, 1217509 除一下得出 0.049 的 FRiP.

Peak calling

这里分别用 Genrich 和 MACS2 进行 Peak calling. Genrich 更加方便，把前面那些过滤等步骤都集合在一起了，同时他可以利用生物学重复样本，计算出显著的开放区域。不需要手动去分开每个样本进行 Peak calling 然后合并。MACS2 是很常用的 Chip-Seq 软件。

Genrich

这里 -j 表示 ATAC 模式，一定要指定。多个生物学重复样本用逗号分隔。-c 参数是前文介绍的实验技术对照样本，一般 ATAC-Seq 没有，所以不指定。
Genrich 集合了移除 PCR 重复、线粒体区域、Blacklist 区域等功能，-r 移除 PCR 重复，-e chrM 移除线粒体区域，-E bed 移除 bed 文件区域。

Genrich -t ${BamDir}/KO1_ATAC.bam,${BamDir}/KO2_ATAC.bam -o ${GenrichDir}/KO.narrowPeak -f ${GenrichDir}/KO_pq.bed -k ${GenrichDir}/KO_pileup_p.bed -b ${GenrichDir}/KO_reads.bed -r -e chrM -E ${Blacklist} -m 30 -j

Genrich -t ${BamDir}/WT1_ATAC.bam,${BamDir}/WT2_ATAC.bam -o ${GenrichDir}/WT.narrowPeak -f ${GenrichDir}/WT_pq.bed -k ${GenrichDir}/WT_pileup_p.bed -b ${GenrichDir}/WT_reads.bed -r -e chrM -E ${Blacklist} -m 30 -j

参数 -m 过滤比对 MAPQ 值，Bowtie2 比对往往取 30 以上认为是好的比对，不同软件 MAPQ 值设定不相同，不可以等同使用。

-b 生成的 fragments 文件需要用 tabix 建立索引才能给 IGV 使用。

bedtools sort -i KO_reads.bed | bgzip > KO_reads_sorted.bed.gz
tabix -p bed KO_reads_sorted.bed.gz

bedtools sort -i WT_reads.bed | bgzip > WT_reads_sorted.bed.gz
tabix -p bed WT_reads_sorted.bed.gz

参数 -o 输出 narrowPeak 格式结果文件，表示每个峰区域的范围和强度。可以在开头添加 track type=narrowPeak 后用 USUC genome browser 打开。narrowPeak 格式也是一种 bed 格式。

chr1    9946    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180045  180182  peak_1  1000    .       403.740967      6.181179        -1      56
chr1    180465  181043  peak_2  1000    .       4381.033203     16.855917       -1      330

参数 -f 输出的类 bedgraph 格式文件保存汇总的 p,q 值。如果有样本重复，会给每个重复的 p 值和合并后的 p 值，以及是否显著差异。
两样本重复的举例

chr     start   end     -log(p)_0       -log(p)_1       -log(p)_comb    signif
chr1    0       9944    0.000000        0.000000        0.000000
chr1    9944    9945    0.252055        0.574477        0.363661
chr1    9945    9946    1.027748        1.432146        1.636151
chr1    9946    9947    1.394456        2.120604        2.556318        *
chr1    9947    9950    2.183553        2.826582        3.911967        *

参数 -k 输出类 bedgraph 格式文件，会用一行注释开始表示是哪个样本，背景对照是什么。把这个文件修改成 bedgraph 格式，就可以用 genome browser 去可视化。下面是注释行。

$ grep "experimental file" KO_pileup_p.bed
# experimental file: /Example/KO1_ATAC.bam; control file: NA
# experimental file: /Example/KO2_ATAC.bam; control file: NA

下面是开头 5 行内容。

# experimental file: /Example/KO1_ATAC.bam; control file: NA
chr     start   end     experimental    control -log(p)
chr1    0       9944    0.000000        0.900041        0.000000
chr1    9944    9945    0.500000        0.900041        0.252055
chr1    9945    9946    2.000000        0.900041        1.027748

下面是第五列背景数值总结。

$ awk '{print$5}' KO_pileup_p.bed | sort | uniq
0.000000
0.866656
0.900041
control

$ awk '{print$5}' WT_pileup_p.bed | sort | uniq
0.000000
2.114438
2.238139
control

MACS2

不像 Genrich 直接采用生物学重复数据合并分析，MACS2 进行一个个样本分析，后续再合并峰区域。

ATAC-seq 要用 --nomodel 参数，-g 要根据自己测序物种选择。

macs2 callpeak -t ${BamDir}/${Sample}_Shift.bam -n ${Sample} --outdir ${MACS2Dir} \ 
-f BAMPE -g hs --nomodel --keep-dup all --cutoff-analysis --bdg

MACS2 peak calling 结果里 \*_peaks.xls 是可以用 Excel 打开的 peak 结果文件，注意这里坐标以 1 开始而不是 0. \*_peaks.narrowPeak 是 narrowPeak 格式结果；*_summits.bed 是包含顶峰（Peak summit）的位置信息的 BED 格式文件，适合用于 motif 的分析。

对于有重复的实验来说，也许想要将得到的 peak 进行合并，可以用 bedtools merge 命令。下面的命令顺便去除了不想要的部分结果，只保留常染色体和性染色体。

cat ${Macs2Dir}/WT1_ATAC_peaks.narrowPeak ${Macs2Dir}/WT2_ATAC_peaks.narrowPeak \ 
| awk -v FS="\t" -v OFS="\t" 'length($1)<=5' | sort -k1,1 -k2,2n > ${Macs2Dir}/WT_SortPeaks.bed
bedtools merge -i ${Macs2Dir}/WT_SortPeaks.bed  > ${Macs2Dir}/WT_MergePeaks.bed

这样合并等于说取重复的并集，是没办法的简单化方法，个人不太建议。不过如果使用 DiffBind 进行差异的峰分析，他也会采用这样的合并方法。所以像 Genrich 在 peak calling 阶段直接考虑实验重复是更为合适的。

差异分析

差异分析可以采用 Peak 结果直接用 bedtools 简单化处理，得到多组的峰共同或特有区域。个人觉得下面的命令还是多调整 -f 的值试试，默认的也许不适合。

# 取共同交集
bedtools intersect -f 0.2 -F 0.2 -e -a ${GenrichDir}/KO.narrowPeak -b ${GenrichDir}/WT.narrowPeak > ${GenrichDir}/CommonPeak.bed
# KO 特异，但是有重叠则移除整条记录
bedtools subtract -A -f 0.3 -a ${GenrichDir}/KO.narrowPeak -b ${GenrichDir}/WT.narrowPeak > ${GenrichDir}/KO_Specific.bed
# WT 特异，有重叠则移除整条记录
bedtools subtract -A -f 0.3 -a ${GenrichDir}/WT.narrowPeak -b ${GenrichDir}/KO.narrowPeak > ${GenrichDir}/WT_Specific.bed

下面分别是 CommonPeak.bed 输出和原来 KO.narrowPeak 的峰区域数据。除了区域位置信息，其余信息会用 -a 的文件信息。也可以把后面那些列移除掉，免得制造混乱。

# Common
chr1    9957    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180691  181041  peak_2  1000    .       4381.033203     16.855917       -1      330
chr1    199594  199987  peak_3  1000    .       487.511292      5.329252    

# KO 
chr1    9946    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180045  180182  peak_1  1000    .       403.740967      6.181179        -1      56
chr1    180465  181043  peak_2  1000    .       4381.033203     16.855917       -1      330

直接取峰区域的交集差集未免过于简单化处理，另一个方法是用 R 包 DiffBind. DiffBind 是借鉴了 RNA-seq 差异分析。首先将所有样本的峰区域合并，合并方法是用 bedtools merge 所以是所有的样本峰区域的并集。有了这个合并的峰区域，那么可以对每个样本计算里面每个峰区域的 Reads 数，这就像是 RNA-seq 分析每个基因/转录本的 Reads 数一样，然后调用 DESeq2/edgeR 来分析有显著差异的峰区域。
DiffBind 我觉得主要存在两个问题，一是用取所有峰区域并集作为总的峰区域，二是无法确定 RNA-seq 分析软件的假设前提是否在 ATAC-seq 成立，比如说 RNA-seq 是不会大量低丰度基因的，但是 ATAC-seq 是会大量低丰度的峰 Reads.

ATAC-Seq 与 RNA-Seq 相比存在更多低丰度数据

总而言之方法还不算完美，只能自己拿数据去试了，结合课题看是否合适。
DiffBind 教程： 《ATAC-Seq 差异分析 DiffBind》

注释

没有注释就只有峰位置信息，对研究人员来说不太可用。注释可以用 ChIPseeker 包在 R 完成，还可以顺便出一些图。ChIPseeker 的使用参照官方文档。一般来说只看常染色体和性染色体，其余的 contig 可以移除。

awk -v FS="\t" -v OFS="\t" 'length($1)<=5' PeakResult > PeakFilter

同时 ChIPseeker 还可以做通路富集分析。

Motif finding

只依赖于对峰区域的注释是不够直接的，分析转录因子结合的 Motif 能够更直接反映基因转录调控的变化。Motif 分析可以用 HOMMER 或 MEME-Chip 完成。
Motif 分析教程：《ATAC-Seq Motif 富集分析》

TF footprints

另一种转录因子相关的分析叫 TF footprint 分析。转录因子结合到开放染色质导致 Tn5 无法剪切该区域，会在覆盖度高的开放区域（Peak）中出现小范围的覆盖度下降，因此这个下降的小区域称为 TF footprint.

Motif 附近 Reads 分布

可视化

ChIPseeker 有不少总体峰区域分析的可视化图片，这里不再做介绍。下面介绍 2 种生成 Reads 覆盖信号 bw 文件的方法，然后用 pyGenomeTracks 生成特定区域的 Reads 信号图。

手动转换到 bw 文件

警告： 这种方法不适用于建库深度有差异的情况，看看就好。具体解释在《Genrich 的类 bedgraph 文件到 bigwig 文件不合理处》文章。

Genrich -k 参数输出文件是类 Bedgraph 格式的，手动转换到 Bedgraph 就可以用 UCSC genome broswer 浏览了。
我这里实验有重复，所以多个 Bam 文件给 Genrich 同时输入，此时 -k 产生的文件会是 2 个 Bam 文件数据，用一个注释行隔开。所以首先我需要把这 2 数据分开，如果是每个样本单独跑 Genrich 那跳过此步。

# 查看注释行的位置
$ grep -n "#" KO_pileup_p.bed
1:# experimental file: /Example/KO1_ATAC.bam; control file: NA
54091955:# experimental file: /Example/KO2_ATAC.bam; control file: NA

# 分别用 head, tail 取需要的数据
$ head -n 54091954 KO_pileup_p.bed > KO1_pileup_p.bed
$ tail -n +54091955 KO_pileup_p.bed > KO2_pileup_p.bed

# 检查一下行数对不对
$ wc -l KO*pileup_p.bed
   54091954 KO1_pileup_p.bed
   52204177 KO2_pileup_p.bed
  106296131 KO_pileup_p.bed
  212592262 total
  
# WT 样本进行同样处理

分割后将前 2 行即注释行和表头移除，移除后用 cut 命令取前 4 列，输出到 TAB 分隔文件。顺便进行排序。

cut -f 1,2,3,4 KO1_pileup_p.bed | LC_COLLATE=C sort -k1,1 -k2,2n > KO1_Cut.bed

# 其余文件相同处理

然后在文件开头添加一行 "track type=bedGraph", 最后用 bedGraphToBigWig 转换到 bigwig 文件。bedGraphToBigWig 在 Index of /admin/exe/linux.x86_64 下载。

bedGraphToBigWig ${GenrichDir}/${Sample}_Cut.bed ${ChromSize} ${GenrichDir}/${Sample}.bw

其中 ChromSize 用参考基因组和 faToTwoBit, twoBitInfo 两个工具生成，工具也在上面链接下载。

faToTwoBit GRCh38.primary_assembly.genome.fa GRCh38.primary_assembly.2bit
twoBitInfo GRCh38.primary_assembly.2bit GRCh38.primary_assembly.tab

deeptools bamCoverage

deeptools bamCoverage 命令将 Bam 文件转换成 bigwig 文件，免去了自己从 peak calling 结果里各种手动转换文件格式过程。推荐使用这个方法。
这里 --binSize 默认 50 对于 ATAC-seq 我觉得需要设置小一些，--effectiveGenomeSize 的设置值可以在 Effective Genome Size — deepTools 3.4.3 documentation 查看。

bamCoverage --binSize 10 --blackListFileName ${BlackList} --numberOfProcessors 5 --effectiveGenomeSize 2747877777 --normalizeUsing RPGC --outFileFormat bigwig -b ${BamDir}/${Sample}_ShiftSort.bam -o ${BwDir}/${Sample}.bw

其中 RPGC 计算方法：
number of reads per bin / scaling factor for 1x average coverage
Scale factor 计算方法：
(total number of mapped reads * fragment length) / effective genome size

pyGenomeTracks

先从 bw/bed 文件生成 pyGenomeTracks 可用的 tracks 文件。

make_tracks_file --trackFiles ${bwDir}/KO1_ATAC.bw ${bwDir}/KO2_ATAC.bw ${bwDir}/WT1_ATAC.bw ${bwDir}/WT2_ATAC.bw -o ${TrackDir}/Track1.ini

生成的 ini 文件里面包含了画图配置，所以需要修改画图配置的，手动修改这个文件。参数表在 https://pygenometracks.readthedocs.io/en/latest/content/possible-parameters.html 查看。
下面这个配置（不展示重复部分）

[x-axis]

[spacer]
height = 1

[KO1_ATAC]
file = /Example/Vis/KO1_ATAC.bw

title = KO1
height = 2

color = #F56F08
min_value = 0
max_value = 25
number_of_bins = 2000
nans_to_zeros = true
summary_method = mean
show_data_range = true
file_type = bigwig

[gtf]
file = /Example/GRCh38/GENCODE/gencode.v35.annotation.gtf
height = 3
title = Gencode Annotation
file_type = gtf
prefered_name = gene_name
merge_transcripts = true

生成如下图片

一个特定区域 Reads 覆盖

EnrichedHeatmap

EnrichedHeatmap 用于生成一些区域的信号富集图，这个图用 ChIPseeker 和 deeptools 也可以生成，依据自己的喜好选择。
《EnrichedHeatmap 教程》

图来源于 Yin, Senlin 等论文

参考资料

ATAC-Seq data analysis
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Brouilette, Scott, et al. "A simple and novel method for RNA‐seq library preparation of single cell cDNA analysis by hyperactive Tn5 transposase." Developmental Dynamics 241.10 (2012): 1584-1590.
Gontarz, Paul, et al. "Comparison of differential accessibility analysis strategies for ATAC-seq data." Scientific reports 10.1 (2020): 1-13.
Yin, Senlin, et al. "Transcriptomic and open chromatin atlas of high-resolution anatomical regions in the rhesus macaque brain." Nature communications 11.1 (2020): 1-13.
ATAC-seq Guidelines - Harvard FAS Informatics
Buenrostro, Jason D., et al. "ATAC‐seq: a method for assaying chromatin accessibility genome‐wide." Current protocols in molecular biology 109.1 (2015): 21-29.
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Buenrostro, Jason D., et al. "Transposition of native chromatin for multimodal regulatory analysis and personal epigenomics." Nature methods 10.12 (2013): 1213.
Li, Zhijian, et al. "Identification of transcription factor binding sites using ATAC-seq." Genome biology 20.1 (2019): 45.
ATAC-seq 150bp reads