仓库源文，站点原文

ATAC-Seq 技术是"Assay for Transposase-Accessible Chromatin with high-throughput Sequencing"的缩写。 ATAC-Seq 方法依赖于使用高活性转座酶 Tn5 的下一代测序（NGS）文库的构建。将 NGS 接头连接到转座酶上，该转座酶可以使染色质断裂并同时将这些接头整合到开放的染色质区域中。构建的文库通过 NGS 测序，并使用生物信息学分析具有可及或可访问染色质的基因组区域。

<!--more-->

ATAC-seq 概述

虽然 ATAC-seq 实验方法操作相对简单且结果稳定，但是专门为 ATAC-seq 测序数据开发的生物信息学分析软件却非常少。由于 ATAC-seq 和 ChIP-seq 数据的相似性较高，ChIP-seq 分析使用的软件一般也可用于 ATAC-seq 的分析，但是使用 ChIP-seq 软件分析得到的 ATAC-seq 结果尚未得到系统性的评估。

上游的预处理部分与 RNA-seq 的差不多，可以参考我之前的文章：https://yuanj.top/p/rna-seq-analysis-learning/

版本信息

一些名词解释

软件安装

只需要使用 conda 就可以安装所有需要的软件，主要使用的软件有以下一些：

先创建 conda 虚拟环境，安装所需要的软件。

数据获取与预处理

从 NCBI SRA 数据库下载 SRR_Acc_List.txt 文件：

后面的分析都可以基于这个文件进行批处理。

使用 prefetch 工具批量下载测序数据，下载比较慢，就用 nohup 挂在后台下载，自己可以做点别的事：

#!/bin/bash
nohup prefetch --option-file SRR_Acc_List.txt &

下载后的数据是放在以 SRR 号命名的文件夹里面的，我们用一个小脚本把它全部移动到一个 sra 文件夹便于处理：

mkdir -p sra
cat ./SRR_Acc_List.txt | while read id; do
    mv -f "${id}/${id}.sra" ./sra/
    rm -rf "${id}"
done

刚刚下载好的数据是 sra 格式的，使用 sratools 将其拆分并输出到 fastqgz 文件夹：

mkdir fastqgz
nohup fastq-dump --gzip --split-files ./sra/*.sra --outdir ./fastqgz &

参考基因组及注释文件

植物的我一般在 Ensembl Plants 下载，用 wget 或 curl 都可以，内存不大，这组数据是水稻的，所以就从 Ensembl Plants 下载水稻的基因组和注释文件：

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-57/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.57.gff3.gz

然后解压：

gzip -d Oryza_sativa.IRGSP-1.0.57.gff3.gz
gzip -d Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

FastQc 质控

常用参数

fastqc [-o output dir] [–(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

进行质量检测

直接进行批处理：

mkdir -p fastqc_report
fastqc ./fastqgz/*.fastq.gz -o ./fastqc_report

当然，数据比较多的时候还是挂在后台批处理然后等着就行：

nohup fastqc ./fastqgz/*.fastq.gz -o ./fastqc_report &
# 等待 fastqc 运行完毕
multiqc ./fastqc_report # 将多个质量检测报告合并

质量检测报告查看和解读之前的文章有，这里就不赘述了。

Trim Galore 过滤低质量序列

常用参数

trim_galore [options] <filename(s)>

过滤低质量序列

使用一个批处理对所有数据进行处理：

mkdir -p clean
cat ./SRR_Acc_List.txt | while read id; do 
    trim_galore \
    -q 20 \
    --length 36 \
    --max_n 3 \
    --stringency 3 \
    --fastqc \
    --paired \
    -o ./clean/ \
    ./fastqgz/${id}_1.fastq.gz ./fastqgz/${id}_2.fastq.gz
done

完成后文件会输出到 clean 文件夹。

bwa 比对到参考基因组

常用参数

bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>

bwa mem [options] ref.fa reads.fq [mates.fq]

进行比对

使用前面下载的水稻基因组建立索引：

bwa index -a bwtsw Oryza_sativa.IRGSP-1.0.dna.toplevel.fa

批量进行比对：

mkdir -p compared
cat ./SRR_Acc_List.txt | while read id; do
    bwa mem \
        -v 3 \
        -t 4 \
        Oryza_sativa.IRGSP-1.0.dna.toplevel.fa \
        ./clean/${id}_1_val_1.fq.gz \
        ./clean/${id}_2_val_2.fq.gz \
        -o ./compared/${id}.sam
done

比对后的结果会输出到 compared 文件夹。

samtools 排序压缩

这个就很常见了，几乎是组学分析必用的软件，之前也详细介绍过了，不再赘述，使用下列批处理将所有数据按坐标排序、压缩并输出到 sorted 文件夹：

mkdir -p sorted
cat ./SRR_Acc_List.txt | while read id ;do
    samtools \
    sort \
    -@ 5 \
    ./compared/${id}.sam \
    -o ./sorted/${id}.bam
done

接着需要给每个 bam 文件创建索引，后面会用到：

cat ./SRR_Acc_List.txt | while read id;
do
    samtools index \
    -@ 4 \
    ./sorted/${id}.bam
done

macs2 进行 Peak Calling

这是非常关键的一步，ATAC-seq 分析的核心结果都要从这里出，使用的是 macs2。

常用参数

macs2 callpeak -t treatment.bam -c control.bam -f BAM -g "effective genome size" -n "output file name"

Peak Calling

水稻基因组大小是 3.7e+8，这个基因组的大小是很关键的，如果指定错误或者未指定都会影响后面的分析。使用下列批处理进行分析：

mkdir -p peaks
cat ./SRR_Acc_List.txt | while read id;do
    macs2 callpeak \
    -t ./sorted/${id}.bam \
    -f BAMPE \
    --nomodel \
    -g 3.7e+8 \
    -n ${id} \
    --outdir ./peaks/
done

这里如果不知道基因组的大小的话可以通过faCount来计算：

conda install -c bioconda ucsc-facount
faCount /path/to/genome.fa > size.txt

完成后每个样本会输出几个文件：

deeptools 计算基因组区段 reads 覆盖度

deeptools 包含了很多的工具：

这里我们主要使用的是 bamCoverage 和 computeMatrix。

bamCoverage

bamCoverage 利用测序数据比对结果转换为基因组区域 reads 覆盖度结果。可以自行设定覆盖度计算的窗口大小 (bin)。

bamCoverage -b reads.bam -o coverage.bw

必须的参数：

标准化参数：

可选的参数很多，可以参考 deeptools 的手册。

此处我使用的命令是：

mkdir -p bigwig
cat ./SRR_Acc_List.txt | while read id;do
    bamCoverage \
    -b ./sorted/${id}.bam \
    -o ./bigwig/${id}.bw \
    --normalizeUsing RPKM
done

这里就需要我们前面给 bam 文件创建的 samtools 索引，不然的话会报错。完成之后会输出一系列 bw 文件，这个文件可以导入到 IGV 进行可视化：

computeMatrix

computeMatrix tool_mode [-h] [--version]  ...

必须的参数：

输出参数：

此处我对单个文件进行批处理计算：

cat ./SRR_Acc_List.txt | while read id;do
    computeMatrix reference-point \
    --referencePoint TSS \
    --missingDataAsZero \
    -b 3000 -a 3000 \
    -R ./peaks/${id}_peaks.narrowPeak \
    -S ./bigwig/${id}.bw  \
    -o ./deeptools_matrix/${id}.TSS.gz
done

绘制热图与折线图：

cat ./SRR_Acc_list.txt | while read id;do
    plotHeatmap -m ./deeptools_matrix/${id}.TSS.gz -out ./deeptools_matrix/${id}.TSS.png
done

再同时对多个样本进行计算，将输出结果放在一个文件：

computeMatrix reference-point \
--referencePoint TSS \
--missingDataAsZero \
-b 3000 -a 3000 \
-R ./peaks/SRR28122456_peaks.narrowPeak \
-S ./bigwig/SRR28122455.bw ./bigwig/SRR28122456.bw \
-o ./deeptools_matrix/TSS.gz

计算完成之后绘制多个样本的折线图在同一个图内：

plotProfile --dpi 720 -m ./deeptools_matrix/TSS.gz -out ./deeptools_matrix/TSS.pdf --plotFileFormat pdf --perGroup

这个图左边是前面绘制的两个热图+折线图，右边是前面绘制多个样本折线图在同一图内。

ChIPseeker 可视化 peak

R 包 ChIPseeker 主要是为了能对 ChIP-seq 数据进行注释与可视化，主要对 peak 位置及 peak 邻近基因的注释。然而，在之后对 ChIPseeker 的应用中，发现它不局限于 ChIP-seq，可用于其他的 peak（如 ATAC-seq，DNase-seq 等富集得到的）注释。该包功能强大之处还是在于可视化上。

安装 ChIPseeker 包：

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ChIPseeker")

所需的文件有两个，一个是 BED 格式的文件，至少得有染色体名字、染色体起始位点和染色体终止位点，其它信息如 name，score，strand 等可有可无。这里直接用前面 macs2 输出的 _peaks.narrowPeak 文件即可。另一个是有注释信息的 TxDb 对象，Bioconductor 包提供了 30 个 TxDb 包，包含了很多物种，如人，老鼠等。当所研究的物种没有已有的 TxDb 时，可通过 GenomicFeatures 包使用基因组注释文件进行制作：

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GenomicFeatures")

这里我使用前面下载的水稻注释文件进行制作：

library("GenomicFeatures")

txdb <- makeTxDbFromGFF("Oryza_sativa.IRGSP-1.0.57.gff3")

然后就可以读入 peak 文件进行可视化：

library("ChIPseeker")
# 为了这里方便分析，将两个 narrowPeak 文件按照样本信息进行重命名了
# 假设已将 SRR28122455_peaks.narrowPeak 和 SRR28122456_peaks.narrowPeak 重命名
# 为 Mock_peaks.narrowPeak 和 RSV_peaks.narrowPeak
Mock <- readPeakFile('./peaks/Mock_peaks.narrowPeak')
RSV <- readPeakFile('./peaks/RSV_peaks.narrowPeak')

Mock_peakAnno <- annotatePeak(Mock, tssRegion =c(-3000, 3000), TxDb = txdb)
RSV_peakAnno <- annotatePeak(RSV, tssRegion =c(-3000, 3000), TxDb = txdb)

plotAnnoPie(Mock_peakAnno)
plotAnnoPie(RSV_peakAnno)

也可以将多个样本绘制在同一个图：

peaks <- list(Mock=Mock, RSV=RSV)
peakAnnoList <- lapply(peaks, annotatePeak, tssRegion = c(-3000, 3000), TxDb = txdb, annoDb = "org.Oryza.sativa.eg.db") # 如果有对应的物种注释包，可以加上
plotAnnoBar(peakAnnoList)

在染色体上可视化 peaks 的位置

这个使用 TBtools 来进行，先制作一个染色体长度文件：

1   43270923
2   35937250
3   36413819
4   35502694
5   29958434
6   31248787
7   29697621
8   28443022
9   23012720
10  23207287
11  29021106
12  27531856

然后使用 TBtools 的 circos 绘图工具，只把染色体长度文件导入：

修改 macs2 输出文件中的 narrowPeak 文件，只保留染色体、开始位置、终止位置和 q-value/p-value 值四列即可，如：

1   1074    1075    3.63838
1   84253   84254   0.48511
1   104155  104156  0.78204
1   104518  104519  0.65582
1   104876  104877  0.92383

将文件导入到 TBtools：

然后改一下颜色、字体什么的进行一下简单的美化：

到这里大致的流程基本上就完成了，还有一些分析就是视情况才做的了。

参考资料