????哈喽大家好,好久没有跟新啦,小编最近面临博士三年级毕业的压力,所以公众号都没时间打理了,不过偶尔还是能看到有人关注我的公众号或者留言,还是很受鼓舞哒~所以还是想抽时间将之前学的ATAC-seq数据分析的代码分享给大家,以下是正文——
? ? ATAC-seq(Assay for Transposase Accessible Chromatin with high-throughput sequencing)是一种利用转座酶研究染色质可接近性的测序技术,利用DNA转座酶Tn5切割开放的DNA区域结合高通量测序研究染色质的开放状态。与传统的MNase-seq以及DNase-seq相比,其具有可重复性强,实验步骤简单,需要的实验样本量少等优点,因而被广泛应用。ATAC-seq建库原理如下图所示:
前期准备:
#conda create -n atac python=2创建环境
#conda activate atac激活环境
#conda install bedtools deeptools macs2 bowtie bowtie2 samtools
分析流程:
1. 质控
采用FASTQC查看测序数据质量
#fastqc -o FASTQC/ -t 8 Control_R1.fastq.gz Control_R2.fastq.gz Treated_R1.fastq.gz Treated _R2.fastq.gz
#multiqc ./
2. 过滤
采用Cutadapt对测序文件进行过滤,目的包括:去除测序引物及接头、去除reads两端低质量碱基、去除N碱基过多的reads、去除截短后单端reads长度小于75bp的reads。
#cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -q 30 -m 75 --trim-n --report=minimal -o Control_out_R1.fastq.gz -p Control_out_R2.fastq.gz Control_R1.fastq.gz Control_R2.fastq.gz
3.比对及排序
使用Bowtie2将clean data与参考基因组进行比对并采用samtools对BAM文件进行排序。其中-X代表最大插入片段,宽泛的插入范围为10-1000bp,一般一个核小体147 180,大片段很稀有。
#bowtie2 -p 10 -X 1000 -x bowtie_index -1 Control_out_R1.fastq -2 Control_out_R2.fastq |samtools sort -O bam -@ 5 -o Control.bam
4. 采用sambamba及samtools去除PCR重复以及线粒体基因,去除低质量序列
#sambamba markdup -r Control.bam Control.sambamba.rmdup.bam去除PCR重复
#samtools view -h -f 2 -q 30 Control.sambamba.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o -> Control.last.bam?去除线粒体重复,去除低质量序列
5. samtools对测序深度、覆盖度、比对率、重复率等进行统计
#samtools index Control.sambamba.rmdup.bam
#samtools flagstat Control.sambamba.rmdup.bam > Control.rmdup.stat
#cat Control.rmdup.stat?查看stat内容
6. 采用bedtools将bam文件转换为bed文件
#bedtools bamtobed -i Control.last.bam >Control.last.bed
#less Control.last.bed?查看bed文件
#bedtools intersect -a ../peaks/Control.bed -b Control_summits.bed?计算插入片段长度
7.采用macs2进行call peaks
# macs2 callpeak -t Control.last.bed -g hs --nomodel --shift -100 --extsize 200 -n Control --outdir ../peaks
# wc Control_peaks.narrowPeak查看样本peaks
# wc -l *bed?查看所有样本的peaks数
8. 采用deeptools及IGV进行可视化
对bam文件进行归一化之后可进行IGV可视化,同时将bam文件转化为bw文件
#ls *last.bam |while read id; do
> nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.last.bw &
> done |
怎么样,是不是很简单呢?快来试一试吧!