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

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调用的mecat2pwmecat2cns, 与之相关的配置如下

  • 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。

参考资料

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 193,968评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,682评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,254评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,074评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 60,964评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,055评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,484评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,170评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,433评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,512评论 2 308
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,296评论 1 325
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,184评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,545评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,880评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,150评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,437评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,630评论 2 335