Python and RDKit to extract sub-structures in a SMILES

Question

I am not a chemist, I come from a computer science background. I am involved in a cheminformatics project. I have a list of molecules in SMILES format, which I want to extract sub-structures from.

I found some interesting sub-structures that I want to look for in a given molecule in the paper https://www.mdpi.com/1420-3049/21/1/75 published in 2015.

In the abstract it is mentioned that: We first uncovered three important criteria closely related to drug-likeness, namely:

• (1) the best numbers of aromatic and non-aromatic rings are 2 and 1, respectively;

• (2) the best functional groups of candidate drugs are usually -OH, -COOR and -COOH in turn, but not -CONHOH, -SH, -CHO and -SO3H. In addition, the -F functional group is beneficial to CNS drugs, and -NH2 functional group is beneficial to anti-infective drugs and anti-cancer drugs;

• (3) the best R-value intervals of candidate drugs are in the range of 0.05–0.50 (preferably 0.10–0.35), and R-value of the candidate CNS drugs should be as small as possible in this interval

• also, Acyclic groups as mentioned in this paper (figure 4), https://bmcchem.biomedcentral.com/articles/10.1186/s13065-021-00737-2

My question is, is there a method in RDKit (or other Python libraries) that extracts these criteria from smiles?

Answer 1

I chose the example in Figure 1 in the article for your first three points.

from rdkit import Chem
from rdkit.Chem import Draw, Descriptors, rdqueries

m = Chem.MolFromSmiles('NC(=O)C1(CCN(CCCN2C3=C(CCC4=C2C=C(Cl)C=C4)C=CC=C3)CC1)N1CCCCC1')
Draw.MolToImage(m, size=(400,200))

print('Aromatic ring count =', Descriptors.NumAromaticRings(m))
Aromatic ring count = 2

print('Non-aromatic ring count =', Descriptors.NumAliphaticRings(m))
Non-aromatic ring count = 3

You have to define SMARTS for all functional groups.

fg = Chem.MolFromSmarts('C(=O)[NX3;H2]') # SMARTS for -CONH2
print('Functional group:', len(m.GetSubstructMatches(fg)), '-CONH2')
Functional group: 1 -CONH2

If I see it correctly, the R value is (heavyatoms - carbons) / heavyatoms.

heavyatoms = Descriptors.HeavyAtomCount(m)
q = rdqueries.AtomNumEqualsQueryAtom(6) # 6 for carbon
carbons = len(m.GetAtomsMatchingQuery(q))
r = (heavyatoms-carbons)/heavyatoms
print('R value =', round(r, 2))
R value = 0.18

What you call Acyclic groups just means, that there are no rings in molecule.

from rdkit.Chem.Scaffolds import MurckoScaffold

m1 = Chem.MolFromSmiles('CCC')
core = MurckoScaffold.GetScaffoldForMol(m1)
s = Chem.MolToSmiles(core)
if len(s) == 0:
    print('No ring in the molecule')
No ring in the molecule