Skip to content

Commit 88db5cd

Browse files
andreasprlicclaude
andreasprlic
and
claude
committed
Backport selenocysteine support (#420/#710) from main to 1.5.x
Cherry-pick of merge commit 8a515ce (PR #814) with conflict resolution. Also adds src/hgvs/utils/position.py (needed by variantmapper imports). Key functional change: c_to_p now uses reference_data.translation_table (auto-detected) instead of the caller-supplied translation_table, enabling correct handling of selenocysteine transcripts. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 5f27045 commit 88db5cd

6 files changed

Lines changed: 361 additions & 65 deletions

File tree

‎src/hgvs/utils/position.py‎

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
"""Shared position handling functions for HGVS."""
2+
3+
import hgvs.location
4+
5+
6+
def get_start_end(
7+
var, outer_confidence=True
8+
) -> tuple[
9+
hgvs.location.SimplePosition | hgvs.location.BaseOffsetPosition,
10+
hgvs.location.SimplePosition | hgvs.location.BaseOffsetPosition,
11+
]:
12+
"""Get start and end positions from a variant or interval.
13+
14+
This function handles all position types (SimplePosition, BaseOffsetPosition,
15+
Interval, BaseOffsetInterval) and returns the appropriate start and end positions.
16+
It can be expected that the returned positions have a base and an uncertain property.
17+
18+
By default we return the outer confidence positions. However, if that position
19+
does not have a base, we return the inner confidence positions.
20+
21+
TODO: add a new optional parameter that allows to define the strictness of the returned positions.
22+
The current behavior is more alike to an "auto" mode, since we might fall back to the inner confidence positions
23+
if the outer confidence positions do not have a base. A potential "strict" mode would only return the outer confidence positions,
24+
and raise an error if the outer confidence positions do not have a base.
25+
26+
Args:
27+
var: A variant object with posedit.pos attribute, or an Interval object
28+
outer_confidence: If True, return the outer confidence positions, otherwise return the inner confidence positions
29+
30+
Returns:
31+
tuple: (start_position, end_position) where positions can be SimplePosition or BaseOffsetPosition
32+
"""
33+
34+
# Handle Interval objects directly
35+
if isinstance(var, hgvs.location.Interval):
36+
s = var.start
37+
e = var.end
38+
if not isinstance(
39+
s, (hgvs.location.SimplePosition, hgvs.location.BaseOffsetPosition)
40+
):
41+
if outer_confidence and not s.start.base is not None:
42+
s = s.start
43+
else:
44+
s = s.end
45+
if not isinstance(
46+
e, (hgvs.location.SimplePosition, hgvs.location.BaseOffsetPosition)
47+
):
48+
if outer_confidence and not e.end.base is not None:
49+
e = e.end
50+
else:
51+
e = e.start
52+
return s, e
53+
54+
# Handle position objects directly
55+
if isinstance(
56+
var, (hgvs.location.SimplePosition, hgvs.location.BaseOffsetPosition)
57+
):
58+
return var, var
59+
60+
# if there is no posedit, return None, all steps below this would fail
61+
if var.posedit is None:
62+
return None, None
63+
64+
# Handle variants with posedit
65+
pos = var.posedit.pos
66+
67+
if isinstance(pos, hgvs.location.SimplePosition):
68+
return pos, pos
69+
elif isinstance(pos, hgvs.location.BaseOffsetInterval):
70+
return pos.start, pos.end
71+
elif isinstance(pos, hgvs.location.Interval):
72+
s = pos.start
73+
e = pos.end
74+
75+
if isinstance(s, hgvs.location.AAPosition) and isinstance(
76+
e, hgvs.location.AAPosition
77+
):
78+
s = s.base
79+
e = e.base
80+
return s, e
81+
82+
if not isinstance(
83+
s, (hgvs.location.SimplePosition, hgvs.location.BaseOffsetPosition)
84+
):
85+
orig_s = s
86+
87+
if outer_confidence and s.start.base is not None:
88+
s = s.start
89+
else:
90+
s = s.end
91+
if s.is_uncertain:
92+
s = orig_s.end
93+
94+
if not isinstance(
95+
e, (hgvs.location.SimplePosition, hgvs.location.BaseOffsetPosition)
96+
):
97+
orig_e = e
98+
if outer_confidence and e.end.base is not None:
99+
e = e.end
100+
else:
101+
e = e.start
102+
if e.is_uncertain:
103+
e = orig_e.start
104+
return s, e
105+
else: # BaseOffsetPosition
106+
return pos, pos
107+
108+
109+
def get_start_end_interbase(
110+
pos: hgvs.location.BaseOffsetPosition | hgvs.location.Interval,
111+
outer_confidence=True,
112+
) -> tuple[int, int]:
113+
"""Get start and end integer positions from a SequenceVariant.
114+
115+
This function extracts integer positions from a SequenceVariant, handling uncertain positions
116+
and intervals. For uncertain positions, it uses the more confident boundary.
117+
118+
Args:
119+
var: A SequenceVariant object
120+
outer_confidence: If True, return the outer confidence positions, otherwise return the inner confidence positions
121+
122+
Returns:
123+
tuple: (start_int, end_int) where both are integer positions (0-based)
124+
"""
125+
126+
# Handle start position
127+
if pos.start.uncertain:
128+
# For uncertain start, use the more confident boundary
129+
if isinstance(pos.start, hgvs.location.Interval):
130+
if outer_confidence and pos.start.start.base is not None:
131+
seq_start = pos.start.start.base - 1
132+
else:
133+
seq_start = pos.start.end.base - 1
134+
else:
135+
seq_start = pos.start.base - 1
136+
else:
137+
# For certain start, use the start boundary
138+
if isinstance(pos.start, hgvs.location.Interval):
139+
if outer_confidence:
140+
seq_start = pos.start.start.base - 1
141+
else:
142+
seq_start = pos.start.end.base - 1
143+
else:
144+
seq_start = pos.start.base - 1
145+
146+
# Handle end position
147+
if pos.end.uncertain:
148+
# For uncertain end, use the more confident boundary
149+
if isinstance(pos.end, hgvs.location.Interval):
150+
if outer_confidence and pos.end.end.base is not None:
151+
seq_end = pos.end.end.base
152+
else:
153+
seq_end = pos.end.start.base
154+
else:
155+
seq_end = pos.end.base
156+
else:
157+
# For certain end, use the end boundary
158+
if isinstance(pos.end, hgvs.location.Interval):
159+
if outer_confidence:
160+
seq_end = pos.end.end.base
161+
else:
162+
seq_end = pos.end.start.base
163+
else:
164+
seq_end = pos.end.base
165+
166+
return seq_start, seq_end

‎src/hgvs/utils/reftranscriptdata.py‎

Lines changed: 30 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,14 @@
1-
from bioutils.sequences import translate_cds, TranslationTable
1+
import logging
2+
3+
from bioutils.sequences import TranslationTable, translate_cds
24

35
from hgvs.exceptions import HGVSDataNotAvailableError
46

7+
# Constants for mitochondrial transcript handling
8+
_REMAINDER_TWO = 2
9+
10+
_logger = logging.getLogger(__name__)
11+
512

613
class RefTranscriptData:
714
def __init__(self, hdp, tx_ac, pro_ac, translation_table=TranslationTable.standard):
@@ -10,9 +17,8 @@ def __init__(self, hdp, tx_ac, pro_ac, translation_table=TranslationTable.standa
1017
tx_seq = hdp.get_seq(tx_ac)
1118

1219
if tx_info is None or tx_seq is None:
13-
raise HGVSDataNotAvailableError(
14-
"Missing transcript data for accession: {}".format(tx_ac)
15-
)
20+
msg = f"Missing transcript data for accession: {tx_ac}"
21+
raise HGVSDataNotAvailableError(msg)
1622

1723
# use 1-based hgvs coords
1824
cds_start = tx_info["cds_start_i"] + 1
@@ -24,16 +30,13 @@ def __init__(self, hdp, tx_ac, pro_ac, translation_table=TranslationTable.standa
2430
if len(tx_seq_to_translate) % 3 == 1 and tx_seq_to_translate[-1] == "T":
2531
tx_seq = tx_seq[:cds_stop] + "AA" + tx_seq[cds_stop:]
2632
tx_seq_to_translate += "AA"
27-
if len(tx_seq_to_translate) % 3 == 2 and tx_seq_to_translate[-2:] == "TA":
33+
if len(tx_seq_to_translate) % 3 == _REMAINDER_TWO and tx_seq_to_translate[-2:] == "TA":
2834
tx_seq = tx_seq[:cds_stop] + "A" + tx_seq[cds_stop:]
2935
tx_seq_to_translate += "A"
3036
# coding sequences that are not divisable by 3 are not yet supported
3137
if len(tx_seq_to_translate) % 3 != 0:
32-
raise NotImplementedError(
33-
"Transcript {} is not supported because its sequence length of {} is not divisible by 3.".format(
34-
tx_ac, len(tx_seq_to_translate)
35-
)
36-
)
38+
msg = f"Transcript {tx_ac} is not supported because its sequence length of {len(tx_seq_to_translate)} is not divisible by 3."
39+
raise NotImplementedError(msg)
3740

3841
protein_seq = translate_cds(tx_seq_to_translate, translation_table=translation_table)
3942

@@ -42,6 +45,22 @@ def __init__(self, hdp, tx_ac, pro_ac, translation_table=TranslationTable.standa
4245
# TODO: drop get_acs_for_protein_seq; use known mapping or digest (wo/pro ac inference)
4346
pro_ac = hdp.get_pro_ac_for_tx_ac(tx_ac) or hdp.get_acs_for_protein_seq(protein_seq)[0]
4447

48+
# Auto-detect selenocysteine: check if reference protein sequence contains 'U'
49+
# If it does, re-translate using the selenocysteine translation table
50+
actual_translation_table = translation_table
51+
try:
52+
ref_protein_seq = hdp.get_seq(pro_ac) if pro_ac else None
53+
if ref_protein_seq and "U" in ref_protein_seq:
54+
# Re-translate with selenocysteine table
55+
actual_translation_table = TranslationTable.selenocysteine
56+
protein_seq = translate_cds(
57+
tx_seq_to_translate, translation_table=actual_translation_table
58+
)
59+
except Exception as e:
60+
# If we can't get reference sequence, use what we have
61+
# This is expected for some accessions that may not be available
62+
_logger.debug("Could not fetch reference protein sequence for %s: %s", pro_ac, e)
63+
4564
exon_start_positions = [-tx_info["cds_start_i"]]
4665
exon_end_positions = [exon_start_positions[0] + tx_info["lengths"][0]]
4766
for exon_length in tx_info["lengths"][1:]:
@@ -53,5 +72,6 @@ def __init__(self, hdp, tx_ac, pro_ac, translation_table=TranslationTable.standa
5372
self.cds_start = cds_start
5473
self.cds_stop = cds_stop
5574
self.protein_accession = pro_ac
75+
self.translation_table = actual_translation_table
5676
self.exon_start_positions = exon_start_positions
5777
self.exon_end_positions = exon_end_positions

0 commit comments

Comments
 (0)