| Previous changeset 6:55ac04db36aa (2020-07-28) |
|
Commit message:
"planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/chemicaltoolbox/sucos commit 05dc325ce687441e5d3bdbdedcc0e3529cd5e070" |
|
modified:
sucos.py sucos_cluster.py sucos_max.py utils.py |
| b |
| diff -r 55ac04db36aa -r 9b48456a96fe sucos.py --- a/sucos.py Tue Jul 28 08:47:57 2020 -0400 +++ b/sucos.py Wed Apr 14 09:30:28 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 55ac04db36aa -r 9b48456a96fe sucos_cluster.py --- a/sucos_cluster.py Tue Jul 28 08:47:57 2020 -0400 +++ b/sucos_cluster.py Wed Apr 14 09:30:28 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 55ac04db36aa -r 9b48456a96fe sucos_max.py --- a/sucos_max.py Tue Jul 28 08:47:57 2020 -0400 +++ b/sucos_max.py Wed Apr 14 09:30:28 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 55ac04db36aa -r 9b48456a96fe utils.py --- a/utils.py Tue Jul 28 08:47:57 2020 -0400 +++ b/utils.py Wed Apr 14 09:30:28 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 |