Develop a module that converts between coordinate systems

[juhtornr]

@mfiume commented on Thu Nov 03 2016

It would be nice to install into the Beacon Network a module that translates from:

• dbSNP IDs
• HGVS coordinates
• LRG
• etc.
  into the genomic coordinate system so that we can issue queries to Beacons regardless of input.

---

@antbro commented on Thu Nov 03 2016

Thats a knurly challenge
But one we're working hard on here in Leicester (led by Prof Raymond
Dalgleish)
I cc him in, so he can summarise the interface and functionality they
have build (and are now optimising) on top of Reece Harts hgvs python
library and Universal Transcript Archive (UTA)
FYI... they were key in getting the several thousand EDS mutations from
HGVS names to genome coordinates, for the Vancouver demo
Cheers
Tony

Marc Fiume wrote:

It would be nice to install into the Beacon Network a module that
translates from:

* dbSNP IDs
* HGVS coordinates
* LRG
* etc.
  into the genomic coordinate system so that we can issue queries
  to Beacons regardless of input.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
ga4gh-beacon/specification#61, or mute the thread
https://github.com/notifications/unsubscribe-auth/AI_EVJE222jGGHkhBIeg7W1IbnCYGu-8ks5q6fFjgaJpZM4KobqZ.

---

@antbro commented on Thu Nov 03 2016

Thats a knurly challenge
But one we're working hard on here in Leicester (led by Prof Raymond
Dalgleish)
I cc him in, so he can summarise the interface and functionality they
have build (and are now optimising) on top of Reece Harts hgvs python
library and Universal Transcript Archive (UTA)
FYI... they were key in getting the several thousand EDS mutations from
HGVS names to genome coordinates, for the Vancouver demo
Cheers
Tony

Marc Fiume wrote:

It would be nice to install into the Beacon Network a module that
translates from:

* dbSNP IDs
* HGVS coordinates
* LRG
* etc.
  into the genomic coordinate system so that we can issue queries
  to Beacons regardless of input.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
ga4gh-beacon/specification#61, or mute the thread
https://github.com/notifications/unsubscribe-auth/AI_EVJE222jGGHkhBIeg7W1IbnCYGu-8ks5q6fFjgaJpZM4KobqZ.

---

@skeenan commented on Mon Nov 07 2016

This can be achieved using the Ensembl REST services

Thanks to @andrewyatz for this info

dbSNP IDs

http://rest.ensembl.org/variation/human/rs56116432?content-type=application/json

Requires knowing the species but brings back things like

{
      "location": "9:133256042-133256042",
      "assembly_name": "GRCh38",
      "end": 133256042,
      "seq_region_name": "9",
      "strand": 1,
      "coord_system": "chromosome",
      "allele_string": "C/T",
      "start": 133256042
    },

HGVS coordinates

http://rest.ensembl.org/vep/human/hgvs/AGT:c.803T%3EC?content-type=application/json

    "assembly_name": "GRCh38",
    "end": 230710048,
    "seq_region_name": "1",
    "strand": -1,
    "id": "AGT:c.803T>C",
    "allele_string": "T/C",
    "most_severe_consequence": "missense_variant",
    "start": 230710048

LRG

For LRG, you can get Reference genomic coordinate through the EBI REST API, e.g.:

GRCh37: http://www.ebi.ac.uk/ebisearch/ws/rest/lrg/entry/LRG_1?fields=chr_name_grch37,chr_start_grch37,chr_end_grch37,chr_strand_grch37&format=json

GRCh38: http://www.ebi.ac.uk/ebisearch/ws/rest/lrg/entry/LRG_1?fields=chr_name_grch38,chr_start_grch38,chr_end_grch38,chr_strand_grch38&format=json

Here is the list of the other LRG "fields" available
http://www.ebi.ac.uk/ebisearch/ws/rest/lrg

---

@antbro commented on Thu Nov 10 2016

Hi All
Raymond Dalgleish has summarised the issues and a good solution very
nicely/thoroughly below, and asked me to forward this to the group.
[seems he cannot post to the list directly]
He really knows this stuff inside-out, and is creating solutions
Cheers
Tony

-------- Original Message --------
Subject: Re: [ga4gh/beacon-team] Develop a module that converts between
coordinate systems (#61)
Date: Thu, 10 Nov 2016 16:22:09 +0000
From: Dalgleish, Raymond W.M. (Prof.) [email protected]
To: Brookes, Anthony J. (Prof.) [email protected]

Dear all,

First, apologies for the length of this response. However, it's important to fully explore these issues and to draw attention to points which may not have been fully appreciated. I would encourage you to read beyond the first few lines so that you can see in full what we propose.

Short Response

I think that the problem boils down to something rather straightforward. There is essentially a need to perform the following data transformation:

variant description in the context of a reference sequence that's not chromosomal => HGVS-compliant variant description in the context of a chromosomal reference sequence

This transformation can be carried out accurately using VariantValidator (https://variantvalidator.org/). VariantValidator was conceived originally as an alternative to Mutalyzer but we have extended its capabilities extensively in recent weeks and an API should be ready for full testing in a day or two. VariantValidator will be able to create VCF files from HGVS-compliant descriptions and is fully aware of both 5´ and 3´ normalisation (shuffling)

These newly developed features are not yet available on our live server, but are in testing. However, the API basically works. The following example shows a chromosomal variant being mapped to a specific transcript of the TSC2 gene:

pjf9@login188:~> curl http://127.0.0.1:5000/GRCh37/'Chr16:2099572TC>T'/'NM_001318829.1'
[
{
"HGVS_RefSeqGene_description": "NG_005895.1:g.5269del",
"HGVS_genomic_description": "NC_000016.9:g.2099575del",
"HGVS_predicted_protein_consequence": "NP_001305758.1:p.?",
"HGVS_transcript_description": "NM_001318829.1:c.-9-826del",
"gene_symbol": "TSC2",
"intronic_variation_validated_against": "NC_000016.9",
"sumbitted_variant": "Chr16:2099572TC>T",
"validation_warnings": "Non-HGVS compliant variant description Chr16:2099572TC>T automapped to NC_000016.9:g.2099572TC>T: Non HGVS compliant variant description NC_000016.9:g.2099572TC>T automapped to NC_000016.9:g.2099575delC"
}
]

Please notice especially that the API produces "validation_warnings".

We humbly propose that the VariantValidator API may be more accurate than other possible solutions. Please read the full response for the details.

Full Response

Let's begin by summarising the history of this thread. The original query came from Marc Fiume who wrote:

It would be nice to install into the Beacon Network a module that
translates from:

 * dbSNP IDs
 * HGVS coordinates
 * LRG
 * etc.
   into the genomic coordinate system so that we can issue queries
   to Beacons regardless of input.

The thinking behind the request is entirely reasonable. It's desirable that the Beacon Network should be able to accept variants queries which are submitted in different valid formats. The key issue here is that queries must be unambiguous and adherence to standards is essential to avoid erroneous results being returned from a Beacon request. A dbSNP ID should presents no problem as an ID uniquely identifies a variant in the genome. However, the ID is entirely "divorced" from the actual variant to which it's assigned. Transposing two digits of an ID entirely changes the variant to which it refers and there is no simple way (any way?) to trap that an input error has occurred. An HGVS-complaint variant description is somewhat less prone to this type of error because the description can, at least, be parsed to ensure that it is valid (even if it is incorrect).

I'm uncertain what Marc is asking with respect to "LRG". LRG is simply a sophisticated sequence reference system which can be used as an alternative to RefSeq or Ensembl reference sequences. Indeed, every sequence (genomic, transcript& amino acid) within an LRG maps precisely onto corresponding RefSeq or Ensembl reference sequences. It ought to be relatively trivial to adapt existing tools to handle LRG sequences if support is currently not implemented.

Stephen Keenan then contributed that Marc's needs can be met using Ensembl REST services. The part that needs most attention is his suggestion with respect to HGVS coordinates:

  |http://rest.ensembl.org/vep/human/hgvs/AGT:c.803T%3EC?content-type=application/json
     >       {
     >            "assembly_name": "GRCh38",
     >            "end": 230710048,
     >            "seq_region_name": "1",
     >            "strand": -1,
     >            "id": "AGT:c.803T>C",
     >            "allele_string": "T/C",
     >            "most_severe_consequence": "missense_variant",
     >            "start": 230710048
     >            }

The key issue here is that the variant description provided as input into VEP (AGT:c.803T>C ) is not HGVS compliant.

Although VEP can analyse variant descriptions in this format, use of a gene symbol as a proxy for an accession number and version is not allowed in the HGVS nomenclature. As it happens, the AGT gene has a single known RefSeq transcript variant (splice variant), NM_000029.3. A syntactically correct and valid HGVS description would be NM_000029.3:c.803T>C. In addition, the following descriptions are also valid, NC_000001.11:g.230710048A>G (chromosome), NG_008836.1:g.9543T>C (RefSeqGene) and the predicted effect on the protein NP_000020.1:p.(Met268Thr).

The implicit suggestion that the format gene_symbol:variant_descrption is syntactically correct sends a misleading message to those who are not fully acquainted with the HGVS nomenclature.

VEP appears to be aware that there is a single transcript for the AGT gene and maps the variant appropriately to that transcript. However, what would happen for a hypothetical gene that has two coding transcripts for each of which a single variant description is valid, but for which the altered bases of the two transcripts map to different genome locations?

This can be illustrated by the COL5A1 gene. There are two coding transcripts for the COL5A1 gene (NM_000093.4 and NM_001278074.1) which encode mutually exclusive alternative versions of exon 64 (64A and 64B respectively). The variant descriptions NM_000093.4:c.5071A>G and NM_001278074.1:c.5071A>G are both valid and syntactically correct.

Variant NM_000093.4:c.5071A>G maps to NC_000009.12:g.134829979A>G and to NM_001278074.1:c.5068-129A>G.

Variant NM_001278074.1:c.5071A>G maps to NC_000009.12:g.134830111A>G and to NM_000093.4:c.5136+67A>G.

COL5A1:c.5071A>G when inputted into VEP correctly identifies the variant NM_000093.4:c.5071A>G, maps it to NC_000009.12:g.134829979A>G and also to NM_001278074.1:c.5068-129A>G. However, VEP seems to be unaware that the variant description c.5071A>G might be valid in the context of NM_001278074.1.

Additional code could be added to VEP to correct this anomalous behaviour, but it would be better that inputs such as COL5A1: c.5071A>G be rejected by VEP as not being HGVS-compliant.

There is a further issue here. Formally, the variant description NM_001278074.1:c.5068-129A>G is invalid. It's impossible in the context of NM_001278074.1 to verify that the base at intron positon -129 relative to base 5068 is an A. If there is the desire to indicate clearly that the variant is intronic, the correct description is NG_008030.1(NM_001278074.1):c.5068-129A>G. VariantValidator (https://variantvalidator.org/) informs the user that NG_008030.1(NM_001278074.1):c.5068-129A>G is correct if the inputted query is NM_001278074.1:c.5068-129A>G.

The Ensembl APIs do not, as far as we are aware, correct for mismatches between the genomic sequence and corresponding transcript sequences at equivalent positions. For example, NC_000001.10:g.216219781A>C (GRCh37) is the genomic-level variant description for a variant in the USH2A gene. However, the corresponding valid transcript-level variant description is NM_206933.2:c.6317C>G. The reason for the apparent mismatch is that the allele represented in the chromosome differs from that in the RefSeq transcript record (NM_206933.2).

VEP produces the outputs "WARNING: Unable to parse HGVS notation 'USH2A:c.6317C>G'" and "WARNING: Unable to parse HGVS notation 'NM_206933.2:c.6317C>G'". The latter warning suggests that VEP maps the variant position first to the genome to find the genomic base, but does not verify that c.6317 really is a C base in NM_206933.2. VariantValidator correctly handles this mismatch: see attached screenshot.

Single-base mismatches are not uncommon between corresponding reference sequences. They arise primarily because of community requests that the major allele be represented in transcript sequences which are, de facto, the primary reference sequences for the reporting of sequence variants in many genes, especially in diagnostic laboratories.

We are currently beta-testing a range of VariantValidator tools that can take VCF files and generate validated HGVS-compliant genomic variant descriptions. Additionally, the files can be filtered for variants that map to specific genes, and further still, the user can additionally specify a particular transcript to which the genomic position should be mapped.

From a VCF call, we are then able to output validated:

1. Genomic HGVS descriptions
2. Transcript-level descriptions (to a single user-specified transcript, a range of user-specified transcripts, or all transcripts overlapping a specified gene)
3. HGVS-compliant descriptions of the predicted effect on the protein(s).
4. Other relevant descriptions e.g. RefSeqGene-based variant descriptions.

Although these tools were developed as batch analysis tools for mainstream genetics users, they could easily be integrated into two specifically designed APIs. One which takes VCF calls (beacons) and maps to HGVS compliant descriptions, and one which takes HGVS descriptions and maps to VCF (beacons).

Best wishes,

Raymond

---

@fcunningham commented on Fri Nov 11 2016

Hi Raymond,

Thanks for your detailed response. We have written the VEP so that it can process 'AGT:c.803T>C' and to output correct HGVS. As endorsed by the HGVS rules, (http://varnomen.hgvs.org/recommendations/DNA/variant/substitution/), in a manuscript it is correct HGVS to use this shortened version e.g. AGT:c.803T>C (given certain conditions). Therefore, this format is a much sought after and useful use case. Given this I don't think it stands that we are misleading those who are not fully acquainted with the HGVS nomenclature.

We currently cannot elegantly handle the RefSeq transcripts when they differ from the genome but we do warn when this happens are we are working to improve this. Note that we correctly handle GENCODE transcripts which are simple to use for NGS workflows when aligning to the reference genome. We are working with RefSeq to use their alignments, rather than risk potential discrepancies in calculating our own.

The inclusion of genomic sequence accession for intronic annotation is not something we have come across previously, so we will look into it.

Hope that helps clarify things.

Best wishes,

Fiona

---

@antbro commented on Thu Nov 24 2016

Please see below, provided on behalf of Raymond who cannot presently post to this list...

Hi Fiona,

As things stand, VEP will not be able to perform the necessary conversions required by this work package. As we previously pointed out, it is unaware of discrepancy between transcript sequence and genome sequence. We have noticed that it cannot parse HGVS c. inversions which we have submitted. It also does not provide user-friendly warnings that enable users to correct their errors.

We now have the tool specifically requested in the title of this work package: i.e. VariantValidator can take an HGVS description, validate it and return a valid VCF. In fact, it can also be used to validate and re-format errant VCF calls.

It's unsafe in the context of clinical decision making to express variant descriptions in the format gene_symbol:variant_description and VariantValidator guards against the use of non-compliant descriptions.

Regards

Raymond

P.S. I'll contact you directly and address your individual points in more detail.

---

From: Fiona Cunningham [[email protected]]
Sent: 11 November 2016 18:48
To: ga4gh/beacon-team
Cc: Brookes, Anthony J. (Prof.); Comment
Subject: Re: [ga4gh/beacon-team] Develop a module that converts between coordinate systems (#61)

Hi Raymond,

Thanks for your detailed response. We have written the VEP so that it can process 'AGT:c.803T>C' and to output correct HGVS. As endorsed by the HGVS rules, (http://varnomen.hgvs.org/recommendations/DNA/variant/substitution/), in a manuscript it is correct HGVS to use this shortened version e.g. AGT:c.803T>C (given certain conditions). Therefore, this format is a much sought after and useful use case. Given this I don't think it stands that we are misleading those who are not fully acquainted with the HGVS nomenclature.

We currently cannot elegantly handle the RefSeq transcripts when they differ from the genome but we do warn when this happens are we are working to improve this. Note that we correctly handle GENCODE transcripts which are simple to use for NGS workflows when aligning to the reference genome. We are working with RefSeq to use their alignments, rather than risk potential discrepancies in calculating our own.

The inclusion of genomic sequence accession for intronic annotation is not something we have come across previously, so we will look into it.

Hope that helps clarify things.

Best wishes,

Fiona

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHubga4gh-beacon/specification#61 (comment), or mute the threadhttps://github.com/notifications/unsubscribe-auth/AI_EVFle88XCYdwPoNqqCWlIM8TWHYV0ks5q9LiBgaJpZM4KobqZ.