Repository 'sucos_docking_scoring'
hg clone https://toolshed.g2.bx.psu.edu/repos/bgruening/sucos_docking_scoring

Changeset 6:4f1896782f7c (2021-04-14)
Previous changeset 5:fe318c648502 (2020-07-28)
Commit message:
"planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/chemicaltoolbox/sucos commit 05dc325ce687441e5d3bdbdedcc0e3529cd5e070"
modified:
sucos.py
sucos.xml
sucos_cluster.py
sucos_max.py
utils.py
b
diff -r fe318c648502 -r 4f1896782f7c sucos.py
--- a/sucos.py Tue Jul 28 08:48:35 2020 -0400
+++ b/sucos.py Wed Apr 14 09:31:11 2021 +0000
[
b'@@ -9,26 +9,39 @@\n """\n \n from __future__ import print_function\n-import argparse, os, sys, gzip\n+\n+import argparse\n+import os\n+\n import numpy as np\n-from rdkit import Chem, rdBase, RDConfig\n+import utils\n+from rdkit import Chem, RDConfig\n from rdkit.Chem import AllChem, rdShapeHelpers\n from rdkit.Chem.FeatMaps import FeatMaps\n-import utils\n \n-\n-### start function definitions #########################################\n+# start function definitions #########################################\n \n # Setting up the features to use in FeatureMap\n-fdef = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir, \'BaseFeatures.fdef\'))\n+fdef = AllChem.BuildFeatureFactory(\n+    os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")\n+)\n \n fmParams = {}\n for k in fdef.GetFeatureFamilies():\n     fparams = FeatMaps.FeatMapParams()\n     fmParams[k] = fparams\n \n-keep = (\'Donor\', \'Acceptor\', \'NegIonizable\', \'PosIonizable\', \'ZnBinder\',\n-        \'Aromatic\', \'Hydrophobe\', \'LumpedHydrophobe\')\n+keep = (\n+    "Donor",\n+    "Acceptor",\n+    "NegIonizable",\n+    "PosIonizable",\n+    "ZnBinder",\n+    "Aromatic",\n+    "Hydrophobe",\n+    "LumpedHydrophobe",\n+)\n+\n \n def filterFeature(f):\n     result = f.GetFamily() in keep\n@@ -37,6 +50,7 @@\n         utils.log("Filtered out feature type", f.GetFamily())\n     return result\n \n+\n def getRawFeatures(mol):\n \n     rawFeats = fdef.GetFeaturesForMol(mol)\n@@ -44,7 +58,10 @@\n     filtered = list(filter(filterFeature, rawFeats))\n     return filtered\n \n-def get_FeatureMapScore(small_feats, large_feats, tani=False, score_mode=FeatMaps.FeatMapScoreMode.All):\n+\n+def get_FeatureMapScore(\n+    small_feats, large_feats, tani=False, score_mode=FeatMaps.FeatMapScoreMode.All\n+):\n     """\n     Generate the feature map score.\n \n@@ -58,7 +75,10 @@\n     for rawFeats in [small_feats, large_feats]:\n         # filter that list down to only include the ones we\'re interested in\n         featLists.append(rawFeats)\n-    fms = [FeatMaps.FeatMap(feats=x, weights=[1] * len(x), params=fmParams) for x in featLists]\n+    fms = [\n+        FeatMaps.FeatMap(feats=x, weights=[1] * len(x), params=fmParams)\n+        for x in featLists\n+    ]\n     # set the score mode\n     fms[0].scoreMode = score_mode\n \n@@ -69,24 +89,35 @@\n             B = len(featLists[1])\n             if B != fms[1].GetNumFeatures():\n                 utils.log("Why isn\'t B equal to number of features...?!")\n-            tani_score = float(c) / (A+B-c)\n+            tani_score = float(c) / (A + B - c)\n             return tani_score\n         else:\n-            fm_score = fms[0].ScoreFeats(featLists[1]) / min(fms[0].GetNumFeatures(), len(featLists[1]))\n+            fm_score = fms[0].ScoreFeats(featLists[1]) / min(\n+                fms[0].GetNumFeatures(), len(featLists[1])\n+            )\n             return fm_score\n     except ZeroDivisionError:\n         utils.log("ZeroDivisionError")\n         return 0\n \n     if tani:\n-        tani_score = float(c) / (A+B-c)\n+        tani_score = float(c) / (A + B - c)\n         return tani_score\n     else:\n-        fm_score = fms[0].ScoreFeats(featLists[1]) / min(fms[0].GetNumFeatures(), len(featLists[1]))\n+        fm_score = fms[0].ScoreFeats(featLists[1]) / min(\n+            fms[0].GetNumFeatures(), len(featLists[1])\n+        )\n         return fm_score\n \n \n-def get_SucosScore(ref_mol, query_mol, tani=False, ref_features=None, query_features=None, score_mode=FeatMaps.FeatMapScoreMode.All):\n+def get_SucosScore(\n+    ref_mol,\n+    query_mol,\n+    tani=False,\n+    ref_features=None,\n+    query_features=None,\n+    score_mode=FeatMaps.FeatMapScoreMode.All,\n+):\n     """\n     This is the key function that calculates the SuCOS scores and is expected to be called from other modules.\n     To improve performance you can pre-calculate the features and pass them in as optional parameters to avoid having\n@@ -107,29 +138,41 @@\n         query_features = getRawFeatures(query_mol)\n \n     fm_score = get_FeatureMapScore(ref_features, query_features, tani, score_mode)\n-  '..b'core", fm_score)\n             if tani:\n@@ -155,9 +204,9 @@\n                 mol.SetDoubleProp("SuCOS_Protrude_Score", val3)\n             utils.log("Scores:", sucos_score, fm_score, val3)\n             writer.write(mol)\n-            total +=1\n+            total += 1\n         except ValueError as e:\n-            errors +=1\n+            errors += 1\n             utils.log("Molecule", count, "failed to score:", e.message)\n \n     input_file.close()\n@@ -165,41 +214,73 @@\n     writer.close()\n     output_file.close()\n \n-    utils.log("Completed.", total, "processed, ", count, "succeeded, ", errors, "errors")\n+    utils.log(\n+        "Completed.", total, "processed, ", count, "succeeded, ", errors, "errors"\n+    )\n+\n \n def parse_score_mode(value):\n-    if value == None or value == \'all\':\n+    if value is None or value == "all":\n         return FeatMaps.FeatMapScoreMode.All\n-    elif value == \'closest\':\n+    elif value == "closest":\n         return FeatMaps.FeatMapScoreMode.Closest\n-    elif value == \'best\':\n+    elif value == "best":\n         return FeatMaps.FeatMapScoreMode.Best\n     else:\n         raise ValueError(value + " is not a valid scoring mode option")\n \n \n-### start main execution #########################################\n+# start main execution #########################################\n+\n \n def main():\n \n-    parser = argparse.ArgumentParser(description=\'SuCOS with RDKit\')\n-    parser.add_argument(\'-i\', \'--input\', help=\'Input file in SDF format. Can be gzipped (*.gz).\')\n-    parser.add_argument(\'-r\', \'--refmol\', help=\'Molecule to compare against in Molfile (.mol) or SDF (.sdf) format\')\n-    parser.add_argument(\'--refmol-format\', help="Format for the reference molecule (mol or sdf). " +\n-                                                "Only needed if files don\'t have the expected extensions")\n-    parser.add_argument(\'--refmolidx\', help=\'Reference molecule index in SD file if not the first\', type=int, default=1)\n-    parser.add_argument(\'-o\', \'--output\', help=\'Output file in SDF format. Can be gzipped (*.gz).\')\n-    parser.add_argument(\'--tanimoto\', action=\'store_true\', help=\'Include Tanimoto distance in score\')\n-    parser.add_argument(\'--score_mode\', choices=[\'all\', \'closest\', \'best\'],\n-                        help="choose the scoring mode for the feature map, default is \'all\'.")\n+    parser = argparse.ArgumentParser(description="SuCOS with RDKit")\n+    parser.add_argument(\n+        "-i", "--input", help="Input file in SDF format. Can be gzipped (*.gz)."\n+    )\n+    parser.add_argument(\n+        "-r",\n+        "--refmol",\n+        help="Molecule to compare against in Molfile (.mol) or SDF (.sdf) format",\n+    )\n+    parser.add_argument(\n+        "--refmol-format",\n+        help="Format for the reference molecule (mol or sdf). "\n+        + "Only needed if files don\'t have the expected extensions",\n+    )\n+    parser.add_argument(\n+        "--refmolidx",\n+        help="Reference molecule index in SD file if not the first",\n+        type=int,\n+        default=1,\n+    )\n+    parser.add_argument(\n+        "-o", "--output", help="Output file in SDF format. Can be gzipped (*.gz)."\n+    )\n+    parser.add_argument(\n+        "--tanimoto", action="store_true", help="Include Tanimoto distance in score"\n+    )\n+    parser.add_argument(\n+        "--score_mode",\n+        choices=["all", "closest", "best"],\n+        help="choose the scoring mode for the feature map, default is \'all\'.",\n+    )\n \n     args = parser.parse_args()\n     utils.log("SuCOS Args: ", args)\n \n     score_mode = parse_score_mode(args.score_mode)\n \n-    process(args.refmol, args.input, args.output, refmol_index=args.refmolidx,\n-            refmol_format=args.refmol_format, tani=args.tanimoto, score_mode=score_mode)\n+    process(\n+        args.refmol,\n+        args.input,\n+        args.output,\n+        refmol_index=args.refmolidx,\n+        refmol_format=args.refmol_format,\n+        tani=args.tanimoto,\n+        score_mode=score_mode,\n+    )\n \n \n if __name__ == "__main__":\n'
b
diff -r fe318c648502 -r 4f1896782f7c sucos.xml
--- a/sucos.xml Tue Jul 28 08:48:35 2020 -0400
+++ b/sucos.xml Wed Apr 14 09:31:11 2021 +0000
[
@@ -2,7 +2,7 @@
     <description>- compare shape and feature overlap of docked ligand poses to a reference molecule</description>
     <macros>
         <import>sucos_macros.xml</import>
-        <token name="@GALAXY_VERSION@">0</token>
+        <token name="@GALAXY_VERSION@">1</token>
     </macros>
     <expand macro="requirements"/>
     <command detect_errors="exit_code"><![CDATA[
b
diff -r fe318c648502 -r 4f1896782f7c sucos_cluster.py
--- a/sucos_cluster.py Tue Jul 28 08:48:35 2020 -0400
+++ b/sucos_cluster.py Wed Apr 14 09:31:11 2021 +0000
[
@@ -10,15 +10,16 @@
 GitHub: https://github.com/susanhleung/SuCOS
 Publication: https://doi.org/10.26434/chemrxiv.8100203.v1
 """
+import argparse
 
-import sucos, utils
-import argparse, gzip
-from rdkit import Chem
 import numpy as np
 import pandas as pd
-from scipy.cluster.hierarchy import linkage, fcluster
+import sucos
+import utils
+from rdkit import Chem
+from scipy.cluster.hierarchy import fcluster, linkage
 
-### start main execution #########################################
+# start main execution #########################################
 
 
 def calc_distance_matrix(mols):
@@ -44,13 +45,17 @@
             if tuple1[0] == tuple2[0]:
                 tmp.append(0.0)
             else:
-                #utils.log("Calculating SuCOS between", mol1, mol2)
-                sucos_score, fm_score, tani_score = sucos.get_SucosScore(tuple1[0], tuple2[0],
-                    tani=True, ref_features=tuple1[1], query_features=tuple2[1])
+                # utils.log("Calculating SuCOS between", mol1, mol2)
+                sucos_score, fm_score, tani_score = sucos.get_SucosScore(
+                    tuple1[0],
+                    tuple2[0],
+                    tani=True,
+                    ref_features=tuple1[1],
+                    query_features=tuple2[1],
+                )
                 tmp.append(1.0 - sucos_score)
         matrix.append(tmp)
 
-
     return matrix
 
 
@@ -64,24 +69,25 @@
 
     indexes = [x for x in range(0, len(matrix))]
     cols = [x for x in range(0, len(matrix[0]))]
-    #utils.log("indexes", indexes)
-    #utils.log("cols", cols)
+    # utils.log("indexes", indexes)
+    # utils.log("cols", cols)
     df = pd.DataFrame(matrix, columns=cols, index=indexes)
     utils.log("DataFrame:", df.shape)
-    #utils.log(df)
+    # utils.log(df)
     indices = np.triu_indices(df.shape[0], k=1)
-    #utils.log("Indices:", indices)
+    # utils.log("Indices:", indices)
     t = np.array(df)[indices]
-    Z = linkage(t, 'average')
+    Z = linkage(t, "average")
     lig_clusters = []
-    cluster_arr = fcluster(Z, t=threshold, criterion='distance')
+    cluster_arr = fcluster(Z, t=threshold, criterion="distance")
     for i in range(np.amax(cluster_arr)):
-        clus = df.columns[np.argwhere(cluster_arr==i+1)]
+        clus = df.columns[np.argwhere(cluster_arr == i + 1)]
         lig_clusters.append([x[0] for x in clus.tolist()])
 
     utils.log("Clusters", lig_clusters)
     return lig_clusters
 
+
 def write_clusters_to_sdfs(mols, clusters, basename, gzip=False):
     """
     Write the molecules to SDF files, 1 file for each cluster.
@@ -99,7 +105,9 @@
         filename = basename + str(i) + ".sdf"
         if gzip:
             filename += ".gz"
-        utils.log("Writing ", len(cluster), "molecules in cluster", i, "to file", filename)
+        utils.log(
+            "Writing ", len(cluster), "molecules in cluster", i, "to file", filename
+        )
         output_file = utils.open_file_for_writing(filename)
         writer = Chem.SDWriter(output_file)
         for index in cluster:
@@ -110,14 +118,26 @@
         output_file.close()
 
 
-
 def main():
-    parser = argparse.ArgumentParser(description='Clustering with SuCOS and RDKit')
-    parser.add_argument('-i', '--input', help='Input file in SDF format. Can be gzipped (*.gz).')
-    parser.add_argument('-o', '--output', default="cluster", help="Base name for output files in SDF format. " +
-                                               "e.g. if value is 'output' then files like output1.sdf, output2.sdf will be created")
-    parser.add_argument('--gzip', action='store_true', help='Gzip the outputs generating files like output1.sdf.gz, output2.sdf.gz')
-    parser.add_argument('-t', '--threshold', type=float, default=0.8, help='Clustering threshold')
+    parser = argparse.ArgumentParser(description="Clustering with SuCOS and RDKit")
+    parser.add_argument(
+        "-i", "--input", help="Input file in SDF format. Can be gzipped (*.gz)."
+    )
+    parser.add_argument(
+        "-o",
+        "--output",
+        default="cluster",
+        help="Base name for output files in SDF format. "
+        + "e.g. if value is 'output' then files like output1.sdf, output2.sdf will be created",
+    )
+    parser.add_argument(
+        "--gzip",
+        action="store_true",
+        help="Gzip the outputs generating files like output1.sdf.gz, output2.sdf.gz",
+    )
+    parser.add_argument(
+        "-t", "--threshold", type=float, default=0.8, help="Clustering threshold"
+    )
 
     args = parser.parse_args()
     utils.log("SuCOS Cluster Args: ", args)
@@ -131,4 +151,4 @@
 
 
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
b
diff -r fe318c648502 -r 4f1896782f7c sucos_max.py
--- a/sucos_max.py Tue Jul 28 08:48:35 2020 -0400
+++ b/sucos_max.py Wed Apr 14 09:31:11 2021 +0000
[
@@ -34,12 +34,17 @@
 Publication: https://doi.org/10.26434/chemrxiv.8100203.v1
 """
 
-import sucos, utils
-import argparse, gzip, os
+import argparse
+import os
+
+import sucos
+import utils
 from rdkit import Chem
 
 
-def process(inputfilename, clusterfilenames, outputfilename, filter_value, filter_field):
+def process(
+    inputfilename, clusterfilenames, outputfilename, filter_value, filter_field
+):
     all_clusters = {}
     for filename in clusterfilenames:
         cluster = []
@@ -49,13 +54,20 @@
         for mol in suppl:
             i += 1
             if not mol:
-                utils.log("WARNING: failed to generate molecule", i, "in cluster", filename)
+                utils.log(
+                    "WARNING: failed to generate molecule", i, "in cluster", filename
+                )
                 continue
             try:
                 features = sucos.getRawFeatures(mol)
                 cluster.append((mol, features))
-            except:
-                utils.log("WARNING: failed to generate features for molecule", i, "in cluster", filename)
+            except Exception:
+                utils.log(
+                    "WARNING: failed to generate features for molecule",
+                    i,
+                    "in cluster",
+                    filename,
+                )
 
         cluster_file.close()
         all_clusters[filename] = cluster
@@ -75,8 +87,10 @@
             continue
         try:
             query_features = sucos.getRawFeatures(mol)
-        except:
-            utils.log("WARNING: failed to generate features for molecule", mol_num, "in input")
+        except Exception:
+            utils.log(
+                "WARNING: failed to generate features for molecule", mol_num, "in input"
+            )
             continue
         scores_max = [0, 0, 0]
         scores_cum = [0, 0, 0]
@@ -89,9 +103,13 @@
                 ref_features = entry[1]
                 index += 1
                 comparisons += 1
-                sucos_score, fm_score, vol_score = sucos.get_SucosScore(hit, mol,
-                                                                        tani=False, ref_features=ref_features,
-                                                                        query_features=query_features)
+                sucos_score, fm_score, vol_score = sucos.get_SucosScore(
+                    hit,
+                    mol,
+                    tani=False,
+                    ref_features=ref_features,
+                    query_features=query_features,
+                )
 
                 if sucos_score > scores_max[0]:
                     scores_max[0] = sucos_score
@@ -104,11 +122,14 @@
                 scores_cum[1] += fm_score
                 scores_cum[2] += vol_score
 
-
         # utils.log("Max SuCOS:", scores[0], "FM:", scores[1], "P:", scores[2],"File:", cluster_file_name_only, "Index:", cluster_index)
         mol.SetDoubleProp("Max_SuCOS_Score", scores_max[0] if scores_max[0] > 0 else 0)
-        mol.SetDoubleProp("Max_SuCOS_FeatureMap_Score", scores_max[1] if scores_max[1] > 0 else 0)
-        mol.SetDoubleProp("Max_SuCOS_Protrude_Score", scores_max[2] if scores_max[2] > 0 else 0)
+        mol.SetDoubleProp(
+            "Max_SuCOS_FeatureMap_Score", scores_max[1] if scores_max[1] > 0 else 0
+        )
+        mol.SetDoubleProp(
+            "Max_SuCOS_Protrude_Score", scores_max[2] if scores_max[2] > 0 else 0
+        )
 
         if cluster_name:
             cluster_file_name_only = cluster_name.split(os.sep)[-1]
@@ -117,8 +138,12 @@
 
         # utils.log("Cum SuCOS:", scores[0], "FM:", scores[1], "P:", scores[2])
         mol.SetDoubleProp("Cum_SuCOS_Score", scores_cum[0] if scores_cum[0] > 0 else 0)
-        mol.SetDoubleProp("Cum_SuCOS_FeatureMap_Score", scores_cum[1] if scores_cum[1] > 0 else 0)
-        mol.SetDoubleProp("Cum_SuCOS_Protrude_Score", scores_cum[2] if scores_cum[2] > 0 else 0)
+        mol.SetDoubleProp(
+            "Cum_SuCOS_FeatureMap_Score", scores_cum[1] if scores_cum[1] > 0 else 0
+        )
+        mol.SetDoubleProp(
+            "Cum_SuCOS_Protrude_Score", scores_cum[2] if scores_cum[2] > 0 else 0
+        )
 
         if filter_value and filter_field:
             if mol.HasProp(filter_field):
@@ -136,20 +161,35 @@
     utils.log("Completed", comparisons, "comparisons")
 
 
-### start main execution #########################################
+# start main execution #########################################
+
 
 def main():
-    parser = argparse.ArgumentParser(description='Max SuCOS scores with RDKit')
-    parser.add_argument('-i', '--input', help='Input file to score in SDF format. Can be gzipped (*.gz).')
-    parser.add_argument('-o', '--output', help='Output file in SDF format. Can be gzipped (*.gz).')
-    parser.add_argument('clusters', nargs='*', help="One or more SDF files with the clustered hits")
-    parser.add_argument('--filter-value', type=float, help='Filter out values with scores less than this.')
-    parser.add_argument('--filter-field', help='Field to use to filter values.')
+    parser = argparse.ArgumentParser(description="Max SuCOS scores with RDKit")
+    parser.add_argument(
+        "-i",
+        "--input",
+        help="Input file to score in SDF format. Can be gzipped (*.gz).",
+    )
+    parser.add_argument(
+        "-o", "--output", help="Output file in SDF format. Can be gzipped (*.gz)."
+    )
+    parser.add_argument(
+        "clusters", nargs="*", help="One or more SDF files with the clustered hits"
+    )
+    parser.add_argument(
+        "--filter-value",
+        type=float,
+        help="Filter out values with scores less than this.",
+    )
+    parser.add_argument("--filter-field", help="Field to use to filter values.")
 
     args = parser.parse_args()
     utils.log("Max SuCOS Args: ", args)
 
-    process(args.input, args.clusters, args.output, args.filter_value, args.filter_field)
+    process(
+        args.input, args.clusters, args.output, args.filter_value, args.filter_field
+    )
 
 
 if __name__ == "__main__":
b
diff -r fe318c648502 -r 4f1896782f7c utils.py
--- a/utils.py Tue Jul 28 08:48:35 2020 -0400
+++ b/utils.py Wed Apr 14 09:31:11 2021 +0000
b
@@ -4,40 +4,54 @@
 """
 
 from __future__ import print_function
-import sys, gzip
+
+import gzip
+import sys
+
 from rdkit import Chem
 
+
 def log(*args, **kwargs):
-    """Log output to STDERR
-    """
+    """Log output to STDERR"""
     print(*args, file=sys.stderr, **kwargs)
 
+
 def open_file_for_reading(filename):
     """Open the file gunzipping it if it ends with .gz."""
-    if filename.lower().endswith('.gz'):
-        return gzip.open(filename, 'rb')
+    if filename.lower().endswith(".gz"):
+        return gzip.open(filename, "rb")
     else:
-        return open(filename, 'rb')
+        return open(filename, "rb")
+
 
 def open_file_for_writing(filename):
-    if filename.lower().endswith('.gz'):
-        return gzip.open(filename, 'at')
+    if filename.lower().endswith(".gz"):
+        return gzip.open(filename, "at")
     else:
-        return open(filename, 'w+')
+        return open(filename, "w+")
+
 
 def read_single_molecule(filename, index=1, format=None):
     """Read a single molecule as a RDKit Mol object. This can come from a file in molfile or SDF format.
     If SDF then you can also specify an index of the molecule that is read (default is the first)
     """
     mol = None
-    if format == 'mol' or filename.lower().endswith('.mol') or filename.lower().endswith('.mol.gz'):
+    if (
+        format == "mol"
+        or filename.lower().endswith(".mol")
+        or filename.lower().endswith(".mol.gz")
+    ):
         file = open_file_for_reading(filename)
         mol = Chem.MolFromMolBlock(file.read())
         file.close()
-    elif format == 'sdf' or filename.lower().endswith('.sdf') or filename.lower().endswith('.sdf.gz'):
+    elif (
+        format == "sdf"
+        or filename.lower().endswith(".sdf")
+        or filename.lower().endswith(".sdf.gz")
+    ):
         file = open_file_for_reading(filename)
         supplier = Chem.ForwardSDMolSupplier(file)
-        for i in range(0,index):
+        for i in range(0, index):
             if supplier.atEnd():
                 break
             mol = next(supplier)
@@ -46,4 +60,4 @@
     if not mol:
         raise ValueError("Unable to read molecule")
 
-    return mol
\ No newline at end of file
+    return mol