How do I use local UTA server AND local seqrepo (or local fasta) together?

[sodiumnitrate]

I want to use the local UTA server I already set up for looking up transcripts, but don't want to use E-utils to fetch chromosomes from NCBI. I already have local fasta. I tried setting up and using seqrepo as well. When I use seqrepo, it looks for the transcripts in the seqrepo as well, which I don't want.

The following uses e-utils:

hdp = uta.connect(url)
hp = hgvs.parser.Parser()
am = hgvs.assemblymapper.AssemblyMapper(
    hdp,
    assembly_name='GRCh37',
    alt_aln_method="splign",
    replace_reference=True
)

gdot = hp.parse_hgvs_variant('NC_000015.10:g.123123_123125del')
norm = Normalizer(hdp, alt_aln_method="splign", validate=False)
res = norm.normalize(gdot)

If I set:

import os
os.environ["HGVS_SEQREPO_DIR"] = "/path/to/seqrepo"

then it ignores the UTA server altogether.

I feel like this should be doable, but I haven't been able to find a working solution. How do I do this? (version 1.5.7)

[korikuzma]

Hi,

Below, I set HGVS_SEQREPO_DIR environment variable to my local SeqRepo and set url to my local UTA.

import os
os.environ["HGVS_SEQREPO_DIR"] = "/usr/local/share/seqrepo/2024-12-20"

from hgvs.dataproviders.uta import connect
from hgvs.parser import Parser
from hgvs.assemblymapper import AssemblyMapper
from hgvs.normalizer import Normalizer

url = "postgresql://anonymous@localhost:5432/uta/uta_20241220"
hdp = connect(url)
am = AssemblyMapper(
    hdp,
    assembly_name='GRCh37',
    alt_aln_method="splign",
    replace_reference=True
)

hp = Parser()
gdot = hp.parse_hgvs_variant('NC_000001.10:g.150550916_150550920delGACAAinsCAATACC')
norm = Normalizer(hdp, alt_aln_method="splign", validate=False)
res = norm.normalize(gdot)

and can see the following in the logs:

DEBUG:hgvs.dataproviders.uta:connecting to postgresql://anonymous@localhost:5432/uta/uta_20241220...
INFO:biocommons.seqrepo.fastadir.fastadir:File descriptor caching disabled
INFO:hgvs.dataproviders.seqfetcher:Fetching sequences with SeqRepo (/usr/local/share/seqrepo/2024-12-20)
INFO:hgvs.dataproviders.uta:connected to postgresql://anonymous@localhost:5432/uta/uta_20241220...
DEBUG:hgvs.parser:Exposed 33 rules (c_hgvs_position, c_interval, c_pos, c_posedit, c_variant, g_hgvs_position, g_interval, g_pos, g_posedit, g_variant, hgvs_position, hgvs_variant, m_hgvs_position, m_interval, m_pos, m_posedit, m_variant, n_hgvs_position, n_interval, n_pos, n_posedit, n_variant, p_hgvs_position, p_interval, p_pos, p_posedit, p_posedit_special, p_variant, r_hgvs_position, r_interval, r_pos, r_posedit, r_variant)
DEBUG:biocommons.seqrepo.seqaliasdb.seqaliasdb:Executing: select seqalias_id, seq_id, alias, added, is_current, namespace from seqalias where (alias = ?) and (is_current = 1) order by seq_id, namespace, alias with params ['NC_000001.10']
DEBUG:biocommons.seqrepo.fastadir.fastadir:Opening for reading uncached: [/usr/local/share/seqrepo/2024-12-20/sequences/2016/0824/051321/1472015601.985206.fa.bgz](https://file+.vscode-resource.vscode-cdn.net/usr/local/share/seqrepo/2024-12-20/sequences/2016/0824/051321/1472015601.985206.fa.bgz)
DEBUG:biocommons.seqrepo.fastadir.fastadir:Opening for reading uncached: [/usr/local/share/seqrepo/2024-12-20/sequences/2016/0824/051321/1472015601.985206.fa.bgz](https://file+.vscode-resource.vscode-cdn.net/usr/local/share/seqrepo/2024-12-20/sequences/2016/0824/051321/1472015601.985206.fa.bgz)

So this approach does seem to be connecting to my local UTA and SeqRepo.

Would you be able to provide your logs?

[sodiumnitrate]

Hi,

Thank you so much for providing the code. This seems to try to fetch transcript info from seqrepo. I set mine up with only two chromosomes and no transcripts, so fails with an hgvs.exceptions.HGVSDataNotAvailableError.

I added the following line:

cdot = am.g_to_c(gdot, transcript_accession)

and it gives the error. If you set up the seqrepo the default way (i.e. pulling all the data from the default repository), I suspect this should not produce an error on your end. The transcript info would be in seqrepo.

Your logs only show connection messages to UTA -- I'm not sure if an explicit fetching message should also appear.

[korikuzma]

@sodiumnitrate Correct, when I add:

transcript_accession = am.relevant_transcripts(gdot)[0]
cdot = am.g_to_c(gdot, transcript_accession)
cdot

I get the following:

SequenceVariant(ac=NM_001197320.2, type=c, posedit=277_281delinsGGTATTG, gene=None)