1. 准備文件：
- ref.fa
- ref.gtf或者gff3，最好是gtf3，可將gff3轉化為gtf
- sample.vcf
2. 用gff3ToGenePred與gtfToGenePred工具將gtf或gff3文件轉化為reference_refGene.txt （軟件來自http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/）
gtfToGenePred.dms -genePredExt ref.gtf SP_refGene.txt &
gtf：
SpoScf_00032 maker exon 12508 13665 . + . transcript_id "Spo06120"; gene_id "Spo06120";
SpoScf_00032 maker exon 14070 17062 . + . transcript_id "Spo06120"; gene_id "Spo06120";
SpoScf_00032 maker exon 17626 17899 . + . transcript_id "Spo06120"; gene_id "Spo06120";
SpoScf_00032 maker exon 17979 18066 . + . transcript_id "Spo06120"; gene_id "Spo06120";
3. 將ref.fa文件轉化為SP_refGeneMrna.fa
1 perl retrieve_seq_from_fasta.pl --format refGene --seqfile ref.fa SP_refGene.txt Sp_refGeneMrna.fa
4. 再將vcf文件轉化為annovar格式
1 perl convert2annovar.pl -includeinfo -allsample -withfreq -format vcf4 sample.VCF >sample.avinput 2 3 4 5 6 ## 7 --includeinfo: 輸出文件含有特定額外的信息 8 --allsample： 多樣本的vcf，輸出多個樣本的結果 9 --withfreq： 輸出文件包含頻率信息 10 --format： 輸入文件格式
5. 用table_annovar.pl進行注釋（可一次性完成三種類型的注釋, 本次只有基於基因)
1 perl ../table_annovar.pl test.avinput sp/ --buildver SP --outfile myanno --protocol refGene --operation g 2 3 ##參數 4 sp: 含有SP_refGeneMrna.fa和SP_refGene.txt文件夾 5 --buildver: 基因組建立的版本 6 --outfile： 輸出文件前綴 7 --protocol： 逗號分隔的注釋流程，代表庫的名字 8 --operation： g（gene），r（region），f（filter）
最終得到兩個注釋文件文件和一個log文件exonic_variant_function和variant_function
關注下方公眾號可獲得更多精彩
