From f3121da73a3675dd84a82e56624b4f9b2e961d7a Mon Sep 17 00:00:00 2001 From: ecrum19 Date: Thu, 8 Oct 2026 18:19:35 +0200 Subject: [PATCH] Link symbolic-allele and breakend records to the genes their REF span overlaps The interval join skipped every record with a symbolic ALT or a breakend, though such a record writes its REF, anchor base included, like any other. A skipped record has no gene link, so a vcfp:LinkedSelector never selects it, and a prohibition on a gene panel does not reach it. vcf-rdfizer-testing experiment 17's arm 2 found this with the record-level check that found the `*` gap in v3.3.0. The 104 high-coverage 1000 Genomes genomes carry structural-variant calls (, , ), and the rule withholding the 28 cancer-predisposition genes from research released the ones inside those genes: 110 records to a disease-specific research view and 58 to a general-research view, where the bcftools baseline withheld them. Carrier lists agreed throughout; only the record counts showed it. Every record whose REF is DNA bases now keys its REF span, whatever its ALT. INFO/END and breakend mates are still not read, so a structural variant links to the genes its anchor lies in, not to every gene it affects; the docs say so and point to a vcfp:RegionSelector where that matters. The SPDI join is unchanged: symbolic alleles and breakends have no allele sequence to express. Co-Authored-By: Claude Opus 5.5 --- docs/datalinking.md | 12 ++++++---- docs/limitations.md | 18 ++++++++------ docs/policy-demonstrator.md | 3 ++- test/test_linking_unit.py | 24 ++++++++++++++----- .../linkers/ensembl-genes-grch38/README.md | 7 +++--- vcf_rdfizer_data/linkers/gene-demo/README.md | 4 ++-- vcf_rdfizer_linking/runner.py | 14 +++++------ 7 files changed, 51 insertions(+), 31 deletions(-) diff --git a/docs/datalinking.md b/docs/datalinking.md index 5fea7b7..77de9bd 100644 --- a/docs/datalinking.md +++ b/docs/datalinking.md @@ -132,12 +132,14 @@ Chromosome names must match exactly (`1` and `chr1` are different) unless the manifest declares `vcfl:contigAliases`: a digest-pinned sequence map (§4, allele identity) through which both the GFF3 seqids and the records' contigs resolve to an accession. `ensembl-genes-grch38` uses one, so Ensembl's `17` meets `chr17`. -There is no liftover. Symbolic alleles and -breakends are skipped and counted in `skipped_records`; a spanning-deletion -`*` keeps the record's REF span, so it is linked (it was skipped up to v3.3.0). -The runner does not interpret INFO/END or confidence ranges. +There is no liftover. A record keys its REF span whatever its ALT: a +spanning-deletion `*` (skipped up to v3.3.0), a symbolic allele or a breakend +(skipped up to v3.3.1) is linked by the REF it writes. Only a record whose REF +is not DNA bases is skipped and counted in `skipped_records`. +The runner does not interpret INFO/END, breakend mates or confidence ranges. It uses explicit DNA REF spans, including their anchor base, rather than an -inferred biological affected region. +inferred biological affected region, so a structural variant links to the genes +its anchor lies in, not to every gene it affects. References may be HTTPS URLs, local `file:` URLs, or relative IRI references such as `` resolved beside `linker.ttl`. An absolute URL string is diff --git a/docs/limitations.md b/docs/limitations.md index f3b35d9..b3cc196 100644 --- a/docs/limitations.md +++ b/docs/limitations.md @@ -255,13 +255,17 @@ the remaining proposal. The base graph remains a faithful, separate artifact. **A rule over links reaches only linked records.** `vcfp:LinkedSelector` selects records through their calls' links, so a record a linker could not key -is never selected: a prohibition on a gene panel does not reach a symbolic -allele or breakend inside the panel, because the gene linkers key explicit REF -spans only. The link report counts these records per linker -(`skipped_records`). Check it before relying on such a prohibition, or add a -`vcfp:RegionSelector` over the same coordinates. Up to v3.3.0 the gene linkers -also skipped `*` alleles, which released 13 records in cancer-predisposition -genes to a research view of one whole genome; v3.3.1 links them. +is never selected. The gene linkers key every record by its REF span, so a +structural variant is selected by the genes its anchor lies in, not by every +gene its END or breakend mate reaches: a prohibition on a gene panel does not +reach a deletion anchored upstream of the panel. Add a `vcfp:RegionSelector` +over the same coordinates where that matters. The link report counts the +records a linker could not key per linker (`skipped_records`). Up to v3.3.0 the +gene linkers skipped `*` alleles, which released 13 records in +cancer-predisposition genes to a research view of one whole genome; up to +v3.3.1 they skipped symbolic alleles and breakends, which released structural +variants in those genes from 104 cohort genomes to research views: 110 to a +disease-specific one and 58 to a general one. **No disclosure control.** Conversion is all-or-nothing: every sample, every genotype, every header line and every free-text `Description` goes into the diff --git a/docs/policy-demonstrator.md b/docs/policy-demonstrator.md index dee8d56..8306395 100644 --- a/docs/policy-demonstrator.md +++ b/docs/policy-demonstrator.md @@ -122,7 +122,8 @@ It fails closed: a graph with no `?predicate` triple at all, because the link graph was left out, is refused rather than evaluated as selecting nothing. It does not fail closed per record: a record the linker could not key (the link report's `skipped_records`) is never selected, so a prohibition on a -panel does not reach it. +panel does not reach it. The gene linkers key every record by its REF span, so +a structural variant is selected by the genes its anchor lies in. Selector types are read from the profile files, and also from the policy file itself, so a policy can bring its own (§8). diff --git a/test/test_linking_unit.py b/test/test_linking_unit.py index 3169b6c..ebde6f1 100644 --- a/test/test_linking_unit.py +++ b/test/test_linking_unit.py @@ -148,12 +148,13 @@ def test_interval_closed_boundaries_ref_span_and_exact_chromosome(self): self.assertEqual(list(index.overlaps(LinkKey(chrom="chr1", start=100, end=100))), []) row = next(r for r in read_vcf(EXAMPLE) if isinstance(r, Record)) self.assertEqual(keys_for(replace(row, pos="99", ref="AT"), self.interval), [LinkKey(chrom="1", start=99, end=100)]) - # `*` and `.` keep the REF span; END, symbolic alleles and breakends do not. - for alt in ("*", "G,*", "."): + # Every ALT keeps the REF span: `*`, `.`, symbolic alleles and breakends + # alike. END and breakend mates are not read. + for alt in ("*", "G,*", ".", "", "", "N]2:20]", "G,", "*,"): self.assertEqual(keys_for(replace(row, pos="99", ref="AT", alt=alt), self.interval), [LinkKey(chrom="1", start=99, end=100)]) - for alt in ("", "N]2:20]", "G,", "*,"): - self.assertEqual(keys_for(replace(row, alt=alt), self.interval), []) + # Only a REF that is not bases keys nothing. + self.assertEqual(keys_for(replace(row, ref="."), self.interval), []) def test_nested_gff_features_are_not_missed(self): gff = self.root / "nested.gff3" @@ -696,10 +697,21 @@ def test_a_spanning_deletion_allele_links_to_the_gene_it_lies_in(self): # NB72462M has 13 `*` records in cancer genes; unlinked, a rule withholding # those genes from research released them. linked = self.link(self.manifest(), allele_record("chr17", 150, "AT", "*", row=1), - allele_record("chr17", 199, "C", "T,*", row=2), - allele_record("chr17", 150, "A", "", row=3)) + allele_record("chr17", 199, "C", "T,*", row=2)) self.assertEqual(linked, {"1": {"https://example.org/GENE_A"}, "2": {"https://example.org/GENE_A"}}) + def test_a_symbolic_allele_or_breakend_links_by_its_anchor(self): + # Experiment 17's arm 2 released 110 structural variants in cancer genes + # (, , ) to a research view: unkeyed, no rule over + # gene links reached them. The anchor links; END is not read, so a + # anchored outside the gene does not. + linked = self.link(self.manifest(), allele_record("chr17", 150, "A", "", row=1), + allele_record("chr17", 150, "A", "", row=2), + allele_record("chr17", 150, "A", "A]chr2:20]", row=3), + allele_record("chr17", 250, "A", "", row=4)) + self.assertEqual(linked, {"1": {"https://example.org/GENE_A"}, "2": {"https://example.org/GENE_A"}, + "3": {"https://example.org/GENE_A"}}) + def test_an_unmapped_contig_still_matches_by_name(self): linked = self.link(self.manifest(), allele_record("chrUn_x", 10, "A", "G")) self.assertEqual(linked, {"1": {"https://example.org/GENE_U"}}) diff --git a/vcf_rdfizer_data/linkers/ensembl-genes-grch38/README.md b/vcf_rdfizer_data/linkers/ensembl-genes-grch38/README.md index 440f3fb..0995201 100644 --- a/vcf_rdfizer_data/linkers/ensembl-genes-grch38/README.md +++ b/vcf_rdfizer_data/linkers/ensembl-genes-grch38/README.md @@ -10,9 +10,10 @@ Contig names are resolved through the GRCh38 sequence map `17` and a VCF's `chr17` meet. Only `gene` features count (21,581 in release 116, including all protein-coding genes); `ncRNA_gene` and pseudogenes do not. -A record's span is its REF span whatever its ALT, `*` included. Symbolic -alleles and breakends, whose extent is END or a mate, are not linked; the -report counts them in `skipped_records`. +A record's span is its REF span whatever its ALT: `*`, symbolic alleles and +breakends included. END and breakend mates are not read, so a structural +variant links to the genes its anchor lies in. Only a record whose REF is not +DNA bases is skipped; the report counts it in `skipped_records`. ```bash vcf-rdfizer-link run -i sample.vcf --link ensembl-genes-grch38 -o sample.links.nt diff --git a/vcf_rdfizer_data/linkers/gene-demo/README.md b/vcf_rdfizer_data/linkers/gene-demo/README.md index 3823914..f27c8ad 100644 --- a/vcf_rdfizer_data/linkers/gene-demo/README.md +++ b/vcf_rdfizer_data/linkers/gene-demo/README.md @@ -13,8 +13,8 @@ vcf-rdfizer-link dry-run my-gene-linker -i sample.vcf --limit 100 Intervals are 1-based closed; keys cover POS through POS + len(REF) - 1. The example contains overlapping genes A/B and a gene C on chromosome 2. Exons are -ignored. Chromosome names match exactly. A `*` ALT keeps the REF span; -symbolic alleles and breakends are skipped. +ignored. Chromosome names match exactly. Every ALT keeps the REF span: `*`, +symbolic alleles and breakends link by their anchor. For a real reference, change the plug-in ID, reference URL, actual SHA-256, assembly and object template. Point `idAttribute` at the GFF3 attribute your diff --git a/vcf_rdfizer_linking/runner.py b/vcf_rdfizer_linking/runner.py index 88528d5..32eb053 100644 --- a/vcf_rdfizer_linking/runner.py +++ b/vcf_rdfizer_linking/runner.py @@ -108,13 +108,13 @@ def keys_for(record, manifest): if allele: keys.append(LinkKey(chrom=record.chrom, start=allele[0], token=f"{allele[1]}:{allele[2]}")) return keys - # POS and GFF3 intervals are both 1-based closed. Explicit REF spans only; - # END, symbolic alleles and breakends require a different coordinate policy. - # `*` (an allele spanning an upstream deletion) and `.` leave the REF span - # as written, so they key like any other ALT: skipping them would let a - # rule that selects records by their genes miss records inside those genes. - if not re.fullmatch(r"[ACGTNacgtn]+", record.ref) or any( - not re.fullmatch(r"[ACGTNacgtn]+|[.*]", alt) for alt in record.alt.split(",")): + # POS and GFF3 intervals are both 1-based closed. A record keys its REF span + # whatever its ALT: `*`, `.`, symbolic alleles and breakends all leave the + # REF, anchor base included, as written, and skipping them would let a rule + # that selects records by their genes miss records inside those genes. + # INFO/END and breakend mates are not read, so a structural variant links by + # its anchor, not by the whole extent it affects. + if not re.fullmatch(r"[ACGTNacgtn]+", record.ref): return [] try: start = int(record.pos)