注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

云之南

风声,雨声,读书声,声声入耳;家事,国事,天下事,事事关心

 
 
 

日志

 
 
关于我

专业背景:计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

网易考拉推荐

基因预测流程  

2013-08-28 10:40:10|  分类: 生信分析软件 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |


http://blog.sina.com.cn/s/blog_5542c24a0100uadb.html

在读下面的文字之前,先推荐一篇NRG上面的文章,去读一下,那么下面的文字不用看就可以了。

A beginner’s guide to eukaryotic genome annotation

在基因注释之前,先要mask基因组,这个mask其实还是有技巧的,因为在植物里面很多repeat也存在于基因内部,甚至在exon里面也不少转座子 的残留,简单粗暴的把所有repeat给mask掉了会丢失不少信息,不过我没有深究,我只是简单在abinitio阶段使用mask过的基因组,而在同 源搜索和EST比对阶段则不对基因组mask.
1.abinitio工具-augustus,glimmerHMM,genemark.HMM
我觉得没必要用太多的denovo工具,就挑了几个有禾本科matrix的工具,其实可能一个augustus就够了。
如果一个基因组没有近缘种已测完,并且很少EST序列,那么abinitio应该还是比较有意义的。
2.同源预测工具
2.1 exonerate,用它,主要是速度快,而且结果还算不错,比其他一些简单的工具要强。我用了所有禾本科的蛋白来做evidence,运气好的时候,我抢到接近500个CPU时,36小时可以跑完。
2.2 genewise,这个其实很麻烦,本来我不想用的,但是在不用它之前,我注释的结果很不理想,比ensembl的结果差了一截,而ensembl注释的核心工具就是genewise,这工具也是他们相关的小组开发的。
要用它,得自己写不少代码,我懒,所以不想写,不过最后避免不了还是得写。
首先,wise很慢,无敌的慢,你要用wise直接对基因组进行注释绝对是疯了。
要解决这个问题,就得先用蛋白搜基因组,找到可能的目标区域,把这个区域截取出来,再对目标区域进行wise。
2.3 目标区域搜索及HSP linker
搜索目标区域看起来简单,其实还是一件麻烦事,你用tblastn搜,但是blast是个local aligment,一个蛋白会有无数的hsp,如果每个hsp你都去延伸,再截取基因组,一样没效率没意义。
所以,我需要一个hsp linker之类的工具,就是根据hsp之间的关系把那些在基因组上顺序排列并且相隔不太远的hsp link在一起,之后再去做延伸。
hsp link是有一些工具可以用,比如你可以用blat,blat的结果可以用它的一个小工具link起来,也可以用bioperl的模块来寻找hsp clusters,或者wublast,使用一个参数-link也可以达到目标。
但是,我仔细检查了一下结果,发现这些工具都不行,结果很垃圾,一个1kb的蛋白最后匹配个100kb的基因组区域,甚至1Mb,根本没法用。
也萌生过自己写的想法,但是要达到比较好的效果,我的能力不够,根据我的观察,就算计算机出身的,如果对算法不够精通,开发的工具很多都没什么用-可以得到结果,但是结果很不好。
最后,我找到了一个工具,叫genblastA,非常符合我的需求,我觉得是现阶段最好的hsp linker了,当然这个工具主要是针对基因注释的。
这个工具很简单,用它搜完基因组,把结果提取出来,整理成tab,就可以进行基因组目标区域截取了。
目标区域截取我是这样的,根据hsp在protein和基因组上start,end,向前后各自延伸3*(protein start-1),3*(protein length-protein end),在此基础上再延伸1kb,因为我的基因组还是有近缘种做reference的,所以1kb应该够了。
目标区域拿到之后就可以做wise了。
wise提供了gff格式的输出,还不错。

3.EST工具
EST工具无疑问,最好的工具就是PASA,虽然慢,但是质量很高。
PASA的文档很详细,参见http://pasa.sourceforge.net
按部就班就可以了。

4.整合工具
当你用各类工具完成了注释后,一个问题是,每一个基因组区域,你都会获得大量redundant的基因结构注释,到底哪一个注释才是最可靠的?所以,我们需要一个整合工具。有些基因组注释过程比较简单粗暴,就是把所有的gene model都翻译了,然后搜已知蛋白数据库,留下最长的,同源性最好的gene model,或者用更复杂一点的筛选组合,例如:
The best predictions for each locus are selected using multiple positive factors: C-score, peptide
homology coverage, transcript splicing and coverage scores by EST assembly alignments, and
one negative factor: overlap with repeats.
另一条路线就是统计打分,先手工注释一些基因,然后把所有的自动注释结果跟手工注释比较,给各个工具打个分,最后用这个打分矩阵扩展到所有的基因上。常用的有glean,jigsaw。
glean,这个工具应该还是不错的,它不需要training,之前用于果蝇12个物种基因组的注释,但是我不会用,文档就半页。。。
jigsaw需要training,我没时间做手工注释。
我用evidence modeler的原因比较直接,它是pasa的开发者做的另一个工具,文档详细,使用简单,而且它自带很多格式转换代码,最重要的是,这个开发者之前也是 做禾本科的,所以虽然这个工具也需要training,但是它也可以人工打分,虽然人工给的不是最好,但也是次好,至少它是根据禾本科的注释推出来的一个 分数。
http://evidencemodeler.sourceforge.net/

5.流程
abinitio,exonerate,genblasta,pasa可以同时跑,各不相干,
genblastA跑完后再做wise.
所有工具跑出来的结果全部转换成gff3格式
用evidence modeler提取最好的gene model
这个gene model再丢个pasa,进行UTR延伸,可变剪接注释。

全部流程跑完4-10天,得看运气,超算紧张时候,100个CPU我都抢不到,宽松时候500多个CPU都有。

6.并行化
以上的所有工具都是单线程的工具,不过大部分工具的输入都是序列文件,只需要分割文件,每个小文件单独一个进程就可以了。wise麻烦一点,因为有几十万 的wise比对,如果每个都做一个进程提交,效率非常低,CY同学给我提供了一个他写的脚本,我感觉对我用处不大,最后用了个折衷的办法,把每一个 wise的命令行都写到一个文本里,然后分割文本,每个小文本作为一个进程直接提交,就相当于一个进程循环执行500条左右的wise运算了。
lsf这个集群管理系统比较诡异,如果进程采用循环提交,甚至使用job array的方式提交,都会出现进程丢失,我搞不清什么原因。
我现在只发现一个提交模式可以避免进程丢失,就是写一个提交进程的规范脚本,然后提交脚本,例如:
#!/bin/bash
# enable your environment, which will use .bashrc configuration in your home directory
#BSUB -L /bin/bash
# the name of your job showing on the queue system
#BSUB -J Augustus[1-200] 这里是你分割的文件的总数
# the queue that you will use, the example here use the queue called "normal"
# please use bqueus command to check the available queues
#BSUB -q Q_Serial
# the system output and error message output, %J will show as your jobID
#BSUB -o "./abinitio/aug_tmp/%J.%I.out" 屏幕输出,程序运行信息
#BSUB -e "./abinitio/aug_tmp/%J.%I.err" 屏幕输出,错误信息
#the CPU number that you will collect (Attention: each node has 2 CPU)
#BSUB -n 1
#enter your working directory, change to your own dir
./abinitio
#Finally, Start the program:命令行
augustus --gff3=on --species=maize --outfile=./abinitio/aug_out/gla.masked.f2.$LSB_JOBINDEX.aug_gff3 ./genomic_chunk/gla.masked.f2.$LSB_JOBINDEX.seq

7.最后筛选
注释出来的gene model还是会有很多repeat,还有很多是假基因,还有一些ORF太短,可能是错误注释,这一步筛选比较简单,先砍掉编码蛋白<50AA的gene model,然后所有蛋白或者gene model的序列对repeat数据库blast,把那些repeat序列占据超过20% CDS区域的基因扔掉,再把中间有提前终止密码子的扔掉,基本上就完工了。

我的流程只是针对大规模基因组注释的,特别是刚测出的基因组,要精细注释,这个流程还是很粗糙的,影响最大的步骤可能是最后一步的整合工具,没有有效的training可能结果并不够好。
  评论这张
 
阅读(1241)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016