|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +# Copyright 2018 Informatics Matters Ltd. |
| 4 | +# |
| 5 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 6 | +# you may not use this file except in compliance with the License. |
| 7 | +# You may obtain a copy of the License at |
| 8 | +# |
| 9 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 10 | +# |
| 11 | +# Unless required by applicable law or agreed to in writing, software |
| 12 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | +# See the License for the specific language governing permissions and |
| 15 | +# limitations under the License. |
| 16 | + |
| 17 | +import argparse, os |
| 18 | + |
| 19 | +from rdkit import Chem, rdBase, RDConfig |
| 20 | +from rdkit.Chem import AllChem, rdShapeHelpers |
| 21 | +from rdkit.Chem.FeatMaps import FeatMaps |
| 22 | + |
| 23 | +from pipelines_utils import parameter_utils, utils |
| 24 | +from pipelines_utils_rdkit import rdkit_utils |
| 25 | + |
| 26 | +### start field name defintions ######################################### |
| 27 | + |
| 28 | +field_SuCOS_Score = "SuCOS_Score" |
| 29 | + |
| 30 | +### start function definitions ######################################### |
| 31 | + |
| 32 | +################################################# |
| 33 | +#### Setting up the features to use in FeatureMap |
| 34 | +fdef = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir, 'BaseFeatures.fdef')) |
| 35 | +# keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic') |
| 36 | + |
| 37 | +fmParams = {} |
| 38 | +for k in fdef.GetFeatureFamilies(): |
| 39 | + fparams = FeatMaps.FeatMapParams() |
| 40 | + fmParams[k] = fparams |
| 41 | + |
| 42 | +keep = ('Donor', 'Acceptor', 'NegIonizable', 'PosIonizable', 'ZnBinder', |
| 43 | + 'Aromatic', 'Hydrophobe', 'LumpedHydrophobe') |
| 44 | + |
| 45 | +def get_FeatureMapScore(ref_mol, query_mol): |
| 46 | + featLists = [] |
| 47 | + for m in [ref_mol, query_mol]: |
| 48 | + rawFeats = fdef.GetFeaturesForMol(m) |
| 49 | + # filter that list down to only include the ones we're interested in |
| 50 | + featLists.append([f for f in rawFeats if f.GetFamily() in keep]) |
| 51 | + fms = [FeatMaps.FeatMap(feats=x, weights=[1] * len(x), params=fmParams) for x in featLists] |
| 52 | + #utils.log("Calc:", str(fms[0].ScoreFeats(featLists[1])), "/ float(min(", str(fms[0].GetNumFeatures()), str(len(featLists[1])), "))") |
| 53 | + fm_score = fms[0].ScoreFeats(featLists[1]) / float(min(fms[0].GetNumFeatures(), len(featLists[1]))) |
| 54 | + return fm_score |
| 55 | + |
| 56 | +def get_SucosScore(ref_mol, query_mol, field_name): |
| 57 | + fm_score = get_FeatureMapScore(ref_mol, query_mol) |
| 58 | + #utils.log("FeatureMapScore:", str(fm_score)) |
| 59 | + protrude_dist = rdShapeHelpers.ShapeProtrudeDist(ref_mol, query_mol, allowReordering=False) |
| 60 | + #utils.log("ProtrudeDistance:", str(protrude_dist)) |
| 61 | + #utils.log("Sucos calc: 0.5 *", str(fm_score), "+ 0.5 * (1.0 -", protrude_dist, ")") |
| 62 | + score = 0.5 * fm_score + 0.5 * (1.0 - protrude_dist) |
| 63 | + #utils.log("SucosScore:", str(score)) |
| 64 | + query_mol.SetDoubleProp(field_name, score) |
| 65 | + return score |
| 66 | + |
| 67 | +### start main execution ######################################### |
| 68 | + |
| 69 | +def main(): |
| 70 | + |
| 71 | + parser = argparse.ArgumentParser(description='SuCOS with RDKit') |
| 72 | + parser.add_argument('--target', help='molecule to compare against') |
| 73 | + parser.add_argument('--targetidx', help="Target molecule index in SD file if not the first", type=int, default=1) |
| 74 | + parameter_utils.add_default_io_args(parser) |
| 75 | + |
| 76 | + args = parser.parse_args() |
| 77 | + utils.log("SuCOS Args: ", args) |
| 78 | + |
| 79 | + # TODO - handle molecules with multiple fragments |
| 80 | + |
| 81 | + ref_mol = rdkit_utils.read_single_molecule(args.target, index=args.targetidx) |
| 82 | + utils.log("Reference mol has", str(ref_mol.GetNumHeavyAtoms()), "heavy atoms") |
| 83 | + |
| 84 | + source = "sucos.py" |
| 85 | + datasetMetaProps = {"source":source, "description": "SuCOS using RDKit " + rdBase.rdkitVersion} |
| 86 | + clsMappings = { "SuCOS_score": "java.lang.Float" } |
| 87 | + fieldMetaProps = [ |
| 88 | + {"fieldName":field_SuCOS_Score, "values": {"source":source, "description":"SuCOS score"}} |
| 89 | + ] |
| 90 | + |
| 91 | + input,output,suppl,writer,output_base = rdkit_utils.\ |
| 92 | + default_open_input_output(args.input, args.informat, args.output, |
| 93 | + 'sucos', args.outformat, |
| 94 | + valueClassMappings=clsMappings, |
| 95 | + datasetMetaProps=datasetMetaProps, |
| 96 | + fieldMetaProps=fieldMetaProps) |
| 97 | + |
| 98 | + count = 0 |
| 99 | + total = 0 |
| 100 | + errors = 0 |
| 101 | + for mol in suppl: |
| 102 | + count +=1 |
| 103 | + if mol is None: |
| 104 | + continue |
| 105 | + #utils.log("Mol has", str(mol.GetNumHeavyAtoms()), "heavy atoms") |
| 106 | + try: |
| 107 | + fm_score = get_SucosScore(ref_mol, mol, field_SuCOS_Score) |
| 108 | + utils.log("Score:", str(fm_score)) |
| 109 | + writer.write(mol) |
| 110 | + total +=1 |
| 111 | + except ValueError as e: |
| 112 | + errors +=1 |
| 113 | + utils.log("Molecule", count, "failed to score:", e.message) |
| 114 | + |
| 115 | + input.close() |
| 116 | + writer.flush() |
| 117 | + writer.close() |
| 118 | + output.close() |
| 119 | + |
| 120 | + if args.meta: |
| 121 | + utils.write_metrics(output_base, {'__InputCount__':count, '__OutputCount__':total, '__ErrorCount__':errors, 'RDKitSuCOS':total}) |
| 122 | + |
| 123 | +if __name__ == "__main__": |
| 124 | + main() |
| 125 | + |
0 commit comments