GENEWISE 的使用

1. GeneWise 简介

Genewise主要用于将蛋白质序列和DNA序列进行比对,从而对DNA序列上的编码区进行预测。

2. GeneWise 安装

$ wget http://www.ebi.ac.uk/~birney/wise2/wise2.4.1.tar.gz
$ tar zxf wise2.4.1.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/src

$ yum install *glib*
$ find ./ -name makefile | xargs sed -i 's/glib-config/pkg-config --libs glib-2.0/'
$ export C_INCLUDE_PATH=/usr/include/glib-2.0/:/usr/lib64/glib-2.0/include/:$C_INCLUDE_PATH
$ perl -p -i -e 's/getline/get_line/g' ./HMMer2/sqio.c
$ perl -p -i -e 's/isnumber/isdigit/' models/phasemodel.c
$ make all
$ export WISECONFIGDIR=/opt/biosoft/wise2.4.1/wisecfg/
$ make test
$ echo 'PATH=$PATH:/opt/biosoft/wise2.4.1/src/bin/' >> ~/.bashrc
$ echo 'export WISECONFIGDIR=/opt/biosoft/wise2.4.1/wisecfg/' >> ~/.bashrc 
$ source ~/.bashrc

3. GeneWise的使用

在GeneWise的安装目录下,有一个wise2.pdf文件,阐述了详细的genewise的使用方法。其软件最常用的命令是genewise。该命令的常用示例:

genewise protein.fasta dna.fasta -both -gff

程序输入的蛋白质序列和DNA序列分别是2个fasta文件。这两个fasta文件中仅有第一条序列是有效的,genewise仅对其中的2个第一条序列进行比对。以上示例对dna序列的正负链都进行cds预测,并将gff格式结果文件输出到标准输出。

genewise的常用参数:

-trev
    仅对负义链进行cds预测。
-tfor
    仅对正义链进行cds预测。该参数是默认值。
-both
    对负链都进行cds预测。
-genes
    给出gene结构的结果,非常简单的exon信息结果。默认情况下仅输出适合人类阅读的比对结果。
-gff
    给出gff格式的结果。
-cdna
    给出cdna序列。
-pep
    给出cds翻译出的蛋白质序列。
-splice [model/flat]
    使用的split site是model(默认值)或GT/AG。
-help
    给出帮助信息。
-version
    给出版本信息。
-silent
    标准错误输出不输出messages信息。
-quiet
    标准错误输出不输出report/info信息。

genewise的运行原理简述:

1. genewise的算法:21:93算法是genewise的基础算法。该算法简单讲就是 Match-Insert-Delete,在蛋白质序列和DNA序列比对后能准确划定intron边界。算法将intron分成5部分:5'端splice site、中间intron主体、富含CT区域、连接区、3'端splice site。根据蛋白质序列和DNA序列的比对结果算出Intron部分,从而将DNA序列的CDS区分成了Match、Insert和Delete 3部分,再对这3部分进行蛋白质翻译或移码翻译,从而划定intron边界,得到CDS信息。
2. 6:23算法则是2:93算法的简单版本,也是软件的默认设置。和2:93算法相比,6:23算法的intron没有第3和第4部分(富含CT区域、连接区)。同时,6:23算法更适合于DNA序列中没有屏蔽重复或introns序列比较怪异的情况。使用该算法的时候,-intron参数的值得是tied(也是该参数默认的值),否则会得到错误的很长的intron结果。
3. 若是算法后面带个 L 字样,则表示适用于进行输入的蛋白质序列是 HMM 模型。此外, 还有其它的一些算法,可以参考wise2.pdf文件。
4. genewise对基因进行预测后,有一个得分。该得分 = log2(预测模型的可能性/随机结果的可能性) 。因此,0表示该结果是个随机的结果,不可靠的。根据软件作者的经验,得分高于35的结果是非常可靠的;得分25~35的结果是可信的;得分18~25的结果可能仅适用于某些蛋白质家族;得分低于15的是不可信的。

4. GeneWise的高级使用

用临近物种的protein序列对基因组进行homolog gene预测的时候,需要通过blast将proteins序列和基因组序列进行比对,再提取基因组的目标基因区域和最佳结果protein进行genewise分析。因此,需要自己写一些程序进行并行化的genewise计算,从而达到对全基因组大数据的分析。Genewise软件提供了一支程序/opt/biosoft/wise2.4.1/src/perl/scripts/blastwise.pl程序能进行该项处理(我没有用过该程序,我自己写了想要的程序)。

使用 ICORN 进行基因组核苷酸的修正

1. ICORN 简介

ICORN (Iterative Correction Of Reference Nucleotides), 能通过将 reads 比对到基因组,从而修正 SNP 和 小于 3bp 的 INDEL 位点。现在出了新的版本 ICORN2 。
ICORN 官网: http://icorn.sourceforge.net/
ICORN 参考文献: Otto T D, Sanders M, Berriman M, et al. Iterative Correction of Reference Nucleotides (iCORN) using second generation sequencing technology[J]. Bioinformatics, 2010, 26(14): 1704-1707.

2. ICORN 原理

ICORN 的原理如下图所示。ICORN2 调用 SMALT 将 reads map 到基因组序列上;然后 Call SNPs 和 <=3bp INDELs;根据其位点的覆盖度情况决定基因组上该位点的核苷酸类型。最后通过多轮迭代直到不能 call 到新的 variants 为止。 http://icorn.sourceforge.net/workflow.jpg

3. ICORN 的下载和安装

ICORN 的使用会调用第三方的软件: GATK、SMALT、samtools 和 SNP-o-AMTIC 。

$ wget ftp://ftp.sanger.ac.uk/pub4/resources/software/pagit/ICORN2/icorn2.V0.95.tgz
$ tar zxf icorn2.V0.95.tgz -C /opt/biosoft/
$ /opt/biosoft/ICORN2/icorn2.sh --help

4. ICORN 的使用

4.1 使用前准备

准备的输入文件是: readroot_1.fastq、readroot_2.fastq、genome.fasta,将这3个文件放置于工作目录中。
同时,需要给一些环境变量赋值:

设定程序所在的目录
$ export ICORN2_HOME=/opt/biosoft/ICORN2/
设定运行的线程数
$ export ICORN2_THREADS=24
设定输出信息的多或少,对 debug 有用
$ export ICORN2_VERBOSE=2

4.2 运行软件 $ /opt/biosoft/ICORN2/icorn2.sh readroot 350 genome.fasta 1 3

以上命令第 2 个参数是数据的 insert size; 1 和 3 代表迭代的起始和终止,表示迭代 3 次。如果需要继续迭代,则设置起始为 4 。上述程序的结果文件为 ICORN2.Query.contigs.fa.4 。在 ICORN2_3 文件夹中也含有该 fasta 文件,同时程序也生成了一个 gff 文件。作者推荐将这两个文件载入到 artemis 基因组浏览器中进行查看。

5. 思考

有些人在有参考基因组的状况下做了重测序,总是想要得到其重测序的基因组结果。可以使用 ICORN 对参考基因组进行修正,即得到了对应其品种的基因组序列。

使用 Edena 进行基因组组装

1. 简介

Edena 是一个基于 overlaps–graph-based 的 de novo assembler。仅能很好地使用 Illumina 数据进行基因组组装。软件能使用 direct-reverse (paired-ends) 和 reverse-direct (mate-pairs) 的 Illumina 数据。同时要求 reads 的长度必须一致,否则所有 reads 会从 3′ 端开始截短到最短的 reads 的长度。

2. Edena 下载和安装

$ wget http://www.genomic.ch/edena/EdenaV3.131028.tar.gz
$ tar zxf EdenaV3.131028.tar.gz -C /opt/biosoft
$ cd /opt/biosoft/EdenaV3.131028/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/EdenaV3.131028/bin/' >> ~/.bashrc
$ source ~/.bashrc

3. Edena 的使用

软件直接使用 Edena 命令即可。包含两个步骤: overlapping 和 assembling。

3.1 overlapping

使用例子:

$ edena -nThreads 24 -DRpairs frag1_1.fq frag1_2.fq frag2_1.fq frag2_2.fq\
 -RDpairs jump1_1.fq jump1_2.fq jump2_1.fq jump2_2.fq

最后的结果文件为 out.ovl

3.2 Assembling

使用例子:

$ edena -e out.ovl

最终的结果文件为 out_contigs.fasta 。

使用 CISA 进行多个 genome Assemblies 的整合

1. CISA 简介

CISA (Contig Integrator for Sequence Assembly). 主要用于细菌基因组(小基因组)的整合。该软件使用简单,算法较好。
CISA 软件官网:http://sb.nhri.org.tw/CISA/en/CISA
CISA参考文献:Lin S H, Liao Y C. CISA: contig integrator for sequence assembly of bacterial genomes[J]. PloS one, 2013, 8(3): e60843.

2. CISA 算法

CISA分4步进行,如下图所示:
CISA算法图

2.1 鉴定代表性的 contigs 并对其进行延伸

首先,将所有 Assemblies 的 contigs 进行合并,并按长度进行排序。然后调用 NUCmer 软件按从长到短的顺序对序列逐个进行基因组比对。通过此方法来得到序列在各个 Assemblies 之间的重叠状况,并将结果写入到 CISA1/explained.txt(序列重叠信息) 和 CISA1/Extend_info(序列延伸信息) 文件中。
设定的阈值为:非代表性序列与代表性序列的比对结果达到 >95% 覆盖度 和 >95% 的一致性; 或非代表性序列超出了代表性序列的两端时,比对结果 >80% 覆盖度 和 >95% 的一致性。
通过这种方法,最终得到延长的代表性序列,输出到 R1_Contigs.fa 文件中。
如图所示:对 3 个 Assemblies 进行 CISA 分析。首先找到所有 contigs 序列中最长的 contig 序列(蓝色线条表示)进行全基因组比对,得到了第 1 轮结果;然后使用剩下的序列(+号后面的序列)挑出最长的 contig 序列(蓝色线条表示),对剩下的序列使用 NUCmer 进行比对,得到第 2 轮结果; 继续使用剩下的序列挑出最长的 contig 序列(黑色线条表示),最为代表性序列进行比对;最终得到了 4 条(蓝蓝黑红)代表性序列,其中的空心椭圆和实心箭头分别代表延伸的位置和方向。

2.2 对 misassembled 的 contigs 进行去除或截断, 并对 contigs两端不确定区域进行截短

将 R1_Contigs.fa 的序列使用 NUCmer 比对到所有的 contigs 上。通过比对结果,从而对 misassembled 的 contigs 进行去除或截断, 或对 contigs两端不确定区域进行截短。
如图所示:
Phase(2)左图表示 1 条代表性序列比对到了 2 个不同的 contigs 的中间,则表明该代表性序列是错误的组装,则去除该代表性序列,同时引进相应的能比对到该代表性序列的 contigs。这些信息在 CISA2/Remove_Info 文件中。
Phase(2)右图表示 1 条代表性序列中间区域没有对应的 contigs 能匹配,则表明该序列是错误的组装,则将该区域进行去除,截断该代表性序列。有关Gaps 的信息在 CISA2/Gaps 文件中,若 Gap 长度大于其代表性序列长度的 95%,则去除该序列。
此外,对代表性序列首尾不确定性区域的信息在 CISA2/clip_info 和 CISA2/clip_out 文件中。
最终,对 misassembled 的 contigs 进行处理后,得到结果文件 R2_Contigs.fa。该文件的序列总长度偏大。

2.3 通过 Overlap 信息将 contigs 融合

如果 2 个 contigs 有 end-to-end 重叠,并且重叠区对较小 contigs 长度的比例 >30%, 则将其进行融合。使用多轮 blastn 进行迭代来融合 contigs 序列。同时也进行重复序列鉴定。
融合信息写入到每轮文件夹下的 Extend_info 文件中,同时每轮生成 temp.fa 文件,作为下一轮 blastn 的数据库构建的输入文件。
重复序列信息文件写入到每轮文件夹下的 Repeat_Region.txt 文件中。

2.4 通过 small overlap 进行 contigs 融合

上一步骤将 contigs 进行融合需要 >30% 覆盖度的 overlap,同时也鉴定了重复序列。如图所示:Phase3中,黑色与红色线条重叠区比例大,被整合到一起,该序列和其他两条蓝色序列的overlap比较small,不能进行融合。但是鉴定了其重复序列(绿色块表示)后,通过本步骤,若重复区域长度小于其 small overlap,则能将 contigs 进行融合。
此步骤和上一步骤类似,使用 blastn 进行多轮运算。同时将结果写入到每轮文件夹的 Extend_info 和 temp.fa 文件中。
最终的输出文件为 CISA.fa (由 cisa.config 配置文件指定的名称)

3. CISA 的下载与安装

$ wget http://sb.nhri.org.tw/CISA/upload/en/2014/3/CISA_20140304-05194132.tar
$ tar xf CISA_20140304-05194132.tar -C /opt/biosoft

4. CISA 的使用

4.1 创建 merge.config 配置文件

配置文件如下:

count=3 
data=abyss.fasta,title=abyss  
data=velvet.fasta,title=velvet
data=allpathslg.fasta,title=allpathslg
Master_file=merge_contigs.fa
min_length=100 (default:100)
Gap=11

4.2 Merge.py

$ /opt/biosoft/CISA1.3/Merge.py merge.config

根据上面 merge.config 配置文件信息,对 3 个 Assemblies 的 contigs 序列进行提取,并整合到 merge_contigs.fasta 文件中,同时输出各个 Assemblies 的 contigs size。

4.3 创建 cisa.config 配置文件

配置文件内容如下:

genome=47814562 (填写Assemblies中最大的 contigs size)
infile=merge_contigs.fa
outfile=cisa.contig.fa
nucmer=/opt/biosoft/MUMmer3.23/nucmer
R2_Gap=0.95 (default:0.95)
CISA=/opt/biosoft/CISA1.3
makeblastdb=/opt/biosoft/ncbi-blast/bin/makeblastdb
blastn=/opt/biosoft/ncbi-blast/bin/blastn

4.4 CISA.py

$ /opt/biosoft/CISA1.3/CISA.py cisa.config

启动主程序进行运算。基因组越大,Assemblies 数目越多,耗时越多(呈几何级数增加)。

5. 思考

CISA 应该要修改下,在第 2 步进行并行化运算,以提高速度。
可以分批进行两两融合,直到最终融合。多个两两融合进行并行化运算。
进行大基因组的融合时,CISA 不适用了。

多个 genome Assemblies 的整合的经验

进行多个 genome Assemblies 整合的软件有很多:GAA, GAM, GAM-NGS, GARM, CISA, minimus2, Reconciliator, MAIA 等。

我对其中的一些进行了尝试,结果下:

1. GAA:软件很快就给出了结果,但是结果明显是错的。我不知道是不是用错了,但是软件运行过程中没有报错。

2. GAM-NGS: 软件安装可能会出现一些问题。有时候能出结果, 将 master 和 slave 反过来运行却会有报错。同时使用 bowtie2 的结果作为输入会报错。该程序安装和运行太过容易出错。 没有使用 MUMer ,虽然说是能节约时间,其实生成 sorted 的 bam 文件也极为耗时的。 最后,对于我的数据,该软件的整合效果太差。

3. GARM:软件的 pipeline 明显没写好,直接报错,连使用其 example data 都狂报错。 一看下载的程序是数年前的,估计作者也没修改其程序了。

4. Reconciliator: 从其发表文章所述的链接,完全下载不到该软件。

5. MAIA: 貌似需要使用 matlab,没装 matlab 切需要用钱购买的软件,这个就麻烦了。

6. minimus2: 使用貌似需要一些其他格式的文件,貌似比较麻烦,同时基因组整合的效果听说比较差。我个人就不去费精力验证了。

7. GAM: 貌似 GAM-NGS 是在该软件基础上做的,那么这个估计也好不到哪儿去。

8. CISA: 该软件结果还不错,使用起来也非常简单。唯一的缺点是运算极度耗时,特别是运行第2步的时候,对于较大的基因组,使用此软件几乎是不可能的事情。推荐使用此软件。

使用 L_RNA_scaffolder 进行转录组序列对基因组序列的 rescaffolding

1. L_RNA_scaffolder简介

L_RNA_scaffolder能利用长的转录子序列来辅助基因组的scaffolding拼接。这些长的转录子序列可以是454/sanger/Ion_Torrent的测序结果,也可以是illumina测序的de novo组装结果。
在使用此软件进行scaffolding之前,需要将这些长的转录子序列使用BLAT软件map到基因组上,得到PSL(without heading)文件。将此PSL文件和基因组的初步组装结果文件输入到L_RNA_scaffolder,很快能得到结果文件。

2. L_RNA_scaffolder的安装

$ wget http://www.fishbrowser.org/software/L_RNA_scaffolder/downloads/L_RNA_scaffolder.tar.gz
$ tar zxf L_RNA_scaffolder.tar.gz
$ cd L_RNA_scaffolder
运行软件的主程序,不加参数给出帮助信息
$ ./L_RNA_scaffolder.sh

3. L_RNA_scaffolder的使用

先使用blat,将转录子序列比对到基因组
$ blat genome.fasta transcripts.fasta -noHead out.psl
再使用L_RNA_scaffolder.sh得到最终的基因组文件
$ mkdir L_RNA_scaffolder_output
$ /opt/biosoft/L_RNA_scaffolder/L_RNA_scaffolder.sh -d /opt/biosoft/L_RNA_scaffolder/ -i out.psl -j genome.fasta -o ./L_RNA_scaffolder_output
最终的结果文件为./L_RNA_scaffolder_output/L_RNA_scaffolder.fasta

使用 musket 进行 reads correction

1. musket 简介

Musket is a an efficient multistage k-mer spectrum based corrector for illumina short read data。
其发表文献为:Liu Y, Schröder J, Schmidt B. Musket: a multistage k-mer spectrum-based error corrector for Illumina sequence data[J]. Bioinformatics, 2013, 29(3): 308-315.
该软件使用简单,能进行多线程运行,速度快。

2. musket 下载与安装

$ wget http://cznic.dl.sourceforge.net/project/musket/musket-1.1.tar.bz
$ tar zxf musket-1.1.tar.bz -C /opt/biosoft
$ cd /opt/biosoft/musket-1.1/
$ make
$ echo 'PATH=$PATH:/opt/biosoft/musket-1.1/' >> ~/.bashrc
$ source ~/.bashrc

在软件编译前,可以修改 Makefile 中第二行内容 “MACROS = -DMAX_SEQ_LENGTH=200 -DMAX_KMER_SIZE=28”。此行表示,软件允许的最大 read 长度为 200, 最大的 kmer 为 28。可以根据自己的需求修改此值。

3. musket 参数

-k int unit
此参数后有 2 个值,前者为 k-mer 的大小,默认值为 21。后者是 distinct k-mers 的数目,如果是全基因组测序的数据,则可以理解为基因组的大小。当然,由于 read 中可能存在错误,因此 distinct k-mers 的数目比基因组要大。如果内存够用,则我个人推荐设置为基因组大小的 2 倍。该值的设置主要用于平衡内存的消耗,对结果无影响。
-o str
此参数设置出处文件的文件名。设置此参数后,则不管输入文件有多少,只会输出一个文件。此参数和下一个参数 "-omulti" 是互斥的,仅能设置其中一个。
-omulti str
此参数设置出处文件的文件名。设置此参数后,每个输入文件都会对应一个输出文件。
-p int
程序能使用的线程数。默认值为 2。该值要设置为 >= 2。
-zlib int
设置输出的文件是否压缩,默认值为 0, 即不压缩。
-maxtrim int
允许 trim 掉的最大碱基数。该参数默认值为 0, 即表示不进行 trim。若此参数不为 0, 则程序在对一个 read 进行 error correction 后,再对该 read 进行 trim。trim 的方法为:程序寻找此 read 中没有错误碱基的子序列,然后取出最长的子序列,若最长子序列长度满足此参数的设定,则进行 trim,若不满足参数设定,则不进行 trim。
-inorder
此参数是输出文件的 reads 顺序和输入文件一致。在进行 paired reads 的修正时,一定要加上该参数。
-lowercase
是输出文件中的被修正的碱基为小写字母。默认值为 0 。
-multik bool
设置是否可以使用多个 k-mer sizes 来进行 error correction。默认值为 0。当设置为 1 时, 则可以设置多个 k-mer 值,如 "-k 21 536870912 -k 23 268435456"。

4. 常用的一些例子

$ musket -k 21 536870912 -p 12 reads.fa reads2.fa -o output.fa
$ musket -p 12 reads.fa reads2.fa -o output.fa
$ musket reads.fa reads2.fa -o output.fa -inorder

$ musket lib1_reads1.fastq lib1_reads2.fastq lib2_reads1.fastq lib2_reads2.fastq -omulti output -inorder
上述命令生成 4 个文件 output.0 output.1 output.2 output.3 按顺序分别对应前 4 个输入文件。

使用 SSPACE 进行 scaffoding

SSPACE 能利用 paired reads 的比对结果,将 contigs 或 scaffolds 连接成 scaffolds。其参考文献:Boetzer M, Henkel C V, Jansen H J, et al. Scaffolding pre-assembled contigs using SSPACE[J]. Bioinformatics, 2011, 27(4): 578-579.

1. 安装 SSPACE

软件下载页面:http://www.baseclear.com/lab-products/bioinformatics-tools/sspace-standard/

$ tar zxf SSPACE-STANDARD-3.0_linux-x86_64.tar.gz
$ ./SSPACE-STANDARD-3.0_linux-x86_64/SSPACE_Standard_v3.0.pl

解压缩软件包后,运行软件文件夹中的 perl 程序即可运行 SSPACE。软件主目录下包含一些软件使用说明和示例等,其中 README 文件描述得非常详细。

2. SSPACE 使用方法

2.1 library 文件

首先要建立一个描述 library 信息的文本文件,例如:

Lib1 bwa file1.1.fasta file1.2.fasta 400 0.25 FR
Lib1 bowtie file2.1.fasta file2.2.fasta 400 0.25 FR
Lib2 bwasw file3.1.fastq file3.2.fastq 4000 0.5 RF
Lib2 TAB file4.tab 4000 0.5 RF
Lib3 TAB file5.tab 10000 0.5 RF
unpaired bowtie unpaired_reads1.fasta
unpaired bwasw unpaired_longreads1.gz

此 library 文件由多列组成,列与列之间由 1 个 空格 或 tab 分隔,各列意义如下:

第 1 列: library 名称。程序运行过程中产生的临时文件以此来命名; 多个行可以拥有同一个 library 名称,则其具有相同的 library 设置和不同的数据文件; 同时,libraries 必须按 insert size 来排序,inert size 最小的必须放到第一行,这是因为进行 scaffold 构建时,按此文件提供的 libraries 的顺序来输入数据的; unpaired reads, 则第一列是 ‘unpaired’。
第 2 列: 将 reads 比对到基因组上所使用的软件名, 可以为 bowtie 、 bwa 和 bwasw 等; 如果输入的数据是 reads 比对过后的 tab 格式结果,则此列为 “TAB”。
第 3,4 列: Fasta 或 Fastq 格式的双末端测序文件,并且文件中成对的 paired reads 必须在两个文件中并处于相同的行号上,同时,软件读取数据与序列的 headers 无关。如果是 unpaired reads,则仅需要第 3 列,为 tab 格式的 reads mapping 结果,过后详述。
第 5,6 列:第 5 列为 insert size 的期望值; 第 6 列为 insert size 允许的最小偏差。 比如,这两列值分别为 4000 和 0.5,则 insert size 在 2000-6000 之间的 pairs 才是有效 pairs。
第 7 列:paired-reads 的方向,有 FF,FR,RF 或 RR 几种选项。

2.2 程序参数

-l 输入的 library 文件
-s 输入的 Fasta 文件
-x 是否对 contigs 进行延长。其值可以为 0 或 1。 1 表示进行延伸,0 表示不延伸。默认值为 0。

延伸参数:
-m 进行延伸时,read 和基因组序列最小的 overlap。此值越大,则结果越准确,同时耗内存越少。推荐此值接近最长的 read 的长度。比如,对于 26 bp 长度的 reads, 该值适合设为 32~35。 默认此值为 32 。此值取值范围为 15~50 。软件运行时,将 unmapped reads 全部打断成 m+1 长度的序列,这些序列用于进行 contigs 的延伸。
-o 进行延伸时,延伸 1 个碱基需要的最小 reads 数。此值越大,则结果越准确。默认值为 20 。
-r 进行延伸时,延伸 1 个碱基,此碱基在所有匹配的 reads 中的最小比例。此值越大,则结果越准确。默认值为 0.9 。

Scaffolding 参数:
-k 将两个 contigs 连接成 scaffold 时,需要的最小的 reads pairs 数目。默认值为 5 。
-a 将两个 contigs 连接成 scaffold 时,这两个 contigs 之间的连接数 与 其和其它 contigs 的连接数之间的最小比值。此值越大,则结果越准确。默认值为 0.70
-n 在 scaffold 中,将两个邻近的 contigs 合并到一起需要的最小的 overlap。默认值为 15。
-z 进行 scaffolding 时,允许的最小的 contig 长度。低于此长度的 contig 将不能用于进行 scaffold 组装。默认值为 0 。较长的 contigs 产生的 scaffolds 比较可信; 而小于 100bp 的 contigs 容易是重复序列。

bowtie 比对参数:
-g 使用 bowtie 进行比对时,允许的最大 gaps 数。默认值为 0

其它参数:
-T 设定运行的线程数。默认值为 1。
-b 输出文件夹名及文件夹内的文件前缀。
-S 当程序正在运行时,跳过读取 reads 的阶段。和 -b 参数结合使用,则可以同时运行多个 SSPACE 程序,对每个程序设置不同的参数,这样能较快得到较好的结果。
-v verbose mode
-p 生成可供可视化的 .dot 文件。

2.3 其它工具

SSPACE 提供了一些其它比较有用的小工具:

estimate_insert_size.pl 用于计算 insert size。此程序计算的结果有些问题。
fastq_qualitytrim_pairs.pl 对 reads pairs 进行质量控制的程序。

sam_bam2tab.pl 将 bam sam 文件转换为 tab 格式的程序。

GAGE-B 进行微生物基因组组装软件的评估

参考文献:Magoc T, Pabinger S, Canzar S, et al. GAGE-B: an evaluation of genome assemblers for bacterial organisms[J]. Bioinformatics, 2013, 29(14): 1718-1725.

评价了 8 个基因组组装软件: ABySS CABOG MIRA MaSuRCA SGA SOAPdenovo SPAdes Velvet。
若这 8 个软件组装效果全部一致,则值全部为 1 。

对 7 个物种的组装结果, 根据 N50, contig 的组装性能:
ABySS 0.79
CABOG 0.55
MIRA 0.85
MaSuRCA 2.08
SGA 0.24
SOAPdenovo 1.10
SPAdes 1.91
Velvet 0.48

对 5 个物种的组装结果, 根据 N50, scaffold 的组装性能:
ABySS 0.76
CABOG 0.81
MIRA 0.85
MaSuRCA 1.53
SGA 0.31
SOAPdenovo 0.81
SPAdes 1.20
Velvet 0.72

可以看出,从基因组的连续性考虑,最佳的细菌基因组组装软件是:
MaSuRCA > SPAdes 。仅有这两个软件能高于平均水平, 特别是 MaSuRCA 显著比其它软件要好。

文中没有进行组装结果的准确性分析。但是提到 SOAPdenvo 可以产生最大的 scaffolds,但是比其它软件的结果有根多的 local errors。

由于数据问题,没有对 ALLPATHS-LG 进行分析。