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的注释。

发表评论

邮箱地址不会被公开。 必填项已用*标注