ChIP-seq数据分析实战训练（二）

四、过滤

raw data -> fastqc

###过滤前先看一下fastqc质检报告确定是否需要过滤以及参数设置
cd /public/workspace/fangwen/learn/chip-seq/qc/
ls ../srr/raw/*gz|xargs fastqc -t 10 -o  ./raw/

使用trim_galore软件进行质控，过滤

###单端测序
cd /public/workspace/fangwen/learn/chip-seq/clean
analysis_dir=/public/workspace/fangwen/learn/chip-seq/
bin_trim_galore="trim_galore"    #可不要
ls ../raw/*gz | while read fq1;
do 
nohup $bin_trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $analysis_dir/clean  $fq1 & 
done

过滤后的文件如下：

(chip) [fangwen@server clean]$ ls -lh *.gz
-rw-rw-r-- 1 fangwen fangwen 325M Jul 14 20:43 Control_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 405M Jul 14 20:44 Control_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 419M Jul 14 20:44 H2Aub1_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 638M Jul 14 20:45 H2Aub1_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 723M Jul 14 20:46 H3K36me3_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 647M Jul 14 20:46 H3K36me3_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 200M Jul 14 20:43 Ring1B_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 200M Jul 14 20:43 Ring1B_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 719M Jul 14 20:47 Ring1B_3_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 368M Jul 14 20:44 RNAPII_8WG16_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 633M Jul 14 20:46 RNAPII_8WG16_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 716M Jul 14 20:47 RNAPII_S2P_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 258M Jul 14 20:43 RNAPII_S2P_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 420M Jul 14 20:44 RNAPII_S2P_3_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 553M Jul 14 20:45 RNAPII_S5P_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 573M Jul 14 20:45 RNAPII_S5P_2_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 190M Jul 14 20:42 RNAPII_S5PRepeat_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 373M Jul 14 20:44 RNAPII_S7P_1_trimmed.fq.gz
-rw-rw-r-- 1 fangwen fangwen 279M Jul 14 20:43 RNAPII_S7P_2_trimmed.fq.gz

clean data -> fastqc

cd /public/workspace/fangwen/learn/chip-seq/qc/
ls ../clean/*gz|xargs fastqc -t 10 -o  ./clean/
multiqc ./

五、进行比对（bowtie2）

参考基因组及注释文件下载

bowtie2官方manual

http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml

###直接下载小鼠索引文件mm10
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
###单端测序数据的比对：
cd /public/workspace/fangwen/learn/chip-seq/align/
#一定要搞清楚自己的索引文件bowtie2_index在什么地方！！！
ls ../clean/*gz |while read id;
do 
file=$(basename $id )
sample=${file%%.*}
#可选择性的将trimmed去掉
echo $file $sample
## 比对过程3分钟一个样本
bowtie2  -p 5  -x  /public/workspace/fangwen/learn/chip-seq/referece/mm10 -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done 

Notes:

1.现有比对工具在做mapping之前，都需要下载对应物种的参考基因组做index，而如何选择合适的参考基因组是一件非常重要的事情。
2.用bowtie2进行比对和统计比对率,需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件
现有的参考基因组存储网站三个：
ENSEMBL
UCSC
NCBI

对bam文件进行QC

cd /public/workspace/fangwen/learn/chip-seq/align/
ls  *.bam  |xargs -i samtools index {} 
ls  *.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done

qc后的结果如下：

(chip) [fangwen@server align]$ ls -lh *.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 Control_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 Control_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 H2Aub1_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 H2Aub1_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 H3K36me3_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 H3K36me3_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 Ring1B_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 Ring1B_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 Ring1B_3_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_8WG16_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 RNAPII_8WG16_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 RNAPII_S2P_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_S2P_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_S2P_3_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 RNAPII_S5P_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 388 Jul 15 13:40 RNAPII_S5P_2_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_S5PRepeat_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_S7P_1_trimmed.stat
-rw-rw-r-- 1 fangwen fangwen 386 Jul 15 13:40 RNAPII_S7P_2_trimmed.stat
(chip) [fangwen@server align]$ cat Control_1_trimmed.stat
8449874 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7438540 + 0 mapped (88.03% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

参考资料：

https://www.jianshu.com/p/f84ffba2ec1e
https://www.jianshu.com/p/1384173c353b