alignmcs

alignmcs.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+import argparse, sys, rdkit, collections
+from rdkit.Chem import AllChem as Chem
+from rdkit.Chem import rdMolAlign
+from rdkit.Chem.rdShape import Align
+from rdkit.Chem.rdMolAlign import AlignMol
+from rdkit.Chem import rdMolAlign
+from rdkit.Chem import rdFMCS
+from rdkit.Chem import rdMolTransforms
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Align query molecules to a reference molecule using their maximum common substructure. May produce multiple alignments per an input.')
+    parser.add_argument('ref',help='Reference SDF file for comparison')
+    parser.add_argument('query',help="SDF file to compare with")
+    parser.add_argument('out',help='Output SDF')
+    parser.add_argument('-s','--strict',action='store_true',help='Use strict atom/bond matching')
+    args = parser.parse_args()
+
+    #load reference file molecules
+    refmols = [mol for mol in Chem.SDMolSupplier(args.ref)]
+
+    if len(refmols) > 1:
+        print("Only using first reference molecule")
+    refmol = refmols[0]
+
+    out = Chem.SDWriter(args.out)
+    for mol in Chem.SDMolSupplier(args.query):
+        if args.strict:
+            mcs = rdFMCS.FindMCS([refmol,mol])
+        else:
+            mcs = rdFMCS.FindMCS([refmol,mol],atomCompare=rdFMCS.AtomCompare.CompareAny,bondCompare=rdFMCS.BondCompare.CompareAny)
+
+        submol = Chem.MolFromSmarts(mcs.smartsString)
+
+        for refmatch in refmol.GetSubstructMatches(submol):
+            for qmatch in mol.GetSubstructMatches(submol):
+                AlignMol(mol, refmol, atomMap=list(zip(qmatch,refmatch)))
+                out.write(mol)

rdconf.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/python3
 
 import sys,string,argparse
 from rdkit.Chem import AllChem as Chem
@@ -49,10 +49,10 @@ def getRMS(mol, c1,c2):
 output = args[1]
 smifile = open(input)
 if options.verbose:
-    print "Generating a maximum of",options.maxconfs,"per a mol"
+    print("Generating a maximum of",options.maxconfs,"per a mol")
 
 if options.etkdg and not Chem.ETKDG:
-    print "ETKDB does not appear to be implemented.  Please upgrade RDKit."
+    print("ETKDB does not appear to be implemented.  Please upgrade RDKit.")
     sys.exit(1)
 
 split = os.path.splitext(output)
@@ -69,23 +69,23 @@ def getRMS(mol, c1,c2):
     sdwriter = Chem.SDWriter(outf)
 
 if sdwriter is None:
-    print "Could not open ".output
+    print("Could not open ".output)
     sys.exit(-1)
 
 for line in smifile:
     toks = line.split()
     smi = toks[0]
-    name = string.join(toks[1:])    
+    name = ' '.join(toks[1:])    
 
     pieces = smi.split('.')
     if len(pieces) > 1:
         smi = max(pieces, key=len) #take largest component by length
-        print "Taking largest component: %s\t%s" % (smi,name)
+        print("Taking largest component: %s\t%s" % (smi,name))
 
     mol = Chem.MolFromSmiles(smi)
     if mol is not None:
         if options.verbose:
-            print smi
+            print(smi)
         try:
             Chem.SanitizeMol(mol)
             mol = Chem.AddHs(mol)
@@ -96,7 +96,7 @@ def getRMS(mol, c1,c2):
             else:
                 cids = Chem.EmbedMultipleConfs(mol, int(options.sample*options.maxconfs),randomSeed=options.seed)
             if options.verbose:
-                print len(cids),"conformers found"
+                print(len(cids),"conformers found")
             cenergy = []            
             for conf in cids:
                 #not passing confID only minimizes the first conformer
@@ -110,7 +110,7 @@ def getRMS(mol, c1,c2):
                     converged = not Chem.UFFOptimizeMolecule(mol,confId=conf)
                     cenergy.append(Chem.UFFGetMoleculeForceField(mol,confId=conf).CalcEnergy())
                 if options.verbose:
-                    print "Convergence of conformer",conf,converged
+                    print("Convergence of conformer",conf,converged)
 
             mol = Chem.RemoveHs(mol)
             sortedcids = sorted(cids,key = lambda cid: cenergy[cid])
@@ -133,7 +133,7 @@ def getRMS(mol, c1,c2):
                         break
                     #check rmsd
                     passed = True
-                    for seenconf in written.iterkeys():
+                    for seenconf in written.keys():
                         rms = getRMS(mol,seenconf,conf) 
                         if(rms < options.rms) or (options.energy > 0 and cenergy[conf]-mine > options.energy):
                             passed = False
@@ -144,9 +144,9 @@ def getRMS(mol, c1,c2):
         except (KeyboardInterrupt, SystemExit):
             raise                
         except Exception as e:
-            print "Exception",e
+            print("Exception",e)
     else:
-        print "ERROR:",smi
+        print("ERROR:",smi)
 
 sdwriter.close()
 outf.close()