1 #!/bin/env python 2 # 3 # File: RDKitPerformRGroupDecomposition.py 4 # Author: Manish Sud <msud@san.rr.com> 5 # 6 # Copyright (C) 2026 Manish Sud. All rights reserved. 7 # 8 # The functionality available in this script is implemented using RDKit, an 9 # open source toolkit for cheminformatics developed by Greg Landrum. 10 # 11 # This file is part of MayaChemTools. 12 # 13 # MayaChemTools is free software; you can redistribute it and/or modify it under 14 # the terms of the GNU Lesser General Public License as published by the Free 15 # Software Foundation; either version 3 of the License, or (at your option) any 16 # later version. 17 # 18 # MayaChemTools is distributed in the hope that it will be useful, but without 19 # any warranty; without even the implied warranty of merchantability of fitness 20 # for a particular purpose. See the GNU Lesser General Public License for more 21 # details. 22 # 23 # You should have received a copy of the GNU Lesser General Public License 24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or 25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330, 26 # Boston, MA, 02111-1307, USA. 27 # 28 29 from __future__ import print_function 30 31 import os 32 import sys 33 import time 34 import re 35 36 # RDKit imports... 37 try: 38 from rdkit import rdBase 39 from rdkit import Chem 40 from rdkit.Chem import AllChem 41 from rdkit.Chem import rdRGroupDecomposition as rgd 42 from rdkit.Chem import rdFMCS 43 except ImportError as ErrMsg: 44 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg) 45 sys.stderr.write("Check/update your RDKit environment and try again.\n\n") 46 sys.exit(1) 47 48 # MayaChemTools imports... 49 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python")) 50 try: 51 from docopt import docopt 52 import MiscUtil 53 import RDKitUtil 54 except ImportError as ErrMsg: 55 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg) 56 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n") 57 sys.exit(1) 58 59 ScriptName = os.path.basename(sys.argv[0]) 60 Options = {} 61 OptionsInfo = {} 62 63 64 def main(): 65 """Start execution of the script.""" 66 67 MiscUtil.PrintInfo( 68 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 69 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()) 70 ) 71 72 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 73 74 # Retrieve command line arguments and options... 75 RetrieveOptions() 76 77 # Process and validate command line arguments and options... 78 ProcessOptions() 79 80 # Perform actions required by the script... 81 PerformRGroupDecomposition() 82 83 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 84 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 85 86 87 def PerformRGroupDecomposition(): 88 """Perform R group decomposition.""" 89 90 # Retrieve molecules... 91 Mols = RetrieveMolecules() 92 93 # Identify R groups and write them out... 94 RGroups, UnmatchedMolIndices = RetrieveRGroups(Mols) 95 WriteRGroups(Mols, RGroups, UnmatchedMolIndices) 96 97 98 def RetrieveRGroups(Mols): 99 """Retrieve R groups.""" 100 101 CoreMols = SetupCoreScaffolds(Mols) 102 DecompositionParams = SetupRGroupDecompositionParams() 103 RGroupDecompositionObject = rgd.RGroupDecomposition(CoreMols, DecompositionParams) 104 105 MiscUtil.PrintInfo("\nPerforming R group decomposition...") 106 107 UnmatchedMolIndices = [] 108 for MolIndex, Mol in enumerate(Mols): 109 Status = RGroupDecompositionObject.Add(Mol) 110 if Status < 0: 111 UnmatchedMolIndices.append(MolIndex) 112 113 if not RGroupDecompositionObject.Process(): 114 MiscUtil.PrintWarning("R group decomposition failed to match any molecule to core scaffold(s)...") 115 116 RGroups = RGroupDecompositionObject.GetRGroupsAsColumns(asSmiles=True) 117 118 return (RGroups, UnmatchedMolIndices) 119 120 121 def SetupCoreScaffolds(Mols): 122 """Setup core scaffold molecules(s).""" 123 124 if re.match("^(BySMARTS|BySMILES)$", OptionsInfo["CoreScaffold"], re.I): 125 return SetupCoreScaffoldsBySMARTSOrSMILES() 126 elif re.match("^ByMCS$", OptionsInfo["CoreScaffold"], re.I): 127 return SetupCoreScaffoldsByMCS(Mols) 128 else: 129 MiscUtil.PrintError( 130 'The value, %s, specified for "-c, --coreScaffold" option is not supported.' 131 % (OptionsInfo["CoreScaffold"]) 132 ) 133 134 135 def SetupCoreScaffoldsBySMARTSOrSMILES(): 136 """Setup core scaffold molecules(s) using specified SMARTS or SMILES.""" 137 138 BySMARTS = True if re.match("^BySMARTS$", OptionsInfo["CoreScaffold"], re.I) else False 139 CoreScaffoldList = OptionsInfo["SMARTSOrSMILESCoreScaffoldList"] 140 141 if BySMARTS: 142 MiscUtil.PrintInfo( 143 "\nSetting up core scaffold(s) using SMARTS...\nSMARTS core scaffold(s): %s" % " ".join(CoreScaffoldList) 144 ) 145 else: 146 MiscUtil.PrintInfo( 147 "\nSetting up core scaffold(s) using SMILES...\nSMILES core scaffold(s): %s" % " ".join(CoreScaffoldList) 148 ) 149 150 CoreMols = [] 151 for Core in CoreScaffoldList: 152 if BySMARTS: 153 CoreMol = Chem.MolFromSmarts(Core) 154 else: 155 CoreMol = Chem.MolFromSmiles(Core) 156 if CoreMol is None: 157 MiscUtil.PrintError("Failed to generate mol for core scaffold: %s" % (Core)) 158 CoreMols.append(CoreMol) 159 160 return CoreMols 161 162 163 def SetupCoreScaffoldsByMCS(Mols): 164 """Setup core scaffold molecule using MCS.""" 165 166 MiscUtil.PrintInfo("\nSetting up core scaffold using MCS...") 167 168 MCSParams = OptionsInfo["MCSParams"] 169 170 CoreMols = [] 171 172 MCSResultObject = rdFMCS.FindMCS( 173 Mols, 174 maximizeBonds=MCSParams["MaximizeBonds"], 175 threshold=MCSParams["Threshold"], 176 timeout=MCSParams["TimeOut"], 177 verbose=MCSParams["Verbose"], 178 matchValences=MCSParams["MatchValences"], 179 ringMatchesRingOnly=MCSParams["RingMatchesRingOnly"], 180 completeRingsOnly=MCSParams["CompleteRingsOnly"], 181 matchChiralTag=MCSParams["MatchChiralTag"], 182 atomCompare=MCSParams["AtomCompare"], 183 bondCompare=MCSParams["BondCompare"], 184 seedSmarts=MCSParams["SeedSMARTS"], 185 ) 186 187 if MCSResultObject.canceled: 188 MiscUtil.PrintError( 189 'MCS failed to identify a core scaffold. Specify a different set of parameters using "-m, --mcsParams" option and try again.' 190 ) 191 192 CoreNumAtoms = MCSResultObject.numAtoms 193 CoreNumBonds = MCSResultObject.numBonds 194 SMARTSCore = MCSResultObject.smartsString 195 196 if not len(SMARTSCore): 197 MiscUtil.PrintError( 198 'MCS failed to identify a core scaffold. Specify a different set of parameters using "-m, --mcsParams" option and try again.' 199 ) 200 201 MiscUtil.PrintInfo( 202 "SMARTS core scaffold: %s\nNumber of atoms in core scaffold: %s\nNumber of bonds in core scaffold: %s" 203 % (SMARTSCore, CoreNumAtoms, CoreNumBonds) 204 ) 205 206 if CoreNumAtoms < MCSParams["MinNumAtoms"]: 207 MiscUtil.PrintError( 208 'Number of atoms, %d, in core scaffold identified by MCS is less than, %d, as specified by "minNumAtoms" parameter in "-m, --mcsParams" option.' 209 % (CoreNumAtoms, MCSParams["MinNumAtoms"]) 210 ) 211 212 if CoreNumBonds < MCSParams["MinNumBonds"]: 213 MiscUtil.PrintError( 214 'Number of bonds, %d, in core scaffold identified by MCS is less than, %d, as specified by "minNumBonds" parameter in "-m, --mcsParams" option.' 215 % (CoreNumBonds, MCSParams["MinNumBonds"]) 216 ) 217 218 CoreMol = Chem.MolFromSmarts(SMARTSCore) 219 CoreMols.append(CoreMol) 220 221 return CoreMols 222 223 224 def SetupRGroupDecompositionParams(): 225 """Setup R group decomposition parameters.""" 226 227 DecompositionParams = rgd.RGroupDecompositionParameters() 228 229 DecompositionParams.alignment = OptionsInfo["DecompositionParams"]["RGroupCoreAlignment"] 230 DecompositionParams.chunkSize = OptionsInfo["DecompositionParams"]["chunkSize"] 231 DecompositionParams.matchingStrategy = OptionsInfo["DecompositionParams"]["RGroupMatching"] 232 DecompositionParams.onlyMatchAtRGroups = OptionsInfo["DecompositionParams"]["matchOnlyAtRGroups"] 233 DecompositionParams.removeAllHydrogenRGroups = OptionsInfo["DecompositionParams"]["removeHydrogenOnlyGroups"] 234 DecompositionParams.removeHydrogensPostMatch = OptionsInfo["DecompositionParams"]["removeHydrogensPostMatch"] 235 236 return DecompositionParams 237 238 239 def WriteRGroups(Mols, RGroups, UnmatchedMolIndices): 240 """Write out R groups.""" 241 242 Outfile = OptionsInfo["Outfile"] 243 UnmatchedOutfile = OptionsInfo["UnmatchedOutfile"] 244 RemoveUnmatchedMode = OptionsInfo["RemoveUnmatchedMode"] 245 246 Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"] 247 248 TextOutFileMode = OptionsInfo["TextOutFileMode"] 249 TextOutFileDelim = OptionsInfo["TextOutFileDelim"] 250 Quote = OptionsInfo["TextOutQuote"] 251 252 SMILESIsomeric = OptionsInfo["OutfileParams"]["SMILESIsomeric"] 253 SMILESKekulize = OptionsInfo["OutfileParams"]["SMILESKekulize"] 254 255 # Setup writers... 256 Writer = None 257 UnmatchedWriter = None 258 if TextOutFileMode: 259 Writer = open(Outfile, "w") 260 if RemoveUnmatchedMode: 261 UnmatchedWriter = open(UnmatchedOutfile, "w") 262 else: 263 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"]) 264 if RemoveUnmatchedMode: 265 UnmatchedWriter = RDKitUtil.MoleculesWriter(UnmatchedOutfile, **OptionsInfo["OutfileParams"]) 266 267 if Writer is None: 268 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile) 269 if RemoveUnmatchedMode: 270 if UnmatchedWriter is None: 271 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % UnmatchedOutfile) 272 273 if RemoveUnmatchedMode: 274 MiscUtil.PrintInfo("\nGenerating files: %s %s..." % (Outfile, UnmatchedOutfile)) 275 else: 276 MiscUtil.PrintInfo("\nGenerating file %s..." % Outfile) 277 278 # Set up data keys and labels for core and R groups... 279 RGroupsDataKeys = [] 280 RGroupsDataLabels = [] 281 CoreDataLabelPresent = False 282 for DataLabel in sorted(RGroups): 283 if re.match("^Core", DataLabel, re.I): 284 CoreDataLabelPresent = True 285 RGroupsDataKeys.append(DataLabel) 286 RGroupsDataLabels.append("SMILES%s" % DataLabel) 287 elif re.match("^R", DataLabel, re.I): 288 RGroupsDataKeys.append(DataLabel) 289 RGroupsDataLabels.append("SMILES%s" % DataLabel) 290 else: 291 MiscUtil.PrintWarning( 292 "Ignoring unknown R group data label, %s, found during R group decomposition..." % DataLabel 293 ) 294 295 if CoreDataLabelPresent: 296 RGroupsCategoriesCount = len(RGroupsDataLabels) - 1 297 else: 298 RGroupsCategoriesCount = len(RGroupsDataLabels) 299 RGroupsMolUnmatchedCount = len(UnmatchedMolIndices) 300 RGroupsMolMatchedCount = len(Mols) - RGroupsMolUnmatchedCount 301 302 # Wite out headers for a text file... 303 if TextOutFileMode: 304 LineWords = ["SMILES", "Name"] 305 306 if RemoveUnmatchedMode: 307 Line = MiscUtil.JoinWords(LineWords, TextOutFileDelim, Quote) 308 UnmatchedWriter.write("%s\n" % Line) 309 310 LineWords.extend(RGroupsDataLabels) 311 Line = MiscUtil.JoinWords(LineWords, TextOutFileDelim, Quote) 312 Writer.write("%s\n" % Line) 313 314 MolCount = 0 315 RGroupsResultIndex = -1 316 317 for MolIndex, Mol in enumerate(Mols): 318 MolCount += 1 319 320 UnmatchedMol = False 321 if MolIndex in UnmatchedMolIndices: 322 UnmatchedMol = True 323 324 if UnmatchedMol: 325 RGroupsDataSMILES = [""] * len(RGroupsDataKeys) 326 else: 327 RGroupsResultIndex += 1 328 RGroupsDataSMILES = [RGroups[RGroupsDataKey][RGroupsResultIndex] for RGroupsDataKey in RGroupsDataKeys] 329 330 if TextOutFileMode: 331 # Write out text file including SMILES file... 332 MolSMILES = Chem.MolToSmiles(Mol, isomericSmiles=SMILESIsomeric, kekuleSmiles=SMILESKekulize) 333 MolName = RDKitUtil.GetMolName(Mol, MolCount) 334 LineWords = [MolSMILES, MolName] 335 336 if UnmatchedMol and RemoveUnmatchedMode: 337 Line = MiscUtil.JoinWords(LineWords, TextOutFileDelim, Quote) 338 UnmatchedWriter.write("%s\n" % Line) 339 else: 340 LineWords.extend(RGroupsDataSMILES) 341 Line = MiscUtil.JoinWords(LineWords, TextOutFileDelim, Quote) 342 Writer.write("%s\n" % Line) 343 else: 344 # Write out SD file... 345 if Compute2DCoords: 346 AllChem.Compute2DCoords(Mol) 347 348 if UnmatchedMol and RemoveUnmatchedMode: 349 UnmatchedWriter.write(Mol) 350 else: 351 for Name, Value in zip(RGroupsDataLabels, RGroupsDataSMILES): 352 Mol.SetProp(Name, Value) 353 Writer.write(Mol) 354 355 if Writer is not None: 356 Writer.close() 357 if UnmatchedWriter is not None: 358 UnmatchedWriter.close() 359 360 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount) 361 MiscUtil.PrintInfo("Number of R group categories: %d" % RGroupsCategoriesCount) 362 MiscUtil.PrintInfo("Number of matched molecules containing core scaffold(s): %d" % RGroupsMolMatchedCount) 363 MiscUtil.PrintInfo("Number of unmatched molecules containing no core scaffold(s): %d" % RGroupsMolUnmatchedCount) 364 365 366 def RetrieveMolecules(): 367 """Retrieve molecules.""" 368 369 Infile = OptionsInfo["Infile"] 370 371 # Read molecules... 372 MiscUtil.PrintInfo("\nReading file %s..." % Infile) 373 OptionsInfo["InfileParams"]["AllowEmptyMols"] = False 374 ValidMols, MolCount, ValidMolCount = RDKitUtil.ReadAndValidateMolecules(Infile, **OptionsInfo["InfileParams"]) 375 376 MiscUtil.PrintInfo("Total number of molecules: %d" % MolCount) 377 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount) 378 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount)) 379 380 return ValidMols 381 382 383 def ProcessDecompositionParameters(): 384 """Set up and process decomposition parameters.""" 385 386 SetupDecompositionParameters() 387 ProcessSpecifiedDecompositionParameters() 388 389 390 def SetupDecompositionParameters(): 391 """Set up default decomposition parameters.""" 392 393 OptionsInfo["DecompositionParams"] = { 394 "RGroupCoreAlignment": rgd.RGroupCoreAlignment.MCS, 395 "RGroupMatching": rgd.RGroupMatching.GreedyChunks, 396 "chunkSize": 5, 397 "matchOnlyAtRGroups": False, 398 "removeHydrogenOnlyGroups": True, 399 "removeHydrogensPostMatch": False, 400 } 401 402 403 def ProcessSpecifiedDecompositionParameters(): 404 """Process specified decomposition parameters.""" 405 406 if re.match("^auto$", OptionsInfo["SpecifiedDecompositionParams"], re.I): 407 # Nothing to process... 408 return 409 410 # Parse specified parameters... 411 DecompositionParams = re.sub(" ", "", OptionsInfo["SpecifiedDecompositionParams"]) 412 if not DecompositionParams: 413 MiscUtil.PrintError( 414 'No valid parameter name and value pairs specified using "-d, --decompositionParams" option.' 415 ) 416 417 DecompositionParamsWords = DecompositionParams.split(",") 418 if len(DecompositionParamsWords) % 2: 419 MiscUtil.PrintError( 420 'The number of comma delimited paramater names and values, %d, specified using "-d, --decompositionParams" option must be an even number.' 421 % (len(DecompositionParamsWords)) 422 ) 423 424 # Setup canonical parameter names... 425 ValidParamNames = [] 426 CanonicalParamNamesMap = {} 427 for ParamName in sorted(OptionsInfo["DecompositionParams"]): 428 ValidParamNames.append(ParamName) 429 CanonicalParamNamesMap[ParamName.lower()] = ParamName 430 431 # Validate and set paramater names and value... 432 for Index in range(0, len(DecompositionParamsWords), 2): 433 Name = DecompositionParamsWords[Index] 434 Value = DecompositionParamsWords[Index + 1] 435 436 CanonicalName = Name.lower() 437 if CanonicalName not in CanonicalParamNamesMap: 438 MiscUtil.PrintError( 439 'The parameter name, %s, specified using "-d, --decompositionParams" option is not a valid name. Supported parameter names: %s' 440 % (Name, " ".join(ValidParamNames)) 441 ) 442 443 ParamName = CanonicalParamNamesMap[CanonicalName] 444 if re.match("^RGroupCoreAlignment$", ParamName, re.I): 445 if re.match("^MCS$", Value, re.I): 446 ParamValue = rgd.RGroupCoreAlignment.MCS 447 elif re.match("^None$", Value, re.I): 448 ParamValue = rgd.RGroupCoreAlignment.names["None"] 449 else: 450 MiscUtil.PrintError( 451 'The parameter value, %s, specified using "-d, --decompositionParams" option for parameter, %s, is not a valid value. Supported values: MCS None' 452 % (Value, Name) 453 ) 454 elif re.match("^RGroupMatching$", ParamName, re.I): 455 if re.match("^Greedy$", Value, re.I): 456 ParamValue = rgd.RGroupMatching.Greedy 457 elif re.match("^GreedyChunks$", Value, re.I): 458 ParamValue = rgd.RGroupMatching.GreedyChunks 459 elif re.match("^Exhaustive$", Value, re.I): 460 ParamValue = rgd.RGroupMatching.Exhaustive 461 else: 462 MiscUtil.PrintError( 463 'The parameter value, %s, specified using "-d, --decompositionParams" option for parameter, %s, is not a valid value. Supported values: Greedy GreedyChunks Exhaustive' 464 % (Value, Name) 465 ) 466 elif re.match("^chunkSize$", ParamName, re.I): 467 Value = int(Value) 468 if Value <= 0: 469 MiscUtil.PrintError( 470 'The parameter value, %s, specified using "-d, --decompositionParams" option for parameter, %s, is not a valid value. Supported values: > 0' 471 % (Value, Name) 472 ) 473 ParamValue = Value 474 else: 475 if not re.match("^(Yes|No|True|False)$", Value, re.I): 476 MiscUtil.PrintError( 477 'The parameter value, %s, specified using "-d, --decompositionParams" option for parameter, %s, is not a valid value. Supported values: Yes No True False' 478 % (Value, Name) 479 ) 480 ParamValue = False 481 if re.match("^(Yes|True)$", Value, re.I): 482 ParamValue = True 483 484 # Set value... 485 OptionsInfo["DecompositionParams"][ParamName] = ParamValue 486 487 488 def ProcessMCSParameters(): 489 """Set up and process MCS parameters.""" 490 491 SetupMCSParameters() 492 ProcessSpecifiedMCSParameters() 493 494 495 def SetupMCSParameters(): 496 """Set up default MCS parameters.""" 497 498 OptionsInfo["MCSParams"] = { 499 "MaximizeBonds": True, 500 "Threshold": 0.9, 501 "TimeOut": 3600, 502 "Verbose": False, 503 "MatchValences": True, 504 "MatchChiralTag": False, 505 "RingMatchesRingOnly": True, 506 "CompleteRingsOnly": True, 507 "AtomCompare": rdFMCS.AtomCompare.CompareElements, 508 "BondCompare": rdFMCS.BondCompare.CompareOrder, 509 "SeedSMARTS": "", 510 "MinNumAtoms": 1, 511 "MinNumBonds": 0, 512 } 513 514 515 def ProcessSpecifiedMCSParameters(): 516 """Process specified MCS parameters.""" 517 518 if re.match("^auto$", OptionsInfo["SpecifiedMCSParams"], re.I): 519 # Nothing to process... 520 return 521 522 # Parse specified parameters... 523 MCSParams = re.sub(" ", "", OptionsInfo["SpecifiedMCSParams"]) 524 if not MCSParams: 525 MiscUtil.PrintError('No valid parameter name and value pairs specified using "-m, --mcsParams" option.') 526 527 MCSParamsWords = MCSParams.split(",") 528 if len(MCSParamsWords) % 2: 529 MiscUtil.PrintError( 530 'The number of comma delimited paramater names and values, %d, specified using "-m, --mcsParams" option must be an even number.' 531 % (len(MCSParamsWords)) 532 ) 533 534 # Setup canonical parameter names... 535 ValidParamNames = [] 536 CanonicalParamNamesMap = {} 537 for ParamName in sorted(OptionsInfo["MCSParams"]): 538 ValidParamNames.append(ParamName) 539 CanonicalParamNamesMap[ParamName.lower()] = ParamName 540 541 # Validate and set paramater names and value... 542 for Index in range(0, len(MCSParamsWords), 2): 543 Name = MCSParamsWords[Index] 544 Value = MCSParamsWords[Index + 1] 545 546 CanonicalName = Name.lower() 547 if CanonicalName not in CanonicalParamNamesMap: 548 MiscUtil.PrintError( 549 'The parameter name, %s, specified using "-m, --mcsParams" option is not a valid name. Supported parameter names: %s' 550 % (Name, " ".join(ValidParamNames)) 551 ) 552 553 ParamName = CanonicalParamNamesMap[CanonicalName] 554 if re.match("^Threshold$", ParamName, re.I): 555 Value = float(Value) 556 if Value <= 0.0 or Value > 1.0: 557 MiscUtil.PrintError( 558 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: > 0 and <= 1.0' 559 % (Value, Name) 560 ) 561 ParamValue = Value 562 elif re.match("^Timeout$", ParamName, re.I): 563 Value = int(Value) 564 if Value <= 0: 565 MiscUtil.PrintError( 566 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: > 0' 567 % (Value, Name) 568 ) 569 ParamValue = Value 570 elif re.match("^MinNumAtoms$", ParamName, re.I): 571 Value = int(Value) 572 if Value < 1: 573 MiscUtil.PrintError( 574 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: >= 1' 575 % (Value, Name) 576 ) 577 ParamValue = Value 578 elif re.match("^MinNumBonds$", ParamName, re.I): 579 Value = int(Value) 580 if Value < 0: 581 MiscUtil.PrintError( 582 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: >=0 ' 583 % (Value, Name) 584 ) 585 ParamValue = Value 586 elif re.match("^AtomCompare$", ParamName, re.I): 587 if re.match("^CompareAny$", Value, re.I): 588 ParamValue = rdFMCS.AtomCompare.CompareAny 589 elif re.match("^CompareElements$", Value, re.I): 590 ParamValue = Chem.rdFMCS.AtomCompare.CompareElements 591 elif re.match("^CompareIsotopes$", Value, re.I): 592 ParamValue = Chem.rdFMCS.AtomCompare.CompareIsotopes 593 else: 594 MiscUtil.PrintError( 595 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: CompareAny CompareElements CompareIsotopes' 596 % (Value, Name) 597 ) 598 elif re.match("^BondCompare$", ParamName, re.I): 599 if re.match("^CompareAny$", Value, re.I): 600 ParamValue = Chem.rdFMCS.BondCompare.CompareAny 601 elif re.match("^CompareOrder$", Value, re.I): 602 ParamValue = rdFMCS.BondCompare.CompareOrder 603 elif re.match("^CompareOrderExact$", Value, re.I): 604 ParamValue = rdFMCS.BondCompare.CompareOrderExact 605 else: 606 MiscUtil.PrintError( 607 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: CompareAny CompareOrder CompareOrderExact' 608 % (Value, Name) 609 ) 610 elif re.match("^SeedSMARTS$", ParamName, re.I): 611 if not len(Value): 612 MiscUtil.PrintError( 613 'The parameter value specified using "-m, --mcsParams" option for parameter, %s, is empty. ' 614 % (Name) 615 ) 616 ParamValue = Value 617 else: 618 if not re.match("^(Yes|No|True|False)$", Value, re.I): 619 MiscUtil.PrintError( 620 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: Yes No True False' 621 % (Value, Name) 622 ) 623 ParamValue = False 624 if re.match("^(Yes|True)$", Value, re.I): 625 ParamValue = True 626 627 # Set value... 628 OptionsInfo["MCSParams"][ParamName] = ParamValue 629 630 631 def ProcessOptions(): 632 """Process and validate command line arguments and options.""" 633 634 MiscUtil.PrintInfo("Processing options...") 635 636 # Validate options... 637 ValidateOptions() 638 639 OptionsInfo["CoreScaffold"] = Options["--coreScaffold"] 640 641 OptionsInfo["Infile"] = Options["--infile"] 642 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters( 643 "--infileParams", Options["--infileParams"], Options["--infile"] 644 ) 645 646 OptionsInfo["Outfile"] = Options["--outfile"] 647 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters( 648 "--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"] 649 ) 650 651 TextOutFileMode = False 652 TextOutFileDelim = "" 653 654 if MiscUtil.CheckFileExt(Options["--outfile"], "csv"): 655 TextOutFileMode = True 656 TextOutFileDelim = "," 657 elif MiscUtil.CheckFileExt(Options["--outfile"], "tsv txt"): 658 TextOutFileMode = True 659 TextOutFileDelim = "\t" 660 661 OptionsInfo["TextOutFileMode"] = TextOutFileMode 662 OptionsInfo["TextOutFileDelim"] = TextOutFileDelim 663 664 TextOutQuote = False 665 if re.match("^auto$", Options["--quote"], re.I): 666 if MiscUtil.CheckFileExt(Options["--outfile"], "csv"): 667 TextOutQuote = True 668 else: 669 if re.match("^yes$", Options["--quote"], re.I): 670 TextOutQuote = True 671 OptionsInfo["TextOutQuote"] = TextOutQuote 672 673 OptionsInfo["Overwrite"] = Options["--overwrite"] 674 675 RemoveUnmatchedMode = False 676 UnmatchedOutfile = None 677 if re.match("^yes$", Options["--removeUnmatched"], re.I): 678 RemoveUnmatchedMode = True 679 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Outfile"]) 680 UnmatchedOutfile = "%sUnmatched.%s" % (FileName, FileExt) 681 OptionsInfo["RemoveUnmatchedMode"] = RemoveUnmatchedMode 682 OptionsInfo["UnmatchedOutfile"] = UnmatchedOutfile 683 684 OptionsInfo["SpecifiedDecompositionParams"] = Options["--decompositionParams"] 685 ProcessDecompositionParameters() 686 687 OptionsInfo["SpecifiedMCSParams"] = Options["--mcsParams"] 688 ProcessMCSParameters() 689 690 SMARTSOrSMILESCoreScaffold = "" 691 SMARTSOrSMILESCoreScaffoldList = [] 692 if not re.match("^none$", Options["--smartsOrSmilesCoreScaffold"], re.I) or len( 693 Options["--smartsOrSmilesCoreScaffold"] 694 ): 695 if re.match("^(BySMARTS|BySMILES)$", Options["--coreScaffold"], re.I): 696 SMARTSOrSMILESCoreScaffold = re.sub(" ", "", Options["--smartsOrSmilesCoreScaffold"]) 697 if not SMARTSOrSMILESCoreScaffold: 698 MiscUtil.PrintError( 699 'A non empty value must be specified for "-s, --smartsOrSmilesCoreScaffold" during %s value of "-c, --coreScaffold" option ' 700 % (Options["--coreScaffold"]) 701 ) 702 SMARTSOrSMILESCoreScaffoldList = SMARTSOrSMILESCoreScaffold.split(",") 703 OptionsInfo["SMARTSOrSMILESCoreScaffold"] = SMARTSOrSMILESCoreScaffold 704 OptionsInfo["SMARTSOrSMILESCoreScaffoldList"] = SMARTSOrSMILESCoreScaffoldList 705 706 707 def RetrieveOptions(): 708 """Retrieve command line arguments and options.""" 709 710 # Get options... 711 global Options 712 Options = docopt(_docoptUsage_) 713 714 # Set current working directory to the specified directory... 715 WorkingDir = Options["--workingdir"] 716 if WorkingDir: 717 os.chdir(WorkingDir) 718 719 # Handle examples option... 720 if "--examples" in Options and Options["--examples"]: 721 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 722 sys.exit(0) 723 724 725 def ValidateOptions(): 726 """Validate option values.""" 727 728 MiscUtil.ValidateOptionTextValue("-c, --coreScaffold", Options["--coreScaffold"], "ByMCS BySMARTS BySMILES") 729 730 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"]) 731 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi csv tsv txt") 732 733 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd csv tsv txt") 734 MiscUtil.ValidateOptionsOutputFileOverwrite( 735 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 736 ) 737 MiscUtil.ValidateOptionsDistinctFileNames( 738 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"] 739 ) 740 741 if re.match("^none$", Options["--smartsOrSmilesCoreScaffold"], re.I) or not len( 742 Options["--smartsOrSmilesCoreScaffold"] 743 ): 744 if re.match("^(BySMARTS|BySMILES)$", Options["--coreScaffold"], re.I): 745 MiscUtil.PrintError( 746 'A non empty value must be specified for "-s, --smartsOrSmilesCoreScaffold" during %s value of "-c, --coreScaffold" option ' 747 % (Options["--coreScaffold"]) 748 ) 749 else: 750 if not re.match("^(BySMARTS|BySMILES)$", Options["--coreScaffold"], re.I): 751 MiscUtil.PrintError( 752 '%s value of "-s, --smartsOrSmilesCoreScaffold" is not allowed during %s value of "-c, --coreScaffold" option ' 753 % (Options["--smartsOrSmilesCoreScaffold"], Options["--coreScaffold"]) 754 ) 755 756 MiscUtil.ValidateOptionTextValue("-q, --quote", Options["--quote"], "yes no auto") 757 MiscUtil.ValidateOptionTextValue("-r, --removeUnmatched", Options["--removeUnmatched"], "yes no") 758 759 760 # Setup a usage string for docopt... 761 _docoptUsage_ = """ 762 RDKitPerformRGroupDecomposition.py - Perform R group decomposition analysis 763 764 Usage: 765 RDKitPerformRGroupDecomposition.py [--coreScaffold <ByMCS, BySMARTS or BySMILES>] 766 [--decompositionParams <Name,Value,...>] 767 [--infileParams <Name,Value,...>] [--mcsParams <Name,Value,...>] 768 [--outfileParams <Name,Value,...>] [--overwrite] [--quote <yes or no>] 769 [--removeUnmatched <yes or no>] [--smartsOrSmilesCoreScaffold <text>] 770 [-w <dir>] -i <infile> -o <outfile> 771 RDKitPerformRGroupDecomposition.py -h | --help | -e | --examples 772 773 Description: 774 Perform R group decomposition for a set of molecules in a series containing 775 a common core scaffold. The core scaffold is identified by a SMARTS string, 776 SMILES string, or using maximum common substructure (MCS) search. 777 Multiple core scaffolds may be specified using SMARTS or SMILES strings for 778 set of molecules corresponding to multiple series. 779 780 The core scaffolds along with appropriate R groups are written out as SMILES 781 strings to a SD or text file. The unmatched molecules without any specified 782 core scaffold are written to a different output file. 783 784 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi, 785 .txt, .csv, .tsv) 786 787 The supported output file formats are: SD (.sdf, .sd), CSV/TSV (.csv, .tsv, .txt) 788 789 Options: 790 -c, --coreScaffold <ByMCS, BySMARTS or BySMILES> [default: ByMCS] 791 Specify a core scaffold for a set of molecules in a series. The core scaffold 792 is identified by an explicit SMARTS string, SMILES string, or using maximum 793 common substructure (MCS) search. Multiple core scaffolds may be specified 794 using SMARTS or SMILES strings for set of molecules corresponding to multiple 795 series. 796 -d, --decompositionParams <Name,Value,...> [default: auto] 797 Parameter values to use during R group decomposition for a series of molecules. 798 In general, it is a comma delimited list of parameter name and value pairs. The 799 supported parameter names along with their default values are shown below: 800 801 RGroupCoreAlignment,MCS, RGroupMatching,GreedyChunks,chunkSize,5, 802 matchOnlyAtRGroups,no,removeHydrogenOnlyGroups,yes, 803 removeHydrogensPostMatch,no 804 805 A brief description of each supported parameter taken from RDKit documentation, 806 along with their possible values, is as follows. 807 808 RGroupCoreAlignment - Mapping of core labels: 809 810 MCS - Map core labels to each other using MCS 811 None - No mapping 812 813 RGroupMatching: Greedy, GreedyChunks, Exhaustive 814 815 matchOnlyAtRGroups - Allow R group decomposition only at specified R groups. 816 Possible values: yes, no. 817 818 removeHydrogenOnlyGroups - Remove all R groups that only have hydrogens. 819 Possible values: yes, no. 820 821 removeHydrogensPostMatch - Remove all hydrogens from the output molecules. 822 Possible values: yes, no. 823 -e, --examples 824 Print examples. 825 -h, --help 826 Print this help message. 827 -i, --infile <infile> 828 Input file name. 829 --infileParams <Name,Value,...> [default: auto] 830 A comma delimited list of parameter name and value pairs for reading 831 molecules from files. The supported parameter names for different file 832 formats, along with their default values, are shown below: 833 834 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes 835 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space, 836 smilesTitleLine,auto,sanitize,yes 837 838 Possible values for smilesDelimiter: space, comma or tab. 839 -m, --mcsParams <Name,Value,...> [default: auto] 840 Parameter values to use for identifying a maximum common substructure 841 (MCS) in a series of molecules. In general, it is a comma delimited list of 842 parameter name and value pairs. The supported parameter names along with 843 their default values are shown below: 844 845 atomCompare,CompareElements,bondCompare,CompareOrder, 846 maximizeBonds,yes,matchValences,yes,matchChiralTag,no, 847 minNumAtoms,1,minNumBonds,0,ringMatchesRingOnly,yes, 848 completeRingsOnly,yes,threshold,1.0,timeOut,3600,seedSMARTS,none 849 850 Possible values for atomCompare: CompareAny, CompareElements, 851 CompareIsotopes. Possible values for bondCompare: CompareAny, 852 CompareOrder, CompareOrderExact. 853 854 A brief description of MCS parameters taken from RDKit documentation is 855 as follows: 856 857 atomCompare - Controls match between two atoms 858 bondCompare - Controls match between two bonds 859 maximizeBonds - Maximize number of bonds instead of atoms 860 matchValences - Include atom valences in the MCS match 861 matchChiralTag - Include atom chirality in the MCS match 862 minNumAtoms - Minimum number of atoms in the MCS match 863 minNumBonds - Minimum number of bonds in the MCS match 864 ringMatchesRingOnly - Ring bonds only match other ring bonds 865 completeRingsOnly - Partial rings not allowed during the match 866 threshold - Fraction of the dataset that must contain the MCS 867 seedSMARTS - SMARTS string as the seed of the MCS 868 timeout - Timeout for the MCS calculation in seconds 869 870 -o, --outfile <outfile> 871 Output file name. 872 --outfileParams <Name,Value,...> [default: auto] 873 A comma delimited list of parameter name and value pairs for writing 874 molecules to files. The supported parameter names for different file 875 formats, along with their default values, are shown below: 876 877 SD: compute2DCoords,auto,kekulize,yes,forceV3000,no 878 SMILES: smilesKekulize,no,smilesIsomeric,yes 879 880 Default value for compute2DCoords: yes for SMILES input file; no for all other 881 file types. The kekulize and smilesIsomeric parameters are also used during 882 generation of SMILES strings for CSV/TSV files. 883 --overwrite 884 Overwrite existing files. 885 -q, --quote <yes or no> [default: auto] 886 Quote SMILES strings and molecule names before writing them out to text 887 files. Possible values: yes or no. Default: yes for CSV (.csv) text files; no for 888 TSV (.tsv) and TXT (.txt) text files. 889 -r, --removeUnmatched <yes or no> [default: no] 890 Remove unmatched molecules containing no specified core scaffold from the 891 output file and write them to a different output file. 892 -s, --smartsOrSmilesCoreScaffold <text> [default: none] 893 SMARTS or SMILES string to use for core scaffold during 'SMARTS' or 'SMILES' 894 value of '-c, --coreScaffold' option. Multiple core scaffolds may be specified using a 895 comma delimited set of SMARTS or SMILES strings. 896 -w, --workingdir <dir> 897 Location of working directory which defaults to the current directory. 898 899 Examples: 900 To perform R group decomposition for a set of molecules in a series using MCS 901 to identify a core scaffold and write out a CSV file containing R groups, type: 902 903 % RDKitPerformRGroupDecomposition.py -i SampleSeriesD3R.smi 904 -o SampleSeriesD3ROut.csv 905 906 To perform R group decomposition for a set of molecules in a series using a 907 specified core scaffold and write out a SD file containing R groups, type: 908 909 % RDKitPerformRGroupDecomposition.py -c BySMARTS 910 -s "Nc1nccc(-c2cnc(CNCc3ccccc3)c2)n1" -i SampleSeriesD3R.smi 911 -o SampleSeriesD3ROut.sdf 912 913 To perform R group decomposition for a set of molecules in a series using MCS 914 to identify a core scaffold and write out CSV files containing matched and 915 unmatched molecules without quoting values, type: 916 917 % RDKitPerformRGroupDecomposition.py -c ByMCS -r yes -q no 918 -i SampleSeriesD3R.sdf -o SampleSeriesD3ROut.csv 919 920 To perform R group decomposition for a set of molecules in multiple series using 921 specified core scaffolds and write out a TSV file containing R groups, type: 922 923 % RDKitPerformRGroupDecomposition.py -c BySMARTS 924 -s "Nc1nccc(-c2cnc(CNCc3ccccc3)c2)n1,[#6]-[#6]1:[#6]:[#6]:[#6]:[#6]: 925 [#6]:1" -i SampleMultipleSeriesD3R.smi -o 926 SampleMultipleSeriesD3ROut.tsv 927 928 To perform R group decomposition for a set of molecules in a CSV SMILES file, 929 SMILES strings in olumn 1, name in column 2, and write out a CSV file containing 930 R groups, type: 931 932 % RDKitPerformRGroupDecomposition.py --infileParams 933 "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1, 934 smilesNameColumn,2" --outfileParams "compute2DCoords,yes" 935 -i SampleSeriesD3R.smi -o SampleSeriesD3ROut.csv 936 937 Author: 938 Manish Sud(msud@san.rr.com) 939 940 See also: 941 RDKitConvertFileFormat.py, RDKitSearchFunctionalGroups.py, RDKitSearchSMARTS.py 942 943 Copyright: 944 Copyright (C) 2026 Manish Sud. All rights reserved. 945 946 The functionality available in this script is implemented using RDKit, an 947 open source toolkit for cheminformatics developed by Greg Landrum. 948 949 This file is part of MayaChemTools. 950 951 MayaChemTools is free software; you can redistribute it and/or modify it under 952 the terms of the GNU Lesser General Public License as published by the Free 953 Software Foundation; either version 3 of the License, or (at your option) any 954 later version. 955 956 """ 957 958 if __name__ == "__main__": 959 main()
