二代测序的数据的分析——质量控制

2017-08-11 本文已影响0人科研虾

质量控制

测序质量检查-FastQC

Fastqc
Fastqc website (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

质量控制的测序质量检测是通过FastQC软件实现。fastqc可以不设置任何参数运行，这样会直接在当前目录下生成一个质量报告的压缩文件和文件夹，报告是网页格式。也可以设置输出目录和是否解压缩(--noextract)，默认设置会解压缩。命令如下：

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

其中--noextract命令是不解压缩输出文件。-t参数是指定使用线程数，fastqc似乎并不是并行运算，而是通过线程数同时执行多个程序，比如线程数指定为4，并不是用4个进程去跑一个文件，而是同时跑4个文件，不过4个线程速度提高很大，个人测试感觉10倍速度于2个线程。-q为屏蔽进程信息并只输出错误信息，-f参数为指定输入文件格式(有bam, sam, fastq可选)

fastqc的结果在v0.11.5版下共有12项。

Basic Statistics
包含文件名(Filename)、文件类型、总序列数(Total Sequences)、序列长度(Sequence length)这些基本信息。
Per base sequence quality
序列每个碱基的平均质量，越高越好，需要注意会有部分序列在开头部分质量差，所以根据这个图在做质控时选择两端都去低质量或只去3'末端。
Per tile sequence quality
新版本增加的功能，不太清楚是干啥的。
Per sequence quality scores
序列平均得分的数量。右侧越高越好，也就是大多数序列质量都得分在30以上。均值低于27(也就是错误率大于0.002)记为警告，均值低于20(错误率0.01)记为不合格。
Per base sequence content
显示碱基比例。正常情况下碱基比例应该差不多。AT与CG差异大于10%记为警告，大于20%记为不合格。
Per sequence GC content
每条序列GC含量百分比与模式化的正态分布GC含量相比较，超过15%记为警告，超过30%为不合格。
Per base N content
每个碱基位置N（未测到或不确定碱基）的比例
Sequence Length Distribution
各序列长度占比
Sequence Duplication Levels
重要序列重复级别
理想情况下所有序列应该是被随机打断了，测序后理论上不太会出现完全相同的序列。但由于PCR扩增或者污染等原因会造成一些重复序列，这里显示重复序列的比例。为了节省内存和时间，fastqc仅抽取了前100,000条序列，并只要超过75bp的序列就会被截断到50bp来分析。fastqc的说明文档对此进行了说明，结果不影响整体结果的反应。
所提供的图像有两条线来反应不同重复级别，蓝线表示整个序列集合中重复序列的分布，红线表示去除重复序列后的结果。这是新版本提供的功能，v0.10版本只有重复序列的程度。
一个好的结果应该是红线蓝线最左侧越大越好。通常会在红色线中看到一个高重复的峰，但同时在蓝色线上消失，这表明重复序列并不显著。如果在蓝色的线中有峰值，说明在大量不同的高度重复序列，这可能是污染或严重的技术重复。
如果非唯一序列数超过总数的50%就会被软件判断为不合格。
Overrepresented sequence
过表达序列。一般认为打断的序列只有很少部分会重复，如果一段序列重复达到总值的0.1%记为警告，超过1%记为不合格。

Trimming and Quality Filtering

根据结果去接头（adapter）、引物（Primary）尾巴（Poly-A）等。必须要去的是接头。常用的软件有cutadapt、trim_galore等等。一般用cutadapt，很多去接头软件的底层其实也是调用cutadapt。

Cutadapt

眼科中心服务器cutadapt 1.9.1版本安装在c0，c10节点上，需要提交到这两个节点才可以运行，否则很多节点用的是1.4.1，老版本的问题是功能有限，尤其是对于双端数据不支持（如-A参数）。cutadapt官网对于Illumina接头去除的说明如下：

If you have reads containing Illumina TruSeq adapters, follow these steps.

Single-end reads as well as the first reads of paired-end data need to be trimmed with A + the “TruSeq Indexed Adapter”. Use only the prefix of the adapter sequence that is common to all Indexed Adapter sequences:
cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -o trimmed.fastq.gz reads.fastq.gz
If you have paired-end data, trim also read 2 with the reverse complement of the “TruSeq Universal Adapter”. The full command-line looks as follows:
cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC  -o trimmed.1.fastq.gz -p trimmed.2.fastq.gz  reads.1.fastq.gz reads.2.fastq.gz
The adapter sequences can be found in the document Illumina TruSeq Adapters De-Mystified.

因此单端数据只需要用-a参数去掉“AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC”就可以了。

按照推荐我双端数据(Pair-End)的命令如下：

#PBS -N cutadapt-HBV 
#PBS -j oe 
#PBS -l nodes=c0:ppn=1 
#PBS -l walltime=5000:00:00 
#PBS -q low  

# set up some parameter

outputdir="/public/users/chentingting/Zoubin/HBV/QC"  

cd /public/users/chentingting/Zoubin/HBV  

cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC -q 30 -m 20 --trim-n -O 10 -o $outputdir/Sample_capsule-1.R1.trimmed.fastq.gz -p $outputdir/Sample_capsule-1.R2.trimmed.fastq.gz Sample_capsule-1.R1.fastq.gz Sample_capsule-1.R2.fastq.gz

其中的参数说明:
-a 序列 正向接头序列，单端测序只用这个。
-A 序列 反向接头序列，双端情况下设置。
-q 数字 表示最低质量值，在去接头前先将低于此数值的bases去除。如果只设置一个数值则从3'末端去除，如果用逗号分割两个数值则先去5'末端后去3'末端。一般设为30。

-q [5'CUTOFF,]3'CUTOFF, --quality-cutoff=[5'CUTOFF,]3'CUTOFF

Trim low-quality bases from 5' and/or 3' ends of each read before adapter removal. Applied to both reads if data is paired. If one value is given, only the 3' end is trimmed. If two comma-separated cutoffs are given, the 5' end is trimmed with the first cutoff, the 3' end with the second.

-m 数字 表示trim后最短bp低于此数的reads被抛弃，一般设为20。

-M 数字 表示长于此数字的reads被抛弃，默认值不限制。

--max-n=COUNT 抛弃有太多N的reads。COUNT如果设置为整数，就是按N的绝对个数来处理；如果设置为小数（0到1之间），就按每条reads中N的百分比来处理。

-O 数字 表示adapt和序列比对最少overlap的值，高于此值就认为是接头并修剪，默认是3，个人设置至少到5。

-o 目录 Read1的输出路径

-p 目录 Read2的输出路径

根据fastqc的报告，如果是RNA数据尾巴较多的情况，最好再去一次PolyA尾巴，少就不用了。

Trim Galore!

Trim Galore 合并了FastQC和Cutadapt到一个程序中。它的优势在于它可以根据FastQC分析的个体质量对每个reads进行修剪。同时可以设置程序对剪切后的序列用FastQC生成一个统计信息。对双端序列支持也很好。

选项

--phred33 使用ASCII+33质量得分作为Phred得分标准。默认选项
--fastqc 当剪切结束后用默认选项对结果文件进行fastqc分析
--fastqc_args "<ARGS>" 对fastQC设置参数。
--paired 设置双端序列
--dont_gzip 输出文件不压缩
--gzip 压缩输出文件，如果输入文件是压缩文件则自动压缩。
--length <INT> 设置低于数值长度的reads就抛弃掉，默认值20.
-q/--quality <INT> 切除质量得分低于设置值的序列。默认值20.
-a/--adapter <STRING> -a参数后为切除的接头序列
-a2/--adapter2 <STRING> 双端序列中read2所切除的接头序列，需配合'--paired'参数。
--rrbs 这是Trim Galore最独特的功能，也是我一开始找到使用这个软件的原因：特异性处理MspI所处理的RRBS样本数据(识别位点：CCGG)

示例：

trim_galore --phred33 --fastqc --fastqc_args "--noextract"  --paired --retain_unpaired    --dont_gzip  Sample_keratopathy-31R1.fastq.gz Sample_keratopathy-31R2.fastq.gz