Link symbolic-allele and breakend records to the genes their REF span overlaps

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 `<genes.gff3>` resolved beside `linker.ttl`. An absolute URL string is

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

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).

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,*", ".", "<DEL>", "<INS:ME:ALU>", "N]2:20]", "G,<DEL>", "*,<DEL>"):
             self.assertEqual(keys_for(replace(row, pos="99", ref="AT", alt=alt), self.interval),
                              [LinkKey(chrom="1", start=99, end=100)])
-        for alt in ("<DEL>", "N]2:20]", "G,<DEL>", "*,<DEL>"):
-            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", "<DEL>", 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
+        # (<INS:ME:ALU>, <DUP>, <INS>) to a research view: unkeyed, no rule over
+        # gene links reached them. The anchor links; END is not read, so a <DUP>
+        # anchored outside the gene does not.
+        linked = self.link(self.manifest(), allele_record("chr17", 150, "A", "<DEL>", row=1),
+                           allele_record("chr17", 150, "A", "<INS:ME:ALU>", row=2),
+                           allele_record("chr17", 150, "A", "A]chr2:20]", row=3),
+                           allele_record("chr17", 250, "A", "<DUP>", 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"}})

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

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

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)