From f992b46e209df1c0320647938491f2602f5d5e78 Mon Sep 17 00:00:00 2001 From: Clare72 Date: Thu, 24 Sep 2026 14:52:49 +0100 Subject: [PATCH] prefer stocks matching full combination for splits --- ...t_expression_pattern_individual_queries.py | 70 ++++++-- src/test/test_flybase_stocks.py | 104 ++++++++++- src/vfbquery/flybase_stocks.py | 169 +++++++++++++++++- src/vfbquery/vfb_queries.py | 44 ++++- 4 files changed, 363 insertions(+), 24 deletions(-) diff --git a/src/test/test_expression_pattern_individual_queries.py b/src/test/test_expression_pattern_individual_queries.py index def192e..eadfb36 100644 --- a/src/test/test_expression_pattern_individual_queries.py +++ b/src/test/test_expression_pattern_individual_queries.py @@ -106,32 +106,74 @@ def test_instance_stock_query_matches_its_driver_live(self): """Data-health (scheduled only): the same against live SOLR.""" self.assertEqual(self._anchors_or_skip(self.EP_INDIVIDUAL), ["FBtp0060056"]) - def test_split_class_offers_a_stock_query_per_hemidriver(self): - """Render code (PR-blocking): FindStocks is generated per hemidriver from a - complete split fixture (its ``has_hemidriver`` relationships). Deterministic - — isolates the generation code from live-data state. See TESTING.md - 'Fixture vs live (data_health)'.""" - self.assertEqual( - _stock_anchors(_fixture_term_info("split_VFBexp_FBtp0129935FBtp0129968")), - ["FBtp0129935", "FBtp0129968"]) + SPLIT_COMBO = "FBco0001890" # FlyBase combination of the two + + def _split_fixture_anchors(self, combos): + """Parse the split fixture with the FlyBase combination lookup stubbed, + keeping the render test hermetic (no Chado).""" + import vfbquery.flybase_stocks as fbs + real = fbs.combinations_for_constructs + fbs.combinations_for_constructs = combos + try: + return _stock_anchors( + _fixture_term_info("split_VFBexp_FBtp0129935FBtp0129968")) + finally: + fbs.combinations_for_constructs = real + + def test_split_class_offers_one_stock_query_on_its_combination(self): + """Render code (PR-blocking): a split's two ``has_hemidriver`` constructs + collapse onto the FlyBase combination they build, so the menu offers a + single FindStocks (exact combination first, hemidriver fallback) rather + than one per half. See TESTING.md 'Fixture vs live (data_health)'.""" + seen = [] + + def combos(construct_ids): + seen.append(sorted(construct_ids)) + return [self.SPLIT_COMBO] + + self.assertEqual(self._split_fixture_anchors(combos), [self.SPLIT_COMBO]) + self.assertEqual(seen, [["FBtp0129935", "FBtp0129968"]]) + + def test_split_combination_stock_label_names_the_combination(self): + """Render code (PR-blocking): the label keeps its usual wording and names + the FBco in brackets, where it used to name each FBtp.""" + import vfbquery.flybase_stocks as fbs + real = fbs.combinations_for_constructs + fbs.combinations_for_constructs = lambda ids: [self.SPLIT_COMBO] + try: + ti = _fixture_term_info("split_VFBexp_FBtp0129935FBtp0129968") + finally: + fbs.combinations_for_constructs = real + labels = [q["label"] for q in ti["Queries"] if q.get("query") == "FindStocks"] + self.assertEqual(labels, [f"Find fly stocks for {ti['Name']} ({self.SPLIT_COMBO})"]) + + def test_split_without_curated_combination_keeps_a_query_per_hemidriver(self): + """Render code (PR-blocking): no FlyBase combination -> per-hemidriver queries.""" + self.assertEqual(self._split_fixture_anchors(lambda ids: []), + ["FBtp0129935", "FBtp0129968"]) + + def test_split_combination_lookup_failure_keeps_a_query_per_hemidriver(self): + """Render code (PR-blocking): a Chado failure must not drop the stock queries.""" + def combos(ids): + raise RuntimeError("chado down") + self.assertEqual(self._split_fixture_anchors(combos), + ["FBtp0129935", "FBtp0129968"]) @pytest.mark.data_health - def test_split_class_offers_a_stock_query_per_hemidriver_live(self): + def test_split_class_offers_one_stock_query_on_its_combination_live(self): """Data-health (scheduled only): the same against the live SOLR document, so an incomplete production split document (no ``Expression_pattern`` type / no ``has_hemidriver``) is caught. Deselected on PRs via ``-m 'not data_health'``.""" - self.assertEqual(self._anchors_or_skip(self.SPLIT_CLASS), - ["FBtp0129935", "FBtp0129968"]) + self.assertEqual(self._anchors_or_skip(self.SPLIT_CLASS), [self.SPLIT_COMBO]) @pytest.mark.data_health - def test_split_instance_inherits_hemidriver_stock_queries(self): + def test_split_instance_inherits_combination_stock_query(self): # No own driver edge: the features come from the pattern class it # instantiates, which `_stock_features_via_parent_pattern` fetches with a # LIVE `_load_term_info(parent)` call. A fixture for the individual alone # cannot make this hermetic (the parent lookup still hits SOLR), so this # stays a live (data_health) test rather than a PR-blocking fixture one. - self.assertEqual(self._anchors_or_skip(self.SPLIT_INDIVIDUAL), - ["FBtp0129935", "FBtp0129968"]) + self.assertEqual(self._anchors_or_skip(self.SPLIT_INDIVIDUAL), [self.SPLIT_COMBO]) class TestExpressionPatternIndividualQueries(unittest.TestCase): diff --git a/src/test/test_flybase_stocks.py b/src/test/test_flybase_stocks.py index 2501edf..7bae36b 100644 --- a/src/test/test_flybase_stocks.py +++ b/src/test/test_flybase_stocks.py @@ -1,7 +1,14 @@ """Tests for flybase_stocks module — entity resolution and stock discovery.""" +import pandas as pd import pytest -from vfbquery.flybase_stocks import resolve_entity, find_stocks +from vfbquery.flybase_stocks import ( + COMBO_MATCH_EXACT, + COMBO_MATCH_HEMIDRIVER, + COMBO_MATCH_OTHER_COMBO, + resolve_entity, + find_stocks, +) from vfbquery.vfb_queries import ( _flybase_report_url, _md_link, @@ -231,6 +238,101 @@ def test_nonexistent_combination(self): stocks = find_stocks("FBco9999999") assert stocks == [] + @pytest.mark.integration + def test_exact_combination_hides_hemidriver_stocks(self): + # FBco0001000 = R50C12-p65.AD ∩ VT049369-GAL4.DBD; 607102 carries both, + # so the single-hemidriver stocks (89605, 72840, 603713) are dropped. + stocks = find_stocks("FBco0001000") + assert [s["stock_number"] for s in stocks] == ["607102"] + assert stocks[0]["match"] == COMBO_MATCH_EXACT + + @pytest.mark.integration + def test_fallback_is_per_hemidriver(self): + # FBco0002000 has no exact stock. Its R20A03-GAL4.DBD has stocks of its + # own, so its other-combination stocks are hidden; its R82F03-p65.AD has + # none, so it falls back to the stocks pairing it with other DBDs. + stocks = {s["stock_number"]: s["match"] for s in find_stocks("FBco0002000")} + assert stocks["75808"] == COMBO_MATCH_HEMIDRIVER # R20A03 DBD alone + assert stocks["602487"] == COMBO_MATCH_OTHER_COMBO # R82F03 AD + R48G01 DBD + assert "602978" not in stocks # R20A03 DBD + R52B07 AD + + @pytest.mark.integration + def test_uncurated_partner_hemidriver_is_detected(self): + # 88099 pairs R24H08-GAL4.DBD with VT058427-p65.AD, which is in no FBco, + # so it must be recognised by the split-system tool it encodes. + import vfbquery.flybase_stocks as fbs + from vfbquery.flybase_db import get_connection + conn = get_connection(statement_timeout_ms=60000) + try: + df = fbs._run_query(conn, fbs._STOCK_HEMIDRIVERS_SQL, + {"stock_ids": ["FBst0088099"]}) + finally: + conn.close() + assert "FBal0331043" in set(df["allele_id"]) + + @pytest.mark.integration + def test_match_column_only_for_combinations(self): + combo = get_flybase_stocks("FBco0001000", return_dataframe=False) + assert combo["headers"]["match"]["type"] == "text" + assert combo["rows"][0]["match"] == COMBO_MATCH_EXACT + allele = get_flybase_stocks("FBal0034227", return_dataframe=False, limit=3) + assert "match" not in allele["headers"] + + +class TestCombinationRanking: + """Ranking of combination stocks, with the chado queries stubbed out.""" + + AD, DBD, OTHER_DBD = "FBal_AD", "FBal_DBD", "FBal_DBD2" + + def _stocks(self, monkeypatch, per_allele, hemidrivers): + import vfbquery.flybase_stocks as fbs + + def fake_run_query(conn, sql, params): + if sql is fbs._COMBO_COMPONENTS_SQL: + return pd.DataFrame({"allele_name": ["AD", "DBD"], + "allele_id": [self.AD, self.DBD]}) + if sql is fbs._STOCK_HEMIDRIVERS_SQL: + return pd.DataFrame(hemidrivers, columns=["stock_id", "allele_id"]) + raise AssertionError("unexpected query") + + def fake_allele(conn, allele_id, collection_filter=None): + return pd.DataFrame( + [{"stock_id": sid, "stock_number": sid[-1], "genotype": "", + "collection": "BDSC"} for sid in per_allele.get(allele_id, [])], + columns=["stock_id", "stock_number", "genotype", "collection"]) + + monkeypatch.setattr(fbs, "_run_query", fake_run_query) + monkeypatch.setattr(fbs, "_find_stocks_allele", fake_allele) + df = fbs._find_stocks_combination(None, "FBco_test") + return list(zip(df["stock_id"], df["match"])) + + def test_exact_combination_hides_lower_tiers(self, monkeypatch): + got = self._stocks( + monkeypatch, + per_allele={self.AD: ["st1", "st2", "st3"], self.DBD: ["st3", "st4"]}, + hemidrivers=[("st1", self.AD), ("st1", self.OTHER_DBD), ("st2", self.AD), + ("st3", self.AD), ("st3", self.DBD), ("st4", self.DBD)]) + assert got == [("st3", COMBO_MATCH_EXACT)] + + def test_hemidriver_alone_hides_its_other_combinations(self, monkeypatch): + # AD has a stock of its own (st2), so st1 is hidden; DBD has only st5. + got = self._stocks( + monkeypatch, + per_allele={self.AD: ["st1", "st2"], self.DBD: ["st4", "st5"]}, + hemidrivers=[("st1", self.AD), ("st1", self.OTHER_DBD), ("st2", self.AD), + ("st4", self.DBD), ("st5", self.DBD)]) + assert got == [("st2", COMBO_MATCH_HEMIDRIVER), + ("st4", COMBO_MATCH_HEMIDRIVER), + ("st5", COMBO_MATCH_HEMIDRIVER)] + + def test_hemidriver_without_own_stock_falls_back_to_other_combinations(self, monkeypatch): + # DBD has a stock of its own; AD only appears paired with another DBD. + got = self._stocks( + monkeypatch, + per_allele={self.AD: ["st1"], self.DBD: ["st2"]}, + hemidrivers=[("st1", self.AD), ("st1", self.OTHER_DBD), ("st2", self.DBD)]) + assert got == [("st2", COMBO_MATCH_HEMIDRIVER), ("st1", COMBO_MATCH_OTHER_COMBO)] + class TestFindStocksTableSchema: CONSTRUCT_WITH_STOCKS = "FBtp0000352" # P{GawB} diff --git a/src/vfbquery/flybase_stocks.py b/src/vfbquery/flybase_stocks.py index 4e50cad..bb4fa63 100644 --- a/src/vfbquery/flybase_stocks.py +++ b/src/vfbquery/flybase_stocks.py @@ -392,6 +392,74 @@ def _resolve_entity_impl(conn, name_or_id): ORDER BY a.uniquename """ +# Split system component alleles (hemidrivers) carried by each stock. A genotype +# usually lists the FBti insertion rather than the allele, so reach the allele +# the same three ways _ALLELE_STOCKS_SQL does: held directly, via the construct +# the insertion was produced by, or via an associated_with insertion. An allele +# counts as a hemidriver when some combination is partially_produced_by it, or +# when it encodes a split-system tool (GAL4(DBD)::Zip-, p65(AD)::Zip+, ...). +# The tool test matters: many hemidrivers (e.g. Hsap\RELA[AD.VT058427]) sit in +# stocks without ever being curated into a combination. The tool set is read +# from the combinations themselves rather than hard-coded. +_STOCK_HEMIDRIVERS_SQL = """ +WITH hemi_tools AS ( + SELECT DISTINCT et.object_id AS tool_fid + FROM feature_relationship ppb + JOIN cvterm c ON ppb.type_id = c.cvterm_id AND c.name = 'partially_produced_by' + JOIN feature_relationship et ON ppb.object_id = et.subject_id + JOIN cvterm ec ON et.type_id = ec.cvterm_id AND ec.name = 'encodes_tool' +), gf AS ( + SELECT DISTINCT s.uniquename AS stock_id, fg.feature_id + FROM stock s + JOIN stock_genotype sg ON s.stock_id = sg.stock_id + JOIN feature_genotype fg ON sg.genotype_id = fg.genotype_id + WHERE s.uniquename = ANY(%(stock_ids)s) +), carried AS ( + SELECT gf.stock_id, gf.feature_id AS allele_fid + FROM gf + + UNION + + SELECT gf.stock_id, fr2.subject_id + FROM gf + JOIN feature_relationship fr1 ON gf.feature_id = fr1.subject_id + JOIN cvterm c1 ON fr1.type_id = c1.cvterm_id AND c1.name = 'producedby' + JOIN feature_relationship fr2 ON fr1.object_id = fr2.object_id + + UNION + + SELECT gf.stock_id, fr.subject_id + FROM gf + JOIN feature_relationship fr ON gf.feature_id = fr.object_id + JOIN cvterm c ON fr.type_id = c.cvterm_id AND c.name = 'associated_with' +) +SELECT DISTINCT carried.stock_id, a.uniquename AS allele_id +FROM carried +JOIN feature a ON carried.allele_fid = a.feature_id AND a.is_obsolete = false +WHERE EXISTS ( + SELECT 1 + FROM feature_relationship ppb + JOIN cvterm c ON ppb.type_id = c.cvterm_id AND c.name = 'partially_produced_by' + WHERE ppb.object_id = a.feature_id +) OR EXISTS ( + SELECT 1 + FROM feature_relationship et + JOIN cvterm ec ON et.type_id = ec.cvterm_id AND ec.name = 'encodes_tool' + JOIN hemi_tools ht ON et.object_id = ht.tool_fid + WHERE et.subject_id = a.feature_id +) +""" + +# Stock match tiers for a split system combination, best first. +COMBO_MATCH_EXACT = "Exact combination" +COMBO_MATCH_HEMIDRIVER = "Hemidriver alone" +COMBO_MATCH_OTHER_COMBO = "Hemidriver in other combination" +_COMBO_MATCH_RANK = { + COMBO_MATCH_EXACT: 0, + COMBO_MATCH_HEMIDRIVER: 1, + COMBO_MATCH_OTHER_COMBO: 2, +} + def _add_collection_filter(sql, params, collection_filter, use_where=False): """Add optional collection filter to a stock query.""" @@ -447,7 +515,19 @@ def _find_stocks_construct(conn, construct_id, collection_filter=None): def _find_stocks_combination(conn, combo_id, collection_filter=None): - """Find stocks for a split system combination via its component alleles.""" + """Find stocks for a split system combination via its component alleles. + + Each stock is given a ``match`` tier, and only the best available tier is + returned: + 1. Exact combination — the stock carries every component hemidriver. + If any exist, nothing else is returned. + 2. Hemidriver alone — one component, with no other hemidriver. + 3. Hemidriver in other combination — one component, paired with a + hemidriver that is not part of this combination. + Tiers 2 and 3 are chosen per hemidriver: a hemidriver with a stock of its + own shows only those, while one without falls back to the stocks that + pair it with something else. + """ components = _run_query(conn, _COMBO_COMPONENTS_SQL, {"combo_id": combo_id}) if components.empty: return pd.DataFrame() @@ -463,14 +543,95 @@ def _find_stocks_combination(conn, combo_id, collection_filter=None): if not frames: return pd.DataFrame() + hits = pd.concat(frames, ignore_index=True) + component_ids = set(components["allele_id"]) + + # A stock reached through more than one component carries all of them. + stocks = ( + hits.groupby("stock_id", sort=False) + .agg( + stock_number=("stock_number", "first"), + genotype=("genotype", "first"), + collection=("collection", "first"), + component=("component", lambda s: "; ".join(sorted(set(s)))), + component_id=("component_id", lambda s: "; ".join(sorted(set(s)))), + n_components=("component_id", "nunique"), + ) + .reset_index() + ) + + hemidrivers = _run_query( + conn, _STOCK_HEMIDRIVERS_SQL, {"stock_ids": list(stocks["stock_id"])}) + other_combo_stocks = set() + if not hemidrivers.empty: + other = hemidrivers[~hemidrivers["allele_id"].isin(component_ids)] + other_combo_stocks = set(other["stock_id"]) + + def _match(row): + if row["n_components"] >= len(component_ids): + return COMBO_MATCH_EXACT + if row["stock_id"] in other_combo_stocks: + return COMBO_MATCH_OTHER_COMBO + return COMBO_MATCH_HEMIDRIVER + + stocks["match"] = stocks.apply(_match, axis=1) + stocks["_rank"] = stocks["match"].map(_COMBO_MATCH_RANK) + + # Keep only the best tier. Non-exact stocks carry a single component, so + # grouping on it applies the fallback to each hemidriver independently. + exact = stocks["match"] == COMBO_MATCH_EXACT + if exact.any(): + stocks = stocks[exact] + else: + best = stocks.groupby("component_id")["_rank"].transform("min") + stocks = stocks[stocks["_rank"] == best] + return ( - pd.concat(frames, ignore_index=True) - .drop_duplicates(subset=["stock_id"]) - .sort_values(["collection", "stock_number"]) + stocks.sort_values(["_rank", "component", "collection", "stock_number"], + na_position="last") + .drop(columns=["_rank", "n_components"]) .reset_index(drop=True) ) +# Split system combinations whose hemidrivers are made from all of the given +# constructs. VFB models a split expression pattern by its two hemidriver +# constructs (FBtp) only, so this is how a pattern reaches its FBco. +_COMBOS_FOR_CONSTRUCTS_SQL = """ +SELECT combo.uniquename AS combo_id +FROM feature tp +JOIN feature_relationship ta ON tp.feature_id = ta.object_id +JOIN cvterm tc ON ta.type_id = tc.cvterm_id + AND tc.name IN ('associated_with', 'derived_tp_assoc_alleles') +JOIN feature_relationship ppb ON ta.subject_id = ppb.object_id +JOIN cvterm pc ON ppb.type_id = pc.cvterm_id AND pc.name = 'partially_produced_by' +JOIN feature combo ON ppb.subject_id = combo.feature_id AND combo.is_obsolete = false +WHERE tp.uniquename = ANY(%(construct_ids)s) +GROUP BY combo.uniquename +HAVING count(DISTINCT tp.uniquename) = %(n)s +ORDER BY combo.uniquename +""" + + +def combinations_for_constructs(construct_ids): + """Return the FBco ids of split system combinations built from all of + ``construct_ids`` (the hemidriver FBtp constructs of a split pattern). + + :param construct_ids: FBtp ids, e.g. ``["FBtp0099469", "FBtp0099561"]`` + :return: list of FBco uniquenames (usually one; empty if none is curated) + """ + construct_ids = sorted(set(construct_ids)) + if len(construct_ids) < 2: + return [] + conn = get_connection(statement_timeout_ms=30000) + try: + df = _run_query(conn, _COMBOS_FOR_CONSTRUCTS_SQL, + {"construct_ids": construct_ids, "n": len(construct_ids)}) + finally: + conn.close() + return [] if df.empty else list(df["combo_id"]) + + def _find_stock_details(conn, stock_id): """Look up details for a specific stock ID.""" return _run_query(conn, _STOCK_DETAILS_SQL, {"stock_id": stock_id}) diff --git a/src/vfbquery/vfb_queries.py b/src/vfbquery/vfb_queries.py index af80324..ee36982 100644 --- a/src/vfbquery/vfb_queries.py +++ b/src/vfbquery/vfb_queries.py @@ -1517,11 +1517,14 @@ def term_info_parse_object(results, short_form): # Split-GAL4 image instances carry no driver edge of their own; # follow the cached parent link to the pattern class's features. ep_feature_ids = _stock_features_via_parent_pattern(vfbTerm) + ep_feature_ids = _split_combination_ids(ep_feature_ids) multi = len(ep_feature_ids) > 1 for fb_id in ep_feature_ids: - # Disambiguate the label only when a pattern drives several - # features (a split), so single-feature patterns keep a clean name. - stock_name = f"{termInfo['Name']} ({fb_id})" if multi else termInfo["Name"] + # Name the feature only for a split (its FBco, or each hemidriver + # when no combination is curated), so single-driver patterns keep + # a clean label. + split = multi or fb_id.startswith("FBco") + stock_name = f"{termInfo['Name']} ({fb_id})" if split else termInfo["Name"] q = FindStocks_to_schema(stock_name, {"short_form": fb_id}) queries.append(q) @@ -2712,6 +2715,26 @@ def _stock_features_from_relationships(vfbTerm): return feature_ids +def _split_combination_ids(feature_ids): + """Collapse a split pattern's hemidriver constructs onto its combination. + + A split pattern reaches FindStocks as its two hemidriver FBtp constructs, + which gave one stock query per half. The FlyBase combination (FBco) built + from them ranks stocks carrying both halves first and falls back to each + hemidriver only when needed, so offer that single query instead. Anything + other than a construct pair with a curated combination — or a FlyBase + lookup failure — keeps the per-feature queries. + """ + if len(feature_ids) < 2 or not all(f.startswith("FBtp") for f in feature_ids): + return feature_ids + from .flybase_stocks import combinations_for_constructs + try: + return combinations_for_constructs(feature_ids) or feature_ids + except Exception as e: + print(f"Could not resolve split combination for {feature_ids}: {e}") + return feature_ids + + def _stock_features_via_parent_pattern(vfbTerm): """Feature IDs for an expression-pattern instance with no driver edge. @@ -5349,12 +5372,18 @@ def get_flybase_stocks(short_form: str, return_dataframe=True, limit: int = -1): # entry where that centre has one, and the collection name resolves to the # centre's homepage. `id` stays the bare FBst — it is the row's selection # id, not a rendered cell. + # Split system combination (FBco) stocks come back ranked exact + # combination -> hemidriver alone -> hemidriver in another combination; + # show that tier so the ordering reads as deliberate. Other feature types + # carry no `match`, and keep their four-column table. + has_match = any(s.get('match') for s in stocks) + rows = [] for s in stocks: stock_id = s.get('stock_id', '') or '' stock_number = s.get('stock_number', '') or '' collection = s.get('collection', '') or '' - rows.append({ + row = { # Hidden identity column (the FBst the row is about). Without a # `selection_id`-typed column the website consumes the first data # column as the row identity and hides it — which dropped Stock ID @@ -5366,7 +5395,10 @@ def get_flybase_stocks(short_form: str, return_dataframe=True, limit: int = -1): 'genotype': s.get('genotype', ''), 'collection': _md_link( collection, homepages.get(collection, {}).get('homepage_url')), - }) + } + if has_match: + row['match'] = s.get('match', '') or '' + rows.append(row) total_count = len(rows) if limit != -1: @@ -5382,6 +5414,8 @@ def get_flybase_stocks(short_form: str, return_dataframe=True, limit: int = -1): 'genotype': {'title': 'Genotype', 'type': 'text', 'order': 2}, 'collection': {'title': 'Collection', 'type': 'markdown', 'order': 3}, } + if has_match: + headers['match'] = {'title': 'Match', 'type': 'text', 'order': 4} return {'headers': headers, 'rows': rows, 'count': total_count}