Identifying rings in system with fused rings.

[thesonofmoon]

Hi there,

Currently I am trying to write a code that is able to identify rings and determine how large the rings are and whether they are aromatic or not. What this means is that I want the output of the code to be for one the non-aromatic rings, if they are fused or not, and how many aromatic rings there are and if they are fused give their trivial name like naphthalene or anthracene.

print("Non-aromatic rings found:", unique_ring_matches)
print("Fused rings found:", fused_ring_matches)
print("Aromatic rings found:", aromatic_ring_matches)

I truly appreciate your help!

Here my initial ideas, so if you want to modify it so it works for example for anthracene or Linagliptin

#finding non aromatic rings, doesn't work for annulated rings
from rdkit import Chem

def find_unique_non_aromatic_rings(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return "Invalid SMILES string. Unable to parse molecule."

    # Define SMARTS patterns to recognize rings of specific sizes
    ring_smarts = {
        '3_membered_ring': '[r3]',
        '4_membered_ring': '[r4]',
        '5_membered_ring': '[r5]',
        '6_membered_ring': '[r6]',
        '7_membered_ring': '[r7]',
        '8_membered_ring': '[r8]',
        '9_membered_ring': '[r9]',
        '10_membered_ring': '[r10]',
    }

    # Define a list to store the chemical groups found in the SMILES
    chemical_groups = []

    # Find matches of each ring pattern in the molecule
    for name, pattern in ring_smarts.items():
        ring_size = int(name.split('_')[0])  # Extract the ring size from the name
        ring_matches = mol.GetSubstructMatches(Chem.MolFromSmarts(pattern))
        non_aromatic_rings = [match for match in ring_matches if not any(mol.GetAtomWithIdx(idx).GetIsAromatic() for idx in match)] #exclude all aromatic rings
        unique_rings_count = len(non_aromatic_rings) // ring_size  # Correctly count unique rings
        if unique_rings_count > 0:
            chemical_groups.extend([name] * unique_rings_count)

    # Return the list of identified chemical groups
    return chemical_groups

# Example alpha-Hopan
smiles_string = 'CC(C)[C@H]1CC[C@]2([C@@H]1CC[C@@]3([C@@H]2CC[C@H]4[C@]3(CC[C@@H]5[C@@]4(CCCC5(C)C)C)C)C)C'
unique_ring_matches = find_unique_non_aromatic_rings(smiles_string)
aromatic_ring_matches = find_aromatic_ring(smiles_string)
print("Non-aromatic rings found:", unique_ring_matches)
print("Aromatic rings found:", aromatic_ring_matches)

# Problems Linagliptin
smiles_string = 'CC1=C(C(=O)C2=CC=CC=C2C1=O)C/C=C(\C)/CCC[C@H](C)CCC[C@H](C)CCCC(C)C'
unique_ring_matches = find_unique_non_aromatic_rings(smiles_string)
aromatic_ring_matches = find_aromatic_ring(smiles_string)
print("Non-aromatic rings found:", unique_ring_matches)
print("Aromatic rings found:", aromatic_ring_matches)

[bp-kelley]

It might be easier to use the mol.GetRingInfo and just prune the molecule

from rdkit import Chem

def prune_rings(mol, aromatic):
    atoms_to_keep = set()
    for atom_ring in mol.GetRingInfo().AtomRings():
        for atom_idx in atom_ring:
            if mol.GetAtomWithIdx(atom_idx).GetIsAromatic() == aromatic:
                if aromatic == False: # if any atom is non aromatic keep the whole ring
                    for atom_idx in atom_ring:
                        atoms_to_keep.add(atom_idx)
                    break
                else:
                    atoms_to_keep.add(atom_idx)

    mol.BeginBatchEdit()                        
    for atom in mol.GetAtoms():
        if atom.GetIdx() not in atoms_to_keep:
            mol.RemoveAtom(atom.GetIdx())
    mol.CommitBatchEdit()

    return Chem.GetMolFrags(mol, asMols=True)

smiles_strings = ('CC(C)[C@H]1CC[C@]2([C@@H]1CC[C@@]3([C@@H]2CC[C@H]4[C@]3(CC[C@@H]5[C@@]4(CCCC5(C)C)C)C)C)C',
                  'CC1=C(C(=O)C2=CC=CC=C2C1=O)C/C=C(\C)/CCC[C@H](C)CCC[C@H](C)CCCC(C)C')
for smiles_string in smiles_strings:
    print("smiles", smiles_string)
    m = Chem.MolFromSmiles(smiles_string)
    
    aromatic = Chem.RWMol(m)
    frags = prune_rings(aromatic, True)
    print("aromatic", Chem.MolToSmiles(aromatic))
    print("frags", [Chem.MolToSmiles(x) for x in frags])
    
    nonaromatic = Chem.RWMol(m)
    frags = prune_rings(nonaromatic, False)
    print("non aromatic", Chem.MolToSmiles(nonaromatic))
    print("frags", [Chem.MolToSmiles(x) for x in frags])

output

smiles CC(C)[C@H]1CC[C@]2([C@@H]1CC[C@@]3([C@@H]2CC[C@H]4[C@]3(CC[C@@H]5[C@@]4(CCCC5(C)C)C)C)C)C
aromatic 
frags []
non aromatic C1CC[C@@H]2CCC3C4CC[C@@H]5CCCC5[C@H]4CC[C@@H]3C2C1
frags ['[CH]1CC[C]2[C@H]3CC[C@H]4[C](CC[C@H]5CCCC[C]54)[C]3CC[C@H]12']

smiles CC1=C(C(=O)C2=CC=CC=C2C1=O)C/C=C(\C)/CCC[C@H](C)CCC[C@H](C)CCCC(C)C
aromatic c1ccccc1
frags ['c1ccccc1']
non aromatic C1=CCccC1
frags ['C1=CCC=CC1']

You have one case where a non aromatic is fused to an aromatic ring, I'm assuming you want to report this as non-aromatic.

[thesonofmoon]

Thank you so much for your help! Yes that is correct!