基因结构注释（2）：同源注释

同源预测(homology prediction)利用近缘物种已知基因进行序列比对，找到同源序列。然后在同源序列的基础上，根据基因信号如剪切信号、基因起始和终止密码子对基因结构进行预测.

在同源预测上，目前看到的大部分基因组文章都是基于TBLASTN + GeneWise,但是目前GeneWise已经不在维护，在本次推文中将使用GenomeThreader进行同源注释。

确定同源物种蛋白序列

首先选择同源物种，在本次分析中，我使用Saccharomyces_cerevisiae、Laccaria_bicolor、Amanita_thiersii、Pleurotus_pulmonarius、Pterula_gracilis进行同源注释

cat Saccharomyces_cerevisiae.fa\       
    Laccaria_bicolor.fa \
    Amanita_thiersii \
    Pleurotus_pulmonarius \
    Pterula_gracilis >all.pep.fa

然后使用TBLASTN确定匹配到基因组的蛋白序列

makeblastdb -in pudorinus.fa -parse_seqids -dbtype nucl -out index/pu&
## tblastn比对
nohup tblastn -query all.pep.fa -out pu.blast -db index/pu -outfmt 6 -evalue 1e-10 -num_threads 8 -qcov_hsp_perc 50.0 -num_alignments 5 &
## 提取匹配到基因组的蛋白序列
awk '{print $1}' pu.blast >pu.list
sort pu.list| unique >pi.ho.list
seqkit seq  all.pep.fa   -w   0   >  all.fa
vi sh.sh
   cat list.ru | while read line
   do
   grep "$line" -A 1 all.fa >$line.1
   done
   cat *.1 > pudorinus.homo.fa
   rm -rf *.1

使用gth进行同源注释

gth -genomic pudorinus.fa -protein pudorinus.homo.fa -intermediate -gff3out > pudorinus.gff

删除一些无用的信息可以看到基本已经注释完成，但是exon之类的没有ID，在后续识别会有问题

grep -v "^#" pudorinus.gth.gff | grep -v "prime_cis_splice_site" | awk -F ";" '{print$1}'>pudorinus.homo.gff
less pudorinus.homo.gff

使用脚本处理下(谢谢课题组师姐帮忙写的脚本(#^.#))

#!/lustre/home/guohan_lab/local/python-3.6/bin/python3
#liyumei
#./change-name_lym.py proteinprediction.gff proteinprediction.gff proteinprediction_r.gff

import sys
import re
fout = open(sys.argv[3],'w')

ref_dict={}
with open(sys.argv[1]) as gff:
        for line in gff:
                line_s = line.strip().split('\t')
                if 'gene' == line_s[2]:
                        ref_g = re.split('=',line_s[8])
                        ref_gene = ref_g[1]
                        ref_dict[ref_gene]=[]
                else:
                        pos = line_s[3]
                        ref_dict[ref_gene].append(pos)

with open(sys.argv[2]) as gff_r:
        for eachline in gff_r:
                i = eachline.strip().split('\t')
                info = re.split('=',i[8])
                name = info[1]
                ref_set = []
                for n in range(0,8):
                        ref_set.append(i[n])
                ref_list = '\t'.join(ref_set)
                ref_set1 = ['CDS' if x == 'exon' else x for x in ref_set]
                ref_list1 = '\t'.join(ref_set1)
                if 'gene' == i[2]:
                        fout.write(eachline)
                elif 'exon' == i[2]:
                        if name in ref_dict:
                                vs = ref_dict[name]
                                len_vs = len(vs)
                                for a in range(0,len_vs):
                                        if vs[a] == i[3]:
                                                b = vs.index(vs[a])
                                                fout.write('%s\tID=%s.t%d;Parent=%s\n'%(ref_list,name,b+1,name))
                                                fout.write('%s\tID=%s.c%d;Parent=%s\n'%(ref_list1,name,b+1,name))
fout.close()

python3 ../annotion/change-gff_lym.py pudorinus.homo.gff pudorinus.homo.gff prediction.pudorinus.gff

至此同源注释已经结束