使用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 (<>) {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。

匿名进行回复 取消回复

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据