diff --git a/claude/maps/commands.md b/claude/maps/commands.md index 51a8fd86c..80329dbb2 100644 --- a/claude/maps/commands.md +++ b/claude/maps/commands.md @@ -144,7 +144,7 @@ Generated by `vg map commands` (do not edit; run `scripts/vg map commands` after | `fix_migrate_flags_to_imported_allele_info` | classification | (helper module, not a command) (guessed) | classification, flags, snpdb | | `fix_nbk_citations` | classification | | classification | | `fix_orphanet_revalidation` | classification | | classification | -| `fix_rematch_release_symbols_to_genes` | genes | | genes | +| `fix_rematch_release_symbols_to_genes` | genes | Resyncs ReleaseGeneSymbolGene (the symbol -> gene matches for a GeneAnnotationRelease) with what | genes | | `fix_tag_colors_collection_permissions` | snpdb | | snpdb | | `fix_truncated_annotation_runs` | annotation | Find AnnotationRuns whose VEP output was silently truncated, and mark them ERROR so they can be | annotation | | `fix_variant_annotation_add_hgvs_g` | annotation | | annotation | diff --git a/claude/plans/1669_release_gene_matcher_alias_chaining_plan.md b/claude/plans/1669_release_gene_matcher_alias_chaining_plan.md deleted file mode 100644 index 14893d338..000000000 --- a/claude/plans/1669_release_gene_matcher_alias_chaining_plan.md +++ /dev/null @@ -1,269 +0,0 @@ -# #1669 — ReleaseGeneMatcher: single-hop alias matching - -Written by Claude Fable 5 (claude-fable-5), 2026-08-31 - -[#1669](https://github.com/SACGF/variantgrid/issues/1669): `ReleaseGeneMatcher` walks the gene symbol -alias table transitively, so an alias string shared by two unrelated genes bridges them. Reported from -prod (SACGF/variantgrid_sapath#426) as a hypogonadotropic hypogonadism panel returning a **PDCD2** -variant via `MT-TS2 → RP8 → PDCD2`. - -The fix is one idea: **a symbol matches a release gene directly, or through exactly one -`GeneSymbolAlias` row that touches a symbol the release knows.** Everything else in this plan is -cleaning up the rows the old traversal wrote. - ---- - -## 1. Diagnosis (confirmed against the code and the local DB) - -`genes/gene_matching.py`, `ReleaseGeneMatcher`: - -* `aliases_dict` (line 165) builds `alias_graph` with every `GeneSymbolAlias` row hung under **both** - its `alias` and its `gene_symbol_id` (lines 173–176), then starts a walk from both ends of every row. -* `_aliases()` (line 186) recurses through that graph until it reaches any symbol in `self.genes` - (symbols the release actually has), then assigns those gene ids to **every symbol on the path** - (lines 196–211). Path length is unbounded; only loops are guarded. - -Local DB reproduces the prod rows exactly — `MT-TS2 → 5134 (PDCD2)` in releases 3, 4, 7, 8 and -`MT-TS2 → ENSG05220024333` in release 5. - -Two further findings that shape the cleanup: - -1. **`match_info` cannot identify chained rows.** `symbol_match_path` is a dict keyed by symbol, so - when the walk re-enters a symbol (`RP8 → MT-TS2 → RP8 → PDCD2`, which the loop guard allows because - the root symbol is never added to `visited_symbols`) the earlier hop is overwritten and the stored - text collapses to a single hop. Release 3's `MT-TS2` row reads `RP8 is an alias for PDCD2 (HGNC)`, - so the issue's `LIKE '%), %'` audit query undercounts. -2. **Rows also go stale when the alias table changes.** `genes/cached_web_resource/refseq.py` deletes - and re-imports NCBI aliases, but `ReleaseGeneSymbolGene` rows derived from the deleted aliases - survive (`match_symbols_to_genes` only inserts). Release 1 has `ADAM12-OT1 → ADAM12` citing - `ADAM12-OT1 is an alias for ADAM12 (NCBI)`; the only alias rows for `ADAM12-OT1` today are - `CAR10` and `FLJ31066` (HGNC). - -So the cleanup recomputes what the fixed matcher would produce and diffs against the table, rather -than pattern-matching `match_info`. - -### What the prototype single-hop matcher says about existing rows - -Alias/gene-version derived rows (`match_info IS NOT NULL`) per release on the local DB, classified by -the shortest explanation the single-hop rules give them: - -| release | rows | gene-version | forward hop only | backward hop only | multi-hop / stale | -|---|---|---|---|---|---| -| 1 GRCh38 RefSeq 110 | 1110 | 0 | 129 | 111 | 870 | -| 2 GRCh37 RefSeq 105 | 2012 | 82 | 177 | 113 | 1640 | -| 3 GRCh37 RefSeq 105 | 1502 | 91 | 180 | 113 | 1118 | -| 4 GRCh38 RefSeq 2023 | 1288 | 0 | 299 | 110 | 879 | -| 5 T2T Ensembl 2022-06 | 2249 | 0 | 171 | **1110** | 968 | -| 7 GRCh38 RefSeq 2024-08 | 1275 | 0 | 315 | 112 | 848 | -| 8 GRCh38 RefSeq 2025-08 | 606 | 0 | 128 | 125 | 353 | - -The "backward hop only" column is why the hop stays bidirectional (§2). The last column is what the -resync command (§4) removes. - ---- - -## 2. Design: one hop, either direction - -`GeneSymbolAlias(alias=A, gene_symbol=S)` reads "A is an alias for S". For a queried symbol `Q` that -is **absent** from the release: - -* **Forward** — `Q == A`, `S` in release → `Q` matches `S`'s genes. Gene list says `KAL1`, release - has `ANOS1`. This is the common "list uses an old name" case. -* **Backward** — `Q == S`, `A` in release → `Q` matches `A`'s genes. Gene list says `ANOS1`, release - still calls the gene `KAL1`. HGNC is refreshed far more often than gene annotation releases, so a - current-symbol gene list against an older release hits this constantly — 1,110 rows in the T2T - Ensembl 2022 release alone (`ABTB3`, `ACE2-DT`, `ACTMAP`, `ADISSP`, …), and ~110 per RefSeq release. - -One hop is enough for renames: HGNC's `prev_symbol` lists *every* previous symbol of a gene, so a -gene renamed twice has direct alias rows from both old names. Chaining is what lets an alias string -shared by two genes (`RP8`, `ALP`, `CD`, `ADMR`) act as a bridge, and dropping it stops -`MT-TS2 → RP8 → PDCD2`, `ATHS → ALP → {SLPI, ATRNL1, CCL27, PDLIM3, ASRGL1}`, -`CELIAC2 → CD → {NOD2, CTLA4}`, `ACKR5 → ADMR → GPR182`. - -Precedence, unchanged from today's `_get_gene_id_and_match_info_for_symbol`: - -1. Symbol is in the release → its genes, `match_info=None`, aliases never consulted. -2. Otherwise the union of gene-version matches (`Gene v0/GRCh38`, from `_get_genes_dict`) and - single-hop alias matches. Where both name the same gene the gene-version text wins (it is set - first; alias entries use `setdefault`). - -A symbol can still legitimately map to several genes (e.g. `RP8` forward-matches both `MT-TS2` and -`PDCD2`; only `PDCD2` is in the release, so `RP8` → `PDCD2` is correct and kept). - -The residual risk of the backward hop — HGNC's informal `alias_symbol` occasionally being another -gene's approved symbol (`AURKAIP1` lists `AIP`) — only bites when the queried symbol is absent from -the release *and* the alias is present as a different gene. Distinguishing `prev_symbol` from -`alias_symbol` is [#1668](https://github.com/SACGF/variantgrid/issues/1668) item 2; once that -provenance exists, the backward hop should accept `prev_symbol` rows only. `ReleaseGeneMatcher` -becomes a consumer of #1668's shared resolver at that point. - ---- - -## 3. `genes/gene_matching.py` - -Replace `aliases_dict` and delete `_aliases`: - -```python -@cached_property -def aliases_dict(self) -> dict[str, dict]: - """ Upper-case symbol -> {gene_id: match_info} for symbols the release lacks, reached either - from an older GeneVersion's symbol or by a single GeneSymbolAlias hop (in either direction) - to a symbol the release has. One hop is deliberate: HGNC lists every previous symbol of a - gene, so renames never need chaining, and chaining lets an alias string shared by two - unrelated genes bridge them (#1669). """ - genes_dict = self._get_genes_dict() - alias_qs = GeneSymbolAlias.objects.all() - for gsa in alias_qs: - alias = gsa.alias.upper() - symbol = gsa.gene_symbol_id.upper() - if alias == symbol: - continue - for query, target in ((alias, symbol), (symbol, alias)): - if query in self.genes: - continue - for gene_id in self.genes.get(target, []): - genes_dict[query].setdefault(gene_id, gsa.match_info) - return genes_dict -``` - -Notes for the implementer: - -* `self.genes` keys are `Upper("gene_symbol")`; `_get_genes_dict` keys are `gene_symbol_id.upper()`. - Keep everything upper-cased on the Python side — `alias` has `case_insensitive` collation in - Postgres but that does nothing for dict lookups. -* Loading all alias rows replaces the old `gene_symbol__releasegenesymbol__release=self.release` - filter. With one hop both ends are checked against `self.genes` directly, so the pre-filter buys - nothing, and removing it means `aliases_dict` (a `cached_property`) no longer depends on which - `ReleaseGeneSymbol` rows existed at the moment it was first computed. -* `_get_genes_dict`, `_get_gene_id_and_match_info_for_symbol`, `match_symbols_to_genes` and - everything below are unchanged. `ReleaseGeneSymbolGene.match_info` keeps its existing format - (`"X is an alias for Y (HGNC)"` / `"Gene v0/GRCh38"`), so `gene_grid.js` keeps rendering it. - ---- - -## 4. Resync command: `fix_rematch_release_symbols_to_genes` - -Extend the existing command (`genes/management/commands/fix_rematch_release_symbols_to_genes.py`) -from add-only to a full resync of the derived table, per release: - -1. `release_gene_symbols = list(gar.releasegenesymbol_set.all())` -2. `expected = gm._get_gene_id_and_match_info_for_symbol(rgs.gene_symbol_id for rgs in ...)` — - a dict `gene_symbol_id -> [(gene_id, match_info)]`. -3. Load existing rows: `ReleaseGeneSymbolGene.objects.filter(release_gene_symbol__release=gar) - .values_list("pk", "release_gene_symbol__gene_symbol_id", "gene_id", "match_info")`. -4. Diff on `(gene_symbol_id, gene_id)`: - * existing pair absent from `expected` → **delete** (`filter(pk__in=...)` in batches of 2000) - * pair present in both with different `match_info` → **update** (`bulk_update`) - * pair only in `expected` → **insert** (reuse `gm.match_symbols_to_genes(release_gene_symbols)` - after the deletes; it already does `ignore_conflicts=True`) -5. Print per-release counts: deleted / updated / inserted, and symbols left with no gene (as now). - At `--verbosity 2` also print each deleted `(symbol, gene_id, match_info)` so an operator can - eyeball what went. - -Add `--dry-run` (report the diff, make no writes) so it can be run on prod for inspection first. - -Cache: `GeneAnnotationRelease.genes_for_symbols` is a plain queryset with no Redis caching, and the -analysis gene-list nodes query `ReleaseGeneSymbolGene` at run time, so the change takes effect on -the next node load. Existing `NodeCount` figures for already-loaded gene-list nodes reflect the old -rows until the node re-runs; that is the normal behaviour after any gene-data change and needs no -extra invalidation. - ---- - -## 5. Migration: `genes/migrations/0089_one_off_resync_release_gene_symbol_genes.py` - -Follow `genes/migrations/0065_one_off_fix_make_panel_app_gene_lists_public.py`: - -```python -def _has_alias_derived_release_gene_matches(apps): - ReleaseGeneSymbolGene = apps.get_model("genes", "ReleaseGeneSymbolGene") - return ReleaseGeneSymbolGene.objects.filter(match_info__contains="is an alias for").exists() - -operations = [ - ManualOperation(task_id=ManualOperation.task_id_manage(["fix_rematch_release_symbols_to_genes"]), - note="Remove chained-alias gene matches (e.g. MT-TS2 -> RP8 -> PDCD2) and resync " - "ReleaseGeneSymbolGene to single-hop matching (#1669)", - test=_has_alias_derived_release_gene_matches), -] -``` - -Dependency: `("genes", "0088_one_off_stamp_existing_pfam_domains_imported")`. - ---- - -## 6. Tests — `genes/tests/test_gene_matching.py`, new `TestReleaseGeneMatcher` - -Build one small release in `setUpTestData`. `annotation/tests/test_data_fake_genes.py` has -`_create_fake_gene_version(genome_build, gene_id, symbol, consortium)`; the release needs a -`GeneAnnotationImport`, a `GeneAnnotationRelease(version, annotation_consortium, genome_build, -gene_annotation_import)` and a `ReleaseGeneVersion(release, gene_version)` per gene — a local helper -`_release_gene(symbol, gene_id)` wrapping those keeps the fixture readable. - -Release genes: `PDCD2`=`5134`, `ANOS1`=`3730`, `OLDNAME`=`111` (a gene the release still carries -under its previous symbol). Aliases (all HGNC): `RP8→MT-TS2`, `RP8→PDCD2`, `KAL1→ANOS1`, -`OLDNAME→NEWNAME`. Also create `GeneSymbol` rows for `MT-TS2`, `KAL1`, `NEWNAME`, `RP8`. - -Call `gm = ReleaseGeneMatcher(release)` and assert on -`gm._get_gene_id_and_match_info_for_symbol([...])`: - -| test | query | expected | -|---|---|---| -| `test_direct_symbol_wins` | `PDCD2` | `[("5134", None)]` | -| `test_forward_alias_hop` | `KAL1` | `[("3730", "KAL1 is an alias for ANOS1 (HGNC)")]` | -| `test_backward_alias_hop` | `NEWNAME` | `[("111", "OLDNAME is an alias for NEWNAME (HGNC)")]` | -| `test_shared_alias_does_not_bridge` | `MT-TS2` | nothing (`RP8` is absent from the release, so the walk stops) | -| `test_alias_matches_only_release_target` | `RP8` | `[("5134", "RP8 is an alias for PDCD2 (HGNC)")]` | - -Plus one command test: create `ReleaseGeneSymbol`s for `PDCD2`, `KAL1`, `MT-TS2`, plant a -`ReleaseGeneSymbolGene(MT-TS2 → 5134, "RP8 is an alias for PDCD2 (HGNC)")` and a direct -`PDCD2 → 5134`, run `call_command("fix_rematch_release_symbols_to_genes")`, and assert the planted -chained row is gone, the direct row remains, and `KAL1 → 3730` was inserted. Run the same with -`--dry-run` and assert nothing changed. - -Run with `python3 manage.py test --keepdb genes.tests.test_gene_matching`, then the wider -`genes.tests` and `annotation.tests.test_gene_level_annotation` (it builds `ReleaseGeneSymbolGene` -rows by hand and exercises `genes_for_symbol`). - ---- - -## 7. Changelog - -`variantgrid/templates/default_templates/changelog.html`, current release block: - -```html -
  • #1669 - Gene list matching follows a single alias hop; symbols sharing an alias string (e.g. MT-TS2 / PDCD2 via "RP8") are no longer linked
    • -``` - ---- - -## 8. Surfacing the substitution to the user - -GeneGrid already shows `Matched : ` for alias-matched symbols -(`variantgrid/static_files/default_static/js/gene_grid.js:384`). Showing the same on the gene list -page is #1668 item 4 and stays with that issue. - ---- - -## 9. Immediate prod workaround (before deploy) - -```sql -DELETE FROM genes_releasegenesymbolgene rgsg -USING genes_releasegenesymbol rgs -WHERE rgs.id = rgsg.release_gene_symbol_id - AND rgs.gene_symbol_id = 'MT-TS2' - AND rgsg.match_info LIKE '%PDCD2%'; -``` - -The full cleanup is §4 run via the §5 migration task. - ---- - -## Files - -| file | change | -|---|---| -| `genes/gene_matching.py` | `aliases_dict` rewritten to single hop; `_aliases` removed | -| `genes/management/commands/fix_rematch_release_symbols_to_genes.py` | add-only → full resync, `--dry-run` | -| `genes/migrations/0089_one_off_resync_release_gene_symbol_genes.py` | `ManualOperation` for the resync | -| `genes/tests/test_gene_matching.py` | `TestReleaseGeneMatcher` | -| `variantgrid/templates/default_templates/changelog.html` | entry | diff --git a/claude/plans/tso500_overall_plan.md b/claude/plans/tso500_overall_plan.md index 4097570f8..40f652a6b 100644 --- a/claude/plans/tso500_overall_plan.md +++ b/claude/plans/tso500_overall_plan.md @@ -404,8 +404,9 @@ one, so they want to stay stable from the first client. Phase 6's gene-symbol item rests on that assumption holding. Phase 5's fusion parser goes through the same resolver for `SEPT14` → `SEPTIN14`, so one check against a real database covers both — and if the alias is missing, a fusion partner still imports, just under a local `GENE:` id rather than its HGNC one. -[`1669_release_gene_matcher_alias_chaining_plan.md`](1669_release_gene_matcher_alias_chaining_plan.md) -changes alias resolution to single-hop, so do the check after that lands. +[#1669](https://github.com/SACGF/variantgrid/issues/1669) made `ReleaseGeneMatcher` single-hop; that is a +different resolver from this one, but do the check against a database that has had +`fix_rematch_release_symbols_to_genes` run. **Clients send a build's own name (`GRCh37`), not an alias.** `GenomeBuild.get_name_or_alias("hg19")` raises `MultipleObjectsReturned` rather than `DoesNotExist`, so a declared build that will not resolve diff --git a/claude/research/genes.md b/claude/research/genes.md index 6aef86115..b6a1227c1 100644 --- a/claude/research/genes.md +++ b/claude/research/genes.md @@ -81,8 +81,11 @@ A symbol does not name a gene; it names a gene *in a release*. The same symbol h RefSeq ids over time (and TAZ became TAFAZZIN between builds), so `genes/gene_matching.py:ReleaseGeneMatcher` writes `genes/models/models_gene_annotation_release.py:ReleaseGeneSymbol` / `ReleaseGeneSymbolGene` rows per release: first a direct hit on the release's own GeneVersion symbols, then `genes/gene_matching.py:ReleaseGeneMatcher.aliases_dict`, -which walks the alias graph (`ReleaseGeneMatcher._aliases`, loop-guarded) and the symbols other builds' GeneVersions -gave the same gene, recording the path as `match_info` so the gene list grid can show why. Readers never touch the +which takes a single `GeneSymbolAlias` hop in either direction (the symbol is the alias, or the alias is of the symbol) +plus the symbols other builds' GeneVersions gave the same gene, recording the hop as `match_info` so the gene list grid +can show why. The hop count is the whole point: chaining hops let an alias string shared by two unrelated genes bridge +them, which is how an MT-TS2 gene list matched PDCD2 via "RP8" (#1669), and HGNC lists every previous symbol of a gene +so renames never need more than one. Readers never touch the matcher: `genes/models/models_gene_annotation_release.py:GeneAnnotationRelease.genes_for_symbols` and `genes/models/models_gene_list.py:GeneList.get_genes` read the cached rows, and the release always comes from the VAV (`GeneAnnotationRelease.get_for_latest_annotation_versions_for_builds`). Matching is triggered whenever symbols @@ -247,8 +250,10 @@ dicts on first use and `GeneSymbolMatcher.create_gene_list_gene_symbols` re-matc a loop creating lists should share one matcher. Inserting `GeneListGeneSymbol` rows any other way leaves symbols with no release rows and the list matches nothing in analyses until `genes/management/commands/rematch_unmatched_gene_list_symbols.py:Command` runs; -`genes/management/commands/fix_rematch_release_symbols_to_genes.py:Command` re-runs the alias walk for symbols that -have a release row but no gene. +`genes/management/commands/fix_rematch_release_symbols_to_genes.py:Command` resyncs the whole derived table for every +release - it inserts, updates `match_info` and deletes matches the current rules no longer make, which matching itself +never does (`match_symbols_to_genes` only inserts), so rows left behind by an alias re-import or an older matcher need +it. Run it `--dry-run` first. Genes prefixed `unknown_` are legacy placeholders from pre-GFF imports; `genes/management/commands/fix_fake_genes.py:Command` re-points their transcripts where another version names the gene and `Gene.delete_orphaned_fake_genes` removes the diff --git a/genes/CLAUDE.md b/genes/CLAUDE.md index 956c78aea..5908a6e94 100644 --- a/genes/CLAUDE.md +++ b/genes/CLAUDE.md @@ -32,6 +32,8 @@ Gotchas: - `PanelAppPanel.cache_valid` expires after `settings.PANEL_APP_CACHE_DAYS` (models/models_panel_app.py:PanelAppPanel.cache_valid); panel_app.py:get_panel_app_local_cache re-fetches from the live API when stale, so tests must not depend on it. - GeneCoverageCollection is a partitioned model (models/models_gene_coverage.py:GeneCoverageCollection); delete via the model so partitions are dropped. - gene_matching.py:GeneSymbolMatcher and gene_matching.py:ReleaseGeneMatcher cache whole-table dicts on first use; build one per import, not per symbol. +- gene_matching.py:ReleaseGeneMatcher takes exactly one GeneSymbolAlias hop (either direction) - chaining hops lets an alias string shared by two unrelated genes bridge them (#1669). Keep it single-hop. +- Matching only ever inserts ReleaseGeneSymbolGene rows, so a rematch can't remove a match that's since become wrong; `manage.py fix_rematch_release_symbols_to_genes` (`--dry-run` first) is the full resync that also updates and deletes. - `` and `` have no HGVS at all - neither a ranged form nor an explicit expansion - so hgvs/hgvs_matcher.py:HGVSMatcher raises hgvs/hgvs_converter.py:HGVSNoRepresentationException before any converter runs, and classification records it as `ResolvedVariantInfo.error` rather than a Rollbar bug. Tests: - annotation/tests/test_data_fake_genes.py:create_fake_transcript_version builds Gene/GeneVersion/Transcript/TranscriptVersion (RUNX1, ENST00000300305.7) for a build; `create_gata2_transcript_version` / `create_pten_transcript_version` add RefSeq examples. diff --git a/genes/gene_matching.py b/genes/gene_matching.py index 1653f738d..128845cb5 100644 --- a/genes/gene_matching.py +++ b/genes/gene_matching.py @@ -170,70 +170,33 @@ def _get_genes_dict(self): @cached_property def aliases_dict(self) -> dict[str, dict]: - """ Get symbols from other GeneVersions that match genes from our release """ - genes_dict = self._get_genes_dict() - - # Gene Symbol alias - qs = GeneSymbolAlias.objects.filter(gene_symbol__releasegenesymbol__release=self.release) - gene_symboli_alias_list = [gsa for gsa in qs if gsa.alias != gsa.gene_symbol_id] - - alias_graph = defaultdict(list) - for gsa in gene_symboli_alias_list: - alias_graph[gsa.alias].append(gsa) - alias_graph[gsa.gene_symbol_id].append(gsa) - - for gene_symbol_alias in gene_symboli_alias_list: - for gene_symbol in [gene_symbol_alias.alias, gene_symbol_alias.gene_symbol_id]: - symbol_match_path = {gene_symbol: gene_symbol_alias.match_info} + """ Upper case symbol -> {gene_id: match_info} for symbols the release doesn't have, reached either + from an older GeneVersion's symbol or by a single GeneSymbolAlias hop (in either direction) to a + symbol the release does have. - self._aliases(alias_graph, genes_dict, gene_symbol, symbol_match_path) - - return genes_dict + One hop is deliberate: HGNC lists every previous symbol of a gene, so renames never need + chaining, while chaining lets an alias string shared by two unrelated genes bridge them - e.g. + 'RP8' is an alias of both MT-TS2 and PDCD2, which used to match MT-TS2 gene lists to PDCD2 (#1669) """ + genes_dict = self._get_genes_dict() - def _aliases(self, alias_graph, genes_dict, gene_symbol, symbol_match_path, visited_symbols=None): - # print(f"_aliases(alias_graph, genes_dict, {gene_symbol} - ({symbol_match_path})") - # Keep track of visited symbols to detect loops in graph - if visited_symbols is None: - visited_symbols = set() - else: - if gene_symbol in visited_symbols: # Stop unnecessary descent - return - visited_symbols.add(gene_symbol) - - if gene_id_list := self.genes.get(gene_symbol): - # Only need to build up match path for where we are - symbol_total_path = {} - match_paths = list(symbol_match_path.values()) - for i, symbol in enumerate(symbol_match_path): - mp = ", ".join(match_paths[i+1:]) - symbol_total_path[symbol] = mp - - for symbol in symbol_match_path: - if symbol != gene_symbol: # No point putting actual one in - match_info = symbol_total_path[symbol] + # Ordered so the match_info a symbol/gene pair reached by several alias rows gets is stable across runs + alias_values = GeneSymbolAlias.objects.order_by("pk").values_list("alias", "gene_symbol_id", "source") + for alias, symbol, source in alias_values.iterator(): + uc_alias = alias.upper() + uc_symbol = symbol.upper() + if uc_alias == uc_symbol: + continue + + # "alias is an alias for symbol" - a gene list can use either end and mean the other + for query, target_symbol in ((uc_alias, uc_symbol), (uc_symbol, uc_alias)): + if query in self.genes: + continue # Release has the symbol itself, so aliases are never consulted for it + if gene_id_list := self.genes.get(target_symbol): + match_info = GeneSymbolAlias(alias=alias, gene_symbol_id=symbol, source=source).match_info for gene_id in gene_id_list: - if existing_match_info := genes_dict[symbol].get(gene_id): - if len(match_info) >= len(existing_match_info): - continue # Don't override with longer - genes_dict[symbol][gene_id] = match_info - else: - if aliases_list := alias_graph.get(gene_symbol): - original_gene_symbol = gene_symbol - - for gene_symbol_alias in aliases_list: - # We may have looked it up via alias or gene symbol - use the other one + genes_dict[query].setdefault(gene_id, match_info) - if gene_symbol_alias.gene_symbol_id == original_gene_symbol: - gene_symbol = gene_symbol_alias.alias - else: - gene_symbol = gene_symbol_alias.gene_symbol_id - - # Make a copy for recursion - child_symbol_match_path = symbol_match_path.copy() - child_symbol_match_path[gene_symbol] = gene_symbol_alias.match_info - - self._aliases(alias_graph, genes_dict, gene_symbol, child_symbol_match_path, - visited_symbols=visited_symbols) + return genes_dict def _get_gene_id_and_match_info_for_symbol(self, gene_symbols) -> dict[str, list]: gene_symbol_gene_id_and_match_info = defaultdict(list) # list items = (gene_id, match_info) diff --git a/genes/management/commands/fix_rematch_release_symbols_to_genes.py b/genes/management/commands/fix_rematch_release_symbols_to_genes.py index 1bb17d309..2f910a2f4 100644 --- a/genes/management/commands/fix_rematch_release_symbols_to_genes.py +++ b/genes/management/commands/fix_rematch_release_symbols_to_genes.py @@ -1,24 +1,65 @@ +""" +Resyncs ReleaseGeneSymbolGene (the symbol -> gene matches for a GeneAnnotationRelease) with what +ReleaseGeneMatcher produces today: inserts missing matches, updates changed match_info and deletes +matches the current rules no longer make. + +The delete is what makes this more than a rematch: matching only ever inserted, so rows written by +the old multi-hop alias traversal (#1669) and rows derived from alias rows since deleted by a +re-import survive a rematch. Run with --dry-run first to see the diff. +""" from django.core.management import BaseCommand from genes.gene_matching import ReleaseGeneMatcher from genes.models import GeneAnnotationRelease, ReleaseGeneSymbol, ReleaseGeneSymbolGene +BATCH_SIZE = 2000 + class Command(BaseCommand): category = "one-off" + def add_arguments(self, parser): + parser.add_argument('--dry-run', action='store_true', help="Report the diff without writing anything") + def handle(self, *args, **options): - for gar in GeneAnnotationRelease.objects.all(): - no_match_qs = ReleaseGeneSymbol.objects.filter(release=gar, releasegenesymbolgene__isnull=True) + dry_run = options["dry_run"] + verbosity = options["verbosity"] - qs = ReleaseGeneSymbolGene.objects.filter(release_gene_symbol__release=gar) - num_genes_original = qs.count() - num_no_match_original = no_match_qs.count() - print(f"{gar} - symbols w/o gene: {num_no_match_original}") + for gar in GeneAnnotationRelease.objects.all(): + release_gene_symbols = list(gar.releasegenesymbol_set.all()) gm = ReleaseGeneMatcher(gar) - release_gene_symbols = gar.releasegenesymbol_set.all() - gm.match_symbols_to_genes(release_gene_symbols) - num_genes = qs.count() + expected = gm._get_gene_id_and_match_info_for_symbol(rgs.gene_symbol_id for rgs in release_gene_symbols) + expected_matches = {} # (gene_symbol_id, gene_id) -> match_info + for gene_symbol_id, gene_id_and_match_info in expected.items(): + for gene_id, match_info in gene_id_and_match_info: + expected_matches[(gene_symbol_id, gene_id)] = match_info - print(f"{gar} - matched {num_genes - num_genes_original} genes") - print(f"{gar} - {no_match_qs.count() - num_no_match_original} less symbols w/o genes") + existing_qs = ReleaseGeneSymbolGene.objects.filter(release_gene_symbol__release=gar) + existing = existing_qs.values_list("pk", "release_gene_symbol__gene_symbol_id", "gene_id", "match_info") + + delete_pks = [] + update_records = [] + for pk, gene_symbol_id, gene_id, match_info in existing: + key = (gene_symbol_id, gene_id) + if key in expected_matches: + expected_match_info = expected_matches.pop(key) # What's left over needs inserting + if expected_match_info != match_info: + update_records.append(ReleaseGeneSymbolGene(pk=pk, match_info=expected_match_info)) + else: + delete_pks.append(pk) + if verbosity >= 2: + print(f"{gar} - delete {gene_symbol_id} -> {gene_id} ({match_info})") + + num_insert = len(expected_matches) + print(f"{gar} - delete: {len(delete_pks)}, update: {len(update_records)}, insert: {num_insert}") + + if not dry_run: + for i in range(0, len(delete_pks), BATCH_SIZE): + ReleaseGeneSymbolGene.objects.filter(pk__in=delete_pks[i:i + BATCH_SIZE]).delete() + if update_records: + ReleaseGeneSymbolGene.objects.bulk_update(update_records, ["match_info"], batch_size=BATCH_SIZE) + if num_insert: + gm.match_symbols_to_genes(release_gene_symbols) + + no_match_qs = ReleaseGeneSymbol.objects.filter(release=gar, releasegenesymbolgene__isnull=True) + print(f"{gar} - symbols w/o gene: {no_match_qs.count()}") diff --git a/genes/migrations/0089_one_off_resync_release_gene_symbol_genes.py b/genes/migrations/0089_one_off_resync_release_gene_symbol_genes.py new file mode 100644 index 000000000..f5c12eddb --- /dev/null +++ b/genes/migrations/0089_one_off_resync_release_gene_symbol_genes.py @@ -0,0 +1,22 @@ +from django.db import migrations + +from manual.operations.manual_operations import ManualOperation + + +def _has_alias_derived_release_gene_matches(apps): + ReleaseGeneSymbolGene = apps.get_model("genes", "ReleaseGeneSymbolGene") + return ReleaseGeneSymbolGene.objects.filter(match_info__contains="is an alias for").exists() + + +class Migration(migrations.Migration): + + dependencies = [ + ('genes', '0088_one_off_stamp_existing_pfam_domains_imported'), + ] + + operations = [ + ManualOperation(task_id=ManualOperation.task_id_manage(["fix_rematch_release_symbols_to_genes"]), + note="Remove chained-alias gene matches (e.g. MT-TS2 -> RP8 -> PDCD2) and resync " + "ReleaseGeneSymbolGene to single-hop matching (#1669)", + test=_has_alias_derived_release_gene_matches), + ] diff --git a/genes/tests/test_gene_matching.py b/genes/tests/test_gene_matching.py index 9325a750f..448d3a745 100644 --- a/genes/tests/test_gene_matching.py +++ b/genes/tests/test_gene_matching.py @@ -1,14 +1,30 @@ from django.contrib.auth.models import User +from django.core.management import call_command from django.test.testcases import TestCase from genes.gene_matching import ( MAX_GENE_SYMBOL_LENGTH, GeneSymbolMatcher, HGNCMatcher, + ReleaseGeneMatcher, tokenize_gene_symbols, ) -from genes.models import HGNC, GeneList, GeneSymbol, GeneSymbolAlias, HGNCImport -from genes.models_enums import GeneSymbolAliasSource, HGNCStatus +from genes.models import ( + HGNC, + Gene, + GeneAnnotationImport, + GeneAnnotationRelease, + GeneList, + GeneSymbol, + GeneSymbolAlias, + GeneVersion, + HGNCImport, + ReleaseGeneSymbol, + ReleaseGeneSymbolGene, + ReleaseGeneVersion, +) +from genes.models_enums import AnnotationConsortium, GeneSymbolAliasSource, HGNCStatus +from snpdb.models import GenomeBuild class TestGeneMatching(TestCase): @@ -130,3 +146,96 @@ def test_create_gene_list_filters_oversized_from_prefiltered_list(self): self.assertIn("oversized", gene_list.error_message.lower()) saved_names = set(gene_list.genelistgenesymbol_set.values_list("original_name", flat=True)) self.assertEqual({"BRCA1"}, saved_names) + + +class ReleaseGeneMatcherTestCase(TestCase): + """ Base fixture: a release with PDCD2/ANOS1/OLDNAME, and HGNC aliases where "RP8" is shared by + MT-TS2 (absent from the release) and PDCD2. """ + + @classmethod + def setUpTestData(cls): + genome_build = GenomeBuild.get_name_or_alias("GRCh38") + cls.gene_annotation_import = GeneAnnotationImport.objects.create( + url="fake", genome_build=genome_build, annotation_consortium=AnnotationConsortium.REFSEQ) + cls.release = GeneAnnotationRelease.objects.create( + version="test_1669", annotation_consortium=AnnotationConsortium.REFSEQ, + genome_build=genome_build, gene_annotation_import=cls.gene_annotation_import) + + def _release_gene(symbol, gene_id): + gene_symbol = GeneSymbol.objects.get_or_create(symbol=symbol)[0] + gene = Gene.objects.create(identifier=gene_id, annotation_consortium=AnnotationConsortium.REFSEQ) + gene_version = GeneVersion.objects.create(gene=gene, gene_symbol=gene_symbol, version=1, + genome_build=genome_build, + import_source=cls.gene_annotation_import) + ReleaseGeneVersion.objects.create(release=cls.release, gene_version=gene_version) + return gene + + _release_gene("PDCD2", "5134") + _release_gene("ANOS1", "3730") + _release_gene("OLDNAME", "111") # Release still carries this gene under its previous symbol + + def _alias(alias, symbol): + GeneSymbol.objects.get_or_create(symbol=alias) + GeneSymbol.objects.get_or_create(symbol=symbol) + GeneSymbolAlias.objects.create(alias=alias, gene_symbol_id=symbol, + source=GeneSymbolAliasSource.HGNC) + + # HGNC has RP8 as an alias for both MT-TS2 (not in release) and PDCD2 (in release) + _alias("RP8", "MT-TS2") + _alias("RP8", "PDCD2") + _alias("KAL1", "ANOS1") + _alias("OLDNAME", "NEWNAME") + + +class TestReleaseGeneMatcher(ReleaseGeneMatcherTestCase): + """ A symbol reaches a release's genes directly, or through exactly one GeneSymbolAlias hop (either + direction). Chaining hops let an alias string shared by two unrelated genes bridge them (#1669). """ + + def _match(self, symbol): + gm = ReleaseGeneMatcher(self.release) + return gm._get_gene_id_and_match_info_for_symbol([symbol])[symbol] + + def test_direct_symbol_wins(self): + self.assertEqual([("5134", None)], self._match("PDCD2")) + + def test_forward_alias_hop(self): + self.assertEqual([("3730", "KAL1 is an alias for ANOS1 (HGNC)")], self._match("KAL1")) + + def test_backward_alias_hop(self): + self.assertEqual([("111", "OLDNAME is an alias for NEWNAME (HGNC)")], self._match("NEWNAME")) + + def test_shared_alias_does_not_bridge(self): + self.assertEqual([], self._match("MT-TS2")) + + def test_alias_matches_only_release_target(self): + self.assertEqual([("5134", "RP8 is an alias for PDCD2 (HGNC)")], self._match("RP8")) + + +class TestFixRematchReleaseSymbolsToGenes(ReleaseGeneMatcherTestCase): + """ The resync command has to delete matches the current rules no longer make - a rematch only inserts. """ + + def setUp(self): + release_gene_symbols = {} + for symbol in ["PDCD2", "KAL1", "MT-TS2"]: + release_gene_symbols[symbol] = ReleaseGeneSymbol.objects.create(release=self.release, + gene_symbol_id=symbol) + # What the old multi-hop traversal wrote: MT-TS2 -> RP8 -> PDCD2 + ReleaseGeneSymbolGene.objects.create(release_gene_symbol=release_gene_symbols["MT-TS2"], gene_id="5134", + match_info="RP8 is an alias for PDCD2 (HGNC)") + ReleaseGeneSymbolGene.objects.create(release_gene_symbol=release_gene_symbols["PDCD2"], gene_id="5134") + + def _matched_genes(self, symbol): + return set(ReleaseGeneSymbolGene.objects.filter(release_gene_symbol__release=self.release, + release_gene_symbol__gene_symbol_id=symbol) + .values_list("gene_id", flat=True)) + + def test_resync_deletes_chained_inserts_missing_and_keeps_direct(self): + call_command("fix_rematch_release_symbols_to_genes") + self.assertEqual(set(), self._matched_genes("MT-TS2")) + self.assertEqual({"5134"}, self._matched_genes("PDCD2")) + self.assertEqual({"3730"}, self._matched_genes("KAL1")) + + def test_dry_run_writes_nothing(self): + call_command("fix_rematch_release_symbols_to_genes", dry_run=True) + self.assertEqual({"5134"}, self._matched_genes("MT-TS2")) + self.assertEqual(set(), self._matched_genes("KAL1")) diff --git a/variantgrid/templates/default_templates/changelog.html b/variantgrid/templates/default_templates/changelog.html index ed7179c70..f83d11ab1 100644 --- a/variantgrid/templates/default_templates/changelog.html +++ b/variantgrid/templates/default_templates/changelog.html @@ -150,6 +150,7 @@

    Genes

  • PanelApp panels are matched on HGNC ID rather than the gene symbol from a dated snapshot, and deleted panels are handled
  • Gene Lists - faster grid, copy/paste genes in, and a graphs tab showing gene locations
  • A mitochondrial 'm.' HGVS on a non-MT reference now warns rather than silently resolving as 'g.'
    • +
  • Gene Lists - a symbol only matches genes through a single alias hop, so genes sharing an alias string (eg MT-TS2 / PDCD2 via "RP8") are no longer linked

    • Genome Builds

    @@ -294,6 +295,10 @@

    Liftover

    • #1273 - Liftover page - retry failed liftovers per tool (eg all failed bcftools +liftover)
    +

    Genes

    +
      +
    • #1669 - Gene list matching follows a single alias hop; symbols sharing an alias string (eg MT-TS2 / PDCD2 via "RP8") are no longer linked
      • +