????哈嘍大家好，好久沒有跟新啦，小編最近面臨博士三年級畢業(yè)的壓力，所以公眾號都沒時間打理了，不過偶爾還是能看到有人關(guān)注我的公眾號或者留言，還是很受鼓舞噠~所以還是想抽時間將之前學(xué)的ATAC-seq數(shù)據(jù)分析的代碼分享給大家，以下是正文——
? ? ATAC-seq（Assay for Transposase Accessible Chromatin with high-throughput sequencing）是一種利用轉(zhuǎn)座酶研究染色質(zhì)可接近性的測序技術(shù)，利用DNA轉(zhuǎn)座酶Tn5切割開放的DNA區(qū)域結(jié)合高通量測序研究染色質(zhì)的開放狀態(tài)。與傳統(tǒng)的MNase-seq以及DNase-seq相比，其具有可重復(fù)性強(qiáng)，實(shí)驗步驟簡單，需要的實(shí)驗樣本量少等優(yōu)點(diǎn)，因而被廣泛應(yīng)用。ATAC-seq建庫原理如下圖所示：
前期準(zhǔn)備：
#conda create -n atac python=2創(chuàng)建環(huán)境
#conda activate atac激活環(huán)境
#conda install bedtools deeptools macs2 bowtie bowtie2 samtools
分析流程：
1. 質(zhì)控
采用FASTQC查看測序數(shù)據(jù)質(zhì)量
#fastqc -o FASTQC/ -t 8 Control_R1.fastq.gz Control_R2.fastq.gz Treated_R1.fastq.gz Treated _R2.fastq.gz
#multiqc ./
2. 過濾
采用Cutadapt對測序文件進(jìn)行過濾，目的包括：去除測序引物及接頭、去除reads兩端低質(zhì)量堿基、去除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與參考基因組進(jìn)行比對并采用samtools對BAM文件進(jìn)行排序。其中-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重復(fù)以及線粒體基因，去除低質(zhì)量序列
#sambamba markdup -r Control.bam Control.sambamba.rmdup.bam去除PCR重復(fù)
#samtools view -h -f 2 -q 30 Control.sambamba.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o -> Control.last.bam?去除線粒體重復(fù)，去除低質(zhì)量序列
5. samtools對測序深度、覆蓋度、比對率、重復(fù)率等進(jìn)行統(tǒng)計
#samtools index Control.sambamba.rmdup.bam
#samtools flagstat Control.sambamba.rmdup.bam > Control.rmdup.stat
#cat Control.rmdup.stat?查看stat內(nèi)容
6. 采用bedtools將bam文件轉(zhuǎn)換為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進(jìn)行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數(shù)
8. 采用deeptools及IGV進(jìn)行可視化
對bam文件進(jìn)行歸一化之后可進(jìn)行IGV可視化，同時將bam文件轉(zhuǎn)化為bw文件
#ls *last.bam |while read id; do
> nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.last.bw &
> done |
怎么樣，是不是很簡單呢？快來試一試吧！
