diff --git a/checks.sh b/checks.sh new file mode 100755 index 0000000..696eb08 --- /dev/null +++ b/checks.sh @@ -0,0 +1,49 @@ +#!/usr/bin/env bash +# Run exactly what CI runs, and say what each step did. +# +# Two failures made this worth having. A commit went out with six mypy +# errors because the command was piped and the pipeline's exit code was the +# pipe's, not mypy's. And PR #290 failed CI lint while `mypy src` passed +# locally -- CI runs bare `mypy`, which covers 159 files rather than 152, +# and the extra seven include the tests. +# +# So: no pipes around the commands, every exit code printed, one verdict at +# the end. Mirrors .github/workflows/ci.yml; if that changes, change this. +set -uo pipefail +cd "$(dirname "$0")" + +PY=.venv/bin/python +[ -x "$PY" ] || PY=$(command -v python3) + +declare -a names=() codes=() +step() { + local name="$1"; shift + "$@" >/tmp/checks-$$.log 2>&1 + local code=$? + names+=("$name"); codes+=("$code") + if [ "$code" -eq 0 ]; then + printf ' %-16s ok\n' "$name" + else + printf ' %-16s FAILED (exit %d)\n' "$name" "$code" + sed 's/^/ /' /tmp/checks-$$.log | tail -25 + fi +} + +echo "checks:" +step "ruff check" "$PY" -m ruff check . +step "ruff format" "$PY" -m ruff format --check . +step "mypy" "$PY" -m mypy +step "pytest" "$PY" -m pytest +rm -f /tmp/checks-$$.log + +failed=0 +for i in "${!codes[@]}"; do + [ "${codes[$i]}" -ne 0 ] && failed=$((failed + 1)) +done + +if [ "$failed" -eq 0 ]; then + echo "all ${#names[@]} checks passed" + exit 0 +fi +echo "$failed of ${#names[@]} checks FAILED" +exit 1 diff --git a/specs/012-run-gene-analysis/spec.md b/specs/012-run-gene-analysis/spec.md new file mode 100644 index 0000000..8a1d2b7 --- /dev/null +++ b/specs/012-run-gene-analysis/spec.md @@ -0,0 +1,155 @@ +# Feature Specification: Run a gene set analysis from the chat + +**Feature Branch**: `012-run-gene-analysis` + +**Created**: 2026-09-21 + +**Status**: Draft + +**Input**: User description: "Run a gene set analysis from an uploaded expression matrix in the chat, and return the pathway results as a downloadable file" + +## Context + +Reactome runs gene set analysis at `gsa.reactome.org` — PADOG, Camera +("similar to the classical GSEA algorithm"), ssGSEA, terapadog. The chatbot +cannot start one. The gap was found the hard way: asked to "run a GSEA with +my list of genes", the chatbot replied that Reactome could not, and offered +`fgsea` and a YouTube tutorial. + +Five catalogue tools were added to reactome-mcp in September (methods, data +types, dataset search, examples, sources). They let the model *describe* the +service. None of them runs anything. + +**Everything below was measured against the live service on 2026-09-21**, not +read off the specification. + +### What the service actually requires + +`POST /analysis` takes `datasets[].data`: the whole expression matrix, inline, +as a tab-delimited string. It is a required field and there is no +by-reference variant. + +| | measured | +|---|---| +| matrix, 16-sample example dataset | **1.2 MB** | +| submitted payload | **1.5 MB** | +| result | **2.0 MB** — 2,679 pathways × 9 columns | +| run time | minutes (PADOG, 1,000 permutations) | + +A full round trip was completed: load → summary → submit → poll → result. +Reactome release 97, columns `Pathway, Name, Direction, FDR, PValue, NGenes, +MeanAbsT0, MeanWeightT0, av_foldchange`, plus a Pathway Browser link carrying +the analysis token. + +### Two routes in, and only one needs an upload + +`POST /data/load/{resourceId}` loads a **public** dataset by identifier — +Expression Atlas, Single Cell Expression Atlas, GREIN, GEO, and the bundled +examples. `GET /data/summary/{id}` then returns the sample IDs and factors +(measured: 16 samples; `condition` = MCM/MOCK, `cell.type` = PBMCB/TIBC, +`patient.id` = P1–P4), which is exactly what a user needs to choose a +comparison, and is small enough to show them. + +The user's own data has no such route. It must be uploaded. + +### Where the data goes, and where the decision is + +Uploading is **not** the boundary that matters. Reactome is not OICR: the +service runs on AWS operated by the Reactome project, and reactome.org's own +analysis page already accepts exactly these files. A user uploading their +matrix here is doing what they can already do on the website. Users know what +they are uploading. + +**The boundary this feature moves is the results reaching OpenAI**, and only +if the user asks for a summary. That already has a mechanism and a warning, +and this feature reuses both rather than inventing either. + +## User Scenarios & Testing *(mandatory)* + +### User Story 1 — Analyse a public dataset (Priority: P1) + +A user says "run a gene set analysis on GSE12345 comparing treated with +control". The chatbot loads the dataset by identifier, shows the sample +groups it found, runs the analysis, and returns the significant pathways plus +a downloadable table. + +**Why this priority**: it is the whole feature minus the upload, it needs no +file handling and no disk, and it is the case a user with no data of their own +can still reach. If only this ships, the chatbot can run real analyses. + +**Acceptance**: given a public dataset identifier and a named comparison, the +chat reports the pathway count, shows the top pathways with FDR, offers the +full table as a file, and links the Pathway Browser view. + +### User Story 2 — Analyse an uploaded matrix (Priority: P2) + +A user uploads their own expression matrix through the chat UI, names the +comparison, and gets the same output. + +**Why this priority**: it is what was asked for, and it is P2 only because it +depends on the submit-and-poll machinery that Story 1 already builds. + +**Acceptance**: a matrix uploaded through the UI reaches the service and +produces results, and the file is deleted from the container once submitted. + +### User Story 3 — Results without a model summary (Priority: P1) + +A user who declines the summary still gets the pathway table and the Pathway +Browser link. + +**Why this priority**: same priority as Story 1 deliberately. The analysis is +the product; the summary is a convenience. A user must never have to send +results to a third party in order to see their own results. + +**Acceptance**: with the summary declined, the file and the link are present +and nothing has been sent to the model. + +## Requirements *(mandatory)* + +- **FR-001** The matrix MUST NOT pass through the model's context or an MCP + tool call. At 1.2 MB it would bury the answer and cost more than the + analysis. The chatbot submits it server-side. +- **FR-002** The result table MUST NOT be sent to the model whole. 2,679 rows + is the measured size of one small run; the model receives a bounded top-N. +- **FR-003** Results sent to the model MUST go through the existing + `src/analysis/disclosure.py` allow-list, **extended for ReactomeGSA's own + free-text fields**. The current `NEVER_SENT` names `fileName`, `sampleName` + and `columnNames`, which are the Analysis Service's. GSA carries the same + hazard under different names — `datasets[].name` is chosen by the user, and + `design.samples` is their column headers, which in the measured example were + `patient.id` values. An allow-list is wrong only by omission; these must be + omitted deliberately, not by luck. +- **FR-004** The existing warning MUST be shown before any result reaches the + model, and declining MUST still yield the file (Story 3). +- **FR-005** A successful `POST /analysis` MUST be treated as *accepted*, never + as *succeeded*. Measured: a submission returned HTTP 200 and then failed with + `CONNECTION_FORCED - broker forced connection closure` while the service was + mid-upgrade. The failure was visible only through `GET /status`. +- **FR-006** Uploads MUST be capped well below Chainlit's default. The config + ships `max_size_mb = 500`; the host has **4.8 GB free of 88 GB (95% used)**, + and `~/update-beta-chat.sh` already needs 6 GB to deploy. A handful of + default-sized uploads takes down the chat and blocks the next deploy. ~20 MB + covers a real matrix. +- **FR-007** An uploaded file MUST be deleted once submitted, whether the + analysis succeeds or fails. +- **FR-008** The run takes minutes, so the flow MUST NOT be shaped like a + streamed answer. The user is told it started and told again when it finishes. + +## Out of scope + +- Hosting any of this in reactome-mcp. The MCP is public, stateless and was + deliberately narrowed this month; file handling and long-running jobs belong + in the chat application. The five catalogue tools stay where they are. +- Single-cell clustering parameters (`k` for Single Cell Expression Atlas). +- Report generation: `GET /report_status/{id}` returned 404 for a completed + analysis, so whatever produces reports is not reachable this way and needs + its own investigation. + +## Open questions + +- Which method to default to. PADOG is the service's own default and is what + was measured; Camera is faster. Needs a decision informed by run time on a + realistic dataset, not by preference. +- Whether an analysis should survive a chat session ending. Results live + behind an analysis ID at the service, so resuming is possible; whether it is + wanted is a product question. diff --git a/src/gsa/__init__.py b/src/gsa/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/gsa/client.py b/src/gsa/client.py new file mode 100644 index 0000000..5531830 --- /dev/null +++ b/src/gsa/client.py @@ -0,0 +1,306 @@ +"""Talk to ReactomeGSA, the service that runs gene set analysis. + +A different service from the Analysis Service, doing a different thing: + + AnalysisService over-representation over a LIST of identifiers + ReactomeGSA gene set analysis over an EXPRESSION MATRIX + +Everything asserted here was measured against the live service on +2026-09-21, because the swagger is silent on the parts that matter. + +**The matrix is always inline.** `POST /analysis` requires +`datasets[].data`, the whole tab-delimited matrix as a string. There is no +by-reference variant, not even for a dataset the service itself just loaded. +Measured: 1.2 MB for the 16-sample melanoma example, a 1.5 MB submission, +and a 2.0 MB result holding 2,679 pathways. Nothing of that size may reach +the model or cross an MCP tool call, so this module exists to keep it on the +server. + +**A public dataset needs no upload.** `POST /data/load/{resourceId}` takes an +identifier -- Expression Atlas, Single Cell Expression Atlas, GREIN, GEO, or +the bundled examples -- and `GET /data/summary/{id}` then returns the sample +IDs and factors a user needs in order to name a comparison. + +**A 200 from /analysis means accepted, not succeeded.** Measured: a +submission returned 200 and then failed with `CONNECTION_FORCED - broker +forced connection closure` while the service was updating its Reactome +version. The failure was visible only through `GET /status`. Treat the +submission as a receipt and the status as the truth. + +**No browser user-agent here, unlike `analysis/client.py`.** That one needs +one because reactome.org sits behind automation blocking that answers a +library client with 403. `gsa.reactome.org` serves hypercorn directly and +answers `python-httpx` with 200 -- verified rather than assumed, because the +opposite mistake (testing only with curl, which is exempt from that block) +is how a sibling service was published and believed to work. +""" + +import os +import re +from dataclasses import dataclass +from typing import Any + +import httpx + +from util.logging import logging + +logger = logging.getLogger(__name__) + +BASE_URL_ENV = "REACTOME_GSA_URL" +DEFAULT_BASE_URL = "https://gsa.reactome.org/0.1" + +#: Submitting carries the matrix and the service answers only once it is +#: queued, so it needs far longer than a catalogue read. +TIMEOUT_SECONDS = 30.0 +SUBMIT_TIMEOUT_SECONDS = 180.0 + +#: Identifiers are interpolated into URL paths. A dataset ID comes from the +#: model, which means it comes from the user, so it is constrained to what +#: real ones look like -- `EXAMPLE_MEL_RNA`, `GSE12345`, `E-MTAB-2770`. +#: Anything else is rejected rather than escaped: a value containing `/` +#: addresses a different endpoint, which is a bypass rather than a 404. +_IDENTIFIER = re.compile(r"^[A-Za-z0-9][A-Za-z0-9_.-]{0,127}$") + + +class GsaError(RuntimeError): + """The service refused, or answered in a shape this code cannot use.""" + + +class GsaNotReadyError(GsaError): + """The analysis has not finished, so there is no result yet. + + `GET /result` answers 406 for this, which is an expected state on the + happy path rather than a fault. Collapsing it into a generic failure + would make "still running" indistinguishable from "broken", and the + caller polls on exactly that distinction. + """ + + +def _identifier_from(payload: Any, kind: str) -> str: + """Read an ID out of a response body, refusing anything else. + + Both POSTs answer with the bare ID as a quoted JSON string. `str()` on + an unexpected shape -- an error object, a list -- produces a plausible + string like `{'detail': ...}` that is then used as a path segment and + fails somewhere far from here, as a 404 that looks like a missing + analysis rather than a malformed response. + """ + if not isinstance(payload, str): + raise GsaError( + f"expected a {kind}, got {type(payload).__name__}: {payload!r:.120}" + ) + value = payload.strip() + if not _IDENTIFIER.match(value): + raise GsaError(f"{kind} is not a valid identifier: {value!r:.120}") + return value + + +def base_url() -> str: + return os.environ.get(BASE_URL_ENV, DEFAULT_BASE_URL).rstrip("/") + + +def _checked(kind: str, value: str) -> str: + if not _IDENTIFIER.match(value): + raise GsaError(f"{kind} is not a valid identifier: {value!r}") + return value + + +@dataclass(frozen=True) +class LoadingStatus: + """Progress of `POST /data/load`, which is not instant.""" + + status: str + description: str + completed: float + dataset_id: str | None + + @property + def finished(self) -> bool: + return self.status in {"complete", "failed"} + + +@dataclass(frozen=True) +class DatasetSummary: + """What a user needs in order to name a comparison. + + Small enough to show: the measured example is 16 samples and three + factors. The matrix it describes is 1.2 MB and is not here. + """ + + dataset_id: str + title: str + type: str + samples: list[str] + #: factor name -> that factor's value for each sample, in sample order + factors: dict[str, list[str]] + + def groups(self, factor: str) -> list[str]: + return sorted(set(self.factors.get(factor, []))) + + +@dataclass(frozen=True) +class AnalysisStatus: + status: str + description: str + completed: float + + @property + def finished(self) -> bool: + return self.status in {"complete", "failed"} + + @property + def failed(self) -> bool: + return self.status == "failed" + + +class GsaClient: + """Thin async client. Holds no state about a running analysis.""" + + def __init__(self, url: str | None = None) -> None: + self._base = (url or base_url()).rstrip("/") + + async def _get(self, path: str, *, timeout: float = TIMEOUT_SECONDS) -> Any: + async with httpx.AsyncClient(timeout=timeout) as client: + response = await client.get(f"{self._base}{path}") + if response.status_code == 406: + raise GsaNotReadyError(f"GET {path}: analysis is not complete") + if response.status_code != 200: + raise GsaError(f"GET {path} returned {response.status_code}") + return response.json() + + async def methods(self) -> list[dict[str, Any]]: + result = await self._get("/methods") + return list(result) if isinstance(result, list) else [] + + async def load_public_dataset(self, resource_id: str, dataset_id: str) -> str: + """Ask the service to load a public dataset. Returns a loading ID. + + The body is a list of name/value parameters rather than an object -- + the shape the service wants, confirmed by a successful load. + """ + _checked("resource", resource_id) + _checked("dataset", dataset_id) + body = [{"name": "dataset_id", "value": dataset_id}] + async with httpx.AsyncClient(timeout=TIMEOUT_SECONDS) as client: + response = await client.post( + f"{self._base}/data/load/{resource_id}", json=body + ) + if response.status_code not in (200, 202): + raise GsaError( + f"loading {dataset_id} from {resource_id} returned " + f"{response.status_code}: {response.text[:200]}" + ) + return _identifier_from(response.json(), "loading id") + + async def loading_status(self, loading_id: str) -> LoadingStatus: + data = await self._get(f"/data/status/{_checked('loading id', loading_id)}") + return LoadingStatus( + status=str(data.get("status", "")), + description=str(data.get("description", "")), + completed=float(data.get("completed") or 0.0), + dataset_id=data.get("dataset_id"), + ) + + async def dataset_summary(self, dataset_id: str) -> DatasetSummary: + data = await self._get(f"/data/summary/{_checked('dataset', dataset_id)}") + factors = { + str(entry.get("name")): [str(v) for v in entry.get("values", [])] + for entry in data.get("sample_metadata") or [] + if entry.get("name") + } + return DatasetSummary( + dataset_id=dataset_id, + title=str(data.get("title") or dataset_id), + type=str(data.get("type") or ""), + samples=[str(s) for s in data.get("sample_ids") or []], + factors=factors, + ) + + async def download_matrix(self, dataset_id: str) -> str: + """The expression matrix, as a tab-delimited string. + + Measured at 1.2 MB for a 16-sample dataset. It exists only to be + handed straight back to `submit`; it must not be logged, returned to + a caller that could show it, or put in a prompt. `format=expr` is + required -- omitting it is a 400, and `tsv` is not one of the + accepted values. + """ + url = f"{self._base}/data/download/{_checked('dataset', dataset_id)}" + async with httpx.AsyncClient(timeout=SUBMIT_TIMEOUT_SECONDS) as client: + response = await client.get(url, params={"format": "expr"}) + if response.status_code != 200: + raise GsaError(f"downloading {dataset_id} returned {response.status_code}") + logger.info( + "gsa matrix downloaded", + extra={"dataset": dataset_id, "bytes": len(response.content)}, + ) + return response.text + + async def submit( + self, + *, + method: str, + dataset_name: str, + dataset_type: str, + matrix: str, + samples: list[str], + analysis_group: list[str], + group1: str, + group2: str, + ) -> str: + """Start an analysis. Returns an analysis ID. + + The ID is a *receipt*. The analysis may still fail, and it will say + so through `analysis_status`, not here. + """ + body = { + "methodName": method, + "datasets": [ + { + "name": dataset_name, + "type": dataset_type, + "data": matrix, + "design": { + "samples": samples, + "analysisGroup": analysis_group, + "comparison": {"group1": group1, "group2": group2}, + }, + } + ], + } + async with httpx.AsyncClient(timeout=SUBMIT_TIMEOUT_SECONDS) as client: + response = await client.post(f"{self._base}/analysis", json=body) + if response.status_code != 200: + raise GsaError( + f"submitting returned {response.status_code}: {response.text[:200]}" + ) + analysis_id = _identifier_from(response.json(), "analysis id") + logger.info( + "gsa analysis submitted", + extra={ + "analysis": analysis_id, + "method": method, + "matrix_bytes": len(matrix), + }, + ) + return analysis_id + + async def analysis_status(self, analysis_id: str) -> AnalysisStatus: + data = await self._get(f"/status/{_checked('analysis', analysis_id)}") + return AnalysisStatus( + status=str(data.get("status", "")), + description=str(data.get("description", "")), + completed=float(data.get("completed") or 0.0), + ) + + async def result(self, analysis_id: str) -> dict[str, Any]: + """The finished result: ~2 MB, 2,679 pathways in the measured run. + + Returned whole because the caller needs the table to write a file. + It must be bounded before any of it reaches a model. + """ + data = await self._get( + f"/result/{_checked('analysis', analysis_id)}", + timeout=SUBMIT_TIMEOUT_SECONDS, + ) + return data if isinstance(data, dict) else {} diff --git a/src/gsa/results.py b/src/gsa/results.py new file mode 100644 index 0000000..e836513 --- /dev/null +++ b/src/gsa/results.py @@ -0,0 +1,229 @@ +"""Turn a 2 MB ReactomeGSA result into a small summary and a file. + +Two consumers with opposite needs. The **file** wants everything: 2,679 +pathways in the measured run, every column, so a researcher can sort and +filter it themselves. The **model** wants almost nothing: a bounded handful +of the most significant rows, because 2,679 rows would bury the answer, cost +more than the analysis, and say nothing the top twenty do not. + +So nothing here returns "the result". Each function names its audience. + +What may reach the model is an **allow-list**, for the reason +`analysis/disclosure.py` gives: a denial list is wrong by default the moment +the service adds a field, an allow-list is wrong only by omission. The +fields excluded here are not hypothetical -- a GSA submission carries the +user's own dataset name and their column headers, and the measured example's +headers were patient identifiers (`P1`..`P4` under a `patient.id` factor). +""" + +from dataclasses import dataclass, field +from typing import Any + +#: Columns of the pathway table that may be shown to a model. Everything +#: else -- present or future -- is excluded. `Pathway` and `Name` identify a +#: Reactome pathway, the rest are statistics over it; none is user content. +MODEL_COLUMNS: tuple[str, ...] = ( + "Pathway", + "Name", + "Direction", + "FDR", + "PValue", + "NGenes", +) + +#: Never sent to a model, under any circumstances, and named so a test can +#: assert on them. +#: +#: These are ReactomeGSA's own free-text carriers, and they are *different +#: field names* from the Analysis Service's `fileName` / `sampleName` / +#: `columnNames`. Reusing that list without adding these would have looked +#: like protection and provided none: +#: +#: datasets[].name chosen by the user, e.g. `smith_lab_unpublished` +#: design.samples the user's column headers +#: fold_changes one row per gene, and the columns are the samples +#: mappings the user's own row identifiers, mapped to UniProt -- +#: 8,035 entries and 499 KB in the measured run, and +#: the field nobody would have thought to exclude +NEVER_SENT: tuple[str, ...] = ("fold_changes", "design", "samples", "mappings") + +#: A result can hold thousands of pathways. This bounds what is described, +#: not what is saved. +DEFAULT_TOP = 20 + + +@dataclass(frozen=True) +class Pathway: + stable_id: str + name: str + direction: str + fdr: float + p_value: float + gene_count: int + + +@dataclass(frozen=True) +class GsaResult: + """A parsed result. `pathways` is the whole table; bound it before use.""" + + method: str + release: str + dataset_name: str + pathways: list[Pathway] + #: The service's own Pathway Browser view. Safe to show: it is a URL the + #: service minted, and it is how a user sees the result properly. + browser_links: list[tuple[str, str]] + #: The raw table, kept verbatim for the file so a researcher gets every + #: column rather than the six a model is allowed. + #: + #: `repr=False` is load-bearing. A dataclass's generated `__repr__` + #: includes every field, and this one is ~500 KB in a real run, so a + #: `logger.debug("%s", result)`, an exception context, or a failing + #: test's output would print the entire pathway table into a log. The + #: care taken over what reaches a *model* is wasted if the same content + #: reaches a *log file* by default. + raw_table: str = field(repr=False) + + @property + def significant(self) -> list[Pathway]: + return [p for p in self.pathways if p.fdr < 0.05] + + +def _as_float(value: str) -> float: + try: + return float(value) + except (TypeError, ValueError): + # A non-numeric FDR must not sort as "most significant". The + # service has been seen returning "NA". + return 1.0 + + +def parse(result: dict[str, Any]) -> GsaResult: + """Parse the service's result. Does no bounding -- see `for_model`.""" + datasets = result.get("results") or [] + first = datasets[0] if datasets else {} + table = first.get("pathways") + table = table if isinstance(table, str) else "" + + pathways: list[Pathway] = [] + lines = [line for line in table.splitlines() if line.strip()] + if lines: + header = lines[0].split("\t") + index = {name: i for i, name in enumerate(header)} + + def cell(row: list[str], column: str) -> str: + position = index.get(column) + return row[position] if position is not None and position < len(row) else "" + + for line in lines[1:]: + row = line.split("\t") + try: + genes = int(float(cell(row, "NGenes") or 0)) + except ValueError: + genes = 0 + pathways.append( + Pathway( + stable_id=cell(row, "Pathway"), + name=cell(row, "Name"), + direction=cell(row, "Direction"), + fdr=_as_float(cell(row, "FDR")), + p_value=_as_float(cell(row, "PValue")), + gene_count=genes, + ) + ) + + links = [ + (str(link.get("name") or "Reactome"), str(link.get("url"))) + for link in result.get("reactome_links") or [] + if link.get("url") + ] + + return GsaResult( + # `method_name`, not `methodName`. The swagger's AnalysisResult + # definition says the latter; the service returns the former, so + # this field read None until a fixture built from a real response + # showed it. Same defect as the ten tools reactome-mcp fixed in + # September: a field path asserted rather than verified. + method=str(result.get("method_name") or result.get("methodName") or ""), + release=str(result.get("release") or ""), + dataset_name=str(first.get("name") or ""), + pathways=pathways, + browser_links=links, + raw_table=table, + ) + + +def for_model(result: GsaResult, *, top: int = DEFAULT_TOP) -> dict[str, Any]: + """The bounded, allow-listed view that may go into a prompt. + + Counts are exact and stated as such, so a summary can say "2,679 + pathways, 412 significant" without having seen 2,679 of anything. + """ + ranked = sorted(result.pathways, key=lambda p: (p.fdr, p.p_value))[:top] + + # A result with no pathways is not an analysis that found nothing -- it + # is an analysis whose table did not arrive, which is what a failed or + # half-written result looks like. Reporting `0` with + # `counts_are_exact: True` invites a model to tell the user their data + # contained no enriched pathways, which is a confident answer to a + # question nobody managed to ask. Say "no result" instead. + if not result.pathways: + return { + "release": result.release, + "no_result": True, + "counts_are_exact": False, + "pathway_count": 0, + "significant_count": 0, + "showing": 0, + "top_pathways": [], + } + + return { + "release": result.release, + "no_result": False, + "pathway_count": len(result.pathways), + "significant_count": len(result.significant), + "counts_are_exact": True, + "showing": len(ranked), + "top_pathways": [ + { + "stId": p.stable_id, + "name": p.name, + "direction": p.direction, + "fdr": p.fdr, + "genes": p.gene_count, + } + for p in ranked + ], + } + + +def for_user(result: GsaResult) -> list[tuple[str, str]]: + """Links for the chat to show the user. **Not for the model.** + + The Pathway Browser URL embeds the analysis token: + + https://reactome.org/PathwayBrowser/#/DTAB=AN&ANALYSIS=MjAyNj... + + Anyone holding that token can fetch the whole result back from the + service -- including `mappings`, which is the user's own gene + identifiers, and `fold_changes`, whose columns are their samples. So + the link is a *capability*, not a citation. + + `for_model` used to include it. Stripping user content from the payload + and then handing over a key that retrieves it is not protection; it is + the same content by a longer route. The user should absolutely see this + link -- it is how they view their own result properly -- so the chat + renders it and the model never receives it. + """ + return list(result.browser_links) + + +def as_tsv(result: GsaResult) -> str: + """The whole table, for the file the user downloads. + + Verbatim rather than reassembled from `Pathway` objects: the parsed form + keeps six columns and the service returned nine, and a researcher asked + for their results should get their results. + """ + return result.raw_table diff --git a/tests/gsa/result_fixture.json b/tests/gsa/result_fixture.json new file mode 100644 index 0000000..0cbd413 --- /dev/null +++ b/tests/gsa/result_fixture.json @@ -0,0 +1,33 @@ +{ + "release": "97", + "method_name": "padog", + "results": [ + { + "name": "melanoma", + "pathways": "Pathway\tName\tDirection\tFDR\tPValue\tNGenes\tMeanAbsT0\tMeanWeightT0\tav_foldchange\nR-HSA-176417\tPhosphorylation of Emi1\tDown\t1e-05\t0.001\t6\t5.8688267847137405\t5.0890032000411685\t-0.823247774951112\nR-HSA-69481\tG2/M Checkpoints\tDown\t1e-05\t0.006\t132\t6.223467132558967\t2.530308576586947\t-0.30499714859719396\nR-HSA-68881\tMitotic Metaphase/Anaphase Transition\tDown\t0.001\t0.008\t2\t4.045413186464432\t3.892355916421484\t-0.9074415514449177\nR-HSA-9839923\tDengue Virus Infection\tDown\t0.001\t0.003\t327\t4.906607633776362\t3.0380303853953463\t-0.157866405311507\nR-HSA-3700989\tTranscriptional Regulation by TP53\tDown\t0.001\t0.008\t342\t3.7316284909879025\t2.0425949256767058\t-0.04056302612315259\nR-HSA-9609690\tHCMV Early Events\tDown\t0.001\t0.014\t83\t4.520669578004483\t0.32774621918271796\t-0.3689377510074243", + "fold_changes": "GENE\tPATIENT_001_TUMOUR\nENSG00000000419\t-0.8231\n" + } + ], + "reactome_links": [ + { + "url": "https://reactome.org/PathwayBrowser/#/DTAB=AN&ANALYSIS=MjAyNjA5MjExOTMyNTlfNjExMTE%3D", + "name": "Gene Set Analysis Summary", + "token": "MjAyNjA5MjExOTMyNTlfNjExMTE%3D", + "description": "Overview over all submitted datasets showing significantly and non-significantly up- and down- regulated pathways" + } + ], + "mappings": [ + { + "identifier": "ENSG00000131089", + "mapped_to": [ + "O43307" + ] + }, + { + "identifier": "PATIENT_001_GENE", + "mapped_to": [ + "P00001" + ] + } + ] +} \ No newline at end of file diff --git a/tests/gsa/test_gsa_client.py b/tests/gsa/test_gsa_client.py new file mode 100644 index 0000000..5441224 --- /dev/null +++ b/tests/gsa/test_gsa_client.py @@ -0,0 +1,72 @@ +"""The client's handling of shapes the service actually produces. + +No network here. Everything asserted was seen from the live service on +2026-09-21 and is reproduced as a stub, so the suite keeps testing it after +the service moves on. +""" + +import pytest + +from gsa import client as gsa_client + + +@pytest.mark.parametrize( + "value", + ["../etc/passwd", "a/b", "GSE1 2345", "", "x" * 200, "tok;drop"], +) +def test_hostile_identifiers_are_refused(value: str) -> None: + # These are interpolated into a URL path. A value containing `/` + # addresses a different endpoint, which is a bypass, not a 404. + with pytest.raises(gsa_client.GsaError): + gsa_client._checked("dataset", value) + + +@pytest.mark.parametrize("value", ["EXAMPLE_MEL_RNA", "GSE12345", "E-MTAB-2770"]) +def test_real_identifiers_are_accepted(value: str) -> None: + assert gsa_client._checked("dataset", value) == value + + +def test_an_id_response_must_be_a_string() -> None: + # The service answers a submission with a bare quoted ID. `str()` on an + # error object yields a plausible-looking path segment that fails far + # away, as a 404 that reads like a missing analysis. + with pytest.raises(gsa_client.GsaError, match="expected a analysis id"): + gsa_client._identifier_from({"detail": "Bad Request"}, "analysis id") + + +def test_a_real_id_response_is_accepted() -> None: + assert ( + gsa_client._identifier_from( + "d8f825e8-b5f2-11f1-9de2-863bc094cc6c", "analysis id" + ) + == "d8f825e8-b5f2-11f1-9de2-863bc094cc6c" + ) + + +def test_status_finished_covers_failure_not_just_success() -> None: + # Measured: a submission returned 200 and then reached + # status=failed with "CONNECTION_FORCED - broker forced connection + # closure". A poll loop that waits for "complete" alone never exits. + running = gsa_client.AnalysisStatus("running", "Permutation 900 / 1000", 0.6) + failed = gsa_client.AnalysisStatus("failed", "Failed to analyse dataset", 1.0) + done = gsa_client.AnalysisStatus("complete", "Analysis done", 1.0) + + assert not running.finished + assert failed.finished + assert failed.failed + assert done.finished + assert not done.failed + + +def test_summary_exposes_the_groups_a_user_must_choose_between() -> None: + # The real shape: sample_metadata is a list of named factors, each with + # one value per sample. + summary = gsa_client.DatasetSummary( + dataset_id="EXAMPLE_MEL_RNA", + title="Melanoma RNA-seq example", + type="rnaseq_counts", + samples=[f"S{i}" for i in range(4)], + factors={"condition": ["MCM", "MOCK", "MCM", "MOCK"]}, + ) + assert summary.groups("condition") == ["MCM", "MOCK"] + assert summary.groups("no-such-factor") == [] diff --git a/tests/gsa/test_gsa_results.py b/tests/gsa/test_gsa_results.py new file mode 100644 index 0000000..f857c8d --- /dev/null +++ b/tests/gsa/test_gsa_results.py @@ -0,0 +1,193 @@ +"""What reaches the model, and what only reaches the file. + +The fixture is a real ReactomeGSA response, trimmed: the same keys, the same +column header, real rows, and `mappings` cut from 8,035 entries to two that +keep its hazard. Trimmed rather than invented, because the two defects found +while writing this code were both invisible to an invented one -- the result +key is `method_name` where the swagger says `methodName`, and `mappings` +carries the user's own row identifiers. +""" + +import json +from pathlib import Path +from typing import Any + +import pytest + +from gsa import results as gsa_results + +FIXTURE = Path(__file__).parent / "result_fixture.json" + + +@pytest.fixture +def raw() -> dict[str, Any]: + loaded: dict[str, Any] = json.loads(FIXTURE.read_text()) + return loaded + + +@pytest.fixture +def parsed(raw: dict[str, Any]) -> gsa_results.GsaResult: + return gsa_results.parse(raw) + + +def test_parses_the_service_shape(parsed: gsa_results.GsaResult) -> None: + assert parsed.release == "97" + # The assertion that would have failed against the swagger's field name. + assert parsed.method == "padog" + assert len(parsed.pathways) == 6 + first = parsed.pathways[0] + assert first.stable_id.startswith("R-HSA-") + assert first.name + assert first.direction in {"Up", "Down"} + assert 0.0 <= first.fdr <= 1.0 + + +def test_browser_link_survives(parsed: gsa_results.GsaResult) -> None: + # The user's real way of seeing the result. Losing it silently would + # leave them with a table and no picture. + assert parsed.browser_links + assert all(url.startswith("http") for _, url in parsed.browser_links) + + +def test_model_view_is_bounded(parsed: gsa_results.GsaResult) -> None: + view = gsa_results.for_model(parsed, top=3) + assert view["showing"] == 3 + assert view["pathway_count"] == 6 + assert view["counts_are_exact"] is True + # Bounded on what is shown, exact on what is counted: a summary may say + # how many there were without having seen them. + assert len(view["top_pathways"]) == 3 + + +def test_model_view_is_ranked_by_significance(parsed: gsa_results.GsaResult) -> None: + view = gsa_results.for_model(parsed, top=6) + fdrs = [p["fdr"] for p in view["top_pathways"]] + assert fdrs == sorted(fdrs) + + +def test_a_non_numeric_fdr_does_not_rank_first(raw: dict[str, Any]) -> None: + # Measured behaviour elsewhere in Reactome: "NA" appears in numeric + # columns. float() raises, and a default of 0.0 would sort it to the top + # as the most significant result in the analysis. + lines = raw["results"][0]["pathways"].split("\n") + header = lines[0].split("\t") + row = lines[1].split("\t") + row[header.index("FDR")] = "NA" + lines[1] = "\t".join(row) + raw["results"][0]["pathways"] = "\n".join(lines) + + parsed = gsa_results.parse(raw) + ranked = gsa_results.for_model(parsed, top=6)["top_pathways"] + assert ranked[0]["fdr"] < 1.0 + assert ranked[-1]["fdr"] == 1.0 + + +@pytest.mark.parametrize("field", gsa_results.NEVER_SENT) +def test_never_sent_fields_are_absent_by_name( + parsed: gsa_results.GsaResult, field: str +) -> None: + assert field not in json.dumps(gsa_results.for_model(parsed)) + + +def test_no_user_content_reaches_the_model(raw: dict[str, Any]) -> None: + """The test that matters: a marker planted in every user-supplied place + must not appear anywhere in the model's view. + + Asserting on field *names* is not enough -- a field could be copied into + a differently named one. This asserts on the values. + """ + marker = "SMITH-LAB-UNPUBLISHED-2026" + raw["results"][0]["name"] = marker + raw["results"][0]["fold_changes"] = f"GENE\t{marker}\nENSG1\t1.0\n" + raw["mappings"] = [{"identifier": marker, "mapped_to": ["P00001"]}] + + view = json.dumps(gsa_results.for_model(gsa_results.parse(raw))) + assert marker not in view + + +def test_the_marker_test_can_fail(raw: dict[str, Any]) -> None: + """The control for the test above. + + An absence assertion proves nothing until the same construction has been + shown capable of producing a presence. If the marker cannot reach the + view even when deliberately placed in an allow-listed field, the test + above is passing vacuously. + """ + marker = "SMITH-LAB-UNPUBLISHED-2026" + lines = raw["results"][0]["pathways"].split("\n") + header = lines[0].split("\t") + row = lines[1].split("\t") + row[header.index("Name")] = marker + lines[1] = "\t".join(row) + raw["results"][0]["pathways"] = "\n".join(lines) + + assert marker in json.dumps(gsa_results.for_model(gsa_results.parse(raw))) + + +def test_the_file_keeps_every_column(parsed: gsa_results.GsaResult) -> None: + # The model gets six columns; the researcher asked for their results and + # the service returned nine. + header = gsa_results.as_tsv(parsed).splitlines()[0].split("\t") + assert set(gsa_results.MODEL_COLUMNS).issubset(header) + assert {"MeanAbsT0", "MeanWeightT0", "av_foldchange"}.issubset(header) + + +def test_the_result_does_not_print_its_table(parsed: gsa_results.GsaResult) -> None: + """A dataclass repr includes every field by default. + + The table is ~500 KB in a real run, so one `logger.debug("%s", result)` + or one exception context would put the whole thing in a log. Care over + what reaches a model is wasted if the same content reaches a log file + for free. + """ + printed = repr(parsed) + assert "raw_table" not in printed + assert "MeanWeightT0" not in printed # a column only the raw table has + + +def test_an_empty_result_is_not_reported_as_an_exact_zero() -> None: + """The failure mode that reads like an answer. + + A result whose table did not arrive has no pathways. Reporting that as + `pathway_count: 0, counts_are_exact: True` invites a summary saying the + user's data contained no enriched pathways -- a confident answer to a + question nobody managed to ask. + """ + empty = gsa_results.parse({"release": "97", "method_name": "padog", "results": []}) + view = gsa_results.for_model(empty) + + assert view["no_result"] is True + assert view["counts_are_exact"] is False + + +def test_a_real_result_is_not_flagged_as_missing(parsed: gsa_results.GsaResult) -> None: + # The control: without it, `no_result: True` on everything would pass + # the test above and break the feature. + view = gsa_results.for_model(parsed) + assert view["no_result"] is False + assert view["counts_are_exact"] is True + assert view["pathway_count"] == 6 + + +def test_the_analysis_token_never_reaches_the_model( + parsed: gsa_results.GsaResult, +) -> None: + """The Pathway Browser link is a capability, not a citation. + + Its URL embeds the analysis token, and anyone holding that token can + fetch the whole result back -- including `mappings`, the user's own gene + identifiers. Stripping user content from the payload and then sending a + key that retrieves it is the same disclosure by a longer route. + """ + view = json.dumps(gsa_results.for_model(parsed)) + assert "ANALYSIS=" not in view + assert "PathwayBrowser" not in view + assert "reactome.org" not in view + + +def test_the_user_still_gets_the_link(parsed: gsa_results.GsaResult) -> None: + # Withholding it from the model must not withhold it from the person. + # It is how they see their own result properly. + links = gsa_results.for_user(parsed) + assert links + assert any("PathwayBrowser" in url for _, url in links)