大家好这里是有时候季更，有时候月更的基因组学研究生。上一期记录了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.bedcat 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命名为countcount$std <- apply(count,1,sd) # 计算标准差count.variable <- count[which(count$std>1),] #直接截取标准差大于1的Regionpheatmap(count.variable[,1:6],cluster_cols = F,cluster_rows = T,show_rownames = F,scale = "row")# 可以根据 pheatmap 参数调整图片，这里选取前6列作图差异分析可以用DESeq2包，注释可以用clusterProfiler, Chipseeker包。今天就记录到这里吧，欢迎后台回复和交流~
欢迎关注基因组学研究生
