标签归档:QC

fastq压缩之后的gzip文件大小与样本数据量

在测序的时候,我们先拿到的是样本fastq压缩后的gzip文件,这个时候可能最关心的是数据量够不够,那么fastq.gz文件大小和测序数据量有什么关系呢。
我用Miseq测序数据(gz文件200M左右),Hiseq panel(gz文件50M左右)和WES测序数据(gz文件4G左右)进行了简单的分析。有意思的地方是,虽然R1和R2的数据量是一样的,解压出来的文件大小是一样的,但R2的gzip文件总比R1大。不管是Miseq还是Hiseq的panel测序,压缩后的R2均大于R1文件,且文件越小,差异越大。

因为R1和R2数据量相同的原因,我只看R1的真实文件和gz文件大小与数据量之间的关系。
数据量=FASTQ文件行数/4*151/1000/1000 单位为M
真实文件大小估计=FASTQ文件行数/4*357/1024/1024 单位为M,预测值,差别不大,因为FASTQ文件中每四行357个字符(和平台和设置有关系),每个字符1byte。
GZ文件大小通过ll -h查看
因为FASTQ文件是规范的,每四行字符基本一致,所以FASTQ真实文件大小和数据量成正比。比如我前面提到的每四行有357个字符,其中序列只占151个字符,也就是说FASTQ文件大小大概是测序量的357/151≈2.3倍多。但因为FASTQ文件为文本文件,占用空间较大,所以一般将FASTQ文件压缩成gzip格式文件。

用gzip文件大小除以估计的真实文件大小,得到压缩比,发现压缩比和测序平台有关系,Miseq的压缩比在0.35左右,Hiseq平台中的panel测序在0.24左右,WES在0.2左右。可以发现WES的数据的压缩效率最高,猜测不同平台在压缩过程中设置的压缩比例不同。Hiseq测序平台的gzip文件压缩性能大于Miseq平台。

以Miseq的gzip文件大小和数据量作图,可以很直观的发现gzip文件大小和数据量之间的线性关系。如图,斜率是1.245,截距可以忽略不计,因为纵坐标很大。也就是说R1文件的数据量约等于gzip文件的1.245倍。

从另外一个角度计算,R1文件的数据量=gzip文件大小/压缩效率/2.3,则因为Miseq的压缩效率在0.35,Hiseq在0.23左右,所以Miseq R1文件的数据量大概是gzip文件的1.242倍(1/0.35/2.3),和上图拟合的结果很相近,Hiseq约在1.89(1/0.23/2.3)倍。真实的数据量其实是R1文件数据量的2倍(R2文件中的数据量和R1文件相同)。

结论:不同平台和不同测序方法之间的fastq文件压缩比例并不一致,但同一种方法和平台内的压缩比例是相近的,因而可以根据gzip文件的大小推测出测序的数据量。
我手头的数据中,Miseq的测序量约是R1 gzip文件大小的2.49倍,Hiseq的测序量约是R1 gzip文件大小的3.78倍。不同平台和方法之间,这个值是不固定的,可以根据手头的数据计算一下这个比例,就能迅速的估算出样本的数据量。

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

测序数据的预处理

测序得到的原始测序序列,里面含有低质量的reads。低质量的reads可能因为flowcell上的cluters不是有单一DNA扩增而来,或者几个cluters混为一起等。还有测序仪在前几个和几个cycle测序质量不好,需要关注一条read的前后几个碱基的质量。如果质量非常不好,测出来的碱基可能为N(无法确定碱基类型)。

此外,原始reads中还包含测序接头等序列。如果一个文库的平常插入长度为450bp的话,不一定每个插入长度都为450bp,如果个别分子插入长度为100bp,双端配对150bp测序的,会将该片断测穿,配对的reads会多包含50bp的index或者SP等序列。为了保证信息分析质量,需要对下机的raw reads 进行精细过滤,得到clean reads,后续分析都基于clean reads进行。

为了提高下一步的比对质量,此时数据预处理的过程主要包括:
• 去掉接头,去掉开头和结尾几个碱基中质量不好的碱基
• 滑窗扫描,检查是否有好几个连续碱基质量不好的情况
• 丢弃过短的read
• 去接头
• 去掉前端碱基质量低于一定值的碱基
• 去掉后端碱基质量低于一定值的碱基
• 以4bp为窗口滑窗扫描read,如果4个碱基平均质量低于15,则截断
• 丢弃序列长度小于36bp的reads
继续阅读

Asymmetric trimmomatic output with paired-end sequencing reads

运行trimmomatic的默认参数

java -jar trimmomatic.jar PE -threads 16 -phred33 sample_1.fastq sample_2.fastq sample_trimmed_paired_1.fastq.gz sample_trimmed_unpaired_1.fastq.gz sample_trimmed_paired_2.fastq.gz sample_trimmed_unpaired_2.fastq.gz ILLUMINACLIP:adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

输出文件中

sample_trimmed_unpaired_1.fastq.gz是sample_trimmed_unpaired_2.fastq.gz的十多倍,异常的大。1文件(forward reads)中unpaired reads非常多,显著多余2文件(reverse reads)中的unpaired reads。

原因是

After read-though has been detected by palindrome mode, and the adapter sequence removed, the reverse read contains the same sequence information as the forward read, albeit in reverse complement. For this reason, the default behaviour is to entirely drop the reverse read.

reverse read被丢弃之后,导致forward read中unpaired read非常多,出现上述非对称结果。

ILLUMINCLIP的最后两个参数可以额外的,其中minAdpterLength的默认参数为8,keepBothReads的参数为FALSE
ILLUMINACLIP::::::

设置这两个额外的参数后,输出文件中的unpaired 文件大小趋于正常。

By specifying “true” for this parameter, the reverse read will also be retained, which may be useful e.g. if the downstream tools cannot handle a combination of paired and unpaired reads.

PS:

trimmmomatic默认认为下游处理工具可以处理paired和unpaired reads,所以在palindrome mode(回文模式?)中丢弃了相关reads,而日常工作中,我们常常只利用paired-end read,所以keepBothReads参数最好设置为TRUE。优化后的参数见下:

java -jar trimmomatic.jar PE -threads 16 -phred33 sample_1.fastq sample_2.fastq sample_trimmed_paired_1.fastq.gz sample_trimmed_unpaired_1.fastq.gz sample_trimmed_paired_2.fastq.gz sample_trimmed_unpaired_2.fastq.gz ILLUMINACLIP:adapters/TruSeq3-PE-2.fa:2:30:10:8:TRUE LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

REF

http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf

http://seqanswers.com/forums/showthread.php?t=38068

http://www.usadellab.org/cms/?page=trimmomatic

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
####################################################################

高通量质控统计软件包Rqc–Quality Control Tool for High-Throughput Sequencing Data

在得到下机数据之后,首先会对下机数据进行质控,并对碱基质量碱基分布等进行统计和绘图,得到质控报告。大家常用的质控软件有Faxtx tookit,fastqc,Trimmomatic 等,但出图要么是自己根据统计结果自己画图要么是用fastqc的图。个人觉得fastqc的图有点土,于是找到了Rqc包,顾名思义R quality control。发现国内没有人介绍该包,在此向广大人民群众推荐该包。

Rqc包的介绍网址在http://www.bioconductor.org/packages/release/bioc/html/Rqc.html

有两种安装方式,
一种是通过bioconductor
source("http://bioconductor.org/biocLite.R")
# 如果下载包的网速非常慢,请更换镜像chooseBioCmirror()
biocLite("Rqc")

另一种方式是通过github
install.packages("devtools")
library(devtools)
install_github("labbcb/Rqc")

Rqc包自带测试数据用于测试Rqc:
library(Rqc)
folder rqc(path = folder, pattern = “.fastq.gz”)

结果输出如下:

read质量分布Per Read Mean Quality Distribution of Files

继续阅读