NGS数据分析实践：05. 测序数据的基本质控 [1] - FastQC
NGS数据分析实践：05. 测序数据的基本质控 [1] - FastQC 目录
前言
1. FastQC
1.1 帮助信息及运行代码
1.2 报告解读
1.3 小结
文接上篇：NGS数据分析实践：04. 准备测序数据
前言
目前，以illumina为首的NGS测序，一般均采用边合成边测序的技术。碱基的合成依靠的是化学反应，这使得碱基链不断地从5’端一直往3’端合成并延伸下去。但，随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，因此，越到后面碱基合成的错误率就会越高，这也是为何当前NGS测序读长普遍偏短的一个原因。有时测序仪在刚开始进行合成反应的时候也会由于反应还不够稳定，同样会带来质量值的波动，不过这个波动一般都在高质量值区域（如下图3）。
一般我们可以从如下几个方面来分析测序数据质量：
read各个位置的碱基质量值分布 (Per base sequence quality)
碱基的总体质量值分布 (Per sequence quality scores)
read各个位置上碱基分布比例 (Per base sequence content)
GC含量分布 (Per sequence GC content)
read各位置的N含量 (Per base N content)
read是否还包含测序的接头序列 (Adapter Content)
read重复率，这个是实验的扩增过程所引入的 (Sequence Duplication Levels)
常用的测序质量分析工具为FastQC和MultiQC。
1. FastQC
FastQC是一款基于Java的软件，一般在linux环境下使用命令行运行，它可以快速多线程地对测序数据进行质量评估（Quality Control）。在做read质量值分析的时候，FastQC并不单独查看具体某一条read中碱基的质量值，而是将Fastq文件中所有的read数据都综合起来一起分析。
1.1 帮助信息及运行代码
fastqc --help
# 基本格式# fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
# 主要是包括前面的各种选项和最后面的可以加入N个文件
# -o --outdir FastQC生成的报告文件的储存路径，生成的报告的文件名是根据输入来定的
# --extract 生成的报告默认会打包成1个压缩文件，使用这个参数是让程序不打包
# -t --threads 选择程序运行的线程数，每个线程会占用250MB内存，越多越快咯
# -c --contaminants 污染物选项，输入的是一个文件，格式是Name [Tab] Sequence，里面是可能的污染序列，如果有这个选项，FastQC会在计算时候评估污染的情况，并在统计的时候进行分析，一般用不到
# -a --adapters 也是输入一个文件，文件的格式Name [Tab] Sequence，储存的是测序的adpater序列信息，如果不输入，目前版本的FastQC就按照通用引物来评估序列时候有adapter的残留
# -q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况。
代码如下：
dat_batch1=/filepath
ls $dat_batch1/Cleandata/ | while read file; do fastqc $file -t 10 -o fastqc_out_Clean/ ; done
# ls $dat_batch1/Rawdata/ | while read file; do fastqc $file -t 10 -o fastqc_out_Raw/ ; done
# or
mkdir fastqc_out
ls *gz |xargs fastqc -t 10 -o fastqc_out # --extract
cd fastqc_out
ls *zip|while read id;do unzip $id; done #把所有压缩包批量解压开。每个测序数据都含十几项统计结果和可视化的图片
1.2 报告解读
使用浏览器打开后缀是html的文件，就是图表化的fastqc报告。
(1) Summary
结果分为：绿色表示“PASS（通过）”，红色表示“FAIL（未通过）”，黄色表示“WARN（警告，不太好）”。
图1：summary - 整体情况
从页面左侧的的summary中可以看出有哪些选项没有通过，上图可以看出此数据的测序质量尚可。
(2) Basic Statics
图2：Basic Statics - 基本统计信息
可以看出数据的序列数量，测序平台以及GC含量等相关信息：
A.Encoding：Sanger / Illumina 1.9 (Phred33质量值体系)；
B.Sequence length：3-150 (来自HiSeq X10的测序结果，read长度达150bp，应该算是比较长的二代测序序列了)；
C.%GC：45 (人类基因组的GC含量一般在40%左右)。
D.Total sequences: reads数量（可译为读段）。
(3) Per base sequence quality
可以看到read上的每一个位置都有一个黄色的箱线图表示在该位置上所有碱基的质量分布情况。除了最后一个碱基之外，其他的碱基质量值都基本都在大于30，而且波动较小，说明质量很稳定，这其实是一个比较高质量的结果。而且可以看到图中质量值的分布都在绿色背景（代表高质量）的区域。
图3：Per base sequence quality - 每个位置的碱基的质量情况
A.此图中的横轴是测序序列第1个碱基到第150个碱基；纵轴是质量得分（Q值），Q = -10*log10（P error），即20表示错误率为1%，30表示错误率为0.1%。
B.图中每1个boxplot，都是该位置的所有序列的测序质量的一个统计，上面的bar是90%分位数，下面的bar是10%分位数，箱子的中间的横线是50%分位数，箱子的上边是75%分位数，下边是25%分位数。
C.图中蓝色的细线是各个位置的平均值的连线，一般要求此图中的所有位置的10%分位数大于20，也就是我们常说的Q20过滤。
D.WARN (对应summary的⚠️) 出现的情况：任何碱基质量低于10，或者是任何中位数低于25；
E.FAIL (对应summary的❌)：任何碱基质量低于5，或者是任何中位数低于20。
(4) Per tile sequence quality
这个图显示了各个tile的序列质量情况。
简单介绍几个名词：测序反应的载体/容器称flowcell；测序反应的平行泳道称为lane，lane是试剂添加、洗脱等过程的发生位置；每次荧光扫描的位置称为tile，肉眼是看不到的。1个flowcell有8个lane，1个lane包含2列（swath），每1列有60个tiles，每个tile会种下不同的cluster，每个tile在一次循环中会拍照4次（每个碱基一次）。相关概念详见：测序的世界 https://www.jianshu.com/p/101c14c3a1d2。
图4：Per tile sequence quality - 每个tile测序的情况
A.横轴和之前一样，代表150个碱基的每个不同位置；纵轴是tile的Index编号。
B.这个图主要是为了防止在测序过程中，某些tile受到不可控因素的影响而出现测序质量偏低。
C.蓝色代表测序质量很高，暖色代表测序质量不高，如果某些tile出现暖色，可以在后续分析中把该tile测序的结果全部都去除。
(5) Per sequence quality scores
这个图可以看出各个序列质量的分布情况，可看到绝大部分序列质量都在30以上，质量很好。
图5：Per sequence quality scores - 每条序列的测序质量统计
A.假如我测的1条序列长度为150bp，那么这150个位置每个位置Q值的平均值就是这条reads的质量值。
B.该图横轴是17-41，表示Q值；纵轴是每个值对应的reads数目。测序结果主要集中在高分中，证明测序质量良好。
C.当峰值小于27（错误率0.2%）时报“WRAN”，当峰值小于20（错误率1%）是报“FAIL”。
(6) Per base sequence content
对所有reads的每个位置，统计ATGC四种碱基的分布。这个图可以看出每条序列中各个位置的平均碱基比例，如果出现AT或GC分离的情况说明这个数据有问题，需要处理。由于AT配对，CG配对，假如测序过程是比较随机的话（随机意味着好），那么在每个位置上A和T比例应该差不多，C和G的比例也应该差不多。如下图所示，两者之间即使有偏差也不应该太大，最好平均在1%以内，除非有合理的原因，比如某些特定的捕获测序所致的偏差过高。
图6：Per base sequence content - 每个位置上的碱基的比例分布
A.横轴是位置1 - 150 bp；纵轴是百分比。
B.图中四条线代表A T C G在每个位置平均含量。理论上来说，A和T应该相等，G和C应该相等，但是一般测序的时候，刚开始测序仪状态不稳定，很可能出现严重分离的情况。像这种情况，即使测序的得分很高，也需要cut开始部分的序列信息，一般像这种情况，会cut前面5-10bp。
C.正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该是平行且接近的。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示有overrepresented sequence 的污染，当所有的碱基比例一致地表现出bias时，即四条线平行但分开，往往代表文库有bias（建库过程或本身特点），或者是测序中的系统误差。当任一位置的A/T比例与G/C比例相差超过10%，报“WRAN”，当任一位置的A/T比例与G/C比例相差超过20%，报“FAIL”。
(7) Per sequence GC content
GC含量指的是G和C这两中碱基占总碱基的比例。二代测序平台或多或少都存在一定的测序偏向性，我们可以通过查看这个值来协助判断测序过程是否足够随机。对于人类来说，我们基因组的GC含量一般在40%左右。因此，如果发现GC含量的图谱明显偏离这个值那么说明测序过程存在较高的序列偏向性，结果就是基因组中某些特定区域被反复测序的几率高于平均水平，除了覆盖度会有偏离之后，将会影响下游的变异检测和CNV分析。
图7：Per sequence GC content - read的GC含量的频率分布图
A.横轴是0 - 100%；纵轴是每条序列GC含量对应的数量。
B.蓝色的线是程序根据经验分布给出的理论值，红色是真实值，两个应该比较接近才比较好。
C.当红色的线出现双峰，基本肯定是混入了其他物种的DNA序列。
D.偏离理论分布的reads超过15%时，报“WRAN”， 偏离理论分布的reads超过30%时，报“FAIL”。
(8) Per base N content
序列中各个位点的N含量，越小越好。正常情况下N的比例很小，所以图上常常看到一条直线。当任一位置的N的比例超过5%，报“WRAN”， 当任一位置的N的比例超过20%，报“FAIL”。
图8：Per base N content - N含量分布图
N在测序数据中一般是不应该出现的，如果出现则意味着，测序的光学信号无法被清晰分辨，如果这种情况多的话，往往意味着测序系统或者测序试剂的错误。
(9) Sequence Length Distribution
序列测序长度统计，从图中可以看出序列的平均长度为150。
图9：Sequence Length Distribution - 序列测序长度分布
A.每次测序仪测出来的长度在理论上应该是完全相等的，但是总会有一些偏差。比如此图中，150bp是主要的，但是还是有少量的149和151bp的长度，不过数量比较少，不影响后续分析。
B.当测序的长度不同时，如果很严重，则表明测序仪在此次测序过程中产生的数据不可信。
C.当reads长度不一致时报“WRAN”，当有长度为0的read时报“FAIL”。
(10) Sequence Duplication Levels
统计序列完全一样的reads的频率。测序深度越高，越容易产生一定程度的duplication，这是正常现象，但如果duplication的程度越高，就提示可能有bias的存在。
图10：Sequence Duplication Levels - read重复的频率分布
A.sequences duplication是指在测序前建库PCR过程中导致的一些序列扩增次数过多导致的。若重复较高则需要进行处理这些dup。
B.横轴代表 reads 的重复次数 ( 1 表示 unique 的序列，2 表示有 2 条完全相同的 reads …)，大于 10 次重复后则按不同的重复次数合并显示。纵坐标表示各重复次数下的 reads 数占总 reads 的百分比。
C.蓝线展示所有 reads 的重复情况，红线表示在去掉重复以后，原重复水平下的 reads 占去重后 reads 总数的百分比。
D.需要注意的是，这个重复水平和建库方式有关。如果是 GBS 建库，由于存在 PCR 扩增的步骤，观察到的重复水平会比较高。
E.注：fastqc中用fastq数据的前100,000条reads统计其在全部数据中的重复情况。当非unique reads占总数的比例大于20%时，报“WRAN”， 当非unique reads占总数的比例大于50%时，报“FAIL”。
(11) Overrepresented sequences
Overrepresented sequences 显示同一条 read 出现次数超过总测序 reads 数的 0.1 % 的统计情况。
图11：Overrepresented sequences
A.正常文库内序列的多样性水平很高，不会有同一条 read 大量出现的情况，这部分结果会把大量出现的 reads 列出来，并给出可能来源。
B.在这部分分析中任何长度大于 75 bp 的序列都会被截成 50 bp。
C.当发现超过总reads数0.1%的reads时报“WRAN”， 当发现超过总reads数1%的reads时报“FAIL”。
(12) Adapter Content
Adapter Content 显示 reads 中的接头含量，并显示可能的来源。
图12：Adapter Content
A.横轴为碱基位置，纵轴为含有接头序列的比例。
B.如接头的含量随着碱基的位置增大而逐渐升高，则表示 reads 中含有接头，这部分接头会影响后续的分析，需要截掉 reads 中的接头序列或者将含有接头的 reads 完全删除。
C.如果任何重复 read 超过总 reads 数的 5 % 则报 Warn， 超过总 reads 数的 10 % 则报 Fail。
1.3 小结
原始测序数据经过fastqc质检后，很少有全部通过的情况，一般都会出现一些warning，因此没有全部通过并不意味着不能进行后续分析。但是前提条件是几个关键指标不能太差，根据经验，一般需要重点关注的主要是 Per base sequence quality、Per base sequence content和Adapter Content。
其中，如果Per base sequence quality太差的话，说明数据的质量远没有达到符合要求的Q30或着Q20的比例，这样测到的reads很多碱基是不可信的，对下游的分析结果影响比较大。
如果Per base sequence content的结果中出现很大异常的话(比如碱基的曲线出现明显的波动)，很可能提示原始下机数据中出现了很大比例的重复reads，这些重复的reads虽然本身的测序质量可能没有问题，但是有可能导致最终可用于分析的clean reads 大大减少，需要引起注意。
如果Adapter Content参数曲线中，出现很大比例的adapter（接头）序列的话，一般需要先根据接头序列先去掉接头序列再进行分析的。否则可能会影响后续的比对分析结果。
QC的结果的评判还需要结合具体的项目，测序平台以及分析目的这些因素，不同的因素导致的判断也是不同的。比如是RNA-seq 还是DNA-seq ，是否是捕获测序，是否是酶切后测序，是否是多重PCR测序等等。
参考阅读：