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?

  • $\begingroup$ Are you asking about only substructure search, or all of the 4 criteria? Because I believe that you would need some kind of numerical analysis to get those criteria. $\endgroup$
    – S R Maiti
    Feb 21, 2022 at 21:42
  • $\begingroup$ No, each substructure/criteria alone. Not all of them at once. $\endgroup$
    – mac179
    Feb 22, 2022 at 9:54

1 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.