input in HGVSp and HGVSc gives unexpected REF allele in reversed strand genes

[FedericaBrando]

Describe the issue

I'm not sure if it's an issue or I am not fully getting the rational behind this behaviour, so I would like to understand more the results I get.

When inputting HGVSp (protein sequence) and HGVSc (coding sequence) in VEP of genes that are on the reversed strand I get opposite nucleotides in the REF and ALT. I was expecting the REF to be standard in giving the 5'->3' in the output. But this doesn't seem the case.

Additional information

For example, I have the mutation: NRAS:p.Gly13Asp . If I ran VEP on this mutation in HGVS format and in genomic coordinates, this is what I get as result:

Uploaded_variants	NRAS:p.Gly13Asp	1_114716123_G/A	1_114716123_C/T
Location	1:114716123	1:114716123	1:114716123
REF_ALLELE	G	G	C
Allele	A	A	T
Gene	ENSG00000213281	ENSG00000213281	ENSG00000213281
Consequence	missense_variant	missense_variant	missense_variant
SYMBOL	NRAS	NRAS	NRAS
HGNC_ID	HGNC:7989	HGNC:7989	HGNC:7989
ENSP	ENSP00000358548	ENSP00000358548	ENSP00000358548
CANONICAL	YES	YES	YES
Amino_acids	G/D	A/V	G/D
Protein_position	13	13	13
Feature	ENST00000369535	ENST00000369535	ENST00000369535
EXON	2/7	2/7	2/7
INTRON	-	-	-
Codons	gGt/gAt	gCt/gTt	gGt/gAt
STRAND	-1	-1	-1
HGVSp	ENSP00000358548.4:p.Gly13Asp	ENSP00000358548.4:p.Ala13Val	ENSP00000358548.4:p.Gly13Asp
HGVSc	ENST00000369535.5:c.38G>A	ENST00000369535.5:c.38G>T	ENST00000369535.5:c.38G>A
HGVSg	-	-	-
MANE_SELECT	NM_002524.5	NM_002524.5	NM_002524.5

Even if the gene is on the reverse strand I am expecting the REF to be C, I do not fully understand this results. Could I have an explanation?

System

I am using docker://ensemblorg/ensembl-vep:release_111.0

Full VEP command line

singularity exec docker://ensemblorg/ensembl-vep:release_111.0 vep  -i /tmp/tmps5xkl4ld_input.tsv -o STDOUT   --cache             --cache_version 111   --dir_cache /home/user/.vep  --assembly GRCh38 --force_overwrite   --offline --hgvs  --tab  --fields Uploaded_variation,Location,REF_ALLELE,Allele,Gene,Consequence,SYMBOL,HGNC_ID,ENSP,CANONICAL,Amino_acids,Protein_position,Feature,EXON,INTRON,Codons,STRAND,HGVSp,HGVSc,HGVSg,MANE_SELECT  --show_ref_allele   --symbol             --canonical             --protein             --numbers             --no_stats             --fork 4             --mane             --warning_file tmp/preprocessing/vep/vep_warnings.txt         

Full error message

This is the only warning that is actually not relevant I guess

WARNING: Possible invalid use of gene or protein identifier 'NRAS' as HGVS reference; most likely transcript will be selected

[likhitha-surapaneni]

Hi @FedericaBrando , sorry to hear that you are facing issues. Can you please upload the input file used (tmps5xkl4ld_input.tsv), if possible? This will help us in our investigation.

[FedericaBrando]

Here's the inputs:

• This is the protein format: tmpwqnkshru_input.txt
• This is the GTF tmpywb1f6s_input.txt

Now that I see the the input of the GTF (we are programmatically running VEP and we build the input on the fly) - the discrepancy of the second column it's probably due to the STRAND columns that is assigning + to both of the variants.
Hence If I assign:

1	114716123	114716123	G/A	-	1_114716123_G/A
1	114716123	114716123	C/T	+	1_114716123_C/T

I get the expected results. My conclusion is that HGVSp are mapped to the CDS on the strand where the gene is for the REF allele.

Another thing that I did is to run the online VEP tool for NM_002524.5:p.Gly13Asp which gives as as ALT a T (differently from above). While on VEP 111 gives an A - with the command reported above. Has anything changed between the two versions? (VEP online run

---

Although - This leaves me a doubt - for example in this run I am querying

17 7675235 7675235 T/C +
17 7675235 7675235 A/G -

And I get the REF allele in the first case T and in the second an A. Which is in line with what I just commented above.

Although the clinical significance is different. Actually - it is only present in the case where the REF is T. ClinVAR result

This reflects if we input the HGVS for that variant, which is: NP_000537.3:p.Tyr126Cys. In fact the VEP run maps as REF allele the A, leaving out the ClinVAR annotation and potentially a very important clinical annotation (run here).

[likhitha-surapaneni]

Hi @FedericaBrando ,
Thank you for the detailed explanation.

Regarding the first question for why NM_002524.5:p.Gly13Asp gives Allele column as T. The Allele column here represents the variant allele used to calculate the consequence. It may or may not be matching alternate allele corresponding to the reference allele column. It is worth downloading the VCF output to check if the results are as expected.

For the ClinVar result showing different results with different representations, this is due to the allele matching algorithm used. We have plans to improve this in the coming releases.

Hope this answers your question.

[FedericaBrando]

I see! Thanks for the clarification. Although I do understand the different use of Allele and result may differ I still don't quite understand why in VEP online I have the Allele mapped to the forward STRAND and - for the same HGVSp - I get the Allele on the reverse STRAND in vep command-line.

Another thing that I did is to run the online VEP tool for NM_002524.5:p.Gly13Asp which gives as as ALT a T (differently from above). While on VEP 111 gives an A - with the command reported above. Has anything changed between the two versions? (VEP online run

Is there any reason why this is happening?

---

From a batch analysis of HGVSp ran in VEP it seems that if a mutation (specifically an SNV) is in a gene on the reverse strand then VEP cmd-line will always map REF and Allele to the reverse STRAD. Do you think it would be safe to assume this and perform a mapping to the forward one from the tsv format output by VEP consistently?