【ChIP-seq 实战】四、fastq数据的质量控制和比对

作者: 佳奥 | 来源:发表于2022-08-10 17:03 被阅读0次

【ChIP-seq 实战】四、fastq数据的质量控制和比对
ChIP-Seq分析
RNA-Seq数据分析—fastp v0.23.1
Biostar handbook学习笔记五-基因组测序技术原理简
RNA-Seq数据分析——原始数据质量控制（QC）
ChIP-seq数据分析实战训练（二）
生信文件格式fastqc
临床生信实操（centos7）1
RNA-Seq数据分析—cutadapt质量控制
RNA-seq分析：从fastq到差异表达基因

这里是佳奥，继续处理数据吧。

1 质量控制（trim_galore）

##本次数据是单端测序数据

##进入到raw目录
(chipseq) root 21:24:01 /home/kaoku/chipseq/mouse_project/raw
$ ls
Control_1.fastq.gz  H3K36me3_1.fastq.gz  Ring1B_3.fastq.gz        RNAPII_S2P_2.fastq.gz  RNAPII_S5PRepeat_1.fastq.gz
Control_2.fastq.gz  H3K36me3_2.fastq.gz  RNAPII_8WG16_1.fastq.gz  RNAPII_S2P_3.fastq.gz  RNAPII_S7P_1.fastq.gz
H2Aub1_1.fastq.gz   Ring1B_1.fastq.gz    RNAPII_8WG16_2.fastq.gz  RNAPII_S5P_1.fastq.gz  RNAPII_S7P_2.fastq.gz
H2Aub1_2.fastq.gz   Ring1B_2.fastq.gz    RNAPII_S2P_1.fastq.gz    RNAPII_S5P_2.fastq.gz

##首先进行fastqc，进入qc目录，qc文件在raw目录，结果会输入到qc目录
ls ../raw/*gz|xargs fastqc -t 10 -o  ./

(chipseq) root 21:43:49 /home/kaoku/chipseq/mouse_project/qc
$ ls
Control_1_fastqc.html  H3K36me3_1_fastqc.html  Ring1B_3_fastqc.html        RNAPII_S2P_2_fastqc.html  RNAPII_S5PRepeat_1_fastqc.html
Control_1_fastqc.zip   H3K36me3_1_fastqc.zip   Ring1B_3_fastqc.zip         RNAPII_S2P_2_fastqc.zip   RNAPII_S5PRepeat_1_fastqc.zip
Control_2_fastqc.html  H3K36me3_2_fastqc.html  RNAPII_8WG16_1_fastqc.html  RNAPII_S2P_3_fastqc.html  RNAPII_S7P_1_fastqc.html
Control_2_fastqc.zip   H3K36me3_2_fastqc.zip   RNAPII_8WG16_1_fastqc.zip   RNAPII_S2P_3_fastqc.zip   RNAPII_S7P_1_fastqc.zip
H2Aub1_1_fastqc.html   Ring1B_1_fastqc.html    RNAPII_8WG16_2_fastqc.html  RNAPII_S5P_1_fastqc.html  RNAPII_S7P_2_fastqc.html
H2Aub1_1_fastqc.zip    Ring1B_1_fastqc.zip     RNAPII_8WG16_2_fastqc.zip   RNAPII_S5P_1_fastqc.zip   RNAPII_S7P_2_fastqc.zip
H2Aub1_2_fastqc.html   Ring1B_2_fastqc.html    RNAPII_S2P_1_fastqc.html    RNAPII_S5P_2_fastqc.html
H2Aub1_2_fastqc.zip    Ring1B_2_fastqc.zip     RNAPII_S2P_1_fastqc.zip     RNAPII_S5P_2_fastqc.zip

##multiqc看一下
(chipseq) root 22:03:24 /home/kaoku/chipseq/mouse_project/qc
$ multiqc ./

有点窒息，这个结果

QQ截图20220809220656.png

所以需要进一步过滤。

##选择trim_galore软件进行过滤
analysis_dir=/home/kaoku/chipseq/mouse_project
ls /home/kaoku/chipseq/mouse_project/raw/*.gz | while read fq1;
do 
trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $analysis_dir/clean  $fq1 & 
done

##过滤结果
Total reads processed:               5,861,669
Reads with adapters:                    80,435 (1.4%)
Reads written (passing filters):     5,861,669 (100.0%)

Total basepairs processed:   211,020,084 bp
Quality-trimmed:              25,712,707 bp (12.2%)
Total written (filtered):    183,863,864 bp (87.1%)

##把原本qc目录清空，把过滤后的数据再次进行qc
ls ../clean/*gz| xargs fastqc -t 10 -o  ./
multiqc ./

QQ截图20220809223626.png

发现过滤效果不错，可以进行下一步的比对了。

2 建立索引（bowtie2）

用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件：

下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10

2.1 官方索引

##索引大小为3.2GB， 可以直接下载索引文件，代码如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip

##这个方法，简便，省事，但是我下载不下来啊orz
mm10.zip        3%[=>                                                       ] 104.95M  33.0KB/s    剩余 19h 53m

2.2 自建索引

##使用浏览器下载或者wget
wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz 
tar -zxvf chromFa.tar.gz 
cat *.fa > mm10.fa

##构建索引
bowtie2-build mm10.fa mm10

##这个方法很好，但是8G内存不够，提示out of memory。。。orz

2.3 遇事不决浏览器

##如果wget下载过慢，自建索引内存不足的话，回到官网
http://ccb.jhu.edu/software.shtml

##选择Software&Data，bowtie进入网页，选择侧边的bowtie2，可以查看到对应软件的索引文件
##复制该链接到浏览器或wget即可下载
https://genome-idx.s3.amazonaws.com/bt/mm10.zip

(chipseq) root 13:47:52 /home/kaoku/reference/mouse
$ ls -lh
总用量 6.7G
-rwxr-xr-x 1 root  root  2.7K  5月  3  2012 make_mm10.sh
-rw-r--r-- 1 root  root  848M  5月  2  2012 mm10.1.bt2
-rw-r--r-- 1 root  root  633M  5月  2  2012 mm10.2.bt2
-rw-r--r-- 1 root  root  6.0K  5月  2  2012 mm10.3.bt2
-rw-r--r-- 1 root  root  633M  5月  2  2012 mm10.4.bt2
-rw-r--r-- 1 root  root  848M  5月  3  2012 mm10.rev.1.bt2
-rw-r--r-- 1 root  root  633M  5月  3  2012 mm10.rev.2.bt2
-rw-r--r-- 1 kaoku kaoku 3.2G  8月 10 13:43 mm10.zip

3 序列比对（bowtie2）

##明确软件安装的的位置和参考文件的位置

##从过滤后clean目录的文件开始，进入align目录
(chipseq) root 14:50:51 /home/kaoku/chipseq/mouse_project
$ cd align/

##比对流程
bowtie2_index=/home/kaoku/reference/mouse/mm10

ls ../clean/*gz | while read id;
do 
file=$(basename $id )
sample=${file%%.*}
echo $file $sample

bowtie2  -p 5  -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done 

##大约4-5分钟一个结果，耐心等待（8G内存）

##约1.2小时，跑完全部bam文件
(chipseq) root 16:10:27 /home/kaoku/chipseq/mouse_project/align
$ ls -lh
总用量 9.4G
-rw-r--r-- 1 root root 369M  8月 10 15:00 Control_1_trimmed.bam
-rw-r--r-- 1 root root 463M  8月 10 15:03 Control_2_trimmed.bam
-rw-r--r-- 1 root root 496M  8月 10 15:07 H2Aub1_1_trimmed.bam
-rw-r--r-- 1 root root 759M  8月 10 15:12 H2Aub1_2_trimmed.bam
-rw-r--r-- 1 root root 824M  8月 10 15:19 H3K36me3_1_trimmed.bam
-rw-r--r-- 1 root root 764M  8月 10 15:23 H3K36me3_2_trimmed.bam
-rw-r--r-- 1 root root 235M  8月 10 15:24 Ring1B_1_trimmed.bam
-rw-r--r-- 1 root root 235M  8月 10 15:26 Ring1B_2_trimmed.bam
-rw-r--r-- 1 root root 751M  8月 10 15:33 Ring1B_3_trimmed.bam
-rw-r--r-- 1 root root 421M  8月 10 15:36 RNAPII_8WG16_1_trimmed.bam
-rw-r--r-- 1 root root 679M  8月 10 15:42 RNAPII_8WG16_2_trimmed.bam
-rw-r--r-- 1 root root 716M  8月 10 15:50 RNAPII_S2P_1_trimmed.bam
-rw-r--r-- 1 root root 298M  8月 10 15:52 RNAPII_S2P_2_trimmed.bam
-rw-r--r-- 1 root root 478M  8月 10 15:56 RNAPII_S2P_3_trimmed.bam
-rw-r--r-- 1 root root 589M  8月 10 16:00 RNAPII_S5P_1_trimmed.bam
-rw-r--r-- 1 root root 603M  8月 10 16:05 RNAPII_S5P_2_trimmed.bam
-rw-r--r-- 1 root root 215M  8月 10 16:06 RNAPII_S5PRepeat_1_trimmed.bam
-rw-r--r-- 1 root root 411M  8月 10 16:08 RNAPII_S7P_1_trimmed.bam
-rw-r--r-- 1 root root 304M  8月 10 16:10 RNAPII_S7P_2_trimmed.bam

##对bam文件进行QC
ls  *.bam  | xargs -i samtools index {} 
ls  *.bam  | while read id ;do ( samtools flagstat $id > $(basename $id ".bam").stat & ); done

##查看结果
(chipseq) root 16:53:09 /home/kaoku/chipseq/mouse_project/align
$ ls -lh
总用量 9.5G
-rw-r--r-- 1 root root 369M  8月 10 15:00 Control_1_trimmed.bam
-rw-r--r-- 1 root root 3.6M  8月 10 16:13 Control_1_trimmed.bam.bai
-rw-r--r-- 1 root root  473  8月 10 16:15 Control_1_trimmed.stat
-rw-r--r-- 1 root root 463M  8月 10 15:03 Control_2_trimmed.bam
-rw-r--r-- 1 root root 3.3M  8月 10 16:13 Control_2_trimmed.bam.bai
-rw-r--r-- 1 root root  473  8月 10 16:15 Control_2_trimmed.stat
-rw-r--r-- 1 root root 496M  8月 10 15:07 H2Aub1_1_trimmed.bam
-rw-r--r-- 1 root root 3.0M  8月 10 16:13 H2Aub1_1_trimmed.bam.bai
-rw-r--r-- 1 root root  473  8月 10 16:15 H2Aub1_1_trimmed.stat
-rw-r--r-- 1 root root 759M  8月 10 15:12 H2Aub1_2_trimmed.bam
-rw-r--r-- 1 root root 2.6M  8月 10 16:14 H2Aub1_2_trimmed.bam.bai
-rw-r--r-- 1 root root  477  8月 10 16:15 H2Aub1_2_trimmed.stat
-rw-r--r-- 1 root root 824M  8月 10 15:19 H3K36me3_1_trimmed.bam
-rw-r--r-- 1 root root 2.8M  8月 10 16:14 H3K36me3_1_trimmed.bam.bai
-rw-r--r-- 1 root root  477  8月 10 16:15 H3K36me3_1_trimmed.stat
-rw-r--r-- 1 root root 764M  8月 10 15:23 H3K36me3_2_trimmed.bam
-rw-r--r-- 1 root root 3.0M  8月 10 16:14 H3K36me3_2_trimmed.bam.bai
-rw-r--r-- 1 root root  477  8月 10 16:15 H3K36me3_2_trimmed.stat
略

##查看比对结果
ls *.stat | while read id;  do echo $id; cat $id | awk 'NR'==7; done > result1
sed 'N; s/\n/ /' result1

Control_1_trimmed.stat 7438540 + 0 mapped (88.03% : N/A)
Control_2_trimmed.stat 7221780 + 0 mapped (86.40% : N/A)
H2Aub1_1_trimmed.stat 8969576 + 0 mapped (97.40% : N/A)
H2Aub1_2_trimmed.stat 13229915 + 0 mapped (97.53% : N/A)
H3K36me3_1_trimmed.stat 11737311 + 0 mapped (98.89% : N/A)
H3K36me3_2_trimmed.stat 13414318 + 0 mapped (98.45% : N/A)
Ring1B_1_trimmed.stat 4634240 + 0 mapped (93.59% : N/A)
Ring1B_2_trimmed.stat 4646919 + 0 mapped (93.85% : N/A)
Ring1B_3_trimmed.stat 22937484 + 0 mapped (95.34% : N/A)
RNAPII_8WG16_1_trimmed.stat 7472851 + 0 mapped (96.16% : N/A)
RNAPII_8WG16_2_trimmed.stat 20697874 + 0 mapped (95.50% : N/A)
RNAPII_S2P_1_trimmed.stat 25018794 + 0 mapped (97.26% : N/A)
RNAPII_S2P_2_trimmed.stat 6112834 + 0 mapped (95.00% : N/A)
RNAPII_S2P_3_trimmed.stat 8675513 + 0 mapped (96.99% : N/A)
RNAPII_S5P_1_trimmed.stat 11801469 + 0 mapped (97.53% : N/A)
RNAPII_S5P_2_trimmed.stat 12182274 + 0 mapped (98.17% : N/A)
RNAPII_S5PRepeat_1_trimmed.stat 4163762 + 0 mapped (82.81% : N/A)
RNAPII_S7P_1_trimmed.stat 6386270 + 0 mapped (80.90% : N/A)
RNAPII_S7P_2_trimmed.stat 5971177 + 0 mapped (82.66% : N/A)

##可以看到比对结果还是挺好的

获得bam文件后，下一步就是判断是否去除PCR重复。

那么我们下一篇再见！

【ChIP-seq 实战】四、fastq数据的质量控制和比对
这里是佳奥，继续处理数据吧。 1 质量控制（trim_galore）有点窒息，这个结果所以需要进一步过滤。发...
ChIP-Seq分析
植物转录因子 ChIP-Seq 实战系列 - 质量控制篇 https://mp.weixin.qq.com/s?_...
RNA-Seq数据分析—fastp v0.23.1
用于执行FASTQ数据的质量控制，支持单端和双端短读数据，作为一体化的 FASTQ预处理器，fastp提供了诸...
Biostar handbook学习笔记五-基因组测序技术原理简
测序技术及原理比较测序质量控制： FASTQ文件中测序Reads需要与指定的参考基因组进行序列比对，定位cDNA...
RNA-Seq数据分析——原始数据质量控制（QC）
获得转录组数据（.fastq文件）后的第一步就是对原始数据的质量控制。目的质量控制的目的是全面查看原始数据的质...
ChIP-seq数据分析实战训练（二）
ChIP-seq数据分析实战训练（二）四、过滤 raw data -> fastqc 使用trim_galore...
生信文件格式fastqc
资料推荐生信菜鸟团的浅谈FastQ和FastA格式,以及测序数据质量控制之FastQC 生信技能书论坛的blat...
临床生信实操（centos7）1
1. 原始数据（fastq） → 质量控制（fastqc）操作及路径：在fastqc路径下双击打开可以修改脚本，...
RNA-Seq数据分析—cutadapt质量控制
可对fastq数据进行质量控制的软件有很多，今天我们先学习cutadapt。安装官方帮助文档常用参数我处理...
RNA-seq分析：从fastq到差异表达基因
RNA-seq的数据分析是比较简单基础的分析，大概流程就是处理下机的fastq数据（trimmomatic），比对...