什么？基因注释文件出错？
往下看！
有时候在做进化树+结构域这种组合图时，会发现偶尔有个别基因的结构域与同一基因家族的其他成员相比，会出现“缺胳膊短腿”的情况。
比如下面这种情况：
找了半天也没找到合适的例子，那就采用CJ在推文中用到的GRAS基因家族。由于我没下载菠萝基因组信息，所以就在拟南芥试试看。
花了几分钟，搞出来这个组合图。大致就是下面这个样子。但具体是不是我想要描述的问题，现在还不知道。
接下来试试看。从下图可以看出，AT5G67411.1这个基因的GRAS结构域明显比其他成员要短半截，对应到motif，应该是N端和C端都缺了一段。
那这是由什么原因造成的呢？
（当然，这个现象不是个别，是普遍）
大致就是两个原因：
要么就是注释文件出错了！导致该基因两端缺失。
要么就是它本身就是一个假基因，不应该被列入后续的分析。
更多内容请关注微信公众号：今日之森。
下面我们具体看看，到底是注释文件错了还是这个基因本身就是假基因呢？
如果注释文件错了，那我们就需要矫正它！如果本身就是假基因，那就需要舍弃它！
1、简单粗暴从gff文件获取AT5G67411.1的位置信息。
2、提取AT5G67411.1上下游5000bp内的序列。
3、提交序列至softberry网站
http://www.softberry.com/berry.phtml?topic=fgenesh&group=programs&subgroup=gfind
打开这个网站还比较费劲。
结果如下：
我们直接打开PDF报告，可以看到具体的编码情况。
但是仔细看起来还是比较头疼。
4、我们直接把刚才的mRNA序列依次提交到NCBI CDD-search
最终发现，将这段长1179bp的序列提交后，结果如下：
然后我们再把上述序列提交至Pfam。结果一致：
这个结果有什么用呢？其实还是让人有点小惊喜。
仔细看会发现，第一张图中AT5G67411.1的GRAS结构域只有约200aa。而其他成员均在350-370aa。
现在我们已经得到了AT5G67411.1基因的新序列，至少补全了GRAS结构域。所以还是有所改善的。
那么既然得到了新的序列，下一步该怎么办呢？
也许这个结果只是一个坑，如果不是坑，那我们需要矫正注释文件。
那这里如果有坑，会有几个坑呢？
坑一：
刚才我们选取了基因本身的序列及上下游各5000bp，那么刚才预测得到的长1179 bp的基因是在AT5G67411.1基础上向5‘端延伸了几百bp呢？还是只是另外一个基因呢？
坑二：
这里假设1179bp的新序列是在AT5G67411.1基础上向5‘端延伸了几百bp后补全了，那么有没有和其他基因有overlap呢?
带着这些问题，那我们继续。
先不着急对注释文件进行注释。我们先来看看刚才选取的序列有没有overlap。其实这一步在刚提取序列的时候就该考虑，既然刚才没有考虑。那现在就来看看。
果不其然，前后都有overlap。
那怎么办呢？那就不要在上下游取5000了，缩短一点再试试。
那就取刚刚不overlap的区间。AT5G67411.1 26897680 26903865
结果如下：
看到这个结果就放心了。说明刚才我们选取的长度为1179的序列是没问题的。
好，到这里，就算对AT5G67411.1这个基因的GRAS结构域进行了完善。
我们再来比对一下完善之前和完善之后的序列。
可以看到，我们对5'和3'都进行了补充。
到这里，我们的分析就告一段落了。基本可以确定这个基因的注释信息有些问题，那我们发现之后就需要对原始的注释文件进行矫正。
既然写到这里了，那就一气呵成。
下面开始对注释文件进行矫正。
最后只需要将矫正后的注释文件粘贴进原始gff文件即可。
题外话：
说到gff文件，大家都感到陌生。其实到底是什么，我们打开看看就知道了。
比如：
很清楚。应该不需要解释吧。
