软件2 —— fastqc和multiqc

一、 基本介绍

FastQC是一款基于Java的软件，可以快速多线程地对测序数据（基因组测序、ChIP-Seq、RNA-Seq、BS-Seq等）进行质量评估（Quality Control）。

MultiQC将多个样本的生物信息学分析结果汇总到一份报告中。它在给定的目录中搜索analysis logs，并汇编一个HTML报告。这是一个通用工具，非常适合汇总众多生物信息学工具的输出。运行时可以提供一个或多个目录以扫描分析结果。

二、 背景知识

(1) FastA格式

FASTA是存储有顺序的序列数据的文件的后缀，最后一个字母A，其实就是Alignment的意思。

>gi|187608668|ref|NM001043364.2| Bombyx mori moricin (Mor), mRNA
AAACCGCGCAGTTATTTAAAATATGAATATTTTAAAACTTTTTGTGGCAATGTCTCTGGTGTCATGTAGTACAGCCGCTCC

• 头信息：Fasta格式首先以大于号 > 开头，接着是序列的标识符，然后是序列的描述信息（FASTA的头信息并没有被严格地限制）。

• 序列信息：换行后是序列信息，序列中允许空格，换行，空行，直到下一个大于号，表示该序列的结束。

所有来源于NCBI的序列都有一个唯一的gi号“gi|gi_identifier”。gi号后面是序列的标识符，标识符由序列来源标识、序列标识（如接收号、名称等）等几部分组成，它们之间用“|”隔开，如果某项缺失，可以留空，但是“|”不能省略。“ref|NM001043364.2|”表示序列来源于NCBI的参考序列库，接收号为“NM_001043364.2”。

(2) FastQ格式

FastQ格式存储了生物序列以及相应的质量评价。最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。FASTQ文件中每个序列通常有四行：

• 第一行：必须以“@”开头，后面跟着唯一的序列ID标识符，然后跟着可选的序列描述内容，标识符与描述内容用空格分开。

• 第二行：序列字符（核酸为[AGCTN]，蛋白为氨基酸字符）。

• 第三行：必须以“+”开头，后面跟着可选的ID标识符和可选的描述内容。如果“+”后面有内容，该内容必须与第一行“@”后的内容相同。

• 第四行：碱基质量字符，每个字符对应第二行相应位置碱基或氨基酸的质量。该字符（例如ASCII字符?，对应值为63）可以按一定规则转换为碱基质量得分（63-33=30，Q30，即-10*lg0.001），碱基质量得分可以反映该碱基的错误率（例如0.001）。这一行的字符数与第二行中的字符数必须相同。

第一行的信息例如：

@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133

@：开始的标记符号
ST-E00126:128:HJFLHCCXX：测序仪唯一的设备名称
2：lane的编号
1101：tile的坐标
7405：在tile中的X坐标
1133：在tile中的Y坐标

碱基的质量值（Phred quality score）是测序错误率的对数（10为底数）乘以-10（并取整），即Q = -10 log10 P。P代表该碱基被测序错误的概率，p_error的值和测序时的多个因素有关，体现为测序图像数据点的清晰程度，并由测序过程中的base calling 算法计算出来。如果该碱基测序出错的概率为0.001，则Q应该为30，称为Q30。Q20和Q30的比例常常被我们用来评价某次测序结果的好坏，比例越高就越好。为了格式存储以及处理时的方便，这个数字被直接转换成了ASCII码，并与第二行的read序列构成一一对应的关系。小于33的ASCII码值所表示的都是不可见字符，因此加上一个固定的整数（33或64，现在一般都是使用Phred33体系，illumina测序1.8版本以上加33，以下加64）。那么Q30即30+33=63对应的ASCII码为“?”。一般地，碱基质量从0-40，既ASCII码为从 “!”（0+33）到“I”（40+33）。

三、 使用方法

(1) 用法和参数

fastqc  *.fq.gz  -t 4  [-o outputdir]  [--(no)extract]  [-f fastq|bam|sam]  [-c contaminantfile]
fastqc  *.fq.gz  -t 4  -o outputdir  -f fastq

*.fq.gz：可以输入一个或多个fq文件
-t / --threads：线程数
-o / --outdir：用于指定FastQC报告的输出目录，它是不能自动新建目录的。如果不指定特定的目录，那么FastQC的结果会默认输出到输入文件*.fq.gz的同一个目录下。生成的报告的文件名是根据输入来定的，它输出结果有两个，一个html和一个.zip压缩包。
-f：读入文件格式，会自动识别，可省略。
--extract：生成的报告默认会打包成1个压缩文件，使用这个参数是让程序不打包。
-c / --contaminants：污染物选项，输入的是一个文件，格式是Name [Tab] Sequence，里面是可能的污染序列。如果有这个选项，FastQC会在计算时候评估污染的情况，并在统计的时候进行分析，一般用不到。
-a / --adapters：输入一个文件，文件的格式Name [Tab] Sequence，储存的是测序的adpater序列信息。如果不输入，目前版本的FastQC就按照通用引物来评估序列adapter的残留。
-q / --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况。

(2) 举个例子

fastqc  testdata/WTF5.Input_1.fastq.gz  -t 4  -o ./01_fastqc  -f fastq  #单个文件
fastqc  testdata/*.fastq.gz  -t 4  -o ./01_fastqc  -f fastq  2>>log/fastqc_log.txt  &  #多个文件

(3) 使用multiqc整合报告

multiqc  ./01_fastqc  -o ./01_multiqc  2>>log/multiqc_log.txt

四、 fastqc结果解读

(1) Basic Statistics ——基本统计

 Filename：文件名。
 File type：文件类型，一般是Conventional base calls。
 Encoding：测序平台的版本和相应的编码版本号，用于计算Phred反推error P时用，一般是Sanger / Illumina 1.9。
 Total Sequences：输入的reads的数量，一般有几十M。
 Total Bases：输入的碱基数量，一般有几G。
 Sequence length：测序长度，例如50、100、150，也可能是类似20-50的一个范围。
 %GC：GC含量，表示整体序列的GC含量。这个数值一般是物种特异的，比如人类细胞就是42%左右。如果测序原始数据的GC含量远远偏离这个比例，说明测序数据存在一定偏好性。

(2) Per base sequence quality ——沿read每个碱基的质量值统计

横轴代表在read中的位置，纵轴代表quality。箱线图中间红色表示中位数，黄色区间是25%-75%，触须是10%-90%区间，蓝线是平均数。若任一位置的下四分位数低于10或中位数低于25，报"WARN"；若任一位置的下四分位数低于5或中位数低于20，报"FAIL"。目前边合成边测序的技术，碱基的合成依靠的是化学反应，这使得碱基链可以不断地从5'端一直往3'端合成并延伸下去。但随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，这就会带来一个问题——越到后面碱基合成的错误率就会越高。有时候测序仪在刚开始进行合成反应的时候也会由于反应还不够稳定，同样会带来质量值的波动，不过这个波动一般都在高质量值区域。

(3) Per tile sequence quality ——沿read每个tile测序情况

横轴表示碱基位置，纵轴表示tail的index编号。这个图主要是为了防止在测序过程中某些tail受到不可控因素的影响而出现测序质量偏低。蓝色表示测序质量很高，暖色表示测序质量不高。当某些tail出现暖色，在后续的分析中应把该tile测序结果全部去除。这一模块是检查reads中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色表示低于平均偏离度，越红则说明偏离平均质量方差越多，也就是说质量越差。如果出现质量问题可能是短暂的，如有气泡产生，也可能是长期的，如在某一小孔中存在残骸。有时分析结果会没有这一项。

(4) Per sequence quality scores ——reads的quality均值的分布

横轴为quality，纵轴是reads数目。当出现右下图的情况时，我们就会知道有一部分reads具有比较差的质量。当峰值小于27（错误率0.2%）时报"WARN"，当峰值小于20（错误率1%）时报"FAIL"。只要大部分quality都高于20，那么就比较正常。对于二代测序，最好是达到Q20的碱基要在95%以上（最差不低于90%），Q30要求大于85%（最差也不要低于80%）。

(5) Per base sequence content ——沿read四种碱基的分布

对所有reads的每一个位置，统计ATCG四种碱基的分布。横轴为位置，纵轴为百分比。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织（右下图），往往提示我们有overrepresented sequence的污染。当所有位置的碱基比例一致的表现出bias时，即四条线平行但分开，往往代表文库有bias（建库过程或本身特点），或者是测序中的系统误差。在碱基含量分布图中，前几个碱基可能会出现较大波动，这是由于随机引物扩增偏差原因造成的（左下图）。当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"。

(6) Per sequence GC content ——reads的平均GC含量的分布

红线是实际情况，蓝线是理论分布（正态分布，均值不一定在50%，而是由平均GC含量推断的；对于人类来说，我们基因组的GC含量一般在40%左右）。曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"。

(7) Per base N content ——沿read每个位置N的比例

当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”。正常情况下N的比例是很小的，所以图上常常看到一条直线，但放大Y轴之后会发现还是有N的存在，这不算问题。当Y轴在0%-100%的范围内也能看到“鼓包”时，说明测序系统出了问题。当任意位置的N的比例超过5%，报"WARN"；当任意位置的N的比例超过20%，报"FAIL"。

(8) Sequence Length Distribution ——reads长度的分布

当reads长度不一致时报"WARN"；当有长度为0的read时报"FAIL"。

(9) Sequence Duplication Levels ——重复序列的水平

统计序列完全一样的reads的频率。横坐标是duplication的次数，纵坐标是reads的百分数。蓝线代表在total reads中的重复情况，红色线代表Deduplicated reads的重复情况（新版只有%total reads了）。这里的deduplicated reads并不是代表只出现一次的reads，它是代表所有不同的reads（去重之后的reads）。比如说有20条reads，其中有10条出现1次，有5条分别出现2次，这里有15条reads是各不相同的，则total reads为20，deduplicated reads为15。表头的deduplicated percentage = deduplicated reads / total reads * 100%。下面的蓝色折线的意思是total reads中出现1次、2次、3次...的reads分别占total reads的比例；红色折线的意思是deduplicated reads中出现1次、2次、3次...的reads分别占deduplicated reads的比例。

例如，有20条reads，10条出现1次，1条出现10次，则产生下图。此时Total reads=20，deduplicated reads=11，deduplicated percentage=11/20=55%。当出现次数为1时，%total sequences=10/20=50%，%deplicated sequences=10/11=91%；当出现次数为10时，%total sequences=10/20=50%，%deplicated sequences=1/11=9%。

在RNA-seq数据的评估结果中，deduplicated reads有可能低于50%，即有一半可能存在重复。但在RNA-seq数据中，可能存在某些转录本较短，表达水平很高，导致在随机打断后被重复抽到测序的概率大大增加。因此对于RNA-seq数据来说，重复性较高的序列不一定就是PCR重复。例如左下图，观察两条折线知道，大部分序列的重复性在2-9次，它们极可能覆盖了一些“normal”基因；部分序列的重复性在10-100次，它们可能覆盖在rRNA或一些重复序列上；只有极少数序列的重复度超过100。因此，总体而言，该RNA-Seq数据可以代表一个多样性较丰富的文库。

而右下图则代表有高度的重复。首先该数据的deduplicated reads非常低，观察折线发现当出现次数为1时，%Deduplicated sequences为60%，而%Total sequences则低于5%；当出现次数大于1k时，%Deduplicated sequences在1%左右，而%Total sequences则在20%以上。这些数据说明了有极少数序列重复出现了上1千次，导致整个文库的重复率特别高，这少数部分的序列应该是PCR重复。

(9) Overrepresented sequences ——过表达的序列

如果有某个序列大量出现，就叫做over-represented。当发现（前200,000条）超过总reads数0.1%的over-represented reads时报"WARN"，当发现超过总reads数1%的over-represented reads时报"FAIL"。

(10) Adapter Content ——沿read接头的含量

横轴表示碱基位置，纵轴表示百分比。当fastqc分析时没有选择参数-a adapter list时，默认使用图例中的几种通用adapter序列进行统计。若有adapter残留，后续必须去接头。