如何在5分钟内拿到细胞器基因组的进化树
前记-HomoBlock
雷神之锤的我,之前设计了一个HomBlocks流程,可以很方便地帮助大家构建细胞器基因组的比对序列,再用来建树。大家如果不明所以的话,可以搜索查看基迪奥论坛上我的其他的帖子,就知道这个流程具体怎么运行,最终用来干什么了。
HomBlocks也发表在了Genomics上,地址:http://www.sciencedirect.com/sci ... i/S088875431730068X
当时发帖的时候,号称一天之内拿到细胞器基因组树。自然,在序列数目不多的时候,HomBlocks的运行还是比较快的,一天之内拿到比对好的序列,用MEGA跑一下NJ树还是没有问题的。但是随着序列数目的增加、物种亲缘性的疏远、序列长度差异越大或者清一色使用叶绿体基因组,会导致软件运行的时间还是比较久。比如使用70个叶绿体序列,来跑HomBlocks可能会花上几天时间。
因为Mauve的比对只能单线程运行,所以速度上没有办法再提高。
虽然HomBlocks的依然能跑出结果,但是让人等待的过程也实在是煎熬。什么时候能跑出来?跑到哪里了?程序是不是中断了?
再虽然,HomBlocks相比传统手动方法(手动挑取基因序列进行比对再合并序列)在效率上和简便性上,已经算得上很不错啦。起码你只需要等待,无须手动,也不怕出来的序列有问题。
但是!就真没有更快的方法了吗?
升级-BLAST2OGMSA
雷神の锤告诉你,还真有!我给你搞出来了!现在的方法号称:
“五!分!钟!内!拿!到!细!胞!器!基!因!组!的!进!化!树!”。
我们这个新的方法(也算不上什么流程),叫做BLAST2OGMSA。
这个新方法的文章题目我也想好了——BLAST2OGMSA: an convenient way to construct whole organelle genomes alignment to facilitate phylogeny analysis.
大家看到这个题目大概也了解了,为什么我们能在5分钟之内拿到这么个东西了。
比对的话,什么最快?mafft?clustalw?muscle? No, No, No, No。是BLAST!
不然NT数据库那么大,你怎么能够那么快速地在网页上检索序列呢?
那我们需要自行构建序列的库,然后两两比对,再解析blast结果吗?
不需要!
正常的思路的确是这样,利用blast的快速比对(虽然在比对精度上不如其他正统序列比对软件),来解析blast结果。但是,还是挺麻烦的。首先要把blast结果保存成xml格式,然后再提取序列,序列再定位,再提取.............然后你会发现,对位还是有问题,还不如直接用mafft来比呢。
但BLAST2OGMSA直接避开了这么个过程,直接使用ncbi网页版BLAST工具就行了。
BLAST2OGMSA就是一个脚本,只需要你提供两个文件就行。
那这两个文件从哪来?
操作-超快
那顺便说一下BLAST2OGMSA的使用方法:
首先,将你需要构建系统发生树的所有序列,合并在一个文件里。
打开ncbi的BLAST工具, 点开Align two or more sequences选项。
image
将你准备好的序列,随便选一个,放入第一栏;将所有序列选中,放入第二栏,然后比对啦。
Blast.png
在蹦出来的结果页面,点开MSA View项。
MSA view.png
在蹦出来的MSA View界面,将blast比对出来的结果download下来(这个是最初的结果,但不能用,我们需要从中提取有效信息),举个例子我们保存成X.aln文件。
x.aln.png
在原来的结果界面,在下部分的结果部分,我们选中所有序列,选择download,文件名就默认为seqdump.txt就好。
seqdump.txt.png
seqdump2.txt.png
我们有了X.aln和seqdump.txt两个文件,然后运行BLAST2OGMSA:
perl BLAST2OGMSA.pl -method=Gblocks X.aln seqdump.txt output.fasta
然后产生的output.fasta就是我们可以用来建树的序列啦!然后简单跑一下NJ:
nj.png
结果-完美
还记得我第一篇帖子中的36个啮齿类动物线粒体的系统发生树吗?结构一致哦,节点的支持度都OK!这个过程花了多长时间?5分钟都没有吧!
BLAST2OGMSA也跟HomBlocks提供了4中序列修剪方法:
Gblock
trimAl
BMGE
Noisy
所以大家用的时候,可以按照自己的需要进行选择。
BLAST2OGMSA的用法也实在很简单,没有别的参数,只需要按顺序把文件列好就行了,可以参考例子中的命令,如果不知道怎么用,那么直接perl BLAST2OGMSA.pl查看一下就行了。
软件下载与使用
BLAST2OGMSA的地址:https://github.com/fenghen360/BLAST2OGMSA
下载解压之后,需要将bin目录下的所有文件可执行化:chmod 755 * 就行。
刚有空写出来,readme都没来得及详细写,大家忽略readme就行,我有空再好好写写。
引用
大家引用的话,就先这么写就好:
the genome-wide alignment of XXX genomes was constructed by BLAST2OGMSA ([https://github.com/fenghen360/BLAST2OGMSA) (https://github.com/fenghen360/BLAST2OGMSA), resulting in XXXXX characters of each species, which including all PCGs and rRNA genes.
后记-联系作者
有三个方式可以找到作者
bioinformatics{*}中国(276151571) - 还剩几个空位
bioinformatics{*}中国大分舵(744366744)
Omicshare论坛