<p>I have some SARS-CoV-2 sequencing data that I'm trying to annotate with SnpEff, however SnpEff doesn't appear to be recognizing multiple adjacent SNPs within the same codon and correctly calling them as a single MNP (multiple-nucleotide polymorphism).</p>
<p>I'm using the GenBank reference MN908947.3. Here are the relevant lines of the annotated VCF file</p>
<pre><code>MN908947.3  28280   .   G   C   2911.06 PASS    AC=2;AF=1.00;AN=2;DP=65;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.56;SOR=8.380;ANN=C|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.7G&gt;C|p.Asp3His|7/1260|7/1260|3/419||,C|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1278G&gt;C|||||1278|,C|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2896G&gt;C|||||2896|,C|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2060G&gt;C|||||2060|,C|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1808G&gt;C|||||1808|,C|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1089G&gt;C|||||1089|,C|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*893G&gt;C|||||893|,C|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*521G&gt;C|||||521|,C|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*21G&gt;C|||||21| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2925,196,0
MN908947.3  28281   .   A   T   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.91;SOR=8.380;ANN=T|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.8A&gt;T|p.Asp3Val|8/1260|8/1260|3/419||,T|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1277A&gt;T|||||1277|,T|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2897A&gt;T|||||2897|,T|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2061A&gt;T|||||2061|,T|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1809A&gt;T|||||1809|,T|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1090A&gt;T|||||1090|,T|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*894A&gt;T|||||894|,T|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*522A&gt;T|||||522|,T|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*22A&gt;T|||||22| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0
MN908947.3  28282   .   T   A   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.71;SOR=8.380;ANN=A|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.9T&gt;A|p.Asp3Glu|9/1260|9/1260|3/419||,A|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1276T&gt;A|||||1276|,A|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2898T&gt;A|||||2898|,A|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2062T&gt;A|||||2062|,A|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1810T&gt;A|||||1810|,A|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1091T&gt;A|||||1091|,A|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*895T&gt;A|||||895|,A|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*523T&gt;A|||||523|,A|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*23T&gt;A|||||23| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0
</code></pre>
<p>These are the positions corresponding to codon 3 on the N gene, they should be annotated as a single MNP, p.Asp3Leu, but istead they've been annotated seperately, as three individual SNPs.</p>
<p>I'm using the latest version of SnpEff (v5), and according to the manual SnpEff should have this functionality, what am I doing wrong?</p>
<p>You can see an example of how it's supposed to work <a href="https://pcingola.github.io/SnpEff/se_inputoutput/" rel="nofollow noreferrer">here</a>.</p>

CC BY-SA 4.0
2021-04-09T12:57:43.817Z
Are they expressed as a single variant in the input file? Does the manual say snpEFF can detect MNPs or only that it can handle them?
terdon
CC BY-SA 4.0
2021-04-09T16:19:08.290Z
@terdon no, in the input file they're expressed as seperate SNPs. Based on the example in the link I provided, my take away from the SnpEff manual was that adjacent SNPs on the same codon would be grouped as an MNP, however there is little elaboration so I can't say for sure what is required for SnpEff to call adjacent SNPs as an MNP. I should note that since the variant caller does not have access to the ORFs, I don't think it would be able to call MNPs, so it has to be snpEff that does it.
Set
CC BY-SA 4.0
2021-04-09T16:23:53.733Z
The ORFs are not relevant here. In order to call a MNP, the variant caller needs to find all alt residues supported _by the same read_. Given that in your example, the reference allele has 0 support in all three, it seems very likely that the caller would find all three in the same read, but I don't see how snpEFF can be expected to _change_ the list of variants you give it. If you join these yourself, snpEFF will likely deal with it, but I don't see anything in the manual that suggests it can merge them for you.
terdon
CC BY-SA 4.0
2021-04-09T18:02:04.663Z
@terdon The ORFs are absolutely relevant here. The purpose of calling adjacent SNPs as a single MNP is so that the amino acid mutation is correctly annotated, thus you need to know where the breaks are between each codon, and cannot be deduced by the variant caller just from the ref. genome. Again, if you click on the link that I provided, and take a look at the example in the SnpEff manual, you will see how SnpEff annotates separate adjacent SNPs as a single MNP.
Set
CC BY-SA 4.0
2021-04-09T18:02:11.650Z
There is of course the problem of needing to determine phase (in the case diploid organisms), or even in the case of haploid organisms, that, as you say, SnpEff cannot in general determine from the VCF file whether adjacent SNPs always occur together on the same read. However I can find no info in the manual about how SnpEff handles these types of ambiguities, and my attempts at modifying the INFO field in the VCF file to indicate ploidy has not fixed the issue.
Set
CC BY-SA 4.0
2021-04-09T18:07:31.327Z
No, the ORFs are not relevant because the variant calling has no concept of genes, let alone translated products or codons. The calling is done using the reference genome and the reads. Of course the ORfs are relevant when you try to figure out the effect of your variants, but they do not affect and cannot affect whether something should be considered one variant or many. That is a function of the reads alone. All I see in the manual is an example of MNP, not a case where the input VCF had them separate and snpEFF joined them. I may be missing it, of course, but that's all I can find.
terdon
CC BY-SA 4.0
2021-04-09T18:09:14.603Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/122832/discussion-between-terdon-and-thoth).
terdon
<sequence-annotation><variant-calling><snpeff>
SnpEff not correctly annotating multiple-nucleotide polymorphisms
# SnpEff not correctly annotating multiple-nucleotide polymorphisms

(Asked by: username_0 on 2021-04-08)

I have some SARS-CoV-2 sequencing data that I'm trying to annotate with SnpEff, however SnpEff doesn't appear to be recognizing multiple adjacent SNPs within the same codon and correctly calling them as a single MNP (multiple-nucleotide polymorphism).
I'm using the GenBank reference MN908947.3. Here are the relevant lines of the annotated VCF file
MN908947.3  28280   .   G   C   2911.06 PASS    AC=2;AF=1.00;AN=2;DP=65;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.56;SOR=8.380;ANN=C|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.7G>C|p.Asp3His|7/1260|7/1260|3/419||,C|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1278G>C|||||1278|,C|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2896G>C|||||2896|,C|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2060G>C|||||2060|,C|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1808G>C|||||1808|,C|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1089G>C|||||1089|,C|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*893G>C|||||893|,C|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*521G>C|||||521|,C|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*21G>C|||||21| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2925,196,0
MN908947.3  28281   .   A   T   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.91;SOR=8.380;ANN=T|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.8A>T|p.Asp3Val|8/1260|8/1260|3/419||,T|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1277A>T|||||1277|,T|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2897A>T|||||2897|,T|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2061A>T|||||2061|,T|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1809A>T|||||1809|,T|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1090A>T|||||1090|,T|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*894A>T|||||894|,T|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*522A>T|||||522|,T|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*22A>T|||||22| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0
MN908947.3  28282   .   T   A   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.71;SOR=8.380;ANN=A|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.9T>A|p.Asp3Glu|9/1260|9/1260|3/419||,A|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1276T>A|||||1276|,A|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2898T>A|||||2898|,A|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2062T>A|||||2062|,A|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1810T>A|||||1810|,A|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1091T>A|||||1091|,A|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*895T>A|||||895|,A|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*523T>A|||||523|,A|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*23T>A|||||23| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0

These are the positions corresponding to codon 3 on the N gene, they should be annotated as a single MNP, p.Asp3Leu, but istead they've been annotated seperately, as three individual SNPs.
I'm using the latest version of SnpEff (v5), and according to the manual SnpEff should have this functionality, what am I doing wrong?
You can see an example of how it's supposed to work here.

**Tags:** sequence-annotation, variant-calling, snpeff
> **Comments:**
> **username_1**: Are they expressed as a single variant in the input file? Does the manual say snpEFF can detect MNPs or only that it can handle them?  
> **username_0**: @terdon no, in the input file they're expressed as seperate SNPs. Based on the example in the link I provided, my take away from the SnpEff manual was that adjacent SNPs on the same codon would be grouped as an MNP, however there is little elaboration so I can't say for sure what is required for SnpEff to call adjacent SNPs as an MNP. I should note that since the variant caller does not have access to the ORFs, I don't think it would be able to call MNPs, so it has to be snpEff that does it.  
> **username_1**: The ORFs are not relevant here. In order to call a MNP, the variant caller needs to find all alt residues supported _by the same read_. Given that in your example, the reference allele has 0 support in all three, it seems very likely that the caller would find all three in the same read, but I don't see how snpEFF can be expected to _change_ the list of variants you give it. If you join these yourself, snpEFF will likely deal with it, but I don't see anything in the manual that suggests it can merge them for you.  
> **username_0**: @terdon The ORFs are absolutely relevant here. The purpose of calling adjacent SNPs as a single MNP is so that the amino acid mutation is correctly annotated, thus you need to know where the breaks are between each codon, and cannot be deduced by the variant caller just from the ref. genome. Again, if you click on the link that I provided, and take a look at the example in the SnpEff manual, you will see how SnpEff annotates separate adjacent SNPs as a single MNP.  
> **username_0**: There is of course the problem of needing to determine phase (in the case diploid organisms), or even in the case of haploid organisms, that, as you say, SnpEff cannot in general determine from the VCF file whether adjacent SNPs always occur together on the same read. However I can find no info in the manual about how SnpEff handles these types of ambiguities, and my attempts at modifying the INFO field in the VCF file to indicate ploidy has not fixed the issue.  
> **username_1**: No, the ORFs are not relevant because the variant calling has no concept of genes, let alone translated products or codons. The calling is done using the reference genome and the reads. Of course the ORfs are relevant when you try to figure out the effect of your variants, but they do not affect and cannot affect whether something should be considered one variant or many. That is a function of the reads alone. All I see in the manual is an example of MNP, not a case where the input VCF had them separate and snpEFF joined them. I may be missing it, of course, but that's all I can find.  
> **username_1**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/122832/discussion-between-terdon-and-thoth).  

---

## 0 answers

