ChIP-seq数据分析课程学习笔记之 测序数据质量控制和比对

其中中国医科大的“小高”同学给大家带来的就是ChIP-seq数据分析实战视频课程的配套笔记，希望可以帮助大家更好的吸收消化课程内容！

首先视频免费共享在B站：【生信技能树】Chip-seq测序数据分析ChIP-SEQ实战演练的素材：链接：https://share.weiyun.com/53CwQ8B 密码：ju3rrh， 包括一些公司PPT，综述以及文献 ChIP-SEQ 实战演练的思维导图：文档链接：https://mubu.com/doc/11taEb9ZYg 密码：wk29

接下来是根据生信技能树课程、思维导图、PPT和综述文献整理的笔记。前面的笔记是：ChIP-seq数据分析课程学习笔记之fastq测序数据的准备

8-fastq数据的质量控制

前面我们已经下载好了全部的fastq数据，接下来就

使用trim_galore软件进行质控

先走一波QC

cd /data/gy/project/epi/qc
## 相对目录需要理解
ls ../raw/*gz|xargs fastqc -t 10 -o  ./

FastQC结果解读

官网上有个相关的视频教程：Using FastQC to check the quality of high throughput sequence，链接：https://www.youtube.com/watch?v=bz93ReOv87Y。（我上传了百度云链接: https://pan.baidu.com/s/1El3ZCli9uVC7u-58tj3HBw 密码: 2wj2）帮助文档在：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/。

左侧可以看到大致的结果：绿色PASS、黄色Warning、红色Failure

1、Basic Statistics

包含基本信息：测序平台的版本和相应的编码版本号、reads的数量、长度、GC含量等信息。

2、Per base sequence quality

image-20210707160448002

所有碱基的测序质量的统计结果。

此图横轴是测序序列1-36个碱基。纵轴是质量得分，20表示0.01的错误率，30表示0.001。

箱线图中红线表示中位数，蓝线是平均数的连线。

Warning 报警：任一碱基质量<10，或者是任一中位数<25 Failure 报错：任一碱基质量<5，或者是任一中位数<20

3、Per tile sequence quality

Per tile sequence quality

各个tile的序列质量情况。此图横轴是测序序列1-36个碱基。纵轴是tile的编号。

蓝色表示测序质量很好，颜色越红越不好，后续分析可以把这个tile结果去除。

4、Per sequence quality scores

Per sequence quality scores

各个序列质量的分布。横图表示Q值（测序质量quality），纵轴是read数目。绝大多数序列质量在30以上可以说质量很好。

5、Per base sequence content

Per base sequence content

每条序列中各个位置的平均碱基比例，如出现AT或GC分离的情况说明这个数据有问题，需要处理。

6、Per sequence GC content

Per sequence GC content

序列平均GC含量分布图。横轴是GC比例0 - 100%；纵轴是每条序列GC含量对应的数量。

蓝色：理论值。红色：真实值。两条越接近约好。

7、Per base N content

Per base N content

序列中各个位点的N含量（测序机器识别不出到底是何种碱基时，会产生N），越小越好。

8、Sequence Length Distribution

Sequence Length Distribution

序列测序长度统计。每次测序仪测出来的长度在理论上应该是完全相等的，但是实际总会有一些偏差。此图中，36bp是是主要的，但是还是有少量的35和37bp。当测序的长度严重不同时，表明此次测序不成功。

9、Sequence Duplication Levels

Sequence Duplication Levels

统计序列的重复度。低重复水平可能表明目标序列的覆盖率非常高，但高重复水平可能表明存在某种富集偏倚，例如 PCR over amplification。

10、Overrepresented sequences

Overrepresented sequences

某个序列大量出现是over-represented。

11、Adapter Content

Adapter Content

这个图表示序列中两端adapter的情况。一般测序仪下机之前会去除接头序列。

运行MultiQC

友情提示MultiQC安装可能会遇到坑：

依赖python2.7+, 3.4+ 或者 3.5+

#pip安装
pip install git+https://github.com/ewels/MultiQC.git
# conda安装
conda install -y bioconda multiqc

直接指定MultiQC要分析的文件路径即可，若数据在当前目录下输入multiqc ./即可。

并不是很好

GC含量还可以

以上为简要介绍，更多精彩内容可以参加**生信入门课**

这个时候选择trim_galore软件进行过滤，过滤条件：测序得到的原始序列含有接头序列或低质量序列，为了保证信息分析的准确性， 需要对原始数据进行质量控制，得到高质量序列(即Clean Reads)，原始序列质量控制的标准为:

①去除含接头的reads;

②过滤去除低质量值数据，确保数据质量;

③去除含有N(无法确定碱基信息)的比例大于5%的reads;

trim_galore常用参数

单端测序数据的代码如下；

cd /data/gy/project/epi/clean
#改成自己的
analysis_dir=/data/gy/project/epi
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 
#再QC一下
ls ../clean/*gz|xargs fastqc -t 10 -o  ./
multiqc ./

9-比对及质控

使用bowtie2进行比对

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

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

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

单端测序数据的比对代码如下：

cd /data/gy/project/epi/align
## 相对目录需要理解
bin_bowtie2='/root/miniconda3/envs/chipseq/bin/bowtie2'
bin_bowtie2=bowtie2
bowtie2_index=/data/gy/public/reference/index/bowtie/mm10
## 一定要搞清楚自己的bowtie2软件安装在哪里，以及自己的索引文件在什么地方！！！
ls ../clean/*gz |while read id;
do 
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
$bin_bowtie2  -p 5  -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done 

对bam文件进行QC

cd /data/gy/project/epi/align
ls  *.bam  |xargs -i samtools index {} 
ls  *.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done

比对成功率都挺好的：

Control_1_trimmed.stat:7438540 + 0 mapped (88.03% : N/A)
Control_2_trimmed.stat:7221781 + 0 mapped (86.40% : N/A)
H2Aub1_1_trimmed.stat:8969578 + 0 mapped (97.40% : N/A)
H2Aub1_2_trimmed.stat:13229916 + 0 mapped (97.53% : N/A)
H3K36me3_1_trimmed.stat:11737310 + 0 mapped (98.89% : N/A)
Ring1B_1_trimmed.stat:4634240 + 0 mapped (93.59% : N/A)
Ring1B_2_trimmed.stat:4646919 + 0 mapped (93.85% : 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:8675514 + 0 mapped (96.99% : N/A)
RNAPII_S5P_2_trimmed.stat:12182274 + 0 mapped (98.17% : N/A)
RNAPII_S5PRepeat_1_trimmed.stat:4163763 + 0 mapped (82.81% : N/A)
RNAPII_S7P_1_trimmed.stat:6386269 + 0 mapped (80.90% : N/A)
RNAPII_S7P_2_trimmed.stat:5971178 + 0 mapped (82.66% : N/A)

题外话：真的是 大神一句话，菜鸟跑半年。这两节加起来20分钟的课程，感觉自己踩的坑大概比这篇笔记的字数还多。经验教训：多读文献，多做实战，shell脚本还要好好看。