0_based numbering: the initial element of a sequence is assigned the index 0; 1_based numbering: the initial element of a sequence is assigned the index 1; BAM/BED data are 0-based and SAM data are 1-based. An example shows below: The screenshot came from a sam file. The sequence (10.4|823|29663|59.364374575|+|61.9047619048|73.3|1) was mapped to mm10, starting from 3259346 in chr10. As sam file is 1_base, which means first locus 3259346 in chr10 is the first nucleotide “A” of the sequence. However, after I converted this to bed file using bamToBed, the iniial locus of the region converts from 3259346 to 3259345. And as Bed file is 0-based, so the second locus in chr10 is the first nucleotide “A”. To double check, I extracted the the sequence of this region (chr10:3259345-3259372) in bed file as shown below. Comparing with the sequence (AAGGGGCTGGACTTGCATGCCATGGAT) in the sam file, we can know that the sequence indeed start from the second locus.
http://blog.sina.com.cn/s/blog_5542c24a0100uadb.html 基因组注释-基因注释流程 (2011-09-16 06:44:11) 转载 ▼ 标签： 基因注释 杂谈 分类： NGS 在基因注释之前，先要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 这里是你分割的文件的总数 # 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可能结果并不够好。
标签: 基因组注释
