It's all well and good to have a bunch of variant sites and genotypes, but how can you know whether they might be biologically relevant? That's where functional annotation comes in; this is the process of adding annotations to variant sites that tell you whether they fall within a coding or regulatory region, whether the mutation is sense or nonsense, etc. The GATK does not include any tool that perform functional annotation as such, but it does include conversion tools that allow you to interface with the popular annotation program SnpEff.
Our testing has shown that not all combinations of snpEff/database versions produce high-quality results. Be sure to read this document completely to familiarize yourself with our recommended best practices BEFORE running snpEff.
Until recently we were using an in-house annotation tool for genomic annotation, but the burden of keeping the database current and our lack of ability to annotate indels has led us to employ the use of a third-party tool instead. After reviewing many external tools (including annoVar, VAT, and Oncotator), we decided that SnpEff best meets our needs as it accepts VCF files as input, can annotate a full exome callset (including indels) in seconds, and provides continually-updated transcript databases. We have implemented support in the GATK for parsing the output from the SnpEff tool and annotating VCFs with the information provided in it.
Download the SnpEff core program. If you want to be able to run VariantAnnotator on the SnpEff output, you'll need to download a version of SnpEff that VariantAnnotator supports from this page (currently supported versions are listed below). If you just want the most recent version of SnpEff and don't plan to run VariantAnnotator on its output, you can get it from here.
After unzipping the core program, open the file snpEff.config in a text editor, and change the "database_repository" line to the following:
database_repository = http://sourceforge.net/projects/snpeff/files/databases/
Then, download one or more databases using SnpEff's built-in download command:
java -jar snpEff.jar download GRCh37.64
You can find a list of available databases here. The human genome databases have GRCh or hg in their names. You can also download the databases directly from the SnpEff website, if you prefer.
The download command by default puts the databases into a subdirectory called data within the directory containing the SnpEff jar file. If you want the databases in a different directory, you'll need to edit the
data_dir entry in the file
snpEff.config to point to the correct directory.
Run SnpEff on the file containing your variants, and redirect its output to a file. SnpEff supports many input file formats including VCF 4.1, BED, and SAM pileup. Full details and command-line options can be found on the SnpEff home page.
If you want to take advantage of SnpEff integration in the GATK, you'll need to run SnpEff version **2.0.5*. Note: newer versions are currently unsupported by the GATK, as we haven't yet had the reources to test it.
These best practices are based on our analysis of various snpEff/database versions as described in detail in the Analysis of SnpEff Annotations Across Versions section below.
We recommend using only the GRCh37.64 database with SnpEff 2.0.5. The more recent GRCh37.65 database produces many false-positive Missense annotations due to a regression in the ENSEMBL Release 65 GTF file used to build the database. This regression has been acknowledged by ENSEMBL and is supposedly fixed as of 1-30-2012; however as we have not yet tested the fixed version of the database we continue to recommend using only GRCh37.64 for now.
We recommend always running with
-onlyCoding true with human databases (eg., the GRCh37.* databases). Setting
-onlyCoding false causes snpEff to report all transcripts as if they were coding (even if they're not), which can lead to nonsensical results. The
-onlyCoding false option should only be used with databases that lack protein coding information.
Do not trust annotations from versions of snpEff prior to 2.0.4. Older versions of snpEff (such as 2.0.2) produced many incorrect annotations due to the presence of a certain number of nonsensical transcripts in the underlying ENSEMBL databases. Newer versions of snpEff filter out such transcripts.
See our analysis of the SNP annotations produced by snpEff across various snpEff/database versions here.
Both snpEff 2.0.2 + GRCh37.63 and snpEff 2.0.5 + GRCh37.65 produce an abnormally high Missense:Silent ratio, with elevated levels of Missense mutations across the entire spectrum of allele counts. They also have a relatively low (~70%) level of concordance with the 1000G Gencode annotations when it comes to Silent mutations. This suggests that these combinations of snpEff/database versions incorrectly annotate many Silent mutations as Missense.
snpEff 2.0.4 RC3 + GRCh37.64 and snpEff 2.0.5 + GRCh37.64 produce a Missense:Silent ratio in line with expectations, and have a very high (~97%-99%) level of concordance with the 1000G Gencode annotations across all categories.
See our comparison of SNP annotations produced using the GRCh37.64 and GRCh37.65 databases with snpEff 2.0.5 here
The GRCh37.64 database gives good results on the condition that you run snpEff with the
-onlyCoding true option. The
-onlyCoding false option causes snpEff to mark all transcripts as coding, and so produces many false-positive Missense annotations.
The GRCh37.65 database gives results that are as poor as those you get with the
-onlyCoding false option on the GRCh37.64 database. This is due to a regression in the ENSEMBL release 65 GTF file used to build snpEff's GRCh37.65 database. The regression has been acknowledged by ENSEMBL and is due to be fixed shortly.
See our analysis of the INDEL annotations produced by snpEff across snpEff/database versions here
Below is an example of how to run SnpEff version 2.0.5 with a VCF input file and have it write its output in VCF format as well. Notice that you need to explicitly specify the database you want to use (in this case, GRCh37.64). This database must be present in a directory of the same name within the
data_dir as defined in
java -Xmx4G -jar snpEff.jar eff -v -onlyCoding true -i vcf -o vcf GRCh37.64 1000G.exomes.vcf > snpEff_output.vcf
In this mode, SnpEff aggregates all effects associated with each variant record together into a single INFO field annotation with the key EFF. The general format is:
EFF=Effect1(Information about Effect1),Effect2(Information about Effect2),etc.
And here is the precise layout with all the subfields:
It's also possible to get SnpEff to output in a (non-VCF) text format with one Effect per line. See the SnpEff home page for full details.
Once you have a SnpEff output VCF file, you can use the VariantAnnotator walker to add SnpEff annotations based on that output to the input file you ran SnpEff on.
There are two different options for doing this:
NOTE: This option works only with supported SnpEff versions as explained above. VariantAnnotator run as described below will refuse to parse SnpEff output files produced by other versions of the tool, or which lack a SnpEff version number in their header.
The default behavior when you run VariantAnnotator on a SnpEff output file is to parse the complete set of effects resulting from the current variant, select the most biologically-significant effect, and add annotations for just that effect to the INFO field of the VCF record for the current variant. This is the mode we plan to use in our Production Data-Processing Pipeline.
When selecting the most biologically-significant effect associated with the current variant, VariantAnnotator does the following:
Prioritizes the effects according to the categories (in order of decreasing precedence) "High-Impact", "Moderate-Impact", "Low-Impact", and "Modifier", and always selects one of the effects from the highest-priority category. For example, if there are three moderate-impact effects and two high-impact effects resulting from the current variant, the annotator will choose one of the high-impact effects and add annotations based on it. See below for a full list of the effects arranged by category.
Within each category, ties are broken using the functional class of each effect (in order of precedence: NONSENSE, MISSENSE, SILENT, or NONE). For example, if there is both a NON_SYNONYMOUS_CODING (MODERATE-impact, MISSENSE) and a CODON_CHANGE (MODERATE-impact, NONE) effect associated with the current variant, the annotator will select the NON_SYNONYMOUS_CODING effect. This is to allow for more accurate counts of the total number of sites with NONSENSE/MISSENSE/SILENT mutations. See below for a description of the functional classes SnpEff associates with the various effects.
Effects that are within a non-coding region are always considered lower-impact than effects that are within a coding region.
java -jar dist/GenomeAnalysisTK.jar \ -T VariantAnnotator \ -R /humgen/1kg/reference/human_g1k_v37.fasta \ -A SnpEff \ --variant 1000G.exomes.vcf \ (file to annotate) --snpEffFile snpEff_output.vcf \ (SnpEff VCF output file generated by running SnpEff on the file to annotate) -L 1000G.exomes.vcf \ -o out.vcf
VariantAnnotator adds some or all of the following INFO field annotations to each variant record:
SNPEFF_EFFECT- The highest-impact effect resulting from the current variant (or one of the highest-impact effects, if there is a tie)
SNPEFF_IMPACT- Impact of the highest-impact effect resulting from the current variant (
SNPEFF_FUNCTIONAL_CLASS- Functional class of the highest-impact effect resulting from the current variant (
SNPEFF_CODON_CHANGE- Old/New codon for the highest-impact effect resulting from the current variant
SNPEFF_AMINO_ACID_CHANGE- Old/New amino acid for the highest-impact effect resulting from the current variant
SNPEFF_GENE_NAME- Gene name for the highest-impact effect resulting from the current variant
SNPEFF_GENE_BIOTYPE- Gene biotype for the highest-impact effect resulting from the current variant
SNPEFF_TRANSCRIPT_ID- Transcript ID for the highest-impact effect resulting from the current variant
SNPEFF_EXON_ID- Exon ID for the highest-impact effect resulting from the current variant
Example VCF records annotated using SnpEff and VariantAnnotator:
1 874779 . C T 279.94 . AC=1;AF=0.0032;AN=310;BaseQRankSum=-1.800;DP=3371;Dels=0.00;FS=0.000;HRun=0;HaplotypeScore=1.4493;InbreedingCoeff=-0.0045; MQ=54.49;MQ0=10;MQRankSum=0.982;QD=13.33;ReadPosRankSum=-0.060;SB=-120.09;SNPEFF_AMINO_ACID_CHANGE=G215;SNPEFF_CODON_CHANGE=ggC/ggT; SNPEFF_EFFECT=SYNONYMOUS_CODING;SNPEFF_EXON_ID=exon_1_874655_874840;SNPEFF_FUNCTIONAL_CLASS=SILENT;SNPEFF_GENE_BIOTYPE=protein_coding;SNPEFF_GENE_NAME=SAMD11; SNPEFF_IMPACT=LOW;SNPEFF_TRANSCRIPT_ID=ENST00000342066 1 874816 . C CT 2527.52 . AC=15;AF=0.0484;AN=310;BaseQRankSum=-11.876;DP=4718;FS=48.575;HRun=1;HaplotypeScore=91.9147;InbreedingCoeff=-0.0520; MQ=53.37;MQ0=6;MQRankSum=-1.388;QD=5.92;ReadPosRankSum=-1.932;SB=-741.06;SNPEFF_EFFECT=FRAME_SHIFT;SNPEFF_EXON_ID=exon_1_874655_874840; SNPEFF_FUNCTIONAL_CLASS=NONE;SNPEFF_GENE_BIOTYPE=protein_coding;SNPEFF_GENE_NAME=SAMD11;SNPEFF_IMPACT=HIGH;SNPEFF_TRANSCRIPT_ID=ENST00000342066
VariantAnnotator also has the ability to take the EFF field from the SnpEff VCF output file containing all the effects aggregated together and copy it verbatim into the VCF to annotate.
Here's an example of how to do this:
java -jar dist/GenomeAnalysisTK.jar \ -T VariantAnnotator \ -R /humgen/1kg/reference/human_g1k_v37.fasta \ -E resource.EFF \ --variant 1000G.exomes.vcf \ (file to annotate) --resource snpEff_output.vcf \ (SnpEff VCF output file generated by running SnpEff on the file to annotate) -L 1000G.exomes.vcf \ -o out.vcf
Of course, in this case you can also use the VCF output by SnpEff directly, but if you are using VariantAnnotator for other purposes anyway the above might be useful.
Below are the possible genomic effects recognized by SnpEff, grouped by biological impact. Full descriptions of each effect are available on this page.
SnpEff assigns a functional class to certain effects, in addition to an impact:
NONSENSE: assigned to point mutations that result in the creation of a new stop codon
MISSENSE: assigned to point mutations that result in an amino acid change, but not a new stop codon
SILENT: assigned to point mutations that result in a codon change, but not an amino acid change or new stop codon
NONE: assigned to all effects that don't fall into any of the above categories (including all events larger than a point mutation)
The GATK prioritizes effects with functional classes over effects of equal impact that lack a functional class when selecting the most significant effect in VariantAnnotator. This is to enable accurate counts of NONSENSE/MISSENSE/SILENT sites.
I am planning to use snpEff to map my variant data (SNPs & Indels) to genes data.
My query is - Can snpEff handle vcfs with both Indels as well as SNPs? (the snpEff site does indicate that it is possible, just want to confirm the same.) - Secondly, all the mappings that are done, are they done per SNP basis or per codon basis? For instance if you have two SNPs next to each other, their combined effects may change the codon (and possibly the amino acid it is coding for). Does the tool do this?
Thanks for your input!
Hi, We use GATK with snpEff to annotate variants of our WES experiments. We use VariantAnnotator to keep only the effect with highest impact. We've realized that for some experiments it could be useful for us to customize priorities assignment (i.e. if a variant is categorized both as DOWNSTREAM for a gene and INTRON for the adjacent gene, the effect picked by default is DOWNSTREAM while we'd like to keep INTRON, etc.) Is there a way to modify VariantAnnotator default behaviour?
I have the genomes of several isolates of a parasite, and I would like to investigate synonymous/non-synonymous substitution for identifying potential antigens, as well as SNPs genome-wide and I am wondering how well BWA/GATK are suited for this purpose. I've been told that BWA is only very good with sequences <2% divergent, and some of the antigens in this specie are known to be >20% divergent. However, I also know that GATK does local realignments of indels. So I would specifically like to know - is BWA/GATK good for looking at substitutions/SNPs in highly variable genes, and if not which other alignment tools are compatible and appropriate for this purpose?
Hi GATK Team,
I have been experimenting with newer versions of snpEff and would like to incorporate a newer version into my local GATK build. However, when I update the jar to the snpeff 3.0 release (by updating the ivy specification and copying the jar in into the proper settings/repository/... directory) I get an internal scalar compiler error:
scala.compile.public: [mkdir] Created dir: /hpc/users/lindem03/packages/gatk-mssm/dev/build/scala/classes [echo] Building Scala... [scalac] Compiling 87 source files to /hpc/users/lindem03/packages/gatk-mssm/dev/build/scala/classes [scalac] Compiling 107 source files to /hpc/users/lindem03/packages/gatk-mssm/dev/build/scala/classes [scalac] Exception in thread "main" java.lang.AssertionError: assertion failed: List(object Byte, object Byte) [scalac] at scla.tools.nsc.symtab.Symbols$Symbol.suchThat(Symbols.scala:1056) [scalac] at scala.tools.nsc.symtab.Symbols$Symbol.companionModule0(Symbols.scala:1271) [scalac] at scala.tools.nsc.symtab.Symbols$Symbol.companionModule(Symbols.scala:1281) [scalac] at scala.tools.nsc.symtab.Symbols$Symbol.linkedClassOfClass(Symbols.scala:1302) [scalac] at scala.tools.nsc.symtab.Definitions$definitions$.addModuleMethod$1(Definitions.scala:711) [scalac] at scala.tools.nsc.symtab.Definitions$definitions$.initValueClasses(Definitions.scala:714) [scalac] at scala.tools.nsc.symtab.Definitions$definitions$.init(Definitions.scala:791) [scalac] at scala.tools.nsc.Global$Run.<init>(Global.scala:604) [scalac] at scala.tools.nsc.Main$.process(Main.scala:105) [scalac] at scala.tools.nsc.Main$.main(Main.scala:120) [scalac] at scala.tools.nsc.Main.main(Main.scala)
Any thoughts? Do I need to modify the snpEff jar in some way to get it play nicely with GATK? That error is so opaque I am not sure where to start debugging. This is against a "clean" directory (I ran
ant clean) with java 1.6.0_30
$ java -version java version "1.6.0_30" Java(TM) SE Runtime Environment (build 1.6.0_30-b12) Java HotSpot(TM) 64-Bit Server VM (build 20.5-b03, mixed mode)
and GATK 1.6-24-gdc14575.
Broad recommends using snpEff to add annotations to VCF files created by GATK. This gives annotations about the effect of a given variant: is it in a coding region? Does it cause a frameshift? What transcripts are impacted? etc. However, snpEff does not provide other annotations you might want, such as 1000 genomes minor allele frequency, SIFT scores, phyloP conservation scores, and so on. I've previously used annovar to get those sorts of things, and that worked well enough, though I did not find it to be especially user-friendly.
So my question is, what other ways have users found of getting this sort of annotation information? I'm interested specifically in human exomes, but I am sure other users reading this Ask the Community post will be interested in answers for other organisms as well. I'm looking for recommendations on what's quick, simple, easy to use, and has been used successfully with VCFs produced by GATK. I'm open to answers in the form of other software tools or sources of raw data that I can easily manipulate on my own.
Thanks in advance.
Hi there, I've done with a run of HaplotypeCaller on my samples. I'm now analysing everything with snpEff, although I'm doing this "outside" GATK. I had to stop the analysis because a huge number of errors, all dealing with indels, such as:
Error while processing VCF entry (line 580649) : chr21 26718345 . TAATCCTGAGTTTAA TATCCTAAATGTTTAC 943.26 […] java.lang.RuntimeException: Insertion '-A+AT' does not start with '+'. This should never happen! chr21 35260360 . CATAACAGTTCAT AGAGACAGAG 425.22 […] java.lang.RuntimeException: Deletion '+G-TTC' does not start with '-'. This should never happen!
Of course, this is a snpEff error, nevertheless the Indel format looks quite different from what I've ever seen. Consider the first line above: shouldn't it be like
chr21 26718345 . AT T 943.26 […]
(I can't resolve the second right now). Any hint is appreciated at this point. I'm writing to snpEff developer for the same reason...