太玄学了,仅供参考和日后借鉴。
运行命令:
###生成refgene
~/tools/gtfToGenePred -genePredExt Oryza_sativa.IRGSP-1.0.50.gtf gcf_refGene.txt # *_refGene.txt,*
###fa
~/tools/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile IRGSP-1.0_genome.fasta gcf_refGene.txt --outfile gcf_refGeneMrna.fa
### 转annovar
~/tools/annovar/convert2annovar.pl -format vcf4 sample.vcf > sample.annovar
###注释命令
$> /public/home/fengting/tools/annovar/annotate_variation.pl -buildver gcf -geneanno -outfile sample.anno sv.annovar ricedb/
我看了下,只注释了后三条染色体,可能是染色体编号的原因,把前9条染色体名改一下:
sed -i "s/chr0/chr/g" sv.annovar
![](https://img.haomeiwen.com/i25274977/fd219fa7ffb79121.png)
接下来是提取相应区间的变异信息:
cat pav-gene-pos.txt | while read chr start end id
do
##awk '{if($3 == "'"$chr"'" && $3 == "gene" && $4 >= "'"${start}"'" && $5 <= "'"${end}"'")print } ' sample.anno.variant_function > tmp1.txt
awk '{if($3 == "'"$chr"'" && $4 >= "'"${start}"'" && $5 <= "'"${end}"'")print } ' sample.anno.variant_function > tmp1.txt
cat tmp1.txt |while read line1
do
echo $line1 >> res/"${id}".csv
done
##awk '$1 == "1" && $3 == "gene" && $4 >= 10000 && $5 <= 200000 {print $0} ' Arabidopsis_thaliana.TAIR10.31.gff3 > out.gff
done
![](https://img.haomeiwen.com/i25274977/ab0985206e754efd.png)
借鉴,防遗忘笔记:
annovar基因(水稻)注释 - 知乎 (zhihu.com)
用annovar对vcf(SNP&INDEL)文件进行注释 - 简书 (jianshu.com)
网友评论