标签归档:Sequencing

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
#####################################################################

GC偏好

GC偏好

测序中的GC偏好指的是基因组上GC含量在50%左右的区域更容易被测到,产生的reads更多,这些区域的覆盖度更高,在高GC或者低GC区域,不容易被测到,产生较少的reads,这些区域的覆盖度更少。用基因组单位长度的bin中的GC含量作为横坐标,覆盖度作为纵坐标作图,可以明显的看到该趋势。这种趋势在100kb为单位的bin中依然存在。

如图A中可以看出随着GC含量的增加,counts是先增加后减少,bin的大小为10kb。图C可以看出大部分片断的GC含量0.4到0.6之间。

GC偏好也存在其他地方,比如基因编码区内密码子的最后一位,C碱基往往占优势;基因的长度和GC含量成相关性;Aquifex aeolicus 的基因组整体GC含量是43%,而核糖体RNA操纵子的GC含量是65%。

如图,鸡(Gallus_gallus-5.0)基因组的GC含量与基因密度之间的散点图和拟合曲线,相关性非常明显。

影响

举个例子,1)在检测拷贝数的时候,GC含量低或者高的区域,其覆盖度小于GC含量中等的,但不意味着仅仅根据测序的覆盖度,就认为GC含量中等的拷贝数比高/低GC含量区域的高。2)在做RNA测序分析的时候,GC含量高/低的区域reads数少,并不一定说明这个基因的表达量低。3)在做基因组拼接的时候,因为GC偏好的存在,高/低GC含量的区域被测的少,这些区域的拼接难度就较大。
继续阅读

测序中加入Phix的作用

测序建库的时候,会加入一定比例的Phix,那么Phix文库有什么作用呢,我转了两篇文章,方便大家理解。Phix文库最主要的目的1)是调节碱基平衡,改善测序仪的空间校正,便于后期提高base calling的准确性,2)由于Phix序列已知基因组较小,在测序的过程中Illumina的测序仪就开始将测的read与phix基因组进行比较,预估测序指标。我也遇到过,Illumina工程师在维护测序仪时,用Phix文库测试。转载内容详见下文

继续阅读

bismark DNA甲基化测序比对-bisulfite-seq

bismark调用bowtie2进行比对,调用samtools生成bam文件,因此在运行bismark之前,需要安装bowtie2和samtools
请注意,fastq文件要进行质控,比如去掉低质量的reads,去掉adaptor等,可以看本文最下方推荐的PPT,本文不介绍,此外本本只介绍到序列比对,后续的统计分析没有介绍,有兴趣的朋友可以关注swDMR和methykit工具包。

安装bismark

wget http://www.bioinformatics.babraham.ac.uk/projects/bismark/bismark_v0.16.1.tar.gz
tar -xvzf bismark_v0.16.1.tar.gz 

生成转换后的基因组

# --path_to_bowtie后面跟的是文件夹
# --verbose 输出log信息
# ./ref 文件夹中有一个基因组fasta文件
# --bowtie2指明用的是bowtie2
./bismark_v0.16.1/bismark_genome_preparation --path_to_bowtie /home/zzx/bowtie2-2.2.9/ --bowtie2 --verbose ./ref/

因为在重亚硫酸盐甲基化测序中,因为未甲基化的C会变为T,在正链表现为C–>T,但是在负链有C变为T,转换为正链时,即为G–>A,所以基因组需要进行两种转化,才能用于比对。在基因组目录下产生Bisulfite_Genome目录,有CT_conversion和GA_conversion文件夹,这两个文件夹包含转换后的fasta文件和bowtie2建立的索引bt2文件。

fastq中的BS转换后的read与转换的参考基因组比对,得到在参考基因组中的位置,再与原始的参考基因组比较,确定methylate call

bismark比对

--non_bs_mm     Optionally outputs an extra column specifying the number of non-bisulfite mismatches a read during the alignment step. 
--nucleotide_coverage   该选项 会生成nucleotide_stats.txt文件。creates a '...nucleotide_stats.txt' file that is also automatically detected by bismark2report and incorporated into the HTML report.
--multicore   Sets the number of parallel instances of Bismark to be run concurrently.  一个实例已经将1个核分配给bismark,两个或者4个核分配给bowtie2,一个分给samtools
./bismark_v0.16.1/bismark --genome_folder ./ref --bowtie2 --nucleotide_coverage --non_bs_mm --basename test --temp_dir ./tmp --samtools_path ./samtools-1.2/ --path_to_bowtie ./bowtie2-2.2.9/ -1 1.fastq -2 2.fastq --output_dir ./

继续阅读

转录组分析新工具流程–HISAT2-stringtie-ballgown

一,HISAT2(Hierarchical Indexing for Spliced Alignment of Transcripts2)

HISAT2是一个对比对RNA-seq reads的快速灵敏的spliced alignment工具,HISAT2支持DNA和RNA比对。针对reads覆盖多个外显子,HISAT其包含两种索引:1,global FM索引,代表整个基因组,2,许许多多的local FM索引,每个索引代表~56,000bp,~55,000个local索引覆盖整个基因组。HISAT基于Bowtie2来处理大多数FM索引。心的索引scheme叫做 Hierarchical Graph FM index (HGFM)。HISAT也支持indel和paire end模式,并且支持多线程和SRA Toolkit。HISAT2官网提到HISAT2是HISAT和TopHat2继承,建议大家从HISAT和TopHat2迁移到HISAT2上来。

HISAT2的优势如图(图片来自HISAT2文章)

继续阅读