Merge pull request #499 from Steinbeck-Lab/development

feat: substructure highlighting in depict

Dockerfile

@@ -1,10 +1,8 @@
 FROM continuumio/miniconda3:24.1.2-0 AS cheminf-python-ms
 
-# Set environment variables
-ENV PYTHON_VERSION=3.10 \
-    RDKIT_VERSION=2023.09.4 \
-    OPENBABEL_VERSION=v3.1.1 \
-    JAVA_HOME=/usr/lib/jvm/java-11-openjdk-amd64/
+ENV PYTHON_VERSION=3.10
+ENV RDKIT_VERSION=2023.09.4
+ENV OPENBABEL_VERSION=v3.1.1
 
 # Install runtime dependencies
 RUN apt-get update && \
@@ -44,4 +42,4 @@ RUN pip3 install --no-cache-dir chembl_structure_pipeline --no-deps
 
 COPY ./app /code/app
 
-CMD ["uvicorn", "app.main:app", "--host", "0.0.0.0", "--port", "80", "--workers", "8"]
+CMD ["uvicorn", "app.main:app", "--host", "0.0.0.0", "--port", "80", "--workers", "2"]

app/modules/depiction.py

@@ -15,8 +15,10 @@ def get_cdk_depiction(
     molecule: any,
     molSize=(512, 512),
     rotate=0,
+    kekulize=True,
     CIP=True,
     unicolor=False,
+    highlight="",
 ):
     """This function takes the user input SMILES and Depicts it.
 
@@ -28,6 +30,7 @@ def get_cdk_depiction(
     Returns:
         image (SVG): CDK Structure Depiction as an SVG image.
     """
+    print(unicolor)
 
     cdk_base = "org.openscience.cdk"
     StandardGenerator = JClass(
@@ -36,6 +39,10 @@ def get_cdk_depiction(
     Color = JClass("java.awt.Color")
     UniColor = JClass(cdk_base + ".renderer.color.UniColor")
     CDK2DAtomColors = JClass(cdk_base + ".renderer.color.CDK2DAtomColors")()
+    Kekulization = JClass(cdk_base + ".aromaticity.Kekulization")
+    SmartsPattern = JClass(cdk_base + ".smarts.SmartsPattern")
+    SCOB = JClass(cdk_base + ".silent.SilentChemObjectBuilder")
+
     if unicolor:
         DepictionGenerator = (
             JClass(cdk_base + ".depict.DepictionGenerator")()
@@ -63,6 +70,12 @@ def get_cdk_depiction(
 
     if SDGMol:
         # Rotate molecule
+        if kekulize:
+            try:
+                Kekulization.kekulize(SDGMol)
+            except Exception as e:
+                print(e + "Can't Kekulize molecule")
+
         point = JClass(
             cdk_base + ".geometry.GeometryTools",
         ).get2DCenter(SDGMol)
@@ -72,6 +85,16 @@ def get_cdk_depiction(
             (rotate * JClass("java.lang.Math").PI / 180.0),
         )
 
+        if highlight and highlight.strip():
+            tmpPattern = SmartsPattern.create(highlight, SCOB.getInstance())
+            SmartsPattern.prepare(SDGMol)
+            tmpMappings = tmpPattern.matchAll(SDGMol)
+            tmpSubstructures = tmpMappings.toSubstructures()
+            lightBlue = Color(173, 216, 230)
+            DepictionGenerator = DepictionGenerator.withHighlight(
+                tmpSubstructures, lightBlue
+            ).withOuterGlowHighlight()
+
         mol_imageSVG = (
             DepictionGenerator.depict(
                 SDGMol,
@@ -89,32 +112,68 @@ def get_cdk_depiction(
         return "Error reading SMILES string, check again."
 
 
-def get_rdkit_depiction(molecule: any, molSize=(512, 512), rotate=0, kekulize=True):
-    """This function takes the user input SMILES and Canonicalize it using the.
-
-    RDKit.
+def get_rdkit_depiction(
+    molecule: Chem.Mol,
+    mol_size=(512, 512),
+    rotate=0,
+    kekulize=True,
+    CIP=False,
+    unicolor=False,
+    highlight: str = "",
+) -> str:
+    """
+    Generate a 2D depiction of the input molecule using RDKit.
 
     Args:
         molecule (Chem.Mol): RDKit molecule object.
+        mol_size (tuple, optional): Size of the output image. Defaults to (512, 512).
+        rotate (int, optional): Rotation angle of the molecule. Defaults to 0.
+        kekulize (bool, optional): Whether to kekulize the molecule. Defaults to True.
+        CIP (bool, optional): Whether to assign CIP stereochemistry. Defaults to False.
+        unicolor (bool, optional): Whether to use a unicolor palette. Defaults to False.
+        highlight (str, optional): SMARTS pattern to highlight atoms/bonds. Defaults to empty.
 
     Returns:
-        image (SVG): RDKit Structure Depiction as an SVG image.
+        str: RDKit Structure Depiction as an SVG image.
     """
-    if molecule:
-        mc = Chem.Mol(molecule.ToBinary())
-        if kekulize:
-            try:
-                Chem.Kekulize(mc)
-            except Exception as e:
-                print(e, flush=True)
-                mc = Chem.Mol(molecule.ToBinary())
-        if not mc.GetNumConformers():
-            rdDepictor.Compute2DCoords(mc)
-        drawer = rdMolDraw2D.MolDraw2DSVG(molSize[0], molSize[1])
-        drawer.drawOptions().rotate = rotate
-        drawer.DrawMolecule(mc)
-        drawer.FinishDrawing()
-        svg = drawer.GetDrawingText()
-        return svg.replace("svg:", "")
+
+    mc = Chem.Mol(molecule.ToBinary())
+
+    if kekulize:
+        try:
+            Chem.Kekulize(mc)
+        except Chem.KekulizeException:
+            mc = Chem.Mol(molecule.ToBinary())
+
+    if not mc.GetNumConformers():
+        rdDepictor.Compute2DCoords(mc)
+
+    if CIP:
+        Chem.AssignStereochemistry(mc, force=True, cleanIt=True)
+
+    drawer = rdMolDraw2D.MolDraw2DSVG(mol_size[0], mol_size[1])
+    drawer.drawOptions().rotate = rotate
+    drawer.drawOptions().addStereoAnnotation = CIP
+
+    if unicolor:
+        drawer.drawOptions().useBWAtomPalette()
+
+    if highlight:
+        patt = Chem.MolFromSmarts(highlight)
+        if patt:
+            hit_ats = mc.GetSubstructMatch(patt)
+            hit_bonds = [
+                mc.GetBondBetweenAtoms(at1, at2).GetIdx()
+                for at1, at2 in zip(hit_ats[:-1], hit_ats[1:])
+            ]
+            rdMolDraw2D.PrepareAndDrawMolecule(
+                drawer, mc, highlightAtoms=hit_ats, highlightBonds=hit_bonds
+            )
+        else:
+            drawer.DrawMolecule(mc)
     else:
-        return "Error reading SMILES string, check again."
+        drawer.DrawMolecule(mc)
+
+    drawer.FinishDrawing()
+    svg = drawer.GetDrawingText()
+    return svg.replace("svg:", "")

app/modules/toolkits/rdkit_wrapper.py

@@ -58,7 +58,7 @@ def get_MolVolume(molecule: any) -> float:
         float: The volume of the molecule.
     """
     molecule = Chem.AddHs(molecule)
-    AllChem.EmbedMolecule(molecule)
+    AllChem.EmbedMolecule(molecule, useRandomCoords=True)
     volume = AllChem.ComputeMolVolume(molecule, gridSpacing=0.2)
     return volume
 

app/routers/depict.py

@@ -108,16 +108,21 @@ async def depict_2d_molecule(
         title="Rotate",
         description="The rotation angle of the molecule in degrees.",
     ),
-    CIP: Optional[bool] = Query(
+    CIP: bool = Query(
         False,
         title="CIP",
         description="Whether to include Cahn-Ingold-Prelog (CIP) stereochemistry information.",
     ),
-    unicolor: Optional[bool] = Query(
+    unicolor: bool = Query(
         False,
         title="Unicolor",
         description="Whether to use a single colour for the molecule.",
     ),
+    highlight: Optional[str] = Query(
+        "COSN",
+        title="Substructure",
+        description="SMARTS pattern to highlight atoms/bonds.",
+    ),
 ):
     """Generates a 2D depiction of a molecule using CDK or RDKit with the given.
 
@@ -132,6 +137,7 @@ async def depict_2d_molecule(
     - **rotate**: (int, optional): The rotation angle of the molecule in degrees. Defaults to 0.
     - CIP (bool, optional): Whether to include Cahn-Ingold-Prelog (CIP) stereochemistry information. Defaults to False.
     - unicolor (bool, optional): Whether to use a single colour for the molecule. Defaults to False.
+    - highlight (Optional[str], optional): SMARTS pattern to highlight atoms/bonds. Defaults to "COSN".
 
     Returns:
         Response: An HTTP response containing the generated image in SVG+xml format.
@@ -155,12 +161,15 @@ async def depict_2d_molecule(
                 mol,
                 [width, height],
                 rotate,
-                CIP,
-                unicolor,
+                CIP=CIP,
+                unicolor=unicolor,
+                highlight=highlight,
             )
         elif toolkit == "rdkit":
             mol = parse_input(smiles, "rdkit", False)
-            depiction = get_rdkit_depiction(mol, [width, height], rotate)
+            depiction = get_rdkit_depiction(
+                mol, [width, height], rotate, unicolor=unicolor, highlight=highlight
+            )
         else:
             raise HTTPException(
                 status_code=422,

tests/test_functions.py

@@ -88,7 +88,7 @@ def test_get_rdkit_depiction_rotate(test_RDKit_Mol):
 
 
 def test_get_rdkit_depiction_size(test_RDKit_Mol):
-    svg = get_rdkit_depiction(test_RDKit_Mol, molSize=(512, 512))
+    svg = get_rdkit_depiction(test_RDKit_Mol, mol_size=(512, 512))
     assert isinstance(svg, str)
     assert "svg" in svg
     assert "Error" not in svg
@@ -158,7 +158,7 @@ def test_all_rdkit_descriptors(test_smiles_descriptors):
         0,
         1.0,
         0,
-        62.15,
+        62.12,
     )
     assert expected_result == descriptors
 
@@ -210,7 +210,7 @@ def test_all_combined_descriptors(test_smiles_descriptors):
         "Formal Charge": (0, 0),
         "FractionCSP3": (1.0, 1.0),
         "Number of Minimal Rings": (0, 0),
-        "Van der Waals Volume": (62.11, 60.444412578400105),
+        "Van der Waals Volume": (62.01, 60.444412578400105),
     }
     assert expected_result == descriptors