Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions tests/test_frameshift_cterminus_regression.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,30 @@
yielded 73 aa ending ``...F-R-K-K-Q-N``.
"""

import pytest
from pyensembl import cached_release

from varcode import Variant
from varcode.effects import FrameShift

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
Expand Down
14 changes: 12 additions & 2 deletions tests/test_timings.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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, \
Expand Down
87 changes: 61 additions & 26 deletions varcode/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -52,33 +65,55 @@ 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)
# 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:
# 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 not overlapping:
continue
variants.append(variant)
else:
return VariantCollection(variants)
Expand Down
Loading