使用DBG2OLC对二、三代混合数据进行基因组组装

1. DBG2OLC软件简介

DBG2OLC能利用二代和三代混合数据组装大基因组。其文章于2016年发表在Scientific Reports上。

2. DBG2OLC软件下载与安装

使用git下载DBG2OLC软件

$ cd /opt/biosoft/
$ git clone https://github.com/yechengxi/DBG2OLC.git
$ cd /opt/biosoft/DBG2OLC
按照说明中对软件进行编译,编译出的3个可执行程序全部都是DBG2OLC命令
$ g++ -O3 -o SparseAssebmler DBG2OLC.cpp
$ g++ -O3 -o DBG2OLC *.cpp
$ g++ -O3 -o Sparc *.cpp
直接拷贝作者编译好的程序即可
$ chmod 755 compiled/*
$ cp compiled/* .
$ echo 'PATH=$PATH:/opt/biosoft/DBG2OLC' >> ~/.bashrc
$ source ~/.bashrc

DBG2OLC程序第三步需要blasr, sparc/pbdagcon软件支持。其中sparc在DBG2OLC安装文件夹下。
安装blasr

下载BLASR
$ git clone https://github.com/PacificBiosciences/blasr.git /opt/biosoft/blasr
$ cd /opt/biosoft/blasr/
下载libcpp和pbbam两个submodules
$ make update-submodule

blasr编译需要hdf5支持,从hdf5官网下载适合centos6的二进制包并安装 
$ wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/bin/linux-centos6-x86_64-gcc447/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz
$ tar zxf hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz -C /opt/sysoft/

可以使用cmake, pitchfork和make三种方式对BLASR进行编译,以下使用常规的make方法进行编译,需要高版本gcc支持
对BLASR进行configure
$ ./configure.py --shared --sub --no-pbbam HDF5_INCLUDE=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/include/ HDF5_LIB=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/
对submodules进行configure
$ make configure-submodule
对submodules进行make
$ make build-submodule -j 4
对BLASR进行make
$ make blasr -j 4
对其它工具,例如pls2fasta, loadPulses, sawriter等进行编译,其结果文件在utils文件夹中
$ make -j 4

可选步骤:手动将有用的命令和库文件放置到指定的地方
blasr的正常运行需要依赖libcpp里面三个库文件和hdf5软件中的库文件
$ mkdir bin lib
$ cp blasr bin/
$ find utils -maxdepth 1 -perm 775 -exec cp {} bin/ \;
$ cp ./libcpp/pbdata/libpbdata.so ./libcpp/hdf/libpbihdf.so ./libcpp/alignment/libblasr.so lib/
$ echo 'export LD_LIBRARY_PATH=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/:/opt/biosoft/blasr/lib/:$LD_LIBRARY_PATH
PATH=/opt/biosoft/blasr/bin/:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

若DBG2OLC流程第三步选择使用pbdagcon进行运算,则需要安装pbdagcon软件

pbdagcon软件的编译需要高版本gcc支持
$ git clone https://github.com/PacificBiosciences/pbdagcon.git /opt/biosoft/pbdagcon
$ cd /opt/biosoft/pbdagcon
$ ./configure.py --boost --gtest --sub --no-pbbam
$ make init-submodule
$ make -j 4
$ make check
$ mkdir bin
$ cp src/cpp/dazcon src/cpp/pbdagcon bin/
$ echo 'PATH=/opt/biosoft/pbdagcon/bin:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

3. 程序运行

使用DBG2OLC软件利用二代和三代数据混合的基因组组装,其运行流程分3步。

3.1 使用SparseAssembler利用二代数据进行DBG组装

首先,利用Illumina小片段文库数据使用SparseAssembler命令组装出contigs序列。此外,也可以使用其他基因组组装软件组装出contigs序列后,直接跳到DBG2OLC的第二个步骤。值得注意的是:输入到第二步骤的contigs必须是没有经过repeat resolving的原始序列;绝大部分基因组组装软件为了获得更完整连续的contigs序列,都牺牲了contigs的准确性,其结果不能用于DBG2OLC软件的第二步,否则最终结果会很差;作者推荐可以直接用于第二步的其它contig组装软件有Platanus和Meraculous。
一般情况下,输入到SparseAssembler命令中~50x的Illumina数据,能获得较好的contigs结果。
SparseAssembler命令参数:

常用参数:
LD 
    是否载入k-mer graph。第一次运行SparseAssembler命令的时候,该参数的值必须是0;若为了使用SparseAssembler得到更好的contigs结果,则需要调整NodeCovTh和EdgeCovTh参数;调整这些参数的时候,不需要再次计算k-mer graph,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置使用DBG方法计算时的Kmer大小,支持的Kmer大小为15-127。
g 
  number of skipped intermediate k-mers, support 1-25.该参数软件示例中使用的值是15。
NodeCovTh 
  设置k-mers覆盖度阈值,去除覆盖度较低的k-mers。该值设定范围为0-16,默认值为1。
EdgeCovTh 
  设置link覆盖度阈值,去除覆盖度较低的links。该值设定范围为0-16,默认值为0。
GS 
  设置一个基因组大小的值。该参数用于决定预先占用的内存量。推荐设置得比基因组大,例如设置为2倍基因组大小。
f 
  输入单端测序数据的路径。输入文件可以是fasta或fastq文件。若有多个输入文件,则使用多个f参数。
i1  i2 
  输入inward paired-end数据。若有多组paired-end数据,则多次使用i1/i2参数。

其它参数:
o1  o2 
  输入outward paired-end数据。
i1_mp  i2_mp 
  输入插入片段长度>10kb的inward paired-end数据。
o1_mp  o2_mp 
  输入插入片段长度>10kb的outward paired-end数据。
PathCovTh 
  设置path覆盖度阈值,去除覆盖度较低的paths。该值设定范围为0-100。根据经验,不推荐添加该参数。
TrimLen 
  将所有过长的序列截短到此指定的长度。
TrimN 
  若read中的碱基N数目超过此设定的值,则去除该read数据。
TrimQual 
  从尾部截短质量低于此值的碱基。
QualBase 
  设置Fastq文件中最低碱基质量对应的ASCII码符号。默认值是'!',等同于Pred33。
Denoise 
  设置是否对reads进行修正。默认值是0,表示不对reads进行修正。
H 
  混合模式。默认值是0,表示对reads的尾部进行截短,来保证高质量的数据进行reads修正。
CovTh 
  覆盖度 < 此设定值的k-mer会被检测,从而被校正。若该参数值设置为0,则软件会自动计算该值。
CorrTh 
  覆盖度 >= 此设定值的k-mer可以用来对reads做校正。若该参数值设置为0,则软件会自动计算该值。

SparseAssembler运行示例:

对某物种Illumina小片段文库测序的PE150bp数据使用trimmomatic质控,再使用FindErrors进行修正,再运行SparseAssembler:
$ SparseAssembler LD 0 k 95 g 15 NodeCovTh 1 EdgeCovTh 0 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq
$ cp Contigs.txt Contigs.txt.00
增大NodeCovTh和EdgeCovTh参数后,再次运行SparseAssembler,并比较两次结果。第二次运行较第一次运行,耗时少了很多很多。
$ SparseAssembler LD 1 k 95 g 15 NodeCovTh 2 EdgeCovTh 1 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq

SparseAssembler在当前目录下生成了18个文件结果,其中Contigs.txt文件是Fasta格式的Contigs序列。
运行SparseAssembler的注意事项:

1. SparseAssembler只可以简单地对Fastq文件进行质控和错误修正。推荐使用其它软件进行reads质控和修正,以获得更好的结果。
2. 参数k设置了k-mer的大小,该参数的值对结果的影响较大。若基因组较小,推荐设置多个k-mer值进行多次计算,从而选择最优k-mer值。个人经验,PE150bp数据的最优的k-mer值约为91~99。
3. 选定了k-mer大小后,使用默认的NodeCovTh和EdgeCovTh参数(默认参数一般能得到很好的结果)运行一遍SparseAssembler。然后尝试增大NodeCovTh和EdgeCovTh参数,设置LD 1参数再次运行SparseAssembler,以获得最优的Contigs结果。
4. 可能是先使用了最小的NodeCovTh和EdgeCovTh参数做运算后,才能再次使用更大的参数进行运算。
5. SparseAssembler虽然也有输入大片段文库数据的参数和Scaffolding参数,但是不推荐输入大片段文库数据进行Scaffolding操作,没太大意义。
6. 虽然SparseAssembler命令的文件输入方式有多种,若是仅进行contigs组装,没有利用到paired信息,因此使用i1 i2参数输入文件和使用f参数输入文件的结果是一模一样的。

3.2 使用DBG2OLC找Contigs序列和Pacbio reads的Overlap并进行Layout

DBG2OLC通过比较contigs和Pacbio reads之间的overlap,将contigs序列定位到Pacbio reads上,将DBG的contigs结果运用到OLC算法中。
DBG2OLC命令参数:

主要参数:
LD 
  是否载入compressed reads information。第一次运行DBG2OLC命令的时候,该参数的值必须是0;若为了得到更好的结果,则需要调整其它参数;调整这些参数的时候,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置k-mer大小。k-mer用来比较contig和pacbio read之间的重叠,而不是用于基因组组装,推荐设置为 17 即可。
AdaptiveTh 
  若contig和pacbio read之间匹配的k-mers数目 < AdaptiveTh * contig长度,则认为contig和pacbio read没有重叠。推荐设置为0.001-0.02。
KmerCovTh 
  若contig和pacbio read之间匹配k-mers的覆盖度 < KmerCovTh,则认为contig和pacbio read没有重叠。推荐设置为2-10。
MinOverlap 
  两条Pacbio read之间匹配的k-mers数目 < MinOverlap,则认为它们之间没有重叠。推荐设置为10-150。
RemoveChimera 
  去除嵌合体Pacbio reads。若Pacbio数据覆盖度大于10x,推荐设置该参数为 1 。
Contigs 
  输入contigs序列文件路径
f 
  输入Pacbio测序Fasta/Fastq文件路径。

其它参数:
MinLen 
  设置能用于计算的最小Pacbio reads长度。
ChimeraTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。
ContigTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。

DBG2OLC运行示例:

$ DBG2OLC LD 0 k 17 AdaptiveTh 0.001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq
$ DBG2OLC LD 1 k 17 AdaptiveTh 0.005 KmerCovTh 3 MinOverlap 30 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq

DBG2OLC的结果文件很多,其主要结果文件是backbone_raw.fasta和DBG2OLC_Consensus_info.txt,是第三步的输入文件。
运行DBG2OLC的注意事项:

1. AdaptiveTh, KmerCovTh和minOverlap这3个计算Overlap的参数对结果的影响最大。对于10x/20x PacBio数据:KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01;对于50x-100x PacBio数据:KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02。
2. 不推荐对Pacbio数据就行修正后再运行DBG2OLC。可以比较使用修正前的数据用于DBG2OLC的结果,一般情况下使用未修正的Pacbio数据能获得更好的结果。此外,DBG2OLC运行过程中有一步多序列比对模块来进行错误修正。
3. 可能是先使用了最小的AdaptiveTh, KmerCovTh和minOverlap参数做运算后,才能再次使用更大的参数进行运算。

3.3 Call consensus

本步骤需要使用/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本,来调用blasr和consensus模块Sparc(可以考虑使用pbdagcon作为consensus模块,但DBG2OLC没有提供相应的脚本)进行运算。

先将/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本拷贝到当前目录
$ cp /opt/biosoft/DBG2OLC/utility/*.sh /opt/biosoft/DBG2OLC/utility/*.py ./
若使用了最新版本的blasr软件,其参数书写方法有一个中划线变成了两个中划线,因此需要修改.sh文件中blasr命令的参数书写方法。
此外,也需要修改.sh文件中Sparc命令路径,或者将Sparc命令也拷贝到当前目录。

将Contigs序列和Pacbio reads数据合并成一个文件ctg_pb.fasta
$ cp Contigs.txt ctg_pb.fasta
$ perl -e 'while (<>) {print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {print; $_ = <>; print; <>; <>;}' Pacbio_Cell02.fastq >> ctg_pb.fasta

运行脚本程序split_and_run_sparc.sh
$ ./split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta ./ 2 > cns_log.txt
结果会输出到 ./ 目录下。最后的结果文件是final_assembly.fasta。

AGP格式简单说明

AGP文件为NCBI数据上传要求的标准格式,用来描述小片段序列(比如contig)如何构成大片段序列(比如scaffold和chromosome)。详细的说明文档请见:http://www.ncbi.nlm.nih.gov/projects/genome/assembly/agp/AGP_Specification.shtml
AGP文件有9列,分别是:

1. 大片段的序列名(object)
2. 大片段起始(object_begin)
3. 大片段结束(object_end)
4. 该段序列在大片段上的编号(part_number)
    一般一个大片段由多个小片段和gap组成。此处则为这些小片段和gap在大片段上的编号。
5. 该段序列的类型(component_type)
    常用的是W、N和U。W表示WGS contig;N表示指定大小的gap;U表示不明确长度的gap,一般用100bp长度。
6. 小片段的ID或gap长度(component_id or gap_length)
    如果第5列不为N或U,则此列为小片段的ID。
    如果第5列是N或U,则此列为gap的长度。如果第5列为U,则此列值必须为100。
7. 小片段起始或gap类型(component_begin or gap_type)
    如果第5列是N或U,则此列表示gap的类型。常用的值是scaffold,表示是scaffold内2个contigs之间的gap。其它值有:contig,2个contig序列之间的unspanned gap,这样的gap由于没有证据表明有gap,应该要打断大片段序列;centromere,表示中心粒的gap;short_arm,a gap inserted at the start of an acrocentric chromosome;heterochromatin,a gap inserted for an especially large region of heterochromatic sequence;telomere,a gap inserted for the telomere;repeat,an unresolvable repeat。
8. 小片段结束或gap是否被连接(component_end or linkage)
    如果第5列是N或U,则此列一般的值为yes,表示有证据表明临近的2个小片段是相连的。
9. 小片段方向或gap的连接方法(orientation or linkage_evidence)
    如果第5列不为N或U,则此列为小片段的方向。其常见的值为 +、-或?。
    如果第5列是N或U,则此列表明临近的2个小片段能连接的证据类型。其用的值是paired-ends,表明成对的reads将小片段连接起来。其它值有:na,第8列值为no的时候使用;align_genus,比对到同属的参考基因组而连接;align_xgenus,比对到其它属的参考基因组而连接;align_trnscpt,比对到同样物种的转录子序列上;within_clone,gap两边的序列来自与同一个clone,但是gap没有paired-ends跨越,因此这种连接两边小片段无法确定方向和顺序;clone_contig,linkage is provided by a clone contig in the tiling path (TPF);map,根据连锁图,光学图等方法确定的连接;strobe,根据PacBio序列得到的连接;unspecified。如果有多中证据,则可以写上多种证据,之间用分号分割。

例子:
Scaffold from component (WGS)
Chromosome from scaffold (WGS)

使用 abyss 进行 scaffolding

1. ABySS 进行 scaffolding 的目的与优点

目的: 对其它基因组 denovo 的 assembly 结果,使用 abyss 再进行一次 scaffolding。
优点: 可以使用 RNA-seq 的转录子数据进行基因组的辅助组装。

2. ABySS 进行 scaffolding 的命令行

输入文件: assembly.fasta, 2000.1.fastq, 2000.2.fastq, 5000.1.fastq, 5000.2.fastq。
输入文件是基因组的组装结果,和 3 对 mate-paired Illumina 数据。

2.1 对 assembly.fasta 进行序列改名

去除序列之间的换行
fasta_no_blank.pl assembly.fasta > 11; mv 11 assembly.fasta
给序列按顺序重命名
perl -e '$num = 0; while (<>) {if (/^>/) { s/>(.*)/>$num/; print; $num ++; } else { print } }' assembly.fasta > ledodes-6.fa

2.2 将 mate-paired 数据比对到基因组序列上

根据比对结果,得到 mate-paired library 的 insertSize 信息(以 .hist 为后缀的文件)和 序列之间的连接、距离与顺序信息 (以 .dist.dot 为后缀的 graph 文件)。

abyss-map -j24 -l87 3000.1.fastq 3000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp1-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp1-6.dist.dot mp1-6.hist
abyss-map -j24 -l87 8000.1.fastq 8000.2.fastq ledodes-6.fa \
  |abyss-fixmate -l87 -h mp2-6.hist \
  |sort -snk3 -k4 \
  |DistanceEst --dot -j24 -k87 -l87 -s200 -n10 -o mp2-6.dist.dot mp2-6.hist

以上命令行中参数:

-j24
    使用 24 个线程运行
-l87
    使用的 kmer值 为 87
-s200
    sedd contigs的最小长度为 200bp
-n10
    所允许连接两条序列的最小的pairs的数目

2.3 进行 scaffolding

abyss-scaffold -k87 -s200 -n5 -g ledodes-6.path.dot ledodes-6.fa mp1-6.dist.dot mp2-6.dist.dot > ledodes-6.path
PathConsensus -k87 -p0.9 -s ledodes-7.fa -g ledodes-7.adj -o ledodes-7.path ledodes-6.fa ledodes-6.fa ledodes-6.path
cat ledodes-6.fa ledodes-7.fa \
  | MergeContigs -k87 -o ledodes-8.fa - ledodes-7.adj ledodes-7.path
ln -sf ledodes-8.fa ledodes-scaffolds.fa
PathOverlap --overlap --dot -k87 ledodes-7.adj ledodes-7.path > ledodes-8.dot

2.4 使用转录子序列进行 rescaffolding

bwa index ledodes-8.fa
bwa mem -a -t2 -S -P -k87 ledodes-8.fa transcripts.fasta \
  |gzip > long1-8.sam.gz
abyss-longseqdist -k87 long1-8.sam.gz \
  |grep -v "l=" >long1-8.dist.dot
abyss-scaffold -k87 -s200 -n1 -g ledodes-8.path.dot ledodes-8.dot long1-8.dist.dot > ledodes-8.path
PathConsensus -k87 -p0.9 -s ledodes-9.fa -g ledodes-9.adj -o ledodes-9.path ledodes-8.fa ledodes-8.dot ledodes-8.path
cat ledodes-8.fa ledodes-9.fa \
  | MergeContigs -k87 -o ledodes-10.fa - ledodes-9.adj ledodes-9.path
ln -sf ledodes-10.fa ledodes-long-scaffs.fa

使用 GCE 进行基因组大小评估

1. GCE 简介

GCE(Genome Characteristics Estimation) 是华大基因用于基因组评估的软件,其参考文献为:Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects。下载地址:ftp://ftp.genomics.org.cn/pub/gce

GCE 软件包中主要包含 kmer_freq_hash 和 gce 两支程序。前者用于进行 kmer 的频数统计,后者在前者的结果上进行基因组大小的准确估算。

2. GCE 下载和安装

$ wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
$ tar zxf gce-1.0.0.tar.gz -C /opt/biosoft

3. kmer_freq_hash 的使用

kmer_freq_hash 的常用例子:

$ /opt/biosoft/gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
  -k 21 -l reads.list -a 10 -d 10 -t 24 -i 50000000 -o 0 -p species &> kmer_freq.log

kmer_freq_hash 的常用参数:

-k <17>
    设置 kmer 的大小。该值为 9~27,默认值为 17 。
-l string
    list文本文件,其中每行为一个fastq文件的路径。
-t int
    使用的线程数,默认为 1 。
-i int
    初始的 hash 表大小,默认为 1048576。该值最好设置为 (kmer 的种类数 / 0.75)/ 线程数。如果基因组大小为 100M,测序了 40M 个 reads,reads 的长度为 100bp,测序错误率为 1%,kmer的大小为 21,则kmer的种类数为100M+40M*100*1%*21=940M,若使用24线程,则该参数设置为 i=940M/0.75/24=52222222。
-p string
    设置输出文件的前缀。
-o int
    是否输出 k-mer 序列。1: yes, 0: no,默认为 1 。推荐选 0 以节约运行时间。
-q int
    设置fastq文件的phred格式,默认为 64。该值可以为 33 或 63。
-c double
    设置k-mer最小的精度,该值位于 0~0.99,或为 -1。 -1 表示不对 kmer进行过滤。设置较高的精度,可以用于过滤低质量 kmer。精度是由 phred 格式的碱基质量计算得来的。
-r int
    设置获取 k-mer 使用到的 reads 长度。默认使用 reads 的全长。
-a int
    忽略read前面该长度的碱基。
-d int
    忽略read后面该长度的碱基。
-g int
    设置使用该数目的碱基来获取 k-mers,默认是使用所有的碱基来获取 k-mer。

kmer_freq_hash 的主要结果文件为 species.freq.stat。该文件有 2 列:第1列是kmer重复的次数,第二列是kmer的种类数。该文件有255行,第225行表示kmer重复次数>=255的kmer的总的种类数。该文件作为 gce 的输入文件。
kmer_freq_hash 的输出到屏幕上的信息结果保存到文件 kmer_freq.log 文件中。该文件中有粗略估计基因组的大小。其中的 Kmer_individual_num 数据作为 gce 的输入参数。

4. gce 的使用

gce 的常用例子:

$ /opt/biosoft/gce-1.0.0/gce \
  -f species.freq.stat -c 85 -g 4112118028 -m 1 -D 8 -b 1 > species.table 2> species.log

gce 的常用参数:

-f string
    kmer depth frequency file
-c int
    kmer depth frequency 的主峰对应的 depth。gce 会在该值附近找主峰。
-g int
    总共的 kmer 数。一定要设定该值,否则 gce 会直接使用 -f 指定的文件计算 kmer 的总数。由于默认下该文件中最大的 depth 为 255,因此,软件自己计算的值比真实的值偏小。同时注意该值包含低覆盖度的 kmer。
-M int
    支持最大的 depth 值,默认为 256 。
-m int
    估算模型的选择,离散型(0),连续型(1)。默认为 0,对真实数据推荐选择 1 。
-D int
    precision of expect value,默认为 1。如果选择了 -m 1,推荐设置该值为 8。
-H int
    使用杂合模式(1),不使用杂合模式(0)。默认值为 0。只有明显存在杂合峰的时候,才选择该值为 1 。
-b int
    数据是(1)否(0)有 bias。当 K > 19时,需要设置 -b 1 。

gce 的结果文件为 species.table 和 species.log 。species.log 文件中的主要内容:

raw_peak	now_node	low_kmer	now_kmer	cvg	genome_size	a[1]	b[1]
84	35834245	22073804	4044916750	84.6637	4.83093e+07	0.928318	0.637648

raw_peak: 覆盖度为 84 的 kmer 的种类数最多,为主峰。
now_node: kmer的种类数。
low_kmer: 低覆盖度的 kmer 数。
now_kmer: 去除低覆盖度的 kmer 数,此值 = (-g 参数指定的总 kmer 数) - low_kmer 。
cvg:估算出的平均覆盖度
genome_size:基因组大小,该值 = now_kmer / cvg 。
a[1]: 在基因组上仅出现 1 次的 kmer 之 种类数比例。
b[1]: 在基因组上仅出现 1 次的 kmer 之 数量比例。该值代表着基因组上拷贝数为 1 的序列比例。

如果使用 -H 1 参数,则会得额外得到如下信息:

for hybrid: a[1/2]=0.223671 a1=0.49108
kmer-species heterozygous ratio is about 0.125918

上面结果中,0.125918 是由 a[1/2] 计算出来的。 0.125918 = a[1/2] / ( 2- a[1/2] ) 。
a[1/2]=0.223671 表示在所有的 uniqe kmer 中,有 0.223671 比例的 kmer 属于杂合 kmer 。

此外,有 a[1/2] 和 b[1/2] 的值在最后的统计结果中。重复序列的含量 = 1 - b[1/2] - b[1] 。

则杂合率 = 0.125918 / kmer_size 。 若计算出的杂合率低于 0.2%,个人认为测序数据应该是纯合的。这时候,应该不使用 -H 1 参数。使用 -H 1 参数会对基因组的大小和重复序列含量估算造成影响。

5. 不同杂合率,有无重复序列的 kmer species 和 kmer individuals 图

下图中 a 和 b 是对理想中无重复的基因组在不同杂合率下的曲线图;
下图中 c 和 d 是对有重复的基因组(human)在不同杂合率下的曲线图。
从下图可以参考不同杂合率下的曲线状况。
kerm_pictures

使用FGAP进行补洞

1. FGAP简介

FGAP利用BLAST将contigs序列比对到基因组草图序列上,寻找重叠到gap区间的最优序列,从而进行补洞。其参考文献:Piro, Vitor C., et al. “FGAP: an automated gap closing tool.” BMC research notes 7.1 (2014): 371.

2. FGAP下载和安装

FGAP官网:http://www.bioinfo.ufpr.br/fgap/

$ wget http://sourceforge.net/projects/fgap/files/MCR_LINUX64b.tar.gz/download
$ tar zxf MCR_LINUX64b.tar.gz
$ cd MCR_LINUX64b
$ ./installMCR.sh /opt/biosoft/MCR

$ wget http://sourceforge.net/projects/fgap/files/FGAP_1_7_LINUX64b.tar.gz/download
$ tar zxf FGAP_1_7_LINUX64b.tar.gz -C /opt/biosoft/

2. FGAP的使用

FGAP的简单使用示例:

$ ln -s /opt/biosoft/FGAP_1_7_LINUX64b/sample_data/* .
$ /opt/biosoft/FGAP_1_7_LINUX64b/run_fgap.sh /opt/biosoft/MCR/v717/ \
 -d DRAFT_ecoli_hiseq454.fasta \
 -a "DATASET_ecoli_454.fasta,DATASET_ecoli_hiseq.fasta"

FGAP的使用参数:

-d /--draft-file	Draft genome file [fasta format - Ex: "draft.fasta"]
-a /--datasets-files	List of datasets files to close gaps [fasta format - Ex: "dataset1.fasta,dataset2.fasta"]

-s /--min-score		Min Score (raw) to return results from BLAST (integer) - Default: 25
-e /--max-evalue	Max E-Value to return results from BLAST (float) - Default: 1e-7
-i /--min-identity	Min identity (%) to return results from BLAST (integer [0-100]) - Default: 70

-C /--contig-end-length	Length (bp) of contig ends to perform BLAST alignment (integer) - Default: 300
-T /--edge-trim-length	Length of ignored bases (bp) upstream and downstrem of the gap (integer) - Default: 0
-R /--max-remove-length	Max number of bases (bp) that can be removed (integer) - Default: 500
-I /--max-insert-length	Max number of bases (bp) that can be inserted (integer) - Default: 500

-p /--positive-gap	Enable closing of positive gaps (with insertion) (integer [0-1]) - Default: 1
-z /--zero-gap		Enable closing of zero gaps (without insert any base) (integer [0-1]) - Default: 0
-g /--negative-gap	Enable closing of negative gaps (overlapping contig ends) (integer [0-1]) - Default: 0

-c /--gap-char				Base that represents the gap (char) - Default: "N"
-b /--blast-path			Blast+ package path (only makeblastdb and blastn are needed, version 2.2.28+ or higher) - Default: ""
-l /--blast-alignment-parameters	BLAST alignment parameters (opengap,extendgap,match,mismatch,wordsize) - Default: "1,1,1,-3,15"
-r /--blast-max-results			Max results from BLAST for each query (integer) - Default: 200
-t /--threads				Number of threads (integer) - Default: 1

-m /--more-output	More output files with gap regions after and before gap closing (integer [0-1]) - Default: 0
-o /--output-prefix	Output prefix [File or folder - Ex: "out" or "out/" ] - Default: "output_fgap"
-h /--help		This help message

使用ABACAS进行有参考基因组的contigs连接

1. ABACAS简介

ABACAS 用于在有参考序列的情况下,对contigs序列进行连接。软件官网:http://abacas.sourceforge.net/

2. ABACAS 下载和安装

ABACAS是一支perl程序,其运行需要安装有MUMmer。

$ wget http://sourceforge.net/projects/abacas/files/abacas.1.3.1.pl
$ chmod 755 abacas.1.3.1.pl
$ mv abacas.1.3.1.pl ~/bin/
$ ln -s ~/bin/abacas.1.3.1.pl ~/bin/abacas.pl
$ perl -p -i -e 's#/usr/local/bin/perl#/usr/bin/perl#' ~/bin/abacas.1.3.1.pl

3. ABACAS 的使用

ABACAS 的简单例子:

$ abacas.pl -r ref.fasta -q query.fasta -p nucmer -o out_prefix

ABACAS 的使用参数:

-r string
    输入参考序列。参考序列为fasta格式并仅有1条序列。
-q string
    输入query序列。query序列为fasta格式。
-p nucmer|promer
    该参数可以设为nucmer或promer,表示使用相应的程序进行contig ordering。
-o string
    设定输出文件的前缀。默认为query.fasta_ref.fastq。
-c
    参考序列是一个环。
-m
    使用该参数后,则将定位并定向的query序列按顺序打印到文件 out_prefix.MULTIFASTA.fa 中。
-b
    使用该参数后,则将未能定位的 congtig 序列内容打印到文件 ou_prefix.contigsInbin.fas 中。
-N
    使用该参数后,则将由定位定向的conitg序列连接而成的pseudomolecule的序列打印到文件 out_prefix.NoNs.fasta 中。默认下,输出文件 out_prefix.fasta 中的 pseudomolecule 序列中 congtig 之间的 gap 使用 N 表示,有重叠则加 100个N;该参数则是生成额外的fasta文件,其中contig之间序列无N。
-t
    使用该参数后,则对 ou_prefix.contigsInbin.fas 中的序列运行 blastall,额外生成 blast.out 文件。使用该参数则必须有 -b 参数。
-g out_prefix
    使用该参数后,则将 pseudomolecule 的gap部分在ref上的序列答应到文件 out_prefix.Gaps_onRef 中。
-P
    使用该参数后,则设计用于 close gaps 的引物。需要安装 prime3_core,但我个人加该参数运行失败。
-R
    使用该参数,则不需要运行 mumer,适合于多次运行程序,节约时间。但我个人加该参数运行失败。
-e
    不进行 contig ordering,有利于直接进行引物设计。
-i int
    允许的最小的identity的百分比。默认为 40。
-v int
    允许的最小的 contig 覆盖率。默认为 40。
-V 0|1
    contig比对到ref序列上的位置的数目。默认为1。使用 -V 0 表示将contig序列随机放到1个比对上的位点。
-l int
    不使用低于此长度的contig序列。默认为 100。 
-d
    ABACAS在调用MUmer的时候,默认下使用了--maxmatch参数来增加匹配,使用 -d 参数则会使用MUmer的默认参数,即不会使用--maxmatch参数。

使用 REAPR 进行基因组错误评估

1. REAPR 简介

REAPR(Recognition of Errors in Assemblies using Paired Reads)能利用成对的reads来识别基因组序列中的错误。从而,能将基因组序列从错误的gap处断开或将错误序列使用 Ns 代替。同时,对错误信息进行统计。

2. REAPR下载和安装

REAPR官网:http://www.sanger.ac.uk/resources/software/reapr/
安装 REAPR 需要先安装 R 和 Perl 模块: File::Basename, File::Copy, File::Spec, File::Spec::Link, Getopt::Long, List::Util。

$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.tar.gz
$ tar zxf Reapr_1.0.17.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/Reapr_1.0.17
$ wget ftp://ftp.sanger.ac.uk/pub/resources/software/reapr/Reapr_1.0.17.manual.pdf
$ ./install.sh
$ echo 'PATH=$PATH:/opt/biosoft/Reapr_1.0.17' >> ~/.bashrc
$ source ~/.bashrc

测试Reapr能否正常运行:
$ wget ftp://ftp.sanger.ac.uk/pub4/resources/software/reapr/Reapr_1.0.17.test_data.tar.gz
$ tar zxf Reapr_1.0.17.test_data.tar.gz 
$ ./test.sh

3. REAPR 使用

reapr 最少需要输入基因组fasta文件和mate-pair数据(需要将RF方向转变成FR方向)文件。常用的运行步骤如下:

使用 facheck 命令检查基因组fasta文件,以防fasta文件序列名含有怪异字符。
$ reapr facheck genome.fasta
使用下列命令对fasta文件进行修改。
$ reapr facheck genome.fasta new_genome.fasta

使用SMALT将 mate-pair reads比对到基因组上。
$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam

运行reapr pipeline对基因组进行评估。将结果输出到 outdir 文件夹中。
$ reapr pipeline genome.fasta jumping.bam outdir

reapr 也支持输入paired-end数据(可选)。通过将paired-end数据uniquely并perfecly比对到基因组,从而能计算出基因组上error-free的碱基位点并进行统计。该项目分析对error calling没有影响。值得注意的是: 默认设置下call一个error-free位点需要至少5x的覆盖度;同时输入的paired-end数据的方向要为inward(一般paired-end数据方向即为inward),此外reapr仅支持inward方向,因此输入的jumping数据也要先转换为inward方向。

使用paired-end数据的命令行:

$ reapr smaltmap genome.fasta jumping.1.fastq jumping.2.fastq jumping.bam
$ reapr perfectmap genome.fasta frag.1.fastq frag.2.fastq insert_size perfect_prefix
$ reapr pipeline genome.fasta jumping.bam outdir perfect_prefix

reapr的运行流程:
reapr_pipeline

4. reapr的结果

REAPR 会对基因组每个位点的 FCD(fragment coverage distribution)进行分析。期望的FCD和实际上的上FCD之差为”FCD error”。 FCD error往往代表不正确的scaffold连接、大的插入缺失或contig序列中不正确的连接。为了能更好地计算FCD,特别是gap位点的FCD,mate-pair数据的insert size越大越好。

REAPR对输入的paired reads数据进行检测,并将结果放入到文件夹 outdir/00.Sample/ 下。

gc_vs_cov.lowess.pdf: GC偏好性检测图
insert.in.pdf: FR(inner/inward)类型插入片段的长度检测图
insert.out.pdf: RF类型插入片段的长度检测图
insert.stats.txt: 插入片段长度统计

从错误位点打断后的基因组序列: outdir/04.break.broken assembly.fa 。
如果FCD error区有gap,则从gap处将基因组序列打断;若FCD error区没有gap,则将其用Ns代替。同时将被代替的序列写入到文件 outdir/04.break.broken assembly bin.fa中。

基因组错误统计文件: outdir/05.summary.report.txt。该文件统计位于 contig 中的 FCD error 的数目,和含有 gap 的 FCD error 的数目。

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使结果更明确,但是在数据处理的时候,不易兼容。