基因组组装软件心得

1. SOAPdenovo 和 Velvet组装基因组是否需要对reads进行预处理?

使用soapdenovo组装基因组,需要先对reads进行预处理,去掉质量不好的reads,剪掉3’端和5’端的一些碱基;然后再进行组装。而使用velvet时候,则不需要,可能是由于velvet的cutoff比较好。
这个结论得出方法:将reads进行预处理前后分别用于组装出基因组,然后再使用bowtie2将两种reads比对到两种基因组组装结果,观察reads比对上的百分比得出。reads的预处理方法是去除低质量reads并将100bp的reads前后各减去5bp碱基。
对于soapdenovo:

                      未处理reads比对率         预处理reads比对率
未处理reads基因组          65.84%                   66.10%
预处理reads基因组          66.68%                   66.99%

对于velvet:

                      未处理reads比对率         预处理reads比对率
未处理reads基因组          86.41%                   88.10%
预处理reads基因组          70.09%                   71.62%

通过以上的对比和数据可以看出:1.velvet组装出的基因组的精确性要比soapdenovo高;2.对reads进行预处理对soapdenovo组装的基因组的精确性影响较小,预处理reads后,组装出的基因组精确性稍有提高;3.使用velvet组装基因组时,不对reads进行预处理组装出的基因组精确性高很多。
讨论:1.数据量不够多,碱基覆盖度只有20~25X, 以上结论具有一定的片面性; 2.velvet在对kmer进行cutoff的算法很棒,因而组装出的基因组精确性高。

Velvet的安装和使用

一. Velvet简介

VelvetEBIDaniel Zerbino and Ewan Birney 开发的利用short reads组装基因组的软件。
Velvet作为一种de novo基因组组装软件,可利用的reads包括 Illumina reads 和 454 reads。
引用文献(2008):Velvet: algorithms for de novo short read assembly using de Bruijn graphs.

二. Velvet下载与安装

make ‘CATEGORIES=57’ ‘MAXKMERLENGTH=75’ ‘BIGASSEMBLY=1’ ‘LONGSEQUENCES=1’ ‘OPENMP=1’ ‘BUNDLEDZLIB=1’

Quake的安装和使用

一. Quake简介

Quake是由CBCB(Center for Bioinformatics and Computational Biology)开发的运用于修正序列错误的软件。Quake采用k-mer的错误修正方式,特别适合于Illumina测序的short reads数据,将reads中的错误碱基进行修正,同时,必须要满足碱基的基因组覆盖度要>15X。其文章2010年发表在Genome Biology上。

二. Quake的安装

1. 下载并安装Boost
2. 下载Quake并安装

$ wget http://www.cbcb.umd.edu/software/quake/downloads/quake-0.3.4.tar.gz
$ tar zxf quake-0.3.4.tar.gz
$ cd Quake/src
$ make

3. 安装JELLYISH,并将jellyfish链接到quake

$ ln -s ..../jellyfish ../bin/jellyfish
这一步需要安装jellyfish,然后将其软链接到Quake的bin目录下。或者修改Quake的bin
目录下的quake.py脚本,将jellyfish所在的目录进行修正。

4. 安装R及R的软件包VGAM

$ R
> install.packages("VGAM")
> q(save="no")

5. 使用Quake的bin目录下的quake.py来运行程序

三. quake的使用参数

1. 主要参数:

-r READSF
    Fastq文件名
-f READS_LISTF
    一个文件,其中包含fastq文件名,每行一个文件名;如果是双末端reads,则是一行两
个文件名
-k K
    使用的k-mer的长度。如果基因组大小为G,则k-mer长度选择为: k ~= log(200G)
/log(4)
-p PROC  default: 4
    使用的CPU线程数
-q QUALITY_SCALE
    使用的碱基质量格式,一般是64或33.如果不给出,则软件会自行猜测

2. 计算k-mer参数:

--no_jelly
    使用quake自带的一个简单的程序来进行k-mer计数,而不是使用jellyfish。该程序
是单线程运行,速度相对慢,能满足像微生物基因组这样的小基因组测序,但是不适合于大的基
因组。
--no_count  default: false
    不进行k-mer计数,而其计数结果已经存于在目标文件[readsfile].qcts和[reads
file].cts中了。
--int  default:false
    对kmers以整数的方式进行计数,而不使用碱基质量值。默认情况下是利用到了碱基质量,
计数的称为qmer。
--hash_size=HASH_SIZE
    Jellyfish的参数,用来设置Hash的大小。如果不设置,Quake则会使用k来估计出一
个值来。

3. 覆盖度模型参数:

--no_cut  default: false
    Quake使用k-mer计数来画直方图,通过调用R脚本cov_model_qmer.r来画直方图,并
决定出Coverage cutoff的阈值。默认情况下是最优化的模型,从而得出一个cutoff的Co
verage值,此值将输出到文件cutoff.txt中。
--ratio=RATIO
    确定Coverage cutoff值的时候,该值处对应的qmer错误的可能性是正确的可能性的
倍数(默认是200倍),即正确率不足0.5%。该值越小,则阈值越松。

4. reads的修正参数:

-l MIN_READ
    输出的长度>=此值的修正后的reads。很多reads由于修正或修剪后会较短
--headers
    输出的fastq文件使用的头和原始文件一致,不包含修正信息
--log
    输出一个log文件,里面包含所有的修正日志,包括"碱基质量 位置 新的碱基 旧的碱基"。

四. Quake的结果

使用Quake对双末端测序文件 A.fastq 和 B.fastq 进行reads修正,则产生如下结果文件:

A.cor.fastq  B.cor.fastq
    修正过后的reads文件               
A.cor_single.fastq  B.err_single.fastq
    A文件中修正过后的reads结果;而B文件中的序列则是error reads
A.err_single.fastq  A.cor_single.fastq
    B文件中修正过后的reads结果;而A文件中的序列则是error reads
A.err.fastq  B.err.fastq
    A和B文件中成对的reads都是error reads

五. Quake的correct程序

Quake的bin文件中有一支C++的程序,用于reads的修正。其实Quake.py在运行过程会调用该程序,但是当出现基因组覆盖度较低,而cutoff的值出现异常结果,比如小于1时,可以考虑使用该程序来重新对reads进行修正,取cutoff值为1.0。其使用方法为:

$ correct -f [fastq list file] -k [k-mer size] -c [cutoff] -m [counts file] -p 24

ALLPATHS-LG的使用

一、ALLPATH简介

ALLPATHS-LG是一个基因组组装软件,适合于组装short reads数据,由Computational Research and Development group at the Broad Institute开发。ALLPATHS-LG是现在行业内公认进行基因组De novo组装效果最好的软件。

二. 基础注意事项

1. 不能只使用一个library数据进行组装;
2. 必须有一个"overlapping"的片段文库的paired-reads数据。比如,reads长度~
100bp,插入片段库长度~180bp;
3. 必须有jumping library数据;
4. 基因组组装需要100x或以上基因组覆盖度的碱基,这个覆盖度是指raw reads数据(在
error correction和filtering之前)的覆盖度;
5. 可以使用PacBio数据;
6. 不能使用454数据和Torrent数据。主要是这两者测序太贵,如果什么时候价格降低,有
需求的话,会写出相应的代码来满足要求;
7. 官方提供了测试用数据;
8. 不支持在整个计算机集群上进行运算;
9. 需要消耗的内存峰值大约是1.7bytes每个碱基,即输入10G的碱基数据量,大约需要17
G内存;
10. 对于试探性的参数,比如K,原则上可以调整。但是我们不会自行调整,并也不推荐。AL
LPATHS-LG不像其它De novo一样,Kmer大小的参数K和read大小之间没有直接的联系,
ALLPATHS-LG会在运行过程中运用一系列的K值。

三. ALLPATHS-LG使用方法

1. 基础的使用方法和命令

使用RunAllPathsLG这个命令来运行。虽然有很多参数,但是在没有指导的情况下不要随意使用,使用默认设置即可。其使用方法为:

$ RunAllPathsLG arg1=value1 arg2=value2 ...

参数主要是设置程序辨别的一些目录,在程序的运行过程,会输入相应目录中的数据,将结果输入到指定的目录。一个简单的命令使用例子:

#!/bin/sh

# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000
# ALLPATHS-LG命令的写法与一般的linux参数写法不是很一样。采用 ‘参数=值’ 的方法,并使之成每行一个参数,使用'\'来连接各个参数,这样看起来直观易懂。初始接触的人可能会不适应。

RunAllPathsLG \
 PRE=$PWD\
 REFERENCE_NAME=species.genome\
 DATA_SUBDIR=data\
 RUN=run\
 SUBDIR=test\
 EVALUATION=STANDARD\
 TARGETS=standard\
 OVERWRITE=True\
 MAXPAR=8
 | tee -a assemble.out

2. 详细的参数说明

必须的参数
PRE (String)
    程序运行的根目录,所有的其它目录全在该目录下
REFERENCE_NAME (String)
    参考基因组目录名称,位于PRE目录下。如果有一个参考基因组,可将参考基因组放到该
目录中;若没有,则创建该文件夹用于基因组组装
DATA_SUBDIR (String)
    DATA子目录名称,位于REFERENCE_NAME目录下。程序从该目录中读取数据。
RUN (String)
    运行目录名称,位于DATA_SUBDIR下。程序将生成的中间文件和结果文件存储于该目录
。比如组装结果是一个名为ASSEMBLES的目录,位于该目录下。

部分可选参数:
SUBDIR (String) default: test
    子目录名,在REF/DATA/RUN/ASSEMBLIES目录下创建的存放基因组组装结果的目录
名。
K (int) default: 96
    核心Kmer大小,只有K=96能很好地运行。
EVALUATION (String: {NONE,BASIC,STANDARD,FULL,CHEAT})default:BASIC
    给定一个参考基因组,pipeline能在基因组组装的不同阶段对组装过程和结果进行评估。
    BASIC:基础评估,不需要参考基因组;
    STANDARD:使用参考基因组来运行评估模块;
    FULL:在某些组装模块下打开in-place评估,不会影响组装结果;
    CHEAT:稍微使用参考基因组指导组装,产生更详细的分析,能对组装结果产生小的(好方
向的)改变。
REFERENCE_FASTA (String) default: REF/genome.fasta
    评估中使用的参考基因组。
MAXPAR (int) default: 1 
    有些模块的运行是独立的,不相互依赖,能同时运行。该参数设定能同时运行的模块的最
大数目。由于pipeline中的绝大部分模块都能多线程运行,因此将该值设定大于1,效果不明
显。
THREADS (String) default: max 
    有些模块能多线程程运行,默认使用最大线程数运行。
OVERWRITE (Bool) default: False 
    是否覆盖存在的文件。可以设置该选项为True,在每次运行程序的时候设定RUN参数为
一个新的目录名,则比较好。
TARGETS (vec) default: standard 
    pipeline会生成一系列的文件,不同的文件的生成需要call不同的模块。如果某文件
已经存在了并且是最新的,则跳过相应的模块的运行。本参数指定生成哪些拟定的目标文件(p
seudo targets)。若目标文件没有相应的模块能生成,则会得到报错。
    none:没有拟定的目标文件,仅仅生成指定的目标文件;
    standard:生成组装文件和选定的评估文件;
    full_eval:生成组装文件和额外的评估文件。
TARGETS_REF (String) 
    在ref_dir目录中生成的目标文件。
    多个目标文件的书写方法为: TARGETS_REF="{target1,target2,target3}" 。
TARGETS_DATA (String) 
    在data目录中生成的目标文件。
TARGETS_RUN (String) 
    在run目录中生成的目标文件。
TARGETS_SUBDIR (String)
    在subdir中生成的目标文件。 
FORCE_TARGETS (Bool) default: False
    生成目标文件,即使文件已经存在并且看起来是很新的。

3. 输入文件与目录的准备

两个文库:插入片段长度为180bp和3000bp,illumina测序文件结果为fastq格式。以此为例来准备ALLPATHS-LG运行所需的文件和目录。

(1) 准备 in_groups.csv 和 in_libs.csv 文件。

这两个文件内容由逗号隔开,in_groups.csv文件内容如下:

group_name, library_name, file_name
firest, Illumina_180bp, seq/species_500bp_read?.fastq
second, Illumina_3000bp, seq/species_3000bp_read?.fastq

in_groups.csv文件的解释:

group_name:数据独特的代号,每一份数据有一个代号;
library_name:数据所属文库的名字,体现出该;
filename:数据文件所存放位置。可以为相对位置,文件名可以包含'*'和'?'(但是扩展名
中不能有该符号,因为要根据扩展名识别文件类型),从而代表paired数据。支持的文件类型有
'.bam','fasta','fa','fastq','fq','fastq.gz'和'fq.gz'。

in_libs.csv文件内容如下:

library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end
Illumina_180bp, species, species.genome, fragment, 1, 180, 10, , , inward, 0, 0
Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0

in_libs.csv文件的解释:

library_name:和in_groups.csv中的相匹配;
project_name:project的名字;
organism_name:测序物种的名字;
type:仅仅只是一个信息;
paired:0:Unpaired reads;1:paired reads;
frag_size:小片段文库插入片段长度的均值;
frag_stddev:小片段文库的插入片段长度估算的标准偏差;
insert_size:大片段文库插入片段长度的均值;
insert_stddev:大片段文库插入片段长度估算的标准偏差;
read_orientation:reads的方向,小片段文库为inward,大片段文库为outward;
genomic_start:reads从该位置开始,读入数据,如果不为0,之前的碱基都被剪掉;
genomic_end:reads从该位置开始,停止读入数据,如果不为0,之后的碱基都被剪掉。

(2) 使用PrepareAllPathsInputs.pl来对数据进行转换

ALLPATHS-LG接受的输入数据要求如下:

1. ALLPATHS-LG的输入数据支持小片段文库(fragment library)、大片段文库(jum
ping library)和超大片段文库(long jumping library)。并且前两种文库至少各有
一个才能进行基因组组装。超大片段文库是只插入片段>20kb的文库,其测序方向和小片段文
库一致,为inward。

2. ALLPATHS-LG的输入数据放置在//文件夹下,包含3种文件:碱基文件,质量文件和配
对信息文件
   frag_reads_orig.fastb
   frag_reads_orig.qualb 
   frag_reads_orig.pairs 

   jump_reads_orig.fastb 
   jump_reads_orig.qualb 
   jump_reads_orig.pairs

以下是可选的超大插入片段文库对应的数据文件(非必须):

  long_jump_reads_orig.fastb 
  long_jump_reads_orig.qualb 
  long_jump_reads_orig.pairs

使用PrepareAllPathsInputs.pl来将fastq等格式的测序结果转换成ALLPATHS-LG可接受的文件。以下是该程序的参数:

DATA_DIR
    将转换后的数据文件放到此文件夹下。
PICARD_TOOLS_DIR
    若输入数据为bam格式,则需要用到Picard软件,该参数Picard的路径
IN_GROUPS_CSV
    输入的in_groups.csv文件名
IN_LIBS_CSV
    输入的in_libs.csv文件名
INCLUDE_NON_PF_READS default: 1
    1:包含non-PF reads;0:仅仅只包含PF reads.
PHRED_64 default: 0
    0:碱基质量是ASCII的33到126,一般情况下Illumina数据的最低碱基质量是'B';
1:碱基质量的ASCII码是从64到126,一般情况下Illumina数据的最低碱基质量是'#'。
PLOIDY
    生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型,2表
示双倍体型。
HOSTS
    列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3.
host2,4.host3"表示使用本地机器的2个CPU线程,host2机器的3个CPU线程和host3机
器的4个CPU线程。

以下是不常用的参数,主要用来选择转换的数据量的大小。当测序数据量太多,而只想使用其
中一部分数据的时候,可以用到
FRAG_FRAC
    使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值,则不能同时设定
FRAG_COVERAGE。
JUMP_FRAC
    使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值,则不能同时设定
JUMP_COVERAGE。
LONG_JUMP_FRAC
    使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值,则不能同时
设定LONG_JUMP_COVERAGE。
GENOME_SIZE
    估计的基因组大小,用来计算对应覆盖度所对应的reads数
FRAG_COVERAGE
    所期望的小片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
JUMP_COVERAGE
    所期望的大片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
LONG_JUMP_COVERAGE
    所期望的超大片度库的覆盖度,比如 1. 要求GENOME_SIZE有设定

一个 PrepareAllPathsInputs.pl 的例子如下:

#!/bin/sh

# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000

mkdir -p species.genome/data

# NOTE: The option GENOME_SIZE is OPTIONAL. 
#       It is useful when combined with FRAG_COVERAGE and JUMP_COVERAGE 
#       to downsample data sets.
#       By itself it enables the computation of coverage in the data sets 
#       reported in the last table at the end of the preparation step. 

# NOTE: If your data is in BAM format you must specify the path to your 
#       picard tools bin directory with the option: 
#
#       PICARD_TOOLS_DIR=/your/picard/tools/bin

PrepareAllPathsInputs.pl\
 DATA_DIR=$PWD/species.genome/data\
 PLOIDY=1\
 IN_GROUPS_CSV=in_groups.csv\
 IN_LIBS_CSV=in_libs.csv\
 OVERWRITE=True\
 | tee prepare.out

(3) 使用Fasta2Fastb来转换得到参考基因组的输入数据

将参考基因组的fasta文件改名为genome.fasta并放到PRE/REFRENCE_NAME/文件夹下。使用Fasta2Fastb来将fasta文件转换为其二进制文件,在PRE/REFRENCE_NAME/目录下生成genome.fastb文件。

$ Fasta2Fastb IN=PRE/REFRENCE_NAME/genome.fasta

4. ALLPATHS Cache的使用

PrepareAllPathsInputs.pl脚本实际上是包含两个步骤:
1. 将测序的原始文本数据(fastq等文件)转换成二进制数据(fastb,qualb文件),并将各个group测序数据的二进制数据放置到缓存文件夹:<DATA>/read_cache。
2. 将部分或者全部的二进制数据结合并生成所需要的运用于基因组组装的数据文件。因此,使用ALLPATHS Cache的优点:a. 当加入新的测序数据或进行不同数据量的基因组组装评估时,需要再次进行基因组组装,这时可以仅仅只转换新的测序数据,节约了重新进行数据转换的时间。只需要将二进制数据根据需要来生成ALLPATHS-LG接受的文件; b. 当测序数据的测序质量有phred64和phred33兼有的时候,则需要使用Cache来分别转换数据。

(1). 使用CacheLibs.pl将文库信息读入到Cache

使用方法:

CacheLibs.pl\
  CACHE_DIR=<CACHE_DIR>\
  IN_LIBS_CSV=in_libs.csv\
  ACTION=Add

使用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
IN_LIBS_CSV
    输入的in_libs.csv文件名
ACTION  deault: List
    Add,List 或者 Remove,从Cache中添加,列出或移除文库信息。

(2). 使用CacheGroups.pl将文本数据文件转换成二进制数据文件到Cache

使用方法:

CacheGroups.pl\
  CACHE_DIR=<CACHE_DIR>\
  PICARD_TOOLS_DIR=/opt/picard/bin\
  IN_GROUPS_CVS=in_groups.csv\
  TMP_DIR=/tmp\
  HOSTS='2,3.host2,4.host3'\
  ACTION=Add

使用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
PICARD_TOOLS_DIR
    若输入数据为bam格式,则需要用到Picard软件,该参数为Picard的路径
IN_GROUPS_CSV
    输入的in_groups.csv文件名
PHRED_64  default: 0
    输入的fastq的碱基质量格式,True: 碱基质量格式为PHRED64; False: 碱基质量
格式为PHRED33.
OVERWRITE  default: 0
    是否覆盖已存在的数据文件
TMP_DIR
    临时文件夹的路径,如果数据量够大,则必须要够大
INCLUDE_NON_PF_READS default: 1
    1:包含non-PF reads;0:仅仅只包含PF reads.
HOSTS
    列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3.
host2,4.host3"表示使用本地机器的2个CPU线程,host2机器的3个CPU线程和host3机
器的4个CPU线程。
ACTION  default: List
    Add,List 或者 Remove,从Cache中添加,列出或移除group的数据信息。

(3). 使用CacheToAllPathsInputs.pl来生成ALLPATHS的输入数据。
使用方法:

CacheToAllPathsInputs.pl\
  CACHE_DIR=<CACHE_DIR>\
  DATA_DIR=DATA_DIR\
  GROUPS="{group1,group2}"\
  FRAG_FRAC=100%\
  JUMP_FRAC=100%

常用参数:

CACHE_DIR
    缓存文件夹的绝对路径名
DATA_DIR
    DATA文件夹就绝对路径名
GROUPS
    需要转换的到ALLPATHS的输入数据的groups名。格式为"{group1,group2,group
...}"
IN_GROUPS_CSV
    in_groups.csv文件名,和 GROUPS 参数二选一
FRAG_FRAC
    使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值,则不能同时设定
FRAG_COVERAGE。
JUMP_FRAC
    使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值,则不能同时设定
JUMP_COVERAGE。
LONG_JUMP_FRAC
    使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值,则不能同时
设定LONG_JUMP_COVERAGE。
FRACTIONS
    同时设置上述3个参数。比如"{0.5,30%,100%}"
GENOME_SIZE
    估计的基因组大小,用来计算对应覆盖度所对应的reads数
FRAG_COVERAGE
    所期望的小片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
JUMP_COVERAGE
    所期望的大片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定
LONG_JUMP_COVERAGE
    所期望的超大片度库的覆盖度,比如 1. 要求GENOME_SIZE有设定
COVERAGES
    同时设置上述3个参数,比如"{45,50,2}"
LONG_READ_MIN_LEN  default: 500
    设置被称为long unpaired read的阈值(适用于PacBio reads)
PLOIDY  
    生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型,2表
示双倍体型。如果没有该参数,则不生成ploidy文件

四. 思考题

对一个简单的真菌物种,运用Illumina Hiseq2000平台,构建了180bp,500bp,3000bp长度的3个DNA片段文库,分别进行了测序。获得了相应的fastq数据文件,如何使用ALLPATHS-LG来进行基因组的De novo组装?

1. 安装ALLPATHS-LG到Unix系统服务器上

2. 准备测序的数据文件

测序的数据文件如下,并放置到当前工作目录下的 seq 文件夹中。

180.reads1.fastq    180.reads2.fastq
500.reads1.fastq    500.reads2.fastq
3000.reads1.fastq   3000.reads2.fastq

3. in_libs.csv 和 in_groups.csv 文件的准备
两个文件都放置到当前工作目录下。 in_libs.csv 内容如下:

library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end
Illumina_180bp, species, species.genome, fragment, 1, 180, 20, , , inward, 0, 0
Illumina_500bp, species, species.genome, fragment, 1, 500, 50, , , inward, 0, 0
Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0

in_groups.csv 内容如下:

group_name, library_name, file_name
180, Illumina_180bp, ./seq/180.reads?.fastq
500, Illumina_500bp, ./seq/500.reads?.fastq
3000, Illumina_3000bp, ./seq/3000.reads?.fastq

4. 使用PrepareAllPathsInputs.pl来对数据进行转换

将以下内容写入prepare.sh并运行

#!/bin/sh
ulimit -s 100000
mkdir -p species.genome/data
PrepareAllPathsInputs.pl\
 DATA_DIR=$PWD/species.genome/data\
 PLOIDY=1\
 IN_GROUPS_CSV=in_groups.csv\
 IN_LIBS_CSV=in_libs.csv\
 OVERWRITE=True\
 | tee prepare.out

5. 运行ALLPATHS-LG主程序进行基因组的De novo组装

将以下内容写入到assemble.sh中,并运行

#!/bin/sh
ulimit -s 100000
RunAllPathsLG\
 PRE=$PWD\
 REFERENCE_NAME=species.genome\
 DATA_SUBDIR=data\
 RUN=run\
 SUBDIR=test\
 OVERWRITE=True\
 MAXPAR=8\
 | tee -a assemble.out

6. 结果文件

基因组组装结果文件位于./species.genome/data/run/ASSEMBLIES/test/final.assembly.fasta文件。同一个目录下有一个文件final.assembly.efasta也是组装结果。
efasta格式即“enhanced fasta”。
efasta和fasta的区别:

fasta: ATGTCNTGTCG
efasta:ATGTC{A,T}GTCG

efasta使结果更明确,但是在数据处理的时候,不易兼容。