1 #!/bin/env python 2 # 3 # File: RDKitSearchFunctionalGroups.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 import multiprocessing as mp 36 37 # RDKit imports... 38 try: 39 from rdkit import rdBase 40 from rdkit import Chem 41 from rdkit.Chem import AllChem 42 from rdkit.Chem import FunctionalGroups 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 FunctionalGroupsMap = {} 64 65 66 def main(): 67 """Start execution of the script.""" 68 69 MiscUtil.PrintInfo( 70 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 71 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()) 72 ) 73 74 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 75 76 # Retrieve command line arguments and options... 77 RetrieveOptions() 78 79 # Process and validate command line arguments and options... 80 ProcessOptions() 81 82 # Perform actions required by the script... 83 PerformFunctionalGroupsSearch() 84 85 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 86 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 87 88 89 def PerformFunctionalGroupsSearch(): 90 """Retrieve functional groups information and perform search.""" 91 92 # Process functional groups info... 93 ProcessFunctionalGroupsInfo() 94 95 # Setup pattern mols for functional group SMARTS... 96 GroupsPatternMols = SetupFunctionalGroupsSMARTSPatterns() 97 98 # Setup a molecule reader... 99 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"]) 100 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"]) 101 102 # Set up molecule writers... 103 Writer, GroupOutfilesWriters = SetupMoleculeWriters() 104 105 MolCount, ValidMolCount, RemainingMolCount, GroupsPatternMatchCountList = ProcessMolecules( 106 Mols, GroupsPatternMols, Writer, GroupOutfilesWriters 107 ) 108 109 if Writer is not None: 110 Writer.close() 111 for GroupOutfileWriter in GroupOutfilesWriters: 112 GroupOutfileWriter.close() 113 114 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount) 115 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount) 116 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount)) 117 118 MiscUtil.PrintInfo("\nTotal number of molecules matched against specified match criteria: %d" % RemainingMolCount) 119 120 MiscUtil.PrintInfo("\nNumber of molecuels matched against individual functional groups:") 121 MiscUtil.PrintInfo("FunctionalGroupName,MatchCount") 122 123 for GroupIndex in range(0, len(OptionsInfo["SpecifiedFunctionalGroups"])): 124 GroupName = OptionsInfo["SpecifiedFunctionalGroups"][GroupIndex] 125 if OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"][GroupIndex]: 126 GroupName = "!" + GroupName 127 GroupMatchCount = GroupsPatternMatchCountList[GroupIndex] 128 MiscUtil.PrintInfo("%s,%d" % (GroupName, GroupMatchCount)) 129 130 131 def ProcessMolecules(Mols, GroupsPatternMols, Writer, GroupOutfilesWriters): 132 """Process and search molecules.""" 133 134 if OptionsInfo["MPMode"]: 135 return ProcessMoleculesUsingMultipleProcesses(Mols, GroupsPatternMols, Writer, GroupOutfilesWriters) 136 else: 137 return ProcessMoleculesUsingSingleProcess(Mols, GroupsPatternMols, Writer, GroupOutfilesWriters) 138 139 140 def ProcessMoleculesUsingSingleProcess(Mols, GroupsPatternMols, Writer, GroupOutfilesWriters): 141 """Process and search molecules using a single process.""" 142 143 MiscUtil.PrintInfo("\nSearching functional groups...") 144 145 Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"] 146 CombineMatchResults = OptionsInfo["CombineMatchResults"] 147 SetSMILESMolProps = OptionsInfo["OutfileParams"]["SetSMILESMolProps"] 148 149 GroupsPatternsMatchCountList = [0] * len(OptionsInfo["SpecifiedFunctionalGroups"]) 150 (MolCount, ValidMolCount, RemainingMolCount) = [0] * 3 151 152 FirstMol = True 153 for Mol in Mols: 154 MolCount += 1 155 156 if Mol is None: 157 continue 158 159 if RDKitUtil.IsMolEmpty(Mol): 160 MolName = RDKitUtil.GetMolName(Mol, MolCount) 161 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 162 continue 163 164 ValidMolCount += 1 165 if FirstMol: 166 FirstMol = False 167 if SetSMILESMolProps: 168 if Writer is not None: 169 RDKitUtil.SetWriterMolProps(Writer, Mol) 170 for GroupOutfileWriter in GroupOutfilesWriters: 171 if GroupOutfileWriter is not None: 172 RDKitUtil.SetWriterMolProps(GroupOutfileWriter, Mol) 173 174 # Match molecule against functional group patterns... 175 MolMatched, GroupsPatternMatchStatusList = MatchMolecule(Mol, GroupsPatternMols) 176 177 # Update functional group match count... 178 for GroupIndex, MatchStatus in enumerate(GroupsPatternMatchStatusList): 179 if MatchStatus: 180 GroupsPatternsMatchCountList[GroupIndex] += 1 181 182 if not MolMatched: 183 continue 184 185 RemainingMolCount += 1 186 WriteMolecule( 187 Writer, GroupOutfilesWriters, Mol, Compute2DCoords, CombineMatchResults, GroupsPatternMatchStatusList 188 ) 189 190 return (MolCount, ValidMolCount, RemainingMolCount, GroupsPatternsMatchCountList) 191 192 193 def ProcessMoleculesUsingMultipleProcesses(Mols, GroupsPatternMols, Writer, GroupOutfilesWriters): 194 """Process and search molecules using multiprocessing.""" 195 196 MiscUtil.PrintInfo("\nSearching functional groups using multiprocessing...") 197 198 MPParams = OptionsInfo["MPParams"] 199 Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"] 200 CombineMatchResults = OptionsInfo["CombineMatchResults"] 201 SetSMILESMolProps = OptionsInfo["OutfileParams"]["SetSMILESMolProps"] 202 203 # Setup data for initializing a worker process... 204 MiscUtil.PrintInfo("Encoding options info and functional groups pattern molecules...") 205 OptionsInfo["EncodedGroupPatternMols"] = [ 206 RDKitUtil.MolToBase64EncodedMolString(PatternMol) for PatternMol in GroupsPatternMols 207 ] 208 InitializeWorkerProcessArgs = ( 209 MiscUtil.ObjectToBase64EncodedString(Options), 210 MiscUtil.ObjectToBase64EncodedString(OptionsInfo), 211 MiscUtil.ObjectToBase64EncodedString(FunctionalGroupsMap), 212 ) 213 214 # Setup a encoded mols data iterable for a worker process... 215 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols) 216 217 # Setup process pool along with data initialization for each process... 218 MiscUtil.PrintInfo( 219 "\nConfiguring multiprocessing using %s method..." 220 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()") 221 ) 222 MiscUtil.PrintInfo( 223 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" 224 % ( 225 MPParams["NumProcesses"], 226 MPParams["InputDataMode"], 227 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]), 228 ) 229 ) 230 231 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs) 232 233 # Start processing... 234 if re.match("^Lazy$", MPParams["InputDataMode"], re.I): 235 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 236 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I): 237 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 238 else: 239 MiscUtil.PrintError( 240 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"]) 241 ) 242 243 GroupsPatternsMatchCountList = [0] * len(OptionsInfo["SpecifiedFunctionalGroups"]) 244 245 (MolCount, ValidMolCount, RemainingMolCount) = [0] * 3 246 247 FirstMol = True 248 for Result in Results: 249 MolCount += 1 250 MolIndex, EncodedMol, MolMatched, GroupsPatternMatchStatusList = Result 251 252 if EncodedMol is None: 253 continue 254 ValidMolCount += 1 255 256 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 257 258 if FirstMol: 259 FirstMol = False 260 if SetSMILESMolProps: 261 if Writer is not None: 262 RDKitUtil.SetWriterMolProps(Writer, Mol) 263 for GroupOutfileWriter in GroupOutfilesWriters: 264 if GroupOutfileWriter is not None: 265 RDKitUtil.SetWriterMolProps(GroupOutfileWriter, Mol) 266 267 # Update functional group match count... 268 for GroupIndex, MatchStatus in enumerate(GroupsPatternMatchStatusList): 269 if MatchStatus: 270 GroupsPatternsMatchCountList[GroupIndex] += 1 271 272 if not MolMatched: 273 continue 274 275 RemainingMolCount += 1 276 WriteMolecule( 277 Writer, GroupOutfilesWriters, Mol, Compute2DCoords, CombineMatchResults, GroupsPatternMatchStatusList 278 ) 279 280 return (MolCount, ValidMolCount, RemainingMolCount, GroupsPatternsMatchCountList) 281 282 283 def InitializeWorkerProcess(*EncodedArgs): 284 """Initialize data for a worker process.""" 285 286 global Options, OptionsInfo, FunctionalGroupsMap 287 288 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid()) 289 290 # Decode Options and OptionInfo... 291 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0]) 292 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1]) 293 FunctionalGroupsMap = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[2]) 294 295 # Decode ChEMBLPatternMols... 296 OptionsInfo["GroupPatternMols"] = [ 297 RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) for EncodedMol in OptionsInfo["EncodedGroupPatternMols"] 298 ] 299 300 301 def WorkerProcess(EncodedMolInfo): 302 """Process data for a worker process.""" 303 304 MolIndex, EncodedMol = EncodedMolInfo 305 306 if EncodedMol is None: 307 return [MolIndex, None, False, None] 308 309 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 310 if RDKitUtil.IsMolEmpty(Mol): 311 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1)) 312 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 313 return [MolIndex, None, False, None] 314 315 # Match molecule against functional group patterns... 316 MolMatched, GroupsPatternMatchStatusList = MatchMolecule(Mol, OptionsInfo["GroupPatternMols"]) 317 318 return [MolIndex, EncodedMol, MolMatched, GroupsPatternMatchStatusList] 319 320 321 def MatchMolecule(Mol, GroupsPatternMols): 322 """Search for functional groups in a molecule.""" 323 324 GroupsPatternMatchStatusList = [] 325 326 # Match pattern mols... 327 for GroupIndex in range(0, len(OptionsInfo["SpecifiedFunctionalGroups"])): 328 Status = DoesPatternMolMatch( 329 GroupsPatternMols[GroupIndex], 330 Mol, 331 OptionsInfo["UseChirality"], 332 OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"][GroupIndex], 333 ) 334 GroupsPatternMatchStatusList.append(Status) 335 336 # Match mol against all specified criteria... 337 MolMatched = DoesMolMeetSpecifiedMatchCriteria( 338 GroupsPatternMatchStatusList, OptionsInfo["CombineMatchResults"], OptionsInfo["AndCombineOperatorMode"] 339 ) 340 341 return (MolMatched, GroupsPatternMatchStatusList) 342 343 344 def DoesMolMeetSpecifiedMatchCriteria(GroupsPatternMolsMatchStatus, CombineMatchResults, AndCombineOperatorMode): 345 """Match molecule using specified match criteia.""" 346 347 if CombineMatchResults and AndCombineOperatorMode: 348 # Must match all specified SMARTS 349 Status = True 350 for MatchStatus in GroupsPatternMolsMatchStatus: 351 if not MatchStatus: 352 Status = False 353 break 354 else: 355 # One match is enough... 356 Status = False 357 for MatchStatus in GroupsPatternMolsMatchStatus: 358 if MatchStatus: 359 Status = True 360 break 361 362 return Status 363 364 365 def WriteMolecule( 366 Writer, GroupOutfilesWriters, Mol, Compute2DCoords, CombineMatchResults, GroupsPatternMatchStatusList 367 ): 368 """Write out molecule.""" 369 370 if OptionsInfo["CountMode"]: 371 return 372 373 if Compute2DCoords: 374 AllChem.Compute2DCoords(Mol) 375 376 if CombineMatchResults: 377 Writer.write(Mol) 378 else: 379 for GroupIndex in range(0, len(GroupsPatternMatchStatusList)): 380 if GroupsPatternMatchStatusList[GroupIndex]: 381 GroupOutfilesWriters[GroupIndex].write(Mol) 382 383 384 def SetupMoleculeWriters(): 385 """Set up molecule writers for output files.""" 386 387 Writer = None 388 GroupOutfilesWriters = [] 389 390 if OptionsInfo["CountMode"]: 391 return (Writer, GroupOutfilesWriters) 392 393 Outfile = OptionsInfo["Outfile"] 394 CombineMatchResults = OptionsInfo["CombineMatchResults"] 395 GroupsOutfiles = OptionsInfo["SpecifiedFunctionalGroupsOutfiles"] 396 397 if CombineMatchResults: 398 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"]) 399 if Writer is None: 400 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile) 401 MiscUtil.PrintInfo("Generating file %s..." % Outfile) 402 else: 403 for GroupOutfile in GroupsOutfiles: 404 GroupOutfileWriter = RDKitUtil.MoleculesWriter(GroupOutfile, **OptionsInfo["OutfileParams"]) 405 if GroupOutfileWriter is None: 406 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Writer) 407 GroupOutfilesWriters.append(GroupOutfileWriter) 408 409 GroupsCount = len(GroupsOutfiles) 410 if GroupsCount > 4: 411 MiscUtil.PrintInfo( 412 "Generating %d output files with the following file name format: %s<GroupName>.%s" 413 % (GroupsCount, OptionsInfo["OutfileBasename"], OptionsInfo["OutfileExt"]) 414 ) 415 else: 416 Delmiter = ", " 417 OutfileNames = Delmiter.join(GroupsOutfiles) 418 MiscUtil.PrintInfo("Generating %d output files: %s..." % (GroupsCount, OutfileNames)) 419 420 return (Writer, GroupOutfilesWriters) 421 422 423 def DoesPatternMolMatch(PatternMol, Mol, UseChirality, NegateMatch): 424 """Perform a substructure match for the presence of pattern molecule in a molecule.""" 425 426 MolMatched = Mol.HasSubstructMatch(PatternMol, useChirality=UseChirality) 427 if NegateMatch: 428 if MolMatched: 429 MolMatched = False 430 else: 431 MolMatched = True 432 433 return MolMatched 434 435 436 def ProcessFunctionalGroupsInfo(): 437 """Process functional groups information.""" 438 439 RetrieveFunctionalGroupsInfo() 440 ProcessSpecifiedFunctionalGroups() 441 442 SetupFunctionalGroupsOutputFileNames() 443 444 445 def ProcessSpecifiedFunctionalGroups(): 446 """Process and validate specified functional groups.""" 447 448 OptionsInfo["SpecifiedFunctionalGroups"] = [] 449 OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"] = [] 450 451 if re.match("^All$", OptionsInfo["FunctionalGroups"], re.I): 452 OptionsInfo["SpecifiedFunctionalGroups"] = FunctionalGroupsMap["Names"] 453 OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"] = [False] * len(OptionsInfo["SpecifiedFunctionalGroups"]) 454 return 455 456 # Set up a map of valid group names for checking specified group names... 457 CanonicalGroupNameMap = {} 458 for GroupName in FunctionalGroupsMap["Names"]: 459 CanonicalGroupNameMap[GroupName.lower()] = GroupName 460 461 # Parse and validate specified names... 462 GroupNames = re.sub(" ", "", OptionsInfo["FunctionalGroups"]) 463 if not GroupNames: 464 MiscUtil.PrintError('No functional group name specified for "-f, --functionalGroups" option') 465 466 SpecifiedFunctionalGroups = [] 467 SpecifiedNegateMatchStatus = [] 468 469 for GroupName in GroupNames.split(","): 470 CanonicalGroupName = GroupName.lower() 471 NegateMatchStatus = False 472 if re.match("^!", CanonicalGroupName, re.I): 473 NegateMatchStatus = True 474 CanonicalGroupName = re.sub("^!", "", CanonicalGroupName) 475 if CanonicalGroupName in CanonicalGroupNameMap: 476 SpecifiedFunctionalGroups.append(CanonicalGroupNameMap[CanonicalGroupName]) 477 SpecifiedNegateMatchStatus.append(NegateMatchStatus) 478 else: 479 MiscUtil.PrintWarning( 480 'The functional group name, %s, specified using "-f, --functionalGroups" option is not a valid name.' 481 % (GroupName) 482 ) 483 484 if not len(SpecifiedFunctionalGroups): 485 MiscUtil.PrintError('No valid functional group names specified for "-f, --functionalGroups" option') 486 487 OptionsInfo["SpecifiedFunctionalGroups"] = SpecifiedFunctionalGroups 488 OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"] = SpecifiedNegateMatchStatus 489 490 491 def SetupFunctionalGroupsSMARTSPatterns(): 492 """Setup SMARTS patterns for specified functional groups.""" 493 494 OptionsInfo["SpecifiedFunctionalGroupsSMARTSPatterns"] = [] 495 FunctionalGroupsPatternMols = [] 496 497 for Name in OptionsInfo["SpecifiedFunctionalGroups"]: 498 SMARTSPattern = FunctionalGroupsMap["SMARTSPattern"][Name] 499 PatternMol = Chem.MolFromSmarts(SMARTSPattern) 500 if PatternMol is None: 501 MiscUtil.PrintError("Failed to parse SMARTS pattern, %s, for function group, %s" % (SMARTSPattern, Name)) 502 503 OptionsInfo["SpecifiedFunctionalGroupsSMARTSPatterns"].append(SMARTSPattern) 504 FunctionalGroupsPatternMols.append(PatternMol) 505 506 return FunctionalGroupsPatternMols 507 508 509 def SetupFunctionalGroupsOutputFileNames(): 510 """Setup output file names for specified functional group names.""" 511 512 OptionsInfo["SpecifiedFunctionalGroupsOutfiles"] = [] 513 514 if OptionsInfo["CountMode"]: 515 # No need of any output file... 516 return 517 518 if OptionsInfo["CombineMatchResults"]: 519 # No need of output files for specified functional groups... 520 return 521 522 OutfileBasename = OptionsInfo["OutfileBasename"] 523 OutfileExt = OptionsInfo["OutfileExt"] 524 SpecifiedFunctionalGroupsOutfiles = [] 525 526 GroupsCount = len(OptionsInfo["SpecifiedFunctionalGroups"]) 527 for GroupIndex in range(0, GroupsCount): 528 GroupName = OptionsInfo["SpecifiedFunctionalGroups"][GroupIndex] 529 if OptionsInfo["SpecifiedFunctionalGroupsNegateMatch"][GroupIndex]: 530 GroupName = "Not" + GroupName 531 GroupName = re.sub(r"\.", "", GroupName) 532 533 GroupOutfile = "%s%s.%s" % (OutfileBasename, GroupName, OutfileExt) 534 SpecifiedFunctionalGroupsOutfiles.append(GroupOutfile) 535 536 OptionsInfo["SpecifiedFunctionalGroupsOutfiles"] = SpecifiedFunctionalGroupsOutfiles 537 538 539 def RetrieveFunctionalGroupsInfo(): 540 """Retrieve functional groups information.""" 541 542 MiscUtil.PrintInfo( 543 "\nRetrieving data from default RDKit functional groups hierarchy file Functional_Group_Hierarchy.txt..." 544 ) 545 546 FunctionalGroupNamesFile = OptionsInfo["GroupNamesFile"] 547 FunctionalGroupsNodes = FunctionalGroups.BuildFuncGroupHierarchy(FunctionalGroupNamesFile) 548 549 FunctionalGroupsMap["Names"] = [] 550 FunctionalGroupsMap["SMARTSPattern"] = {} 551 552 RetrieveDataFromFunctionalGroupsHierarchy(FunctionalGroupsNodes) 553 554 if not len(FunctionalGroupsMap["Names"]): 555 MiscUtil.PrintError("Failed to retrieve any functional group names and SMARTS patterns...") 556 557 MiscUtil.PrintInfo( 558 "Total number of functional groups present functional group hierarchy: %d" % (len(FunctionalGroupsMap["Names"])) 559 ) 560 561 562 def RetrieveDataFromFunctionalGroupsHierarchy(FGNodes): 563 """Retrieve functional groups data by recursively visiting functional group nodes.""" 564 565 for FGNode in FGNodes: 566 Name = FGNode.label 567 SMARTSPattern = FGNode.smarts 568 569 if Name in FunctionalGroupsMap["SMARTSPattern"]: 570 MiscUtil.PrintWarning("Ignoring duplicate functional group name: %s..." % Name) 571 else: 572 FunctionalGroupsMap["Names"].append(Name) 573 FunctionalGroupsMap["SMARTSPattern"][Name] = SMARTSPattern 574 575 RetrieveDataFromFunctionalGroupsHierarchy(FGNode.children) 576 577 578 def ListFunctionalGroupsInfo(): 579 """List functional groups information.""" 580 581 MiscUtil.PrintInfo("\nListing available functional groups names and SMARTS patterns...") 582 MiscUtil.PrintInfo("\nFunctionalGroupName\tSMARTSPattern") 583 584 for Name in sorted(FunctionalGroupsMap["Names"]): 585 SMARTSPattern = FunctionalGroupsMap["SMARTSPattern"][Name] 586 MiscUtil.PrintInfo("%s\t%s" % (Name, SMARTSPattern)) 587 588 MiscUtil.PrintInfo("") 589 590 591 def ProcessOptions(): 592 """Process and validate command line arguments and options.""" 593 594 MiscUtil.PrintInfo("Processing options...") 595 596 # Validate options... 597 ValidateOptions() 598 599 OptionsInfo["CombineMatches"] = Options["--combineMatches"] 600 601 OptionsInfo["CombineMatchResults"] = True 602 if re.match("^No$", Options["--combineMatches"], re.I): 603 OptionsInfo["CombineMatchResults"] = False 604 if Options["--outfile"]: 605 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"]) 606 OptionsInfo["OutfileBasename"] = FileName 607 OptionsInfo["OutfileExt"] = FileExt 608 609 OptionsInfo["CombineOperator"] = Options["--combineOperator"] 610 OptionsInfo["AndCombineOperatorMode"] = True 611 if re.match("^or$", Options["--combineOperator"], re.I): 612 OptionsInfo["AndCombineOperatorMode"] = False 613 614 OptionsInfo["GroupNamesFile"] = None 615 if not re.match("^auto$", Options["--groupNamesFile"], re.I): 616 OptionsInfo["GroupNamesFile"] = Options["--groupNamesFile"] 617 618 OptionsInfo["FunctionalGroups"] = Options["--functionalGroups"] 619 620 OptionsInfo["Infile"] = Options["--infile"] 621 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters( 622 "--infileParams", Options["--infileParams"], Options["--infile"] 623 ) 624 625 OptionsInfo["Outfile"] = Options["--outfile"] 626 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters( 627 "--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"] 628 ) 629 630 OptionsInfo["Overwrite"] = Options["--overwrite"] 631 632 OptionsInfo["CountMode"] = False 633 if re.match("^count$", Options["--mode"], re.I): 634 OptionsInfo["CountMode"] = True 635 636 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False 637 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"]) 638 639 OptionsInfo["UseChirality"] = False 640 if re.match("^yes$", Options["--useChirality"], re.I): 641 OptionsInfo["UseChirality"] = True 642 643 644 def RetrieveOptions(): 645 """Retrieve command line arguments and options.""" 646 647 # Get options... 648 global Options 649 Options = docopt(_docoptUsage_) 650 651 # Set current working directory to the specified directory... 652 WorkingDir = Options["--workingdir"] 653 if WorkingDir: 654 os.chdir(WorkingDir) 655 656 # Handle examples option... 657 if "--examples" in Options and Options["--examples"]: 658 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 659 sys.exit(0) 660 661 # Handle listing of functional group information... 662 if Options and Options["--list"]: 663 ProcessListFunctionalGroupsOption() 664 sys.exit(0) 665 666 667 def ProcessListFunctionalGroupsOption(): 668 """Process list functional groups information.""" 669 670 # Validate and process dataFile option for listing functional groups information... 671 OptionsInfo["GroupNamesFile"] = None 672 if not re.match("^auto$", Options["--groupNamesFile"], re.I): 673 MiscUtil.ValidateOptionFilePath("-g, --groupNamesFile", Options["--groupNamesFile"]) 674 OptionsInfo["GroupNamesFile"] = Options["--groupNamesFile"] 675 676 RetrieveFunctionalGroupsInfo() 677 ListFunctionalGroupsInfo() 678 679 680 def ValidateOptions(): 681 """Validate option values.""" 682 683 MiscUtil.ValidateOptionTextValue("-c, --combineMatches", Options["--combineMatches"], "yes no") 684 MiscUtil.ValidateOptionTextValue("--combineOperator", Options["--combineOperator"], "and or") 685 686 if not re.match("^auto$", Options["--groupNamesFile"], re.I): 687 MiscUtil.ValidateOptionFilePath("-g, groupNamesFile", Options["--groupNamesFile"]) 688 689 if re.match("^none$", Options["--functionalGroups"], re.I): 690 MiscUtil.PrintError('The name(s) of functional groups must be specified using "-f, --functionalGroups" option') 691 692 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"]) 693 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd smi txt csv tsv") 694 if Options["--outfile"]: 695 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd smi") 696 MiscUtil.ValidateOptionsOutputFileOverwrite( 697 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 698 ) 699 MiscUtil.ValidateOptionsDistinctFileNames( 700 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"] 701 ) 702 703 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "retrieve count") 704 if re.match("^retrieve$", Options["--mode"], re.I): 705 if not Options["--outfile"]: 706 MiscUtil.PrintError( 707 'The outfile must be specified using "-o, --outfile" during "retrieve" value of "-m, --mode" option' 708 ) 709 710 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no") 711 712 MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no") 713 714 715 # Setup a usage string for docopt... 716 _docoptUsage_ = """ 717 RDKitSearchFunctionalGroups.py - Search for functional groups using SMARTS patterns 718 719 Usage: 720 RDKitSearchFunctionalGroups.py [--combineMatches <yes or no>] [--combineOperator <and or or>] 721 [--groupNamesFile <FileName or auto>] [--infileParams <Name,Value,...>] 722 [--mode <retrieve or count>] [--mp <yes or no>] [--mpParams <Name,Value,...>] 723 [--negate <yes or no>] [--outfileParams <Name,Value,...>] [--overwrite] 724 [--useChirality <yes or no>] [-w <dir>] [-o <outfile>] -i <infile> -f <Name1,Name2,Name3... or All> 725 RDKitSearchFunctionalGroups.py [--groupNamesFile <FileName or auto>] -l | --list 726 RDKitSearchFunctionalGroups.py -h | --help | -e | --examples 727 728 Description: 729 Perform a substructure search in an input file using SMARTS patterns for functional 730 groups and write out the matched molecules to an output file or simply count the 731 number of matches. 732 733 The SMARTS patterns for specified functional group(s) are retrieved from file, 734 Functional_Group_Hierarchy.txt, available in RDKit data directory. 735 736 The names of valid functional groups and hierarchies are dynamically retrieved from the 737 functional groups hierarchy file and are shown below: 738 739 AcidChloride, AcidChloride.Aromatic, AcidChloride.Aliphatic 740 Alcohol, Alcohol.Aromatic, Alcohol.Aliphatic 741 Aldehyde, Aldehyde.Aromatic, Aldehyde.Aliphatic 742 Amine, Amine.Primary, Amine.Primary.Aromatic, Amine.Primary.Aliphatic, 743 Amine.Secondary, Amine.Secondary.Aromatic, Amine.Secondary.Aliphatic 744 Amine.Tertiary, Amine.Tertiary.Aromatic, Amine.Tertiary.Aliphatic 745 Amine.Aromatic, Amine.Aliphatic, Amine.Cyclic 746 Azide, Azide.Aromatic, Azide.Aliphatic 747 BoronicAcid, BoronicAcid.Aromatic, BoronicAcid.Aliphatic 748 CarboxylicAcid, CarboxylicAcid.Aromatic, CarboxylicAcid.Aliphatic, 749 CarboxylicAcid.AlphaAmino 750 Halogen, Halogen.Aromatic, Halogen.Aliphatic 751 Halogen.NotFluorine, Halogen.NotFluorine.Aliphatic, 752 Halogen.NotFluorine.Aromatic 753 Halogen.Bromine, Halogen.Bromine.Aliphatic, Halogen.Bromine.Aromatic, 754 Halogen.Bromine.BromoKetone 755 Isocyanate, Isocyanate.Aromatic, Isocyanate.Aliphatic 756 Nitro, Nitro.Aromatic, Nitro.Aliphatic, 757 SulfonylChloride, SulfonylChloride.Aromatic, SulfonylChloride.Aliphatic 758 TerminalAlkyne 759 760 The supported input file formats are: SD (.sdf, .sd), SMILES (.smi, .csv, .tsv, .txt) 761 762 The supported output file formats are: SD (.sdf, .sd), SMILES (.smi) 763 764 Options: 765 -c, --combineMatches <yes or no> [default: yes] 766 Combine search results for matching SMARTS patterns of specified functional groups 767 against a molecule. Possible values: yes or no. 768 769 The matched molecules are written to a single output file for "yes" value. Otherwise, 770 multiple output files are generated, one for each functional group. The names of 771 these files correspond to a combination of the basename of the specified output file 772 and the name of the functional group. 773 774 No output files are generated during "count" value of "-m, --mode" option. 775 --combineOperator <and or or> [default: and] 776 Logical operator to use for combining match results corresponding to specified 777 functional group names before writing out a single file. This option is ignored 778 during "No" value of "-c, --combineMatches" option. 779 -e, --examples 780 Print examples. 781 -g, --groupNamesFile <FileName or auto> [default: auto] 782 Specify a file name containing data for functional groups hierarchy or use functional 783 group hierarchy file, Functional_Group_Hierarchy.txt, available in RDKit data directory. 784 785 RDKit data format: Name<tab>Smarts<tab>Label<tab>RemovalReaction (optional) 786 787 The format of data in local functional group hierarchy must match format of the 788 data in functional group file available in RDKit data directory. 789 -f, --functionalGroups <Name1,Name2,Name3... or All> [default: none] 790 Functional group names for performing substructure SMARTS search. Possible values: 791 Comma delimited list of valid functional group names or All. The current set of valid 792 functional group names are listed in the description section. 793 794 The match results for multiple functional group names are combined using 'and' 795 operator before writing them out to single file. No merging of match results takes 796 place during generation of individual result files corresponding to fictional group 797 names. 798 799 The functional group name may be started with an exclamation mark to negate 800 the match result for that fictional group. 801 -h, --help 802 Print this help message. 803 -i, --infile <infile> 804 Input file name. 805 --infileParams <Name,Value,...> [default: auto] 806 A comma delimited list of parameter name and value pairs for reading 807 molecules from files. The supported parameter names for different file 808 formats, along with their default values, are shown below: 809 810 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes 811 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space, 812 smilesTitleLine,auto,sanitize,yes 813 814 Possible values for smilesDelimiter: space, comma or tab. 815 -l, --list 816 List functional groups information without performing any search. 817 -m, --mode <retrieve or count> [default: retrieve] 818 Specify whether to retrieve and write out matched molecules to an output 819 file or simply count the number of matches. 820 --mp <yes or no> [default: no] 821 Use multiprocessing. 822 823 By default, input data is retrieved in a lazy manner via mp.Pool.imap() 824 function employing lazy RDKit data iterable. This allows processing of 825 arbitrary large data sets without any additional requirements memory. 826 827 All input data may be optionally loaded into memory by mp.Pool.map() 828 before starting worker processes in a process pool by setting the value 829 of 'inputDataMode' to 'InMemory' in '--mpParams' option. 830 831 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input 832 data mode may adversely impact the performance. The '--mpParams' section 833 provides additional information to tune the value of 'chunkSize'. 834 --mpParams <Name,Value,...> [default: auto] 835 A comma delimited list of parameter name and value pairs to configure 836 multiprocessing. 837 838 The supported parameter names along with their default and possible 839 values are shown below: 840 841 chunkSize, auto 842 inputDataMode, Lazy [ Possible values: InMemory or Lazy ] 843 numProcesses, auto [ Default: mp.cpu_count() ] 844 845 These parameters are used by the following functions to configure and 846 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and 847 mp.Pool.imap(). 848 849 The chunkSize determines chunks of input data passed to each worker 850 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions. 851 The default value of chunkSize is dependent on the value of 'inputDataMode'. 852 853 The mp.Pool.map() function, invoked during 'InMemory' input data mode, 854 automatically converts RDKit data iterable into a list, loads all data into 855 memory, and calculates the default chunkSize using the following method 856 as shown in its code: 857 858 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4) 859 if extra: chunkSize += 1 860 861 For example, the default chunkSize will be 7 for a pool of 4 worker processes 862 and 100 data items. 863 864 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs 865 'lazy' RDKit data iterable to retrieve data as needed, without loading all the 866 data into memory. Consequently, the size of input data is not known a priori. 867 It's not possible to estimate an optimal value for the chunkSize. The default 868 chunkSize is set to 1. 869 870 The default value for the chunkSize during 'Lazy' data mode may adversely 871 impact the performance due to the overhead associated with exchanging 872 small chunks of data. It is generally a good idea to explicitly set chunkSize to 873 a larger value during 'Lazy' input data mode, based on the size of your input 874 data and number of processes in the process pool. 875 876 The mp.Pool.map() function waits for all worker processes to process all 877 the data and return the results. The mp.Pool.imap() function, however, 878 returns the the results obtained from worker processes as soon as the 879 results become available for specified chunks of data. 880 881 The order of data in the results returned by both mp.Pool.map() and 882 mp.Pool.imap() functions always corresponds to the input data. 883 -o, --outfile <outfile> 884 Output file name. 885 --outfileParams <Name,Value,...> [default: auto] 886 A comma delimited list of parameter name and value pairs for writing 887 molecules to files. The supported parameter names for different file 888 formats, along with their default values, are shown below: 889 890 SD: compute2DCoords,auto,kekulize,yes,forceV3000,no 891 SMILES: smilesKekulize,no,smilesDelimiter,space, smilesIsomeric,yes, 892 smilesTitleLine,yes,smilesMolName,yes,smilesMolProps,no 893 894 Default value for compute2DCoords: yes for SMILES input file; no for all other 895 file types. 896 --overwrite 897 Overwrite existing files. 898 -u, --useChirality <yes or no> [default: no] 899 Use stereochemistry information for SMARTS search. 900 -w, --workingdir <dir> 901 Location of working directory which defaults to the current directory. 902 903 Examples: 904 To list names of all available functional groups along with their SMARTS 905 patterns, type: 906 907 % RDKitSearchFunctionalGroups.py -l 908 909 To retrieve molecules containing amine functional group and write out a 910 SMILES file, type: 911 912 % RDKitSearchFunctionalGroups.py -f Amine -i Sample.smi -o SampleOut.smi 913 914 To retrieve molecules containing amine functional group, perform search in 915 multiprocessing mode on all available CPUs without loading all data into 916 memory, and write out a SMILES file, type: 917 918 % RDKitSearchFunctionalGroups.py --mp yes -f Amine -i Sample.smi 919 -o SampleOut.smi 920 921 To retrieve molecules containing amine functional group, perform search in 922 multiprocessing mode on all available CPUs by loading all data into memory, 923 and write out a SMILES file, type: 924 925 % RDKitSearchFunctionalGroups.py --mp yes --mpParams "inputDataMode, 926 InMemory" -f Amine -i Sample.smi -o SampleOut.smi 927 928 To retrieve molecules containing amine functional group, perform search in 929 multiprocessing mode on specific number of CPUs and chunksize without loading 930 all data into memory, and write out a SMILES file, type: 931 932 % RDKitSearchFunctionalGroups.py --mp yes --mpParams "inputDataMode, 933 lazy,numProcesses,4,chunkSize,8" -f Amine -i Sample.smi -o 934 SampleOut.smi 935 936 To retrieve molecules containing amine functional group but not halogens and carboxylic 937 acid functional groups and write out a SMILES file, type: 938 939 % RDKitSearchFunctionalGroups.py -f 'Amine,!Halogen,!CarboxylicAcid' 940 -i Sample.smi -o SampleOut.smi 941 942 To retrieve molecules containing amine, halogens or carboxylic acid functional groups 943 and write out a SMILES file, type: 944 945 % RDKitSearchFunctionalGroups.py -f 'Amine,Halogen,CarboxylicAcid' 946 --combineOperator or -i Sample.smi -o SampleOut.smi 947 948 To retrieve molecules containing amine and carboxylic acid functional groups defined in 949 a local functional groups hierarchy file and write out individual SD files for each 950 funcitonal group, type: 951 952 % RDKitSearchFunctionalGroups.py -f 'Amine,CarboxylicAcid' -i Sample.sdf 953 -g Custom_Functional_Group_Hierarchy.txt --combineMatches No -o SampleOut.sdf 954 955 To count number of all functional groups in molecules without writing out an output 956 files, type: 957 958 % RDKitSearchFunctionalGroups.py -m count -f All --combineMatches no -i Sample.smi 959 960 To retrieve molecule not containing aromatic alcohol and aromatic halogen functional 961 group along with the use of chirality during substructure search and write out individual 962 SMILES files for each functional group, type: 963 964 % RDKitSearchFunctionalGroups.py --combineMatches no -u yes 965 -f '!Alcohol.Aromatic,!Halogen.Aromatic' -i Sample.smi -o SampleOut.smi 966 967 To retrieve molecule containing amine functional group from a CSV SMILES file, 968 SMILES strings in column 1, name in column 2, and write out a SD file, type: 969 970 % RDKitSearchFunctionalGroups.py -f Amine --infileParams 971 "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1, 972 smilesNameColumn,2" --outfileParams "compute2DCoords,yes" 973 -i SampleSMILES.csv -o SampleOut.sdf 974 975 Author: 976 Manish Sud(msud@san.rr.com) 977 978 See also: 979 RDKitConvertFileFormat.py, RDKitFilterPAINS.py, RDKitSearchSMARTS.py 980 981 Copyright: 982 Copyright (C) 2026 Manish Sud. All rights reserved. 983 984 The functionality available in this script is implemented using RDKit, an 985 open source toolkit for cheminformatics developed by Greg Landrum. 986 987 This file is part of MayaChemTools. 988 989 MayaChemTools is free software; you can redistribute it and/or modify it under 990 the terms of the GNU Lesser General Public License as published by the Free 991 Software Foundation; either version 3 of the License, or (at your option) any 992 later version. 993 994 """ 995 996 if __name__ == "__main__": 997 main()
