基因结构注释（2）：同源注释 |
您所在的位置：网站首页 › 基因组注释信息怎么填 › 基因结构注释（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至此同源注释已经结束 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |