lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

三代测序基因组组装工具MECAT2/NECAT

已有 2784 次阅读 2019-10-30 21:31 |个人分类:Software/Pipeline/Script|系统分类:科研笔记| _生信工具, _基因组组装

三代测序基因组组装工具MECAT2/NECAT

 

MECAT是一款优秀的基因组三代组装工具,前段时间发现MECAT已经更新到第二版了。相较于初版,修复了一些bug,并改进了算法,运行效率更高。并且,原版的MECAT处理PacBioNanopore数据时使用同一程序,而改进版细分为了MECAT2(用于PacBio数据)和NECAT(用于Nanopore数据)。然后就顺便试了试。

对于MECAT软件本身的介绍,这里就跳过了,毕竟网上对它的介绍已经有很多了。而且相信大多数同学应该跟我一样,懒得关注算法,好吧最主要的还是根本看不懂啊2.gif;也不是很在乎它和其它组装软件的优劣,网上关于软件的评估测试也有很多,但通常各说各的好,咱们一般也只当参考,实际用的时候还不是多个工具都试下,毕竟结果谁好谁差还真不一定啊,适合自己数据的工具才是好工具。

那么最直接的关注点首先肯定是软件怎么使用了,下文主要就是使用测试,为了方便快速得到结果,使用的细菌小基因组数据(好吧其实本人也是对高等生物的基因组组装产生严重的心理阴影了……)。

测试使用的PacBioNanopore数据下载自SRA

https://www.ncbi.nlm.nih.gov/sra?term=SAMN12698632

我也保存在网盘中:

https://pan.baidu.com/s/1sCtbbvGW8hxe40p1R11DKQ

 


MECAT2


新版MECAT2github链接: https://github.com/xiaochuanle/MECAT2

github中已经提供了详细的文档供我们参考。

 

MECAT2安装


根据github文档,有两种方式安装MECAT2,一是直接下载已经编译好的二进制程序使用,二是源码编译。

1.png

#本人犯懒,直接下载二进制版本
wget https://github.com/xiaochuanle/MECAT2/releases/download/20192026/mecat2_20190226_linuax_amd64.tar.gz
tar xzvf mecat2_20190226_linuax_amd64.tar.gz
 
#之后将主程序添加至环境变量里,例如我的 MECAT2 存放路径在 /home/my/software/MECAT2/
export PATH=/home/my/software/MECAT2/Linux-amd64/bin:$PATH

 

MECAT2组装PacBio测序数据


接下来,开始测试组装基因组的过程。

升级版的MECAT2使用起来非常便捷。组装过程分为两步,第一步为reads纠错,第二步为基因组组装,整理好输入的配置文件运行就可以了。

#查看帮助
mecat.pl

2.png

 

编辑配置文件

github文档中给定了config配置文件的参考样式,我们可仿照它的样式来。

#创建路径,编辑 config 文件
mkdir mecat2 && cd mecat2
vi mecat_cfg.txt
 
#config 格式可参考 https://github.com/xiaochuanle/MECAT2

例如在这里,对于我们的测试物种,config配置文件“mecat_cfg.txt”编辑如下。编辑完成后,关闭保存config配置文件。

PROJECT=Paenibacillus_sp.R4
RAWREADS=/home/my/data/Paenibacillus_sp.R4/SRR10074456_subreads.fastq
GENOME_SIZE=8600000
THREADS=4
MIN_READ_LENGTH=2000
CNS_OVLP_OPTIONS=""
CNS_OPTIONS="-r 0.6 -a 1000 -c 4 -l 2000"
TRIM_OVLP_OPTIONS="-B"
ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400"
FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1"
FSA_ASSEMBLE_OPTIONS=""
USE_GRID=false
GRID_NODE=0
CLEANUP=0
CNS_OUTPUT_COVERAGE=60

PROJECT:作业名称,随便命个名,会在结果路径中生成以其为命名的主目录。

RAWREADSPacBio数据路径,fastqfasta格式,不能是压缩版的。

GENOME_SIZE:预估基因组大小,示例8600000意为8600000bp

THREADS:程序运行的线程数,示例为4

MIN_READ_LENGTHreads纠错(corrected reads)以及修剪(trimmed reads)的最低长度,示例2000意为2000bp

CNS_OVLP_OPTIONS:用于在reads纠错阶段检测候选重叠的选项,参数将传递给子程序mecat2pw执行;详见mecat2pw -h

CNS_OPTIONS:纠错原始reads的参数选项,将传递给子程序 mecat2cns执行;详见mecat2cns -h

TRIM_OVLP_OPTIONS:用于在reads修剪阶段检测候选重叠的选项,参数将传递给子程序v2asmpm执行;详见v2asmpm -h

ASM_OVLP_OPTIONS:用于在reads组装阶段检测候选重叠的选项,参数将传递给子程序v2asmpm.sh执行;详见v2asmpm.sh -h

FSA_OL_FILTER_OPTIONS:用于过滤重叠的选项,参数将传递给子程序fsa_ol_filter执行;详见fsa_ol_filter -h

FSA_ASSEMBLE_OPTIONS:用于组装修剪后reads的选项,参数将传递给子程序v2asmpm执行;详见v2asmpm -h

USE_GRIDGRID_NODE:使用多个计算节点,正常情况下用不到(修改运行线程THREADS就足以了),USE_GRID设置为falseGRID_NODE0

CLEANUP:是否删除中间文件,这里小基因组小数据,设置为0(否)倒是无妨;大基因组还是设置为1(是)吧,不然存储空间会被严重占用。

CNS_OUTPUT_COVERAGE:以预估基因组大小为准,提取N×覆盖度的最长纠错后reads进行修剪和组装;示例为60,也就是会提取60* 8600000bp=516Mbp总长度的reads用于后续修剪和组装;该参数传递给子程序mecat2elr,详见mecat2elr -h

 

程序运行,reads校正和基因组组装

然后按顺序自动执行两步就可以了。

#reads 纠错
mecat.pl correct mecat_cfg.txt
 
#contig 组装
mecat.pl assemble mecat_cfg.txt

对于reads纠错,是三代测序PacBio数据自我纠错的环节;后面基因组组装使用校正后的reads执行三代组装。关于校正和组装过程的算法原理,有兴趣的话还请自行了解。

这个测试用的细菌基因组数据不大,再加上MECAT本身校正算法也高效,因此运行时间也不会很长。多长时间我没注意,4线程运行大概2小时左右。

 

工作路径中生成了结果目录“Paenibacillus_sp.R4”。

其它中间文件略过,我们关注最终的组装结果“Paenibacillus_sp.R4/4-fsa/contigs.fasta

#查看组装得到的基因组 fasta 文件
less Paenibacillus_sp.R4/4-fsa/contigs.fasta

3.png

 注:MECAT2会自动判别组装所得序列的类型供我们参考,circular意为环状序列linear意为“线性序列”。但是记住这些仅供参考,不一定完全就是准的,后续一定要自己想办法核实下。



NECAT


接下来是NECAT,它被开发用于组装Nanopore数据。

NECATgithub链接:https://github.com/xiaochuanle/NECAT/

github中已经提供了详细的文档供我们参考。

 

NECAT安装


对于NECAT,目前作者仅为我们提供了已经编译好的二进制程序使用,暂未公开源码(或许还有待完善吧)。

#下载二进制版本
wget https://github.com/xiaochuanle/NECAT/releases/download/v0.01/necat_20190307_linux_amd64.tar.gz
tar xzvf necat_20190307_linux_amd64.tar.gz
 
#添加环境变量,例如我的 NECAT 存放路径在 /home/my/software/NECAT/
export PATH=/home/my/software/NECAT/Linux-amd64/bin:$PATH

 

NECAT组装Nanopore测序数据


接下来,开始测试组装基因组的过程。

NECAT的使用过程和上述MECAT2方法一致,非常便捷,整理好输入的配置文件运行就可以了。组装过程分为三步,第一步为reads纠错,第二步组装重叠群,第三步桥接重叠群。

#查看帮助
necat.pl

4.png

 

编辑配置文件

同样地,github文档中给定了config配置文件的参考样式,我们可仿照它的样式来。

#创建路径,编辑 config 文件
mkdir necat && cd necat
vi necat_cfg.txt
 
#config 格式可参考 https://github.com/xiaochuanle/NECAT/

例如在这里,对于我们的测试物种,config配置文件“necat_cfg.txt”编辑如下。编辑完成后,关闭保存config配置文件。

PROJECT=Paenibacillus_sp.R4
ONT_READ_LIST=/home/my/test/necat/read.list
GENOME_SIZE=8600000
THREADS=4
MIN_READ_LENGTH=3000
OVLP_SENSITIVE_OPTIONS="-n 500 -z 10 -e 0.5 -j 0 -u 1 -a 1000"
CNS_SENSITIVE_OPTIONS="-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0"
OVLP_FAST_OPTIONS="-n 500 -z 20 -b 2000 -e 0.5 -j 0 -u 1 -a 1000"
CNS_FAST_OPTIONS="-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0"
TRIM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 1 -a 400"
ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400"
NUM_ITER=2
CNS_OUTPUT_COVERAGE=60
CLEANUP=0
USE_GRID=false
GRID_NODE=0
FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1 --coverage=60"
FSA_ASSEMBLE_OPTIONS=""
FSA_CTG_BRIDGE_OPTIONS="--dump --read2ctg_min_identity=80 --read2ctg_min_coverage=4 --read2ctg_max_overhang=500 --read_min_length=5000 --ctg_min_length=1000 --read2ctg_min_aligned_length=5000 --select_branch=best"

NECAT配置文件和MECAT2差不多,相似去理解就可以了。

PROJECT:作业名称,随便命个名,会在结果路径中生成以其为命名的主目录。

ONT_READ_LIST:先将所有的Nanopore数据(解压后的fastqfasta)路径写至一文本文件(多数据时,一个路径占一行),再将该文件的路径指定在这里,read.list内容示例如下(不要问我为什么不直接使用原始的SRR10074455_1.fastq,它实在太大了,我运行过一次直接存储空间满了......然后我就从原始数据中截了一些出来,cut后的数据共取了原fastq文件的前80W行用于后续组装。我也很好奇上传者为什么测那么多Nanopore数据,应该不止拼完成图那么简单?因为只是完成图的话,后面来看只通过Pacbio数据再外加Illumina数据就完全够用了。测甲基化?Pacbio RSII的效果也比Nanopore要好吧。想不通):

/home/my/data/Paenibacillus_sp.R4/SRR10074455_1.cut.fastq

如果是多个fastq文件,逐行排列,例如:

/home/my/data/Paenibacillus_sp.R4/SRR10074455_1.cut1.fastq
/home/my/data/Paenibacillus_sp.R4/SRR10074455_1.cut2.fastq
/home/my/data/Paenibacillus_sp.R4/SRR10074455_1.cut3.fastq

GENOME_SIZE:预估基因组大小,示例8600000意为8600000bp

THREADS:程序运行的线程数,示例为4

MIN_READ_LENGTHreads纠错(corrected reads)以及修剪(trimmed reads)的最低长度,示例3000意为3000bp

OVLP_SENSITIVE_OPTIONSOVLP_FAST_OPTIONS:存在两次reads纠错过程,第一次灵敏第二次快速;两个选项分别在reads纠错阶段检测候选重叠,参数将传递给子程序oc2pmov执行;详见oc2pmov -h

CNS_SENSITIVE_OPTIONSCNS_FAST_OPTIONS:同样地,分别对应两次纠错原始reads的参数选项,传给子程序oc2cns执行;详见oc2cns -h

TRIM_OVLP_OPTIONS:用于在reads修剪阶段检测候选重叠的选项,参数将传递给子程序oc2asmpm执行;详见oc2asmpm -h

ASM_OVLP_OPTIONS:用于在reads组装阶段检测候选重叠的选项,参数将传递给子程序oc2asmpm.sh执行;详见oc2asmpm.sh

CNS_OUTPUT_COVERAGE:以预估基因组大小为准,提取N×覆盖度的最长纠错后reads进行修剪和组装;示例为60,也就是会提取60* 8600000bp=516Mbp总长度的reads用于后续修剪和组装;该参数传递给子程序oc2elr,详见oc2elr -h

CLEANUP,是否删除中间文件,实例使用的小基因组小数据,设置为0(否)倒是无妨;大基因组还是设置为1(是)吧,不然存储空间会被严重占用。

USE_GRIDGRID_NODE:使用多个计算节点,正常情况下用不到(修改运行线程THREADS就足以了),USE_GRID设置为falseGRID_NODE0

FSA_OL_FILTER_OPTIONS:用于过滤重叠的选项,参数将传递给子程序fsa_ol_filter执行;详见fsa_ol_filter -h

FSA_ASSEMBLE_OPTIONS:用于组装修剪后reads的选项,参数将传递给子程序fsa_assemble执行;详见fsa_assemble -h

FSA_CTG_BRIDGE_OPTIONS:用于桥接初步组装后contigs的选项,参数将传递给子程序fsa_ctg_bridge执行;详见fsa_ctg_bridge -h

 

程序运行,reads校正和基因组组装

然后按顺序自动执行三步就可以了(对,比MECAT2多一步桥接过程)。

#reads 纠错
necat.pl correct necat_cfg.txt
 
#contig 组装
mecat.pl assemble necat_cfg.txt
 
#contig 桥接
necat.pl bridge necat_cfg.txt

对于reads纠错,是三代测序Nanopore数据自我纠错的环节;后面基因组组装使用校正后的reads执行三代组装;初步组装后还有个桥接环节,尽管我没搞懂桥接和组装过程的差别在哪里……关于校正和组装过程的算法原理,有兴趣的话还请自行了解。

示例的细菌基因组,使用截过的Nanopore数据,4线程运行也不到3小时。

 

工作路径中生成了结果目录“Paenibacillus_sp.R4”。

其它中间文件略过,我们关注最终的组装结果“Paenibacillus_sp.R4/6-bridge_contigs/bridged_contigs.fasta

#查看组装得到的基因组 fasta 文件
less Paenibacillus_sp.R4/6-bridge_contigs/bridged_contigs.fasta

5.png

 


顺便再整点其它的


如果您只是想知道MECAT2/NECAT怎么安装及使用,看到上面那一步就可以了,如此即得到了初步组装的基因组序列。初步组装后的工作,通常就是评估组装质量、基因组polish等一些列的过程,这些不属于这里的内容,就不继续提了。

 

从这儿开始的下文内容,算是对上文的两个组装结果作一下简单评估,若无兴趣可无需往下看。

细菌三代测序,不搞个完成图出来简直对不起三代测序啊。由于两组测序数据,PacBioNanopore,来自同一菌株,所以首先比较了一下两个组装结果。

先对组装结果做个初步统计,例如我借助QUAST

#运行后查看报告 quast/report.html
quast.py -o quast_mecat2 mecat2/Paenibacillus_sp.R4/4-fsa/contigs.fasta
quast.py -o quast_necat necat/Paenibacillus_sp.R4/6-bridge_contigs/bridged_contigs.fasta

6.png

7.png

发现本次组装,两组数据均只得到了一条contig序列。

既然没有出现第二条contig,可以推断两个结论:

1)该测序的细菌物种中没有质粒,组装所得contig序列就是这唯一的染色体;

2)染色体装配过程中没有出现断点,没有出现染色体碎片。

并且发现这里Nanopore组装的contigPacBio组装的contig要长那么一截(~47kb,对于细菌这种小基因组来说还真不小),不过并不能说明前者效果优于后者,继续往下看。

 

那么,这条contig是否就已经将细菌染色体组装完整了呢?

通过在NCBIEMBL等数据库中查找了几个同属物种基因组,发现都是环状染色体序列,可推断测序物种的染色体应该也是环状的。(其实我尝试着找过已上传的Paenibacillus sp. R4基因组,但是未找见,即使直接使用全基因组序列与NT数据库作BLAST也未匹配到99%以上相似度的基因组;好奇怪,原始测序数据都上传了,基因组怎么可能没上传呢,难道是我找的方式不对?)

8.png

然后,我借助MUMmer,比较了两个组装结果的共线性,根据共线性的结果判断是否成环。

mkdir mummer && cd mummer
ln -s ../mecat2/Paenibacillus_sp.R4/4-fsa/contigs.fasta mecat2.fasta
ln -s ../necat/Paenibacillus_sp.R4/6-bridge_contigs/bridged_contigs.fasta necat.fasta

#MECAT2 的组装 contig 为参考
nucmer --mum -p test1 mecat2.fasta necat.fasta
delta-filter -m test1.delta > test1.filter
show-coords -T -r -l test1.filter > test1.coords
 
#NECAT 的组装 contig 为参考
nucmer --mum -p test2 necat.fasta mecat2.fasta
delta-filter -m test2.delta > test2.filter
show-coords -T -r -l test2.filter > test2.coords

查看“test1.coords”,发现MECAT2组装PacBio数据所得contig,恰好成环(如上文在介绍MECAT2使用的结尾处提到,现在很多三代组装软件都能自动判别“环状序列”这种情况,不过并不是每次都能自动识别,而且识别的也不一定完全就是对的,所以组装后一定要自己核实下,这里的共线性分析方法就是一种核对方法之一)。

9.png

 

查看“test2.coords”,发现NECAT2组装Nanopore数据所得contig虽然也拼接完整了,不过不是恰到好处的那种,而是冗余了那么一段。正如前面统计的基因组长度所示,Nanopore数据所得contigPacBio组装的contig要长那么一截,原因就在这。所以不同于那些大型基因组,对于这种环状小基因组来讲,并不是说组装序列越长就是好的。

10.png

通过blast的方法,找到fasta中首尾一致的序列,在一端将多出来的截掉,就可以了。

#组装结果自身 blastn 比对,寻找首尾重叠区域
mkdir blast && cd blast
ln -s ../necat/Paenibacillus_sp.R4/6-bridge_contigs/bridged_contigs.fasta necat.fasta
makeblastdb -in necat.fasta -dbtype nucl
blastn -db necat.fasta -query necat.fasta -out blast.txt -outfmt 6 -evalue 1e-5

11.png

#去除多余的序列,保留合适的部分,例如只保留 1-8976585 非重叠区即可
#这里使用 samtools 工具获取序列
samtools faidx necat.fasta
samtools faidx necat.fasta 000000F:1-8976585 > necat_cut.fasta
 
#再和 MECAT2 的组装序列共线性分析
nucmer --mum -p test3 necat_cut.fasta ../mummer/mecat2.fasta
delta-filter -m test3.delta > test3.filter
show-coords -T -r -l test3.filter > test3.coords

这次的共线性结果显示刚刚好成环了。

12.png

 

总之到这一步,两种方式组装的细菌染色体序列都环起来了,也就是说已经将染色体拼接完整了。尽管PacBio的组装结果和Nanopore的组装结果仍然有微小差别,这是肯定避免不了的。再后面就是还需基因组polish、起点调整、正负链调整等一些列的环节,使基因组更可靠。这些在这里就不再提了,后面准备再写个关于细菌基因组完成图组装的分享统一细说。

 



http://blog.sciencenet.cn/blog-3406804-1203984.html

上一篇:列举一些次级代谢物基因簇相关的数据库
下一篇:三代测序基因组组装工具NextDenovo

1 阎影

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 17:56

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部