MECAT2:三代高效组装工具

MECAT2: 三代高效组装工具

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的安装难问题了

1
2
3
4
5
mkdir -p ~/opt/biosoft
cd ~/opt/biosoft
git clone https://github.com/xiaochuanle/MECAT2.git
cd MECAT2
make -j 4

唯一要注意的是尽量要用新版的perl,经测试,Perl 5.26 是可以正常运行脚本。

最后加入到环境变量PATH

1
export PATH=~/opt/biosoft/MECAT2/Linux-amd64/bin:$PATH

软件使用

我们以 E. coli 数据集作为案例,介绍MECAT2应该如何使用。

第一步,我们要下载数据

1
2
3
4
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文件

1
mecat.pl config ecoli_config_file.txt

使用vim修改config文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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

第三步: 原始数据纠错

1
mecat.pl correct ecoli_config_file.txt

第四步: 对纠错后的reads进行组装

1
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

最终的组装结果为

1
2
3
4
5
6
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调用的mecat2pwmecat2cns, 与之相关的配置如下

  • CNS_OVLP_OPTIONS="", 在纠错阶段是检测候选重叠的参数, 会传给mecat2pw
  • CNS_OPTIONS="":原始reads纠错参数,会传递给mecat2cns,
  • TRIM_OVLP_OPTIONS="": 在trim阶段,用于检测重叠的参数,会传给v2asmpm

但是我发现v2asmpm没有文档说明,和软件开发者讨论之后,给出的说明如下

1
2
3
4
5
6
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)值如下:

1
2
3
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。

参考资料