????大家好这里是有时候季更，有时候月更的基因组学研究生。上一期记录了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，就可以做很多事情了?？梢宰霾钜旆治觯琾eak注释，与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包。今天就记录到这里吧，欢迎后台回复和交流~
欢迎关注基因组学研究生
