基因组

前言：

    一般测序下机数据会存在含N比例过大、测序质量较低的碱基数占比过高、含有duplication、序列污染等低质量reads,这些不合格的reads会影响后续的分析，所以，我们拿到测序数据首先要了解测序数据的质量情况，具体内容包括含N比例、GC含量、duplication情况、序列长度分布情况、碱基平衡情况等。

    今天，我们将一起通过数据格式和质量体系、数据质控步骤、Fastqc结果解读及异常处理三大模块进行学习。

第一部分 数据格式和质量体系

    Illumina测序的下机数据一般为fastq格式，至于fastq格式的说明我已经在上期《测序技术原理及常用数据格式简介》中有详细描述，在进行数据质控前，我们需要知道数据中第四行质量字符和序列质量值Q值的关系以及Q值与碱基测序错误率的关系。

    Fastq数据中的质量字符并不是和质量值Q值直接对应起来的，fastq数据格式中的质量字符是ASCII值，在Phred+64体系中，ASCII值-64的结果就是Q值，在Phred+33体系中，ASCII值-33的结果就是Q值。在Phred+33体系中，Q = -10log10(P), 碱基质量值与误率的对应关系表如下所示：

即，Q10准确率为90%，Q20准确率为99%，Q30准确率为99.9%，Q40准确率为99.99%，Q50准确率为99.999%。

第二部分 数据质控

    数据质控现在用得最多的是fastqc,我们今天就以它为工具学习如何了解测序数据质量。

Fastqc下载安装

wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip

unzip fastqc_v0.11.5.zip

cd FastQC

chmod +x fastqc

Fastqc评估测序数据质量

Usage:

fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file][-t]

--(no)extract输出的结果不接压，若无此选项，输出的结果为.zip压缩文件。

-f fastq|bam|sam指定输入文件格式，若无此项，则会自动检测。

-c contaminant file指定一个contaminant文件，文件格式为”Name\tSequence”，fastqc会把overrepreseted sequence往这个contaminant文件搜索。

-t线程数

例子：

fastqc *fq.gz –t 4  #目录下所有fq.gz文件进行质控，线程数一般与文件数一致。

第三部分 fastqc结果解读及异常处理

    横轴代表碱基在序列中的位置，纵轴代表Q值，由前面碱基质量值与错误率的关系可知，若某个位置对应的Q值为30，则该处碱基测序准确率为99.9%。

    如Figure1所示，在箱线图中，红色表示中位数，黄色是25%-75%区间，触须是10%-90%区间，蓝线是平均数。若任一位置的下四分位数低于10或中位数低于25，报"WARN"；若任一位置的下四分位数低于5或中位数低于20，报"FAIL"。

    当出现任一位置的下四分位数低于10或中位数低于25或任一位置的下四分位数低于5或中位数低于20时，表示测序数据存在质量不合格的情况，这时我们可以继续观察Sequence Contentacross图、GC Contentacross all base图、N Content across all bases图、Sequences Duplication level stastics图这几个图进一步判断测序数据的不合格之处具体在哪。

    观察Sequence Contentacross图和GC Contentacross all base图的GC含量的线是否平行于X轴，若不平行，则该位置往往有overrepresentedsequence的污染，可能原因建库过程的误差、测序的系统误差或者文库本身特点。

    由N Content across all bases图可知reads中含N碱基的情况，理想状况下是含N量越少越好，在微生物多样性分析中一般是去除含N碱基比例>5%的序列。

    去除含N比例过高的序列，可以用NGS QC tookit。该工具可在http://www.nipgr.res.in/ngsqctoolkit.html处下载，解压缩在QC/RIMINGREADS文件夹中用AmbiguityFiltering.pl脚本去除含N比例过高的序列，具体使用方法可以参考软件压缩包内的manual文件。

    观察Sequences Duplication level stastics图，横坐标是duplication的次数，纵坐标是duplicated reads的数目，若duplication的程度偏高，则可能存在PCR duplication。去除duplication可以通过Samtools、Picard或Iontorrent，其中Samtools只看5’端的起始位置不考虑reads突变；Picard不仅考虑起始位点也会考虑突变情况和质量值，即reads完全一样的才会被当成duplication被去除；Iontorrent则是看5’端的起始位置和3’端adaptor的比对情况，不考虑reads突变。

    去除PCR duplication的具体操作可参考：仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的。

    最后，在进行去低质量reads和接头等预处理步骤后，再次进行fastqc质控，然后用Multic QC即可把多个样品的质控结果汇总在一起，在报告中图片是交互式的，鼠标停留可显示样品名。如下图所示。

    当然，要是不想经历fastqc、NGS QC tookit、cutadapter、trimomatic这些工具这么麻烦，想省事点，直接一步到位完成所有测序数据预处理步骤，有没有这样的神器呢？有的，fastp可以做到。常规测序预处理流程一般包括Fastqc等软件进行质控、cutadapter去除接头、Trimmonatic进行过滤低质量reads、切除长度不合格的序列等步骤，涉及软件较多，比较繁琐，直接用fastp则可以一步到位完成。详情请点击https://github.com/OpenGene/fastp自行了解。