gatk的germline mutation call完之后,得到VQSR过滤之后的VCF文件,之前使用了gvcf模式,这里是多个样本合并起来的的文件。
首先官网下载ANNOVAR软件,需要注册。官网地址
解压后里面有很多.pl文件,这些是用perl写的脚本,可以直接运行,类似于:
perl table_annovar.pl +参数
初级主要用到的有
annotate.pl 下载数据库,注释数据
convert2annovar.pl 将变异文件转化annovar可以使用的文件格式
annotate_variation.pl 一次注释一个数据库
table_annovar.pl 一次可以对多个数据库进行注释
下面下载需要的注释数据库文件
ANNOVAR下载的时候自带了一部分参考基因组注释文件,如hg19_refGene等。通常需要下载其它的数据库,有gene-based annotation,filter-based annotation数据库用户贡献数据库以及第三方数据库,由于我这次没有用到,也没有去了解其它数据库。
这些数据库在ANNOVAR官网都有列出来,我自己常用的有:
ljb26_all,dbnsfp35c,exac03,1000g2015aug,avsnp150,clinvar_20200316,cytoBand
下载格式类似于:
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar ljb26_all humandb/
由于我的是gvcf模式得到的多个样本的vcf文件,首先拆分转换为annovar的输入文件,代码类似于:
convert2annovar.pl -format vcf4 indir/${cohort}.passed.vcf.gz -outfile outdir/${cohort} -allsample
其中在我这里cohort为样本组的名字
然后使用table_annovar.pl同时进行多个数据库的注释,代码类似于:
samples=$(ls ${outdir} | grep ${cohort} )
for sample in $samples
do
echo $sample
id=$(echo $sample | awk -F_ '{print $1}') #这个因样本名而异
echo $id
perl table_annovar.pl $outdir/${sample} \
humandb/ -buildver hg19 \
-out myanno -remove -protocol refGene,cytoBand,exac03,avsnp150,dbnsfp35c \
-operation g,r,f,f,f -nastring . -vcfinput \
#g,r,f,f,f 代表对refGene,cytoBand,exac03,avsnp150,dbnsfp35c这五个数据库分别进行gene-based annotation,region based annotation, filter-based annotation
-polish --outfile $outdir/${id} && echo "$id done"
done
最终会得到id为前缀的注释文件,即id.hg19_multianno.txt
用于下一步分析。