From 7b0e6f805b0ac27357e886a1b392636a5850bfb7 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Mon, 5 Oct 2026 10:42:33 -0400 Subject: [PATCH] Test RNA phasing with complete shared Sid templates --- CHANGELOG.md | 7 ++ tests/README.md | 74 +++++++++++++++ tests/test_osteosarc_rna_phasing.py | 142 ++++++++++++++++++++++++++++ varcode/version.py | 2 +- 4 files changed, 224 insertions(+), 1 deletion(-) create mode 100644 tests/test_osteosarc_rna_phasing.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 1d923ca..f7f6505 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Change Log +## [v10.11.5](https://github.com/openvax/varcode/tree/v10.11.5) (2026-10-05) + +- Test RNA read phasing and molecular phase resolution on real Sid RNA + templates from the verified `openvax-v2` bundle (#528). Cover short-read + and ONT cis evidence, trans evidence, insufficient shared coverage and + fragment thresholds; verify that mates and supplementary records are kept. + ## [v10.11.4](https://github.com/openvax/varcode/tree/v10.11.4) (2026-10-02) - Ship test helpers, bounded fixtures, provenance/checksums and the docs, diff --git a/tests/README.md b/tests/README.md index 7752d43..618f0fd 100644 --- a/tests/README.md +++ b/tests/README.md @@ -220,3 +220,77 @@ python -m pytest -q tests/test_osteosarc_shared_targets.py tests/test_osteosarc_ The first run downloads and verifies the bundle (28 MB) and exports its members into the osteosarc cache (`OSTEOSARC_CACHE`, else the shared OpenVax cache); later runs reuse them offline. + +## Real Sid RNA phasing + +`tests/test_osteosarc_rna_phasing.py` completes the real-read phasing coverage +in [#528](https://github.com/openvax/varcode/issues/528). It uses the published +`openvax-v2` bundle, with manifest SHA256 +`b42dbce529cb5724153f035fca3783a41e673b31610a30851b2ef96d4ba43734`. +The test exports three named members into a pytest temporary directory, +preserving their original alignment records and indexed BAM format: + +| Source label | Target/member suffix | Records | Templates | Supplementary records | +|---|---|---:|---:|---:| +| `2024.06.11.bostongene.align.tcga.protocol.dr32Aligned.sorted` | `.NTF3-chr12-5494381-compound` | 12 | 6 | 0 | +| `IPISRC044_T1_sclrs_ONT.tagged` | `.NTF3-chr12-5494381-compound` | 15 | 15 | 0 | +| `IPISRC044_T1_sclrs_ONT.tagged` | `.EXOC4-chr7-133274996` | 40 | 38 | 2 | + +Each full member name is its source label followed by the suffix. The bundle +recipe and manifest retain the original source URL, assembly and record +identities. The tests compare exported record identities and multiplicities +with both the member manifest and **every record of its selected RG/QNAME +templates in the bundle's source BAM**. The short-read member includes both +mates of all six templates. The EXOC4 member keeps both supplementary records; +the default evidence filter excludes them, while explicitly enabling them +recognizes their reference alleles without adding alternate support. This +checks completeness within the published bundle, without reacquiring or +making claims about records absent from the original bundle. + +All coordinates below are one-based GRCh38, using Ensembl 81. NTF3's +`AG>GT` compound target is split into its two nonoverlapping SNVs. Two +additional companion alleles were observed in these reads: NTF3 +`12:5494466 G>A` (from the source MD/base observations) and EXOC4 +`7:133895694 G>A` (compared with Ensembl transcript sequences). Their reference +bases were checked against Ensembl 81 transcripts; germline or somatic status +is not assumed. + +The frozen fragment counts were independently audited from the original +aligned bases: apply the SAM flag and quality filters, locate each SNV with +aligned reference/query positions, and merge calls by RG/QNAME. A fragment +contributes only when both loci have unambiguous ref/alt observations; missing +coverage and conflicting or other alleles do not count as reference. These +counts are expectations for the evidence API, rather than another phasing +implementation used as an oracle. + +| Member/case | First allele | Second allele | Both alt | First alt only | Second alt only | Both ref | Expected phase | +|---|---|---|---:|---:|---:|---:|---| +| Short-read NTF3 | `12:5494381 A>G` | `12:5494382 G>T` | 3 | 0 | 0 | 3 | Cis | +| ONT NTF3 | `12:5494381 A>G` | `12:5494382 G>T` | 9 | 0 | 0 | 5 | Cis | +| ONT NTF3 companion | `12:5494381 A>G` | `12:5494466 G>A` | 0 | 9 | 3 | 2 | Trans | +| ONT EXOC4 companion | `7:133274996 G>T` | `7:133895694 G>A` | 1 | 0 | 1 | 0 | Unknown | + +The defaults require MAPQ/base quality 20, a five-base read-edge exclusion, +proper pairing for paired reads, and at least two supporting fragments; +duplicates, secondary and supplementary alignments are filtered as evidence. +The ONT cases allow a 5% allele error rate; the short-read case uses the default +1%. EXOC4 has 19 and 3 alternate-supporting templates at the individual loci, +but only one carries both alternate alleles, so expression alone does not +establish phase. Raising the NTF3 fragment minimum to 10 also changes its cis +call to unknown, with nine supporting fragments. Tests exercise both +`RNAReadPhasingSource` and `MolecularPhaseResolver` in either query order. + +These are reproducible RNA observations under the configured evidence model, +not independently validated tumor haplotypes or estimates of population allele +fractions: the bundle uses bounded, allele-balanced selection. The existing +synthetic tests separately cover the germline/subclone rule fixed in #527. + +To reproduce, install the reference data as described in `CONTRIBUTING.md`, then: + +```sh +python -m pip install -e '.[test-data,rna]' +python -m pytest -q tests/test_osteosarc_rna_phasing.py +``` + +The first run can download the verified 28 MB shared bundle; later runs reuse +it offline. No sibling checkout or BAM copies in `tests/data` are required. diff --git a/tests/test_osteosarc_rna_phasing.py b/tests/test_osteosarc_rna_phasing.py new file mode 100644 index 0000000..3846825 --- /dev/null +++ b/tests/test_osteosarc_rna_phasing.py @@ -0,0 +1,142 @@ +"""RNA phasing from complete templates in the published openvax-v2 bundle. + +These are observed allele combinations, not independently validated tumor +haplotypes. See tests/README.md for the sources, manual fragment audit and +filtering policy. No alignments or bases are synthesized for these tests. +""" + +from collections import Counter +import hashlib +import json +from pathlib import Path + +import pytest +from pyensembl import cached_release + +from varcode import MolecularPhaseResolver, RNAReadPhasingSource, Variant + + +SHORT_SOURCE = "2024.06.11.bostongene.align.tcga.protocol.dr32Aligned.sorted" +LONG_SOURCE = "IPISRC044_T1_sclrs_ONT.tagged" +SHORT_NTF3 = SHORT_SOURCE + ".NTF3-chr12-5494381-compound" +LONG_NTF3 = LONG_SOURCE + ".NTF3-chr12-5494381-compound" +LONG_EXOC4 = LONG_SOURCE + ".EXOC4-chr7-133274996" +MEMBERS = (SHORT_NTF3, LONG_NTF3, LONG_EXOC4) +MANIFEST_SHA256 = "b42dbce529cb5724153f035fca3783a41e673b31610a30851b2ef96d4ba43734" + +NTF3 = ("12", 5494381, "A", "G") +NTF3_COMPOUND_PARTNER = ("12", 5494382, "G", "T") +# Read-observed companion alleles; no germline/somatic status is assumed. +NTF3_COMPANION = ("12", 5494466, "G", "A") +EXOC4 = ("7", 133274996, "G", "T") +EXOC4_COMPANION = ("7", 133895694, "G", "A") + + +@pytest.fixture(scope="module") +def shared_phasing_reads(tmp_path_factory): + osteosarc = pytest.importorskip("osteosarc", reason="install .[test-data]") + pytest.importorskip("pysam", reason="install .[rna]") + root = Path(osteosarc.fetch_bundle("openvax-v2")) + manifest_bytes = (root / "manifest.json").read_bytes() + assert hashlib.sha256(manifest_bytes).hexdigest() == MANIFEST_SHA256 + recipe = json.loads((root / "recipe.json").read_text()) + assert recipe["selection"]["supplementary"] is True + compound = recipe["targets"]["NTF3-chr12-5494381-compound"] + assert (compound["assembly"], compound["contig"], compound["position"], + compound["ref"], compound["alt"]) == ("GRCh38", "chr12", 5494381, "AG", "GT") + paths = osteosarc.export_bundle( + "openvax-v2", tmp_path_factory.mktemp("sid-rna-phasing"), members=MEMBERS) + return paths, root, json.loads(manifest_bytes) + + +@pytest.mark.parametrize("member,records,templates,supplementary", [ + (SHORT_NTF3, 12, 6, 0), + (LONG_NTF3, 15, 15, 0), + (LONG_EXOC4, 40, 38, 2), +]) +def test_shared_phasing_inputs_keep_complete_templates( + shared_phasing_reads, member, records, templates, supplementary): + from osteosarc.records import read_records + + paths, root, manifest = shared_phasing_reads + exported = list(read_records(paths[member])) + identities = Counter(record.digest for record in exported) + assert identities == Counter(manifest["members"][member]["records"]) + selected = {record.template for record in exported} + source = manifest["members"][member]["source"] + # Compare to every record of each selected RG/QNAME in the bundle's source, + # rather than refetching a locus and losing mates or distant split pieces. + complete = Counter(record.digest for record in read_records( + root / manifest["sources"][source]["bam"]) if record.template in selected) + assert identities == complete + assert len(exported) == records + assert len(selected) == templates + assert sum(record.read.is_supplementary for record in exported) == supplementary + if member == SHORT_NTF3: + for template in selected: + assert {record.read.flag & 0xc0 for record in exported + if record.template == template} == {0x40, 0x80} + + +@pytest.mark.parametrize("member,alleles,counts,support,phase", [ + (SHORT_NTF3, (NTF3, NTF3_COMPOUND_PARTNER), (3, 0, 0, 3), (3, 3), True), + (LONG_NTF3, (NTF3, NTF3_COMPOUND_PARTNER), (9, 0, 0, 5), (9, 9), True), + (LONG_NTF3, (NTF3, NTF3_COMPANION), (0, 9, 3, 2), (9, 3), False), + (LONG_EXOC4, (EXOC4, EXOC4_COMPANION), (1, 0, 1, 0), (19, 3), None), +]) +def test_real_sid_rna_fragment_evidence( + shared_phasing_reads, member, alleles, counts, support, phase): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) for allele in alleles) + source = RNAReadPhasingSource( + str(paths[member]), variants=variants, + # Q20/MAPQ20, five-base edge exclusion and SAM flag filters remain + # at their defaults. Allow a 5% allele error rate for the ONT data. + phasing_error_rate=0.01 if member == SHORT_NTF3 else 0.05) + try: + assert source._phase_counts(*variants) == counts + assert tuple(source.supports_variant(v) for v in variants) == support + assert all(source.has_evidence(v) for v in variants) + assert source.in_cis(*variants) is phase + assert source.in_cis(*reversed(variants)) is phase + resolver = MolecularPhaseResolver(source) + assert resolver.in_cis(*variants) is phase + for variant, partner in (variants, tuple(reversed(variants))): + assert tuple(resolver.phased_partners(variant)) == ((partner,) if phase is True else ()) + finally: + source.close() + + +def test_real_sid_cis_call_requires_enough_fragments(shared_phasing_reads): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) + for allele in (NTF3, NTF3_COMPOUND_PARTNER)) + source = RNAReadPhasingSource( + str(paths[LONG_NTF3]), variants=variants, min_alt_reads=10, phasing_error_rate=0.05) + try: + assert tuple(source.supports_variant(v) for v in variants) == (9, 9) + assert not any(source.has_evidence(v) for v in variants) + assert source.in_cis(*variants) is None + assert MolecularPhaseResolver(source).in_cis(*variants) is None + finally: + source.close() + + +@pytest.mark.parametrize("skip_supplementary", [True, False]) +def test_real_sid_supplementary_records_obey_the_evidence_filter( + shared_phasing_reads, skip_supplementary): + import pysam + + paths, _, _ = shared_phasing_reads + variant = Variant(*EXOC4, genome=cached_release(81)) + with pysam.AlignmentFile(paths[LONG_EXOC4], "rb") as bam: + supplementary = [read for read in bam if read.is_supplementary] + assert len(supplementary) == 2 + source = RNAReadPhasingSource( + str(paths[LONG_EXOC4]), skip_supplementary=skip_supplementary, phasing_error_rate=0.05) + try: + assert [source._read_allele(read, variant) for read in supplementary] == ( + [None, None] if skip_supplementary else ["ref", "ref"]) + assert source.supports_variant(variant) == 19 + finally: + source.close() diff --git a/varcode/version.py b/varcode/version.py index d2f4cce..6c37e06 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.11.4" +__version__ = "10.11.5"