fixed chuk size for prepare_tether.py

Author: tdudgeon

src/python/pipelines/xchem/prepare_tether.py

@@ -59,9 +59,9 @@ def get_writer(outfile_base):
     global file_count
 
     if not writer:
-        writer = Chem.SDWriter(outfile_base + '_' + f'{file_count:04}' + '.sdf')
+        writer = Chem.SDWriter(outfile_base + '_' + f'{file_count:03}' + '.sdf')
 
-    if write_count < chunk_size:
+    if write_count <= chunk_size:
         # print('  Using existing writer for', id, writer)
         return writer
 
@@ -73,7 +73,7 @@ def get_writer(outfile_base):
     writer.close()
     file_count += 1
 
-    writer = Chem.SDWriter(outfile_base + '_' + f'{file_count:04}' + '.sdf')
+    writer = Chem.SDWriter(outfile_base + '_' + f'{file_count:03}' + '.sdf')
     #print('  Using new writer for', id, writer)
     return writer
 
@@ -361,11 +361,89 @@ def minimize_mol(mol, target, molMatch, targetMatch, algMap, getForceField):
     rms = AlignMol(mol, target, atomMap=algMap)
     mol.SetDoubleProp('EmbedRMS', rms)
 
+def find_best_mcs(hit, mol):
+    best_score = 0
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, ringMatchesRingOnly=True,
+                          atomCompare=rdFMCS.AtomCompare.CompareElements, bondCompare=rdFMCS.BondCompare.CompareOrder)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 1
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, ringMatchesRingOnly=False,
+                          atomCompare=rdFMCS.AtomCompare.CompareElements, bondCompare=rdFMCS.BondCompare.CompareOrder)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 2
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=True,
+                          atomCompare=rdFMCS.AtomCompare.CompareElements, bondCompare=rdFMCS.BondCompare.CompareOrder)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 3
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=False,
+                          atomCompare=rdFMCS.AtomCompare.CompareElements, bondCompare=rdFMCS.BondCompare.CompareOrder)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 4
+
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, ringMatchesRingOnly=True,
+                          atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 5
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, ringMatchesRingOnly=False,
+                          atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 6
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=True,
+                          atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 7
+    mcs = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=False,
+                          atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
+    score = score_mcs(hit, mol, mcs)
+    if score > best_score:
+        best_score = score
+        best_mcs = mcs
+        best_index = 8
+    print("  Best MCS is", best_index, 'with score', best_score)
+    return best_mcs
+
+def score_mcs(hit, mol, mcs):
+    smartsMol = Chem.MolFromSmarts(mcs.smartsString)
+    hitMatches = hit.GetSubstructMatch(smartsMol)
+    molMatches = mol.GetSubstructMatch(smartsMol)
+    # print('  Matches:', hitMatches, molMatches)
+    score = 0
+    for hitMatch, molMatch in zip(hitMatches, molMatches):
+        score += 1
+        if hit.GetAtomWithIdx(hitMatch).GetAtomicNum() == mol.GetAtomWithIdx(molMatch).GetAtomicNum():
+            score += 1
+    if mol.GetAtomWithIdx(molMatch).IsInRing():
+        score += 1
+    print('  Score:', score)
+    return score
 
 def execute(smi, hit_molfile, outfile_base, min_ph=None, max_ph=None, max_inputs=0, max_outputs=0, modulus=0, timout_embed_secs=5):
 
     global write_count
 
+    GetFF=lambda x,confId=-1:AllChem.MMFFGetMoleculeForceField(x,AllChem.MMFFGetMoleculeProperties(x),confId=confId)
+
     hit = Chem.MolFromMolFile(hit_molfile)
 
     num_mols = 0
@@ -395,6 +473,8 @@ def execute(smi, hit_molfile, outfile_base, min_ph=None, max_ph=None, max_inputs
             num_processed += 1
             num_added = 0
 
+            w = get_writer(outfile_base)
+
             print('Processing', num_mols, num_processed, smiles, Chem.MolToSmiles(hit))
 
             try:
@@ -406,14 +486,19 @@ def execute(smi, hit_molfile, outfile_base, min_ph=None, max_ph=None, max_inputs
                 else:
                     enumerated_mols = [mol]
 
-                # mcs0 = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, matchValences=False, ringMatchesRingOnly=True,
-                #                       atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
-                mcs0 = rdFMCS.FindMCS([hit, mol], completeRingsOnly=True, ringMatchesRingOnly=True)
+                mcs0 = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=False,
+                                      atomCompare=rdFMCS.AtomCompare.CompareAny, bondCompare=rdFMCS.BondCompare.CompareAny)
+                # mcs0 = rdFMCS.FindMCS([hit, mol], completeRingsOnly=False, ringMatchesRingOnly=False)
+                # score_mcs(hit, mol, mcs0)
+
+                # mcs0 = find_best_mcs(hit, mol)
+
                 mcsQuery = Chem.MolFromSmarts(mcs0.smartsString)
 
                 count = 0
                 for ligand in enumerated_mols:
                     molh = Chem.AddHs(ligand)
+                    # mol_match_tuple = MultiConstrainedEmbed(molh, queryMol, getForceField=GetFF)
                     mol_match_tuple = multi_constrained_embed(molh, hit, mcsQuery, timout_embed_secs=timout_embed_secs)
                     print(' ', len(mol_match_tuple), 'mols tethered')
 
@@ -429,10 +514,8 @@ def execute(smi, hit_molfile, outfile_base, min_ph=None, max_ph=None, max_inputs
                         t_mol.SetProp('TETHERED ATOMS', tethers)
                         # print(' Tethers: ', tethers)
 
-                        w = get_writer(outfile_base)
                         w.write(t_mol)
                         write_count += 1
-                        print('  Write count =', write_count)
                         num_outputs += 1
                         num_added += 1
 
@@ -455,7 +538,7 @@ def execute(smi, hit_molfile, outfile_base, min_ph=None, max_ph=None, max_inputs
 def main():
     """
     Example usage:
-    python -m pipelines.xchem.prepare_tether --smi ../../data/mpro/Mpro-x0387_0.smi --mol ../../data/mpro/Mpro-x0387_0.mol -o TETHERED --max-inputs 500 --chunk-size 100
+    python -m pipelines.xchem.prepare_tether_2 --smi ../../data/mpro/Mpro-x0387_0.smi --mol ../../data/mpro/Mpro-x0387_0.mol -o TETHERED --max-inputs 500 --chunk-size 100
 
     :return:
     """