Add offline WURCS-to-SMILES glycan translator and wire it into both pipelines

nextflow/bin/get_ec_information.py

@@ -21,7 +21,7 @@
 import argparse
 from pathlib import Path
 from rdkit.ML.Descriptors.MoleculeDescriptors import MolecularDescriptorCalculator
-from utils import get_terminal_record, get_csdb_from_glycoct, get_smiles_from_csdb, process_ec_records
+from utils import get_terminal_record, get_csdb_from_glycoct, get_smiles_from_csdb, get_smiles_from_wurcs_offline, process_ec_records
 from bs4 import BeautifulSoup
 from urllib.parse import quote
 
@@ -646,6 +646,17 @@ def unpack_sets(row):
         new_smiles_values.rename(columns = {"smiles":"descriptor"}, inplace = True)
         smiles_cache = pd.concat([smiles_cache, new_smiles_values], ignore_index = True)
         smiles_cache.to_pickle(f"{args.smiles_cache}")
+
+        #offline fallback for rows the live glycoct/CSDB chain above didn't resolve.
+        #GlyTouCan's gtcid2seqs API no longer returns glycoct at all (confirmed dead,
+        #0% success on real data - see docs/iupac_translator_plan.md), so in practice
+        #this is currently the path that resolves nearly everything that resolves at
+        #all - benchmarked at ~89% on the real cognate-ligand master set. Kept as a
+        #fallback rather than replacing the live chain outright, in case GlyTouCan's
+        #API is fixed in future.
+        missing_smiles_mask = glycan_compounds_df_merged["smiles"].isna() & glycan_compounds_df_merged["wurcs"].notna()
+        glycan_compounds_df_merged.loc[missing_smiles_mask, "smiles"] = glycan_compounds_df_merged.loc[missing_smiles_mask, "wurcs"].apply(get_smiles_from_wurcs_offline)
+
         glycan_compounds_df_merged = glycan_compounds_df_merged.loc[glycan_compounds_df_merged.smiles.isna() == False]
         kegg_reaction_enzyme_df_exploded_gtc = kegg_reaction_enzyme_df_exploded.merge(glycan_compounds_df_merged, left_on="entities", right_on="compound_id", how = "inner")
         PandasTools.AddMoleculeColumnToFrame(kegg_reaction_enzyme_df_exploded_gtc, smilesCol='smiles')

nextflow/bin/process_all_pdb_contacts.py

@@ -3,7 +3,7 @@
 import argparse
 import pandas as pd
 from gemmi import cif
-from utils import process_ec_records, get_updated_enzyme_records, get_scop_domains_info, extract_interpro_domain_annotations, get_pfam_annotations, get_glycoct_from_wurcs, get_csdb_from_glycoct, get_smiles_from_csdb, build_cath_dataframe, parse_cddf, build_g3dsa_dataframe, get_scop2_domains_info
+from utils import process_ec_records, get_updated_enzyme_records, get_scop_domains_info, extract_interpro_domain_annotations, get_pfam_annotations, get_glycoct_from_wurcs, get_csdb_from_glycoct, get_smiles_from_csdb, get_smiles_from_wurcs_offline, build_cath_dataframe, parse_cddf, build_g3dsa_dataframe, get_scop2_domains_info
 import numpy as np
 from Bio.ExPASy import Enzyme as EEnzyme
 import re
@@ -16,15 +16,20 @@ def get_sugar_smiles_from_wurcs(wurcs_list, csdb_linear_cache, smiles_cache, gly
     updated_csdb_cache = []
     updated_smiles_cache = []
     for wurcs in wurcs_list:
-        smiles = None
-        glycoct = get_glycoct_from_wurcs(wurcs, glycoct_cache)
-        updated_glycoct_cache.append({"WURCS": wurcs, "glycoct": glycoct})
-        if not pd.isna(glycoct):
-            csdb = get_csdb_from_glycoct(glycoct, csdb_linear_cache)
-            updated_csdb_cache.append({"glycoct": glycoct, "csdb": csdb})
-            if not pd.isna(csdb):
-                smiles = get_smiles_from_csdb(csdb, smiles_cache)
-                updated_smiles_cache.append({"csdb": csdb, "descriptor" : smiles})
+        # offline translation (glypy -> IUPAC-condensed -> glyles) is the
+        # primary path here - benchmarked at ~96% on real PDB-deposited
+        # glycans, see docs/iupac_translator_plan.md. Only fall back to
+        # the live GlycoSmos/CSDB chain (below) when it fails.
+        smiles = get_smiles_from_wurcs_offline(wurcs)
+        if pd.isna(smiles):
+            glycoct = get_glycoct_from_wurcs(wurcs, glycoct_cache)
+            updated_glycoct_cache.append({"WURCS": wurcs, "glycoct": glycoct})
+            if not pd.isna(glycoct):
+                csdb = get_csdb_from_glycoct(glycoct, csdb_linear_cache)
+                updated_csdb_cache.append({"glycoct": glycoct, "csdb": csdb})
+                if not pd.isna(csdb):
+                    smiles = get_smiles_from_csdb(csdb, smiles_cache)
+                    updated_smiles_cache.append({"csdb": csdb, "descriptor" : smiles})
         sugar_smiles[wurcs] = smiles
     updated_glycoct_cache_df = pd.concat([pd.DataFrame(updated_glycoct_cache, columns = ["WURCS", "glycoct"]), glycoct_cache]).drop_duplicates()
     updated_csdb_cache_df = pd.concat([pd.DataFrame(updated_csdb_cache, columns = ["glycoct", "csdb"]), csdb_linear_cache]).drop_duplicates()

nextflow/bin/tests/test_get_ec_information_glycan_fallback.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+"""
+get_ec_information.py's Context A glycan-resolution block is inline
+script code (not a function), so it can't be called directly in a test.
+This instead replicates its exact offline-fallback logic
+(missing_smiles_mask / .apply(get_smiles_from_wurcs_offline)) against a
+small synthetic dataframe, to validate that pattern in isolation without
+needing the full pipeline's upstream KEGG/Rhea/GlyTouCan machinery.
+
+    python3 nextflow/bin/tests/test_get_ec_information_glycan_fallback.py
+"""
+
+import sys
+import unittest
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+
+from utils import get_smiles_from_wurcs_offline
+
+CHITOBIOSE_WURCS = "WURCS=2.0/1,2,1/[a2122h-1b_1-5_2*NCC/3=O]/1-1/a4-b1"
+
+
+class TestOfflineFallbackMasking(unittest.TestCase):
+
+    def test_only_missing_rows_with_a_wurcs_value_are_backfilled(self):
+        df = pd.DataFrame({
+            "compound_id": ["G1", "G2", "G3", "G4"],
+            # G1: live chain already resolved it -> must not be overwritten
+            # G2: live chain failed, has a wurcs -> should be backfilled offline
+            # G3: live chain failed, no wurcs at all -> stays nan
+            # G4: live chain failed, wurcs is unparseable -> stays nan
+            "smiles": ["C(already resolved)", np.nan, np.nan, np.nan],
+            "wurcs": [CHITOBIOSE_WURCS, CHITOBIOSE_WURCS, np.nan, "not a wurcs string"],
+        })
+
+        # exact logic from get_ec_information.py
+        missing_smiles_mask = df["smiles"].isna() & df["wurcs"].notna()
+        df.loc[missing_smiles_mask, "smiles"] = df.loc[missing_smiles_mask, "wurcs"].apply(get_smiles_from_wurcs_offline)
+
+        self.assertEqual(df.loc[0, "smiles"], "C(already resolved)")
+        self.assertIsInstance(df.loc[1, "smiles"], str)
+        self.assertNotEqual(df.loc[1, "smiles"], "")
+        self.assertTrue(pd.isna(df.loc[2, "smiles"]))
+        self.assertTrue(pd.isna(df.loc[3, "smiles"]))
+
+
+if __name__ == "__main__":
+    unittest.main()

nextflow/bin/tests/test_utils_glycan_offline.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+
+"""
+Small tests for the offline glycan-to-SMILES helper added to utils.py,
+and for how process_all_pdb_contacts.get_sugar_smiles_from_wurcs layers
+it in front of the existing live GlycoSmos/CSDB chain. No network access
+required (the live-chain calls are monkeypatched, not actually made).
+
+    python3 nextflow/bin/tests/test_utils_glycan_offline.py
+"""
+
+import sys
+import unittest
+from pathlib import Path
+from unittest import mock
+
+import numpy as np
+import pandas as pd
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+
+from utils import get_smiles_from_wurcs_offline
+import process_all_pdb_contacts as pac
+
+CHITOBIOSE_WURCS = "WURCS=2.0/1,2,1/[a2122h-1b_1-5_2*NCC/3=O]/1-1/a4-b1"
+
+
+class TestGetSmilesFromWurcsOffline(unittest.TestCase):
+
+    def test_valid_wurcs_returns_smiles(self):
+        smiles = get_smiles_from_wurcs_offline(CHITOBIOSE_WURCS)
+        self.assertIsInstance(smiles, str)
+        self.assertNotEqual(smiles, "")
+
+    def test_none_returns_nan(self):
+        self.assertTrue(pd.isna(get_smiles_from_wurcs_offline(None)))
+
+    def test_nan_returns_nan(self):
+        self.assertTrue(pd.isna(get_smiles_from_wurcs_offline(np.nan)))
+
+    def test_malformed_wurcs_returns_nan_not_exception(self):
+        # a production pipeline can't have a bad WURCS string crash the
+        # whole run - failures must degrade to nan, same contract as the
+        # existing get_glycoct_from_wurcs/get_csdb_from_glycoct/
+        # get_smiles_from_csdb functions in utils.py
+        self.assertTrue(pd.isna(get_smiles_from_wurcs_offline("not a wurcs string")))
+
+
+class TestContextBFallbackWiring(unittest.TestCase):
+    """
+    process_all_pdb_contacts.get_sugar_smiles_from_wurcs should try the
+    offline route first and only fall back to the live chain
+    (get_glycoct_from_wurcs -> get_csdb_from_glycoct -> get_smiles_from_csdb)
+    when the offline route fails - per the Context B recommendation in
+    docs/iupac_translator_plan.md (offline primary, live chain fallback).
+    """
+
+    def _empty_cache(self, columns):
+        return pd.DataFrame(columns=columns)
+
+    def test_offline_success_skips_live_chain_entirely(self):
+        with mock.patch.object(pac, "get_glycoct_from_wurcs") as mock_glycoct:
+            sugar_smiles, *_ = pac.get_sugar_smiles_from_wurcs(
+                [CHITOBIOSE_WURCS],
+                self._empty_cache(["glycoct", "csdb"]),
+                self._empty_cache(["csdb", "descriptor"]),
+                self._empty_cache(["WURCS", "glycoct"]),
+            )
+        self.assertFalse(pd.isna(sugar_smiles[CHITOBIOSE_WURCS]))
+        mock_glycoct.assert_not_called()
+
+    def test_offline_failure_falls_back_to_live_chain(self):
+        bad_wurcs = "not a wurcs string"
+        with mock.patch.object(pac, "get_glycoct_from_wurcs", return_value="fake_glycoct") as mock_glycoct, \
+             mock.patch.object(pac, "get_csdb_from_glycoct", return_value="fake_csdb") as mock_csdb, \
+             mock.patch.object(pac, "get_smiles_from_csdb", return_value="C") as mock_smiles:
+            sugar_smiles, *_ = pac.get_sugar_smiles_from_wurcs(
+                [bad_wurcs],
+                self._empty_cache(["glycoct", "csdb"]),
+                self._empty_cache(["csdb", "descriptor"]),
+                self._empty_cache(["WURCS", "glycoct"]),
+            )
+        mock_glycoct.assert_called_once()
+        mock_csdb.assert_called_once()
+        mock_smiles.assert_called_once()
+        self.assertEqual(sugar_smiles[bad_wurcs], "C")
+
+
+if __name__ == "__main__":
+    unittest.main()

nextflow/bin/tests/test_wurcs_to_iupac.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+"""
+Small, offline-only regression tests for wurcs_to_iupac.translate().
+Locks in the cases validated during development (see
+docs/iupac_translator_plan.md) so future changes can't silently regress
+them. No network access required - run with:
+
+    python3 nextflow/bin/tests/test_wurcs_to_iupac.py
+"""
+
+import sys
+import unittest
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+
+from wurcs_to_iupac import translate, TranslationError, SUBSTITUENT_MAP, _substituent_token
+from glypy.structure.substituent import Substituent
+
+
+class TestTranslate(unittest.TestCase):
+
+    def test_simple_disaccharide(self):
+        # chitobiose core, pulled from real PDB entry 6VXX
+        wurcs = "WURCS=2.0/1,2,1/[a2122h-1b_1-5_2*NCC/3=O]/1-1/a4-b1"
+        self.assertEqual(translate(wurcs), "Glc2NAc(b1-4)Glc2NAc")
+
+    def test_branched_n_glycan_core(self):
+        # Man3GlcNAc2 - the standard N-glycosylation core
+        wurcs = (
+            "WURCS=2.0/3,5,4/[a2122h-1b_1-5_2*NCC/3=O][a1122h-1a_1-5][a1122h-1a_1-5]"
+            "/1-1-2-3-3/a4-b1_b4-c1_c3-d1_c6-e1"
+        )
+        self.assertEqual(
+            translate(wurcs),
+            "Man(a1-6)[Man(a1-3)]Man(a1-4)Glc2NAc(b1-4)Glc2NAc",
+        )
+
+    def test_furanose_ring_suffix(self):
+        # sucrose - regression test for the "Fructofuranose" vs "Fruf"
+        # naming quirk found during PDB-scale testing
+        wurcs = "WURCS=2.0/2,2,1/[ha122h-2b_2-5][a2122h-1a_1-5]/1-2/a2-b1"
+        self.assertEqual(translate(wurcs), "Glc(a1-2)Fruf")
+
+    def test_undefined_anomer_reducing_end(self):
+        # a free reducing end has no fixed anomeric configuration in real
+        # WURCS data - regression test for the anomer-override retry in
+        # _base_sugar_name (this used to raise TranslationError entirely)
+        wurcs = "WURCS=2.0/1,1,1/[a2122h-1x_1-5]/1/"
+        result = translate(wurcs)
+        self.assertTrue(result.startswith("Glc"))
+
+    def test_unsupported_substituent_raises(self):
+        # n_methyl is confirmed-unsupported (composes to the wrong
+        # chemistry in GlyLES, see docs/iupac_translator_plan.md) and
+        # must not be in the map, and must fail loudly rather than guess
+        self.assertNotIn("n_methyl", SUBSTITUENT_MAP)
+        with self.assertRaises(TranslationError):
+            _substituent_token(2, Substituent("n_methyl"))
+
+    def test_malformed_wurcs_does_not_silently_succeed(self):
+        with self.assertRaises(Exception):
+            translate("not a wurcs string")
+
+
+if __name__ == "__main__":
+    unittest.main()

nextflow/bin/utils.py

@@ -6,12 +6,15 @@
 import requests
 from urllib.parse import quote
 from bs4 import BeautifulSoup
-import json 
+import json
 from Bio.ExPASy import Enzyme as EEnzyme
 from pdbeccdutils.helpers.mol_tools import fix_molecule
 from rdkit import Chem
-import gzip 
+import gzip
 import xml.etree.ElementTree as ET
+import signal
+import glyles
+from wurcs_to_iupac import translate as translate_wurcs_to_iupac
 
 #make this function be applicable to extract_pdbe_info script too - need to check if the grouped output at end is appropriate.
 def process_ec_records(enzyme_dat_file, enzyme_class_file):
@@ -149,7 +152,47 @@ def get_smiles_from_csdb(csdb_linear, cache_df):
         else:
             smiles = np.nan
         return smiles
-
+
+def get_smiles_from_wurcs_offline(wurcs, timeout_seconds = 15):
+    """
+    Convert a WURCS glycan descriptor straight to SMILES, entirely offline:
+    glypy (WURCS parsing) -> wurcs_to_iupac.translate (IUPAC-condensed
+    translation) -> glyles (SMILES generation). No network calls, unlike
+    get_glycoct_from_wurcs/get_csdb_from_glycoct/get_smiles_from_csdb above.
+
+    See docs/iupac_translator_plan.md for the full validation history -
+    benchmarked at ~89% success on the real cognate-ligand master set and
+    ~96% on real PDB-deposited glycans, vs. 0% for the live glycoct-based
+    chain above (GlyTouCan's API no longer returns glycoct at all).
+
+    GlyLES is known to hang on malformed input instead of failing fast
+    (an ANTLR grammar issue, not specific to this translator), hence the
+    timeout - returns np.nan rather than blocking indefinitely.
+    """
+    if wurcs is np.nan or wurcs is None:
+        return np.nan
+    try:
+        iupac = translate_wurcs_to_iupac(wurcs)
+    except Exception:
+        return np.nan
+
+    def _timeout_handler(signum, frame):
+        raise TimeoutError("glyles.convert timed out")
+
+    previous_handler = signal.signal(signal.SIGALRM, _timeout_handler)
+    signal.alarm(timeout_seconds)
+    try:
+        result = glyles.convert(glycan = iupac, verbose = None)
+    except Exception:
+        return np.nan
+    finally:
+        signal.alarm(0)
+        signal.signal(signal.SIGALRM, previous_handler)
+
+    if not result or not result[0][1]:
+        return np.nan
+    return result[0][1]
+
 def pdbe_sanitise_smiles(smiles, return_mol = False, return_sanitisation = False):
     """
     Sanitises a smiles string using pdbeccdutils fix_molecule functions and