annovar注释的进阶使用

生信菜鸟团

发布于 2022-05-24 16:30:57

3.3K0

发布于 2022-05-24 16:30:57

文章被收录于专栏：生信菜鸟团

菜鸟团公众号肯定讲过annovar的使用了。比如Nickier的vcf文件的注释及ANNOVAR的使用。

而在使用 ANNOVAR 之前，你应该知道和ANNOVAR 是如何注释 RS ID 的？中，鲍志炜同学也深入介绍了一下annovar的运行机制与注意事项。

在本文中，我主要介绍我在使用annovar时遇到的需求和这些需求我是如何实现的。

两个示例文件为，test.vcf.gz（单样本）和test_multi.vcf.gz（多样本）。

$ bcftools view -H test.vcf.gz
14      23240713        .       T       TAGC    26248.7 PASS    AC=1;AN=2       GT      0/1
6       12296255        .       G       T       34219.8 .       AC=1;AN=2       GT      0/1
22      44322970        .       G       T       11568.8 .       AC=1;AN=2       GT      0/1
$ bcftools view -H test_multi.vcf.gz
14      23240713        .       T       TAGC    26248.7 PASS    AC=3;AN=6       GT      0/1     ./.     0/1     ./.     0/1
6       12296255        .       G       T       34219.8 .       AC=2;AN=2       GT      ./.     ./.     1/1     ./.     ./.
22      44322970        .       G       T       11568.8 .       AC=1;AN=2       GT      ./.     ./.     0/1     ./.     ./.

annovar最常使用的方式是这样的。

$ perl ~/software/annovar/table_annovar.pl test.vcf.gz ~/project/0-reference/3-Annovar/hg19 -buildver hg19 \
-out test -remove -protocol refGene \
-operation g -nastring . -polish -vcfinput

这里，仅使用了refGene数据库。

运行完毕后，程序生成了test.avinput，test.hg19_multianno.txt和test.hg19_multianno.vcf三个文件。

如果我们想查看对应的注释，显然，可以在test.hg19_multianno.vcf里查看。

$ bcftools view -H test.hg19_multianno.vcf
14      23240713        .       T       TAGC    26248.7 PASS    AC=1;AN=2;ANNOVAR_DATE=2020-06-08;Func.refGene=exonic;Gene.refGene=OXA1L;GeneDetail.refGene=.;ExonicFunc.refGene=nonframeshift_insertion;AAChange.refGene=OXA1L:NM_005015:exon10:c.1254_1255insAGC:p.S422_K423insS;ALLELE_END   GT      0/1
6       12296255        .       G       T       34219.8 .       AC=1;AN=2;ANNOVAR_DATE=2020-06-08;Func.refGene=exonic;Gene.refGene=EDN1;GeneDetail.refGene=.;ExonicFunc.refGene=nonsynonymous_SNV;AAChange.refGene=EDN1:NM_001168319:exon5:c.G591T:p.K197N,EDN1:NM_001955:exon5:c.G594T:p.K198N;ALLELE_END      GT      0/1
22      44322970        .       G       T       11568.8 .       AC=1;AN=2;ANNOVAR_DATE=2020-06-08;Func.refGene=exonic;Gene.refGene=PNPLA3;GeneDetail.refGene=.;ExonicFunc.refGene=nonsynonymous_SNV;AAChange.refGene=PNPLA3:NM_025225:exon2:c.G343T:p.G115C;ALLELE_END  GT     0/1

但是vcf文件并不适合读取分析。

相比较而言，test.hg19_multianno.txt，作为一个tsv，是一个更友好的处理选择。

$ cat test.hg19_multianno.txt
Chr     Start   End     Ref     Alt     Func.refGene    Gene.refGene    GeneDetail.refGene      ExonicFunc.refGene      AAChange.refGenOtherinfo1       Otherinfo2      Otherinfo3      Otherinfo4      Otherinfo5      Otherinfo6      Otherinfo7      Otherinfo8      Otherinfo9      Otherinfo10     Otherinfo11     Otherinfo12     Otherinfo13
14      23240713        23240713        -       AGC     exonic  OXA1L   .       nonframeshift insertion OXA1L:NM_005015:exon10:c.1254_1255insAGC:p.S422_K423insS        0.5     26248.7 .       14      23240713        .       T       TAGC    26248.7 PASS    AC=1;AN=2      GT       0/1
6       12296255        12296255        G       T       exonic  EDN1    .       nonsynonymous SNV       EDN1:NM_001168319:exon5:c.G591T:p.K197N,EDN1:NM_001955:exon5:c.G594T:p.K198N    0.5     34219.8 .       6       12296255        .       G       T       34219.8 .      AC=1;AN=2        GT      0/1
22      44322970        44322970        G       T       exonic  PNPLA3  .       nonsynonymous SNV       PNPLA3:NM_025225:exon2:c.G343T:p.G115C  0.5     11568.8 .       22      44322970        .       G       T       11568.8 .       AC=1;AN=2       GT      0/1

不过这个文件不够完美。让我们来改进一下。

突变格式

首先，是test.hg19_multianno.txt中的突变格式问题，14号染色体的T到TAGC的突变，在test.hg19_multianno.txt里居然变成了，-到AGC。虽然我能理解-到AGC和T到TAGC是一样的，但是这样操作之后，当需要将test.hg19_multianno.txt与test.vcf.gz产出的突变信息文件进行匹配的时候，就要大伤一番脑筋。

其实annovar提供了对indel不进行处理的方式，-keepindelref参数。

但是很怪，这个参数不能在table_annovar.pl里直接使用，所以注释分成更规范的两步。

先使用convert2annovar.pl生成test.avinput，这里添加-keepindelref参数，拒绝annovar对indel的修改，

$ perl ~/software/annovar/convert2annovar.pl -keepindelref  -format vcf4 test.vcf.gz -outfile test.avinput
$ cat test.avinput
14      23240713        23240713        T       TAGC    het     26248.7 .
6       12296255        12296255        G       T       het     34219.8 .
22      44322970        44322970        G       T       het     11568.8 .

然后注释，此时注释就不需要-vcfinput参数了。

$ perl ~/software/annovar/table_annovar.pl test.avinput ~/project/0-reference/3-Annovar/hg19 -buildver hg19 \
-out test -remove -protocol refGene \
-operation g -nastring . -polish
$ cat test.hg19_multianno.txt
Chr     Start   End     Ref     Alt     Func.refGene    Gene.refGene    GeneDetail.refGene      ExonicFunc.refGene      AAChange.refGene
14      23240713        23240713        T       TAGC    exonic  OXA1L   .       nonframeshift substitution      OXA1L:NM_005015:exon10:c.1254delinsTAGC:p.S422_K423insS
6       12296255        12296255        G       T       exonic  EDN1    .       nonsynonymous SNV       EDN1:NM_001168319:exon5:c.G591T:p.K197N,EDN1:NM_001955:exon5:c.G594T:p.K198N
22      44322970        44322970        G       T       exonic  PNPLA3  .       nonsynonymous SNV       PNPLA3:NM_025225:exon2:c.G343T:p.G115C

可以了。

多样本

直接使用table_annovar.pl注释多样本vcf时，会产生Otherinfo列非常多的问题。

$ perl ~/software/annovar/table_annovar.pl test_multi.vcf.gz ~/project/0-reference/3-Annovar/hg19 -buildver hg19 \
-out test_multi -remove -protocol refGene \
-operation g -nastring . -polish -vcfinput
$ cat test_multi.hg19_multianno.txt
Chr     Start   End     Ref     Alt     Func.refGene    Gene.refGene    GeneDetail.refGene      ExonicFunc.refGene      AAChange.refGenOtherinfo1       Otherinfo2      Otherinfo3      Otherinfo4      Otherinfo5      Otherinfo6      Otherinfo7      Otherinfo8      Otherinfo9      Otherinfo10     Otherinfo11     Otherinfo12     Otherinfo13     Otherinfo14     Otherinfo15     Otherinfo16     Otherinfo17
14      23240713        23240713        -       AGC     exonic  OXA1L   .       nonframeshift insertion OXA1L:NM_005015:exon10:c.1254_1255insAGC:p.S422_K423insS        0.5     26248.7 .       14      23240713        .       T       TAGC    26248.7 PASS    AC=3;AN=6      GT       0/1     ./.     0/1     ./.     0/1
6       12296255        12296255        G       T       exonic  EDN1    .       nonsynonymous SNV       EDN1:NM_001168319:exon5:c.G591T:p.K197N,EDN1:NM_001955:exon5:c.G594T:p.K198N    1       34219.8 .       6       12296255        .       G       T       34219.8 .      AC=2;AN=2        GT      ./.     ./.     1/1     ./.     ./.
22      44322970        44322970        G       T       exonic  PNPLA3  .       nonsynonymous SNV       PNPLA3:NM_025225:exon2:c.G343T:p.G115C  0.5     11568.8 .       22      44322970        .       G       T       11568.8 .       AC=1;AN=2       GT      ./.     ./.    0/1      ./.     ./.

五个样本已经很长了，如果是成百上千个样本，根本不敢想象。

而使用convert2annovar.pl两步走的话，Otherinfo列被去除，避免了此问题。

但是多样本使用convert2annovar.pl会有一些新的问题。

$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4 test_multi.vcf.gz -outfile test_multi.avinput
WARNING to old ANNOVAR users: this program no longer does line-to-line conversion for multi-sample VCF files. If you want to include all variants in output, use '-format vcf4old' or use '-format vcf4 -allsample -withfreq' instead.
...
$ cat test_multi.avinput
14      23240713        23240713        T       TAGC    het     26248.7 .

在产生的test_multi.avinput中只有一个突变。为什么呢？注意到annovar给出了warning。如果要使用多样本vcf，必须使用-format vcf4old或-format vcf4 -allsample -withfreq。

重新试一下。

$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4old test_multi.vcf.gz -outfile test_multi.avinput
$ cat test_multi.avinput
14      23240713        23240713        -       AGC     het     26248.7
6       12296255        12296255        G       T       unknown 34219.8
22      44322970        44322970        G       T       unknown 11568.8
$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4 -allsample -withfreq test_multi.vcf.gz -outfile test_multi.avinput
$ cat test_multi.avinput
14      23240713        23240713        T       TAGC    0.5     26248.7 .
6       12296255        12296255        G       T       1       34219.8 .
22      44322970        44322970        G       T       0.5     11568.8 .

两者稍有区别，但都输出了所有突变。

不过，这时候，如果还想包含vcf里的一些信息的话。比如-withfilter。

$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4old test_multi.vcf.gz -outfile test_multi.avinput -withfilter
$ cat test_multi.avinput
14      23240713        23240713        -       AGC     het     PASS    26248.7
6       12296255        12296255        G       T       unknown .       34219.8
22      44322970        44322970        G       T       unknown .       11568.8

这-keepindelref怎么失效了。

$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4 -allsample -withfreq test_multi.vcf.gz -outfile test_multi.avinput
$ cat test.avinput
14      23240713        23240713        T       TAGC    0.5     26248.7 .
6       12296255        12296255        G       T       0.5     34219.8 .
22      44322970        44322970        G       T       0.5     11568.8 .

也没有，看来这-withfilter和-keepindelref不能既要还要啊。

添加自有数据库

官方的数据库列表包含了很多常用数据库，但是有时候会需要一些比较特殊的信息。比如亚洲人群的MAF信息。

在找到数据之后，就想利用annovar的注释机制将此信息也添加到注释结果中去。

看一下annovar数据库中gnomAD数据的样子，仿一下。

$ head -2 ~/project/0-reference/3-Annovar/hg19/hg19_gnomad211_exome.txt
#Chr    Start   End     Ref     Alt     AF      AF_popmax       AF_male AF_female       AF_raw  AF_afr  AF_sas  AF_amr  AF_eas  AF_nfe AF_fin   AF_asj  AF_oth  non_topmed_AF_popmax    non_neuro_AF_popmax     non_cancer_AF_popmax    controls_AF_popmax
1       12198   12198   G       C       .       .       .       .       0.0457  .       .       .       .       .       .       .      .

第一行列名，第二行主要是对应的Chr Start End Ref Alt值，剩下的都是AF。应该挺好仿的。

我下载的MAF信息来自于WBBC[1]。

下载各染色体数据，然后合并。

$ bcftools concat WBBC.chr{1..22}.GRCh37.vcf.gz WBBC.chr{X,Y}.GRCh37.vcf.gz -Oz -o ~/WBBC.GRCh37.vcf.gz

normalize，然后增加end信息。

$ hg19=~/project/0-reference/2-GATK/hg19/human_g1k_v37_decoy.fasta
$ bcftools norm -m-both -f ${hg19} WBBC.GRCh37.vcf.gz | bcftools +fill-tags - -Oz -o WBBC.GRCh37.processed.vcf.gz -- -t END
[W::vcf_parse_info] INFO 'NS' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'North_AF' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'North_AN' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'Central_AF' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'Central_AN' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'South_AF' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'South_AN' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'Lingnan_AF' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'Lingnan_AN' is not defined in the header, assuming Type=String
[W::vcf_parse_info] INFO 'RR' is not defined in the header, assuming Type=String
[E::vcf_format] Invalid BCF, the INFO tag id=59 is too large at 1:12772
[flush_buffer] Error: cannot write to -

咦，出错了，大概是因为这些INFO没有定义，导致出错。改一下。

生成header。

$ bcftools view -h WBBC.GRCh37.vcf.gz > header

在header里照葫芦画瓢添加一下。

##INFO=<ID=NS,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=North_AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=North_AN,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=Central_AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=Central_AN,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=South_AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=South_AN,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=Lingnan_AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=Lingnan_AN,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=RR,Number=A,Type=Float,Description="Allele Frequency">

保存。修改vcf的header，然后处理。

$ bcftools reheader -h header WBBC.GRCh37.vcf.gz -o WBBC.GRCh37.reheader.vcf.gz
$ bcftools norm -m-both -f ${hg19} WBBC.GRCh37.reheader.vcf.gz | bcftools +fill-tags -Oz -o WBBC.GRCh37.processed.vcf.gz -- -t END

看一下数据。

$ bcftools view -H WBBC.GRCh37.processed.vcf.gz | head -2
1       12772   .       A       G       68.74   PASS    AC=1;AF=0.000111607;AN=8960;NS=4480;North_AF=0;North_AN=448;Central_AF=0;Central_AN=100;South_AF=0;South_AN=8070;Lingnan_AF=0;Lingnan_AN=126;RR=4479;DP=40710;VQSLOD=1.91;END=12772
1       12795   .       C       T       46.66   PASS    AC=1;AF=0.000111607;AN=8960;NS=4480;North_AF=0;North_AN=448;Central_AF=0;Central_AN=100;South_AF=0.000123916;South_AN=8070;Lingnan_AF=0;Lingnan_AN=126;RR=4479;DP=44411;VQSLOD=2.94;END=12795

最需要的是AF，提取需要的突变信息和AF列。

$ echo -e '#Chr\tStart\tEnd\tRef\tAlt\tWBBC' > WBBC_AF.txt
$ bcftools query -f '%CHROM\t%POS\t%END\t%REF\t%ALT\t%AF\n' WBBC.GRCh37.processed.vcf.gz >> hg19_WBBC_AF.txt
$ head -5 WBBC_AF.txt
#Chr    Start   End     Ref     Alt     WBBC
1       12772   12772   A       G       0.000111607
1       12795   12795   C       T       0.000111607
1       12807   12807   C       T       0.0371652
1       12839   12839   G       C       0.00368304

像模像样的。但是等等，官方数据库好像有个对应的idx文件作为index。自有的这个数据库也得搞一个。

经过仔细查找，找到了制作index的脚本[2]。

下载之后使用。

$ wget https://github.com/WGLab/doc-ANNOVAR/files/6670482/index_annovar.txt -O ~/software/annovar/index_annovar.pl
$ perl ~/software/annovar/index_annovar.pl WBBC_AF.txt --outfile ~/project/0-reference/3-Annovar/hg19/hg19_WBBC.txt

试一下。

$ perl ~/software/annovar/convert2annovar.pl -keepindelref -format vcf4 -allsample -withfreq test_multi.vcf.gz -outfile test_multi.avinput
$ perl ~/software/annovar/table_annovar.pl test.avinput ~/project/0-reference/3-Annovar/hg19 -buildver hg19 \
-out test -remove -protocol WBBC \
-operation f -nastring . -polish
$ cat test.hg19_multianno.txt
Chr     Start   End     Ref     Alt     WBBC
14      23240713        23240713        T       TAGC    0.108147
6       12296255        12296255        G       T       0.279353
22      44322970        44322970        G       T       0.0487723

真不错。