目前用于新一代的测序的主要仪器有Illumina/Solexa的Genome Analyzer、ABI的Solid和Roche的454,它们都能高通量的测序,产生大量的测序结果,接下来就要对序列进行拼接,用于拼接的软件也有很多,比如velvet、soap、abyss、maq等,454的还有专门的newbler。平时用velvet比较多,就简单介绍一下。
velvet对短序列的拼接效果比较好,所以多用于对Illumina等产生的短序列片段进行组装拼接。下面以Illumina的GAII产生的结果为例进行说明。
一、单端测序
单端测序可以直接对fastq格式的原始文件进行处理,首先是用velveth命令建立hash表子集
输入./velveth会出来使用帮助:
Usage:
./velveth directory hash_length {[-file_format][-read_type] filename} [options]
directory : directory name for output files
hash_length : odd integer (if even, it will be decremented) <= 75 (if above, will be reduced)
filename : path to sequence file or – for standard input
File format options:
-fasta
-fastq
-fasta.gz
-fastq.gz
-eland
-gerald
Read type options:
-short
-shortPaired
-short2
-shortPaired2
-long
-longPaired
Options:
-strand_specific : for strand specific transcriptome sequencing data (default: off)
Output:
directory/Roadmaps
directory/Sequences
这一步主要是要确定使用的hash值,hash值必须为奇数,且小于MAXKMERLENGTH,这个值默认为31,但是在安装的时候可以调整。具体的命令可以是:
velveth result 29 -fastq -short sequence.fastq
建立hash子集后用velvetg命令进行组装,输入./velvetg也会出来帮助:
Usage:
./velvetg directory [options]directory : working directory name
Standard options:
-cov_cutoff <floating-point|auto> : removal of low coverage nodes AFTER tour bus or allow the system to infer it
(default: no removal)
-ins_length <integer> : expected distance between two paired end reads (default: no read pairing)
-read_trkg <yes|no> : tracking of short read positions in assembly (default: no tracking)
-min_contig_lgth <integer> : minimum contig length exported to contigs.fa file (default: hash length * 2)
-amos_file <yes|no> : export assembly to AMOS file (default: no export)
-exp_cov <floating point|auto> : expected coverage of unique regions or allow the system to infer it
(default: no long or paired-end read resolution)Advanced options:
-ins_length2 <integer> : expected distance between two paired-end reads in the second short-read dataset (default: no read pairing)
-ins_length_long <integer> : expected distance between two long paired-end reads (default: no read pairing)
-ins_length*_sd <integer> : est. standard deviation of respective dataset (default: 10% of corresponding length)
[replace ‘*’ by nothing, ‘2’ or ‘_long’ as necessary]
-scaffolding <yes|no> : scaffolding of contigs used paired end information (default: on)
-max_branch_length <integer> : maximum length in base pair of bubble (default: 100)
-max_divergence <floating-point>: maximum divergence rate between two branches in a bubble (default: 0.2)
-max_gap_count <integer> : maximum number of gaps allowed in the alignment of the two branches of a bubble (default: 3)
-min_pair_count <integer> : minimum number of paired end connections to justify the scaffolding of two long contigs (default: 10)
-max_coverage <floating point> : removal of high coverage nodes AFTER tour bus (default: no removal)
-long_mult_cutoff <int> : minimum number of long reads required to merge contigs (default: 2)
-unused_reads <yes|no> : export unused reads in UnusedReads.fa file (default: no)Output:
directory/contigs.fa : fasta file of contigs longer than twice hash length
directory/stats.txt : stats file (tab-spaced) useful for determining appropriate coverage cutoff
directory/LastGraph : special formatted file with all the information on the final graph
directory/velvet_asm.afg : (if requested) AMOS compatible assembly file
这里主要用到的参数是-cov_cutoff(覆盖度)这个参数,其他都可不加,其实覆盖度参数也可以省略,但是最好还是加上,增加可靠性。还有-exp_cov这个参数加上后会对重复区域进行处理,也最好加上,后面的数值用auto就可以了。其他的一些参数就根据你自己的喜好来设置。具体的命令可以是:
velvetg result -cov_cutoff 30 -exp_cov auto -min_contig_lgth 100
最后根据出来的n50和max contig长度来判断拼接的效果,为了达到最好的拼接效果,一般要对hash值和覆盖度进行一系列的设置来进行比较。
二、双端测序
双端测序与单端测序相比在运行velveth时将-short参数改为-shortPaired,其他一样。
在运行velvetg时可加入-ins_length 和-ins_length_sd 参数提高拼接的准确性。具体命令可以是:
velvetg result -cov_cutoff 30 -ins_length 300 -ins_length_sd 100 -exp_cov auto -min_contig_lgth 100
其实运行这几个命令很简单,最重要的还是根据你的测序数据选择合适的参数。
你也是学生物信息的??
是的~
您好,有velvet拼接上的问题,能向您请教一下吗?
@insidi:
有问题你可以在这里提出
好的!illumina paired-end 2×100测序,120x genomic coverage!参数该如何设置呢?已使用您上面的参数命令,结果很不理想!short or shortPaired ?exp_cov 该设为120?盼为解惑~
@insidi:
100bp的测序长度最好先过滤一下低质量的碱基,先用shuffleSequences_fastq.pl脚本把两端的reads放到一个文件中,双端测序用shortpaired,exp_cov用auto参数
以前都是用soap现在想试试这个软件请问一下,如果想PE和MP一起拼接的话命令是写成velveth result 29 -fastq -shortPaired file1.fq file2.fq -shortPaired2 file3.fq file4.fq 吗?
@mt: 可以这样的,使用velvetg的时候插入长度也要输入两个
哦,谢谢,还有个问题,如果我用correction来优化序列,然后用merge_pair.pl把能有pair关系的序列输出到merge.pair里,里面的格式是这样的
read1
sequence
read2
sequence
read1
sequence
……
使用上面那种fasta拼接或者按read标签把corr文件里面的质量文件取出来的fastq拼接你觉得那种拼接方式更好啊?能不能用soapdenovo拼出一个效果最好的然后用你新写的那个软件minimus2拼到一起呢?
看你写这个博客很久了,知道您比较资深,想听听您的建议,谢谢
@mt: 你说的correction就是进行低质量过滤吧?这个你可以用我另一篇文章中介绍的这个工具来进行,velvet在拼接的时候不管序列质量的,所以fasta和fastq的结果都是一样的。soap我只试用过一下,感觉拼接效果没有velvet好,所以我不用的。你是要用minimus2把哪几个并到一起?
我想把相同的数据velvet做一次初步拼接,然后soap做一次拼接,然后把两个的合到一起。不知道会不会好一些。
我这两天在做一个真菌的线粒体,因为没有单独测线粒体,所以我将全基因组的454的数据和solexa的数据map到近源菌的线粒体上,然后用最严格的要求取出数据,用454和soap分别拼接,用minimus2合并到一起,出现了序列重复使用情况,有几K的数据重复出现在两个scaffold里面。
@mt: 在程序进行序列拼接过程中出现误差或者错误都是有可能出现的,如果你确定这几K的序列不是重复序列的话,那就可能是拼接过程中有问题,用不同的参数和不同软件都试一下,找出最佳拼接吧~
嗯,你有文章介绍AMOS的其他功能吗?我只是看了你的博客才学会的minimus的用法,听说AMOS功能挺强大的~
@mt: AMOS是有很多功能,我也是看它们网站上的文档的,没看过其他的
不知道你用过什么重测序的拼接软件吗?我找到一个叫maq的,还在看说明,你用过吗?感觉怎么样啊?
@mt: maq对于重测序的比对用得还是比较多的,后面的格式转换之类的也比较方便。另外推荐bowtie和samtools组合
哦,用过samtools的很少的功能,bowtie还没用过,我就想如果能够有那种有公布过的序列拼接的时候,能把序列map出一个框架,这样比较好拼接,但是我现在接触的几个软件,只能单独map或者是denovo没有直接map拼接的,感觉对于已公布的物种再用denovo拼比较费时间~
@mt: 看你物种的变异程度了,如果变异慢的话可以用reference先map一下,如果变异快的话还是要做de novo,不然错误太高。velvet有columbus模块可以map reference的,不过我也没试过
哦,我看maq的介绍说通过把reads map到reference上可以找到插入缺失,我实验一下这几软件吧,重测序好多都是同一株,可能打入了些新的插入缺失,主要想做转录组可能需要重测序,
用新软件有时候可挠头了,有时候他写个参数什么的缩写都不知道是干什么
@mt: 它说的是单碱基的插入缺失吧,maq可以用来找snp的。用新软件之前都要看文档,没办法的
英语学的不好,所以看文档有点费力,但是要是都讲到了我可以慢慢看,但是有的参数就泛泛的解释了一下字面意思而且有的不理解的单词查翻译软件夜查不出来,也不讲有什么作用,这是让我比较头疼的地方。
@mt: 有些东西就是要靠慢慢积累的,你是研究生?
不是啊,我是本科毕业,跨专业就业到生物信息学的,所以很多东西得学~~我是从什么都不懂开始学的~~
请问你一般怎么处理的病毒的长末端重复区啊?
@mt: 没处理过病毒的序列~
能给发一个velvet吗
@xia 最新版下载地址:http://www.ebi.ac.uk/~zerbino/velvet/velvet_1.1.04.tgz
我刚用velvet 1.0.18的,对于velveth中的参数hash_length有点疑问,帮助提示可以输入“m,M,s”格式的参数值,但当我输入的是“21,31,2”时,它只给我生成21到29和之间的结果,而且除了21,其它的结果都是链接21的。当我输入“21,32,2”时,却也只生成31的结果。
@lsq 我没有用过这种设置方式,你可以用contrib里面的VelvetOptimiser试试~
请问细胞器基因组和和基因组做数据拼接组装的时候,参数设置有什么不一样吗?
@rachel 这个没啥区别的~
有个小问题:
Pair-End在500bp以下时,是forward-reverse,超过2k,是forward-forward
比如说:
-ins_length 500
-ins_length2 2000
但velvet是默认forward-reverse的,是不识别后者的。因此需要先将2k插入文库其中一条read反向互补,
问题是,velvet说明书中是转第一个文库,可我认为应该转第二个文库。因为第一个文库是forward,不需要转换,第二个文库也是forward,需要转换成reverse。
根据illumina那张Pair-End测序原理的图示。
比如说有AT……CG这样一条序列,互补链是CG……AT,forward-reverse测得AT……和CG……,或者CG……和AT……
forward-forward测得AT……和……CG,或者CG……和……AT
现在要将forward-forward的其中一条read反向互补,怎么看都应该转第二个read,才是啊。可是velvet说明书中在shuffleSequence中不是这样写的。就困惑了。
除非,forward-reverse测得是……AT和……CG,可是接头序列在read中如果出现的话一般是在右侧的,所以read应该不是adaptor……AT,adaptor……CG的。
@tang 你说的超过2K的insert长度应该是指用mate pair测序得到的结果吧,在mate pair的结果当中,两个read的方向是<=insert=>的,所以为了让它们变成=>insert<=,需要把两个read全部取反向互补序列。
嗯,mate-pair的方向确实是outward,PE是inward。将read align到scaffold上,在同一个scaffold上的pair read,基本上是一正一负的(只是随机查看了少数几个PE),与mate-pair的图示一致。所以我很奇怪,long PE是怎么得到forward-forward的。假如是forward-forward,原理上,应该反转3‘端的read,这类read可能在文库1,也可能在文库2;本来,mapping的结果,预期会得到都是正、或者都是负的pair End,某个文库的read都在下游,那么就转这个文库就好了。可是现在发现,其实2k的文库应该是mate-pair的。
另外,在SOAP拼接时,将2k的插入文库(outward),seq_reverse设为0,拼接质量大大下降,所以还是有影响的。
现在,我想再分别拼三次,将第一个文库反向互补reverse complement,之后拼接,将第二个文库RC之后拼接,两个文库都RC之后拼接,看看拼接质量会不会变化。也许,velvet根本不需要多此一举,不管是inward还是outward,都是同样处理的。
velvet拼接时,确实需要将两个MP文库都进行反向互补,拼接质量有较为明显的上升。
谢谢!
@tang 感谢您的测试反馈!
我在学习用VelvetOptimiser,这一款你用的怎么样,和Velvet相比哪个更好呢?
@小花ingfyq velvetoptimiser只是对velvet的参数进行最优化选择,用的还是velvet进行拼接
Velvet用得还顺手,就是在用VelvetOptimiser时,也不知道是不是我下的版本的问题,
指令是perl ./VelvetOptimiser.pl -s 21 -e 31 -f ‘-shortPaired -fastq ../d.fastaq ‘ -t 8
我观察了一下运行过程,发现程序自动生成auto_data_21到auto_data_31六个文件,到最后删选过后会剩下auto_data_31一个文件,但是运行到最后,auto_data_31文件里面的内容,例如Sequences文件等被运行结果给覆盖了,,这个文件夹还是存在的。所以会报错说can’t open Sequences.知道这个问题的原因吗?
@小花ingfyq 我都是自己用velvet跑一遍所有kmer自己挑选最佳结果的,不用optimiser
想请教一下,我现在想自己写一个程序,选出最佳的Kmer和cut_off的值,我看了Manual,Kmer的选择范围是10~31,cut_off的值是根据exp_cov的值确定的,所以现在我觉得有点一头雾水,我还是个新手,能帮我解答一下不?
@小花ingfyq Kmer的范围没有那么窄,最大的要根据你的reads长度选择,cut_off和exp_cov的值都可以用auto参数,也可以设定cut_off的值,exp_cov还是auto就好了
我现在每个read的长度为97,那我最好设置为多少呢?
@小花ingfyq 这个最好是要试出来的,一般范围可能在41-59之间,跑的时候再往两边多选几个参数看看
我按照你说的,我分别设Kmer长度为21,31,41,45,51,55,exp_cov设置为auto,我想选取最佳的cov_cutoff的值,用R软件
> weighted.hist(data$lgth,data$short1_cov, breaks=0:50)描述出关系,发现这个组装的效果一点都不好,45之前的取值所得的图形的峰值在1左右,很小,而51和55得到stats.txt的数据少,像55的话就只得到5行的数据,画出的图的效果很不好,在2附近只有一根柱(我画的是柱状图)。我很困惑,请再指点,麻烦了
@小花ingfyq 我从来不看stats里面的数据,直接统计每个contigs.fa的各项数据,找N50最大的
您好,我想问一下,我现在有reads1.fq, reads2.fq两个文件,一个是forwad, 一个是reverse, 请问如何操作呢?
@Zheng Zhou 你是paired-end还是mate pair?paired的直接用velvet就好了,mate pair的需要将序列掉个头
怎么根据质量值过滤一些误差比较大的reads啊?
@zyy 你可以使用fastx_toolkit这个工具进行质量过滤
楼主你好!
我用Bowtie拼接大约5万对88bp的双端测试数据,请问对比参数因该如何选?我直接用其中一条的数据build可以么?
还有对边比后的文件应该如何看?我打开发现是可以配对的序列,但我想要拼接之后的长链
@Zorro Van bowtie是一个比对软件,你为什么用它来拼接呢?它比对的对象应该是已知的基因组序列
楼主你好!
请问拼接结果如何阅读?怎样看拼接结果的好坏?楼主在velvet基础上用过Oases拼接转录组么?
@helloJohn 拼接结果好坏通过拼接长度,N50,contig数目等来判断。还没用过Oases,RNA拼接貌似Trinity用的人更多吧
这些结果在哪里看呢,N50 这些?velvet的结果主要是看contig.fa么。
@helloJohn 这个需要用脚本统计得到的contig.fa结果