同源预测(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
至此同源注释已经结束
