Variant annotation tutorial
Code and tutorials for CIHEAM course. Variation annotation and function prediction sessions
Variant annotation
Slicing and dicing of VCF files
VCF files can get very big. Being able to manipulate them and extract relevant information is important. In practice, flat text files such as VCF are currently the basis for further analysis - see the Imputation and Population practical sessions.
The first task involves selecting all called variants of chromosome 18, bgzip-ping them on-the-fly, and then indexing that file with tabix:
tabix -h allbt.vcf.gz 18 | bgzip >BT18.vcf.gz tabix -p vcf BT18.vcf.gz
Annotating VCF with rs-numbers
The VCF file you are working does not yet have the dbSNP identifiers (rs-numbers) annotated. The rs-numbers of all currently known variants as databased in dbSNP can be found in this file:
/path/to/ BT18.vcf.gz
Have a look at what is inside this file:
gunzip -c BT18.vcf.gz | grep -v '^##' | more
or
tabix BT18.vcf.gz 18:1-10000 | more
You can use vcftools to do many things related to VCF files. This includes annotating the VCF based on other interval-file formatted input, including other VCF files. For instance, you can annotate an interval of the VCF on the fly, by selecting that interval using tabix, and then piping it to vcf-annotate, like so:
tabix -h BT18.vcf.gz 18:100000-101000 | vcf-annotate -a BT_incl_cons.18.vcf.gz -c CHROM,FROM,ID
Now make a new VCF file for all the variants called for chromosome 18 using this tool:
tabix -h BT18.vcf.gz 18 | vcf-annotate -a BT_incl_cons.18.vcf.gz -c CHROM,FROM,ID | bgzip >BT18_rsnumbers.vcf.gz tabix -p vcf BT18_rsnumbers.vcf.gz
Extracting variants using BEDtools
You will very often want to extract variants based on interval files, e.g. a certain region on a chromosome. This can be easily achieve using the tabix command described above. However, sometimes you may want to extract a substantial number of intervals. You could use tabix in a loop, either as a shell-script or in a scripting language like Perl or Python. However, BedTools provides many utilities for manipulating interval data, such as BED, VCF, or GTF/GFF formats.
Say you are interested in the variants in all the genes lying between coordinates 1,000,000 and 2,000,000 on chromosome 18. The annotation of genes can be found in a GTF file. The GTF file for the UMD3.1 genome build can be found at Ensembl and downloaded using the command line (NOTE: the compute nodes you are working on do not have internet connectivity. This line of code is just for completeness):
wget ftp://ftp.ensembl.org/pub/release-78/gtf/bos_taurus/Bos_taurus.UMD3.1.78.gtf.gz
The file is already downloaded and can be found here:
/path/to/Bos_taurus.UMD3.1.78.gtf.gz
Have a quick look at what is in this file:
gunzip -c /path/to/Bos_taurus.UMD3.1.78.gtf.gz | more
Familiarise yourself with the GTF format if you have not already. For more information:
http://www.ensembl.org/info/website/upload/gff.html
You can select the appropriate region from the GTF file:
gunzip -c Bos_taurus.UMD3.1.78.gtf.gz | awk '$3=="gene"' | awk '$1==18&&$4>1000000&&$4<2000000' >mygenes.gtf
Subsequently you can intersect the VCF with the GTF file that contains just the genes between coordinates 1 and 2 million:
bedtools intersect -a BT18_rsnumbers.vcf.gz -b mygenes.gtf >intersection_genes.vcf
Have a quick look at the output, and count the number of variants in this VCF (e.g. by using 'wc -l').
Now do the same, but select only the exonic variants in the same interval, from 1 to 2 million.
How many variants do you have now?
Variant Effect Predictor
Various tools exist to make more precise annotations than what you can achieve using Bedtools and a GTF. Among them are Variant Effect Predictor (VEP), snpEff, and Annovar. For the current practical we will use former two.
VEP can be found here: perl /home/formacion/COMUNES/IAMZ/soft/ensembl-tools-release-78/scripts/variant_effect_predictor/variant_effect_predictor.pl
If you invoke it without parameters you will see the following:
Usage:
perl variant_effect_predictor.pl [--cache|--offline|--database] [arguments]
Basic options
=============
--help Display this message and quit
-i | --input_file Input file
-o | --output_file Output file
--force_overwrite Force overwriting of output file
--species [species] Species to use [default: "human"]
--everything Shortcut switch to turn on commonly used options. See web
documentation for details [default: off]
--fork [num_forks] Use forking to improve script runtime
For full option documentation see:
http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html
For the current session we will use the following options:
--species bos_taurus -o BT18.vep # outputfile --fork 4 # use 4 threads --canonical --sift b --coding_only --no_intergenic --offline # using a cache directory --dir /home/formacion/COMUNES/IAMZ/data/CIHEAM/ReferenceGenome/VEP/ # where the cache dir lives --force_overwrite
You can create a file that will provide annotation information for each SNP:
gunzip -c BT18_rsnumbers.vcf.gz | \ perl /home/formacion/COMUNES/IAMZ/soft/ensembl-tools-release-78/scripts/variant_effect_predictor/variant_effect_predictor.pl \ --dir /home/formacion/COMUNES/IAMZ/data/CIHEAM/ReferenceGenome/VEP/ --species bos_taurus -o BT18.vep \ --fork 4 --canonical --sift b --coding_only --no_intergenic --offline --force_overwrite
Have a look inside the output file and study its content. Why are there no SNPs in this list that are intergenic?
The VEP output you generated before has many advantages, one of which is that it is nicely tabular and easy to parse. In fact a demonstration of that will be covered in the 'Provean' practical. For some purposes, however, it might be more convenient to have the annotations directly in the VCF. Furthermore, you may want annotations for each and every SNP, even if e.g. they are intergenic. You could for instance do the following, omitting the '--no_intergenic and --coding_only flags, and in stead adding the --vcf flag:
gunzip -c BT18_rsnumbers.vcf.gz | \ perl /home/formacion/COMUNES/IAMZ/soft/ensembl-tools-release-78/scripts/variant_effect_predictor/variant_effect_predictor.pl \ --dir /home/formacion/COMUNES/IAMZ/data/CIHEAM/ReferenceGenome/VEP/ --species bos_taurus -o BT18.vep.vcf \ --fork 4 --canonical --sift b --offline --force_overwrite --vcf bgzip BT18.vep.vcf tabix -p vcf BT18.vep.vcf.gz
Again, study and interpret your outcome by looking inside the file. Specifically, note WHERE the annotation ended up.
18 53205918 . G A 0.329662 . AB=0;ABP=0;AC=0;AF=0.125;AN=16;AO=2;CIGAR=1X;DP=70;DPB =70;DPRA=1.03279;EPP=7.35324;EPPR=3.13803;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=60;NS=8;NUMALT=1; ODDS=2.54006;PAIRED=1;PAIREDR=0.985294;PAO=0;PQA=0;PQR=0;PRO=0;QA=70;QR=2421;RO=68;RPP=3.0103; RPPR=11.1853;RUN=1;SAF=1;SAP=3.0103;SAR=1;SRF=49;SRP=31.7504;SRR=19;TYPE=snp;technology.ILLUMINA=1; CSQ=A|ENSBTAG00000018834|ENSBTAT00000025071|Transcript|missense_variant|2132|1984|662|E/K|Gag/Aag|||1| YES|deleterious(0.02) GT:DP:RO:QR:AO:QA:GL 0/0:7:7:263:0:0:0,-1.92588,-10 0/0:7:7:238:0:0:0,-1.92376,-10 0/0:9:7:266:2:70:-3.90732,0,-10 0/0:10:10:362:0:0:0,-2.73805,-10 0/0:10:10:345:0:0:0,-2.73789,-10 0/0:8:8:287:0:0:0,-2.19662,-10 0/0:5:5:178:0:0:0,-1.38409,-10 0/0:14:14:482:0:0:0,-3.82032,-10
A further look at the VEP output
Although the VCF now looks even more complex than before, there is a good reason to add the VEP annotation data the way it was done. The 8th column of the VCF provides a lot of information, each field separated by a ';'. In addition, each field contains a field name and a value, delimited by a '='. The VEP annotation has its own field name, or tag: 'CSQ'. Within this field, the values are again structured, but this time delimited by a '|'. What we in fact have now is a nested data structure:
['18', '53205918', '.', 'G', 'A', '0.329662', '.', {'PAIREDR': '0.985294', 'DPB': '70', 'MQMR': '60', 'PQA': '0', 'PAO': '0',
'SAR': '1', 'AB': '0', 'GTI': '0', 'NUMALT': '1', 'RO': '68', 'DP': '70', 'DPRA': '1.03279',
'CSQ': ['A', 'ENSBTAG00000018834', 'ENSBTAT00000025071', 'Transcript', 'missense_variant',
'2132', '1984', '662', 'E/K', 'Gag/Aag', , , '1', 'YES', ['deleterious', '0.02']], 'EPP': '7.35324',
'SRP': '31.7504', 'QR': '2421', 'SAF': '1', 'ODDS': '2.54006', 'SRR': '19', 'SRF': '49',
'technology.ILLUMINA': '1', 'QA': '70', 'EPPR': '3.13803', 'SAP': '3.0103', 'RPPR': '11.1853', 'PQR': '0',
'NS': '8', 'CIGAR': '1X', 'TYPE': 'snp', 'AF': '0.125', 'ABP': '0', 'LEN': '1', 'AC': '0', 'MEANALT': '1',
'PAIRED': '1', 'AO': '2', 'AN': '16', 'MQM': '60', 'RUN': '1', 'PRO': '0', 'RPP': '3.0103'}, 'GT:DP:RO:QR:AO:QA:GL',
'0/0:7:7:263:0:0:0,-1.92588,-10', '0/0:7:7:238:0:0:0,-1.92376,-10', '0/0:9:7:266:2:70:-3.90732,0,-10',
'0/0:10:10:362:0:0:0,-2.73805,-10', '0/0:10:10:345:0:0:0,-2.73789,-10', '0/0:8:8:287:0:0:0,-2.19662,-10',
'0/0:5:5:178:0:0:0,-1.38409,-10', '0/0:14:14:482:0:0:0,-3.82032,-10']
When dealing with deleterious alleles, we might be interested to learn how frequent that allele is in the population we are studying. Our expectation is that genuinly deleterious alleles should be relatively rare in a population. If not they might hint at either a sign of inbreeding, or, perhaps, something that might be deleterious in the wild but advantageous in the context of domestication/breeding. The VCF file holds the (reference) allele frequency information. The 8th column of the VCF file in fact is a bunch of fields concatenated with the ';' delimitor. One of the fields has the tag 'AF=', where AF stands for (Alternative) Allele Frequency. We can expand the number of tab-delimited fields easily by replacing ';' for '\t', and then look for the 11th column. We need to remove the 'AF=' tag though, but if we do that we can further select for instance for allele frequencies between 0 and 15%.
gunzip -c BT18.vep.vcf.gz | grep deleterious | sed 's/;/\t/g' | awk '{gsub("AF=","",$11); print}' | awk '$11<0.15&&$11>0' | wc -l
Analysing the annotations a bit further, you could discover that there are quite a few variants that have high deleterious allele frequencies. Not only do some putative deleterious alleles occur at high frequency, but there is also a number of genes that has quite a few deleterious alleles annotated in the current population sample. Say you would like to make a table that contains all genes that have more than one putative deleterious variant. We would need to isolate both the AF field and the CSQ field, and from the latter isolate the 'gene' field. While this is possible with a fairly simple shell scripting hack, it does require an awful lot of counting to make sure you end up with the right fields. Ideally you would like to have more control by putting the relevant data in a complex data structure. Languages such as Perl and particularly Python allow far more explicit syntax to achieve that, so we will swap out some complex string of shell-commands by a single Python script. That script does only one thing: takes lines from the VCF file, put all fields in a complex data structure (as seen above) and, in a structured way, prints a selection of fields.
What this script will VCF input from a stream:
gunzip -c LeonEnv/BT18.vep.all.vcf.gz | python3 print_fields_from_vcf.py -m vep
Notice that the script itself does in fact take an option: 'vep'. This is done because it can also work for the snpEff annotation, which, highly inconveniently, is slightly different.
and for each line provide output similar to this:
# chrom coord AF Gene consequence sift_prediction sift_score 18 53205918 0.125 ENSBTAG00000018834 missense_variant deleterious 0.02
Obviously, the functionality of the Python script is limited. However, it would be quite easy to provide additional statistics and output options based on the current script, because it makes the data structure so explicit. This means that the Python script provides you with a means for expansion of functionality, while that will be far more tricky using shell commands solely.
#alternative way, shell hack
gunzip -c BT18.vep.vcf.gz | grep deleterious | sed 's/;/ /g' | awk '{gsub("AF=","",$11); print}' | awk '$11<1&&$11>0' | sed 's/ /#_#/g' | sed 's/|/\t/g' | cut -f2 | sort | uniq -c | sed 's/^ \+//' | awk '$1>1' >genes_with_more_than1.txt
Now we make a file that has the gene names, and counts, for all genes that have more than one deleterious allele:
gunzip -c LeonEnv/BT18.vep.all.vcf.gz | python3 print_fields_from_vcf.py -m vep | grep deleterious | awk '$3>0&&$3<1' | cut -f4 | sort | uniq -c | sed 's/^ \+//' | awk '$1>1' >genes_with_more_than1.txt
Take two or three genes that have the highest number of deleterious variants. Manually look them up at the Ensembl website (ensembl.org). Check the 'orthologous' status. What do you observe? How would you interpret the validity of the gene annotation and/or the likelihood of off-site mapping of reads?
This provides us with a strategy to filter a bit more rigorously for genes that may actually contain genuine deleterious alleles. Ideally we would a priory filter for genes that have proper one-to-one orthologous relations with other mammalian species. Time and lack of direct internet connection makes that unfeasible for this practical. However, genes that have only one deleterious variant segregating in the study population at a low frequency would certainly be better candidates than genes that have multiple deleterious variants, as the above example indicates. To filter you can do the following:
gunzip -c LeonEnv/BT18.vep.all.vcf.gz | python3 print_fields_from_vcf.py -m vep | grep deleterious | awk '$3>0&&$3<1' | cut -f4 | sort | uniq -c | sed 's/^ \+//' | awk '$1==1' | awk '{print $2}' >genes_with1.txt
How many genes are there that meet this condition?
snpEff
Another popular variant annotation tool is 'snpEff'. Among the advantages over VEP is that it is applicable to a wider range of species. Furthermore, it is blazingly fast. For chromosome 18, we will also use this tool:
snpEff is Java-based. The path to the snpEff.jar Java archive is available as an environment variable. Do:
echo $SNP /home/formacion/COMUNES/IAMZ/soft/snpEff/snpEff.jar
The snpEff program can be used like this:
java -Xmx4G -jar $snpEff -dataDir ~/snpEff/data -v UMD3.1.78 BT18_rsnumbers.vcf.gz >testresult_snpEff_Bt.txt
The options used here are:
-dataDir ~/snpEff/data # annotation data goes here -v UMD3.1.78 # Genome / version BT18_rsnumbers.vcf.gz # requires an input file, can be gzipped >testresult_snpEff_Bt.txt # results go to STDOUT and are captured in a file by '>'
Study the outcome. Try to filter out 'intergenic' SNPs and other things that are perhaps a bit less relevant (use your skills in 'grep-ping', etc.). What are some of the differences in how the variants are annotated, compared to VEP?
You will notice that the variants that have an effect are labelled 'MODERATE' or 'HIGH'. Try to filter out the ones that are 'HIGH'. How many are there?
You can now easily make tables for number of different consequences as annotated by both these annotation tools:
cat LeonEnv/testresult_snpEff_Bt.txt | python3 cut_line.py -m snpeff | cut -f6 | sort | uniq -c
gunzip -c LeonEnv/BT18.vep.all.vcf.gz | python3 cut_line.py -m vep | cut -f5 | sort | uniq -c
Bonus question: Compare the results by selecting the 'MODERATE' and "HIGH" variants, and compare that with the 'deleterious' variants as annotated (through SIFT) in VEP. You can use, e.g. vcf-compare.
Inferring function
Provean
Provean uses an approach where a protein sequence is compared (Psiblast) to all known protein sequences (so called 'non-redundant database'). Subsequent clustering defines a group of 'similar', and implicitly homologous sequences. An example of in- and output of a pig sequence:
What goes in: Protein sequence (Fasta):
>ENSSSCP00000018263 pep:novel chromosome:Sscrofa10.2:12:6621092:6624938:1 gene:ENSSSCG00000017236 transcript:ENSSSCT00000018765 gene_biotype:protein_coding t transcript_biotype:protein_coding MTPRVGAVWLPSALLLLRVPGCLSLSGPPTAMGTKGGSLSVQCRYEEEYIDDKKYWDKSP CFLSWKHIVETTESAREVRRGRVSIRDDPANLTFTVTLERLTEEDAGTYCCGITAQFSVD PTHEVEVVVFPALGTSRPPSMPGPPTTLPATTWSFVSERETMANNLGKGPASQDPGQHPR SKHPSIRLLLLVFLEVPLFLGMLGAVLWVHRPLRSSESRSVAMDPVPGNTAPSAGWK
and variation, with coordinates pertaining the above protein sequences, the second column the reference amino acid, and the third/last column the alternative amino acid (alternative allele):
235,G,E 5,V,A 22,C,Y 34,T,I 51,D,N 53,K,N 59,S,Y 61,C,R 64,S,L 67,H,P 68,I,T 75,A,V 108,T,K 115,A,T 124,E,D 130,F,Y 133,L,P 142,P,A 155,F,I 158,E,G 186,I,V 21,G,D
Output looks like this:
## PROVEAN scores ## # VARIATION SCORE 235,G,E 0.076 5,V,A 0.287 22,C,Y -8.028 34,T,I -1.932 51,D,N -1.613 53,K,N 1.565 59,S,Y -2.140 61,C,R -5.826 64,S,L -0.437 67,H,P -0.511 68,I,T -2.664 75,A,V -1.061 108,T,K -4.051 115,A,T 1.557 124,E,D -1.587 130,F,Y -0.983 133,L,P 3.203 142,P,A -2.058 155,F,I 0.502 158,E,G -0.077 186,I,V 0.220 21,G,D -6.014
We will do something similar, although with a somewhat less complex example, for the cow data. At the end of the VEP exercise we created a filtered list of genes that have only a single deleterious variant annotated. Among these genes was BLOC1S3 gene, ENSBTAG00000007070.
http://www.ensembl.org/Bos_taurus/Gene/Summary?db=core;g=ENSBTAG00000007070;r=18:53231599-53232201;t=ENSBTAT00000009293
What we need for the Provean analysis is 1) the protein sequence of the gene (or rather, the transcript), and 2) a file with the variation in the same format as listed above. Obviously, if we want to do this for many genes, we will not parse the relevant information by hand. A multifasta file that contains all protein sequences for cow can be found at Ensembl. It is located here on the machine you are working on:
GENE=ENSBTAG00000007070
PROT=`cat ~/Bos_taurus.UMD3.1.pep.all.fa | grep $GENE | awk '{print $1}' | sed 's/>//'`
cat test.vep | grep $GENE | grep missense | awk '{print $10","$11}' | sed 's/\//,/' >$PROT.var
provean.sh -q $PROT.fa -v $PROT.var --save_supporting_set $PROT.sss >$PROT.result.txt 2>$PROT.error;
faOneRecord ~/Bos_taurus.UMD3.1.pep.all.fa $PROT >$PROT.fa
cat ENSBTAG00000007070.fa >>ENSBTAG00000007070.sss.fasta mafft --auto ENSBTAG00000007070.sss.fasta >ENSBTAG00000007070.sss.fasta.out
Transfer the alignment to your local environment and view (seaview?). What do you observe for the 26th amino-acid of the cow protein sequence? In this 'implicit evolutionary context' there is one alternative amino acid that occurs at this site. What do the D amino acid and this evolutionary viable alternative share in terms of properties?
Through the orthology with human, try to find some possible phenotypic consequences of the alternative allele (e.g. through OMIM). Could there be a reason that the variant is at a relatively high frequency in this cattle population?
Polyphen
/lustre/nobackup/WUR/ABGC/shared/public_data_store/polyphen/polyphen-2.2.2/bin/run_pph.pl -s test1033.fa test1033.coord #o_acc o_pos o_aa1 o_aa2 rsid acc pos aa1 aa2 nt1 nt2 prediction based_on effect site region PHAT dScore Score1 Score2 MSAv Nobs Nstruct Nfilt PDB_id PDB_pos PDB_ch ident lengthNormASA SecStr MapReg dVol dProp B-fact H-bonds AveNHet MinDHet AveNInt MinDInt AveNSit MinDSit Transv CodPos CpG MinDJxn PfamHit IdPmax IdPSNP IdQmin ENSSSCG00000001099ENSSSCT00000001195 368 S A ENSSSCG00000001099ENSSSCT00000001195 368 S A benign alignment +0.721 -1.180 -1.901 2 68 12.350 12.350 42.63 ENSSSCG00000001099ENSSSCT00000001195 453 R H ENSSSCG00000001099ENSSSCT00000001195 453 R H benign alignment +0.351 -2.130 -2.481 2 67 33.194 33.194 93.76 ENSSSCG00000001099ENSSSCT00000001195 431 K T ENSSSCG00000001099ENSSSCT00000001195 431 K T possibly damaging alignment +1.705 -1.678 -3.383 2 68 2.382 68.80
What goes in:
test1033.fa >ENSSSCG00000001099ENSSSCT00000001195 MSSIEQTTEILLCLSPAEAANLKEGINFVRNKSTGKDYILFKNKSRLKACKNMCKHQGGL FIKDIEDLNGRSVKCTKHNWKLDVSSMKYINPPGSFCQDELVVEKDEENGVLLLELNPPN PWDSEPRSPEDLAFGEVQITYLTHACMDLKLGDKRMVFDPWLIGPAFARGWWLLHEPPSD WLERLSRADLIYISHMHSDHLSYPTLKKLAERRPDVPIYVGNTERPVFWNLNQSGVQLTN INVVPFGIWQQVDKNLRFMILMDGVHPEMDTCIIVEYKGHKILNTVDCTRPNGGRLPMKV ALMMSDFAGGASGFPMTFSGGKFTEEWKAQFIKTERKKLLNYKARLVKDLQPRIYCPFAG YFVESHPSDKYIKETNIKNDPNELNNLIKKNSEVVTWTPRPGATLDLGRMLKDPTDSKGI VEPPEGTKIYKDSWDFGPYLNILNAAIGDEIFRHSSWIKEYFTWAGFKDYNLVVRMIETD EDFSPLPGGYDYLVDFLDLSFPKERPSREHPYEEIRSRVDVIRHVVKNGLLWDDLYIGFQ TRLQRDPDIYHHLFWNHFQIKLPLTPPDWKSFLMCSG
test1033.coord ENSSSCG00000001099ENSSSCT00000001195 368 S A ENSSSCG00000001099ENSSSCT00000001195 453 R H ENSSSCG00000001099ENSSSCT00000001195 431 K T