标签归档:FASTQ

什么是ubam文件,为什么ubam文件比fastq文件好

ubam是Unmapped BAM Format,是BAM文件的一个变种,里面的read是未经map的。大部分测序供应FASTQ文件,这是最常见的测序分析的起始文件。FASTQ文件的优势是,压缩比比bam文件好,解压速度快。

但与ubam文件相比,FASTQ并不是最理想的:

1)单个文件更容易分析

在双端测序中,有些软件希望配对的reads放在一个文件中,有些希望配对的文件,有些软件直接根据read在文件中的位置判断read是否配对,当然现在FASTQ通过在文件中加入/1和/2来表明配对的read解决这一问题。在生信分析的时候,一个FASTQ文件往往和另一个FASTQ文件关联配对,比如R1和R2,往往会花费更多时间来验证read是否配对。但如果通过bam文件的话,会更简单。只需要在FLAG这个地方加入对应的值,比如77和141,就能指定。单个文件更加简单,能储存更多metadata的信息。

2)BAM文件可以储存更多信息

实际上,正是这个原因,我才开始用ubam文件的。
因为FASTQ文件格式的问题,很多其他信息不能存在FASTQ文件中,而BAM/SAM文件可以储存很多信息,只需加tag或者header就行。如果你想放,可以放用了什么软件,平台,样本标签等。
上图是不利用ubam文件的流程,下图是利用ubam文件的流程。

这个我体会很深,比如我在分析带分子标签UMI的数据时,需要在bam文件中添加read的UMI标签序列,具有相同的UMI标签的序列说明来自同一个分子。如果我将带UMI标签的reads直接比对基因组,UMI标签序列会影响比对质量,而且后续对bam文件的处理提取UMI时,需要考虑CIGAR标签,这个时候就头大了。实际上,我是先将fastq转成ubam文件(picard的FastqToSam工具),然后提取ubam文件中的UMI,生成新的ubam,将ubam转成fastq文件(reads都放在一个文件fastq文件中,bwa能识别其中的配对reads),与参考基因组进行比对。这个时候问题就来了,比对之后生成的bam文件,不包含我提取的UMI标签序列信息,因为这个信息在ubam文件中。不要着急,picard的另外一个工具MergeBamAlignment,可以将ubam文件和bam文件合并,这样ubam文件中的tag信息,比如UMI标签信息,就会补充到新的bam文件中。这样生成的bam文件供下一步数据分析。

3)BAM能实现随机读写

BAM文件采用的是BGZF压缩,可以通过随机读写快速读取文件特定位置的数据。对于未比对的reads(unaligned reads),这个功能不是很吸引人,因为后续序列比对都是顺序比对的,不可能跳着比对。

参考:
https://blastedbio.blogspot.com/2011/10/fastq-must-die-long-live-sambam.html
https://gatkforums.broadinstitute.org/gatk/discussion/5990/what-is-ubam-and-why-is-it-better-than-fastq-for-storing-unmapped-sequence-data
https://gatkforums.broadinstitute.org/gatk/discussion/11694/why-is-converting-from-fastq-to-ubam-nesessary-before-preprocessing
https://gatkforums.broadinstitute.org/gatk/discussion/6484/how-to-generate-an-unmapped-bam-from-fastq-or-aligned-bam

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

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
继续阅读

文件介绍–FASTQ文件格式

在培训部门同事的时候,发现刚开始学生信的人,只是在学如何运行命令,但对自己手头的文件格式和内容却不了解,这对分析的流程的深入理解和研究是非常不好的,所以刚学习的人,应该在等待分析结果的时候,多去了解下文件的内容,程序的大体算法等,这对以后的工作优化是非常有好处的。本文简单介绍一下Fastq的文件格式,希望新手多查文档,多了解自己接触的东西。

文库构建和测序

DNA分子会通过超声波或者酶被打断成几百碱基的小片段,然后在小片段DNA分子的两端添加接头,便于测序和样本区分。当然现在也有转座酶技术,通过转座酶同时实现DNA片断化和加接头和引物的过程。将文库上机,文库中的DNA分子首先与flowcell上lane中的接头结合,通过桥式PCR进行扩增(cluster簇增长),待达到一定量之后,进行便合成边测序。

测序数据的产生

每进行一个cycle,测序仪会合成过程中产生的荧光进行拍照分析,将产生的数据记录在bcl(base calling)文件,测序仪后续会将bcl文件转化成fastq文件,并同时进行demultiplexing。Demultiplexing会将属于同一Index的reads放在同一个fastq文件中,就个过程是Illumina测序仪自动化进行的,如果涉及更复杂参数的demultiplexing,可以下载illumima的bcl2fastq程序,指定相关文件夹即可自己进行demultiplexing。因为在测序上机之前,会在测序仪中设置样本对应的index,这样通过demultiplexing,测序仪对每个样本产生一对fastq文件。Illumina测序仪会自动将拍照形成的bcl文件转换成FASTQ文件。此文件为数据分析的开始。

FASTQ文件格式简介

Paired End(PE)测序,会生成的一对FASTQ文件,分别为R1和R2,正是因为双端测序的产生,使得序列拼接和比对更加准确,因为如果单端的话,只能依靠单端read的长短信息,如果双端的话,能依靠整个插入DNA片断长度的信息,通过判断配对序列的相近距离,确认reads的位置,减少了read比对到多个地方的情况发生。
配对的reads分别在两个FASTQ文件中,其中一个FASTQ文件的内容如下,每四行表示一条read:

继续阅读

BCL文件与BCL2FAFSTQ程序简介

BCL文件

测序产生的原始文件是BCL(binary base call)文件,测序仪在测序的时候,每个cycle都会测量编码不同颜色的通道强度,并确定最有可能的碱基类型。Real Time Analysis (RTA) 软件会将碱基类型和可信度(一个质量分数)。与FASTQ文件不同的是,BCL文件是实时产生,每个cycle的每个tile都会有一个对应文件,文件放在

<run directory>/Data/Intensities/BaseCalls/L<lane>/C<cycle>.1

文件的命名

s_<lane>_<tile>.bcl

bcl2fastq

该文件需要通过Illunima的软件或者第三方分析工具将BCL文件转成FASTQ文件。一般而言,数据下机之后,Illumina测序仪会自动将BCL转成FASTQ文件。有时候,根据实验需要,需要自己手工将BCL文件转成FASTQ文件,比如自己设计的index中含有简并碱基,或者需要调整一下转换的参数等。

Illumina提供bcl2fastq的程序包,共离线处理BCL文件,生成FASTQ文件。

bcl2fastq  -i /paht/to/run/Data/Intensities/BaseCalls/ \
       -o /output/dir --sample-sheet /paht/to/run/SampleSheet.csv \
       -R /paht/to/run/

bcl2fastq文件有很多参数可调,比如在FASTQ中记录read的index(fastq文件中会记录配对的index具体序列,以及会生成额外对应的index文件),可以添加–create-fastq-for-index-reads选项。

如果允许index有mismatch的话,可以通过–barcode-mismatches设置。

–fastq-compression-level可以设置生成的FASTQ.GZ文件的压缩比例

安装如下

unzip bcl2fastq2-v2.17.1.14.tar.zip
tar -xvzf bcl2fastq2-v2.17.1.14.tar.gz
./bcl2fastq/src/configure --prefix=/path/to/install/
make
make install

安装过程中,如果遇到一下问题,请更新gcc版本
问题1
cc1plus: error: unrecognized command line option “-std=c++11”
问题2
undefined reference to `boost::re_detail::perl_matcher
collect2: error: ld returned 1 exit status

软件选项说明如下:
继续阅读