HGVS Nomenclature06/11
Starting sandbox…
Beginner·3 graded tests

UTA and Data Providers

UTA and Data Providers

Core Principle: The Universal Transcript Archive (UTA) provides relational transcript-to-genome exon alignments and gene structure metadata necessary for variant coordinate projection and validation.

Why a database is needed

Converting NM_000546.6:c.215C>G to genomic coordinates requires knowing:

  • where that transcript aligns on the genome, exon by exon
  • which strand it is on
  • where the CDS starts and ends within it
  • whether the alignment has indels relative to the genome (transcripts and the reference genome genuinely disagree in places)

None of that is derivable from the variant string. It has to be looked up. UTA is that lookup.

Connecting

import hgvs.dataproviders.uta

hdp = hgvs.dataproviders.uta.connect()          # public remote instance by default
hdp = hgvs.dataproviders.uta.connect(
    "postgresql://anonymous@localhost:5432/uta/uta_20210129"
)

Or via the environment, which is the better pattern:

export UTA_DB_URL=postgresql://anonymous@localhost:5432/uta/uta_20210129

The URL's final path segment (uta_20210129) is the schema name, which is the dataset version. Different schemas are different snapshots of RefSeq/Ensembl alignments, and they can give different answers.

Build hdp once and pass it around. Each connect() opens a database connection; doing it per variant is catastrophic for performance.

Remote vs local

The default remote instance (uta.biocommons.org) is a free public service.

Use it for: learning, exploration, a handful of variants, testing.

Do not use it for: any pipeline, any batch job, anything time-sensitive. It is a shared courtesy resource: slow, rate-limited, and its availability is not guaranteed. Running 50,000 variants through it is both painful for you and antisocial.

Running UTA locally

Docker is the practical route:

docker run -d --name uta -p 5432:5432 biocommons/uta:uta_20210129
export UTA_DB_URL=postgresql://anonymous@localhost:5432/uta/uta_20210129

Roughly 5–10 GB on disk. The speed difference is 10–100× on batch work: a job that takes hours remotely takes minutes locally.

For production, restore the SQL dump into your own PostgreSQL instance rather than relying on the container.

What the provider offers

hdp.data_version()                                  # schema version
hdp.schema_version()
hdp.get_tx_info(tx_ac, alt_ac, alt_aln_method)      # CDS bounds, transcript info
hdp.get_tx_exons(tx_ac, alt_ac, alt_aln_method)     # exon structure
hdp.get_tx_for_gene("TP53")                          # all transcripts for a gene
hdp.get_tx_for_region(alt_ac, method, start, end)    # transcripts overlapping a locus
hdp.get_acs_for_protein_seq(seq)
hdp.get_gene_info("TP53")
hdp.get_similar_transcripts(tx_ac)
hdp.get_pro_ac_for_tx_ac(tx_ac)                      # transcript → protein accession

Most of the time you use these indirectly through AssemblyMapper, but they are useful for building transcript-selection logic.

import pandas as pd
tx = pd.DataFrame(hdp.get_tx_for_gene("TP53"),
                  columns=["hgnc", "cds_start_i", "cds_end_i", "tx_ac",
                           "alt_ac", "alt_aln_method"])
tx[tx["alt_ac"].str.startswith("NC_0000")]      # only primary assembly alignments

Alignment methods

[ ]python
out1ms

This code completed without producing output.

Use splign for RefSeq (NM_/NR_) transcripts. Mixing methods for the same transcript can give different genomic coordinates, because the aligners place indels differently in ambiguous regions. Pick one and record it.

Caching for speed and reproducibility

from hgvs.dataproviders.seqfetcher import SeqFetcher
import hgvs.dataproviders.uta

hdp = hgvs.dataproviders.uta.connect(pooling=True)

pooling=True reuses connections: a straightforward win for batch work.

For a fully offline and reproducible setup, biocommons provides caching wrappers, and hgvs's own test suite uses a cached dataset so tests need no network. In a pipeline, wrapping hdp in your own memoisation of get_tx_info/get_tx_exons is cheap and effective, since the same handful of transcripts is queried thousands of times.

Pin your data version

print(hdp.data_version())        # e.g. 'uta_20210129'

Record this in every output. Two facts make it essential:

  1. Normalization results depend on the sequence data, which depends on the UTA snapshot.
  2. Transcript sets change. A transcript present in one snapshot may be suppressed in the next, and CDS boundaries occasionally get corrected.

An analysis that does not record its UTA and SeqRepo versions is not reproducible, no matter how carefully you pinned your Python packages. See Applied - Reproducible Environment.

Common connection errors

Error Cause
HGVSDataNotAvailableError: No transcript definition for ... transcript not in this UTA snapshot: often a version mismatch
psycopg2.OperationalError database unreachable; check UTA_DB_URL and that the container is running
Very slow projection you are on the remote instance, or reconnecting per variant
No alignment for tx X on Y with method splign wrong alignment method, or the transcript is unplaced on that assembly

HGVSDataNotAvailableError with a versioned accession usually means that exact version is absent. Query get_tx_for_gene to see which versions the snapshot has, and decide deliberately whether to fall back to a different version: do not strip the version silently.

Alternative providers

The provider interface is pluggable:

[ ]python
out1ms

This code completed without producing output.

Alternatives exist: including implementations backed by Ensembl, and by SeqRepo-only setups for reference lookups. The openvar/vv_hgvs fork (VariantValidator) uses its own provider. Unless you have a specific need, UTA is the well-trodden path.

Common mistakes

  • Using the public remote instance for a pipeline.
  • Reconnecting per variant.
  • Not recording the UTA version.
  • Mixing alt_aln_method values.
  • Stripping accession versions when a lookup fails.
  • Assuming every transcript is present. UTA covers what RefSeq/Ensembl aligned; novel and suppressed transcripts are not there.
  • Ignoring HGVSDataNotAvailableError instead of reporting which variants could not be processed.

See also

SeqRepo · AssemblyMapper · hgvs easy · biocommons hgvs · hgvs Normalizer · Applied - Reproducible Environment

Exercise

The hgvs package and a UTA instance are both unavailable in this sandbox, so you will stand in a dict for the database and a second dict for the memoisation the lesson recommends. Counting the round-trips yourself is what makes the 10-100x difference concrete, and handling the version mismatch by hand is what stops you from stripping an accession version to make a lookup succeed.

Write serve(text), taking the whole of stdin and returning a list of output lines. Above a line -- stdin gives one tx <transcript> <genomic_accession> <alt_aln_method> line per row the snapshot holds. Below it come queries of the same three fields, separated by spaces. Answer each query in order and emit <transcript> <accession> <method> <status>:

  • hit: this exact triple has been answered from the database before, so no round-trip is needed
  • miss: the triple is in the table; count a round-trip and cache it
  • version_mismatch:<available>: the triple is absent but the table holds the same accession at a different version with the same genomic accession and method; count a round-trip and cache nothing
  • no_data: absent entirely; count a round-trip and cache nothing

Finish with queries=<n> lookups=<m> hit_rate=<r>, where the lookups are the round-trips and the hit rate is the hits divided by the queries, rounded to three decimal places and reported as 0.0 when there are no queries.

exerciseEvery submission runs independently in a fresh namespace.

Test cases · 3

#viainputexpected stdout
1stdintx NM_000546.6 NC_000017.11 splign tx NM_000546.6 NC_000017.11 blat tx NM_000492.4 NC_000007.14 splign --NM_000546.6 NC_000017.11 splign NM_000546.6 NC_000017.11 splign NM_000546.6 NC_000017.11 blat NM_000546.5 NC_000017.11 splign NM_000546.5 NC_000017.11 splign NM_007294.3 NC_000017.11 splign queries=0 lookups=0 hit_rate=0.0
2stdintx NM_000546.6 NC_000017.11 splign tx NM_000546.6 NC_000017.11 blat tx NM_000492.4 NC_000007.14 splign --NM_000492.4 NC_000007.14 splign NM_000492.4 NC_000007.14 blat NM_000492.4 NC_000017.11 splign NM_000492.4 NC_000007.14 splign queries=0 lookups=0 hit_rate=0.0
3stdintx NM_000546.6 NC_000017.11 splign tx NM_000546.6 NC_000017.11 blat tx NM_000492.4 NC_000007.14 splign -- queries=0 lookups=0 hit_rate=0.0

Hints · 3

01Hint

The cache key is the whole (transcript, accession, method) triple; caching on the transcript alone answers the wrong question.

02Hint

Split an accession on its last dot to compare the base against the table's rows.

03Hint

Only successful lookups go in the cache, so a repeated missing transcript costs a round-trip every time.