Support hg19/hg38 builds and extend output with multiple mappings?

[ifokkema]

The variantformatter API method to me seems to look like a lightweight alternative to the variantvalidator API method (but please correct me where I'm wrong). However, the variantvalidator API method supports GRCh37, GRCh38, hg19, hg38, and returns data for all four builds. The variantformatter API method allows for only GRCh37 and GRCh38, but uses only the given one in the output.

I think it would be very useful to also support the other builds here, or at least two builds (also see openvar/variantValidator#126) and to also provide mappings here to the other build(s).

[Peter-J-Freeman]

Hi Ivo, the VariantFormatter tool was designed as a straight-forward genome > transcript > protein tool. The code for supporting hg19 and hg38 could likely be easily integrated from VariantValidator. The issues comes from lifting between builds.

As we discussed previously, we use a transcript based mapping system which ensures that the c. description mapped from say GRCh37 is then mapped onto GRCh38. We cannot, therefore assume that a variant mapped onto all transcripts from GRCh37 will map onto the same genomic position for every single transcript in GRCh38.

The solution would be to provide a t_to_g liftover for each transcript, and VariantValidator has a module in place that will do this. So the coding should not take long, it's the design of the output.

Can you possibly provide a basic map of what the output JSON might look like? This will help me see what you want and inform the design. Please bear in mind the need for t_to_g liftover from each transcript.

[ifokkema]

I've merged part of the VariantValidator output into the VariantFormatter output, and come up with this:

{
  "NC_000015.9:g.40699840C>T": {
    "NC_000015.9:g.40699840C>T": {
      "g_hgvs": "NC_000015.9:g.40699840C>T",
      "genomic_variant_error": null,
      "hgvs_t_and_p": {
        "NM_002225.3": {
          "gap_statement": null,
          "gapped_alignment_warning": null,
          "p_hgvs_slc": "NP_002216.2:p.(R53C)",
          "p_hgvs_tlc": "NP_002216.2:p.(Arg53Cys)",
          "primary_assembly_loci": {
            "grch37": {
              "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
            },
            "grch38": {
              "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
            },
            "hg19": {
              "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
            },
            "hg38": {
              "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
            }
          },
          "t_hgvs": "NM_002225.3:c.157C>T",
          "transcript_variant_error": null
        }
      },
      "p_vcf": "15:40699840:C:T"
    },
    "errors": []
  },
  "metadata": {
    "seqrepo_db": "/local/seqrepo/2018-08-21",
    "uta_schema": "uta_20180821",
    "variantvalidator_hgvs_version": "1.2.5.vv1",
    "variantvalidator_version": "1.0.1.dev53+gc394c83"
  }
}

I immediately notice the inconsistency in naming between g_hgvs and hgvs_genomic_description, which are basically the same fields. There are multiple options here;

• Don't drop g_hgvs (but possibly fix naming inconsistency), but note that user should know the input to realize which build this is.
• Drop g_hgvs if there are transcripts that have mappings back to the genome && create a primary_assembly_loci element within the NC_000015.9:g.40699840C>T subelement to provide genome-to-genome mappings when no transcript is available;
• Combination of the above, always provide an primary_assembly_loci element for the genomic variant as well as for all transcripts (will be slow, I assume).

I think I might be messing with your API design too much, what do you think? I often need the mapper without needing protein prediction, so perhaps we should split this off? Or do you think that passing checkonly=False suffices?

[Peter-J-Freeman]

I think this is the first time the two tools have ever been used in combination. They were never intended to be used in this way hence never having considered the naming conventions. But that does not mean we shouldn't.

I'm pretty sure I can use the VV liftover tool in VF which will create the priomary assembly loci

I find it a little hard to visualise your ideas. Looking at this JSON, I would drop the g_hgvs. It is redundant. The reason to keep the primary assembly loci is that we cannot assume that we will always hit the same genomic position if we lift-over, so I would stick to outputting this on a transcript by transcript basis. I would not put this into the the NC_000015.9:g.40699840C>T sub-element for the same reason. i.e. if we get different GRCh37 or GRCh38 edits, one will be over-written, and we will lose all sense of why this might have arisen. Liftover within genes is transcript based, so we need the transcript info. Hope this makes sense. It is very topical at the moment.

I agree, we need a selected genome build flag, especially if we are providing the vcf. OOPS!

Here is my recommendation

{
  "NC_000015.9:g.40699840C>T": {
    "NC_000015.9:g.40699840C>T": {
      "genomic_variant_error": null,
      "selected_genome_build": "GRCh37"
         "hgvs_t_and_p": {
        "NM_002225.3": {
          "gap_statement": null,
          "gapped_alignment_warning": null,
          "p_hgvs_slc": "NP_002216.2:p.(R53C)",
          "p_hgvs_tlc": "NP_002216.2:p.(Arg53Cys)",
          "primary_assembly_loci": {
            "grch37": {
              "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
            },
            "grch38": {
              "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
            },
            "hg19": {
              "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
            },
            "hg38": {
              "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
            }
          },
          "t_hgvs": "NM_002225.3:c.157C>T",
          "transcript_variant_error": null
        }
      },
    },
    "errors": []
  },
  "metadata": {
    "seqrepo_db": "/local/seqrepo/2018-08-21",
    "uta_schema": "uta_20180821",
    "variantvalidator_hgvs_version": "1.2.5.vv1",
    "THIS WILL BECOME VF INFORMATION"
  }
}

Note - Remove the p_vcf as well and add to the primary assembly loci.

e.g.

      "hg38": {
        "hgvs_genomic_description": "NC_000017.11:g.50198002C>A",
        "vcf": {
          "alt": "A",
          "chr": "chr17",
          "pos": "50198002",
          "ref": "C"
        }

[ifokkema]

In principle, this looks great!
I forgot to remove the p_vcf so thanks for picking that up.

One note:

Looking at this JSON, I would drop the g_hgvs. It is redundant. The reason to keep the primary assembly loci is that we cannot assume that we will always hit the same genomic position if we lift-over, so I would stick to outputting this on a transcript by transcript basis. I would not put this into the NC_000015.9:g.40699840C>T sub-element for the same reason. i.e. if we get different GRCh37 or GRCh38 edits, one will be over-written, and we will lose all sense of why this might have arisen. Liftover within genes is transcript based, so we need the transcript info.

Yes, but, what do you want to do when there are no transcripts mapping to the given variant? Will you then drop the hgvs_t_and_p element and just provide the primary_assembly_loci instead? Like:

{
  "NC_000015.9:g.40699840C>T": {
    "NC_000015.9:g.40699840C>T": {
      "genomic_variant_error": null,
      "selected_genome_build": "GRCh37",
      "primary_assembly_loci": {
        "grch37": {
          "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
        },
        "grch38": {
          "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
        },
        "hg19": {
          "hgvs_genomic_description": "NC_000015.9:g.40699840C>T"
        },
        "hg38": {
          "hgvs_genomic_description": "NC_000015.10:g.40407639C>T"
        }
      }
    },
    "errors": []
  },
  "metadata": {
    "seqrepo_db": "/local/seqrepo/2018-08-21",
    "uta_schema": "uta_20180821",
    "variantvalidator_hgvs_version": "1.2.5.vv1",
    "THIS WILL BECOME VF INFORMATION"
  }
}

The client will need to then check for hgvs_t_and_p OR primary_assembly_loci if the former is absent.

Note - Remove the p_vcf as well and add to the primary assembly loci.

I don't use the VCF fields, but I assume getting them won't take much calculation. Otherwise, I'd like to make those optional as well 😉

[Peter-J-Freeman]

Good point! We cannot assume the key for the nested dict will be hgvs. Will keep the g_hgvs!

[Peter-J-Freeman]

I think that the simplest thing will be to simply update VF to provide this output. Will keep the p_vcf for now as well as same logic applies r.e. g_hgvs. Will add in a a way to halt at transcript level as well

[Peter-J-Freeman]

Actually, we can add the primary_assembly_loci regardless of hitting a transcript. This is slower for intragenic regions so actually I will make this an API method rather than a VF method.

[ifokkema]

Actually, we can add the primary_assembly_loci regardless of hitting a transcript. This is slower for intragenic regions so actually I will make this an API method rather than a VF method.

It would be a requirement for us indeed to get a mapping for intragenic regions. If then no transcripts are present, we need to be able to find the primary_assembly_loci somewhere else in the structure...

[Peter-J-Freeman]

7382c77 closes issue relating to support for hg19 and 38

[Peter-J-Freeman]

bdcbbe9 closes for VF

[Peter-J-Freeman]

openvar/rest_variantValidator@40e1d96 closes for REST API ready for testing. Will spin up on dev server later this afternoon

[Peter-J-Freeman]

Now running on https://www35.lamp.le.ac.uk/ ready for testing. Note @ifokkema. The updated versions of VV and VF need their pytests updating to compensate for known changes. Please let me know what you think

[ifokkema]

Well done, @PeteCausey-Freeman ! I am impressed!! 🤝

For my own understanding;

• Variant with mapping to transcripts, liftover information is in assembly_loci
• Variant without mapping to transcripts, liftover information is in primary_assembly_loci and assembly_loci is null.

Wat is the difference between primary_assembly_loci and assembly_loci? The API allows for arguments True, primary and False. What is the difference between True and primary?