水稻基因组变异注释

太玄学了，仅供参考和日后借鉴。

运行命令：

###生成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 

okk，前9条染色体也出来了

接下来是提取相应区间的变异信息:

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

改成csv文件

借鉴，防遗忘笔记：
annovar基因（水稻）注释 - 知乎 (zhihu.com)
用annovar对vcf（SNP&INDEL）文件进行注释 - 简书 (jianshu.com)