Velvet的安装和使用

1. Velvet的安装

Velvet用于基因组的de novo组装,支持各种原始数据,包括Illumina的short reads和454的long reads。

首先下载velvet的安装包,直接使用make命令来编译,即可获得可执行主程序velveth和velvetg。安装如下:

$ wget http://www.ebi.ac.uk/~zerbino/velvet/velvet_1.2.10.tgz
$ tar zxf velvet_1.2.10.tgz
$ cd velvet_1.2.10
$ make 'CATEGORIES=10' 'MAXKMERLENGTH=57' 'LONGSEQUENCES=1' \
  'OPENMP=1' 'BUNDLEDZLIB=1'

值得注意的是,make后有多种参数,需要对velvet软件进行需要的设置:

CATEGORIES=10: 默认情况下velvet只能输入 2 groups of short reads,此处设置为10. 如果有多个文库或多种类型的原始数据,需要相应增加该值的大小。该值越大,耗内存越大。

MAXKMERLENGTH=57: 最大的Kmer长度,默认情况下该值为 31, 一般de novo组装基因组,31 是不够的,故需要增大该值。组装过程中kmer越大,越耗内存。

BIGASSEMBLY=1: 当超过 2.2G 的reads用于组装基因组的时候,需要设置该值。实际上很少会有如此之多的reads用于基因组组装,可以不用设置该值。设置该值后,会消耗更多的内存。

LONGSEQUENCES=1: 当组装出的contigs长度超过 32kb 长的时候,需要设置该值。会消耗跟多内存。

OPENMP=1:打开多线程运行。需要设置环境变量 OMP_NUM_THREADS 和 OMP_THREAD_LIMIT。 Velvet最多使用 OMP_NUM_THREADS+1 或 OMP_THREAD_LIMIT 个线程。也值哟部分的velvet算法支持多线程运行,不会造成运行时间的线性减少。

BUNDLEDZLIB=1: 默认velvet使用系统自带的zlib,如果系统没有zlib,则需要加入该参数来使用velvet源码包中自带的zlib。

在上述 make 编译完 Velvet 后,继续velvet的使用

2. 使用velveth来准备数据

velveth接受输入的文件,产生一个hash表。生成两个文件:Sequences和Roadmaps。velveth用法为:

$ velveth
$ velveth output_directory hash_length \
  -file_format -read_type filename1 filename2 ... \
  [-file_format -read_type filename1 filename2 ...]

在不带参数情况下,直接运行velveth会给出其帮助文件。而其参数如下:

output_directory:输出文件夹的路径

hash_length: 设置Kmer的大小。该值3点要求:1.必须为奇数;2.必须小于或等于编译velvet时设置的MAXKMERLENGTH值;3.必须小于reads的长度。

file_format: 支持的格式有fasta(默认)、fastq、bam等。

reads_type: short(默认)、shortpaired、short2、shortpaired2 … short10、shortpaired10、long、longPaired。 默认情况下short reads只保留了2个通道。在之前设置中将CATEGORIES值设置为10,则velvet最多支持10种不同类型的short reads数据。long支持长的reads,比如sanger测序数据和454测序数据。 如果reads_type一致,则同一个reads_type下可以有多个文件。

filename1: reads的文件名。如果是成对的reads,则需要将两个reads文件合并成一个文件。Velvet安装文件夹中提供了这样的一个文件。

一个velveth的例子。对5个Illumina测序的结果和一个454测序的结果使用velvet进行组装:

$ velveth output/ 31 \
  -shortPaired -fastq fragment1.reads.fastq \
  -shortPaired2 -fastq fragment2.reads.fastq \
  -shortPaired3 -fastq jumping1.reads.fastq \
  -shortPaired4 -fastq jumping2.reads.fastq \
  -long -fastq 454.fasta

3. 使用velvetg来进行基因组组装

velvetg是vlevet软件的进行de Bruijin图构建和操作的核心。在命令行下直接键入命令velvetg,这样描述:velvetg – de Bruijn graph construction, error removal and repeat resolution。其使用方法如下:

$ velvetg
$ velvetg directory [options]

$ velvetg output -exp_cov auto -cov_cutoff auto \
  -shortMatePaired3 yes -shortMatePaired4 yes \
  -clean yes -scaffolding yes -amos_file yes

velvetg的具体参数如下:

directory 工作的目录名,即为上一步骤velveth中的输出文件夹。

-cov_cutoff <float|auto> default: no removal
设置最低kmer覆盖度的值。默认下会生成很多nodes,而有些nodes很短,覆盖度较低,需要去除掉。auto则自动设置该值,是该值为所有nodes的kmer覆盖度值的median值的1/2。

-exp_cov <float|auto> default: no long or paired-end read resolution。
期望的kmer覆盖度。如果设置了auto,则该值为所有nodes的kmer覆盖度值的median值; 该值设置为auto,则同时自动设置-cov_cutoff为auto。如果对杂合基因组进行组装时,设置auto,却很难进行预测,组装结果肯定不好。auto适用于标准的基因组测序。

-long_cov_cutoff <float> default: no removal
移除低long-read覆盖度的nodes。

-ins_length <interger>
成对的short reads间的distance的期望值,即插入片段的平均长度。

-ins_length* <interger>
*代表第*组shortPaired reads, 和 velveth 中的参数相对应。

-ins_length_long <interger>
和前两个参数一样,代表longPaired reads的插入片段长度。

-ins_length*_sd <interger>
此时’*’代表’nothing、2、3…、long’等,和上面的3个参数相匹配;该值设置插入片段长度的标准差。有关插入片段长度和sd的如果不设置,velvet则会自动计算。velvet是将成对的reads比对到组装出来的nodes上,最后计算出一个insert size 和sd。这样做的话,可能估算的不准确,需要看velvet的运行输出信息,检测其预算是否准确。

-scaffolding <yes|no> defautl: yes
是否要使用paired end信息进行scaffolds组装

-max_branch_length <integer> default: 100
处理气泡(bubble)的参数。默认下,如果两条序列超过了100bp的位点处的kmer不一致,则将这两条序列分开成单独的contigs。

-max_divergence <float> default: 0.2
在一个bubble内两条branches最大的差异率(分母为较长的branch的长度)。

-max_gap_count <integer> default: 3
在一个bubble内两条branches比对结果中,运行的最大gap数。该参数和上两个参数为重要的参数,能很大程度上影响组装结果。如果设置得松点,可分别将值设为200,0.33,5。

-min_pair_count <integer> default: 5
默认,将两个contigs连成scaffold最少需要5对paired reads的验证。如果测序深度较大,则可以提高该值;测序深度低,则降低该值。

-long_mult_cutoff <integer> default: 2
默认下,融合两个contigs需要最少有2个long reads验证。

-max_coverage <float> default: no removal
最高的覆盖度,高于此覆盖度的nodes将被删除。

-coverage_mask <integer> default: 1
contigs最小的置信覆盖度,低于此覆盖度的contigs被maksed。

-shortMatePaired* <yes|no> default: no
如果哪一个shortPaired为mate-pair library测序的结果,则需要指定该参数为yes。

-conserveLong <yes|no> default: no
是否保留含有long reads的序列

-unused_reads <yes|no> default: no
是否输出unused reads, 保存到 UnusedReads.fa 中。

-alignments <yes|no> default: no
是否输出contig比对到参考序列的summary.

-exportFiltered <yes|no> default: no
是否输出由于覆盖度的原因被过滤掉的long nodes。

-clean <yes|no> default: no
是否删除所有的不能用于重新计算的中间文件

-very_clean <yes|no> default: no
是否删除所有的中间文件(删除后不能重新计算)

-min_contig_lgth <integer> default: hash length * 2
输出结果中最小的contigs长度

-amos_file <yes|no> default: no
是否输出AMOS文件

velvetg的默认输出结果

directory/contigs.fa 长度2倍长于kmer的contigs。参数-scaffolding决定生成的该fasta文件是否包含scaffold序列。
directory/stats.txt 用于决定覆盖度cutoff的统计表
directory/PreGraph 初始的de vruijin图
directory/Graph2 最终的de bruijin图 关于该文件中内容的解释,请见velvet PDF manual。
directory/velvet_asm.afg AMOS兼容的组装文件,能用于AMOS基因组组装软件包
directory/Log velvet的运行记录

4. Velvet提供了额外的一些scripts

这些scripts非常有用,位于$velvetHome/contrib/文件夹下。

4.1 estimate-exp_cov

在使用velvetg组装出一个初步结果后,利用结果文件stats.txt来估算出kmer的期望覆盖度。

$ $velvetHome/contrib/estimate-exp_cov/velvet-estimate-exp_cov.pl \
  [options] stats.txt

4.2 observed-insert-length.pl

在使用velvetg组装出一个初步结果后,以结果文件夹为输入,计算insert size和insert sd。

$ $velvetHome/contrib/observed-insert-length.pl/observed-insert-length.pl \
  [options] Velvet_directory

4.3 shuffleSequences

将对称的reads1和reads2文件合并成一个文件,用于velvet的输入。velvet不能将两个文件用于输入。

$ $velvetHome/contrib/shuffleSequences_fasta/shuffleSequences_fastq.pl \
  reads1.fastq reads2.fastq reads.fastq

4.4 VelvetOptimiser

输入该脚本的命令,不加参数,给出详细的帮助信息。该程序用于选取不同的Kmer值,来对原始数据进行组装,同时计算出最优化的组装参数。最后给出最优化方案的组装结果。非常适合于简单基因组的自动化组装操作。其参数设置较为简单。当然需要能在$PATH中有velveth和velvetg两个主要的velvet命令。以下就不详细介绍该命令参数,仅给出一个例子:

$ $velvetHome/contrib/VelvetOptimiser-2.2.4/VelvetOptimiser.pl \
  --s 17 --e 97 --x 2 --a --t 4 --k 'n50' --c 'Lbp' \
  --f '-shortPaired -fastq fragment1.reads.fastq \
  -shortPaired2 -fastq fragment2.reads.fastq \
  -shortPaired3 -fastq jumping1.reads.fastq \
  -shortPaired4 -fastq jumping2.reads.fastq \
  -long -fastq 454.fasta' \
  -o 'shortMatePaired3 yes shortMatePaired4 yes' \
  [-g 40M]

以上命令意义为:使用kmer长度从17到97,每次运行kmer长度+2; –a表示要生成amosfile;使用N50来决定kmer的选取;以’Lbp’,即大于1kb的contigs的总碱基数来决定coverage cutoff的值的选取;–f后为velveth的参数;–o后接velvetg的参数;-g后为预计的基因组大小,用于计算该命令下内存的消耗,然后退出运行,由于该程序会同时运行多个velvet和velveg,消耗的内存巨大。

4.5 AssemblyAssembler

该程序能使用不同的kmer组装出很多Assemblies,然后将这些Assemblies组装融合出一个最终的结果。其使用方法位于该程序的通目录下的README文本文件中。给出一个例子如下:

$ $velvetHome/contrib/AssemblyAssembler1.3/AssemblyAssembler1.3.py \
  -s 17 -e 99 -v $velvetHome -c 5.1 -t 35 -i 300 -m 51 -r y \
  -f "-fastq -shortPaired reads.fastq"

以上命令意义为:使用从17到99的kmer,每种kmer都组装出基因组并最后合并成一个结果。-c为coverage cutoff值;-t设置为只有大于35的kmer时,才进行coverage cutoff;-i表示插入片段长度为300;-m表示期望的kmer覆盖度为51;-r表示结果给出short reads要包含到最后的组装结果中。

可以看出该程序适合于只有一种类型的short reads时,才适合。因为其参数-i只能设置一种插入片段长度大小。

4.6 其它程序

extractContigReads用于提取和指定conting相对应的reads;

afg_handing 用于检查 .afg 文件;

fasta2agp 用于将fasta格式的基因组组装结果(gaps用N表示)转换成AGP文件,用于提交到EMBL或NCBI;

show_repeats 指出Assembly中larger repeated contigs的长度;

read_prepare 和 select_paired 用于NGS数据的过滤

5. 有参考基因组下velvet的组装

转录组De novo组装软件

1. Trinity

Trinity,本博客之前介绍了:Trinity的安装与使用

Citation: Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011 May 15;29(7):644-52. doi: 10.1038/nbt.1883. PubMed PMID: 21572440.

2. Oases

Oases是EMBL-EBI开发出来的,是一种使用short reads(such as Illumina, SOLiD, or 454)进行de novo转录组组装的软件。Oases是在velvet的基础上开发的,其运行需要安装Velvet。

Citation:M.H. Schulz, D.R. Zerbino, M. Vingron and Ewan Birney. Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics, 2012. DOI: 10.1093/bioinformatics/bts094.

上文献中写道:It was tested on human and mouse RNA-seq data and is shown to improve significantly on the transABySS and Trinity de novo transcriptome assemblers.

3. EBARDenovo

EBARDenovo, which stands for Extension, Bridging And Repeat-sensing Denovo. This software package (with patent pending) is free of charge for academic use only.

Citation:Hsueh-Ting Chu,William W. L. Hsiao,Jen-Chih Chen,Tze-Jung Yeh,Mong-Hsun Tsai,Han Lin,Yen-Wenn Liu,Sheng-An Lee,Chaur-Chin Chen,Theresa T. H. Tsao,and Cheng-Yan Kao。EBARDenovo: highly accurate de novo assembly of RNA-Seq with efficient chimera-detection Bioinformatics (2013) 29 (8): 1004-1010 first published online March 1, 2013 doi:10.1093/bioinformatics/btt092

上文献中写道:our algorithm is the most accurate among the examined programs, including de Bruijn graph assemblers, Trinity and Oases.

4. Trans-ABySS

Trans-ABySS: Analyze ABySS multi-k-assembled shotgun transcriptome data.

Trans-ABySS 1.4.4, Released Oct 09, 2012, Supports both transcriptome and genome assemblies. A more robust pipeline and improvement to the 3 analysis suites (fusions, indels, and splicing).
Citation:Robertson G,et al. De novo assembly and analysis of RNA-seq data. Nat. Methods 2010;7:909-912.

Variant calling 软件

1. SNP和INDEL caller

1.1 GATK

1.2 SAMtools

1.3 DISCOVAR

DISCOVAR的官方博客:http://www.broadinstitute.org/software/discovar/blog/。该软件和ALLPATHS-LG是同一个组织开发的,是在ALLPATHS-LG的基础上开发的,用于基因组的组装和Variants的发掘。2013.3 公布出了该软件。

1.4 PyroHMMvar

PyroHMMvar: a sensitive and accurate method to call short INDELs and SNPs for Ion Torrent and 454 data. 2013.8.28发表在bioinformatics上。