「三代组装」Pacbio组装后如何用自身数据进行polish(更新版)

之前那我由于需要对PacBio的组装结果进行polish,于是写了「三代组装」Pacbio组装后如何用自身数据进行polish。最近发现自己又有了需求,于是重新回顾了我之前写的这篇文章,但是在实践的时候,却发现新版本pb-assembly(0.0.8)没有之前的pbalignvariantCaller这两个命令,也就是意味着我之前的文章随着版本更新而失效了。

此时我就面临着两个选择,一个用conda安装之前的版本,逃避可耻但有用。另一个就是探索新版本下的代码,勇敢面对惨烈的人生。由于时间不是特别的紧张,所以我选择了第二条路。

通过我一波检索,我发现第二条路是一条没有人走过的路。https://github.com/PacificBiosciences/pb-assembly只讲了falcon组装后,使用falcon-unzip进行polish,而没有介绍使用其他软件(如Canu, Mecat2) 应该怎么做。https://github.com/PacificBiosciences/GenomicConsensus提到的quiver命令并不存在于我们pb-assembly中。

我虽然把pb-assembly下的命令都看了个遍,但是也没有想到哪些命令可以组合在一起用来polish。但是我突然想到,既然falcon-unzip里面会有polish,那么的它运行日志里面一定会有它如何进行polish蛛丝马迹。经过一波的分析和追踪,终于被我找到了相关的代码。修改之后如下

set -e
REF=$1
BAM=$2
THREADS=80

# align
pbmm2 align --sort  $REF $BAM aln.bam
pbindex aln.bam
# arrow
gcpp --algorithm=arrow -x 5 -X 120 -q 0 -j $THREADS -r $REF  aln.bam -o cns.fa,cns.fq,cns.vcf

新的脚本不再用BLASR进行比对,而是拥抱了minimap2,提高了比对速度,并且避免了之前BLASR可能会出现的内存泄露问题。使用gcpp替代了之前variantCaller,但是参数几乎没有变动,你可以认为就是改了一个名字而已。

另外这是一个简化的脚本,如果是多个bam文件,则需要分别比对然后合并。原来的falcon流程会将参考基因组按照染色体拆分然后并行处理,但是我们的脚本设置了更高的线程数,最后结果速度也差不多。

©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

友情链接更多精彩内容