diff --git a/database/scripts/common.py b/database/scripts/common.py index cd609fd..2101524 100644 --- a/database/scripts/common.py +++ b/database/scripts/common.py @@ -17,7 +17,7 @@ from retromol.model.result import Result from retromol.model.rules import MatchingRule, RuleSet from retromol_database.duckdb import FINGERPRINT_SIZE -from retromol_fingerprint.fingerprint import Fingerprinter, Vocabulary +from retromol_fingerprint.fingerprint import TOKEN_LINK, Fingerprinter, Vocabulary # Silences RDKit's kekulization/valence/etc. warnings in *this* (single) process -- # every pipeline script imports common, so this alone covers create_db.py, @@ -69,7 +69,16 @@ def build_fingerprint_context(ruleset: RuleSet) -> tuple[dict[str, MatchingRule] def per_monomer_tokens(name: str, name_to_rule: dict[str, MatchingRule]) -> list[str]: - """Fingerprinting token list for one compound primary-sequence block (mirrors discovery.py's `_per_monomer_tokens`).""" + """Fingerprinting token list for one compound primary sequence block (mirrors discovery.py's `_per_monomer_tokens`). + + TOKEN_LINK is not a building block -- it just joins two merged paths -- and is + filtered out by callers before fingerprinting (see load_compounds.py), so it's + never actually looked up here. Handled explicitly anyway so a stray call can't + silently fall through to the empty-tokens/TOKEN_UNK path instead. + """ + if name == TOKEN_LINK: + return [TOKEN_LINK] + if name in PK_GROUP_TOKENS: return [name, "PK"] @@ -103,53 +112,26 @@ def npatlas_url(npaid: str | None) -> str | None: return NPATLAS_URL_TEMPLATE.format(npaid=npaid) -def primary_sequences_from_result(result: Result, min_length: int = 2) -> tuple[list[list[str]], list[str]]: +def primary_sequence_from_result(result: Result) -> list[str]: """ - Every candidate primary sequence for a parsed compound, read directly off - `result.linear_readout.paths` -- no backbone reconstruction involved, that's a + The single primary sequence for a parsed compound, read directly off + `result.linear_readout` -- no backbone reconstruction involved, that's a display-only concern this pipeline has no use for. `result.linear_readout` is already computed by retromol.pipelines.parsing.run_retromol (it's just a field - on Result), each path is one candidate ordering of monomers through the - molecule, and each becomes its own db entry (see load_compounds.py). An - unidentified node is named "X", the same convention used everywhere else in - RetroMol. - - `result.linear_readout` also includes single-node paths for tailoring events - that don't connect to any chain -- most commonly glycosylation (AssemblyGraph - only keeps C-C/C-N bonds as "connections", so a sugar attached via a glycosidic - C-O-C linkage always ends up disconnected from the main chain) and methylation. - These aren't a "sequence" in any useful sense, so they're excluded from the - returned primary_sequences by the `min_length` floor -- but they're still real, - identified content RetroMol found on this molecule, so their names are returned - separately as `extra_tokens`, meant to be folded into the fingerprint alongside - each primary_sequence (see per_monomer_tokens) without being displayed as part - of it. Dropping them from the fingerprint entirely -- the previous behavior -- - made two compounds differing only in glycosylation pattern (a common source of - real biosynthetic diversity, e.g. erythromycin vs. megalomicin) invisible to - that difference. - - `extra_tokens` is NOT deduplicated: Fingerprinter.encode has no final - normalization step, so each occurrence adds its own weight -- a molecule with - two glycosylation events should end up with roughly twice the glycosylation - weight of one with a single event, not the same weight collapsed to "some - glycosylation happened". Deduping here would throw that count signal away. + on Result); every path through the molecule -- including single-node paths for + tailoring events that don't connect to any chain, e.g. glycosylation/methylation + (AssemblyGraph only keeps C-C/C-N bonds as "connections", so a sugar attached via + a glycosidic C-O-C linkage always ends up disconnected from the main chain) -- + is merged into one sequence (longest first, ties broken lexicographically, joined + by TOKEN_LINK -- see `LinearReadout.primary_sequence`), and that single sequence + becomes this compound's one db entry (see load_compounds.py). Nothing found in + the assembly graph is dropped. An unidentified node is named "X", the same + convention used everywhere else in RetroMol. :param result: a parsed RetroMol Result - :param min_length: drop paths shorter than this from primary_sequences (default 2) - :return: (primary_sequences, extra_tokens) -- one extra_tokens entry per - qualifying tailoring-event path, duplicates and all + :return: the single primary sequence """ - primary_sequences: list[list[str]] = [] - extra_tokens: list[str] = [] - - for path in result.linear_readout.paths: - names = [node.identity.matched_rule.name if node.is_identified else "X" for node in path] - if len(path) >= min_length: - primary_sequences.append(names) - else: - extra_tokens.extend(n for n in names if n != "X") - - return primary_sequences, extra_tokens + return result.linear_readout.primary_sequence() def _init_worker_quiet(ruleset: RuleSet) -> None: diff --git a/database/scripts/load_compounds.py b/database/scripts/load_compounds.py index fe561c2..645ffec 100644 --- a/database/scripts/load_compounds.py +++ b/database/scripts/load_compounds.py @@ -1,18 +1,11 @@ """Steps 4 & 6: turn parsed compound results into database entries. -For each RetroMol Result, every candidate primary sequence -- one per path in -result.linear_readout.paths, read directly off the Result (see -common.primary_sequences_from_result) -- becomes its own "compound" entry. -Ambiguous parses intentionally produce multiple queryable entries, the same -convention the webapp uses for an uploaded compound with more than one candidate -reading. `raw` is always the original input SMILES, the same for every entry a -given compound produces. - -A compound's tailoring events (glycosylation, methylation -- anything that shows up -as its own disconnected single-node path, see common.primary_sequences_from_result) -aren't part of any displayed primary_sequence, but their tokens still get folded -into every one of that compound's fingerprints, so two compounds that only differ -by e.g. an attached sugar aren't fingerprint-identical. +For each RetroMol Result, every path found in result.linear_readout -- including a +compound's tailoring events (glycosylation, methylation -- anything that shows up as +its own disconnected single-node path) -- is merged into the one primary sequence +stored for that compound: longest path first, ties broken lexicographically, joined +by TOKEN_LINK (see common.primary_sequence_from_result). One compound -> one db +entry, always. `raw` is the original input SMILES. """ import argparse @@ -30,10 +23,11 @@ mibig_url, npatlas_url, per_monomer_tokens, - primary_sequences_from_result, + primary_sequence_from_result, ) from retromol.model.result import Result from retromol_database.duckdb import RetroMolDuckDB +from retromol_fingerprint.fingerprint import TOKEN_LINK log = logging.getLogger(__name__) @@ -96,15 +90,15 @@ def run( name = name or result.submission.name or result.submission.inchikey - sequences, extra_tokens = primary_sequences_from_result(result) - extra_tokens_encoded = [per_monomer_tokens(n, name_to_rule) for n in extra_tokens] + names = primary_sequence_from_result(result) - for names in sequences: - if not names: - skipped += 1 - continue - - tokens = [per_monomer_tokens(n, name_to_rule) for n in names] + extra_tokens_encoded + if not names: + skipped += 1 + else: + # TOKEN_LINK only marks where two merged paths join -- it isn't + # a building block, so it's excluded from the fingerprint (an + # empty token list would otherwise silently add TOKEN_UNK mass). + tokens = [per_monomer_tokens(n, name_to_rule) for n in names if n != TOKEN_LINK] fp = fingerprinter.encode(tokens) db.add_entry( diff --git a/gui/src/client/src/components/MotifHoverCard.tsx b/gui/src/client/src/components/MotifHoverCard.tsx index 17bb851..10bdb88 100644 --- a/gui/src/client/src/components/MotifHoverCard.tsx +++ b/gui/src/client/src/components/MotifHoverCard.tsx @@ -3,6 +3,7 @@ import Box from "@mui/material/Box"; import CircularProgress from "@mui/material/CircularProgress"; import Tooltip from "@mui/material/Tooltip"; import Typography from "@mui/material/Typography"; +import WarningAmberIcon from "@mui/icons-material/WarningAmber"; import { useQuery } from "@tanstack/react-query"; import { fetchMotifStructures } from "../features/motifs/api"; import { MotifName } from "./MotifName"; @@ -23,7 +24,11 @@ function MotifHoverContent({ name, hint }: { name: string; hint?: string }) { gcTime: Infinity, }); - const smiles = structuresQuery.data?.[name]; + const smiles = structuresQuery.data?.structures[name]; + // Some names (e.g. "glycosylation") don't identify one specific rule -- they're + // shared by many rules with different structures, so whichever one `smiles` + // shows was picked arbitrarily (see /api/motifStructures's ambiguousNames). + const isAmbiguous = structuresQuery.data?.ambiguousNames.includes(name) ?? false; // Must be unique per *mounted instance*, not derived from `name` alone -- // SmilesDrawerContainer draws into a DOM node it looks up by this id, and the @@ -62,6 +67,15 @@ function MotifHoverContent({ name, hint }: { name: string; hint?: string }) { )} )} + {isAmbiguous && ( + + + + Multiple structures share this name: the one shown is picked arbitrarily, + not necessarily the one that actually matched here. + + + )} diff --git a/gui/src/client/src/components/MotifName.tsx b/gui/src/client/src/components/MotifName.tsx index 03562e8..cff0fae 100644 --- a/gui/src/client/src/components/MotifName.tsx +++ b/gui/src/client/src/components/MotifName.tsx @@ -1,8 +1,20 @@ import React from "react"; import Box from "@mui/material/Box"; +import { isLinkToken } from "../features/reconstruction/types"; // Renders a monomer name like "A2^R" or "C^E2", superscripting the trailing stereo marker. +// The link token (joining two merged primary sequence paths, see LINK_TOKEN) isn't a real +// building block -- render it as a muted glyph instead of a name, everywhere a monomer +// name would otherwise be rendered (sequence chips, alignment grid cells, ...). export function MotifName({ name }: { name: string }) { + if (isLinkToken(name)) { + return ( + + ⋮ + + ); + } + const parts = name.split(/(\^[SREZ])/g); return ( diff --git a/gui/src/client/src/components/workspace/DialogViewItem.tsx b/gui/src/client/src/components/workspace/DialogViewItem.tsx index 475086c..b02a571 100644 --- a/gui/src/client/src/components/workspace/DialogViewItem.tsx +++ b/gui/src/client/src/components/workspace/DialogViewItem.tsx @@ -15,7 +15,7 @@ import { ErrorBoundary } from "../ErrorBoundary"; import { ExportImageButton } from "../ExportImageButton"; import SmilesDrawerContainer from "../SmilesDrawerContainer.js"; import { DrawingAttribution } from "../DrawingAttribution"; -import { PrimarySequenceRows, usePrimarySequenceEditor } from "./PrimarySequenceEditor"; +import { PrimarySequenceOverview, PrimarySequenceRows, usePrimarySequenceEditor } from "./PrimarySequenceEditor"; import { ClusterReadoutRows } from "./ClusterReadoutRows"; type HighlightAtom = [number, string]; @@ -126,6 +126,7 @@ export const DialogViewItem: React.FC = ({ }); const data = reconstructionQuery.data?.reconstructions ?? null; + const dbMatchingSequences = reconstructionQuery.data?.dbMatchingSequences ?? []; const loading = reconstructionQuery.isLoading; const error = reconstructionQuery.error ? (reconstructionQuery.error as Error).message || "Unknown error" @@ -269,6 +270,17 @@ export const DialogViewItem: React.FC = ({ )} + + {dbMatchingSequences.length > 0 && ( + + + + )} + void; + renderAction?: (sequence: PrimarySequenceItem[]) => React.ReactNode; +}) { + return ( + <> + {data.map((reconstruction, idx) => { + const subsequences = splitOnLinkToken(reconstruction.primary_sequence); + return ( + + {/* pb reserves the same space PrimarySequenceChips' own mb below cancels out of + the row's reported height (see that mb's comment) -- without this, the row + as a whole under-reports its true rendered height by that same amount, and + whatever comes right after (a chain row, or anything a caller of this + component places after it) paints over the scrollbar instead of below it. */} + + + + + {/* mb cancels out PrimarySequenceChips' own reserved scrollbar space (see + horizontalScrollSx's pb) so it doesn't throw off alignItems: "center" + against the label chip -- without it, the chip row reads as sitting + above center, offset by however much space is reserved below it. */} + + + + {renderAction && {renderAction(reconstruction.primary_sequence)}} + + + {subsequences.length > 1 && ( + + {subsequences.map((subsequence, subIdx) => ( + + + + + + + + {renderAction && {renderAction(subsequence)}} + + ))} + + )} + + ); + })} + + ); +} + // Renders one row per reconstruction: a label (+ "Edited" chip), and either the // read-only chip row or the live SequenceEditor + per-row revert control. export function PrimarySequenceRows({ @@ -317,7 +412,7 @@ export function PrimarySequenceRows({ minWidth: 0, }} > - {ordered ? `primary sequence ${idx + 1}` : "parsed motifs (unordered)"} + {ordered ? `reconstructed backbone ${idx + 1}` : "parsed motifs (unordered)"} {override && !editing && ( diff --git a/gui/src/client/src/components/workspace/WorkspaceDiscovery.tsx b/gui/src/client/src/components/workspace/WorkspaceDiscovery.tsx index 63db5a1..a4c663e 100644 --- a/gui/src/client/src/components/workspace/WorkspaceDiscovery.tsx +++ b/gui/src/client/src/components/workspace/WorkspaceDiscovery.tsx @@ -44,6 +44,7 @@ import { useNotifications } from "../NotificationProvider"; import { MotifName } from "../MotifName"; import { horizontalScrollSx } from "../../theme/scrollbarSx"; import { SequenceEditor, type SequenceBlock } from "./SequenceEditor"; +import { PrimarySequenceOverview } from "./PrimarySequenceEditor"; import { DialogViewDiscoveryQuery } from "./DialogViewDiscoveryQuery"; import {MinimalIconButton} from "../MinimalIconButton"; @@ -372,11 +373,6 @@ export const WorkspaceDiscovery: React.FC = ({ session, const [selectedItemId, setSelectedItemId] = React.useState(""); const selectedItem = session.items.find((item) => item.id === selectedItemId); const [blocks, setBlocks] = React.useState([]); - // Carried alongside `blocks` from whichever candidate "Use this" was clicked on -- - // only ever non-empty when that candidate was a dbMatchingSequences entry (see - // handlePickNames). Folded into the query's fingerprint on submit without being - // shown as part of the editable sequence itself. - const [extraFingerprintTokens, setExtraFingerprintTokens] = React.useState([]); const [entryType, setEntryType] = React.useState("compound"); const [scoreMode, setScoreMode] = React.useState("longest_sequence"); @@ -420,9 +416,8 @@ export const WorkspaceDiscovery: React.FC = ({ session, enabled: selectedItem?.kind === "cluster", }); - const handlePickNames = (names: string[], extraTokens: string[] = []) => { + const handlePickNames = (names: string[]) => { setBlocks(blocksFromNames(names)); - setExtraFingerprintTokens(extraTokens); }; const maxTopX = Math.max(1, Math.min(n, MAX_TOP_X)); @@ -452,7 +447,6 @@ export const WorkspaceDiscovery: React.FC = ({ session, includeUserUploads: includeUserUploads || onlyUserUploads, onlyUserUploads, queryOriginSmiles, - extraFingerprintTokens: extraFingerprintTokens.length > 0 ? extraFingerprintTokens : undefined, flags: { computeMsa, computeCompare }, }); setSession((prev) => (prev ? { ...prev, items: [...prev.items, item] } : prev)); @@ -607,45 +601,25 @@ export const WorkspaceDiscovery: React.FC = ({ session, {reconstructionQuery.data && reconstructionQuery.data.dbMatchingSequences.length > 0 && ( <> - Database-matching sequences -- query with one of these for a fingerprint - guaranteed comparable to what's actually stored + Database-matching sequence -- query with this for a fingerprint + guaranteed comparable to what's actually stored. Nothing is + filtered out of it (tailoring events like glycosylation are + included), so a compound with more than one biosynthetic chain + can also be queried with just one of them below. - - {reconstructionQuery.data.dbMatchingSequences.map((reconstruction, idx) => ( - - - - - - - + + ( - - ))} + )} + /> )} @@ -673,7 +647,7 @@ export const WorkspaceDiscovery: React.FC = ({ session, }} > - {`primary sequence ${idx+1}`} + {`reconstructed backbone ${idx+1}`} {override && ( diff --git a/gui/src/client/src/components/workspace/WorkspaceItemCard.tsx b/gui/src/client/src/components/workspace/WorkspaceItemCard.tsx index db9de22..fa163ec 100644 --- a/gui/src/client/src/components/workspace/WorkspaceItemCard.tsx +++ b/gui/src/client/src/components/workspace/WorkspaceItemCard.tsx @@ -26,7 +26,7 @@ import { renameSessionItem } from "../../features/session/api"; import { alpha } from "@mui/material/styles"; import type { Theme } from "@mui/material/styles"; import { DialogViewItem } from "./DialogViewItem"; -import { PrimarySequenceRows, usePrimarySequenceEditor } from "./PrimarySequenceEditor"; +import { PrimarySequenceOverview, PrimarySequenceRows, usePrimarySequenceEditor } from "./PrimarySequenceEditor"; import { ClusterReadoutRows } from "./ClusterReadoutRows"; import { reconstructCompound } from "../../features/reconstruction/api"; import { getClusterReadout } from "../../features/clusters/api"; @@ -162,6 +162,7 @@ export const WorkspaceItemCard: React.FC = ({ enabled: expanded && isCompound, }); const reconstructions = reconstructionQuery.data?.reconstructions ?? null; + const dbMatchingSequences = reconstructionQuery.data?.dbMatchingSequences ?? []; const editor = usePrimarySequenceEditor(session, setSession, item, reconstructions); // See DialogViewItem's isUnordered -- an unordered "bag of motifs" result isn't a // sequence, so there's nothing meaningful to drag-and-drop reorder. @@ -493,6 +494,12 @@ export const WorkspaceItemCard: React.FC = ({ )} + {dbMatchingSequences.length > 0 && ( + + + + )} + {reconstructions && reconstructions.length > 0 && ( diff --git a/gui/src/client/src/components/workspace/alignmentSvgExport.ts b/gui/src/client/src/components/workspace/alignmentSvgExport.ts index 824a9c0..2878669 100644 --- a/gui/src/client/src/components/workspace/alignmentSvgExport.ts +++ b/gui/src/client/src/components/workspace/alignmentSvgExport.ts @@ -4,6 +4,8 @@ // (e.g. via html2canvas), so the output stays crisp at any zoom/print size and opens // cleanly in vector editors like Illustrator/Inkscape. +import { LINK_TOKEN } from "../../features/reconstruction/types"; + const FONT_FAMILY = "Arial, Helvetica, sans-serif"; const NAME_FONT_SIZE = 12; const LABEL_FONT_SIZE = 12; @@ -49,7 +51,12 @@ function escapeXml(value: string): string { } // Splits a name like "A2^R" or "C^E2" into plain/superscript runs, matching MotifName's convention. +// The link token renders as the same muted glyph MotifName uses instead of its literal text. function splitStereoMarkers(name: string): { text: string; sup: boolean }[] { + if (name === LINK_TOKEN) { + return [{ text: "⋮", sup: false }]; + } + return name .split(/(\^[SREZ])/g) .filter((part) => part !== "") diff --git a/gui/src/client/src/features/discovery/types.ts b/gui/src/client/src/features/discovery/types.ts index 1420801..caff6de 100644 --- a/gui/src/client/src/features/discovery/types.ts +++ b/gui/src/client/src/features/discovery/types.ts @@ -43,12 +43,6 @@ export const SubmitDiscoveryQueryReqSchema = z.object({ includeUserUploads: z.boolean().optional(), onlyUserUploads: z.boolean().optional(), queryOriginSmiles: z.string().nullable().optional(), - // A compound's tailoring-event names (glycosylation, methylation, ...), carried - // over from whichever dbMatchingSequences candidate the query was built from - // (see features/reconstruction/types.ts's Reconstruction.extra_fingerprint_tokens). - // Folded into the query fingerprint server-side without being displayed as part - // of primarySequence -- see run_discovery_query's extra_fingerprint_tokens. - extraFingerprintTokens: z.array(z.string().min(1)).optional(), flags: z.object({ computeMsa: z.boolean().optional(), computeCompare: z.boolean().optional(), diff --git a/gui/src/client/src/features/motifs/api.ts b/gui/src/client/src/features/motifs/api.ts index c3040c4..6ceaa3f 100644 --- a/gui/src/client/src/features/motifs/api.ts +++ b/gui/src/client/src/features/motifs/api.ts @@ -1,10 +1,10 @@ import { getJson } from "../http"; -import { MotifStructuresRespSchema, type MotifStructures } from "./types"; +import { MotifStructuresRespSchema, type MotifStructuresData } from "./types"; -// The whole name -> SMILES vocabulary, fetched once and cached by the caller -// (see MotifHoverCard) -- it's small and effectively static for the life of -// the server process, same rationale as searchMonomerNames' rule list. -export async function fetchMotifStructures(signal?: AbortSignal): Promise { - const data = await getJson("/api/motifStructures", MotifStructuresRespSchema, signal); - return data.structures; +// The whole name -> SMILES vocabulary (plus which names are ambiguous -- see +// MotifStructuresData), fetched once and cached by the caller (see +// MotifHoverCard) -- it's small and effectively static for the life of the +// server process, same rationale as searchMonomerNames' rule list. +export async function fetchMotifStructures(signal?: AbortSignal): Promise { + return getJson("/api/motifStructures", MotifStructuresRespSchema, signal); } diff --git a/gui/src/client/src/features/motifs/types.ts b/gui/src/client/src/features/motifs/types.ts index 1d4d0b2..9530e62 100644 --- a/gui/src/client/src/features/motifs/types.ts +++ b/gui/src/client/src/features/motifs/types.ts @@ -4,6 +4,12 @@ import { z } from "zod"; // for this name" (unidentified "X" blocks, hand-edited names, PK_GROUP_TOKENS). export const MotifStructuresRespSchema = z.object({ structures: z.record(z.string()), + // Names a matching rule doesn't uniquely identify -- e.g. "glycosylation" alone + // names dozens of rules, one per distinct sugar. `structures[name]` for one of + // these is just whichever rule happened to be registered first, not necessarily + // the one that actually matched a given occurrence of that name. + ambiguousNames: z.array(z.string()).default([]), }); export type MotifStructures = z.output["structures"]; +export type MotifStructuresData = z.output; diff --git a/gui/src/client/src/features/reconstruction/types.ts b/gui/src/client/src/features/reconstruction/types.ts index 559b026..c42b8ff 100644 --- a/gui/src/client/src/features/reconstruction/types.ts +++ b/gui/src/client/src/features/reconstruction/types.ts @@ -1,5 +1,35 @@ import { z } from "zod"; +// Joins two candidate primary sequence paths that couldn't be threaded into a +// single path (e.g. a disconnected sugar, or a branched/cyclic assembly) so a +// compound/BGC always has exactly one primary sequence -- see +// retromol.model.readout.merge_named_paths on the backend. Alignable like any +// other token (see gui/src/server/routes/discovery.py's DiscoveryContext), so it +// can appear in an aligned sequence too, not just an unaligned primary sequence. +export const LINK_TOKEN = ""; +export const isLinkToken = (name: string | null | undefined): boolean => name === LINK_TOKEN; + +// Splits a merged primary sequence back apart on LINK_TOKEN, e.g. so the discovery +// query editor can offer "search with just this chain" alongside "search with the +// whole thing" for a compound whose readout has more than one path (a branched/ +// disconnected assembly, or a main chain plus tailoring events like glycosylation). +// Names are [name, tags] pairs (see PrimarySequenceItemSchema) or plain strings -- +// either way the link token itself is dropped, never included in a subsequence. +export function splitOnLinkToken(sequence: T[]): T[][] { + const nameOf = (item: T): string => (Array.isArray(item) ? item[0] : item); + + const groups: T[][] = [[]]; + for (const item of sequence) { + if (isLinkToken(nameOf(item))) { + groups.push([]); + } else { + groups[groups.length - 1].push(item); + } + } + + return groups.filter((g) => g.length > 0); +} + export const PrimarySequenceItemSchema = z.tuple([z.string(), z.array(z.number())]); export type PrimarySequenceItem = z.output; @@ -17,13 +47,6 @@ export const ReconstructionSchema = z.object({ // single path (e.g. a branched or cyclic assembly). Render as an unordered set, // not a sequence. ordered: z.boolean().default(true), - // Names of this compound's tailoring events (glycosylation, methylation, ...) -- - // real identified content that isn't part of any chain, so it's never displayed - // as part of primary_sequence, but should still count toward a discovery query's - // fingerprint (see SubmitDiscoveryQueryReqSchema.extraFingerprintTokens). Always - // empty for `data` (reconstruct_linear_readout candidates don't carry this); - // only ever populated on `dbMatchingSequences` entries. - extra_fingerprint_tokens: z.array(z.string()).default([]), }); export type Reconstruction = z.output; @@ -44,10 +67,15 @@ export const ReconstructCompoundRespSchema = z.object({ ok: z.boolean().optional(), status: z.string().optional(), data: z.array(ReconstructionSchema).default([]), - // Candidates read directly off result.linear_readout.paths -- the same - // representation the database is actually populated with, unlike `data` above - // (reconstruct_linear_readout's candidates, which apply backbone-reconstruction - // eligibility/orientation filtering that can diverge from what's stored). Query - // with one of these, not `data`, to get a fingerprint comparable to the database. + // The single primary sequence read directly off result.linear_readout -- the + // same representation the database is actually populated with, unlike `data` + // above (reconstruct_linear_readout's candidates, which apply backbone- + // reconstruction eligibility/orientation filtering that can diverge from what's + // stored). Always at most one element. Nothing is filtered out of it -- every + // path found in the assembly graph, tailoring events (glycosylation, + // methylation, ...) included, is merged into it, joined by LINK_TOKEN -- so use + // splitOnLinkToken to pull out just one biosynthetic chain if that's what's + // wanted for a query, rather than the whole molecule. Query with this, not + // `data`, to get a fingerprint comparable to the database. dbMatchingSequences: z.array(ReconstructionSchema).default([]), }); diff --git a/gui/src/server/routes/discovery.py b/gui/src/server/routes/discovery.py index 5aff0c1..d4c1703 100644 --- a/gui/src/server/routes/discovery.py +++ b/gui/src/server/routes/discovery.py @@ -25,12 +25,11 @@ from retromol_alignment.aligner import setup_aligner from retromol_alignment.msa import calculate_msa from retromol_alignment.pairwise import Converter, align -from retromol_alignment.ranking import rerank +from retromol_alignment.ranking import reorder_target_chains from retromol_alignment.scoring import HARDCODED_PK_SCORING, create_tanimoto_scoring_matrix from retromol_antismash.modules import LinearReadout as BgcLinearReadout, bgc_primary_sequence from retromol_database.duckdb import FINGERPRINT_SIZE, Entry, SearchResult -from retromol_fingerprint.fingerprint import TOKEN_UNK, Fingerprinter, Vocabulary -from retromol_synthesis.reconstruction import reconstruct_linear_readout +from retromol_fingerprint.fingerprint import TOKEN_LINK, TOKEN_UNK, Fingerprinter, Vocabulary from routes.database import open_retromol_db from routes.queue import JobStillRunningError, enqueue_and_wait, enqueue_job @@ -142,7 +141,12 @@ def _build_context() -> DiscoveryContext: radius=2, num_bits=2048, stereochemistry=False, - self_score_tokens=[TOKEN_UNK, *PK_GROUP_TOKENS], + # TOKEN_LINK gets self_score like TOKEN_UNK/PK_GROUP_TOKENS -- 1.0 similarity + # to itself, 0 to everything else (including gaps, handled separately by the + # aligner's own gap penalties) -- so it aligns like any other token: two + # merge points line up as a match, a merge point against a real residue or a + # gap scores as a mismatch/gap like any other substitution would. + self_score_tokens=[TOKEN_UNK, TOKEN_LINK, *PK_GROUP_TOKENS], self_score=1.0, hardcoded_scores=HARDCODED_PK_SCORING, ) @@ -157,7 +161,7 @@ def _build_context() -> DiscoveryContext: ruleset=rules, name_to_rule=name_to_rule, # PK_GROUP_TOKENS aren't matching-rule names (name_to_rule.get would miss - # them), but they are valid primary-sequence block names -- both + # them), but they are valid primary sequence block names -- both # _per_monomer_tokens and _normalize_for_alignment special-case them -- so # they belong in the autocomplete list the Sequence Editor searches. rule_names_sorted=sorted({*name_to_rule, *PK_GROUP_TOKENS}), @@ -291,55 +295,54 @@ def _cosine_similarity(a: np.ndarray, b: np.ndarray) -> float: def _build_compound_upload_candidate(item: dict, ctx: DiscoveryContext, query_fp: np.ndarray) -> list[SearchResult]: """ - Build synthetic search candidates from one uploaded (and successfully parsed) - compound -- each of its reconstructed primary sequences becomes its own candidate, - using whatever sequence is currently effective for it (the user's edited override - if one was saved, otherwise the algorithm's own parse). + Build a synthetic search candidate from one uploaded (and successfully parsed) + compound, using the same single merged primary sequence -- every path in the + readout, tailoring events included, nothing filtered -- the persistent database + would store for it (see retromol.model.readout.LinearReadout.primary_sequence / + database/scripts/load_compounds.py) -- not reconstruct_linear_readout's + per-path chemistry candidates, which apply backbone-reconstruction eligibility + filtering that can diverge from what's actually stored (see + DB_MATCHING_SEQUENCE_NOTE in routes/jobs.py). This is what makes an uploaded + compound's search candidate directly comparable to a persisted database entry. :param item: the session item (kind "compound", status "done", with a payload) :param ctx: the discovery context - :param query_fp: the query's fingerprint, for scoring candidates against - :return: synthetic SearchResult candidates, one per reconstructed sequence + :param query_fp: the query's fingerprint, for scoring the candidate against + :return: a single-element list holding the synthetic candidate, or empty if + the readout has no primary sequence at all """ try: - reconstructions = reconstruct_linear_readout(Result.from_dict(item["payload"])) + result = Result.from_dict(item["payload"]) + names = result.linear_readout.primary_sequence() except Exception: - current_app.logger.exception("discovery_query: failed to reconstruct upload item_id=%s", item.get("id")) + current_app.logger.exception("discovery_query: failed to read primary sequence for upload item_id=%s", item.get("id")) return [] - overrides = item.get("editedPrimarySequences") or {} - label = item.get("name") or "Uploaded compound" - - candidates: list[SearchResult] = [] - for idx, reconstruction in enumerate(reconstructions): - override = overrides.get(str(idx)) - effective_sequence = override if override is not None else reconstruction.to_dict()["primary_sequence"] - names = [name for name, _tags in effective_sequence] - if not names: - continue + if not names: + return [] - fp = ctx.fingerprinter.encode([_per_monomer_tokens(name, ctx) for name in names]) - candidates.append( - SearchResult( - entry=Entry( - id=f"{UPLOAD_ENTRY_ID_PREFIX}{item['id']}:{idx}", - name=label if len(reconstructions) == 1 else f"{label} #{idx + 1}", - url=None, - # The uploaded compound's own SMILES -- same molecule for every - # reconstruction candidate derived from it. Carried through as `raw` - # so it lines up with how database compound entries store their - # SMILES (see the module docstring / _build_context), letting the - # Tanimoto compare endpoint treat both origins identically. - raw=item.get("smiles"), - type="compound", - primary_sequence=names, - fingerprint=fp.tolist(), - ), - similarity=_cosine_similarity(query_fp, fp), - ) + label = item.get("name") or "Uploaded compound" + tokens = [_per_monomer_tokens(name, ctx) for name in names if name != TOKEN_LINK] + fp = ctx.fingerprinter.encode(tokens) + + return [ + SearchResult( + entry=Entry( + id=f"{UPLOAD_ENTRY_ID_PREFIX}{item['id']}:0", + name=label, + url=None, + # The uploaded compound's own SMILES, so it lines up with how database + # compound entries store their SMILES (see the module docstring / + # _build_context), letting the Tanimoto compare endpoint treat both + # origins identically. + raw=item.get("smiles"), + type="compound", + primary_sequence=names, + fingerprint=fp.tolist(), + ), + similarity=_cosine_similarity(query_fp, fp), ) - - return candidates + ] def _build_bgc_upload_candidate(item: dict, ctx: DiscoveryContext, query_fp: np.ndarray) -> list[SearchResult]: @@ -435,12 +438,21 @@ def _build_upload_candidates( def _per_monomer_tokens(name: str, ctx: DiscoveryContext) -> list[str]: """ - Build the fingerprinting token list for one primary-sequence block. + Build the fingerprinting token list for one primary sequence block. + + TOKEN_LINK is not a building block -- it just joins two merged paths -- and + callers filter it out of `names` before mapping over this function (an empty + token list would otherwise fall back to TOKEN_UNK and silently add unknown-token + mass for something that isn't unknown, it's just not a block at all). Handled + explicitly anyway so a stray call can't hit that fallback. :param name: the display-layer block name :param ctx: the discovery context :return: tokens for Fingerprinter.encode (empty list -> falls back to the unknown token) """ + if name == TOKEN_LINK: + return [TOKEN_LINK] + if name in PK_GROUP_TOKENS: # Not a matching-rule name -- the group-level pseudonym a BGC's PKS module # resolves to (see PK_GROUP_TOKENS). Still resolvable to fingerprint tokens @@ -459,7 +471,7 @@ def _per_monomer_tokens(name: str, ctx: DiscoveryContext) -> list[str]: @blp_discovery_monomer_names.get("/api/discoveryMonomerNames") def discovery_monomer_names() -> tuple[Response, int]: """ - Autocomplete endpoint for valid primary-sequence block names -- matching rule + Autocomplete endpoint for valid primary sequence block names -- matching rule names plus the PK_GROUP_TOKENS pseudo-names (e.g. "PK_A") used for a reduction-level-only PKS call. @@ -507,16 +519,32 @@ def motif_structures() -> tuple[Response, int]: one is set. Names with no matching rule (unidentified "X" blocks, hand-edited names, PK_GROUP_TOKENS) are simply absent from the map. + A rule *name* isn't unique to one structure -- e.g. "glycosylation" alone + names 26 different rules, one per distinct sugar, since the primary-sequence + display convention only ever shows the tailoring-event name, not which + specific rule matched. `structures` picks one (via `ctx.name_to_rule`, + whichever rule was first registered for that name) rather than showing + nothing, but that pick is arbitrary -- `ambiguousNames` lists every name this + is true for, so the frontend can warn that the depiction shown isn't + necessarily the one that actually matched a given occurrence of that name. + The whole vocabulary is returned in one shot rather than per-name, since it's small (on the order of a few hundred names) and effectively static for the life of the process -- same rationale as rule_names_sorted powering the autocomplete above. - :return: a tuple containing a dictionary with the name -> SMILES map and an HTTP status code + :return: a tuple containing a dictionary with the name -> SMILES map, the + list of ambiguous names, and an HTTP status code """ ctx = get_discovery_context() structures = {name: (rule.display_smiles or rule.smiles) for name, rule in ctx.name_to_rule.items()} - return jsonify({"structures": structures}), 200 + + smiles_by_name: dict[str, set[str]] = {} + for rule in ctx.ruleset.matching_rules: + smiles_by_name.setdefault(rule.name, set()).add(rule.display_smiles or rule.smiles) + ambiguous_names = sorted(name for name, smiles in smiles_by_name.items() if len(smiles) > 1) + + return jsonify({"structures": structures, "ambiguousNames": ambiguous_names}), 200 def run_discovery_query( @@ -528,7 +556,6 @@ def run_discovery_query( only_user_uploads: bool, include_user_uploads: bool, session_id: str | None, - extra_fingerprint_tokens: list[str] | None = None, ) -> tuple[dict, int]: """ Fingerprint a primary sequence, retrieve nearest neighbors from the database, and @@ -538,22 +565,24 @@ def run_discovery_query( lazily via get_discovery_context() the first time it runs in a given worker process, same as any other caller of that function. - :param extra_fingerprint_tokens: names of a compound's tailoring events - (glycosylation, methylation, ...) that don't belong in `primary_sequence` - itself (they're not part of any chain) but should still count toward the - query fingerprint -- mirrors how database/scripts/load_compounds.py folds - the same thing into a stored compound entry's fingerprint. Never touches - alignment_query below: alignment compares against `primary_sequence` only. + `primary_sequence` is expected to already be whatever the caller wants + fingerprinted and aligned -- the full merged sequence (tailoring events and + all, joined by TOKEN_LINK) or a subsequence split off at TOKEN_LINK, e.g. to + search with just one biosynthetic chain rather than the whole molecule. There's + no separate "extra tokens" side channel any more: everything that should count + toward the query lives in `primary_sequence` itself, the same single + representation database/scripts/load_compounds.py stores. + :return: a (response body, HTTP status code) pair """ ctx = get_discovery_context() # Fingerprint is built from the display sequence directly (Fingerprinter.encode # already treats an empty token list as "unknown", matching how unidentified - # blocks were encoded when the database was originally built), plus any - # tailoring-event tokens folded in on top -- see extra_fingerprint_tokens above. - per_monomer_tokens = [_per_monomer_tokens(name, ctx) for name in primary_sequence] - per_monomer_tokens += [_per_monomer_tokens(name, ctx) for name in (extra_fingerprint_tokens or [])] + # blocks were encoded when the database was originally built). TOKEN_LINK is + # excluded -- it just joins merged paths, it isn't a block (see + # database/scripts/load_compounds.py, which excludes it the same way). + per_monomer_tokens = [_per_monomer_tokens(name, ctx) for name in primary_sequence if name != TOKEN_LINK] query_fp = ctx.fingerprinter.encode(per_monomer_tokens) # Alignment needs every token to exist in the scoring-matrix alphabet, so @@ -616,7 +645,17 @@ def run_discovery_query( else: skipped += 1 - reranked = rerank(alignment_query, usable_targets, ctx.aligner, ctx.converter) if usable_targets else [] + # Reorders (and, per chain, reorients) each target's chains to best match the + # query's own chain order -- an order-agnostic optimal assignment over every + # (query_chain, target_chain) pair, since which chain landed in which position + # after merge_named_paths' longest-first/lexicographic sort has nothing to do + # with which chains actually correspond biosynthetically between two different + # compounds. See retromol_alignment.ranking.reorder_target_chains. + reordered = ( + [reorder_target_chains(alignment_query, target, ctx.aligner, ctx.converter) for target in usable_targets] + if usable_targets + else [] + ) # Only "longest_sequence" needs each candidate's own self-alignment score -- # "subsequence" normalizes every candidate by the same constant (self_score), @@ -632,14 +671,12 @@ def rank_key(align_score: float, target_self_score: float | None) -> float: return _normalized_pct(align_score, max(self_score, target_self_score)) return align_score - scored = list(zip(usable_candidates, usable_targets, reranked, target_self_scores)) - scored.sort(key=lambda item: rank_key(item[2][0], item[3]), reverse=True) + scored = list(zip(usable_candidates, reordered, target_self_scores)) + scored.sort(key=lambda item: rank_key(item[1][0], item[2]), reverse=True) top = scored[:top_x] results: list[dict[str, Any]] = [] - for candidate, target, (score, inverted), target_self_score in top: - oriented_target = target[::-1] if inverted else target - + for candidate, (_assignment_score, oriented_target, inverted), target_self_score in top: try: align_score, aligned_query_display, aligned_target_display = align( ctx.aligner, alignment_query, oriented_target, ctx.converter @@ -650,8 +687,10 @@ def rank_key(align_score: float, target_self_score: float | None) -> float: ) continue - # Self-alignment score is invariant under reversing both operands, so the - # score computed against the un-reversed target is still valid here. + # A self-alignment always achieves a perfect 1:1, gap-free match (score = + # sum of each token's own self-similarity), which doesn't depend on token + # order or per-chain orientation -- so target_self_score, computed on the + # original (pre-reorder) target, is still exactly right for oriented_target. denom = max(self_score, target_self_score) if target_self_score is not None else self_score normalized_pct = _normalized_pct(align_score, denom) @@ -728,7 +767,6 @@ def discovery_query() -> tuple[Response, int]: only_user_uploads = bool(payload.get("onlyUserUploads", False)) include_user_uploads = bool(payload.get("includeUserUploads", False)) or only_user_uploads session_id = payload.get("sessionId") - extra_fingerprint_tokens = payload.get("extraFingerprintTokens") if ( not isinstance(primary_sequence, list) @@ -737,12 +775,6 @@ def discovery_query() -> tuple[Response, int]: ): return jsonify({"error": "primarySequence must be a non-empty list of non-empty strings"}), 400 - if extra_fingerprint_tokens is not None and ( - not isinstance(extra_fingerprint_tokens, list) - or not all(isinstance(x, str) and x for x in extra_fingerprint_tokens) - ): - return jsonify({"error": "extraFingerprintTokens must be a list of non-empty strings, if given"}), 400 - if entry_type not in ENTRY_TYPES_FOR_QUERY: return jsonify({"error": f"entryType must be one of {ENTRY_TYPES_FOR_QUERY}"}), 400 @@ -763,7 +795,6 @@ def discovery_query() -> tuple[Response, int]: body, status = enqueue_and_wait( run_discovery_query, primary_sequence, entry_type, n, top_x, score_mode, only_user_uploads, include_user_uploads, session_id, - extra_fingerprint_tokens, ) except JobStillRunningError as e: return jsonify({"error": str(e)}), 503 @@ -1045,7 +1076,6 @@ def mark_processing(it: dict) -> None: settings.get("onlyUserUploads", False), settings.get("includeUserUploads", False), session_id, - settings.get("extraFingerprintTokens"), ) if query_status != 200: raise RuntimeError(query_body.get("error", "Discovery query failed")) @@ -1138,10 +1168,6 @@ def submit_discovery_query() -> tuple[Response, int]: query_origin_smiles = payload.get("queryOriginSmiles") flags = payload.get("flags") or {} name = payload.get("name") or "Discovery query" - # A compound's tailoring events (glycosylation, methylation, ...) -- see - # run_discovery_query's docstring. Optional: absent/empty means "none", not - # an error, since most queries (and every BGC query) won't have any. - extra_fingerprint_tokens = payload.get("extraFingerprintTokens") if not isinstance(session_id, str) or not session_id: return jsonify({"error": "Missing sessionId"}), 400 @@ -1153,12 +1179,6 @@ def submit_discovery_query() -> tuple[Response, int]: ): return jsonify({"error": "primarySequence must be a non-empty list of non-empty strings"}), 400 - if extra_fingerprint_tokens is not None and ( - not isinstance(extra_fingerprint_tokens, list) - or not all(isinstance(x, str) and x for x in extra_fingerprint_tokens) - ): - return jsonify({"error": "extraFingerprintTokens must be a list of non-empty strings, if given"}), 400 - if entry_type not in ENTRY_TYPES_FOR_QUERY: return jsonify({"error": f"entryType must be one of {ENTRY_TYPES_FOR_QUERY}"}), 400 @@ -1202,7 +1222,6 @@ def submit_discovery_query() -> tuple[Response, int]: "includeUserUploads": include_user_uploads, "onlyUserUploads": only_user_uploads, "queryOriginSmiles": query_origin_smiles, - "extraFingerprintTokens": extra_fingerprint_tokens, } item = { diff --git a/gui/src/server/routes/jobs.py b/gui/src/server/routes/jobs.py index 0801f35..69a3c54 100644 --- a/gui/src/server/routes/jobs.py +++ b/gui/src/server/routes/jobs.py @@ -17,8 +17,10 @@ from retromol.chem.tagging import get_tags_mol from retromol.model.submission import Submission from retromol.model.rules import RuleSet +from retromol.model.readout import merge_named_paths from retromol.model.result import Result from retromol.pipelines.parsing import run_retromol +from retromol_fingerprint.fingerprint import TOKEN_LINK from retromol_synthesis.reconstruction import reconstruct_linear_readout from retromol_antismash.io import parse_antismash_gbk, AntiSmashOptions @@ -267,70 +269,62 @@ def submit_compound() -> tuple[Response, int]: DB_MATCHING_SEQUENCE_NOTE = ( "This is the exact primary sequence the database-population pipeline stores for " - "this molecule (read directly off result.linear_readout.paths -- no backbone " - "reconstruction, no eligibility/orientation filtering). Use one of these to " - "search, not a reconstructed candidate above: reconstruct_linear_readout builds " - "a separate identified-only readout to attempt backbone reconstruction, which " - "diverges from what's actually stored for a real fraction of compounds " - "(anything with unidentified or tailoring-only content) -- querying with that " - "instead silently caps how similar a match can ever score, even against the " - "molecule's own database entry." + "this molecule (read directly off result.linear_readout with no backbone " + "reconstruction, no eligibility/orientation filtering, no length filtering: " + "every path found in the assembly graph is included, tailoring events like " + "glycosylation/methylation among them). Every path is merged into this one " + "sequence, longest first, ties broken lexicographically, joined by \"\". " + "The same merge database/scripts/common.py's primary_sequence_from_result " + "applies before storing a compound. Use this to search, not a reconstructed " + "candidate above: reconstruct_linear_readout builds a separate identified-only " + "readout to attempt backbone reconstruction, which diverges from what's actually " + "stored for a real fraction of compounds (anything with unidentified content). " + "Querying with that instead silently caps how similar a match can ever score, " + "even against the molecule's own database entry. To search with just one " + "biosynthetic chain rather than the whole molecule, split this sequence on " + "\"\" and use one of the resulting subsequences instead." ) -def _db_matching_primary_sequences(result: Result, min_length: int = 2) -> list[dict]: +def _db_matching_primary_sequences(result: Result) -> list[dict]: """ - Every candidate primary sequence read directly off result.linear_readout.paths -- - the same representation database/scripts/common.py's primary_sequences_from_result - uses to populate the persistent database (itself reproducing the original + The single primary sequence read directly off result.linear_readout -- the same + representation database/scripts/common.py's primary_sequence_from_result uses to + populate the persistent database (itself reproducing the original db_scripts/create_database.py recipe, recovered from git history at commit - 72a1bbb). Shaped identically to a Reconstruction dict so the frontend can render - and pick from these the same way it does reconstruct_linear_readout's candidates, - just without ever calling that function -- see DB_MATCHING_SEQUENCE_NOTE for why - the two aren't interchangeable as a database query. - - Paths shorter than min_length -- almost always a lone tailoring event - (glycosylation, methylation: AssemblyGraph only keeps C-C/C-N bonds as - "connections", so a sugar attached via a glycosidic C-O-C linkage always ends - up disconnected from the main chain) -- never become their own candidate, same - as database/scripts/common.py. But they're real identified content, so their - names are collected into every candidate's extra_fingerprint_tokens instead: - the frontend carries this alongside primary_sequence (see - features/reconstruction/types.ts) and the query submission flow folds it into - the query fingerprint without displaying it as part of the sequence (see - run_discovery_query's extra_fingerprint_tokens parameter). Not deduplicated, - same reasoning as the database side: two glycosylation events should count for - roughly twice the weight of one. + 72a1bbb). Every path is merged into one sequence, nothing filtered out (see + retromol.model.readout.LinearReadout.primary_sequence / merge_named_paths). + Shaped identically to a Reconstruction dict, in a one-element list, so the + frontend can render it the same way it does reconstruct_linear_readout's + candidates, just without ever calling that function -- see + DB_MATCHING_SEQUENCE_NOTE for why the two aren't interchangeable as a database + query. :param result: a parsed RetroMol Result - :param min_length: drop paths shorter than this from the candidate list itself - :return: one Reconstruction-shaped dict per path meeting min_length + :return: a one-element list holding the single merged Reconstruction-shaped + dict, or empty if the readout has no paths at all """ tagged_input_smiles = mol_to_smiles(result.submission.mol, include_tags=True) - sequences = [] - extra_fingerprint_tokens: list[str] = [] - for path in result.linear_readout.paths: - names = [node.identity.matched_rule.name if node.is_identified else "X" for node in path] - if len(path) < min_length: - extra_fingerprint_tokens.extend(n for n in names if n != "X") - continue - - primary_sequence = [[name, list(get_tags_mol(node.mol))] for name, node in zip(names, path)] - sequences.append({ - "tagged_input_smiles": tagged_input_smiles, - "tagged_backbone_smiles": None, - "primary_sequence": primary_sequence, - "backbone_warning": DB_MATCHING_SEQUENCE_NOTE, - "ordered": True, - }) - - # Every candidate shares the same extra_fingerprint_tokens -- tailoring events - # belong to the whole molecule, not to any one path through it. - for sequence in sequences: - sequence["extra_fingerprint_tokens"] = extra_fingerprint_tokens - - return sequences + named_paths = [ + [[name, list(get_tags_mol(node.mol))] for name, node in zip( + [n.identity.matched_rule.name if n.is_identified else "X" for n in path], path + )] + for path in result.linear_readout.paths + ] + + if not named_paths: + return [] + + primary_sequence = merge_named_paths(named_paths, key=lambda item: item[0], link_item=[TOKEN_LINK, []]) + + return [{ + "tagged_input_smiles": tagged_input_smiles, + "tagged_backbone_smiles": None, + "primary_sequence": primary_sequence, + "backbone_warning": DB_MATCHING_SEQUENCE_NOTE, + "ordered": True, + }] def run_compound_reconstruction(item_payload: dict | None) -> tuple[dict, int]: diff --git a/src/retromol/model/readout.py b/src/retromol/model/readout.py index d2e58e2..ac5ed46 100644 --- a/src/retromol/model/readout.py +++ b/src/retromol/model/readout.py @@ -1,17 +1,72 @@ """Data structures for representing readouts from RetroMol parsing results.""" from dataclasses import dataclass -from typing import Literal +from typing import Callable, Literal, TypeVar from retromol.model.reaction_graph import MolNode, ReactionGraph from retromol.model.assembly_graph import AssemblyGraph from retromol.model.rules import MatchingRule from retromol.chem.mol import encode_mol from retromol.chem.tagging import get_tags_mol +from retromol_fingerprint.fingerprint import TOKEN_LINK ReadoutMode = Literal["leaf_identified", "first_identified"] +T = TypeVar("T") + + +def merge_named_paths(paths: list[list[T]], key: Callable[[T], str], link_item: T) -> list[T]: + """ + Merge candidate primary sequence paths into a single sequence, joined by + `link_item`, so a compound/BGC with more than one disconnected path (e.g. a + branched/cyclic assembly, or a sugar attached via a glycosidic bond that + AssemblyGraph never treats as a "connection") still gets exactly one primary + sequence to store/query/align, instead of one entry per path. + + Longer paths sort first; paths of equal length are ordered lexicographically + by their own token sequence (via `key`), so the merge is deterministic and + reproducible regardless of the order `paths` was produced in. + + :param paths: candidate paths, each a list of items (e.g. names, or (name, tags) pairs) + :param key: extracts the string used for sorting/tie-breaking from one item + :param link_item: the item inserted between two merged paths (e.g. TOKEN_LINK, + or a (TOKEN_LINK, []) pair matching the shape of the other items) + :return: one flat list: every path's items in sorted order, `link_item` between them + """ + ordered = sorted(paths, key=lambda items: (-len(items), [key(item) for item in items])) + + merged: list[T] = [] + for i, items in enumerate(ordered): + if i > 0: + merged.append(link_item) + merged.extend(items) + + return merged + + +def split_named_paths(sequence: list[T], key: Callable[[T], str], link_token: str = TOKEN_LINK) -> list[list[T]]: + """ + Inverse of `merge_named_paths`: split a merged sequence back apart into its + constituent chains, wherever an item's `key` equals `link_token`. The link + items themselves are dropped, never included in a returned chain -- so + `merge_named_paths(split_named_paths(seq, key, tok), key, link_item)` round-trips + `seq` (given a `link_item` whose `key` is `tok`). + + :param sequence: a merged sequence, as returned by `merge_named_paths` + :param key: extracts the string used to test for the link token from one item + :param link_token: the token marking a join point (default: TOKEN_LINK) + :return: the chains, in their original relative order, empty chains dropped + """ + chains: list[list[T]] = [[]] + for item in sequence: + if key(item) == link_token: + chains.append([]) + else: + chains[-1].append(item) + + return [chain for chain in chains if chain] + @dataclass(frozen=True) class LinearReadout: @@ -74,7 +129,31 @@ def from_reaction_graph( paths.append(path) return cls(assembly_graph=a, paths=paths) - + + def primary_sequence(self) -> list[str]: + """ + The single primary sequence for this readout: every path in `self.paths` + -- including single-node ones, e.g. a lone tailoring event like + glycosylation or methylation that AssemblyGraph never connects to the + main chain (it only keeps C-C/C-N bonds as "connections", so e.g. a sugar + attached via a glycosidic C-O-C linkage always ends up disconnected) -- + merged into one sequence via `merge_named_paths` (longest first, ties + broken lexicographically, joined by TOKEN_LINK). Nothing found in the + assembly graph is dropped; a caller that wants just one biosynthetic + chain, without tailoring events mixed in, can split the result back + apart on TOKEN_LINK. + + An unidentified node is named "X", the convention used everywhere else in + RetroMol. + + :return: the single merged primary sequence + """ + named_paths = [ + [node.identity.matched_rule.name if node.is_identified else "X" for node in path] + for path in self.paths + ] + return merge_named_paths(named_paths, key=lambda n: n, link_item=TOKEN_LINK) + def to_dict(self) -> dict: """ Serialize the LinearReadout to a dictionary. diff --git a/src/retromol_alignment/ranking.py b/src/retromol_alignment/ranking.py index a16cf67..991cc9d 100644 --- a/src/retromol_alignment/ranking.py +++ b/src/retromol_alignment/ranking.py @@ -1,24 +1,114 @@ +import numpy as np +from scipy.optimize import linear_sum_assignment + +from retromol.model.readout import split_named_paths +from retromol_fingerprint.fingerprint import TOKEN_LINK from retromol_alignment.aligner import PairwiseAligner from retromol_alignment.pairwise import Converter, align -def rerank( +def reorder_target_chains( query: list[str], - targets: list[list[str]], + target: list[str], aligner: PairwiseAligner, converter: Converter, -) -> list[tuple[float, bool]]: - scores: list[tuple[float, bool]] = [] + link_token: str = TOKEN_LINK, +) -> tuple[float, list[str], bool]: + """ + Reorder and reorient `target`'s chains (split on `link_token`) to best match + `query`'s chains, order-agnostically -- a merged primary sequence's chain + order is an artifact of how `merge_named_paths` sorted its candidate paths + (longest first, then lexicographically), not a biosynthetic correspondence + between two different compounds' chains, so a straight flat alignment (or + the old whole-sequence-only `rerank`) can easily misalign chain N of one + sequence against a structurally unrelated chain M of the other just because + they happened to land at the same position. + + `query`'s own chain order is left untouched -- it's the fixed reference every + target gets reordered against, the same role it already played in the old + whole-sequence rerank (query fixed, target flipped). + + Matching is a 1:1 optimal assignment (Hungarian algorithm) over every + (query_chain, target_chain) pair's own best-of-forward/reverse alignment + score -- exact rather than greedy, which matters here since a greedy pick + of the single best-scoring pair can lock out a better overall assignment + (e.g. two chains that are each other's second-best match, but each other's + globally best match is a chain the other also wants). Chain counts are + small (typically 2-5), so this optimal solve is cheap. + + When the two sides have different chain counts (e.g. a tailoring event + present on one side but not the other), the leftover chains on the larger + side are appended after the matched ones, in their original relative order + -- unmatched on purpose, so they end up aligning against a gap run in the + final flat alignment (penalizing the missing content) instead of silently + disappearing. + + :param query: alignment-normalized query token sequence (may contain link_token) + :param target: alignment-normalized target token sequence (may contain link_token) + :param aligner: the PairwiseAligner to score/align chain pairs with + :param converter: Converter for the aligner's alphabet + :param link_token: the token marking a join point between chains (default: TOKEN_LINK) + :return: (assignment_score, reordered_target, any_chain_reversed) -- + assignment_score is the sum of each matched pair's own alignment score, + informational/for ranking only: the score actually shown/used for a + candidate should come from re-aligning `query` against `reordered_target` + as one flat sequence (see routes/discovery.py's run_discovery_query), so + chain-boundary gaps are scored consistently with every other candidate. + any_chain_reversed is True if any matched chain needed flipping. + """ + identity = lambda item: item # noqa: E731 -- items here are already the plain string tokens + + query_chains = split_named_paths(query, identity, link_token) + target_chains = split_named_paths(target, identity, link_token) + + if not query_chains or not target_chains: + return 0.0, list(target), False + + n_q, n_t = len(query_chains), len(target_chains) + scores = np.zeros((n_q, n_t)) + reversed_flags = np.zeros((n_q, n_t), dtype=bool) + + for i, q_chain in enumerate(query_chains): + for j, t_chain in enumerate(target_chains): + score_fwd, _, _ = align(aligner, q_chain, t_chain, converter) + score_rev, _, _ = align(aligner, q_chain, t_chain[::-1], converter) + if score_rev > score_fwd: + scores[i, j] = score_rev + reversed_flags[i, j] = True + else: + scores[i, j] = score_fwd + + # linear_sum_assignment minimizes cost; negate to maximize alignment score. + # It handles a non-square matrix natively, matching min(n_q, n_t) pairs and + # leaving the rest of the larger side unmatched -- exactly the "leftover + # chains" case described above. + row_ind, col_ind = linear_sum_assignment(-scores) + + # zip(row_ind, col_ind) is already in ascending query-chain order (that's + # how linear_sum_assignment returns it), so the matched chains come out in + # the query's own chain order -- the natural choice, since query is the + # fixed reference. + reordered: list[str] = [] + assignment_score = 0.0 + any_reversed = False - for target in targets: - score1, _, _ = align(aligner, query, target, converter) - score2, _, _ = align(aligner, query, target[::-1], converter) - score = max(score1, score2) + for i, j in zip(row_ind, col_ind): + chain = target_chains[j] + if reversed_flags[i, j]: + chain = chain[::-1] + any_reversed = True - inverted = False - if score2 > score1: - inverted = True + if reordered: + reordered.append(link_token) + reordered.extend(chain) + assignment_score += scores[i, j] - scores.append((score, inverted)) + matched_target_idx = set(col_ind.tolist()) + for j in range(n_t): + if j in matched_target_idx: + continue + if reordered: + reordered.append(link_token) + reordered.extend(target_chains[j]) - return scores \ No newline at end of file + return assignment_score, reordered, any_reversed diff --git a/src/retromol_antismash/modules.py b/src/retromol_antismash/modules.py index 3bf8beb..690a261 100644 --- a/src/retromol_antismash/modules.py +++ b/src/retromol_antismash/modules.py @@ -1099,7 +1099,7 @@ def module_primary_sequence_tokens(module: Module, ruleset: RuleSet) -> tuple[st def bgc_primary_sequence(readout: LinearReadout, ruleset: RuleSet) -> tuple[list[str], list[list[str]]]: """ - Convert a BGC's linear readout into a primary-sequence representation comparable + Convert a BGC's linear readout into a primary sequence representation comparable to a compound's: one display name and one fingerprint token list per module, both drawn from the same matching-rule vocabulary, in biosynthetic order. diff --git a/src/retromol_database/duckdb.py b/src/retromol_database/duckdb.py index 84e4ad9..0e76812 100644 --- a/src/retromol_database/duckdb.py +++ b/src/retromol_database/duckdb.py @@ -7,7 +7,7 @@ import duckdb import numpy as np -from retromol_fingerprint.fingerprint import TOKEN_UNK +from retromol_fingerprint.fingerprint import TOKEN_LINK, TOKEN_UNK ENTRY_TYPES = ("compound", "bgc") EntryType = Literal["compound", "bgc"] @@ -228,8 +228,9 @@ def stats(self) -> DatabaseStats: Building blocks are the tokens in each entry's primary_sequence -- e.g. amino acid names, PK reduction-state groups, or tailoring events like "methylation". - TOKEN_UNK ("") marks a block RetroMol couldn't identify and is excluded - from unique_block_count since it isn't a real building block. + TOKEN_UNK ("") marks a block RetroMol couldn't identify and TOKEN_LINK + ("") just joins two merged paths within one entry's sequence -- neither + is a real building block, so both are excluded from unique_block_count. :return: a DatabaseStats snapshot """ @@ -258,9 +259,9 @@ def stats(self) -> DatabaseStats: """ SELECT count(DISTINCT token) FROM (SELECT unnest(primary_sequence) AS token FROM entries) - WHERE token != ? + WHERE token NOT IN (?, ?) """, - [TOKEN_UNK], + [TOKEN_UNK, TOKEN_LINK], ).fetchone()[0] ) diff --git a/src/retromol_fingerprint/fingerprint.py b/src/retromol_fingerprint/fingerprint.py index 37bfac0..69eeb6b 100644 --- a/src/retromol_fingerprint/fingerprint.py +++ b/src/retromol_fingerprint/fingerprint.py @@ -4,7 +4,12 @@ TOKEN_UNK = "" -SPECIAL_TOKENS = [TOKEN_UNK] +# Joins two candidate primary sequence paths that couldn't be threaded into one +# path (e.g. a disconnected sugar, or a branched/cyclic assembly) into a single +# sequence, so every compound/BGC always has exactly one primary sequence. See +# retromol.model.readout.merge_named_paths. +TOKEN_LINK = "" +SPECIAL_TOKENS = [TOKEN_UNK, TOKEN_LINK] class Vocabulary: