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 (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; 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。
