【2.3】序列拼接-velvet

一、安装

1.下载地址:http://www.ebi.ac.uk/~zerbino/velvet/

2.首先需要安装连个关联软件包zlib1g-dev 和pdflatex

sudo apt-get install zlib1g-dev (我的系统自带这个了)
sudo apt-get install texlive      (这个得下)

3,解压下载的velvet软件包,

tar xvzf  velvet_version.tgz
cd velvet_version
 make color 'CATEGORIES=2' 'MAXKMERLENGTH=127' 'BIGASSEMBLY=1''LONGSEQUENCES=1' 'OPENMP=1'
 (参数怎么修改见下面)
make clean 就是删除,然后再重新安装

其他一些设置:

make color  设置不同的颜色???
make 'MAXKMERLENGTH=127'  #默认的最大kmer为31,这个需要修改
make 'CATEGORIES=57'   #处理的数据的管道数,不同的管道可以用来处理来自不同插入文库或者不同样品(这里面有一个问题:宏基因组算是多个样品吗?),而这个值越大,对内存消耗越大,而我的样品,默认的为2,我怕内存不够,暂且设为2吧
     make 'BIGASSEMBLY=1'    #read的ID保存在一个32字节的字符串中,如果你有大于22亿(2.2G) 的reads,你就需要更多的内存,设置为1,是调用更多的内存。
     make 'LONGSEQUENCES=1'  #read的长度保存于16字节的字符串中,如果长度大于32kb,就意味着要更多的字节来存放字长度,所以可以直接设置为1,调用更多的内存,但是跟上面一样也会消耗很多的内存。
   make 'OPENMP=1'  #可以利用多个cpu,但是速度不会有指数的变化,因为只有一部分可以利用多线程
    make ‘BUNDLEDZLIB=1’  #如果不能使用系统默认的zlib,也可以造常运行。
也可以直接这样,一致性到位:    make color 'CATEGORIES=57' 'MAXKMERLENGTH=127' 'BIGASSEMBLY=1' 'LONGSEQUENCES=1' 'OPENMP=1'

环境变量:

export PATH=$PATH:/sam/velvet/bin

#注意我把安装好后的velveth,velvetg两个程序罗到了新建文件夹/bin里面,这样调用这两个程序的时候,就直接可以在终端的任何地方输入velveth,而不用cd 到程序目录,然后./velveth罗。

也可以这样来改环境变量:

sudo cp velvetg /usr/bin/
sudo cp velveth /usr/bin/

运行环境:

64bit,gcc(GNU Compiler Collection,GNU编译器套装,是自由的类Unix及苹果计算机Mac OS X 操作系统的标准编译器),>12G内存

二、主要工作原理图:

1、Bruijn算法

给定的数据总是存在着冗余信息,因此,算法的使用很重要。当前大多数是用的overlap-consensus-layout算法,它的每一个read都是单独的实体。而Bruijn算法,将分析基于observed words(或者k-mers)。给定k-mers,不管发现多少次,都是唯一的node。此外,他可以无差别的容纳long和short read的混合物。 构建Bruijn图关键在于对所有的reads进行hash。Hash值通过设置参数给定,或者是默认的21。这对所有的序列进行成对比对时会减少时间,一旦所有的reads已经被hash,由k-mer node可以追溯他们的每一个路径。增值覆盖率和创建arcs就属于这种方法。

2、错误的移除

错误移除算法:Tour Bus,对没有错误连接的graph进行移除。这要确保基因组中低覆盖率的唯一的点没有被破坏。因此,这一步在画图之后进行。 这里有两种类型的错误:tips和bubbles。Tips:起源于低质量的read的结尾部分没有任何的重叠。Bubbles:他产生于长read的中间或两个错误的read结尾部分重叠。 首先去除tips,由长度和覆盖率来定义。然后去除bubbles,运用的是和Dijkstra’s算法相似的算法,最后低覆盖率的点会被淘汰。

三、使用

合并两个fastq文件,如果是fasta,则需将下面的代码中fastq改为fasta

/sam/velvet/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl output_forward_paired.fq output_reverse_paired.fq paired.fq

第一步:建立哈希索引

直接对fastq格式的原始文件进行处理,首先是用velveth命令建立hash表子集 输入./velveth会出来使用帮助:

./velveth directory hash_length {[-file_format][-read_type][-separate|-interleaved] filename1 [filename2 ...]} {...} [options]

  directory    : directory name for output files
  hash_length        : EITHER an odd integer (奇数)(if even, it will be decremented) <= 127 (if above, will be reduced)
                         : OR: m,M,s where m and M are odd(奇数) integers (if not, they will be decremented) with m < M <= 127 (if above, will be reduced)
                              and s is a step (even number). Velvet will then hash from k=m to k=M with a step of s
  filename    : path to sequence file or - for standard input
File format options:
  -fasta           -fastq     -raw  -fasta.gz    -fastq.gz    -raw.gz       -sam -bam -fmtAuto
  (Note: -fmtAuto will detect fasta or fastq, and will try the following programs for decompression : gunzip, pbunzip2, bunzip2(默认的是fasta格式)

File layout options for paired reads (only for fasta and fastq formats):
  -interleaved        : File contains paired reads interleaved in the one file (default)
  -separate   : Read 2 separate files for paired reads

Read type options:
  -short        
-shortPaired 
  -short2       独立的插入大小的文库
-shortPaired2  同上
  -long -longPaired   针对sanger和454
  -reference

Options:
  -strand_specific  : for strand specific transcriptome sequencing data (default: off)
                   原来这个可以针对转录组数据啊
  -reuse_Sequences        : reuse Sequences file (or link) already in directory (no need to provide original filenames in this case (default: off)
  -reuse_binary      : reuse binary sequences file (or link) already in directory (no need to provide original filenames in this case (default: off)
  -noHash                        : simply prepare Sequences file, do not hash reads or prepare Roadmaps file (default: off)
  -create_binary            : create binary CnyUnifiedSeq file (default: off)

Synopsis:(大纲)

- Short single end reads:
velveth Assem 29 -short -fastq s_1_sequence.txt

- Paired-end short reads (PE1和PE2是合并的remember to interleave paired reads):
velveth Assem 31 -shortPaired -fasta interleaved.fna
#  Assem为生成文件所在文件夹,31为设置的kmer,一般为奇数,默认的是31,最开始默认最大的也为31,我已经改为127了,-shortPaired双末端模式,interleaved.fna?难道是是它指定的排列顺序?默认的输入文件是fasta格式

- Paired-end short reads (PE1和PE2是分开的)
velveth Assem 31 -shortPaired -fasta -separate left.fa right.fa

- Two channels and some long reads:
velveth Assem 43 -short -fastq unmapped.fna -longPaired -fasta SangerReads.fasta

- Three channels:
velveth Assem 35 -shortPaired -fasta pe_lib1.fasta -shortPaired2 pe_lib2.fasta -short3 se_lib1.fa

Output:
directory/Roadmaps
directory/Sequences
          [Both files are picked up by graph, so please leave them there]

第二步:根据建立的索引来组装

在终端输入./velvetg 查看说明

Usage:
./velvetg directory [options]

  directory                       : working directory name

Standard options:
  -cov_cutoff : removal of low coverage nodes AFTER tour bus or allow the system to infer it
            (default: no removal)
  -ins_length             : expected distance between two paired end reads (default: no read pairing)
  -read_trkg              : tracking of short read positions in assembly (default: no tracking)
                          追踪段read的位置(默认的是关)
  -min_contig_lgth   : minimum contig length exported to contigs.fa file (default: hash length * 2)
  -amos_file              : export assembly to AMOS file (default: no export)
  -exp_cov      : expected coverage of unique regions or allow the system to infer it
            (default: no long or paired-end read resolution)
  -long_cov_cutoff : removal of nodes with low long-read coverage AFTER tour bus
            (default: no removal)

Advanced options:
  -ins_length*                    : expected distance between two paired-end reads in the respective short-read dataset (default: no read pairing)
  -ins_length_long    : expected distance between two long paired-end reads (default: no read pairing)
  -ins_length*_sd     : est. standard deviation of respective dataset (default: 10% of corresponding length)
            [replace '*' by nothing, '2' or '_long' as necessary]
  -scaffolding            : scaffolding of contigs used paired end information (default: on)
  -max_branch_length       : maximum length in base pair of bubble (default: 100)
  -max_divergence : maximum divergence rate between two branches in a bubble (default: 0.2)
  -max_gap_count    : maximum number of gaps allowed in the alignment of the two branches of a bubble (default: 3)
  -min_pair_count    : minimum number of paired end connections to justify the scaffolding of two long contigs (default: 5)
  -max_coverage     : removal of high coverage nodes AFTER tour bus (default: no removal)
  -coverage_mask : minimum coverage required for confident regions of contigs (default: 1)
  -long_mult_cutoff                  : minimum number of long reads required to merge contigs (default: 2)
  -unused_reads                 : export unused reads in UnusedReads.fa file (default: no)
  -alignments            : export a summary of contig alignment to the reference sequences (default: no)
  -exportFiltered       : export the long nodes which were eliminated by the coverage filters (default: no)
  -clean                     : remove all the intermediary files which are useless for recalculation (default : no)
  -very_clean             : remove all the intermediary files (no recalculation possible) (default: no)
  -paired_exp_fraction       : remove all the paired end connections which less than the specified fraction of the expected count (default: 0.1)
  -shortMatePaired*          : for mate-pair libraries, indicate that the library might be contaminated with paired-end reads (default no)
  -conserveLong                 : preserve sequences with long reads in them (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

velvetg Assem -cov_cutoff 30 -ins_length 350 -ins_length_sd 100 -exp_cov auto -min_contig_lgth 500
#-ins_length为Insert size加上reads的长度,我前者为160,后者为91+91,所以后面的参数为350 ;默认的
#-ins_length_sd 100为什么设置这两个参数就有点不明白了,Assem注意和第一步建立的文件夹的对应

四、讨论

1.velvet拼接过程中主要的几个需要调试的参数

主要有三个:Kmer值,exp_cov和cov_cutoff三个。本文分别设置不同的参数值,进行比较。 Kmer值是基于De Bruijn算法的高通量读段拼接软件中最重要的一个参数,Kmer必须为奇数,在执行velveth时设置。在velvet的使用说明中,有一段关于Kmer设置的一般要求,如下: Ck= C(L-K+1)/L L表示读段长度,K表示Kmer值,C表示碱基的覆盖深度,Ck表示kmer的覆盖深度。L值已知,C值通过基因组规模和测序量来估计,当K值设定之后,可算得Ck值。根据经验,Ck值应当大于10,才能较好地完成拼接,而当Ck值大于20时,就浪费了测序深度。Kmer值越大Ck值就越小。当测序深度加大的时候,Kmer值也可以设的更大,有利于拼接。这个公式只能确定Kmer的取值范围,在此范围内选择采用哪个Kmer值,则需要进行Kmer参数值的调试。

   在执行velvetg时,可设定exp_cov和cov_cutoff两个选项。选项cov_cutoff用来过滤低覆盖深度的contigs,选项exp_cov为预期的覆盖深度。这两个参数可设定具体的值,默认值为“auto”。当设为auto时,软件在拼接时自动取值,如果数据的覆盖比较均匀,设为auto,拼接结果会比较好;exp_cov设为auto时,cov_cutoff也将强制为auto,auto表示覆盖深度的阈值为exp_cov的一半。

2.如何使用多个Kmer值进行拼接

If you wish for Velvet to e?ciently test various k-mer lengths without doing redundant computations, you can request all k-mer lengths such that m ≤ k < M with a step of s:

> ./velveth output_directory/ m,M,s (..data files..)
几个小贴士:
1,  velvet容许通道|,*代替任意字符等快捷方式
2,  可以将第一步分成两步,据说可以节约时间,不理解
> ./velveth output_directory/ 21 (..data files..) -noHash
> ./velveth output_directory/ 21 -reuse_Sequences
   3,对于数据量很大,可以压缩包来读取
> ./velveth output_directory/ 21 -create_binary (..data files..)

##Ps:

修改MAXKMERLENGTH修改了半天,我在查看其MAXKMERLENGTH值的时候是cd到我的安装文件,然后velvet,这样能看到他的值,实际上这个值是环境变量的值,而不是我这个安装程序的值,我不知道我系统之前安装过嘛,最后只有sudo apt-get remove –auto-remove velvet,删除这个velet后,在重新安装,重新修改环境变量。通过这个折腾,我深刻了理解了./velveth和velveth这两个调用程序还是有本质的区别的. 最后感谢 QQ好友:生物信息学-小徐 不厌其烦的帮助我找问题,谢谢罗

参考资料:

阿飞与安达博客:http://liuwei441005.blog.163.com/blog/static/135705811201211751914615/

xiaoshijun的个人博客 http://blog.sciencenet.cn/blog-303373-715739.html

微斯人 http://mingkang1217.blog.163.com/blog/static/20352272011620250132/

百度文库 http://wenku.baidu.com/view/2884921dfc4ffe473368abf4.html

软件自带说明书 强烈推荐这个说明书,逻辑和易读性很强

个人公众号,比较懒,很少更新,可以在上面提问题:

更多精彩,请移步公众号阅读:

Sam avatar
About Sam
专注生物信息 专注转化医学