Beginner’s Guide to Genomic Data Analysis: Variant Annotation | SnpEff Tutorial
This is a continuation of the series of articles focused on computational biology for beginners. This is a beginner-level introduction to knowing the effects of variants. For example, checking if our mutations in a specific gene will make a different amino acid or not. Knowing such a small change can help us identify/suspect diseases or any other trait for which we are looking. Additionally, this tutorial also includes a way to handle the CHROMOSOME NOT FOUND ERROR while using SnpEff.
Setting up SnpEff
Go to the official website of SnpEff and download either using the Download button or alternatively copy the link address and download in the terminal.
wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip
unzip snpEff_latest_core.zip
rm snpEff_latest_core.zipAfter using the above commands, you will have SnpEff and SnpSift both setup and ready to use.
Usage
There are various examples on the official website of SnpEff for various use cases. Additionally, there is also complete documentation online. I will be using SnpEff to work on the VCF files generated in my previous article.
First of all, check the installation/setup and see the usage and available commands.
java -jar snpEff/snpEff.jar
# Output
SnpEff version SnpEff 5.2c (build 2024-04-09 12:24), by Pablo Cingolani
Usage: snpEff [command] [options] [files]Run 'java -jar snpEff.jar command' for help on each specific commandAvailable commands:
[eff|ann] : Annotate variants / calculate effects (you can use either 'ann' or 'eff', they mean the same). Default: ann (no command or 'ann').
build : Build a SnpEff database.
...
Generic options:
...
-version : Show version number and exitDatabase options:
-canon : Only use canonical transcripts.
...We will be using the eff/ann command by default. An example for working with sequencing data is given at this link. But we will deal with our data and work with any errors that are faced.
Since we are working with C. elegans so let’s search for any available databases for this organism or its reference genome.
java -jar snpEff/snpEff.jar databases | grep -i WBcel235
#Output
WBcel235.75 Caenorhabditis_elegans [https://snpeff.blob.core.windows.net/databases/v5_2/snpEff_v5_2_WBcel235.75.zip, https://snpeff.blob.core.windows.net/databases/v5_0/snpEff_v5_0_WBcel235.75.zip, https://snpeff.blob.core.windows.net/databases/v5_1/snpEff_v5_1_WBcel235.75.zip]
WBcel235.99 Caenorhabditis_elegans [https://snpeff.blob.core.windows.net/databases/v5_2/snpEff_v5_2_WBcel235.99.zip, https://snpeff.blob.core.windows.net/databases/v5_0/snpEff_v5_0_WBcel235.99.zip, https://snpeff.blob.core.windows.net/databases/v5_1/snpEff_v5_1_WBcel235.99.zip]I am going to use a recent release and download it. The downloaded database will be stored in the data directory under the snpeff directory snpEff/data/.
java -jar snpEff/snpEff.jar download WBcel235.99
After this, we can then just run the command to annotate and see if there are any errors.
java -jar snpEff/snpEff.jar ann WBcel235.99 vars/igv.vcf > annot1.vcf
# view the annotated file
bcftools view -H annot/annot1.vcf | lessIf you view the annotated VCF file, you may encounter that in the ANN tag, it says, Chromosome not found. Like this.
NC_003279.8 229204 . A T 67.9943 PASS DP=15;VDB=0.0511907;SGB=-0.616816;RPBZ=-2.06786;MQBZ=-3.45462;MQSBZ=-0.935627;BQBZ=0;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=5,3,3,3;MQ=36;ANN=T||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 0/1:101,0,198
NC_003279.8 229214 . G A 80.1209 PASS DP=16;VDB=0.0498254;SGB=-0.636426;RPBZ=-1.00633;MQBZ=-3.71086;MQSBZ=-0.412318;BQBZ=-1.34815;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=5,4,4,3;MQ=36;ANN=A||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 0/1:113,0,204
NC_003279.8 5163168 . T TG 53.0361 PASS INDEL;IDV=5;IMF=0.333333;DP=15;VDB=0.915812;SGB=-0.590765;RPBZ=1.35572;MQBZ=1.4344;MQSBZ=1.7186;BQBZ=-1.60696;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=7,3,0,5;MQ=33;ANN=TG||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 0/1:86,0,166
NC_003279.8 10271897 . C G 63.9635 PASS DP=15;VDB=0.211317;SGB=-0.556411;RPBZ=-1.13137;MQBZ=3.60555;MQSBZ=2.28035;BQBZ=0.506226;SCBZ=0;MQ0F=0;AC=2;AN=2;DP4=7,3,0,4;MQ=12;ANN=G||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 1/1:91,3,0
NC_003279.8 10462734 . C T 82.6088 PASS DP=257;VDB=6.18896e-05;SGB=-0.693147;RPBZ=0.518279;MQBZ=-5.56884;MQSBZ=-2.3572;BQBZ=-3.49299;SCBZ=0;MQ0F=0.14786;AC=1;AN=2;DP4=83,23,62,22;MQ=15;ANN=T||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 0/1:116,0,255
NC_003279.8 10462736 . C T 65.5 PASS DP=257;VDB=0.000298199;SGB=-0.693147;RPBZ=0.293845;MQBZ=-5.02518;MQSBZ=-1.32568;BQBZ=-0.707855;SCBZ=0;MQ0F=0.151751;AC=1;AN=2;DP4=90,22,64,15;MQ=15;ANN=T||MODIFIER|||||||||||||ERROR_CHROMOSOME_NOT_FOUND GT:PL 0/1:99,0,255Solution/Debugging
We will deal with this error in three simple steps and then run the command to annotate the variants.
Step 1. Look for the name of Chromosomes in the database downloaded by snpeff. The database is located in the snpeff folder. There are various ways to do so, we can just dump the database to view the names of chromosomes.
java -jar snpEff/snpEff.jar dump WBcel235.99 > WBcel235.99.db.txt vi WBcel235.99.db.txt
We can see the names of chromosomes in the dump output.
# Number of chromosomes : 7
# Chromosomes : Format 'chromo_name size codon_table'
# 'V' 20924180 Standard
# 'X' 17718942 Standard
# 'IV' 17493829 Standard
# 'II' 15279421 Standard
# 'I' 15072434 Standard
# 'III' 13783801 Standard
# 'MtDNA' 13794 Standard
Now we know the names of Chromosomes/Contigs are I, II, III, IV, V, X and MtDNA (which is in fact mitochondrial DNA).
OR
The other method to download the name of chromosomes, if you are automating is to use the bed file.
java -jar snpEff/snpEff.jar dump -bed WBcel235.99 > WBcel235.99.db.bed
more WBcel235.99.db.bed | cut -f1 | sort | uniq
# Outout is the name of chromosomes in the database
I
II
III
IV
MtDNA
V
X
Step 2. Look for the names of chromosomes in our VCF file.
We can check that in multiple ways. To view it from header:
bcftools view -h vars/igv.vcf | grep contig
# output
##contig=<ID=NC_003279.8,length=15072434>
##contig=<ID=NC_003280.10,length=15279421>
##contig=<ID=NC_003281.10,length=13783801>
##contig=<ID=NC_003282.8,length=17493829>
##contig=<ID=NC_003283.11,length=20924180>
##contig=<ID=NC_003284.9,length=17718942>
##contig=<ID=NC_001328.1,length=13794>To view it shortly, if you are automating, you can:
bcftools view -H vars/igv.vcf | cut -f1 | sort | uniq
# output
NC_001328.1
NC_003279.8
NC_003280.10
NC_003281.10
NC_003282.8
NC_003283.11
NC_003284.9Now we know the names of chromosomes that are used in our VCF files. These names are NC_001328.1, NC_003279.8, NC_003280.10, NC_003281.10, NC_003282.8, NC_003283.11 and NC_003284.9.
Step 3. Replace the chromosome names.
Now if you go to the same page from where we downloaded our reference genome at this URL https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000002985.6/ and scroll down to the Chromosomes and Annotation Details, there you can see the different names used and know how can we replace the names.
Now, we will replace the chromosome names using BCFtools annotate. The command to use is — rename-chr and it takes a file as input. The structures of the file can be seen from its online documentation. You can view at this link.
Now, create that file first.
vi names_to_change.txt
# OR
bcftools view -H vars/igv.vcf | cut -f1 | sort | uniq > names_to_change.txt
# and then update the file so that each of the old name is followed by a
# space and then the new name (the name used in snpeff database)# for our case the contents of file should look similar to this as
more names_to_change.txt
# output
NC_001328.1 MtDNA
NC_003279.8 I
NC_003280.10 II
NC_003281.10 III
NC_003282.8 IV
NC_003283.11 V
NC_003284.9 XFinally, replace the names by using the following command.
bcftools annotate --rename-chr names_to_change.txt vars/igv.vcf -Ov -o updated_names_igv.vcfNow if you view the chromosome names, you can see they are now updated in this new file.
bcftools view -H updated_names_igv.vcf | cut -f1 | sort | uniq
# output
I
II
III
IV
MtDNA
V
XBack to Annotating using SnpEff
Now annotate this new file using SnpEff.
java -jar snpEff/snpEff.jar ann WBcel235.99 updated_names_igv.vcf > annot2.vcf
Let’s see if the file is annotated now or not.
bcftools view -H annot2.vcf | less
# output
I 229204 . A T 67.9943 PASS DP=15;VDB=0.0511907;SGB=-0.616816;RPBZ=-2.06786;MQBZ=-3.45462;MQSBZ=-0.935627;BQBZ=0;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=5,3,3,3;MQ=36;ANN=T|downstream_gene_variant|MODIFIER|Y48G1BM.6|WBGene00021671|transcript|Y48G1BM.6|pseudogene||n.*2488T>A|||||2488|,T|intergenic_region|MODIFIER|Y48G1BL.5-Y48G1BM.6|WBGene00021664-WBGene00021671|intergenic_region|WBGene00021664-WBGene00021671|||n.229204A>T|||||| GT:PL 0/1:101,0,198
I 229214 . G A 80.1209 PASS DP=16;VDB=0.0498254;SGB=-0.636426;RPBZ=-1.00633;MQBZ=-3.71086;MQSBZ=-0.412318;BQBZ=-1.34815;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=5,4,4,3;MQ=36;ANN=A|downstream_gene_variant|MODIFIER|Y48G1BM.6|WBGene00021671|transcript|Y48G1BM.6|pseudogene||n.*2478C>T|||||2478|,A|intergenic_region|MODIFIER|Y48G1BL.5-Y48G1BM.6|WBGene00021664-WBGene00021671|intergenic_region|WBGene00021664-WBGene00021671|||n.229214G>A|||||| GT:PL 0/1:113,0,204
I 5163168 . T TG 53.0361 PASS INDEL;IDV=5;IMF=0.333333;DP=15;VDB=0.915812;SGB=-0.590765;RPBZ=1.35572;MQBZ=1.4344;MQSBZ=1.7186;BQBZ=-1.60696;SCBZ=0;MQ0F=0;AC=1;AN=2;DP4=7,3,0,5;MQ=33;ANN=TG|frameshift_variant|HIGH|let-630|WBGene00002804|transcript|Y110A7A.19.1|protein_coding|4/5|c.1576dupC|p.Gln526fs|1579/2097|1576/1995|526/664||,TG|upstream_gene_variant|MODIFIER|ppw-2|WBGene00004094|transcript|Y110A7A.18.1|protein_coding||c.-965dupC|||||938|,TG|upstream_gene_variant|MODIFIER|hcp-6|WBGene00001833|transcript|Y110A7A.1b.1|protein_coding||c.-1989_-1988insG|||||1983|,TG|upstream_gene_variant|MODIFIER|hcp-6|WBGene00001833|transcript|Y110A7A.1a.1|protein_coding||c.-1989_-1988insG|||||1983|,TG|downstream_gene_variant|MODIFIER|trf-2|WBGene00022454|transcript|Y110A7A.2.1|protein_coding||c.*4701_*4702insG|||||4702|;LOF=(let-630|WBGene00002804|1|1.00) GT:PL 0/1:86,0,166Now, we can see the error of Chromosome not found is resolved.
If we want to see the annotations only. You can do so as I am doing to see only the first two annotations.
bcftools view -H annot2.vcf | cut -f8 | cut -d ';' -f 14 | head -2
# output ANN=T|downstream_gene_variant|MODIFIER|Y48G1BM.6|WBGene00021671|transcript|Y48G1BM.6|pseudogene||n.*2488T>A|||||2488|,T|intergenic_region|MODIFIER|Y48G1BL.5-Y48G1BM.6|WBGene00021664-WBGene00021671|intergenic_region|WBGene00021664-WBGene00021671|||n.229204A>T||||||
ANN=A|downstream_gene_variant|MODIFIER|Y48G1BM.6|WBGene00021671|transcript|Y48G1BM.6|pseudogene||n.*2478C>T|||||2478|,A|intergenic_region|MODIFIER|Y48G1BL.5-Y48G1BM.6|WBGene00021664-WBGene00021671|intergenic_region|WBGene00021664-WBGene00021671|||n.229214G>A||||||If you want to further work on or understand about the ANN field or annotations contained in your VCF file refer to the link and learn about the various possible values and how you can filter or view them.
A Smarter Solution
If you are facing the above error and you know you are working with only one species and you have many samples, then instead of doing so every time on each of the files. You can just replace the name of Chromosomes in your reference genome and index it again. Because the names of chromosomes are passed from the reference genome to the mapping SAM/BAM files and hence the VCF files. So if you know that your database has different names than the reference genome, just rename the chromosomes/contigs in your reference genome.
Exploring ANN Field
ANN field represents different fields that you can see listed from the above link. We can explore the effects according to that.
For starting a basic Genomic analysis from scratch, look at the following article.
Twitter: vet_abdullah_n
