Annotating Variants: Introduction

As we've already seen, determining the presence or absence of a variant from NGS data is not trivial. It is software dependent and has inherent trade-offs between sensitivity and specificity. Inevitably, the number of putative variants in a real data set is very large; for example the first samples from the 1000 Genomes project typically found 0.1% variants (3 million variants), approximately 10% of which had never been previously observed (300,000 novel variants per individual.) False positive discovery rates are also typically very high at this stage.

...

Because these tools draw in information from may disparate sources, they can be very difficult to install, configure, use, and maintain. For example, the vcf files from the 1000 Genomes project are arranged in a deep ftp tree by date of data generation. Large genome centers spend significant resources managing these tools.

Pre-packaged programs

Annovar - one of the most powerful yet simple to run variant annotators available

Annovar is a variant annotator. Given a vcf file from an unknown sample and a host of existing data about genes, other known SNPs, gene variants, etc., Annovar will place the discovered variants in context.

...

This next exercise will give you some idea of how Annovar works; we've taken the liberty of writing the bash script annovar_pipe.sh around the existing summarize_annovar.pl wrapper (a wrapper within a wrapper - a common trick) to even further simplify the process for this course.

Exercise:

First, look at the code for our annovar_pipe.sh command.

...

cp $BI/ngs_course/human_variation/*chrom20.GATK* .

ANNOVAR output

Annovar does a ton of work in assessing variants for us (though if you were going for clinical interpretation, you still have a long way to go - compare this to RUNES or CarpeNovo).  It provides all these output files:

...

Everything after the "LJB_GERP++" field in exome_summary came from the original VCF file, so this file REALLY contains everything you need to go on to functional analysis!  This is one of the many reasons I like Annovar.

Scavenger hunts!

grep "frameshift" NA12878.chrom20.GATK.vcf.exome_summary.csv  # this will print all the lines which contains the text "frameshift"
# From the output you can key into the first few columns having the information you are interested in: location-classification, gene, mutation type each separated by commas. This should lead you to think about adding the awk command to print only some columns.
 
grep "frameshift" NA12878.chrom20.GATK.vcf.exome_summary.csv | awk -F"," '{print $2"\t"$3}'  # the -F"," syntax forces it to split on commas
# you will likely notice this data is easier to visualize, and in this case you can probably see what gene is represented multiple times, but why stop there ... lets add the uniq -c command to the pipes to have linux count for us
 
grep "frameshift" NA12878.chrom20.GATK.vcf.exome_summary.csv | awk -F"," '{print $2"\t"$3}' | uniq -c 
# for the number of mutations we have this is sufficient, but for increased numbers of mutations where you may be interested in displaying them in a particular order. This can be done by adding the sort command

grep "frameshift" NA12878.chrom20.GATK.vcf.exome_summary.csv | awk -F"," '{print $2"\t"$3}' | uniq -c | sort -r  # the -r option on the sort command  sorts in reverse order

...

grep "frameshift" NA12878.chrom20.GATK.vcf.exome_summary.csv | awk -F"," '{print $2"\t"$3}' | uniq -c | sort -r

cat *GATK.vcf.exome_summary.csv | grep "exonic" | awk -F"," '{print $2,"\t"$3}' | sort | uniq -c | sort -nr

Other variant annotators:

• http://www.yandell-lab.org/software/vaast.html
• http://www.broadinstitute.org/gatk/gatkdocs/ VariantAnnotator annotations
• http://www.bioconductor.org/help/workflows/variants/
• http://vat.gersteinlab.org/
• http://code.google.com/p/mu2a/

Notes

Variants consist of single base base changes, insertions and deletions, and larger scale structural changes. "Larger scale" is usually defined relative to the capabilities of the technology; for example, a "small indel" usually means "detectable within a single sequence read". In 2009, sequence reads were about 50 bp but in 2011 they were 100 bp.