annovar在肿瘤项目中的使用

最近在做个乳腺癌的项目,用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!

image.png
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 参数后即可。


jianshujietu1.png

我觉得不好用的一点是,它没有将下载失败的原因告诉我们,我猜测是因为文件太大了,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官方文档中也有 详细的介绍 )

最近持续修改更新中,有没有什么问题呀?

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 159,716评论 4 364
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,558评论 1 294
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 109,431评论 0 244
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 44,127评论 0 209
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,511评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,692评论 1 222
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,915评论 2 313
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,664评论 0 202
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,412评论 1 246
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,616评论 2 245
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,105评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,424评论 2 254
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,098评论 3 238
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,096评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,869评论 0 197
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,748评论 2 276
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,641评论 2 271

推荐阅读更多精彩内容