1.在uniprot下載fungi蛋白質(zhì)序列
在uniport數(shù)據(jù)庫輸入taxonomy:fungi AND reviewed:yes查找蛋白質(zhì)相關(guān)數(shù)據(jù)调俘，總共找到33905條序列问芬，下載并解壓严卖。
2.創(chuàng)建本地數(shù)據(jù)庫
使用makeblastdb程序楣黍，對上述FASTA格式的基因組序列進行處理胖秒，建立本地BLAST數(shù)據(jù)庫馆匿。
/opt/BioBuilds/bin/makeblastdb -in GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -input_type fasta -title yjm1338 -dbtype nucl -out yjm1338
3.全基因組的同源基因搜索
使用tblastn程序滞欠，把已知蛋白質(zhì)序列和上述建立的本地BLAST數(shù)據(jù)庫進行比對忙菠。
nohup tblastn -query uniprot-taxonomy%3Afungi+reviewed%3Ayes.fasta -db yjm1338 -out ./blastx_results.outfmt6 -evalue 1e-5 -outfmt 6 -max_target_seqs 1 -num_threads 10 &
nohup tblastn -query uniprot-taxonomy%3Afungi+reviewed%3Ayes.fasta -db yjm1338 -out ./blastx_results.outfmt7 -evalue 1e-5 -outfmt 7 -max_target_seqs 1 -num_threads 10 &
使用blast92gff3.pl和blast2gff.py程序何鸡，分別把結(jié)果轉(zhuǎn)成GFF3格式
Python轉(zhuǎn)化
python2 blast2gff.py -b blastx_results.outfmt6 -t gene -p sp -F >results_py.gff3
去除低質(zhì)量的序列，共找到36367條序列
Perl轉(zhuǎn)化
perl blast92gff3.pl -geneType gene -exonType exon blastx_results.outfmt6 > results_pl.gff3
Python轉(zhuǎn)化的結(jié)果不包含exon,當去除低質(zhì)量的序列后牛欢，序列條數(shù)少于perl轉(zhuǎn)化的結(jié)果骡男。
- 同源搜索結(jié)果的評估
使用gffcompare工具把Perl轉(zhuǎn)化后的結(jié)果與原始GFF格式數(shù)據(jù)進行比較，查看結(jié)果傍睹，并分析它們之間的異同之處隔盛。
./gffcompare -r file/GCA_000977205.2_Sc_YJM1338_v1_genomic.gff -o compare file/outresult6_pl.gff3
查看compare.stats結(jié)果
從compare.stats的結(jié)果看出堿基水平及基因水平結(jié)果比較好，外顯子水平及轉(zhuǎn)錄本水平比較差拾稳，但可信度高吮炕；內(nèi)含子結(jié)果差且可信度不高。從結(jié)果看出访得，同源搜索在預測結(jié)構(gòu)方面并不是很好龙亲。
附錄
關(guān)于Sensitivity和Precision
計算公式如下：
Sensitivity = TP / (TP+FN)
Precision = TP / (TP+FP)
例如，現(xiàn)在從一堆有好有壞的西瓜中挑好瓜悍抑，好壞各50個鳄炉。二狗認為所有的都是好瓜，則TP搜骡、FP=50,則Sensitivity是100%拂盯，但這可信嗎？很明顯瓜中有壞瓜记靡，所以不可信谈竿，此時需要考慮Precision，Precision=50%,而小明認為其中１個好瓜（這一個確實是好瓜）摸吠，99個壞瓜空凸，則,Sensitivity ＝1.96%,Precision＝100%，預測的好瓜特別對蜕便，但數(shù)量少。所以需要綜合考慮Sensitivity和Precision 贩幻。
