From 327d475e665dd5a86f090de905c1d7b33cf55ab6 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 16:10:31 -0400 Subject: [PATCH 1/5] Fix CI: don't require uninstalled Ensembl release 115 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CI (and any env without release 115) was failing three tests because they depend on Ensembl 115, which the openvax/ensembl-data mirror doesn't publish (it tops out at GRCh38.95): * test_frameshift_cterminus_regression.py hardcodes cached_release(115) for the ATM p.F61fs regression from #396. Guard it: skip cleanly at module level when release 115 isn't installed. The same bug CLASS is covered on release 81 by test_annotator_divergence_scenarios.py and test_protein_diff_parity.py, so CI coverage of #396/#397 is unaffected. * test_timings.py used random_variants(), whose genome_name='GRCh38' default resolves via genome_for_reference_name to pyensembl's LATEST release (now 115) rather than an installed one. Add a backward-compatible ensembl= param to random_variants and pin the timing test to cached_release(81). This unbreaks main, which went red at dead346 (#396 merge) — the failures were masked locally on machines that happen to have release 115 cached. Claude-Session: https://claude.ai/code/session_01VNtZRyKZ7u9jiMGPQEbx4c --- tests/test_frameshift_cterminus_regression.py | 17 +++++++++++++++++ tests/test_timings.py | 14 ++++++++++++-- varcode/util.py | 17 +++++++++++++++-- 3 files changed, 44 insertions(+), 4 deletions(-) diff --git a/tests/test_frameshift_cterminus_regression.py b/tests/test_frameshift_cterminus_regression.py index 8bab9d2..0136462 100644 --- a/tests/test_frameshift_cterminus_regression.py +++ b/tests/test_frameshift_cterminus_regression.py @@ -25,6 +25,7 @@ yielded 73 aa ending ``...F-R-K-K-Q-N``. """ +import pytest from pyensembl import cached_release from varcode import Variant @@ -32,6 +33,22 @@ ensembl_grch38 = cached_release(115) +# This regression is pinned against the exact reported variant (ATM p.F61fs +# on MANE Select ENST00000675843), which only exists in newer Ensembl +# releases. Skip cleanly when release 115 isn't installed -- CI's data mirror +# (openvax/ensembl-data) tops out at GRCh38.95, so 115 is unavailable there. +# The same bug CLASS is exercised on an installed release (81) by +# tests/test_annotator_divergence_scenarios.py (CFTR p.L127fs, BRCA1 p.R71fs) +# and tests/test_protein_diff_parity.py, so CI coverage of #396/#397 does not +# depend on release 115. +try: + ensembl_grch38.transcript_by_id("ENST00000675843") +except Exception as _exc: # pyensembl raises if the GTF DB isn't downloaded + pytest.skip( + "Ensembl release 115 not installed (%s); ATM regression covered on " + "release 81 elsewhere." % type(_exc).__name__, + allow_module_level=True) + def _atm_f61fs_effect(annotator): # VCF-style anchored representation of chr11:g.108227882delT diff --git a/tests/test_timings.py b/tests/test_timings.py index 3666889..fc8f1b1 100644 --- a/tests/test_timings.py +++ b/tests/test_timings.py @@ -13,8 +13,16 @@ from __future__ import print_function, division, absolute_import import time +from pyensembl import cached_release + from varcode.util import random_variants +# Pin to an installed release. genome_for_reference_name("GRCh38") (the +# random_variants default) resolves to pyensembl's latest release, which is +# not necessarily downloaded (CI installs 75/81/95 only); release 81 matches +# the rest of the suite. +ensembl_grch38 = cached_release(81) + def _time_variant_annotation(variant_collection): start_t = time.time() effects = variant_collection.effects() @@ -30,12 +38,14 @@ def test_effect_timing( n_warmup_variants=5): warmup_collection = random_variants( n_warmup_variants, - random_seed=None) + random_seed=None, + ensembl=ensembl_grch38) warmup_collection.effects() variant_collection = random_variants( n_variants, - random_seed=random_seed) + random_seed=random_seed, + ensembl=ensembl_grch38) elapsed_t = _time_variant_annotation(variant_collection) print("Elapsed: %0.4f for %d variants" % (elapsed_t, n_variants)) assert elapsed_t / n_variants < 0.1, \ diff --git a/varcode/util.py b/varcode/util.py index 245b604..c2b771e 100644 --- a/varcode/util.py +++ b/varcode/util.py @@ -26,13 +26,26 @@ def random_variants( genome_name="GRCh38", deletions=True, insertions=True, - random_seed=None): + random_seed=None, + ensembl=None): """ Generate a VariantCollection with random variants that overlap at least one complete coding transcript. + + Parameters + ---------- + ensembl : pyensembl.EnsemblRelease, optional + Explicit genome to draw transcripts from. When ``None`` (default) + the genome is resolved from ``genome_name`` via + :func:`genome_for_reference_name`, which picks pyensembl's *latest* + release for that assembly. Pass a pinned release (e.g. + ``cached_release(81)``) when the caller needs deterministic data + that is actually installed — ``genome_for_reference_name`` can + resolve to a release that hasn't been downloaded. """ rng = random.Random(random_seed) - ensembl = genome_for_reference_name(genome_name) + if ensembl is None: + ensembl = genome_for_reference_name(genome_name) if ensembl in _transcript_ids_cache: transcript_ids = _transcript_ids_cache[ensembl] From 5e639fc1c1421cfdacae1d7bfcd3ecb00d72b697 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 16:29:07 -0400 Subject: [PATCH 2/5] random_variants: skip transcripts on contigs Variant rejects Pinning test_timings to release 81 surfaced a latent bug: random_variants draws from all of ensembl.transcript_ids(), which on release 81 includes transcripts on alternate/patch scaffolds (e.g. CHR_HSCHR19LRC_COX1_CTG3_1) that Variant() rejects as non-standard for GRCh38, raising ValueError and failing the whole generator. It was flaky rather than deterministic because the timing test's warmup draw is unseeded. Wrap the per-draw body in try/except ValueError and skip to the next transcript, so an unlucky pick no longer fails generation. The count*100 iteration budget and the terminal 'Unable to generate' error still catch a systematically-broken genome. Verified: 40 seeds x 50 variants on release 81 now all succeed. Claude-Session: https://claude.ai/code/session_01VNtZRyKZ7u9jiMGPQEbx4c --- varcode/util.py | 59 +++++++++++++++++++++++++++++-------------------- 1 file changed, 35 insertions(+), 24 deletions(-) diff --git a/varcode/util.py b/varcode/util.py index c2b771e..cf084eb 100644 --- a/varcode/util.py +++ b/varcode/util.py @@ -65,33 +65,44 @@ def random_variants( if not transcript.complete: continue - exon = rng.choice(transcript.exons) - base1_genomic_position = rng.randint(exon.start, exon.end) - transcript_offset = transcript.spliced_offset(base1_genomic_position) - seq = transcript.sequence + try: + exon = rng.choice(transcript.exons) + base1_genomic_position = rng.randint(exon.start, exon.end) + transcript_offset = transcript.spliced_offset( + base1_genomic_position) + seq = transcript.sequence - ref = str(seq[transcript_offset]) - if transcript.on_backward_strand: - ref = reverse_complement(ref) + ref = str(seq[transcript_offset]) + if transcript.on_backward_strand: + ref = reverse_complement(ref) - alt_nucleotides = [x for x in STANDARD_NUCLEOTIDES if x != ref] + alt_nucleotides = [x for x in STANDARD_NUCLEOTIDES if x != ref] - if insertions: - nucleotide_pairs = [ - x + y - for x in STANDARD_NUCLEOTIDES - for y in STANDARD_NUCLEOTIDES - ] - alt_nucleotides.extend(nucleotide_pairs) - if deletions: - alt_nucleotides.append("") - alt = rng.choice(alt_nucleotides) - variant = Variant( - transcript.contig, - base1_genomic_position, - ref=ref, - alt=alt, - ensembl=ensembl) + if insertions: + nucleotide_pairs = [ + x + y + for x in STANDARD_NUCLEOTIDES + for y in STANDARD_NUCLEOTIDES + ] + alt_nucleotides.extend(nucleotide_pairs) + if deletions: + alt_nucleotides.append("") + alt = rng.choice(alt_nucleotides) + variant = Variant( + transcript.contig, + base1_genomic_position, + ref=ref, + alt=alt, + ensembl=ensembl) + except ValueError: + # Some transcripts live on alternate/patch contigs (e.g. + # 'CHR_HSCHR19LRC_COX1_CTG3_1') that Variant rejects as + # non-standard for the reference, and a few have sequence / + # offset edge cases. Skip and draw another transcript rather + # than failing the whole generator on an unlucky pick — this + # otherwise made the result depend on the (often unseeded) + # draw order. + continue variants.append(variant) else: return VariantCollection(variants) From fc6c1f20445877afb70096f3695dfd81bf72eb45 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 16:40:14 -0400 Subject: [PATCH 3/5] random_variants: filter transcripts to valid contigs (the real fix) The previous commit wrapped Variant construction in try/except, but the 'Invalid contig name' ValueError is raised LAZILY -- from the .transcripts / .genes properties via _check_that_genome_has_contig, not from __init__ -- so construction succeeded and the error only surfaced later in .effects(). That made the earlier fix (and its stress test, which never called .effects()) ineffective; CI kept failing on a different alt scaffold (CHR_HSCHR11_2_CTG1). Filter transcripts up front to those whose contig is in set(ensembl.contigs()) -- the exact validity set _check_that_genome_has_contig uses -- so every returned variant is annotatable. Variant preserves the scaffold contig name (verified), so the filter aligns exactly with the validator. Note contigs() differs by environment (locally it includes alt scaffolds, CI excludes them); because filter and validator both call contigs(), they stay consistent either way. Verified under a simulated-CI contigs() (272 scaffolds excluded): 25 seeds x 30 variants generate + annotate with zero failures and no leaked contigs. Claude-Session: https://claude.ai/code/session_01VNtZRyKZ7u9jiMGPQEbx4c --- varcode/util.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/varcode/util.py b/varcode/util.py index cf084eb..f5b1696 100644 --- a/varcode/util.py +++ b/varcode/util.py @@ -53,6 +53,16 @@ def random_variants( transcript_ids = ensembl.transcript_ids() _transcript_ids_cache[ensembl] = transcript_ids + # Only draw from transcripts on contigs the genome considers valid. + # ``transcript_ids()`` includes transcripts on alternate/patch scaffolds + # (e.g. 'CHR_HSCHR11_2_CTG1') whose contig is NOT in ``genome.contigs()``; + # a Variant built there passes construction but raises + # ``ValueError: Invalid contig name`` lazily, when effect prediction + # accesses ``.transcripts`` (see Variant._check_that_genome_has_contig). + # Mirror that validity set here so we never hand back a variant that + # can't be annotated. + valid_contigs = set(ensembl.contigs()) + variants = [] # we should finish way before this loop is over but just in case @@ -65,6 +75,11 @@ def random_variants( if not transcript.complete: continue + if transcript.contig not in valid_contigs: + # Alternate/patch scaffold — Variant would reject this contig + # during annotation. Skip and draw another. + continue + try: exon = rng.choice(transcript.exons) base1_genomic_position = rng.randint(exon.start, exon.end) From 0e14718d03e2c1365a76ce93263340b3b0016c90 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 16:51:49 -0400 Subject: [PATCH 4/5] random_variants: validate contig against the variant's OWN genome Third time's the charm. The prior filter used set(ensembl.contigs()), but the lazy validator (_check_that_genome_has_contig) checks self.genome.contigs(), and self.genome = infer_genome(ensembl) is NOT guaranteed to be the same object as the passed ensembl -- in CI their contig sets differed, so the filter excluded nothing and CI kept failing on CHR_HSCHR11_2_CTG1. Validate each built variant against variant.genome.contigs() -- the exact same object and call the validator uses -- so the check is consistent with effect prediction by construction, regardless of how infer_genome resolves ensembl or how contigs() varies across environments. valid_contigs is computed once from the first successfully-built variant (all share a genome). Verified by patching the resolved genome's contigs() to exclude its 272 alt scaffolds (simulating CI): 25 seeds x 30 variants now generate AND annotate with zero failures and no leaked contigs. Claude-Session: https://claude.ai/code/session_01VNtZRyKZ7u9jiMGPQEbx4c --- varcode/util.py | 42 ++++++++++++++++++++---------------------- 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/varcode/util.py b/varcode/util.py index f5b1696..d711028 100644 --- a/varcode/util.py +++ b/varcode/util.py @@ -53,18 +53,19 @@ def random_variants( transcript_ids = ensembl.transcript_ids() _transcript_ids_cache[ensembl] = transcript_ids - # Only draw from transcripts on contigs the genome considers valid. - # ``transcript_ids()`` includes transcripts on alternate/patch scaffolds - # (e.g. 'CHR_HSCHR11_2_CTG1') whose contig is NOT in ``genome.contigs()``; - # a Variant built there passes construction but raises - # ``ValueError: Invalid contig name`` lazily, when effect prediction - # accesses ``.transcripts`` (see Variant._check_that_genome_has_contig). - # Mirror that validity set here so we never hand back a variant that - # can't be annotated. - valid_contigs = set(ensembl.contigs()) - variants = [] + # Effect prediction validates a variant's contig against its OWN + # ``self.genome.contigs()`` (see Variant._check_that_genome_has_contig), + # where ``self.genome = infer_genome(ensembl)`` -- NOT necessarily the + # ``ensembl`` object passed here, and whose contig set can exclude + # alternate/patch scaffolds (e.g. 'CHR_HSCHR11_2_CTG1'). A variant built on + # such a scaffold passes construction but raises ``ValueError: Invalid + # contig name`` lazily, when it's annotated. Validate against the variant's + # own genome so every returned variant is annotatable. Computed once from + # the first successfully-built variant (all share the same genome). + valid_contigs = None + # we should finish way before this loop is over but just in case # something is wrong with PyEnsembl we want to avoid an infinite loop for _ in range(count * 100): @@ -75,11 +76,6 @@ def random_variants( if not transcript.complete: continue - if transcript.contig not in valid_contigs: - # Alternate/patch scaffold — Variant would reject this contig - # during annotation. Skip and draw another. - continue - try: exon = rng.choice(transcript.exons) base1_genomic_position = rng.randint(exon.start, exon.end) @@ -110,13 +106,15 @@ def random_variants( alt=alt, ensembl=ensembl) except ValueError: - # Some transcripts live on alternate/patch contigs (e.g. - # 'CHR_HSCHR19LRC_COX1_CTG3_1') that Variant rejects as - # non-standard for the reference, and a few have sequence / - # offset edge cases. Skip and draw another transcript rather - # than failing the whole generator on an unlucky pick — this - # otherwise made the result depend on the (often unseeded) - # draw order. + # A few transcripts have sequence/offset edge cases that make + # Variant construction itself raise; skip and draw another. + continue + if valid_contigs is None: + valid_contigs = set(variant.genome.contigs()) + if variant.contig not in valid_contigs: + # Alternate/patch scaffold whose contig isn't in the genome's + # valid set — annotating it would raise ValueError lazily. + # This is what previously made the (unseeded) timing test flaky. continue variants.append(variant) else: From 9cc481dee300cfb302514a639eed201f8b7d9ba0 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Wed, 8 Jul 2026 17:00:59 -0400 Subject: [PATCH 5/5] random_variants: validate via variant.transcripts (shared-cache exact path) My prior two attempts compared the variant's contig against a locally computed contig set (ensembl.contigs(), then variant.genome.contigs()). Both were inconsistent with the validator, because Variant._check_that_genome_has_contig reads a PROCESS-WIDE cache (_reference_name_to_valid_contig_names, keyed by reference name) that an earlier test populates with a scaffold-free set. So my set could include scaffolds the cached validator set excluded, and CI kept failing on a different alt scaffold each unseeded run. Skip via the exact path the caller uses: after constructing the variant, access variant.transcripts (which runs _check_that_genome_has_contig against the shared cache) inside the try/except. A scaffold variant raises there and is skipped; what survives is guaranteed annotatable by the same cache .effects() will read. Also drop variants that overlap no transcript. Verified by seeding _reference_name_to_valid_contig_names['GRCh38'] with a scaffold-free set (as CI's earlier tests do) and generating + annotating across 20 seeded and 10 unseeded draws: zero failures, no leaked scaffolds. Claude-Session: https://claude.ai/code/session_01VNtZRyKZ7u9jiMGPQEbx4c --- varcode/util.py | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/varcode/util.py b/varcode/util.py index d711028..7521bad 100644 --- a/varcode/util.py +++ b/varcode/util.py @@ -55,17 +55,6 @@ def random_variants( variants = [] - # Effect prediction validates a variant's contig against its OWN - # ``self.genome.contigs()`` (see Variant._check_that_genome_has_contig), - # where ``self.genome = infer_genome(ensembl)`` -- NOT necessarily the - # ``ensembl`` object passed here, and whose contig set can exclude - # alternate/patch scaffolds (e.g. 'CHR_HSCHR11_2_CTG1'). A variant built on - # such a scaffold passes construction but raises ``ValueError: Invalid - # contig name`` lazily, when it's annotated. Validate against the variant's - # own genome so every returned variant is annotatable. Computed once from - # the first successfully-built variant (all share the same genome). - valid_contigs = None - # we should finish way before this loop is over but just in case # something is wrong with PyEnsembl we want to avoid an infinite loop for _ in range(count * 100): @@ -105,16 +94,25 @@ def random_variants( ref=ref, alt=alt, ensembl=ensembl) + # Force the lazy contig validation NOW, through the exact path + # effect prediction uses: Variant.transcripts calls + # _check_that_genome_has_contig, which reads a process-wide + # valid-contig cache keyed by reference name. A variant built on + # an alternate/patch scaffold (e.g. 'CHR_HSCHR6_MHC_MCF_CTG1') + # constructs fine but raises here; skipping now guarantees we + # never return a variant that blows up a later .effects(). + # Going through .transcripts -- not a separately computed contig + # set -- is what makes this consistent with the caller (earlier + # attempts compared against the wrong set and let scaffolds + # through). + overlapping = variant.transcripts except ValueError: - # A few transcripts have sequence/offset edge cases that make - # Variant construction itself raise; skip and draw another. + # Alternate/patch-scaffold contig, or a transcript with a + # sequence/offset edge case: skip and draw another rather than + # failing the whole generator on an unlucky (often unseeded) + # pick. continue - if valid_contigs is None: - valid_contigs = set(variant.genome.contigs()) - if variant.contig not in valid_contigs: - # Alternate/patch scaffold whose contig isn't in the genome's - # valid set — annotating it would raise ValueError lazily. - # This is what previously made the (unseeded) timing test flaky. + if not overlapping: continue variants.append(variant) else: