ATAC-seq分析流程（三）

前边介绍了使用DiffBind进行ATAC-seq数据的差异分析，DiffBind进行差异分析时调用了edgeR或者DESeq2包；相对于DiffBind流程我更习惯普通转录组分析流程，因为两者本质上是一样的，区别在于ATAC-seq要根据结果的peak文件自己制备gtf文件，peak对应的gtf文件和gene是不一样的，不同样本间也可能存在差异。

cat *.narrowPeak > merged.peaks
sort -k1,1 -k2,2n merged.peaks > merged.peaks.sorted
bedtools merge -i merged.peaks.sorted > bedtools.merged.sorted.bed
awk -v var=ATAC_seq '{print $1"\t"var"\texon\t"$2"\t"$3"\t.\t+\t.\tgene_id \"peak_"NR"\"; transcript_id \"peak_"NR"\";"}' bedtools.merged.sorted.bed > bedtools.merged.atac.gtf
featureCounts -T 18 -p -t exon -g gene_id -a bedtools.merged.atac.gtf -o ATAC_counts.txt *last.bam
sed 1d ATAC_counts.txt | cut -f 1,7-15 > rawATAC_counts.matrix

一个peak相当于一个gene，得到的原始矩阵可以使用DESeq2包分析，筛选出共有特有或者差异的peaks等进行peaks的注释。