最近在做个乳腺癌的项目,用GATK call 变异后得到vcf,需要做一些注释,就要用到annovar,
ANNOVAR支持三种不同形式的注释: gene-based, region-based 和filter-based. 这三种注释分别针对于每一个variant的不同方面:基于基因的注释(gene-based annotation)揭示variant与已知基因直接的关系以及对其产生的功能性影响;基于区域的注释(region-based annotation)揭示variant 与不同基因组特定段的关系,例如:它是否落在已知的保守基因组区域;基于过滤子的注释( filter-based annotation )则给出这个variant的一系列信息,如: population frequency in different populations 和various types of variant-deleteriousness prediction scores, 这些可被用来过滤掉一些公共的及 probably(大概,肯定的成分较大,是most likely) nondeleterious variants。
1,annovar下载,http://www.openbioinformatics.org/annovar/annovar_download_form.php
还要注册,填了邮箱,最新版的下载地址就会发送到邮箱,就几个perl程序和自带的几个简单的database,下载很快的。
不过现在填写的邮箱如果是以.com/.co 结尾的,都会被要求购买license后才能使用。sad!
ANNOVAR │ annotate_variation.pl #主程序,功能包括下载数据库,三种不同的注释
│ coding_change.pl #可用来推断蛋白质序列
│ convert2annovar.pl #将多种格式转为.avinput的程序
│ retrieve_seq_from_fasta.pl #用于自行建立其他物种的转录本
│ table_annovar.pl #注释程序,可一次性完成三种类型的注释
│ variants_reduction.pl #可用来更灵活地定制过滤注释流程
| ─example #存放示例文件
│─humandb #人类注释数据库
2,开始先将vcf文件转化成avinput格式,
perl convert2annovar.pl -format vcf4 example.vcf > example.avinput
3,用table_annovar.pl注释'example.avinput'得到注释文件 'example.hg19_multianno.txt',为什么用table_annovar?用the ‘table_annovar.pl’ 来注释variants(可一次性完成三种类型的注释)。允许在同一命令中用输出的特定顺序来对多个注释类型进行自定义选择(custom selection)。
输入下列命令,用之前下载好的注释数据库来注释vcf格式文件中的variants:
perl table_annovar.pl example.avinput humandb/ -buildver hg19 -out example -remove -protocol refGene,cytoBand,popfreq_all_20150413,dbnsfp33a,dbscsnv11,clinvar_20170905,intervar_20170202,cosmic82_coding,cosmic82_noncoding,icgc21,avsnp147 -operation g,r,f,f,f,f,f,f,f,f,f -nastring . -polish
# example.avinput 参考(refers to )输入的vcf文件的名称
# humandb/ 表示你前面下载的数据库的绝对路径(完整路径)、
# -buildver 后面接数据库的bulid,比如hg19/hg18/hg38
# -remove 表示把注释过程中产生的中间文件在结束时删除掉、
# -nastring . 表示没有值的项用点来代替
# -out 选项是指定输出文件的前缀,输出到指定目录时可以跟路径+前缀。
# -vcfinput 表明输入的是vcf文件,那么输出文件也是vcf文件,否则只会生成一个txt文件
# -protocol 表示注释使用的数据库的准确名称,多个用逗号隔开,且要注意顺序
# -operation 表示对应顺序的数据库的类型(g代表gene-based、r代表region-based、f代表filter-based),用逗号隔开,注意顺序
# -csvout 表示最后输出.csv文件
关键步骤( CR ITICAL STEP):
1、确保注释数据库的名称正确并且是按你想要在输出文件中显示的顺序排列的;
2、确保 ‘-operation’指定的注释类型顺序和‘-protocol’指定的数据库顺序是一致的;
3、确保每个protocal名称或注释类型之间只有一个逗号,并且没有空白。
‘final.hg19_multianno.vcf’.输出文件应该是以个VCF格式文件,INFO那列以 ‘key=value’ 形式、 ‘;’分割成几个小区域. eg:‘Func.refGene=intronic;Gene.refGene=SAMD11’. 每个键值对代表一个ANNOVAR注释信息。输出文件可以用为VCF格式文件设计的基因分析软件进一步处理。
‘final.hg19_multianno.txt’. 每一行代表一个variant 。用tab分隔,多余列为加上的注释信息,顺序按 ‘--protocol’ 选项所设定的注释类型argument。
这个时候发现命令里面的数据库没有下载下来,那就一一下载一下:
缺少的数据库有:cytoBand,popfreq_all_20150413,dbnsfp33a,dbscsnv11,clinvar_20170905,intervar_20170202,cosmic82_coding,cosmic82_noncoding,icgc21,avsnp147,
下载使用的是annovar 的命令行方法:
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar clinvar_20170905 humandb/
这里列出了如下命令可供下载的database:http://annovar.openbioinformatics.org/en/latest/user-guide/download/
也可以用avdblist参数查看,perl annotate_variation.pl -buildver hg19 -downdb avdblist -webfrom annovar ./
没有的数据库可以更换webfrom参数:比如cytoBand,就可以不用加 -webfrom 参数,
perl annotate_variation.pl -buildver hg19 -downdb clinvar_20170905 humandb/
--downdb
download annotation databases from UCSC Genome Browser, Ensembl,
1000 Genomes Project, ANNOVAR website or other resources. The
annotation databases are required for functional annotation of
genetic variants.
--webfrom
specify the source of database (ucsc or annovar or URL) in the
downdb operation. By default, files from UCSC Genome Browser
annotation database will be downloaded.
也可以提供url下载,只要将网页链接提供在-webfrom 参数后即可。
我觉得不好用的一点是,它没有将下载失败的原因告诉我们,我猜测是因为文件太大了,5.3G,只好在浏览器上下载。http://www.openbioinformatics.org/annovar/download/hg19_dbnsfp33a.txt.gz
公司网速受限,也下载了大半天,下载后解压,将文件上传到humandb中。
还有个头疼的数据库:
http://annovar.openbioinformatics.org/en/latest/user-guide/filter/#cosmic-annotations
报错:
Can't exec "annotate_variation.pl": Permission denied at $path/table_annovar.pl line 391.
还要改一下权限:
chmod 755 *.pl
最后重新注释:
perl table_annovar.pl example.avinput humandb/ -buildver hg19 -out example -remove -protocol refGene,cytoBand,popfreq_all_20150413,dbnsfp33a,dbscsnv11,clinvar_20170905,intervar_20170202,cosmic82_coding,cosmic82_noncoding,icgc21,avsnp147 -operation g,r,f,f,f,f,f,f,f,f,f -nastring . -polish
比较大的几个数据库跑起来还挺慢的。一天左右,全部跑完,
使用另外一个脚本variant_classifier1.0.py 得到变异的具体分类信息'result.txt',这个脚本也用到了一些数据库:domain_brca.txt iso_brca.txt lof.genes mis_brca.txt nifty_brca.txt pathogenic_brca.txt,这些是自己准备的。将结果整合到一起:
python variant_classifier1.0.py example.hg19_multianno.txt pathogenic_brca.txt iso_brca.txt mis_brca.txt domain_brca.txt nifty_brca.txt lof.genes result.txt
最后,就算注释了这么多的数据库,到最后,能看懂数据库的又有谁呢,应该写一篇关于各种数据库的文章。
使用annotate_variation.pl进行注释:
Gene-based Annotation(基于基因的注释)
基于基因的注释(gene-based annotation)揭示variant与已知基因直接的关系以及对其产生的功能性影响,需要使用 for gene-based 的数据库。
perl annotate_variation.pl -geneanno -dbtype refGene -out ex1 -buildver hg19 ex1.avinput humandb/
geneanno 表示使用基于基因的注释
# -dbtype refGene 表示使用"refGene"数据库
# -out ex1 表示输出文件以ex1为前缀
Region-based Annotation(基于区域的注释)
基于过滤的注释精确匹配查询变异与数据库中的记录:如果它们有相同的染色体,起始位置,结束位置,REF的等位基因和ALT的等位基因,才能认为匹配。基于区域的注释看起来更像一个区域的查询(这个区域也可以是一个单一的位点),在一个数据库中,它不在乎位置的精确匹配,它不在乎核苷酸的识别。
基于区域的注释(region-based annotation)揭示variant与不同基因组特定段的关系,例如:它是否落在已知的保守基因组区域。基于区域的注释的数据库一般由UCSC提供。
perl annotate_variation.pl -regionanno -buildver hg19 -dbtype cytoBand ex1.avinput humandb/
perl annotate_variation.pl -regionanno -buildver hg19 -dbtype gff3 -gff3dbfile tfbs.gff3 ex1.avinput humandb/
Filter-based Annotation(基于过滤的注释)
filter-based和region-based主要的区别是,filter-based针对mutation(核苷酸的变化)而region-based针对染色体上的位置。例如region-based比对chr1:1000-1000而filter-based比对chr1:1000-1000上的A->G。
基于过滤的注释,使用不同的过滤数据库,可以给出这个variant的一系列信息。如在全基因组数据中的变异频率,可使用1000g2015aug、kaviar_20150923等数据库;在全外显组数据中的变异频率,可使用exac03、esp6500siv2等;在孤立的或者低代表人群中的变异频率,可使用ajews等数据库。(在ANNOVAR官方文档中也有 详细的介绍 )
最近持续修改更新中,有没有什么问题呀?