ATAC-seq数据可视化——超简单教程

    大家好这里是有时候季更，有时候月更的基因组学研究生。上一期记录了ATAC数据从Fastq数据到bam文件的处理脚本。

ATAC数据一键处理

基因组学研究生，公众号：基因组学研究生一个不成熟的小脚本，ATAC数据一键预处理从fastq到treated bam

    后来看到有同学后台问我有没有ATAC-seq可视化教程，所以简单写了写，希望能帮到部分人~ 

Ps：那位问我的同学，由于我消息晚了几天看到，所以回复不了了，希望这篇文章你能看到。

ATAC数据可视化

一、bam to peak    

假设我们已经得到了bam文件，接下来首先使用MACS2计算的peak，代码非常简单，如下：

macs2 callpeak -t bamfile1 bamfile2... -n output_file -g mm -f BAM --nomodel  --shift -100 --extsize 200

二、Make Atlas 

    工具：bedtools

    为了使得数据之间具有可比性，我们需要制作一个矩阵，类似于基因表达矩阵(行为Gene,列为样本)。ATAC矩阵则是行为Genome Region（即peak），列为样本。但是由于每个样本的开放Region不一样，所以我们需要制作一个包含所有Region的Region Atlas。

    MACS2将会输出3个文件，其中包括一个后缀为“_summits.bed”的文件，即peak的中心位置，我们使用这个文件制作总的Peak Atlas。

cat A_summits.bed|awk 'BGEIN{OFS="\t"}{print $1,$2-250,$3+250}' > A_summits_250bp.bed

cat A_summits_250bp.bed B_summits_250bp.bed... | awk '{if($2<0){print $1 "\t" 0 "\t" $3} else {print $0}}' |sortBed -i - |bedtools merge -i - > atlas

三、计算各个样本的Signal

    工具：deeptools

    获得Atlas后，我们来计算每个样本在每个Peak内的信号值。

multiBamSummary BED-file --BED ${atlas_file} --bamfiles ${bamfile}*sort.bam -o ${out}all_count.npz --outRawCounts ${out}all_count.txt

    这样我们就得到了一个Genome Accessibility Matrix

四、可视化

    工具：R

    有了matrix，就可以做很多事情了。可以做差异分析，peak注释，与RNA expr联合分析等。

    ATAC热图就很好看，先简单做个热图：

library(pheatmap) # 原始数据是有很多peak的，peak太多作图很慢，我们直接去掉一些没有变化的Region，# 假设你已经read table命名为count

count$std <- apply(count,1,sd) # 计算标准差

count.variable <- count[which(count$std>1),] #直接截取标准差大于1的Region

pheatmap(count.variable[,1:6],cluster_cols = F,cluster_rows = T,show_rownames = F,scale = "row")# 可以根据 pheatmap 参数调整图片，这里选取前6列作图

    差异分析可以用DESeq2包，注释可以用clusterProfiler, Chipseeker包。今天就记录到这里吧，欢迎后台回复和交流~

公众号原文

欢迎关注基因组学研究生