JiaoYuan's Blog

ATAC-Seq 分析流程

2024-05-03

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

ATAC-seq 概述

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

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

版本信息

一些名词解释

软件安装

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

先创建 conda 虚拟环境，安装所需要的软件，可以自行手动安装，也可以直接导入我的 conda 环境：

git clone https://github.com/imjiaoyuan/NGS-analysis.git
cd NGS-analysis/environment
cp .condarc ~/
conda env create --file atac-seq-env.yml   # ATAC-seq

创建完成后激活环境就可以使用了：

conda activate atac-seq

数据获取与预处理

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

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

使用 sratools 进行批量下载，下载比较慢，就用 nohup 挂在后台下载，自己可以做点别的事：

nohup prefetch -O . $(<SRR_Acc_List.txt) &

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

mkdir 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-3 ./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

进行质量检测

直接进行批处理：

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

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

nohup fastqc SRR*.fastq.gz &
multiqc ./fastqc_report # 将多个质量检测报告合并

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

Trim Galore 过滤低质量序列

常用参数

trim_galore [options] <filename(s)>

过滤低质量序列

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

mkdir 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 Oryza_sativa

批量进行比对：

cat ./SRR_Acc_List.txt | while read id;
do
    bwa mem -v 3 -t 4 ./Oryza_sativa/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa ./clean/${id}_1.fastq.gz ./clean/${id}_2.fastq.gz -o ./compared/${id}.sam
done

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

samtools 排序压缩

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

cat ./SRR_Acc_List.txt | while read id ;do
    samtools \
    sort \
    -n -@ 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 [-h] [--version] {callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak} ...

Peak Calling

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

cat ./SRR_Acc_List.txt | while read id;do
    macs2 \
    callpeak \
    -t ./sorted/${id}.bam \
    -f BAM \
    -g 3.6e+8 \
    -n ./peaks/${id} \
done

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

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

deeptools 包含了很多的工具：

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

bamCoverage

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

bamCoverage -b reads.bam -o coverage.bw

必须的参数：

标准化参数：

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

此处我使用的命令是：

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 [-h] [--version]  ...

必须的参数：

输出参数：

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

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

绘制热图与折线图：

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

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

computeMatrix reference-point \
--referencePoint TSS \
--missingDataAsZero \
-b 3000 -a 3000 \
-R SRR28122456.bed \
-S SRR28122455.bw SRR28122456.bw \
-o TSS.gz

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

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

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

Chipseek 可视化 peak

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

安装 ChIPseeker 包：

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

BiocManager::install("ChIPseeker")

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

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

BiocManager::install("GenomicFeatures")
BiocManager::install("txdbmaker")

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

library("GenomicFeatures")

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

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

# 为了这里方便分析，将两个 bed 文件按照样本信息进行重命名了
Mock <- readPeakFile('./peaks/Mock.bed')
RSV <- readPeakFile('./peaks/RSV.bed')

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

plotAnnoPie(Mock_peakAnno)
plotAnnoPie(RSV_peakAnno)

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

peaks <- list(Mock_peakAnno=Mock, RSV_peakAnno=RSV)
peakAnno <- lapply(peaks, annotatePeak, tssRegion = c(-3000, 3000), TxDb = spombe)
plotAnnoBar(peakAnno)

在染色体上可视化 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 输出文件中的 bed 文件，只保留染色体、开始位置、终止位置和值四列即可，如：

1	1074	1075	33.7137
1	84253	84254	2.21489
1	104155	104156	5.22729
1	104518	104519	3.93076
1	104876	104877	6.8855

将文件导入到 TBtools：

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

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

参考资料

This site is open source. Improve this page.