Dear Jerry, When calculating the RMS, try:
rms = AllChem.CalcRMS(probeMol, refMol, map=[list(zip(mv, refMatch))]) Note that you are calculating MCS across all mols, however; you might want to consider pairwise-MCS (i.e. MCS between each docked ligand and the ref ligand individually and then calculate RMS based on it) since you want difference between each docked ligand and reference ligand. I hope this works for you. Best regards, Omar On Sat, Dec 19, 2020 at 2:03 AM Parks, Jerry M. via Rdkit-discuss < rdkit-discuss@lists.sourceforge.net> wrote: > Hello > > > > I’m trying to calculate RMSDs between a reference molecule from an X-ray > structure and a series of similar but somewhat different molecules with a > common scaffold that I have docked to the protein. I don’t want to align > the structures, I just want to measure how different their docked poses are > to the reference structure. All of the molecules have the same core > structure, but have different functional groups. This functionality is > available in the recently developed web server LigRMSD ( > https://ligrmsd.appsbio.utalca.cl), but I need to make some modifications > to the procedure. > > > > Following this thread ( > https://sourceforge.net/p/rdkit/mailman/message/35678724/) and this blog > post (http://rdkit.blogspot.com/2013/12/using-allchemconstrainedembed.html), > I’ve written the following code: > > > > from rdkit import Chem > > from rdkit.Chem import AllChem > > from rdkit.Chem.Draw import IPythonConsole > > from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions > > from rdkit.Chem import rdFMCS > > import glob > > > > # Load a set of mol files. > > mols = [] > > for m in glob.glob('*.mol'): > > tmpmol = Chem.MolFromMolFile(m) > > name_without_suffix = m[:-4] > > tmpmol.SetProp('_Name', name_without_suffix) > > mols.append(tmpmol) > > > > # Draw the molecules. > > Chem.Draw.MolsToGridImage(mols, molsPerRow=2, subImgSize=(350, 350), > legends=[mol.GetProp('_Name') for mol in mols]) > > > > # Find the MCS. > > mcs = rdFMCS.FindMCS(mols,threshold=1.0, ringMatchesRingOnly=True, > completeRingsOnly=True) > > > > # Align everything to the reference molecule. > > > > patt = Chem.MolFromSmarts(mcs.smartsString) > > refMol = mols[0] > > refMatch = refMol.GetSubstructMatch(patt) > > rmsVs = [] > > for probeMol in mols[1:]: > > mv = probeMol.GetSubstructMatch(patt) > > rms = AllChem.AlignMol(probeMol, refMol, atomMap=list(zip(mv, > refMatch)), maxIters=0) > > print(probeMol.GetProp('_Name'), f'{rms:.2f}') > > rmsVs.append(rms) > > > > However, even though I set maxIters=0 in the call to AlignMol, the > molecules have all been translated to a common origin (which is not what I > want). > > > > I tried replacing that line of code with the following call to CalcRMS > (using the same atom mapping): > > > > rms = AllChem.CalcRMS(probeMol, refMol, map=list(zip(mv, refMatch))) > > > > but then I got the following error: > > > > --------------------------------------------------------------------------- > > ValueError Traceback (most recent call > last) > > <ipython-input-7-0faef2f852d9> in <module>() > > * 8* mv = probeMol.GetSubstructMatch(patt) > > * 9* #rms = AllChem.AlignMol(probeMol, refMol, > atomMap=list(zip(mv, refMatch)), maxIters=0) > > ---> 10 rms = AllChem.CalcRMS(probeMol, refMol, map=list(zip(mv, > refMatch))) > > * 11* print(probeMol.GetProp('_Name'), f'{rms:.2f}') > > * 12* rmsVs.append(rms) > > > > ValueError: sequence does not support length query > > > > Any suggestions for fixing the code would be appreciated. > > > > Thank you, > > Jerry Parks > _______________________________________________ > Rdkit-discuss mailing list > Rdkit-discuss@lists.sourceforge.net > https://lists.sourceforge.net/lists/listinfo/rdkit-discuss >
_______________________________________________ Rdkit-discuss mailing list Rdkit-discuss@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
