From 4bd835d6c25d12ad34ccefdc4848204f757829dc Mon Sep 17 00:00:00 2001
From: Alex Rubinsteyn <alex.rubinsteyn@gmail.com>
Date: Wed, 8 Jul 2026 22:05:01 -0400
Subject: [PATCH] Add selenocysteine (Sec/U) regression tests (#42)

#42 (2015) observed varcode treating selenocysteines as stop codons,
producing spurious StopLoss / UTR effects for selenoproteins. That no
longer happens: varcode annotates against Ensembl's reference
protein_sequence (full-length, with U at the recoded UGA) rather than
re-translating the CDS with the standard table.

Verified on GPX1 (ENST00000419783, Sec at residue 49) and locked in with
regression tests covering the three paths:
- substitution AT the Sec codon -> Substitution p.U49C (aa_ref="U"),
  not StopLoss/UTR;
- substitution downstream -> p.R100W with correct numbering (read-through);
- in-frame deletion upstream -> full-length mutant protein that retains U.

Passes under both the fast and protein_diff annotators.

Closes #42
---
 tests/test_selenocysteine.py | 89 ++++++++++++++++++++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 tests/test_selenocysteine.py

diff --git a/tests/test_selenocysteine.py b/tests/test_selenocysteine.py
new file mode 100644
index 0000000..9923cce
--- /dev/null
+++ b/tests/test_selenocysteine.py
@@ -0,0 +1,89 @@
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#         http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Regression tests for selenocysteine (Sec / U) handling (issue #42).
+
+Selenoproteins recode an in-frame UGA -- normally a stop codon -- as
+selenocysteine. Ensembl's reference protein for such transcripts is
+full-length with a ``U`` at the recoded position, and varcode annotates
+against that reference protein_sequence rather than re-translating the
+coding sequence with the standard table. As a result an internal UGA is
+NOT mistaken for a stop, so selenoprotein variants are no longer
+mis-annotated as StopLoss / UTR effects (the original bug in #42).
+
+GPX1 (glutathione peroxidase 1, ENSG00000233276) is the canonical test
+case: a single selenocysteine at residue 49 (transcript ENST00000419783,
+minus strand). Coordinates are GRCh38 and the transcript is pinned to the
+installed release used across the suite.
+"""
+
+from pyensembl import cached_release
+
+from varcode import Variant
+from varcode.effects import Substitution, Deletion
+
+from .common import expect_effect
+
+ensembl_grch38 = cached_release(81)
+
+# GPX1 selenoprotein transcript: minus strand, selenocysteine at residue 49.
+GPX1_TRANSCRIPT_ID = "ENST00000419783"
+
+
+def _gpx1_variant(contig_pos, ref, alt):
+    return Variant("3", contig_pos, ref, alt, ensembl=ensembl_grch38)
+
+
+def test_substitution_at_selenocysteine_reports_U_not_stop():
+    # The selenocysteine codon itself (UGA) is reported as reference amino
+    # acid ``U``, and a change there is an ordinary missense -- NOT a
+    # StopLoss and NOT a UTR effect, which is what #42 originally observed.
+    variant = _gpx1_variant(49358132, "T", "A")  # UGA -> UGU, Sec -> Cys
+    expect_effect(
+        variant,
+        transcript_id=GPX1_TRANSCRIPT_ID,
+        effect_class=Substitution,
+        aa_ref="U",
+        aa_alt="C",
+        short_description="p.U49C")
+
+
+def test_substitution_downstream_of_selenocysteine_keeps_numbering():
+    # A missense well past the selenocysteine keeps correct residue
+    # numbering (100, not truncated at the internal UGA), confirming the
+    # reference protein reads through the Sec.
+    variant = _gpx1_variant(49357702, "G", "A")
+    expect_effect(
+        variant,
+        transcript_id=GPX1_TRANSCRIPT_ID,
+        effect_class=Substitution,
+        aa_ref="R",
+        aa_alt="W",
+        short_description="p.R100W")
+
+
+def test_inframe_deletion_upstream_of_selenocysteine_reads_through():
+    # An in-frame deletion upstream of the selenocysteine exercises the
+    # re-translation path. The mutant protein must stay full length and
+    # retain the ``U`` rather than truncating at the internal UGA.
+    variant = _gpx1_variant(49358144, "CAC", "")  # deletes one upstream codon
+    transcript = ensembl_grch38.transcript_by_id(GPX1_TRANSCRIPT_ID)
+    effect = variant.effect_on_transcript(transcript)
+    assert isinstance(effect, Deletion), \
+        "Expected an in-frame Deletion, got %s" % type(effect).__name__
+    mutant_protein = effect.mutant_protein_sequence
+    reference_protein = transcript.protein_sequence
+    assert "U" in mutant_protein, \
+        "Selenocysteine dropped from mutant protein: %s" % mutant_protein
+    assert len(mutant_protein) == len(reference_protein) - 1, \
+        "Expected a single-residue deletion, ref=%d mutant=%d" % (
+            len(reference_protein), len(mutant_protein))