哈喽大家好，好久没有跟新啦，小编最近面临博士三年级毕业的压力，所以公众号都没时间打理了，不过偶尔还是能看到有人关注我的公众号或者留言，还是很受鼓舞哒~所以还是想抽时间将之前学的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 |
怎么样，是不是很简单呢？快来试一试吧！
网友评论