Inherit native reference DNA and preserve genome identity (10.11.0)

CHANGELOG.md

@@ -1,5 +1,18 @@
 # Change Log
 
+## [v10.11.0](https://github.com/openvax/varcode/tree/v10.11.0) (2026-09-29)
+
+- `Genome(native_genome)` inherits optional native PyEnsembl reference DNA
+  instead of hiding it (#488), enabling intronic/intergenic reference lookups.
+  Construction and rewrapping remain lazy; explicit `fasta=` takes precedence.
+- Native `sequence()` preserves PyEnsembl's errors for uninstalled DNA and
+  invalid intervals. Tiered reference lookups retain their cDNA fallback.
+  Reads never download DNA. `close()` closes only explicit path readers opened
+  by that wrapper, leaving native genomes and borrowed readers usable.
+- Preserve the identity of supplied genome objects during inference (#565).
+  Equal annotation genomes can carry different DNA; caching them by equality
+  could silently select the wrong reference file.
+
 ## [v10.10.0](https://github.com/openvax/varcode/tree/v10.10.0) (2026-09-29)
 
 - Add `reference_completion_hypotheses` for mapped Exacto and other RNA

docs/api_variants.md

@@ -5,6 +5,45 @@ For examples, see [file loading](getting_started.md#load-a-file),
 
 ## Reference genomes
 
+`Genome` inherits reference DNA configured on a PyEnsembl genome (PyEnsembl
+2.11 or later). This makes intronic and intergenic bases available to Varcode's
+reference lookups as well as features such as indel left alignment.
+
+```python
+from pyensembl import EnsemblRelease
+from varcode import Genome
+
+native = EnsemblRelease(81, genome_fasta_path="/data/GRCh38.fa")
+genome = Genome(native)
+bases = genome.reference_range("7", 117_480_000, 117_480_050)
+```
+
+Construction and rewrapping do not open or download native DNA. The first
+read can index a local FASTA or decompress it into PyEnsembl's cache. For a
+remote source, install DNA explicitly before querying:
+
+```python
+native = EnsemblRelease(81, download_genome_fasta=True)
+native.download_genome_fasta()  # explicit, potentially large download
+native.index_genome_fasta()
+genome = Genome(native)
+```
+
+An explicit `Genome(native, fasta="/data/other.fa")` overrides native DNA;
+the supplied reference must match the annotation assembly. Without DNA,
+`reference_base` and `reference_range` fall back to transcript cDNA and return
+`""` for uncovered positions. `sequence` reads only chromosome DNA: configured
+native DNA preserves PyEnsembl's missing-DNA and invalid-interval errors,
+while unconfigured genomes return `""`. Explicit FASTA overrides retain their
+legacy permissive lookup behavior. Coordinates are 1-based inclusive, and
+returned bases are uppercase on the genomic plus strand.
+
+The native PyEnsembl genome owns its reader. Calling `genome.close()` leaves
+native DNA and caller-provided reader objects open; it closes only readers
+that this wrapper opened from an explicit `fasta=` path. Rewrapping borrows an
+explicit reader, so its owning wrapper must stay open. Close the native genome
+or a caller-provided reader only after all its borrowers have finished.
+
 ::: varcode.Genome
 
 ## Variants

docs/transforms.md

@@ -138,7 +138,7 @@ leftmost equivalent position.
 from varcode import load_vcf, Genome
 from varcode.transforms import left_align_indels
 
-# Default (transcript-cDNA coverage only) — exonic indels normalize,
+# Without installed reference DNA — exonic indels normalize,
 # intronic/intergenic pass through unchanged.
 vc = load_vcf("tumor.vcf", genome="GRCh38")
 vc = left_align_indels(vc)
@@ -153,6 +153,8 @@ No `reference` parameter — `left_align_indels` reads bases via the
 genome the variants already carry (see
 [varcode.Genome](api_variants.md#varcode.Genome)). Coverage depends on which
 genome shape was passed.
+Native PyEnsembl reference DNA is also inherited when wrapping a configured
+genome; see [reference genomes](api_variants.md#reference-genomes) for setup.
 
 ### Behavior
 

tests/test_genome_sequence.py

@@ -31,6 +31,7 @@
 """
 
 import copy
+import pickle
 import warnings
 
 import pytest
@@ -152,6 +153,25 @@ def test_genome_rewrap_overrides_fasta(ensembl, cftr_position):
     assert g2.fasta is fasta2
 
 
+@pytest.mark.parametrize("with_override", [False, True])
+def test_restore_legacy_pickled_genome(ensembl, with_override):
+    # The public Genome's pickle state through 10.10.x stored .fasta
+    # directly. Emulate that old state without the new private fields.
+    legacy = Genome.__new__(Genome)
+    legacy.__dict__.update(
+        _ensembl=ensembl,
+        fasta=_DictBackedFasta({"1": "acgt"}) if with_override else None,
+        _missing_reference_warned=False)
+    restored = pickle.loads(pickle.dumps(legacy))
+    assert restored.reference_name == ensembl.reference_name
+    assert restored.sequence("1", 1, 4) == ("ACGT" if with_override else "")
+    assert "FASTA" in repr(restored)
+    assert Genome(restored).sequence("1", 1, 4) == restored.sequence("1", 1, 4)
+    restored.close()
+    # Legacy explicit readers are borrowed, like caller-provided readers.
+    assert restored.sequence("1", 1, 4) == ("ACGT" if with_override else "")
+
+
 def test_genome_with_integer_fasta_raises_typeerror(ensembl):
     with pytest.raises(TypeError, match="pyfaidx-style protocol"):
         Genome(ensembl, fasta=42, verify=False)

tests/test_native_reference_dna.py

@@ -0,0 +1,206 @@
+"""Native PyEnsembl DNA stays available through varcode.Genome (#488)."""
+
+import inspect
+
+import pytest
+from pyensembl import Genome as AnnotationGenome
+
+from varcode import Genome, Variant
+from varcode.genome_sequence import reference_base, reference_range
+from varcode.reference import infer_genome
+
+
+@pytest.fixture
+def native_genome(tmp_path):
+    if "genome_fasta_path_or_url" not in inspect.signature(AnnotationGenome).parameters:
+        pytest.skip("Native reference DNA requires PyEnsembl >= 2.11")
+    fasta = tmp_path / "reference.fa"
+    fasta.write_text(">1\naacgttcaggtaccgattcgaacgt\n")
+    genome = AnnotationGenome(
+        reference_name="GRCh38", annotation_name="native_dna_test",
+        genome_fasta_path_or_url=str(fasta),
+        cache_directory_path=str(tmp_path / "cache"))
+    yield genome
+    genome.close()
+
+
+def test_native_intronic_and_intergenic_queries(native_genome, tmp_path):
+    gtf = tmp_path / "tiny.gtf"
+    attributes = (
+        'gene_id "g1"; gene_name "G1"; transcript_id "t1"; '
+        'transcript_name "T1"; transcript_biotype "lncRNA";')
+    gtf.write_text("".join(
+        '1\ttest\t%s\t%d\t%d\t.\t-\t.\t%s\n' % (
+            feature, start, end, attributes + extra)
+        for feature, start, end, extra in [
+            ("gene", 2, 13, ""), ("transcript", 2, 13, ""),
+            ("exon", 2, 5, ' exon_id "e1"; exon_number "2";'),
+            ("exon", 10, 13, ' exon_id "e2"; exon_number "1";')]))
+    native = AnnotationGenome(
+        reference_name="GRCh38", annotation_name="native_dna_annotation",
+        gtf_path_or_url=str(gtf),
+        genome_fasta_path_or_url=str(tmp_path / "reference.fa"),
+        cache_directory_path=str(tmp_path / "annotated_cache"))
+    try:
+        native.index()
+        wrapped = Genome(native)
+        # This intron belongs to a minus-strand transcript; genomic reads
+        # must still return plus-strand DNA, with inclusive endpoints.
+        transcript = native.transcript_by_id("t1")
+        assert transcript.contains("1", 6, 9)
+        assert all(not exon.overlaps("1", 6, 9) for exon in transcript.exons)
+        assert not native.transcripts_at_locus("1", 16, 20)
+        assert wrapped.sequence("1", 6, 9) == "TCAG"
+        assert wrapped.reference_range("1", 6, 9) == "TCAG"
+        assert reference_base(wrapped, "1", 7) == "C"
+        assert reference_range(wrapped, "1", 16, 20) == "ATTCG"
+        assert wrapped.fasta is native.fasta
+    finally:
+        native.close()
+
+
+def test_construction_inference_and_repr_do_not_open_native_dna(
+        native_genome, monkeypatch):
+    def unexpected_open(self):
+        pytest.fail("Native FASTA opened before a sequence request")
+
+    monkeypatch.setattr(AnnotationGenome, "fasta", property(unexpected_open))
+    wrapped = Genome(native_genome)
+    rewrapped = Genome(wrapped)
+    assert infer_genome(native_genome) == (native_genome, False)
+    assert infer_genome(wrapped) == (wrapped, False)
+    assert Variant("1", 7, "C", "A", genome=rewrapped).genome is rewrapped
+    assert "native FASTA configured" in repr(wrapped)
+    wrapped.close()
+
+
+def test_rewrap_tracks_native_reader_refresh(native_genome):
+    wrapped = Genome(native_genome)
+    rewrapped = Genome(wrapped)
+    original = wrapped.fasta
+    assert original is not None
+    native_genome.clear_cache()
+    assert wrapped.fasta is not original
+    assert rewrapped.fasta is wrapped.fasta is native_genome.fasta
+    assert str(original["1"][5:9]).upper() == "TCAG"
+    original.close()
+
+
+@pytest.mark.parametrize("rewrap", [False, True])
+def test_closing_borrowing_wrapper_keeps_native_reader_usable(native_genome, rewrap):
+    first = Genome(native_genome)
+    second = Genome(first if rewrap else native_genome)
+    reader = native_genome.fasta
+    first.close()
+    first.close()
+    assert str(reader["1"][5:9]).upper() == "TCAG"
+    assert second.sequence("1", 6, 9) == "TCAG"
+    assert native_genome.fasta is reader
+
+
+def test_explicit_override_and_rewrap_precede_native_dna(native_genome, tmp_path):
+    override = tmp_path / "override.fa"
+    override.write_text(">1\n" + "a" * 24 + "\n")
+    wrapped = Genome(native_genome, fasta=override, verify=False)
+    rewrapped = Genome(wrapped)
+    assert rewrapped.fasta is wrapped.fasta
+    assert rewrapped.sequence("1", 6, 9) == "AAAA"
+    assert rewrapped.reference_range("1", 6, 9) == "AAAA"
+    assert native_genome.sequence("1", 6, 9) == "TCAG"
+    reader = wrapped.fasta
+    rewrapped.close()
+    assert not reader.faidx.file.closed
+    assert wrapped.reference_base("1", 6) == "A"
+    wrapped.close()
+    assert reader.faidx.file.closed
+
+
+def test_caller_provided_reader_is_borrowed(native_genome, tmp_path):
+    from pyfaidx import Fasta
+
+    path = tmp_path / "borrowed.fa"
+    path.write_text(">1\ncccc\n")
+    with Fasta(str(path)) as reader:
+        wrapped = Genome(native_genome, fasta=reader, verify=False)
+        wrapped.close()
+        assert wrapped.sequence("1", 1, 4) == "CCCC"
+        assert not reader.faidx.file.closed
+
+
+def test_assigning_fasta_releases_owned_reader(native_genome, tmp_path):
+    path = tmp_path / "override.fa"
+    path.write_text(">1\ncccc\n")
+    wrapped = Genome(native_genome, fasta=path, verify=False)
+    reader = wrapped.fasta
+    wrapped.fasta = None
+    assert reader.faidx.file.closed
+    assert wrapped.sequence("1", 6, 9) == "TCAG"
+
+
+@pytest.mark.parametrize("remote", [False, True])
+def test_uninstalled_native_dna_never_downloads(native_genome, tmp_path, monkeypatch, remote):
+    from pyensembl.genome_fasta import GenomeFasta, MissingGenomeFastaError
+
+    def unexpected_download(*args, **kwargs):
+        pytest.fail("Reference lookup attempted a download")
+
+    monkeypatch.setattr(GenomeFasta, "_download", unexpected_download)
+    source = ("https://example.invalid/reference.fa" if remote
+              else str(tmp_path / "missing.fa"))
+    native = AnnotationGenome(
+        reference_name="GRCh38", annotation_name="missing_native_dna",
+        genome_fasta_path_or_url=source,
+        cache_directory_path=str(tmp_path / "missing_cache"))
+    monkeypatch.setattr(native, "transcripts_at_locus", lambda *args, **kwargs: [])
+    wrapped = Genome(native)
+    assert wrapped.fasta is None
+    assert wrapped.reference_base("1", 7) == ""
+    assert wrapped.reference_range("1", 6, 9) == ""
+    with pytest.raises(MissingGenomeFastaError):
+        wrapped.sequence("1", 6, 9)
+    wrapped.close()
+
+
+def test_equal_native_genomes_keep_distinct_dna(native_genome, tmp_path):
+    path = tmp_path / "other.fa"
+    path.write_text(">1\n" + "c" * 24 + "\n")
+    other = AnnotationGenome(
+        reference_name=native_genome.reference_name,
+        annotation_name=native_genome.annotation_name,
+        genome_fasta_path_or_url=str(path),
+        cache_directory_path=str(tmp_path / "other_cache"))
+    try:
+        assert native_genome == other
+        assert infer_genome(native_genome)[0] is native_genome
+        assert infer_genome(other)[0] is other
+        assert Genome(native_genome).sequence("1", 6, 9) == "TCAG"
+        assert Genome(other).sequence("1", 6, 9) == "CCCC"
+        assert Variant("1", 7, "C", "A", genome=other).genome is other
+    finally:
+        other.close()
+
+
+@pytest.mark.parametrize("start,end", [(0, 3), (3, 2), (20, 30)])
+def test_native_sequence_preserves_coordinate_errors(native_genome, start, end):
+    with pytest.raises(ValueError):
+        Genome(native_genome).sequence("1", start, end)
+
+
+def test_native_sequence_preserves_missing_contig_error(native_genome):
+    with pytest.raises(ValueError, match="Contig"):
+        Genome(native_genome).sequence("missing", 1, 4)
+
+
+def test_old_pyensembl_without_native_api(monkeypatch, tmp_path):
+    # Simulate supported PyEnsembl versions before native DNA was added.
+    monkeypatch.delattr(AnnotationGenome, "fasta", raising=False)
+    monkeypatch.delattr(AnnotationGenome, "requires_genome_fasta", raising=False)
+    native = AnnotationGenome(
+        reference_name="GRCh38", annotation_name="old_api",
+        cache_directory_path=str(tmp_path / "old_cache"))
+    monkeypatch.setattr(native, "transcripts_at_locus", lambda *args, **kwargs: [])
+    wrapped = Genome(native)
+    assert wrapped.fasta is None
+    assert wrapped.sequence("1", 1, 4) == ""
+    assert wrapped.reference_base("1", 1) == ""
+    wrapped.close()