Chem.SanitizeMol removes aromaticity tag from smarts Mol

[HelloJocelynLu]

Describe the bug
When creating a mol object with MolFromSmarts, the aromaticity tag of atoms is removed by Chem.SanitizeMol.

To Reproduce
Before the sanitization, GetIsAromatic() returns True; After Chem.SanitizeMol(mol), GetIsAromatic() returns False

mol = Chem.MolFromSmarts('a1aaaaa1')
for atom in mol.GetAtoms():
    print(atom.GetIdx(), atom.GetIsAromatic())
Chem.SanitizeMol(mol)
print('='*20)
for atom in mol.GetAtoms():
    print(atom.GetIdx(), atom.GetIsAromatic())

I can confirm that this issue does not occur if we remove Chem.SanitizeMol(mol) or if we have mol = Chem.MolFromSmiles('c1ccccc1')

Expected behavior
I would expect GetIsAromatic() should not change before & after sanitization:

0 True
1 True
2 True
3 True
4 True
5 True
====================
0 True
1 True
2 True
3 True
4 True
5 True

Screenshots

Configuration (please complete the following information):

• RDKit version: 2023.03.2
• OS: [e.g. Ubuntu 20.04]: Ubuntu 20.04 and MacOS Ventura 13.3.1 (a)
• Python version (if relevant): python 3.9.7
• Are you using conda? Yes
• If you are using conda, which channel did you install the rdkit from? conda-forge
• If you are not using conda: how did you install the RDKit? NA

Additional context
Add any other context about the problem here.

[bp-kelley]

I think you are asking a bit much from sanitization. Sanitization really requires a physically reasonable molecule.

That being said, the RDKit is largely based on the daylight smiles model, so let's take a look at what you are really doing, first by converting your query to a molecule and then to a smiles string:

>>> from rdkit import Chem
>>> m = Chem.MolFromSmarts("a1aaaaa1")
>>> smi = Chem.MolToSmiles(m)
>>> print(smi)
*1:*:*:*:*:*:1

What this is converted to in smiles-land are dummy atoms connected by aromatic bonds. In SMILES there is no concept of an aromatic dummy atom so this gets stripped. However, if we read the smiles we just made back into a smarts pattern:

>>> pat = Chem.MolFromSmarts(smi)
>>> print(Chem.MolFromSmiles("c1ccccc1").HasSubstructMatch(pat))
True

It still matches aromatic atoms so all is good. I think an important concept in SMARTS which is kind of the opposite in SMILES is that the matches are controlled by what you are explicitly set to, otherwise you still match. I.e. if you don't say you are explicitly aliphatic or aromatic, you will match both.

[HelloJocelynLu]

Hi Brian,

Thank you for your explanation. That's helpful!
Could you please review this discussion? I am interested in preserving certain properties, like aromaticity, in SMATRS. If you have any insights on this topic, I would greatly appreciate it.

Jocelyn

[greglandrum]

I've made this a discussion, because it's not a bug (the RDKit is behaving as intended), but the discussion is good to have.
@HelloJocelynLu sanitization is generally intended to do things like chemistry perception and correctness checking on molecules. This type of operation doesn't generally apply to query molecules like what you from SMARTS.

It may be possible to do what you want another way: why are you trying to sanitize the molecules constructed from SMARTS?

[HelloJocelynLu]

It is a little bit complicated, let me try to explain it.
I was trying to build a query and search fragments that satisfied the query, then replace a substructure with those qualified fragment. Below is a demo of this operation:

mol = Chem.MolFromSmiles("CCc1ccn(C)c1")
query_smarts = template_from_submol(mol, index_to_be_replaced, {4: "n", 5: "n"})    # query_smarts = "a1anna1"
query = Chem.MolFromSmarts(query_smarts)
# vocab = {"c1cn[nH]c1": Chem.rdmol object}
for k,v in vocab.items():
    if v.HasSubstructMatch(query): # v smiles: c1cn[nH]c1
        match = v.GetSubstructMatch(query) # (0,4,3,2,1)
        new_mol = mol_builder(...)

I have a molecule mol looks like:

I want to replace the 5-member aromatic ring (2,3,4,5,7) to have N in position 4 and 5, this gives me query_smarts = "a1anna1". Then I try to replace the 5-member aromatic ring (2,3,4,5,7) to c1cn[nH]c1 with function mol_builder. The content of mol_builder is:

def mol_builder(old_mol, smotif, atom_map):
    edmol = Chem.RWMol(old_mol)
    for aix in atom_map:
        atom = smotif.GetAtomWithIdx(atom_map[aix])
        added_atom = Chem.Atom(atom.GetSymbol())
        added_atom.SetFormalCharge(atom.GetFormalCharge())
        edmol.ReplaceAtom(aix, atom, preserveProps=True)
        
    for bond in old_mol.GetBonds():
        start, end = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
        if start in atom_map and end in atom_map:
            edmol.ReplaceBond(bond.GetIdx(), smotif.GetBondBetweenAtoms(atom_map[start], atom_map[end]))
    new_mol = edmol.GetMol()
    return new_mol

Basically the function just create a rwmol object, and replace substructure atom by atom, bond by bond based on old->new atom mapping. The resulting new_mol is shown below:

So occasionally I have problematic structures (like this) that cannot be sanitized. I am wondering if there is better way to handle/tackle this problem?

[HelloJocelynLu]

Hi Greg,

Thank you for your response. I really appreciate your work in contributing and supporting RDKit community!

I have been tried many different ways to do the replacement and generate valid molecules. The issue mentioned in the title is just one of the challenges I have encountered during my various explorations.
This is the minimum code that can generate old_mol, substructure to replace and new_mol:

mol = Chem.MolFromSmiles("CCc1ccn(C)c1")
index_to_be_replaced = (2,3,4,5,7)
query_smarts = "a1anna1"
query = Chem.MolFromSmarts(query_smarts)
v = Chem.MolFromSmiles("c1c[nH]nc1")
match = v.GetSubstructMatch(query)
new_mol = mol_builder(mol, v, {old:new for old, new in zip(index_to_be_replaced, match)})
new_mol

Here are some potential areas for improvement:

1. Enhancing the method for generating customized query smarts.
2. Improving the process of editing molecules, such as replacing atoms or bonds.
3. Finding an approach to sanitizing molecules (post-processing?)
   I am eager to hear your insights on this issue.