Merge pull request #217 from datamol-io/feat/isomers

Added function to get the number of stereoisomers

Author: zhu0619

datamol/__init__.py

@@ -135,6 +135,7 @@
     "save_df": "datamol.io",
     # isomers
     "enumerate_stereoisomers": "datamol.isomers",
+    "count_stereoisomers": "datamol.isomers",
     "enumerate_tautomers": "datamol.isomers",
     "enumerate_structisomers": "datamol.isomers",
     "canonical_tautomer": "datamol.isomers",
@@ -329,6 +330,7 @@ def __dir__():
     from .io import open_df
 
     from .isomers import enumerate_stereoisomers
+    from .isomers import count_stereoisomers
     from .isomers import enumerate_tautomers
     from .isomers import enumerate_structisomers
     from .isomers import canonical_tautomer

datamol/isomers/__init__.py

@@ -3,5 +3,6 @@
 from ._enumerate import enumerate_stereoisomers
 from ._enumerate import enumerate_tautomers
 from ._enumerate import enumerate_structisomers
+from ._enumerate import count_stereoisomers
 from ._enumerate import remove_stereochemistry
 from ._enumerate import canonical_tautomer

datamol/isomers/_enumerate.py

@@ -5,6 +5,7 @@
 
 from rdkit.Chem.MolStandardize import rdMolStandardize
 from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers
+from rdkit.Chem.EnumerateStereoisomers import GetStereoisomerCount
 from rdkit.Chem.EnumerateStereoisomers import StereoEnumerationOptions
 
 import datamol as dm
@@ -82,7 +83,8 @@ def enumerate_stereoisomers(
 
     # in case any bonds/centers are missing stereo chem flag it here
     Chem.AssignStereochemistry(mol, force=False, flagPossibleStereoCenters=True, cleanIt=clean_it)  # type: ignore
-    Chem.FindPotentialStereoBonds(mol, cleanIt=clean_it)  # type: ignore
+    # lu: do not clean (overwrite bond stereo information) when set `undefined_only=Ture`
+    Chem.FindPotentialStereoBonds(mol, cleanIt=not undefined_only and clean_it)
 
     # set up the options
     stereo_opts = StereoEnumerationOptions(
@@ -117,6 +119,68 @@ def enumerate_stereoisomers(
     return variants
 
 
+def count_stereoisomers(
+    mol: dm.Mol,
+    n_variants: int = 20,
+    undefined_only: bool = False,
+    rationalise: bool = True,
+    timeout_seconds: int = None,
+    clean_it: bool = True,
+    precise: bool = False,
+):
+    """Get the number of possible stereoisomers for a molecule.
+
+    Warning: By default, this function compute an estimtion number based on the stereo bonds which
+    gives an upper bound of possible stereoisomers. By setting `precise=True`, the number is computed
+    by enumrerating the stereoisomers. However, it can be computationnaly intensive.
+
+    Args:
+        mol: The molecule whose state we should enumerate.
+        n_variants: The maximum amount of molecules that should be returned.
+        undefined_only: If we should enumerate all stereocenters and bonds or only those
+            with undefined stereochemistry.
+        rationalise: If we should try to build and rationalise the molecule to ensure it
+            can exist.
+        timeout_seconds: The maximum amount of time to spend on enumeration. None
+            will disable the timeout. Note that the timeout might be inaccurate as a running single variant
+            computation is not stopped when the duration is reached.
+        clean_it: A flag for assigning stereochemistry. If True, it will remove previous stereochemistry
+            markings on the bonds.
+        precise: Whether compute counts by enumerate the stereoisomers using `enumerate_stereoisomers`.
+
+    """
+    if precise:
+        num_variants = len(
+            enumerate_stereoisomers(
+                mol=mol,
+                n_variants=n_variants,
+                undefined_only=undefined_only,
+                rationalise=rationalise,
+                timeout_seconds=timeout_seconds,
+                clean_it=clean_it,
+            )
+        )
+    else:
+        # safety first
+        mol = dm.copy_mol(mol)
+
+        # in case any bonds/centers are missing stereo chem flag it here
+        Chem.AssignStereochemistry(mol, force=False, flagPossibleStereoCenters=True, cleanIt=clean_it)  # type: ignore
+        # lu: do not clean (overwrite bond stereo information) when set `undefined_only=Ture`
+        Chem.FindPotentialStereoBonds(mol, cleanIt=not undefined_only and clean_it)
+
+        # set up the options
+        stereo_opts = StereoEnumerationOptions(
+            tryEmbedding=rationalise,
+            onlyUnassigned=undefined_only,
+            unique=True,
+        )
+
+        num_variants = GetStereoisomerCount(mol, options=stereo_opts)
+
+    return num_variants
+
+
 def enumerate_structisomers(
     mol: dm.Mol,
     n_variants: int = 20,

tests/test_isomers.py

@@ -43,6 +43,14 @@ def test_enumerate_stereo_timeout():
     dm.enumerate_stereoisomers(mol, n_variants=2, timeout_seconds=1)
 
 
+def test_count_stereoisomers():
+    num_isomers_1 = dm.count_stereoisomers(dm.to_mol("CC=CC"), undefined_only=True)
+    num_isomers_2 = dm.count_stereoisomers(dm.to_mol("CC=CC"), undefined_only=False)
+    assert num_isomers_1 == num_isomers_2
+
+    assert dm.count_stereoisomers(dm.to_mol("Br/C=C\\Br"), undefined_only=True) == 1
+
+
 def test_enumerate_structural():
     mol = dm.to_mol("CCCCC")  # pentane has only three structural isomers