三代组装软件Canu使用
Canu简介
Canu是Celera的继任者,能用于组装PacBio和Nanopore两家公司得到的测序结果。
Canu分为三个步骤,纠错,修整和组装,每一步都差不多是如下几个步骤:
- 加载read到read数据库,gkpStore
- 对k-mer进行技术,用于计算序列间的overlap
- 计算overlap
- 加载overlap到overlap数据库,OvlStore
- 根据read和overlap完成特定分析目标
- read纠错时会从overlap中挑选一致性序列替换原始的噪声read
- read修整时会使用overlap确定read哪些区域是高质量区域,哪些区域质量较低需要修整。最后保留单个最高质量的序列块
- 序列组装时根据一致的overlap对序列进行编排(layout), 最后得到contig。
这三步可以分开运行,既可以用Canu纠错后结果作为其他组装软件的输入,也可以将其他软件的纠错结果作为Canu的输入,因此下面分别运行这三步,并介绍重要的参数。
几个全局参数:genomeSize设置预估的基因组大小,这用于让Canu估计测序深度; maxThreads设置运行的最大线程数;rawErrorRate用来设置两个未纠错read之间最大期望差异碱基数;correctedErrorRate则是设置纠错后read之间最大期望差异碱基数,这个参数需要在 组装 时多次调整;minReadLength表示只使用大于阈值的序列,minOverlapLength表示Overlap的最小长度。提高minReadLength可以提高运行速度,增加minOverlapLength可以降低假阳性的overlap。
组装实战
数据下载
数据来自于发表在 Nature Communication 的一篇文章 "High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell"。 这篇文章提供了 Arabidopsis thaliana KBS-Mac-74 的30X短片段文库二代测序、PacBio和Nanopore的三代测序以及Bionano测序数据, 由于拟南芥的基因组被认为是植物中的金标准,因此文章提供的数据适合非常适合用于练习。根据文章提供的项目编号"PRJEB21270", 在European Nucleotide Archive上找到下载地址。
![](/assets/blank.gif)
## PacBio Sequal
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam
## MinION
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116595/fastq/ont.fq.gz
# Illuminia MiSeq
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_1.fq.gz
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_2.fq.gz
下载的PacBio数据以BAM格式存储,可以通过安装PacBio的smrtlink工具套装,使用其中的bam2fasta工具进行转换.
# build index for convert
~/opt/biosoft/smrtlink/smrtcmds/bin/pbindex pb.bam &
# convert bam to fasta
~/opt/biosoft/smrtlink/smrtcmds/bin/bam2fasta -o pb pb.bam &
其实samtools fasta
也可以将bam转成fasta文件,并且不影响之后的组装。
PacBio的smrtlink工具套装大小为1.4G,不但下载速度慢,安装也要手动确认各种我不清楚的选项, 唯一好处就是工具很全。
运行Canu
第一步:纠错。三代测序本身错误率高,使得原始数据充满了噪音。这一步就是通过序列之间的相互比较纠错得到高可信度的碱基。主要调整两个参数
- corOutCoverage: 用于控制多少数据用于纠错。比如说拟南芥是120M基因组,100X测序后得到了12G数据,如果只打算使用最长的6G数据进行纠错,那么参数就要设置为50(120m x 50)。设置一个大于测序深度的数值,例如120,表示使用所有数据。
canu -correct \-p ath -d pb_ath \Threads=10 gnuplotTested=true\genomeSize=120m minReadLength=2000 minOverlapLength=500\corOutCoverage=120 corMinCoverage=2 \-pacbio-raw pb.fasta.gz
可以将上述命令保存到shell脚本中进行运行, nohup bash run_canu.sh 2> correct.log &
.
注: 1.7版本里会没有默认没有安装gnuplot出错,因此gnuplotTested=true 可以跳过检查。
第二步:修整。这一步的目的是为了获取更高质量的序列,移除可疑区域(如残留的SMRTbell接头).
canu -trim \-p ath -d pb_athmaxThreads=20 gnuplotTested=true\genomeSize=120m minReadLength=2000 minOverlapLength=500\-pacbio-corrected ath/pb_ath.correctedReads.fasta.gz
第三步: 组装。在前两步获得高质量的序列后,就可以正式进行组装. 这一步主要调整的就是纠错后的序列的错误率, correctedErrorRate,它会影响utgOvlErrorRate。这一步可以尝试多个参数,因为速度比较块。
# error rate 0.035
canu -assemble \-p ath -d ath-erate-0.035 \maxThreads=20 gnuplotTested=true \genomeSize=120m\correctedErrorRate=0.035 \-pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz
# error rate 0.050
canu -assemble \-p ath -d ath-erate-0.050 \maxThreads=20 gnuplotTested=true \genomeSize=120m\correctedErrorRate=0.050 \-pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz
最后输出文件下的ath.contigs.fasta
就是结果文件。
一些宝贵的建议
Nanopore组装
对于Nanopore数据,使用Canu组装并不是一个非常好的选择,我曾经以一个600多M的物种100X数据进行组装,在120线程,花了整整一个多月的时间,尽管它的组装效果真的是很好。
- rawErrorRate: 从0.144 调整到 0.12 或者更低,速度会提高5到10倍
- readSamplingCoverage, readSamplingBias: 可以抽取部分数据进行纠错,而不是全部数据
- corOutCoverage=30: 默认的30X或者40X的组装效果其实不错
-fast
: 对于1Gb以下的物种,可以加上该参数,会明显提高速度
高杂合物种组装
对于高杂合物种的组装,Canu建议是用 batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50
参数尽量分出两套单倍型,然后对基因组去冗余。
batOptions
表示传递后续的参数给组装软件bogart
, -dg 3 -db3
降低自动确定阈值时的错误率离差(deviation),从而更好的分开单倍型。-dr 1 -ca 500 -cp 50
会影响错误组装的拆分,对于一个模棱两可的contig,如果至少另一条可选路径的overlap长度至少时500bp,或者说另一条可选路径时在长度上和当前最佳路径存在50%的差异,那么就将contig进行拆分。
关于杂合物种组装的讨论,参考https://github.com/marbl/canu/issues/201#issuecomment-233750764
购买SSD避免服务器IO瓶颈
如果你的服务器线程数很多,你在普通的机械硬盘上运行组装,而且你的系统还是CentOS,那么你需要调整一个参数,避免其中一步的IO严重影响服务器性能。
Canu通过两个策略进行并行,bucketizing (‘ovb’ 任务) 和 sorting (‘ovs’ 任务)。 bucketizing会从1-overlap读取输出的overlap,将他们复制一份作为中间文件。sorting一步将这些文件加载到内存中进行排序然后写出到硬盘上。 如果你的overlap输出特别多,那么该步骤将会瞬间挤爆的你的IO.
为了避免悲剧发生,请增加如下参数: ovsMemory=16G ovbConcurrency=15 ovsConcurrency=15
, 也就是降低这两步同时投递的任务数,缓解IO压力。
如何用不同服务器处理同一个任务
overlap这一步时间久,如果并没有服务器集群,而是有多台服务器,可以参考如下方法进行数据并行处理(必须要安装相同的Canu版本)
假如Canu任务中的prefix参数为coc, 用scp进行数据的传递
scp -r coc.seqStore wangjw@10.10.87.132:/data1/wjw/Coc/Canu
scp -r unitigging/1-overlapper/overlap.sh wangjw@10.10.87.132:/data1/wjw/Coc/Canu
scp -r unitigging/0-mercounts/ wangjw@10.10.87.132:/data1/wjw/Coc/Canu
到目标服务器的canu目录下
mkdir unitigging
mkdir -p unitigging/1-overlapper
mv 0-mercounts/ unitigging/
mv overlap.sh unitigging/1-overlapper/
修改overlap.sh
里的bin路径,指定当前服务器的canu路径
之后运行任务即可
cd unitigging/1-overlapper/
./overlap.sh 77 &> 77.log
运行完任务之后,将目前服务器unitigging/1-overlapper/001
目录下的000077.oc 000077.ovb 000077.stats
文件传送回原来的服务器unitigging/1-overlapper/001
下即可。
参考资料
- 官方文档: https://canu.readthedocs.io/en/latest/index.html
三代组装软件Canu使用相关推荐
- 三代组装软件canu学习笔记
三代组装软件canu学习笔记 (2017-08-07 14:17:43) 转载▼ 分类: 三代 1:这个组装软件起源于PBcR包含在Celera Assembler中(http://wgs-ass ...
- Canu|三代组装软件
Canu Canu是基于OLC算法,具有长reads的自纠错和组装功能,是应用最为广泛的三代组装软件(速度相对较慢). Canu相关使用方法来啦!!! 一.Canu的安装 方法一:conda安装 co ...
- 三代数据组装软件canu
三代拼接工具canu.canu这款软件来自于经典的 Celera Assembler.celera最早用于人类全基因组计划,后来随着二代测序数据的普及,这种基于overlap的方法逐渐被基于kmer的 ...
- nextpolish安装_希望组自主三代组装软件NextDenovo最新版本全球学术开源!
图1. NextDenovo V2.0-beta.1上线Github 测序中国2019年10月18日消息,希望组面向全球释放三代测序数据高效纠错.组装软件NextDenovo最新版本V2.0-beta ...
- 三代组装软件miniasm笔记
我们用来练手的文章发表在 Nature Communication ,"High contiguity Arabidopsis thaliana genome assembly with a ...
- canu三代组装软件安装
官方推荐流程: 1. 下载canu-2.2.tar.xz(官方推荐不下source code) 1.1 curl命令下载 curl -L https://github.com/marbl/canu/r ...
- NBT:宏基因组二、三代混合组装软件OPERA-MS
文章目录 宏基因组二.三代测序混合组装软件OPERA-MS 热心肠日报 摘要 主要结果 图1. OPERA-MS工作流程图 图2. 宏基因组数据混合组装基因组评测 图3. 组装虚拟肠道微生物组 图4. ...
- 不同测序数据应用与基因组组装|Canu
介绍 Canu专门组装PacBio或Oxford Nanopore序列.Canu分为三个阶段:校正.修整和装配.校正阶段将提高读取中基数的准确性.微调阶段将微调显示为高质量序列的部分的读取,删除可疑区 ...
- Nature方法 | 三代长读长宏基因组组装软件metaFlye
简介 标题:metaFlye:基于重复图的可拓展长序列宏基因组序列组装 metaFlye: scalable long-read me ...
最新文章
- 线刷一加5t android 9,一加5/5T 氢OS 9.0稳定四版 侧边工具 通知特效 Magisk 极速流畅 简约实用-刷机之家...
- ASP.NET CORE之上传文件夹
- python中threading模块详解及常用方法_Python常用模块功能简介(二)threading
- 印刷点阵字体_印刷术如何确定可读性:衬线与无衬线,以及如何组合字体。
- 运营商取消话费余额有效期后改收闲置费
- 作家如何利用Git更好地完成工作
- 【Flink】Flink on RocksDB 参数调优指南
- 如何让fragment每次都重新加载_带上税控在哪都能开票,如何重新安装开票软件?...
- react 使用rewired_create-react-app 通过 react-app-rewired 添加 webpack 的 alias
- 【2018-2019 ACM-ICPC Pacific Northwest Regional Contest (Div. 1) D】Count The Bits【数位DP】
- 操作系统 chapter1 操作系统概述
- JTAG接口针脚定义及含义
- 四川自考计算机信息管理专业,计算机信息管理专业2019年10月四川自考科目及考试时间[本科]...
- canvas教程15-变形
- MG513P30 12V直流减速电机编码器电线与杜邦线焊接教程
- DW个人网站设计 练习作业
- 亚马逊入驻商用户画像分析
- Sequencer和Sequence
- 手把手的教你安装PyCharm --Pycharm安装详细教程(一)(非常详细,非常....)
- python 拆分excel工作表_用python编写的excel拆分小工具