ATAC-seq(3) -- 比对

bowtie2 比对

其中bowtie2比对加入了-X 2000 参数，是最大插入片段，宽泛的插入片段范围(10-1000bp)

cat >bowtie2.sh
index=~/my_data/ref_data/mouse/bowtie2_index/mm10
ls *.fq.gz | while read id;
do 
file=$(basename $id)
sample=${file%%_*}
bowtie2  -p 5  --very-sensitive -X 2000 -x  $index -1 ${sample}_1_val_1.fq.gz -2 ${sample}_2_val_2.fq.gz | samtools sort  -O bam  -@ 5 -o - > ../bowtie2/${sample}.bam 
samtools index  ../bowtie2/${sample}.bam 
# 比对结束后，需要了解比对结果的情况
samtools flagstat  ../bowtie2/${sample}.bam  >  ../bowtie2/${sample}.stat
#记录mapping的结果，其中第一列是染色体名称，第二列是序列长度，第三列是mapped reads数，第四列是unmapped reads数
samtools idxstats ../bowtie2/${sample}.bam >  ../bowtie2/${sample}_idxstats.txt
done
nohup sh bowtie2.sh &

计算线粒体 Reads 占比

ls *.txt | while read id
do
sample=${id%.*}
mtReads=$(cat  ${sample}.txt | grep 'chrM' | cut -f 3)
totalReads=$(cat  ${sample}.txt | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
done
==> mtDNA Content: 24.10%
==> mtDNA Content: 23.40%
==> mtDNA Content: 17.64%
==> mtDNA Content: 22.70%

去除比对到线粒体上的reads

mkdir ../no_chrM_bam
ls *.bam | while read id
do 
sample=${id%.*}
samtools view  -h  -f 2 -q 30   ${sample}.bam   | grep -v chrM |samtools sort  -O bam  -@ 5 -o - > ../no_chrM_bam/${sample}.bam
samtools index   ../no_chrM_bam/${sample}.bam 
samtools flagstat  ../no_chrM_bam/${sample}.bam > ../no_chrM_bam/${sample}.stat 
done 

计算插入片段长度

检查比对报告，可能 "aligned concordantly >1 times" 比例会比较高，原因一是 --very-sensitive 参数，让比对质量稍低于最佳比对的仍保留，第二是 ATACseq 数据线粒体含量高，线粒体有许多区域序列类似。

前面说了 ATACseq 插入片段应该在约 <100,200,400,600bp 大小有峰，用 gatk CollectInsertSizeMetrics 命令分析比对后插入片段分布。
安装gatk

conda install -c bioconda/label/cf201901 gatk
或者
nohup wget https://github.com/broadinstitute/gatk/releases/download/4.1.3.0/gatk-4.1.3.0.zip &
unzip gatk-4.1.3.0.zip
echo export PATH="/home/chenyuqiao/my_disk/biosoft/gatk-4.1.3.0:$PATH" >> ~/.bashrc
source .bashrc
gatk --help

用 gatk CollectInsertSizeMetrics 命令分析比对后插入片段分布

ls *.bam | while read id;
do
sample=${id%.*}
gatk CollectInsertSizeMetrics -H ../gatk/${sample}_InsertSize.pdf -I ${sample}.bam -O ../gatk/${sample}_InsertSize.txt
done

image.png

（可选步骤）移动 Reads

如下面示意图所示，Tn5 建库时会插入 9bp 的序列，所以比对到正链（+）的 reads 移动 +4 bp, 比对负链（-）的 reads 移动 -5 bp.

image.png

这个移动可以用 deeptools 的 alignmentSieve 命令完成。

alignmentSieve --numberOfProcessors 8 --ATACshift -b \
${bam_dir}/${sample}_nSorted.bam -o ${bam_dir}/${sample}_Shift.bam

一般的分析这点 reads 位移没影响，不需要进行这个步骤。