1 #!/bin/env python 2 # 3 # File: RDKitCalculateMolecularDescriptors.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 import numpy 37 38 # RDKit imports... 39 try: 40 from rdkit import rdBase 41 from rdkit import Chem 42 from rdkit.Chem import AllChem 43 from rdkit.Chem import rdMolDescriptors 44 from rdkit.Chem import Descriptors 45 from rdkit.Chem import Descriptors3D 46 except ImportError as ErrMsg: 47 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg) 48 sys.stderr.write("Check/update your RDKit environment and try again.\n\n") 49 sys.exit(1) 50 51 # MayaChemTools imports... 52 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python")) 53 try: 54 from docopt import docopt 55 import MiscUtil 56 import RDKitUtil 57 except ImportError as ErrMsg: 58 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg) 59 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n") 60 sys.exit(1) 61 62 ScriptName = os.path.basename(sys.argv[0]) 63 Options = {} 64 OptionsInfo = {} 65 66 DescriptorNamesMap = {} 67 68 69 def main(): 70 """Start execution of the script.""" 71 72 MiscUtil.PrintInfo( 73 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 74 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()) 75 ) 76 77 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 78 79 # Retrieve command line arguments and options... 80 RetrieveOptions() 81 82 # Process and validate command line arguments and options... 83 ProcessOptions() 84 85 # Perform actions required by the script... 86 CalculateMolecularDescriptors() 87 88 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 89 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 90 91 92 def CalculateMolecularDescriptors(): 93 """Calculate molecular descriptors.""" 94 95 ProcessMolecularDescriptorsInfo() 96 PerformCalculations() 97 98 99 def ProcessMolecularDescriptorsInfo(): 100 """Process descriptors information.""" 101 102 RetrieveMolecularDescriptorsInfo() 103 ProcessSpecifiedDescriptorNames() 104 105 106 def PerformCalculations(): 107 """Calculate descriptors for a specified list of descriptors.""" 108 109 # Setup a molecule reader... 110 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"]) 111 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"]) 112 113 # Setup a writer... 114 Writer = SetupMoleculeWriter() 115 116 # Process molecules... 117 MolCount, ValidMolCount = ProcessMolecules(Mols, Writer) 118 119 if Writer is not None: 120 Writer.close() 121 122 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount) 123 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount) 124 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount)) 125 126 127 def ProcessMolecules(Mols, Writer): 128 """Process and filter molecules.""" 129 130 if OptionsInfo["MPMode"]: 131 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer) 132 else: 133 return ProcessMoleculesUsingSingleProcess(Mols, Writer) 134 135 136 def ProcessMoleculesUsingSingleProcess(Mols, Writer): 137 """Process molecules and calculate descriptors using a single process.""" 138 139 DescriptorsCount = len(OptionsInfo["SpecifiedDescriptorNames"]) 140 MiscUtil.PrintInfo( 141 "\nCalculating %d molecular %s for each molecule..." 142 % (DescriptorsCount, ("descroptors" if DescriptorsCount > 1 else "descriptor")) 143 ) 144 145 (MolCount, ValidMolCount) = [0] * 2 146 for MolIndex, Mol in enumerate(Mols): 147 MolCount += 1 148 if Mol is None: 149 continue 150 151 if RDKitUtil.IsMolEmpty(Mol): 152 MolName = RDKitUtil.GetMolName(Mol, MolCount) 153 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 154 continue 155 ValidMolCount += 1 156 157 # Calculate and write descriptor values... 158 CalculatedValues = CalculateDescriptorValues(MolIndex, Mol) 159 WriteDescriptorValues(Mol, MolCount, Writer, CalculatedValues) 160 161 return (MolCount, ValidMolCount) 162 163 164 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer): 165 """Process molecules and calculate descriptors using multiprocessing.""" 166 167 DescriptorsCount = len(OptionsInfo["SpecifiedDescriptorNames"]) 168 MiscUtil.PrintInfo( 169 "\nCalculating %d molecular %s for each molecule using multiprocessing......" 170 % (DescriptorsCount, ("descroptors" if DescriptorsCount > 1 else "descriptor")) 171 ) 172 173 MPParams = OptionsInfo["MPParams"] 174 175 # Setup data for initializing a worker process... 176 MiscUtil.PrintInfo("Encoding options info...") 177 InitializeWorkerProcessArgs = ( 178 MiscUtil.ObjectToBase64EncodedString(Options), 179 MiscUtil.ObjectToBase64EncodedString(OptionsInfo), 180 ) 181 182 # Setup a encoded mols data iterable for a worker process... 183 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols) 184 185 # Setup process pool along with data initialization for each process... 186 MiscUtil.PrintInfo( 187 "\nConfiguring multiprocessing using %s method..." 188 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()") 189 ) 190 MiscUtil.PrintInfo( 191 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" 192 % ( 193 MPParams["NumProcesses"], 194 MPParams["InputDataMode"], 195 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]), 196 ) 197 ) 198 199 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs) 200 201 # Start processing... 202 if re.match("^Lazy$", MPParams["InputDataMode"], re.I): 203 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 204 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I): 205 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 206 else: 207 MiscUtil.PrintError( 208 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"]) 209 ) 210 211 (MolCount, ValidMolCount) = [0] * 2 212 for Result in Results: 213 MolCount += 1 214 MolIndex, EncodedMol, CalculatedValues = Result 215 216 if EncodedMol is None: 217 continue 218 ValidMolCount += 1 219 220 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 221 222 # Write descriptor values... 223 WriteDescriptorValues(Mol, MolCount, Writer, CalculatedValues) 224 225 return (MolCount, ValidMolCount) 226 227 228 def InitializeWorkerProcess(*EncodedArgs): 229 """Initialize data for a worker process.""" 230 231 global Options, OptionsInfo 232 233 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid()) 234 235 # Decode Options and OptionInfo... 236 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0]) 237 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1]) 238 239 RetrieveMolecularDescriptorsInfo(PrintInfo=False) 240 241 242 def WorkerProcess(EncodedMolInfo): 243 """Process data for a worker process.""" 244 245 MolIndex, EncodedMol = EncodedMolInfo 246 247 CalculatedValues = [] 248 if EncodedMol is None: 249 return [MolIndex, None, CalculatedValues] 250 251 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 252 if RDKitUtil.IsMolEmpty(Mol): 253 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1)) 254 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 255 return [MolIndex, None, CalculatedValues] 256 257 return [MolIndex, EncodedMol, CalculateDescriptorValues(MolIndex, Mol)] 258 259 260 def CalculateDescriptorValues(MolIndex, Mol): 261 """Calculate descriptor values.""" 262 263 return [ 264 FormatCalculatedValue(CalculateDescriptorValue(MolIndex, Mol, Name), OptionsInfo["Precision"]) 265 for Name in OptionsInfo["SpecifiedDescriptorNames"] 266 ] 267 268 269 def CalculateDescriptorValue(MolIndex, Mol, Name): 270 """Calculate value for a specific descriptor along with handling any calculation failure.""" 271 272 try: 273 Value = DescriptorNamesMap["ComputeFunction"][Name](Mol) 274 except ValueError as ErrMsg: 275 MiscUtil.PrintWarning( 276 "Failed to calculate descriptor %s for molecule %s:\n%s\n" % (Name, (MolIndex + 1), ErrMsg) 277 ) 278 Value = "NA" 279 280 return Value 281 282 283 def WriteDescriptorValues(Mol, MolNum, Writer, CalculatedValues): 284 """Write out calculated descriptor values.""" 285 286 if OptionsInfo["TextOutFileMode"]: 287 LineWords = [] 288 if OptionsInfo["SMILESOut"]: 289 SMILES = Chem.MolToSmiles(Mol, isomericSmiles=True, canonical=True) 290 LineWords.append(SMILES) 291 292 MolName = RDKitUtil.GetMolName(Mol, MolNum) 293 LineWords.append(MolName) 294 LineWords.extend(CalculatedValues) 295 Line = OptionsInfo["TextOutFileDelim"].join(LineWords) 296 Writer.write("%s\n" % Line) 297 else: 298 for Index, Name in enumerate(OptionsInfo["SpecifiedDescriptorNames"]): 299 Mol.SetProp(Name, CalculatedValues[Index]) 300 301 if OptionsInfo["OutfileParams"]["Compute2DCoords"]: 302 AllChem.Compute2DCoords(Mol) 303 304 Writer.write(Mol) 305 306 307 def SetupMoleculeWriter(): 308 """Setup a molecule writer.""" 309 310 if OptionsInfo["TextOutFileMode"]: 311 Writer = open(OptionsInfo["Outfile"], "w") 312 else: 313 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"]) 314 315 if Writer is None: 316 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"]) 317 318 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"]) 319 320 # Wite out headers for a text file... 321 if OptionsInfo["TextOutFileMode"]: 322 LineWords = [] 323 if OptionsInfo["SMILESOut"]: 324 LineWords.append("SMILES") 325 LineWords.append("MolID") 326 LineWords.extend(OptionsInfo["SpecifiedDescriptorNames"]) 327 Line = OptionsInfo["TextOutFileDelim"].join(LineWords) 328 Writer.write("%s\n" % Line) 329 330 return Writer 331 332 333 def FormatCalculatedValue(Value, Precision): 334 """Format calculated value of descriptor based on its type.""" 335 336 if (type(Value) is float) or (type(Value) is numpy.float64): 337 FormattedValue = "%.*f" % (Precision, Value) 338 if not re.search("[1-9]", FormattedValue): 339 FormattedValue = "0.0" 340 elif type(Value) is list: 341 FormattedValue = "%s" % Value 342 FormattedValue = re.sub(r"\[|\]|,", "", FormattedValue) 343 else: 344 FormattedValue = "%s" % Value 345 346 return FormattedValue 347 348 349 def ProcessSpecifiedDescriptorNames(): 350 """Process and validate specified decriptor names.""" 351 352 OptionsInfo["SpecifiedDescriptorNames"] = [] 353 354 if not re.match("^(2D|3D|All|FragmentCountOnly|Specify)$", OptionsInfo["Mode"], re.I): 355 MiscUtil.PrintError('Mode value, %s, using "-m, --mode" option is not a valid value.' % OptionsInfo["Mode"]) 356 357 if re.match("^2D$", OptionsInfo["Mode"], re.I): 358 OptionsInfo["SpecifiedDescriptorNames"] = DescriptorNamesMap["2D"]["Names"] 359 if OptionsInfo["FragmentCount"]: 360 OptionsInfo["SpecifiedDescriptorNames"].extend(DescriptorNamesMap["FragmentCount"]["Names"]) 361 return 362 elif re.match("^3D$", OptionsInfo["Mode"], re.I): 363 OptionsInfo["SpecifiedDescriptorNames"] = DescriptorNamesMap["3D"]["Names"] 364 return 365 elif re.match("^All$", OptionsInfo["Mode"], re.I): 366 OptionsInfo["SpecifiedDescriptorNames"] = DescriptorNamesMap["2D"]["Names"] 367 if OptionsInfo["FragmentCount"]: 368 OptionsInfo["SpecifiedDescriptorNames"].extend(DescriptorNamesMap["FragmentCount"]["Names"]) 369 OptionsInfo["SpecifiedDescriptorNames"].extend(DescriptorNamesMap["3D"]["Names"]) 370 return 371 elif re.match("^FragmentCountOnly$", OptionsInfo["Mode"], re.I): 372 OptionsInfo["SpecifiedDescriptorNames"] = DescriptorNamesMap["FragmentCount"]["Names"] 373 return 374 375 # Set up a canonical descriptor names map for checking specified names... 376 CanonicalNameMap = {} 377 for Name in DescriptorNamesMap["ComputeFunction"]: 378 CanonicalNameMap[Name.lower()] = Name 379 380 # Parse and validate specified names... 381 DescriptorNames = re.sub(" ", "", OptionsInfo["DescriptorNames"]) 382 if not DescriptorNames: 383 MiscUtil.PrintError('No descriptor names specified for "-d, --descriptorNames" option') 384 385 SMILESInfile = MiscUtil.CheckFileExt(Options["--infile"], "smi") 386 Canonical3DNameMap = {} 387 if SMILESInfile: 388 for Name in DescriptorNamesMap["3D"]["Names"]: 389 Canonical3DNameMap[Name.lower()] = Name 390 391 SpecifiedDescriptorNames = [] 392 for Name in DescriptorNames.split(","): 393 CanonicalName = Name.lower() 394 if CanonicalName in CanonicalNameMap: 395 SpecifiedDescriptorNames.append(CanonicalNameMap[CanonicalName]) 396 else: 397 MiscUtil.PrintError( 398 'The descriptor name, %s, specified using "-d, --descriptorNames" option is not a valid name.' % (Name) 399 ) 400 if SMILESInfile: 401 if CanonicalName in Canonical3DNameMap: 402 MiscUtil.PrintError( 403 'The 3D descriptor name, %s, specified using "-d, --descriptorNames" option is not a valid for SMILES input file.' 404 % (Name) 405 ) 406 407 if not len(SpecifiedDescriptorNames): 408 MiscUtil.PrintError('No valid descriptor name specified for "-d, --descriptorNames" option') 409 410 OptionsInfo["SpecifiedDescriptorNames"] = SpecifiedDescriptorNames 411 412 413 def RetrieveMolecularDescriptorsInfo(PrintInfo=True): 414 """Retrieve descriptors information.""" 415 416 if PrintInfo: 417 MiscUtil.PrintInfo("\nRetrieving information for avalible molecular descriptors...") 418 419 # Initialze data for 2D, FragmentCount and 3D descriptors... 420 DescriptorNamesMap["Types"] = ["2D", "FragmentCount", "3D"] 421 DescriptorNamesMap["ComputeFunction"] = {} 422 423 Autocorr2DExclude = OptionsInfo["Autocorr2DExclude"] if "Autocorr2DExclude" in OptionsInfo else False 424 425 for Type in DescriptorNamesMap["Types"]: 426 DescriptorNamesMap[Type] = {} 427 DescriptorNamesMap[Type]["Names"] = [] 428 429 # Setup data for 2D and FragmentCount... 430 DescriptorsInfo = Descriptors.descList 431 for DescriptorInfo in DescriptorsInfo: 432 Name = DescriptorInfo[0] 433 ComputeFunction = DescriptorInfo[1] 434 435 Type = "2D" 436 if re.match("^fr_", Name, re.I): 437 Type = "FragmentCount" 438 elif re.match("^Autocorr2D$", Name, re.I): 439 if Autocorr2DExclude: 440 continue 441 442 if Name in DescriptorNamesMap["ComputeFunction"]: 443 if PrintInfo: 444 MiscUtil.PrintWarning("Ignoring duplicate descriptor name: %s..." % Name) 445 else: 446 DescriptorNamesMap[Type]["Names"].append(Name) 447 DescriptorNamesMap["ComputeFunction"][Name] = ComputeFunction 448 449 # Add new 2D decriptor name to the list... 450 Type = "2D" 451 if not Autocorr2DExclude: 452 try: 453 Name = "Autocorr2D" 454 ComputeFunction = rdMolDescriptors.CalcAUTOCORR2D 455 if Name not in DescriptorNamesMap["ComputeFunction"]: 456 DescriptorNamesMap[Type]["Names"].append(Name) 457 DescriptorNamesMap["ComputeFunction"][Name] = ComputeFunction 458 except AttributeError: 459 if PrintInfo: 460 MiscUtil.PrintInfo("2D descriptor, %s, is not available in your current version of RDKit." % Name) 461 462 # Set data for 3D descriptors... 463 Type = "3D" 464 NameToComputeFunctionMap = { 465 "PMI1": Descriptors3D.PMI1, 466 "PMI2": Descriptors3D.PMI2, 467 "PMI3": Descriptors3D.PMI3, 468 "NPR1": Descriptors3D.NPR1, 469 "NPR2": Descriptors3D.NPR2, 470 "RadiusOfGyration": Descriptors3D.RadiusOfGyration, 471 "InertialShapeFactor": Descriptors3D.InertialShapeFactor, 472 "Eccentricity": Descriptors3D.Eccentricity, 473 "Asphericity": Descriptors3D.Asphericity, 474 "SpherocityIndex": Descriptors3D.SpherocityIndex, 475 } 476 477 for Name in NameToComputeFunctionMap: 478 ComputeFunction = NameToComputeFunctionMap[Name] 479 if Name in DescriptorNamesMap["ComputeFunction"]: 480 if PrintInfo: 481 MiscUtil.PrintWarning("Ignoring duplicate descriptor name: %s..." % Name) 482 else: 483 DescriptorNamesMap[Type]["Names"].append(Name) 484 DescriptorNamesMap["ComputeFunction"][Name] = ComputeFunction 485 486 # Check and add new 3D descriptors not directly available through Descriptors3D module... 487 Type = "3D" 488 AvailableName3DMap = {} 489 for Name in ["Autocorr3D", "RDF", "MORSE", "WHIM", "GETAWAY"]: 490 ComputeFunction = None 491 try: 492 if re.match("^Autocorr3D$", Name, re.I): 493 ComputeFunction = rdMolDescriptors.CalcAUTOCORR3D 494 elif re.match("^RDF$", Name, re.I): 495 ComputeFunction = rdMolDescriptors.CalcRDF 496 elif re.match("^MORSE$", Name, re.I): 497 ComputeFunction = rdMolDescriptors.CalcMORSE 498 elif re.match("^WHIM$", Name, re.I): 499 ComputeFunction = rdMolDescriptors.CalcWHIM 500 elif re.match("^GETAWAY$", Name, re.I): 501 ComputeFunction = rdMolDescriptors.CalcGETAWAY 502 else: 503 ComputeFunction = None 504 except AttributeError: 505 if PrintInfo: 506 MiscUtil.PrintWarning("3D descriptor, %s, is not available in your current version of RDKit" % Name) 507 508 if ComputeFunction is not None: 509 AvailableName3DMap[Name] = ComputeFunction 510 511 for Name in AvailableName3DMap: 512 ComputeFunction = AvailableName3DMap[Name] 513 if Name not in DescriptorNamesMap["ComputeFunction"]: 514 DescriptorNamesMap[Type]["Names"].append(Name) 515 DescriptorNamesMap["ComputeFunction"][Name] = ComputeFunction 516 517 Count = 0 518 TypesCount = [] 519 for Type in DescriptorNamesMap["Types"]: 520 TypeCount = len(DescriptorNamesMap[Type]["Names"]) 521 TypesCount.append(TypeCount) 522 Count += TypeCount 523 524 if not Count: 525 if PrintInfo: 526 MiscUtil.PrintError("Failed to retrieve any molecular descriptors...") 527 528 # Sort descriptor names... 529 for Type in DescriptorNamesMap["Types"]: 530 DescriptorNamesMap[Type]["Names"] = sorted(DescriptorNamesMap[Type]["Names"]) 531 532 if PrintInfo: 533 MiscUtil.PrintInfo("\nTotal number of availble molecular descriptors: %d" % Count) 534 for Index in range(0, len(DescriptorNamesMap["Types"])): 535 Type = DescriptorNamesMap["Types"][Index] 536 TypeCount = TypesCount[Index] 537 if PrintInfo: 538 MiscUtil.PrintInfo("Number of %s molecular descriptors: %d" % (Type, TypeCount)) 539 540 541 def ListMolecularDescriptorsInfo(): 542 """List descriptors information.""" 543 544 MiscUtil.PrintInfo("\nListing information for avalible molecular descriptors...") 545 546 Delimiter = ", " 547 for Type in DescriptorNamesMap["Types"]: 548 Names = DescriptorNamesMap[Type]["Names"] 549 MiscUtil.PrintInfo("\n%s descriptors: %s" % (Type, Delimiter.join(Names))) 550 551 MiscUtil.PrintInfo("") 552 553 554 def ProcessOptions(): 555 """Process and validate command line arguments and options.""" 556 557 MiscUtil.PrintInfo("Processing options...") 558 559 # Validate options... 560 ValidateOptions() 561 562 OptionsInfo["Autocorr2DExclude"] = True 563 if not re.match("^Yes$", Options["--autocorr2DExclude"], re.I): 564 OptionsInfo["Autocorr2DExclude"] = False 565 566 OptionsInfo["FragmentCount"] = True 567 if not re.match("^Yes$", Options["--fragmentCount"], re.I): 568 OptionsInfo["FragmentCount"] = False 569 570 OptionsInfo["DescriptorNames"] = Options["--descriptorNames"] 571 OptionsInfo["Mode"] = Options["--mode"] 572 573 OptionsInfo["Infile"] = Options["--infile"] 574 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters( 575 "--infileParams", Options["--infileParams"], Options["--infile"] 576 ) 577 578 OptionsInfo["Outfile"] = Options["--outfile"] 579 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters( 580 "--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"] 581 ) 582 583 OptionsInfo["Overwrite"] = Options["--overwrite"] 584 585 TextOutFileMode = False 586 TextOutFileDelim = "" 587 if MiscUtil.CheckFileExt(Options["--outfile"], "csv"): 588 TextOutFileMode = True 589 TextOutFileDelim = "," 590 elif MiscUtil.CheckFileExt(Options["--outfile"], "tsv txt"): 591 TextOutFileMode = True 592 TextOutFileDelim = "\t" 593 OptionsInfo["TextOutFileMode"] = TextOutFileMode 594 OptionsInfo["TextOutFileDelim"] = TextOutFileDelim 595 596 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False 597 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"]) 598 599 OptionsInfo["Precision"] = int(Options["--precision"]) 600 601 OptionsInfo["SMILESOut"] = False 602 if re.match("^Yes$", Options["--smilesOut"], re.I): 603 OptionsInfo["SMILESOut"] = True 604 605 606 def RetrieveOptions(): 607 """Retrieve command line arguments and options.""" 608 609 # Get options... 610 global Options 611 Options = docopt(_docoptUsage_) 612 613 # Set current working directory to the specified directory... 614 WorkingDir = Options["--workingdir"] 615 if WorkingDir: 616 os.chdir(WorkingDir) 617 618 # Handle examples option... 619 if "--examples" in Options and Options["--examples"]: 620 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 621 sys.exit(0) 622 623 # Handle listing of descriptor information... 624 if Options["--list"]: 625 ProcessListMolecularDescriptorsOption() 626 sys.exit(0) 627 628 629 def ProcessListMolecularDescriptorsOption(): 630 """Process list descriptors information.""" 631 632 RetrieveMolecularDescriptorsInfo() 633 ListMolecularDescriptorsInfo() 634 635 636 def ValidateOptions(): 637 """Validate option values.""" 638 639 MiscUtil.ValidateOptionTextValue("-a, --autocorr2DExclude", Options["--autocorr2DExclude"], "yes no") 640 MiscUtil.ValidateOptionTextValue("-f, --fragmentCount", Options["--fragmentCount"], "yes no") 641 642 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "2D 3D All FragmentCountOnly Specify") 643 644 if re.match("^Specify$", Options["--mode"], re.I): 645 if re.match("^none$", Options["--descriptorNames"], re.I): 646 MiscUtil.PrintError( 647 'The name(s) of molecular descriptors must be specified using "-d, --descriptorNames" option during "Specify" value of "-m, --mode" option.' 648 ) 649 650 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"]) 651 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi csv tsv txt") 652 653 if re.match("^3D|All$", Options["--mode"], re.I): 654 if MiscUtil.CheckFileExt(Options["--infile"], "smi"): 655 MiscUtil.PrintError( 656 'The input SMILES file, %s, is not valid for "3D or All" value of "-m, --mode" option.' 657 ) 658 659 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd csv tsv txt") 660 MiscUtil.ValidateOptionsOutputFileOverwrite( 661 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 662 ) 663 MiscUtil.ValidateOptionsDistinctFileNames( 664 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"] 665 ) 666 667 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no") 668 669 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0}) 670 MiscUtil.ValidateOptionTextValue("-s, --smilesOut", Options["--smilesOut"], "yes no") 671 672 673 # Setup a usage string for docopt... 674 _docoptUsage_ = """ 675 RDKitCalculateMolecularDescriptors.py - Calculate 2D/3D molecular descriptors 676 677 Usage: 678 RDKitCalculateMolecularDescriptors.py [--autocorr2DExclude <yes or no>] [--fragmentCount <yes or no>] 679 [--descriptorNames <Name1,Name2,...>] [--infileParams <Name,Value,...>] 680 [--mode <2D, 3D, All...>] [--mp <yes or no>] [--mpParams <Name,Value,...>] 681 [--outfileParams <Name,Value,...>] [--overwrite] [--precision <number>] 682 [--smilesOut <yes or no>] [-w <dir>] -i <infile> -o <outfile> 683 RDKitCalculateMolecularDescriptors.py -l | --list 684 RDKitCalculateMolecularDescriptors.py -h | --help | -e | --examples 685 686 Description: 687 Calculate 2D/3D molecular descriptors for molecules and write them out to a SD or 688 CSV/TSV text file. 689 690 The complete list of currently available molecular descriptors may be obtained by 691 using '-l, --list' option. The names of valid 2D, fragment count, and 3D molecular 692 descriptors are shown below: 693 694 2D descriptors: Autocorr2D, AvgIpc, BCUT2D_*, BalabanJ, BertzCT, Chi0, Chi1, Chi0n - Chi4n, 695 Chi0v - Chi4v, EState_VSA1 - EState_VSA11, ExactMolWt, FpDensityMorgan1, FpDensityMorgan2, 696 FpDensityMorgan3, FractionCSP3, HallKierAlpha, HeavyAtomCount, HeavyAtomMolWt, Ipc, 697 Kappa1 - Kappa3, LabuteASA, MaxAbsEStateIndex, MaxAbsPartialCharge, MaxEStateIndex, 698 MaxPartialCharge, MinAbsEStateIndex, MinAbsPartialCharge, MinEStateIndex, MinPartialCharge, 699 MolLogP, MolMR, MolWt, NHOHCount, NOCount, NumAliphaticCarbocycles, NumAliphaticHeterocycles, 700 NumAliphaticRings, NumAromaticCarbocycles, NumAromaticHeterocycles, NumAromaticRings, 701 NumHAcceptors, NumHDonors, NumHeteroatoms, NumRadicalElectrons, NumRotatableBonds, 702 NumSaturatedCarbocycles, NumSaturatedHeterocycles, NumSaturatedRings, NumValenceElectrons, 703 PEOE_VSA1 - PEOE_VSA14, RingCount, SMR_VSA1 - SMR_VSA10, SlogP_VSA1 - SlogP_VSA12, 704 TPSA, VSA_EState1 - VSA_EState10, qed 705 706 FragmentCount 2D descriptors: fr_Al_COO, fr_Al_OH, fr_Al_OH_noTert, fr_ArN, fr_Ar_COO, 707 fr_Ar_N, fr_Ar_NH, fr_Ar_OH, fr_COO, fr_COO2, fr_C_O, fr_C_O_noCOO, fr_C_S, fr_HOCCN, 708 fr_Imine, fr_NH0, fr_NH1, fr_NH2, fr_N_O, fr_Ndealkylation1, fr_Ndealkylation2, fr_Nhpyrrole, 709 fr_SH, fr_aldehyde, fr_alkyl_carbamate, fr_alkyl_halide, fr_allylic_oxid, fr_amide, fr_amidine, 710 fr_aniline, fr_aryl_methyl, fr_azide, fr_azo, fr_barbitur, fr_benzene, fr_benzodiazepine, 711 fr_bicyclic, fr_diazo, fr_dihydropyridine, fr_epoxide, fr_ester, fr_ether, fr_furan, fr_guanido, 712 fr_halogen, fr_hdrzine, fr_hdrzone, fr_imidazole, fr_imide, fr_isocyan, fr_isothiocyan, fr_ketone, 713 fr_ketone_Topliss, fr_lactam, fr_lactone, fr_methoxy, fr_morpholine, fr_nitrile, fr_nitro, 714 fr_nitro_arom, fr_nitro_arom_nonortho, fr_nitroso, fr_oxazole, fr_oxime, fr_para_hydroxylation, 715 fr_phenol, fr_phenol_noOrthoHbond, fr_phos_acid, fr_phos_ester, fr_piperdine, fr_piperzine, 716 fr_priamide, fr_prisulfonamd, fr_pyridine, fr_quatN, fr_sulfide, fr_sulfonamd, fr_sulfone, 717 fr_term_acetylene, fr_tetrazole, fr_thiazole, fr_thiocyan, fr_thiophene, fr_unbrch_alkane, fr_urea 718 719 3D descriptors: Asphericity, Autocorr3D, Eccentricity, GETAWAY, InertialShapeFactor, MORSE, 720 NPR1, NPR2, PMI1, PMI2, PMI3, RDF, RadiusOfGyration, SpherocityIndex, WHIM 721 722 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi, 723 .txt, .csv, .tsv) 724 725 The supported output file formats are: SD File (.sdf, .sd), CSV/TSV (.csv, .tsv, .txt) 726 727 Options: 728 -a, --autocorr2DExclude <yes or no> [default: yes] 729 Exclude Autocorr2D descriptor from the calculation of 2D descriptors. 730 -f, --fragmentCount <yes or no> [default: yes] 731 Include 2D fragment count descriptors during the calculation. These descriptors are 732 counted using SMARTS patterns specified in FragmentDescriptors.csv file distributed 733 with RDKit. This option is only used during '2D' or 'All' value of '-m, --mode' option. 734 -d, --descriptorNames <Name1,Name2,...> [default: none] 735 A comma delimited list of supported molecular descriptor names to calculate. 736 This option is only used during 'Specify' value of '-m, --mode' option. 737 -e, --examples 738 Print examples. 739 -h, --help 740 Print this help message. 741 -i, --infile <infile> 742 Input file name. 743 --infileParams <Name,Value,...> [default: auto] 744 A comma delimited list of parameter name and value pairs for reading 745 molecules from files. The supported parameter names for different file 746 formats, along with their default values, are shown below: 747 748 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes 749 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space, 750 smilesTitleLine,auto,sanitize,yes 751 752 Possible values for smilesDelimiter: space, comma or tab. 753 -l, --list 754 List molecular descriptors without performing any calculations. 755 -m, --mode <2D, 3D, All, FragmentCountOnly, or Specify> [default: 2D] 756 Type of molecular descriptors to calculate. Possible values: 2D, 3D, 757 All or Specify. The name of molecular descriptors must be specified using 758 '-d, --descriptorNames' for 'Specify'. 2D descriptors also include 1D descriptors. 759 The structure of molecules must contain 3D coordinates for the calculation 760 of 3D descriptors. 761 --mp <yes or no> [default: no] 762 Use multiprocessing. 763 764 By default, input data is retrieved in a lazy manner via mp.Pool.imap() 765 function employing lazy RDKit data iterable. This allows processing of 766 arbitrary large data sets without any additional requirements memory. 767 768 All input data may be optionally loaded into memory by mp.Pool.map() 769 before starting worker processes in a process pool by setting the value 770 of 'inputDataMode' to 'InMemory' in '--mpParams' option. 771 772 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input 773 data mode may adversely impact the performance. The '--mpParams' section 774 provides additional information to tune the value of 'chunkSize'. 775 --mpParams <Name,Value,...> [default: auto] 776 A comma delimited list of parameter name and value pairs to configure 777 multiprocessing. 778 779 The supported parameter names along with their default and possible 780 values are shown below: 781 782 chunkSize, auto 783 inputDataMode, Lazy [ Possible values: InMemory or Lazy ] 784 numProcesses, auto [ Default: mp.cpu_count() ] 785 786 These parameters are used by the following functions to configure and 787 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and 788 mp.Pool.imap(). 789 790 The chunkSize determines chunks of input data passed to each worker 791 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions. 792 The default value of chunkSize is dependent on the value of 'inputDataMode'. 793 794 The mp.Pool.map() function, invoked during 'InMemory' input data mode, 795 automatically converts RDKit data iterable into a list, loads all data into 796 memory, and calculates the default chunkSize using the following method 797 as shown in its code: 798 799 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4) 800 if extra: chunkSize += 1 801 802 For example, the default chunkSize will be 7 for a pool of 4 worker processes 803 and 100 data items. 804 805 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs 806 'lazy' RDKit data iterable to retrieve data as needed, without loading all the 807 data into memory. Consequently, the size of input data is not known a priori. 808 It's not possible to estimate an optimal value for the chunkSize. The default 809 chunkSize is set to 1. 810 811 The default value for the chunkSize during 'Lazy' data mode may adversely 812 impact the performance due to the overhead associated with exchanging 813 small chunks of data. It is generally a good idea to explicitly set chunkSize to 814 a larger value during 'Lazy' input data mode, based on the size of your input 815 data and number of processes in the process pool. 816 817 The mp.Pool.map() function waits for all worker processes to process all 818 the data and return the results. The mp.Pool.imap() function, however, 819 returns the the results obtained from worker processes as soon as the 820 results become available for specified chunks of data. 821 822 The order of data in the results returned by both mp.Pool.map() and 823 mp.Pool.imap() functions always corresponds to the input data. 824 -o, --outfile <outfile> 825 Output file name. 826 --outfileParams <Name,Value,...> [default: auto] 827 A comma delimited list of parameter name and value pairs for writing 828 molecules to files. The supported parameter names for different file 829 formats, along with their default values, are shown below: 830 831 SD: compute2DCoords,auto,kekulize,yes,forceV3000,no 832 833 Default value for compute2DCoords: yes for SMILES input file; no for all other 834 file types. 835 -p, --precision <number> [default: 3] 836 Floating point precision for writing the calculated descriptor values. 837 -s, --smilesOut <yes or no> [default: no] 838 Write out SMILES string to CSV/TSV text output file. 839 --overwrite 840 Overwrite existing files. 841 -w, --workingdir <dir> 842 Location of working directory which defaults to the current directory. 843 844 Examples: 845 To compute all available 2D descriptors except Autocorr2D descriptor and 846 write out a CSV file, type: 847 848 % RDKitCalculateMolecularDescriptors.py -i Sample.smi -o SampleOut.csv 849 850 To compute all available 2D descriptors except Autocorr2D descriptor in 851 multiprocessing mode on all available CPUs without loading all data into 852 memory, and write out a CSV file, type: 853 854 % RDKitCalculateMolecularDescriptors.py --mp yes -i Sample.smi 855 -o SampleOut.csv 856 857 To compute all available 2D descriptors except Autocorr2D descriptor in 858 multiprocessing mode on all available CPUs by loading all data into memory, 859 and write out a CSV file, type: 860 861 % RDKitCalculateMolecularDescriptors.py --mp yes --mpParams 862 "inputDataMode,InMemory" -i Sample.smi -o SampleOut.csv 863 864 To compute all available 2D descriptors except Autocorr2D descriptor in 865 multiprocessing mode on specific number of CPUs and chunk size without 866 loading all data into memory, and write out a SDF file, type: 867 868 % RDKitCalculateMolecularDescriptors.py --mp yes --mpParams 869 "inputDataMode,Lazy,numProcesses,4,chunkSize,8" -i Sample.smi 870 -o SampleOut.sdf 871 872 To compute all available 2D descriptors including Autocorr2D descriptor and 873 excluding fragment count descriptors, and write out a TSV file, type: 874 875 % RDKitCalculateMolecularDescriptors.py -m 2D -a no -f no 876 -i Sample.smi -o SampleOut.tsv 877 878 To compute all available 3D descriptors and write out a SD file, type: 879 880 % RDKitCalculateMolecularDescriptors.py -m 3D -i Sample3D.sdf 881 -o Sample3DOut.sdf 882 883 To compute only fragment count 2D descriptors and write out a SD 884 file file, type: 885 886 % RDKitCalculateMolecularDescriptors.py -m FragmentCountOnly 887 -i Sample.sdf -o SampleOut.sdf 888 889 To compute all available 2D and 3D descriptors including fragment count and 890 Autocorr2D and write out a CSV file, type: 891 892 % RDKitCalculateMolecularDescriptors.py -m All -a no -i Sample.sdf 893 -o SampleOut.csv 894 895 To compute a specific set of 2D and 3D descriptors and write out a 896 write out a TSV file, type: 897 898 % RDKitCalculateMolecularDescriptors.py -m specify 899 -d 'MolWt,MolLogP,NHOHCount, NOCount,RadiusOfGyration' 900 -i Sample3D.sdf -o SampleOut.csv 901 902 To compute all available 2D descriptors except Autocorr2D descriptor for 903 molecules in a CSV SMILES file, SMILES strings in column 1, name in 904 column 2, and write out a SD file without calculation of 2D coordinates, type: 905 906 % RDKitCalculateMolecularDescriptors.py --infileParams 907 "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1, 908 smilesNameColumn,2" --outfileParams "compute2DCoords,no" 909 -i SampleSMILES.csv -o SampleOut.sdf 910 911 Author: 912 Manish Sud(msud@san.rr.com) 913 914 See also: 915 RDKitCalculateRMSD.py, RDKitCompareMoleculeShapes.py, RDKitConvertFileFormat.py, 916 RDKitGenerateConformers.py, RDKitPerformMinimization.py 917 918 Copyright: 919 Copyright (C) 2026 Manish Sud. All rights reserved. 920 921 The functionality available in this script is implemented using RDKit, an 922 open source toolkit for cheminformatics developed by Greg Landrum. 923 924 This file is part of MayaChemTools. 925 926 MayaChemTools is free software; you can redistribute it and/or modify it under 927 the terms of the GNU Lesser General Public License as published by the Free 928 Software Foundation; either version 3 of the License, or (at your option) any 929 later version. 930 931 """ 932 933 if __name__ == "__main__": 934 main()
