RNA-Seq数据分析——原始数据质量控制(QC)
获得转录组数据(.fastq文件)后的第一步就是对原始数据的质量控制。
目的
质量控制的目的是全面查看原始数据的质量,内容包括碱基质量评估、GC含量检验、N碱基数量评估、TCGA碱基分布、k-mer数量检验等。
方法
可以于检验fastq文件质量的软件有很多,例如FastQC、fastp、multiQC等。本文主要介绍应用最多的FastQC。
FastQC是一款基于Java的软件,须在linux环境下使用命令行运行,它可以快速多线程地对测序数据进行质量评估(Quality Control),其官网地址为:Babraham Bioinformatics。
安装
FastQC可以使用conda进行安装。在linux环境下运行命令conda install fastqc
即可,运行结果如下图。
运行命令fastqc -h
可检验其是否成功安装,运行结果如下图。
运行
# 运行命令的基本格式
# 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 安静运行模式,一般不选这个选项的时候,程序会实时报告运行的状况。
使用fastqc -o #输出结果全路径 #数据存储全路径/*reads_R1.fq
命令运行案例数据
运行后可获得如下结果。
Fig.3结果解读
Basic Statistics 基本信息
Fig.4报告第一部分既是对质量检测结果的基本信息统计,如上图所示。其中包括:
-
Filename:检测的fastq文件名称;
-
File type:文件类型;
-
Encoding:测序平台的版本和相应的编码版本号;
-
Total Sequence:总reads数;
-
Sequences flagged as poor quality:低质量序列数量;
-
Sequence:测得的序列长度范围;
-
%GC:GC含量。
Per base sequence quality 序列测序质量统计
Fig.5上图显示了检测fastq文件的整体碱基质量分数统计。
-
横轴表示测序文件中所有序列第一个碱基到最后一个碱基,纵轴表示质量得分;
-
红线表示中位数,蓝线代表平均值;
-
柱状表示该位置所有序列的测序质量的统计,柱状(黄色)是25%~75%区间质量分布,error bar(触须)是10%~90%区间质量分布;
-
一般要求所有位置的10%小于20,即最多允许该位置10%的序列低于Q20,即90%的序列的碱基质量都大于Q20,即90%的序列碱基错误率不超过99%。当任何碱基质量低于10,或者任何中位数低于25时报WARN;当任何碱基质量低于5或者任何中位数低于20报FAIL。
Per tile sequence quality 每个tail测序的情况
Fig.6上图展示了每个tail的测序情况。
-
横轴表示每个碱基的位置;
-
纵轴是tail的Index编号;
-
这个图主要是为了防止,在测序过程中,某些tail受到不可控因素的影响而出现测序质量偏低;
-
蓝色代表测序质量很高,暖色代表测序质量不高,如果某些tail出现暖色,可以在后续分析中把该tail测序的结果全部都去除。
Per sequence quality scores 每条序列的测序质量统计
Fig.7对每条序列(reads)的测序质量统计。
-
假如我测的1条序列长度为101bp,那么这101个位置每个位置Q值的平均值就是这条reads的质量值;
-
该图横轴是0-40,表示Q值,即该序列(reads)质量得分;
-
纵轴是每个值对应的reads数目;
-
我们的数据中,测序结果主要集中在高分中,证明测序质量良好。
Per base sequence content 序列各个位置碱基比例分布
Fig.8上图显示了A T C G在每个位置的平均分布情况。
-
横轴表示每个碱基的位置,纵轴表示百分比;
-
图中四条线代表A T C G在每个位置平均含量;
-
理论上来说,A和T应该相等,G和C应该相等,但是一般测序的时候,刚开始测序仪状态不稳定,很可能出现上图的情况。像这种情况,即使测序的得分很高,也需要cut开始部分的序列信息。
Per sequence GC content 序列平均GC分布
Fig.9上图展示了序列平均GC分布。
-
横轴为平均GC含量; 纵轴为每个GC含量对应的序列数量;
-
蓝线为系统计算得到的理论分布;红线为测量值,二者越接近越好;
-
这里不相符可能有两个原因:
- GC可以作为物种特异性根据,这里出现了其他的峰有可能混入了其他物种的DNA;
- 目前二代测序基本都会有序列偏向性(所说的 bias),也就是某些特定区域会被反复测序,以至于高于正常水平,变相说明测序过程不够随机。这种现象会对以后的变异检测以及CNV分析造成影响。
Per base N content N碱基含量分布
Fig.10上图N碱基含量分布
- N碱基是指仪器不能识别的碱基,一般不会出现。但是如果出现并且量还很大,应该就是测序系统或者试剂的问题;
- 任意位置的N的比例超过5%,报"WARN";任意位置的N的比例超过20%,报"FAIL"。
Sequence Length Distribution 序列测序长度统计
Fig.11上图展示了检验文件中序列的长度统计。
-
每次测序仪测出来的长度在理论上应该是完全相等的,但是总会有一些偏差;
-
比如此图中,126-127bp是主要的,但是还是有少量的120-121bp的长度,不过数量比较少,不影响后续分析;
-
当测序的长度不同时,如果很严重,则表明测序仪在此次测序过程中产生的数据不可信。
Sequence Duplication Levels 统计序列完全一样的reads的频率
Fig.12-
谈到NGS数据的duplicated reads(暂且翻译为“重复数据”),我们通常会直观地认为:duplicated reads是在NGS文库构建过程中,由于PCR过度扩增导致同一个模板DNA片段被反复测序多次,得到一模一样的reads;
-
上图中横坐标是duplication的次数;纵坐标是duplicated reads的数目(红线);
-
正常情况下的确,测序深度越高,越容易产生一定程度的duplication。高程度的duplication level,提示我们可能有bias的存在(如建库过程中的PCR duplication)。
Overrepresented sequences 大量重复序列
Fig.13-
Overrepresented sequences是指一条序列的重复数,因为一个转录组中有非常多的转录本,一条序列再怎么多也不太会占整个转录组的一小部分(比如1%),如果出现这种情况,不是这种转录本巨量表达,就是样品被污染。这个模块列出来大于全部转录组1%的reads序列,但是因为用的是前200,000条,所以其实参考意义不大,完全可以忽略。
-
和duplication计算一样,取前200,000进行统计,大于75bp只取50bp;
-
发现超过总reads数0.1%的reads时报”WARN“,当发现超过总reads数1%的reads时报”FAIL“;
Adapter Content 序列Adapter
Fig.14-
此图衡量的是序列中两端adapter的情况
-
如果在当时fastqc分析的时候-a选项没有内容,则默认使用图例中的四种通用adapter序列进行统计
-
本例中adapter都已经去除,如果有adapter序列没有去除干净的情况,在后续分析的时候需要先使用cutadapt软件进行去接头。
接下来就是基于QC结果对数据进行质量控制,我们应用cutadapt来做。