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