Listing all possible structures for a given Markush structure automatically

Question

I have Markush formulas described as mol blocks, v3000, even if that does not matter that much because I could convert it into another format if needed.

I need the possibility to automatically list all possible structures for a given Markush. More specifically I often have aromatic rings with a substituent that can be in any position, and I'd like some automated way of listing the structures for that.

I realize there is this Knime node called "Markush Enumeration" (https://docs.chemaxon.com/display/docs/markush-enumeration.md), but it needs a license.

Is there a 'free' way to easily do it? I am using Python and Knime, if that matters.

Answer 1

Preface: Since you mention familiarity with/access to Python, and I assume you want to process batches from .sdf from a script rather than from a GUI (note, ChemAxon offers a dedicated Markush Editor in addition to a plugin for their sketcher, Marvin (a tutorial)). Thus, I suggest using RDKit. In case you have access to a Linux distribution of Debian (or a relative of Debian like Ubuntu), there is a repackaged version by DebiChem (tracker), too. Then you do not need to install a separate instance of Miniconda/Anaconda (some notes about this on Matter Modeling Stack Exchange).

The following works well within a Jupyter notebook by Linux Debian 12/bookworm (branch testing) offering Python (version 3.9.9) and RDKit (version 2021.09.2). The credit is due to Greg Landrum, who posted the bits and bolts required for this to work and only rearranged/paraphrased for this answer.

Ad rem. In the first cell, initiate the Jupyter notebook:

from rdkit import Chem
from rdkit.Chem.Draw import rdDepictor
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import rdMolEnumerator
import rdkit
print(rdkit.__version__)

Next, load .sdf to process, e.g.

pv2 = Chem.MolFromMolBlock('''
  Mrv2007 06242006032D          

  0  0  0     0  0            999 V3000
M  V30 BEGIN CTAB
M  V30 COUNTS 10 8 0 0 0
M  V30 BEGIN ATOM
M  V30 1 C -1.7083 2.415 0 0
M  V30 2 C -3.042 1.645 0 0
M  V30 3 C -3.042 0.105 0 0
M  V30 4 N -1.7083 -0.665 0 0
M  V30 5 C -0.3747 0.105 0 0
M  V30 6 C -0.3747 1.645 0 0
M  V30 7 * -3.042 0.875 0 0
M  V30 8 F -5.0434 0.875 0 0
M  V30 9 * -1.0415 2.03 0 0
M  V30 10 Cl -1.0415 4.34 0 0
M  V30 END ATOM
M  V30 BEGIN BOND
M  V30 1 1 1 2
M  V30 2 2 2 3
M  V30 3 1 3 4
M  V30 4 2 4 5
M  V30 5 1 5 6
M  V30 6 2 1 6
M  V30 7 1 7 8 ENDPTS=(2 2 3) ATTACH=ANY
M  V30 8 1 9 10 ENDPTS=(2 1 6) ATTACH=ANY
M  V30 END BOND
M  V30 END CTAB
M  END
''')
pv2

to yield

where grey clouds mark where variation of the position of the substituents is anticipated. By subsequent

pv2_bundle = rdMolEnumerator.Enumerate(pv2)
Draw.MolsToGridImage(pv2_bundle)

you launch the variation of the positions, with the results displayed as

If you want, a subsequent list comprehension offers you to access the result e.g., in SMILES strings (export as .sdf (v3000) equally is possible, check RDKit's documentation and cookbook):

listing = []
listing = [rdkit.Chem.MolToSmiles(entry) for entry in pv2_bundle]

for element in listing:
    print(element)

with output of

Fc1cnccc1Cl
Fc1cncc(Cl)c1
Fc1cc(Cl)ccn1
Fc1ccc(Cl)cn1

Again, the above only is rearranging an example by Greg Landrum's blog post detailing out multiple additional equally interesting examples of application.

---

As for the generation of the .sdf (v3000) with «Markush syntax» with a GUI, beside mentioned Marvin by ChemAxon, BIOVIA Draw by Dassault offers this functionality, too. The later bears some similarity to MDL's ISIS Draw the more senior readers of chemistry.se still might remember/use.