MECAT2是三代测序数据PacBio的高效组装工具,是之前MECAT的改进版, 修复了之前的很多bug, 使用基于string graph的fsa
替换了之前的mecat2canu
。
MECAT2由4个模块组成:
-
mecat2pw
: SMART reads快速和准确地配对比对工具 -
mecat2ref
: SMART reads的参考基因组比对工具 -
mecat2cns
: 基于配对的重叠对存在噪音的reads进行纠错 -
fsa
: 基于string graph的组装工具
目前MECAT2只支持PacBio数据,对于Nanopore数据,肖老师他们开发了NECAT
安装
MECAT2软件安装非常简单,不存在上一代MECAT的安装难问题了
mkdir -p ~/opt/biosoft
cd ~/opt/biosoft
git clone https://github.com/xiaochuanle/MECAT2.git
cd MECAT2
make -j 4
export PATH=~/opt/biosoft/MECAT2/Linux-amd64/bin:$PATH
唯一要注意的是尽量要用新版的perl,经测试,Perl 5.26 是可以正常运行脚本。
最后加入到环境变量PATH
中
export PATH=~/opt/biosoft/MECAT2/Linux-amd64/bin:$PATH
软件使用
我们以 E. coli 数据集作为案例,介绍MECAT2应该如何使用。
第一步,我们要下载数据
mkdir -p ~/ecoli
cd ~/ecoli
wget http://gembox.cbcb.umd.edu/mhap/raw/ecoli_filtered.fastq.gz
gunzip ecoli_filtered.fastq.gz
注意: 目前MECAT2还不支持gz压缩文件,如果你直接用gz作为输入,会提示core dumped。
MECAT2不是根据文件名来确定输入的序列格式, 也就是如果你把 ecoli_filtered.fastq 命名为 ecoli_filtered.fasta, 组装也不会出错。
第二步, 准备模板的config文件
mecat.pl config ecoli_config_file.txt
使用vim
修改config文件
PROJECT=ecoli
RAWREADS=ecoli_filtered.fastq
GENOME_SIZE=4800000
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
CLEANUP=0
CNS_OUTPUT_COVERAGE=30
GRID_NODE=0
第三步: 原始数据纠错
mecat.pl correct ecoli_config_file.txt
第四步: 对纠错后的reads进行组装
mecat.pl assemble ecoli_config_file.txt
输出结果:
- 纠错后的reads:
ecoli/1-consensus/cns_reads.fasta
. - 最长30X纠错后用于trimming的reads:
ecoli/1-consensus/cns_final.fasta
. - trimmed reads:
ecoli/2-trim_bases/trimReads.fasta
- 组装的contigs:
ecoli/4-fsa/contigs.fasta
最终的组装结果为
Count: 1
Tatal: 4630437
Max: 4630437
Min: 4630437
N25: 4630437
L25: 1
能够完美的就装出一条序列
配置文件介绍
MECAT2的组装过程的参数都在配置文件中, 也就是之前的ecoli_config_file.txt
。我按照不同部分对这些参数进行介绍
基本参数
这些参数属于改起来不怎么纠结的参数
-
PROJECT=ecoli
: 项目名,之后的所有输出文件都在该项目名的目录下 -
RAWREADS=
: 原始数据的位置,可以是Fasta或者是Fastq, H5格式要先转成Fasta (参考 Input Format). -
GENOME_SIZE=
: 基因组大小,单位是bp,如果是120M,那么就是120000000. -
THREADS=
: 线程数 -
MIN_READ_LENGTH=
: 用于纠错和trim的reads的最低长度. 数据质量好,就长一点,比如说1000到2000 -
USE_GRID=false
: 是否有多个计算节点 -
CLEANUP=0
: 运行结束后删除MECAT2的中间文件, 大基因组的临时文件很大,所以要设置为1. -
GRID_NODE=0
, 当USE_GRID=1
时,设置用到的计算节点数,如果是单节点服务器,不需要设置。
纠错和修剪阶段
纠错阶段(correct stage)和修剪阶段(trimming stage),MECAT2调用的mecat2pw
和mecat2cns
, 与之相关的配置如下
-
CNS_OVLP_OPTIONS=""
, 在纠错阶段是检测候选重叠的参数, 会传给mecat2pw
-
CNS_OPTIONS=""
:原始reads纠错参数,会传递给mecat2cns
, -
TRIM_OVLP_OPTIONS=""
: 在trim阶段,用于检测重叠的参数,会传给v2asmpm
但是我发现v2asmpm
没有文档说明,和软件开发者讨论之后,给出的说明如下
V2asmpm -Pworkpath -Tt -Sx -Ey -B
-Pworkpath 工作目录是workpath
-Tt 用t个CPU线程
-Sx -Ey 这两个参数一起表示只计算第x到第y卷(包括第y卷), 通常工作目录下会有000001.fasta, 000002.fasta这样的分卷
-B 如果有这个选项, 就输出二进制格式的比对结果; 如果没有这个选项, 就输出文本格式的比对结果
因此用默认的-B
就行了。
组装阶段
组装阶段相关的参数如下
-
CNS_OUTPUT_COVERAGE=30
: 选择多少覆盖度的最长纠错后reads进行trim和组装。举个例子,4.8M的 E. coli 30X, 是144MB。 一般选择20~40之间就行, 会传递给mecat2elr
-
ASM_OVLP_OPTIONS=""
: 在组装阶段,用于检测重叠的参数,传给v2asmpm.sh
-
FSA_OL_FILTER_OPTIONS=""
, 过滤重叠的参数,传递给fsa_ol_filter
-
FSA_ASSEMBLE_OPTIONS=""
, 组装trimm reads的参数, 传给v2asmpm
这些参数调整起来都很复杂,只能建议看函数帮助文档了。举个例子FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1"
根据阅读fsa_ol_filter
的帮助发现,-1
表示让程序自己决定
而ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400"
, 根据和软件开发者的讨论,是不需要额外调整,用默认的即可。
不客观的个人体验测评
MECAT2是基于比对纠错的组装工具,组装速度上比同类型的Canu和Falcon都快。仅看组装结果的N50参数,MECAT2结果和Canu, Falcon结果都是差不多的。
考虑到BUSCO完整度,Canu是三者最高。不过用三代数据和二代数据进行几轮Polish后会提高到Canu相同水平。三个软件上的BUSCO(embryophyta_odb10)值如下:
Canu: C:98.8%[S:96.3%,D:2.5%],F:0.5%,M:0.7%,n:1375
MECAT2: C:93.5%[S:91.6%,D:1.9%],F:4.6%,M:1.9%,n:1375
Falcon: C:96.8%[S:95.1%,D:1.7%],F:2.0%,M:1.2%,n:1375
因此,如果你原本你是用Falcon做组装,那么可以改成MECAT2了,效果差不多,速度还更快了。对于越大的物种,在有限的资源下,更推荐用MECAT2。