1 #!/bin/env python
2 #
3 # File: RDKitFilterTorsionStrainEnergyAlerts.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Collaborator: Pat Walters
7 #
8 # Copyright (C) 2026 Manish Sud. All rights reserved.
9 #
10 # This script uses the torsion strain energy library developed by Gu, S.;
11 # Smith, M. S.; Yang, Y.; Irwin, J. J.; Shoichet, B. K. [ Ref 153 ].
12 #
13 # The torsion strain enegy library is based on the Torsion Library jointly
14 # developed by the University of Hamburg, Center for Bioinformatics,
15 # Hamburg, Germany and F. Hoffmann-La-Roche Ltd., Basel, Switzerland.
16 #
17 # The functionality available in this script is implemented using RDKit, an
18 # open source toolkit for cheminformatics developed by Greg Landrum.
19 #
20 # This file is part of MayaChemTools.
21 #
22 # MayaChemTools is free software; you can redistribute it and/or modify it under
23 # the terms of the GNU Lesser General Public License as published by the Free
24 # Software Foundation; either version 3 of the License, or (at your option) any
25 # later version.
26 #
27 # MayaChemTools is distributed in the hope that it will be useful, but without
28 # any warranty; without even the implied warranty of merchantability of fitness
29 # for a particular purpose. See the GNU Lesser General Public License for more
30 # details.
31 #
32 # You should have received a copy of the GNU Lesser General Public License
33 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
34 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
35 # Boston, MA, 02111-1307, USA.
36 #
37
38 from __future__ import print_function
39
40 import os
41 import sys
42 import time
43 import re
44 import glob
45 import multiprocessing as mp
46 import math
47
48 # RDKit imports...
49 try:
50 from rdkit import rdBase
51 from rdkit import Chem
52 except ImportError as ErrMsg:
53 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
54 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
55 sys.exit(1)
56
57 # MayaChemTools imports...
58 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
59 try:
60 from docopt import docopt
61 import MiscUtil
62 import RDKitUtil
63 from TorsionAlerts.TorsionStrainEnergyAlerts import TorsionStrainEnergyAlerts
64 except ImportError as ErrMsg:
65 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
66 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
67 sys.exit(1)
68
69 ScriptName = os.path.basename(sys.argv[0])
70 Options = {}
71 OptionsInfo = {}
72
73
74 def main():
75 """Start execution of the script."""
76
77 MiscUtil.PrintInfo(
78 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
79 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
80 )
81
82 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
83
84 # Retrieve command line arguments and options...
85 RetrieveOptions()
86
87 if Options["--list"]:
88 # Handle listing of torsion library information...
89 ProcessListTorsionLibraryOption()
90 else:
91 # Process and validate command line arguments and options...
92 ProcessOptions()
93
94 # Perform actions required by the script...
95 PerformFiltering()
96
97 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
98 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
99
100
101 def PerformFiltering():
102 """Filter molecules using SMARTS torsion rules in the torsion strain energy
103 library file."""
104
105 # Setup a molecule reader...
106 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
107 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
108
109 MolCount, ValidMolCount, RemainingMolCount, WriteFailedCount = ProcessMolecules(Mols)
110
111 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
112 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
113 MiscUtil.PrintInfo("Number of molecules failed during writing: %d" % WriteFailedCount)
114 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + WriteFailedCount))
115
116 MiscUtil.PrintInfo("\nNumber of remaining molecules: %d" % RemainingMolCount)
117 MiscUtil.PrintInfo("Number of filtered molecules: %d" % (ValidMolCount - RemainingMolCount))
118
119
120 def ProcessMolecules(Mols):
121 """Process and filter molecules."""
122
123 if OptionsInfo["MPMode"]:
124 return ProcessMoleculesUsingMultipleProcesses(Mols)
125 else:
126 return ProcessMoleculesUsingSingleProcess(Mols)
127
128
129 def ProcessMoleculesUsingSingleProcess(Mols):
130 """Process and filter molecules using a single process."""
131
132 # Instantiate torsion strain energy alerts class...
133 TorsionStrainEnergyAlertsHandle = InstantiateTorsionStrainEnergyAlertsClass()
134
135 MiscUtil.PrintInfo("\nFiltering molecules...")
136
137 OutfileFilteredMode = OptionsInfo["OutfileFilteredMode"]
138
139 # Set up writers...
140 OutfilesWriters = SetupOutfilesWriters()
141
142 WriterRemaining = OutfilesWriters["WriterRemaining"]
143 WriterFiltered = OutfilesWriters["WriterFiltered"]
144 WriterAlertSummary = OutfilesWriters["WriterAlertSummary"]
145
146 # Initialize alerts summary info...
147 TorsionAlertsSummaryInfo = InitializeTorsionAlertsSummaryInfo()
148
149 (MolCount, ValidMolCount, RemainingMolCount, WriteFailedCount, FilteredMolWriteCount) = [0] * 5
150 for Mol in Mols:
151 MolCount += 1
152
153 if Mol is None:
154 continue
155
156 if RDKitUtil.IsMolEmpty(Mol):
157 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % RDKitUtil.GetMolName(Mol, MolCount))
158 continue
159
160 # Check for 3D flag...
161 if not Mol.GetConformer().Is3D():
162 MiscUtil.PrintWarning("3D tag is not set. Ignoring molecule: %s\n" % RDKitUtil.GetMolName(Mol, MolCount))
163 continue
164
165 ValidMolCount += 1
166
167 # Identify torsion library alerts for rotatable bonds..
168 RotBondsAlertsStatus, RotBondsAlertsInfo = (
169 TorsionStrainEnergyAlertsHandle.IdentifyTorsionLibraryAlertsForRotatableBonds(Mol)
170 )
171
172 TrackTorsionAlertsSummaryInfo(TorsionAlertsSummaryInfo, RotBondsAlertsInfo)
173
174 # Write out filtered and remaining molecules...
175 WriteStatus = True
176 if RotBondsAlertsStatus:
177 if OutfileFilteredMode:
178 WriteStatus = WriteMolecule(WriterFiltered, Mol, RotBondsAlertsInfo)
179 if WriteStatus:
180 FilteredMolWriteCount += 1
181 else:
182 RemainingMolCount += 1
183 WriteStatus = WriteMolecule(WriterRemaining, Mol, RotBondsAlertsInfo)
184
185 if not WriteStatus:
186 WriteFailedCount += 1
187
188 WriteTorsionAlertsSummaryInfo(WriterAlertSummary, TorsionAlertsSummaryInfo)
189 CloseOutfilesWriters(OutfilesWriters)
190
191 if FilteredMolWriteCount:
192 WriteTorsionAlertsFilteredByRulesInfo(TorsionAlertsSummaryInfo)
193
194 return (MolCount, ValidMolCount, RemainingMolCount, WriteFailedCount)
195
196
197 def ProcessMoleculesUsingMultipleProcesses(Mols):
198 """Process and filter molecules using multiprocessing."""
199
200 MiscUtil.PrintInfo("\nFiltering molecules using multiprocessing...")
201
202 MPParams = OptionsInfo["MPParams"]
203 OutfileFilteredMode = OptionsInfo["OutfileFilteredMode"]
204
205 # Instantiate torsion strain energy alerts class to list torsion library information...
206 InstantiateTorsionStrainEnergyAlertsClass()
207
208 # Set up writers...
209 OutfilesWriters = SetupOutfilesWriters()
210
211 WriterRemaining = OutfilesWriters["WriterRemaining"]
212 WriterFiltered = OutfilesWriters["WriterFiltered"]
213 WriterAlertSummary = OutfilesWriters["WriterAlertSummary"]
214
215 # Initialize alerts summary info...
216 TorsionAlertsSummaryInfo = InitializeTorsionAlertsSummaryInfo()
217
218 # Setup data for initializing a worker process...
219 MiscUtil.PrintInfo("\nEncoding options info and rotatable bond pattern molecule...")
220 InitializeWorkerProcessArgs = (
221 MiscUtil.ObjectToBase64EncodedString(Options),
222 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
223 )
224
225 # Setup a encoded mols data iterable for a worker process...
226 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
227
228 # Setup process pool along with data initialization for each process...
229 MiscUtil.PrintInfo(
230 "\nConfiguring multiprocessing using %s method..."
231 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
232 )
233 MiscUtil.PrintInfo(
234 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
235 % (
236 MPParams["NumProcesses"],
237 MPParams["InputDataMode"],
238 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
239 )
240 )
241
242 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
243
244 # Start processing...
245 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
246 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
247 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
248 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
249 else:
250 MiscUtil.PrintError(
251 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
252 )
253
254 (MolCount, ValidMolCount, RemainingMolCount, WriteFailedCount, FilteredMolWriteCount) = [0] * 5
255 for Result in Results:
256 MolCount += 1
257 MolIndex, EncodedMol, RotBondsAlertsStatus, RotBondsAlertsInfo = Result
258
259 if EncodedMol is None:
260 continue
261 ValidMolCount += 1
262
263 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
264
265 TrackTorsionAlertsSummaryInfo(TorsionAlertsSummaryInfo, RotBondsAlertsInfo)
266
267 # Write out filtered and remaining molecules...
268 WriteStatus = True
269 if RotBondsAlertsStatus:
270 if OutfileFilteredMode:
271 WriteStatus = WriteMolecule(WriterFiltered, Mol, RotBondsAlertsInfo)
272 if WriteStatus:
273 FilteredMolWriteCount += 1
274 else:
275 RemainingMolCount += 1
276 WriteStatus = WriteMolecule(WriterRemaining, Mol, RotBondsAlertsInfo)
277
278 if not WriteStatus:
279 WriteFailedCount += 1
280
281 WriteTorsionAlertsSummaryInfo(WriterAlertSummary, TorsionAlertsSummaryInfo)
282 CloseOutfilesWriters(OutfilesWriters)
283
284 if FilteredMolWriteCount:
285 WriteTorsionAlertsFilteredByRulesInfo(TorsionAlertsSummaryInfo)
286
287 return (MolCount, ValidMolCount, RemainingMolCount, WriteFailedCount)
288
289
290 def InitializeWorkerProcess(*EncodedArgs):
291 """Initialize data for a worker process."""
292
293 global Options, OptionsInfo
294
295 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
296
297 # Decode Options and OptionInfo...
298 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
299 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
300
301 # Instantiate torsion strain energy alerts class...
302 OptionsInfo["TorsionStrainEnergyAlertsHandle"] = InstantiateTorsionStrainEnergyAlertsClass(Quiet=True)
303
304
305 def WorkerProcess(EncodedMolInfo):
306 """Process data for a worker process."""
307
308 MolIndex, EncodedMol = EncodedMolInfo
309
310 if EncodedMol is None:
311 return [MolIndex, None, False, None]
312
313 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
314 if RDKitUtil.IsMolEmpty(Mol):
315 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
316 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
317 return [MolIndex, None, False, None]
318
319 # Check for 3D flag...
320 if not Mol.GetConformer().Is3D():
321 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
322 MiscUtil.PrintWarning("3D tag is not set. Ignoring molecule: %s\n" % MolName)
323 return [MolIndex, None, False, None]
324
325 # Identify torsion library alerts for rotatable bonds..
326 TorsionStrainEnergyAlertsHandle = OptionsInfo["TorsionStrainEnergyAlertsHandle"]
327 RotBondsAlertsStatus, RotBondsAlertsInfo = (
328 TorsionStrainEnergyAlertsHandle.IdentifyTorsionLibraryAlertsForRotatableBonds(Mol)
329 )
330
331 return [MolIndex, EncodedMol, RotBondsAlertsStatus, RotBondsAlertsInfo]
332
333
334 def InitializeTorsionAlertsSummaryInfo():
335 """Initialize torsion alerts summary."""
336
337 if OptionsInfo["CountMode"]:
338 return None
339
340 if not OptionsInfo["TrackAlertsSummaryInfo"]:
341 return None
342
343 TorsionAlertsSummaryInfo = {}
344 TorsionAlertsSummaryInfo["RuleIDs"] = []
345
346 for DataLabel in [
347 "SMARTSToRuleIDs",
348 "RuleSMARTS",
349 "HierarchyClassName",
350 "HierarchySubClassName",
351 "EnergyMethod",
352 "MaxSingleEnergyAlertTypes",
353 "MaxSingleEnergyAlertTypesMolCount",
354 ]:
355 TorsionAlertsSummaryInfo[DataLabel] = {}
356
357 return TorsionAlertsSummaryInfo
358
359
360 def TrackTorsionAlertsSummaryInfo(TorsionAlertsSummaryInfo, RotBondsAlertsInfo):
361 """Track torsion alerts summary information for matched torsion rules in a
362 molecule."""
363
364 if OptionsInfo["CountMode"]:
365 return
366
367 if not OptionsInfo["TrackAlertsSummaryInfo"]:
368 return
369
370 if RotBondsAlertsInfo is None:
371 return
372
373 MolAlertsInfo = {}
374 MolAlertsInfo["RuleIDs"] = []
375 MolAlertsInfo["MaxSingleEnergyAlertTypes"] = {}
376
377 for ID in RotBondsAlertsInfo["IDs"]:
378 if not RotBondsAlertsInfo["MatchStatus"][ID]:
379 continue
380
381 if SkipRotatableBondAlertInfo(ID, RotBondsAlertsInfo):
382 continue
383
384 MaxSingleEnergyAlertType = SetupMaxSingleEnergyAlertStatusValue(
385 RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID]
386 )
387
388 TorsionRuleNodeID = RotBondsAlertsInfo["TorsionRuleNodeID"][ID]
389 TorsionRuleSMARTS = RotBondsAlertsInfo["TorsionRuleSMARTS"][ID]
390
391 # Track data for torsion alert summary information across molecules...
392 if TorsionRuleNodeID not in TorsionAlertsSummaryInfo["RuleSMARTS"]:
393 TorsionAlertsSummaryInfo["RuleIDs"].append(TorsionRuleNodeID)
394 TorsionAlertsSummaryInfo["SMARTSToRuleIDs"][TorsionRuleSMARTS] = TorsionRuleNodeID
395
396 TorsionAlertsSummaryInfo["RuleSMARTS"][TorsionRuleNodeID] = TorsionRuleSMARTS
397 TorsionAlertsSummaryInfo["HierarchyClassName"][TorsionRuleNodeID] = RotBondsAlertsInfo[
398 "HierarchyClassName"
399 ][ID]
400 TorsionAlertsSummaryInfo["HierarchySubClassName"][TorsionRuleNodeID] = RotBondsAlertsInfo[
401 "HierarchySubClassName"
402 ][ID]
403
404 TorsionAlertsSummaryInfo["EnergyMethod"][TorsionRuleNodeID] = RotBondsAlertsInfo["EnergyMethod"][ID]
405
406 # Initialize number of alert types across all molecules...
407 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID] = {}
408
409 # Initialize number of molecules flagged by each alert type...
410 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypesMolCount"][TorsionRuleNodeID] = {}
411
412 if MaxSingleEnergyAlertType not in TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID]:
413 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID][MaxSingleEnergyAlertType] = 0
414 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypesMolCount"][TorsionRuleNodeID][
415 MaxSingleEnergyAlertType
416 ] = 0
417
418 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID][MaxSingleEnergyAlertType] += 1
419
420 # Track data for torsion alert information in a molecule...
421 if TorsionRuleNodeID not in MolAlertsInfo["MaxSingleEnergyAlertTypes"]:
422 MolAlertsInfo["RuleIDs"].append(TorsionRuleNodeID)
423 MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID] = {}
424
425 if MaxSingleEnergyAlertType not in MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID]:
426 MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID][MaxSingleEnergyAlertType] = 0
427 MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID][MaxSingleEnergyAlertType] += 1
428
429 # Track number of molecules flagged by a specific torsion alert...
430 for TorsionRuleNodeID in MolAlertsInfo["RuleIDs"]:
431 for MaxSingleEnergyAlertType in MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID]:
432 if MolAlertsInfo["MaxSingleEnergyAlertTypes"][TorsionRuleNodeID][MaxSingleEnergyAlertType]:
433 TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypesMolCount"][TorsionRuleNodeID][
434 MaxSingleEnergyAlertType
435 ] += 1
436
437
438 def WriteTorsionAlertsSummaryInfo(Writer, TorsionAlertsSummaryInfo):
439 """Write out torsion alerts summary informatio to a CSV file."""
440
441 if OptionsInfo["CountMode"]:
442 return
443
444 if not OptionsInfo["OutfileSummaryMode"]:
445 return
446
447 if len(TorsionAlertsSummaryInfo["RuleIDs"]) == 0:
448 return
449
450 # Write headers...
451 QuoteValues = True
452 Values = [
453 "TorsionRule",
454 "HierarchyClass",
455 "HierarchySubClass",
456 "EnergyMethod",
457 "MaxSingleEnergyTorsionAlertTypes",
458 "MaxSingleEnergyTorsionAlertCount",
459 "MaxSingleEnergyTorsionAlertMolCount",
460 ]
461 Writer.write("%s\n" % MiscUtil.JoinWords(Values, ",", QuoteValues))
462
463 SortedRuleIDs = GetSortedTorsionAlertsSummaryInfoRuleIDs(TorsionAlertsSummaryInfo)
464
465 # Write alerts information...
466 for ID in SortedRuleIDs:
467 # Remove any double quotes in SMARTS...
468 RuleSMARTS = TorsionAlertsSummaryInfo["RuleSMARTS"][ID]
469 RuleSMARTS = re.sub('"', "", RuleSMARTS, flags=re.I)
470
471 HierarchyClassName = TorsionAlertsSummaryInfo["HierarchyClassName"][ID]
472 HierarchySubClassName = TorsionAlertsSummaryInfo["HierarchySubClassName"][ID]
473
474 EnergyMethod = TorsionAlertsSummaryInfo["EnergyMethod"][ID]
475
476 MaxSingleEnergyAlertTypes = []
477 MaxSingleEnergyAlertTypesCount = []
478 MaxSingleEnergyAlertTypesMolCount = []
479 for MaxSingleEnergyAlertType in sorted(TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][ID]):
480 MaxSingleEnergyAlertTypes.append(MaxSingleEnergyAlertType)
481 MaxSingleEnergyAlertTypesCount.append(
482 "%s" % TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][ID][MaxSingleEnergyAlertType]
483 )
484 MaxSingleEnergyAlertTypesMolCount.append(
485 "%s" % TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypesMolCount"][ID][MaxSingleEnergyAlertType]
486 )
487
488 Values = [
489 RuleSMARTS,
490 HierarchyClassName,
491 HierarchySubClassName,
492 EnergyMethod,
493 "%s" % MiscUtil.JoinWords(MaxSingleEnergyAlertTypes, ","),
494 "%s" % (MiscUtil.JoinWords(MaxSingleEnergyAlertTypesCount, ",")),
495 "%s" % (MiscUtil.JoinWords(MaxSingleEnergyAlertTypesMolCount, ",")),
496 ]
497 Writer.write("%s\n" % MiscUtil.JoinWords(Values, ",", QuoteValues))
498
499
500 def GetSortedTorsionAlertsSummaryInfoRuleIDs(TorsionAlertsSummaryInfo):
501 """Sort torsion rule IDs by alert types molecule count in descending order."""
502
503 SortedRuleIDs = []
504
505 RuleIDs = TorsionAlertsSummaryInfo["RuleIDs"]
506 if len(RuleIDs) == 0:
507 return SortedRuleIDs
508
509 # Setup a map from AlertTypesMolCount to IDs for sorting alerts...
510 RuleIDs = TorsionAlertsSummaryInfo["RuleIDs"]
511 MolCountMap = {}
512 for ID in RuleIDs:
513 MolCount = 0
514 for AlertType in sorted(TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypes"][ID]):
515 MolCount += TorsionAlertsSummaryInfo["MaxSingleEnergyAlertTypesMolCount"][ID][AlertType]
516 MolCountMap[ID] = MolCount
517
518 SortedRuleIDs = sorted(RuleIDs, key=lambda ID: MolCountMap[ID], reverse=True)
519
520 return SortedRuleIDs
521
522
523 def WriteTorsionAlertsFilteredByRulesInfo(TorsionAlertsSummaryInfo):
524 """Write out torsion alerts SD files for individual torsion rules."""
525
526 if OptionsInfo["CountMode"]:
527 return
528
529 if not OptionsInfo["OutfilesFilteredByRulesMode"]:
530 return
531
532 if len(TorsionAlertsSummaryInfo["RuleIDs"]) == 0:
533 return
534
535 # Setup a molecule reader for filtered molecules...
536 FilteredMols = RDKitUtil.ReadMolecules(OptionsInfo["OutfileFiltered"], **OptionsInfo["InfileParams"])
537
538 # Get torsion rule IDs for writing out filtered SD files for individual torsion alert rules...
539 TorsionRuleIDs = GetTorsionAlertsFilteredByRuleFilesRuleIDs(TorsionAlertsSummaryInfo)
540
541 # Setup writers...
542 ByRuleOutfilesWriters = SetupByRuleOutfilesWriters(TorsionRuleIDs)
543
544 for Mol in FilteredMols:
545 # Retrieve torsion alerts info...
546 TorsionAlertsInfo = RetrieveTorsionAlertsInfo(Mol, TorsionAlertsSummaryInfo)
547 if TorsionAlertsInfo is None:
548 continue
549
550 for TorsionRuleID in TorsionRuleIDs:
551 if TorsionRuleID not in TorsionAlertsInfo["RuleSMARTS"]:
552 continue
553
554 WriteMoleculeFilteredByRuleID(
555 ByRuleOutfilesWriters[TorsionRuleID], Mol, TorsionRuleID, TorsionAlertsSummaryInfo, TorsionAlertsInfo
556 )
557
558 CloseByRuleOutfilesWriters(ByRuleOutfilesWriters)
559
560
561 def GetTorsionAlertsFilteredByRuleFilesRuleIDs(TorsionAlertsSummaryInfo):
562 """Get torsion rule IDs for writing out individual SD files filtered by torsion alert rules."""
563
564 # Get torsion rule IDs triggering torsion alerts sorted in the order from the most to
565 # the least number of unique molecules...
566 RuleIDs = GetSortedTorsionAlertsSummaryInfoRuleIDs(TorsionAlertsSummaryInfo)
567
568 # Select torsion rule IDs for writing out SD files...
569 if not OptionsInfo["OutfilesFilteredByRulesAllMode"]:
570 MaxRuleIDs = OptionsInfo["OutfilesFilteredByRulesMaxCount"]
571 if MaxRuleIDs < len(RuleIDs):
572 RuleIDs = RuleIDs[0:MaxRuleIDs]
573
574 return RuleIDs
575
576
577 def RetrieveTorsionAlertsInfo(Mol, TorsionAlertsSummaryInfo):
578 """Parse torsion alerts data field value to retrieve alerts information for rotatable bonds."""
579
580 TorsionAlertsLabel = OptionsInfo["SDFieldIDsToLabels"]["TorsionAlertsLabel"]
581 TorsionAlerts = Mol.GetProp(TorsionAlertsLabel) if Mol.HasProp(TorsionAlertsLabel) else None
582
583 if TorsionAlerts is None or len(TorsionAlerts) == 0:
584 return None
585
586 # Initialize for tracking by rule IDs...
587 TorsionAlertsInfo = {}
588 TorsionAlertsInfo["RuleIDs"] = []
589
590 TorsionAlertsInfo["RuleSMARTS"] = {}
591 TorsionAlertsInfo["HierarchyClassName"] = {}
592 TorsionAlertsInfo["HierarchySubClassName"] = {}
593 TorsionAlertsInfo["EnergyMethod"] = {}
594
595 TorsionAlertsInfo["AtomIndices"] = {}
596 TorsionAlertsInfo["TorsionAtomIndices"] = {}
597 TorsionAlertsInfo["TorsionAngle"] = {}
598
599 TorsionAlertsInfo["Energy"] = {}
600 TorsionAlertsInfo["EnergyLowerBound"] = {}
601 TorsionAlertsInfo["EnergyUpperBound"] = {}
602
603 TorsionAlertsInfo["AngleNotObserved"] = {}
604 TorsionAlertsInfo["MaxSingleEnergyAlertType"] = {}
605
606 TorsionAlertsInfo["AnglesNotObservedCount"] = {}
607 TorsionAlertsInfo["MaxSingleEnergyAlertsCount"] = {}
608
609 ValuesDelimiter = OptionsInfo["IntraSetValuesDelim"]
610 TorsionAlertsSetSize = 12
611
612 TorsionAlertsWords = TorsionAlerts.split()
613 if len(TorsionAlertsWords) % TorsionAlertsSetSize:
614 MiscUtil.PrintError(
615 "The number of space delimited values, %s, for TorsionAlerts data field in filtered SD file must be a multiple of %s."
616 % (len(TorsionAlertsWords), TorsionAlertsSetSize)
617 )
618
619 ID = 0
620 for Index in range(0, len(TorsionAlertsWords), TorsionAlertsSetSize):
621 ID += 1
622
623 (
624 RotBondIndices,
625 TorsionIndices,
626 TorsionAngle,
627 Energy,
628 EnergyLowerBound,
629 EnergyUpperBound,
630 HierarchyClass,
631 HierarchySubClass,
632 TorsionRule,
633 EnergyMethod,
634 AngleNotObserved,
635 MaxSingleEnergyAlertType,
636 ) = TorsionAlertsWords[Index : Index + TorsionAlertsSetSize]
637 RotBondIndices = RotBondIndices.split(ValuesDelimiter)
638 TorsionIndices = TorsionIndices.split(ValuesDelimiter)
639
640 if TorsionRule not in TorsionAlertsSummaryInfo["SMARTSToRuleIDs"]:
641 MiscUtil.PrintWarning(
642 "The SMARTS pattern, %s, for TorsionAlerts data field in filtered SD file doesn't map to any torsion rule..."
643 % TorsionRule
644 )
645 continue
646 TorsionRuleNodeID = TorsionAlertsSummaryInfo["SMARTSToRuleIDs"][TorsionRule]
647
648 # Track data for torsion alerts in a molecule...
649 if TorsionRuleNodeID not in TorsionAlertsInfo["RuleSMARTS"]:
650 TorsionAlertsInfo["RuleIDs"].append(TorsionRuleNodeID)
651
652 TorsionAlertsInfo["RuleSMARTS"][TorsionRuleNodeID] = TorsionRule
653 TorsionAlertsInfo["HierarchyClassName"][TorsionRuleNodeID] = HierarchyClass
654 TorsionAlertsInfo["HierarchySubClassName"][TorsionRuleNodeID] = HierarchySubClass
655 TorsionAlertsInfo["EnergyMethod"][TorsionRuleNodeID] = EnergyMethod
656
657 TorsionAlertsInfo["AtomIndices"][TorsionRuleNodeID] = []
658 TorsionAlertsInfo["TorsionAtomIndices"][TorsionRuleNodeID] = []
659 TorsionAlertsInfo["TorsionAngle"][TorsionRuleNodeID] = []
660
661 TorsionAlertsInfo["Energy"][TorsionRuleNodeID] = []
662 TorsionAlertsInfo["EnergyLowerBound"][TorsionRuleNodeID] = []
663 TorsionAlertsInfo["EnergyUpperBound"][TorsionRuleNodeID] = []
664 TorsionAlertsInfo["AngleNotObserved"][TorsionRuleNodeID] = []
665 TorsionAlertsInfo["MaxSingleEnergyAlertType"][TorsionRuleNodeID] = []
666
667 TorsionAlertsInfo["AnglesNotObservedCount"][TorsionRuleNodeID] = 0
668 TorsionAlertsInfo["MaxSingleEnergyAlertsCount"][TorsionRuleNodeID] = 0
669
670 # Track multiple values for a rule ID...
671 TorsionAlertsInfo["AtomIndices"][TorsionRuleNodeID].append(RotBondIndices)
672 TorsionAlertsInfo["TorsionAtomIndices"][TorsionRuleNodeID].append(TorsionIndices)
673 TorsionAlertsInfo["TorsionAngle"][TorsionRuleNodeID].append(TorsionAngle)
674
675 TorsionAlertsInfo["Energy"][TorsionRuleNodeID].append(Energy)
676 TorsionAlertsInfo["EnergyLowerBound"][TorsionRuleNodeID].append(EnergyLowerBound)
677 TorsionAlertsInfo["EnergyUpperBound"][TorsionRuleNodeID].append(EnergyUpperBound)
678 TorsionAlertsInfo["AngleNotObserved"][TorsionRuleNodeID].append(AngleNotObserved)
679
680 TorsionAlertsInfo["MaxSingleEnergyAlertType"][TorsionRuleNodeID].append(MaxSingleEnergyAlertType)
681
682 # Count angles not observer for a rule ID...
683 if AngleNotObserved == "Yes":
684 TorsionAlertsInfo["AnglesNotObservedCount"][TorsionRuleNodeID] += 1
685
686 # Count max single energy alert for a rule ID...
687 if MaxSingleEnergyAlertType == "Yes":
688 TorsionAlertsInfo["MaxSingleEnergyAlertsCount"][TorsionRuleNodeID] += 1
689
690 return TorsionAlertsInfo
691
692
693 def WriteMolecule(Writer, Mol, RotBondsAlertsInfo):
694 """Write out molecule."""
695
696 if OptionsInfo["CountMode"]:
697 return True
698
699 SetupMolPropertiesForAlertsInformation(Mol, RotBondsAlertsInfo)
700
701 try:
702 Writer.write(Mol)
703 except Exception as ErrMsg:
704 MiscUtil.PrintWarning("Failed to write molecule %s:\n%s\n" % (RDKitUtil.GetMolName(Mol), ErrMsg))
705 return False
706
707 return True
708
709
710 def SetupMolPropertiesForAlertsInformation(Mol, RotBondsAlertsInfo):
711 """Setup molecule properties containing alerts information for rotatable bonds."""
712
713 if not OptionsInfo["OutfileAlerts"]:
714 return
715
716 SDFieldIDsToLabels = OptionsInfo["SDFieldIDsToLabels"]
717 Precision = OptionsInfo["Precision"]
718
719 # Setup rotatable bonds count...
720 RotBondsCount = 0
721 if RotBondsAlertsInfo is not None:
722 RotBondsCount = len(RotBondsAlertsInfo["IDs"])
723 Mol.SetProp(SDFieldIDsToLabels["RotBondsCountLabel"], "%s" % RotBondsCount)
724
725 if RotBondsAlertsInfo is not None:
726 # Setup total energy along with lower and upper bounds...
727 Mol.SetProp(
728 SDFieldIDsToLabels["TotalEnergyLabel"],
729 "%s" % SetupEnergyValueForSDField(RotBondsAlertsInfo["TotalEnergy"], Precision),
730 )
731 Mol.SetProp(
732 SDFieldIDsToLabels["TotalEnergyLowerBoundCILabel"],
733 "%s" % SetupEnergyValueForSDField(RotBondsAlertsInfo["TotalEnergyLowerBound"], Precision),
734 )
735 Mol.SetProp(
736 SDFieldIDsToLabels["TotalEnergyUpperBoundCILabel"],
737 "%s" % SetupEnergyValueForSDField(RotBondsAlertsInfo["TotalEnergyUpperBound"], Precision),
738 )
739
740 # Setup max single energy and alert count...
741 if OptionsInfo["MaxSingleEnergyMode"] or OptionsInfo["TotalOrMaxSingleEnergyMode"]:
742 Mol.SetProp(
743 SDFieldIDsToLabels["MaxSingleEnergyLabel"],
744 "%s" % SetupEnergyValueForSDField(RotBondsAlertsInfo["MaxSingleEnergy"], Precision),
745 )
746 Mol.SetProp(
747 SDFieldIDsToLabels["MaxSingleEnergyAlertsCountLabel"],
748 "%s"
749 % (
750 "NA"
751 if RotBondsAlertsInfo["MaxSingleEnergyAlertsCount"] is None
752 else RotBondsAlertsInfo["MaxSingleEnergyAlertsCount"]
753 ),
754 )
755
756 Mol.SetProp(
757 SDFieldIDsToLabels["AnglesNotObservedCountLabel"],
758 "%s"
759 % (
760 "NA"
761 if RotBondsAlertsInfo["AnglesNotObservedCount"] is None
762 else RotBondsAlertsInfo["AnglesNotObservedCount"]
763 ),
764 )
765
766 # Setup alert information for rotatable bonds...
767 AlertsInfoValues = []
768
769 # Delimiter for multiple values corresponding to specific set of information for
770 # a rotatable bond. For example: TorsionAtomIndices
771 ValuesDelim = OptionsInfo["IntraSetValuesDelim"]
772
773 # Delimiter for various values for a rotatable bond...
774 RotBondValuesDelim = OptionsInfo["InterSetValuesDelim"]
775
776 # Delimiter for values corresponding to multiple rotatable bonds...
777 AlertsInfoValuesDelim = OptionsInfo["InterSetValuesDelim"]
778
779 if RotBondsAlertsInfo is not None:
780 for ID in RotBondsAlertsInfo["IDs"]:
781 if not RotBondsAlertsInfo["MatchStatus"][ID]:
782 continue
783
784 if SkipRotatableBondAlertInfo(ID, RotBondsAlertsInfo):
785 continue
786
787 RotBondValues = []
788
789 # Bond atom indices...
790 Values = ["%s" % Value for Value in RotBondsAlertsInfo["AtomIndices"][ID]]
791 RotBondValues.append(ValuesDelim.join(Values))
792
793 # Torsion atom indices...
794 TorsionAtomIndices = SetupTorsionAtomIndicesValues(
795 RotBondsAlertsInfo["TorsionAtomIndices"][ID], ValuesDelim
796 )
797 RotBondValues.append(TorsionAtomIndices)
798
799 # Torsion angle...
800 RotBondValues.append("%.2f" % RotBondsAlertsInfo["TorsionAngle"][ID])
801
802 # Energy along with its lower and upper bound confidence interval...
803 RotBondValues.append(SetupEnergyValueForSDField(RotBondsAlertsInfo["Energy"][ID], Precision))
804 RotBondValues.append(SetupEnergyValueForSDField(RotBondsAlertsInfo["EnergyLowerBound"][ID], Precision))
805 RotBondValues.append(SetupEnergyValueForSDField(RotBondsAlertsInfo["EnergyUpperBound"][ID], Precision))
806
807 # Hierarchy class and subclass names...
808 RotBondValues.append("%s" % RotBondsAlertsInfo["HierarchyClassName"][ID])
809 RotBondValues.append("%s" % RotBondsAlertsInfo["HierarchySubClassName"][ID])
810
811 # Torsion rule SMARTS...
812 RotBondValues.append("%s" % RotBondsAlertsInfo["TorsionRuleSMARTS"][ID])
813
814 # Energy method...
815 RotBondValues.append("%s" % RotBondsAlertsInfo["EnergyMethod"][ID])
816
817 # Angle not observed...
818 RotBondValues.append("%s" % SetupAngleNotObservedValue(RotBondsAlertsInfo["AngleNotObserved"][ID]))
819
820 # Max single energy alert status...
821 RotBondValues.append(
822 "%s" % SetupMaxSingleEnergyAlertStatusValue(RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID])
823 )
824
825 # Track joined values for a rotatable bond...
826 AlertsInfoValues.append("%s" % RotBondValuesDelim.join(RotBondValues))
827
828 if len(AlertsInfoValues):
829 Mol.SetProp(
830 OptionsInfo["SDFieldIDsToLabels"]["TorsionAlertsLabel"],
831 "%s" % ("%s" % AlertsInfoValuesDelim.join(AlertsInfoValues)),
832 )
833
834
835 def WriteMoleculeFilteredByRuleID(Writer, Mol, TorsionRuleID, TorsionAlertsSummaryInfo, TorsionAlertsInfo):
836 """Write out molecule."""
837
838 if OptionsInfo["CountMode"]:
839 return
840
841 SetupMolPropertiesForFilteredByRuleIDAlertsInformation(
842 Mol, TorsionRuleID, TorsionAlertsSummaryInfo, TorsionAlertsInfo
843 )
844
845 Writer.write(Mol)
846
847
848 def SetupMolPropertiesForFilteredByRuleIDAlertsInformation(
849 Mol, TorsionRuleID, TorsionAlertsSummaryInfo, TorsionAlertsInfo
850 ):
851 """Setup molecule properties containing alerts information for torsion alerts
852 fileted by Rule IDs."""
853
854 # Delete torsion alerts information for rotatable bonds...
855 if Mol.HasProp(OptionsInfo["SDFieldIDsToLabels"]["TorsionAlertsLabel"]):
856 Mol.ClearProp(OptionsInfo["SDFieldIDsToLabels"]["TorsionAlertsLabel"])
857
858 # Delimiter for values...
859 IntraSetValuesDelim = OptionsInfo["IntraSetValuesDelim"]
860 InterSetValuesDelim = OptionsInfo["InterSetValuesDelim"]
861
862 # Setup alert rule information...
863 AlertRuleInfoValues = []
864
865 AlertRuleInfoValues.append("%s" % TorsionAlertsInfo["HierarchyClassName"][TorsionRuleID])
866 AlertRuleInfoValues.append("%s" % TorsionAlertsInfo["HierarchySubClassName"][TorsionRuleID])
867
868 AlertRuleInfoValues.append("%s" % TorsionAlertsInfo["RuleSMARTS"][TorsionRuleID])
869 AlertRuleInfoValues.append("%s" % TorsionAlertsInfo["EnergyMethod"][TorsionRuleID])
870
871 Mol.SetProp(
872 OptionsInfo["SDFieldIDsToLabels"]["TorsionRuleLabel"],
873 "%s" % ("%s" % InterSetValuesDelim.join(AlertRuleInfoValues)),
874 )
875
876 # Setup max single energy alert count for torsion rule...
877 Mol.SetProp(
878 OptionsInfo["SDFieldIDsToLabels"]["TorsionRuleMaxSingleEnergyAlertsCountLabel"],
879 "%s" % TorsionAlertsInfo["MaxSingleEnergyAlertsCount"][TorsionRuleID],
880 )
881
882 # Setup angle not observed count for torsion rule...
883 Mol.SetProp(
884 OptionsInfo["SDFieldIDsToLabels"]["TorsionRuleAnglesNotObservedCountLabel"],
885 "%s" % TorsionAlertsInfo["AnglesNotObservedCount"][TorsionRuleID],
886 )
887
888 # Setup torsion rule alerts...
889 # "TorsionRuleAlertsLabel": "TorsionRuleAlerts (RotBondIndices TorsionIndices TorsionAngle Energy EnergyLowerBoundCI EnergyUpperBoundCI EnergyMethod AngleNotObserved MaxSingleEnergyAlert)
890 AlertsInfoValues = []
891 for Index in range(0, len(TorsionAlertsInfo["AtomIndices"][TorsionRuleID])):
892 RotBondInfoValues = []
893
894 # Bond atom indices...
895 Values = ["%s" % Value for Value in TorsionAlertsInfo["AtomIndices"][TorsionRuleID][Index]]
896 RotBondInfoValues.append(IntraSetValuesDelim.join(Values))
897
898 # Torsion atom indices retrieved from the filtered SD file and stored as strings...
899 Values = ["%s" % Value for Value in TorsionAlertsInfo["TorsionAtomIndices"][TorsionRuleID][Index]]
900 RotBondInfoValues.append(IntraSetValuesDelim.join(Values))
901
902 # Torsion angle...
903 RotBondInfoValues.append(TorsionAlertsInfo["TorsionAngle"][TorsionRuleID][Index])
904
905 # Energy and its bounds...
906 RotBondInfoValues.append(TorsionAlertsInfo["Energy"][TorsionRuleID][Index])
907 RotBondInfoValues.append(TorsionAlertsInfo["EnergyLowerBound"][TorsionRuleID][Index])
908 RotBondInfoValues.append(TorsionAlertsInfo["EnergyUpperBound"][TorsionRuleID][Index])
909
910 # Angle not observed......
911 RotBondInfoValues.append(TorsionAlertsInfo["AngleNotObserved"][TorsionRuleID][Index])
912
913 # Max single energy alert type...
914 RotBondInfoValues.append(TorsionAlertsInfo["MaxSingleEnergyAlertType"][TorsionRuleID][Index])
915
916 # Track alerts informaiton...
917 AlertsInfoValues.append("%s" % InterSetValuesDelim.join(RotBondInfoValues))
918
919 Mol.SetProp(
920 OptionsInfo["SDFieldIDsToLabels"]["TorsionRuleAlertsLabel"], "%s" % (InterSetValuesDelim.join(AlertsInfoValues))
921 )
922
923
924 def SkipRotatableBondAlertInfo(ID, RotBondsAlertsInfo):
925 """Skip rotatble bond alert info for a specific bond during writing to output files."""
926
927 if not OptionsInfo["OutfileAlertsOnly"]:
928 return False
929
930 if RotBondsAlertsInfo["RotBondsAlertsStatus"] is None:
931 return True
932
933 Status = False
934 if OptionsInfo["TotalEnergyMode"]:
935 if not RotBondsAlertsInfo["RotBondsAlertsStatus"]:
936 Status = True
937 elif OptionsInfo["MaxSingleEnergyMode"]:
938 if (
939 RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID] is None
940 or not RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID]
941 ):
942 Status = True
943 elif OptionsInfo["TotalOrMaxSingleEnergyMode"]:
944 if not RotBondsAlertsInfo["RotBondsAlertsStatus"]:
945 Status = True
946
947 return Status
948
949
950 def SetupEnergyValueForSDField(Value, Precision):
951 """Setup energy value for SD field."""
952
953 if Value is None or math.isnan(Value) or math.isinf(Value):
954 return "NA"
955
956 return "%.*f" % (Precision, Value)
957
958
959 def SetupAngleNotObservedValue(Value):
960 """Setup angle not observed value."""
961
962 if Value is None:
963 return "NA"
964
965 return "Yes" if Value else "No"
966
967
968 def SetupMaxSingleEnergyAlertStatusValue(Value):
969 """Setup max single energy alert status value."""
970
971 if Value is None:
972 return "NA"
973
974 return "Yes" if Value else "No"
975
976
977 def SetupTorsionAtomIndicesValues(TorsionAtomIndicesList, ValuesDelim):
978 """Setup torsion atom indices value for output files."""
979
980 Values = ["%s" % Value for Value in TorsionAtomIndicesList]
981
982 return ValuesDelim.join(Values)
983
984
985 def SetupOutfilesWriters():
986 """Setup molecule and summary writers."""
987
988 OutfilesWriters = {"WriterRemaining": None, "WriterFiltered": None, "WriterAlertSummary": None}
989
990 # Writers for SD files...
991 WriterRemaining, WriterFiltered = SetupMoleculeWriters()
992 OutfilesWriters["WriterRemaining"] = WriterRemaining
993 OutfilesWriters["WriterFiltered"] = WriterFiltered
994
995 # Writer for alert summary CSV file...
996 WriterAlertSummary = SetupAlertSummaryWriter()
997 OutfilesWriters["WriterAlertSummary"] = WriterAlertSummary
998
999 return OutfilesWriters
1000
1001
1002 def SetupMoleculeWriters():
1003 """Setup molecule writers."""
1004
1005 Writer = None
1006 WriterFiltered = None
1007
1008 if OptionsInfo["CountMode"]:
1009 return (Writer, WriterFiltered)
1010
1011 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
1012 if Writer is None:
1013 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
1014 MiscUtil.PrintInfo("\nGenerating file %s..." % OptionsInfo["Outfile"])
1015
1016 if OptionsInfo["OutfileFilteredMode"]:
1017 WriterFiltered = RDKitUtil.MoleculesWriter(OptionsInfo["OutfileFiltered"], **OptionsInfo["OutfileParams"])
1018 if WriterFiltered is None:
1019 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["OutfileFiltered"])
1020 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["OutfileFiltered"])
1021
1022 return (Writer, WriterFiltered)
1023
1024
1025 def SetupAlertSummaryWriter():
1026 """Setup a alert summary writer."""
1027
1028 Writer = None
1029
1030 if OptionsInfo["CountMode"]:
1031 return Writer
1032
1033 if not OptionsInfo["OutfileSummaryMode"]:
1034 return Writer
1035
1036 Outfile = OptionsInfo["OutfileSummary"]
1037 Writer = open(Outfile, "w")
1038 if Writer is None:
1039 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
1040
1041 MiscUtil.PrintInfo("Generating file %s..." % Outfile)
1042
1043 return Writer
1044
1045
1046 def CloseOutfilesWriters(OutfilesWriters):
1047 """Close outfile writers."""
1048
1049 for WriterType, Writer in OutfilesWriters.items():
1050 if Writer is not None:
1051 Writer.close()
1052
1053
1054 def SetupByRuleOutfilesWriters(RuleIDs):
1055 """Setup by rule outfiles writers."""
1056
1057 # Initialize...
1058 OutfilesWriters = {}
1059 for RuleID in RuleIDs:
1060 OutfilesWriters[RuleID] = None
1061
1062 if OptionsInfo["CountMode"]:
1063 return OutfilesWriters
1064
1065 if not OptionsInfo["OutfilesFilteredByRulesMode"]:
1066 return OutfilesWriters
1067
1068 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1069 OutfilesRoot = "%s_Filtered_TopRule" % FileName
1070 OutfilesExt = "sdf"
1071
1072 MsgTxt = (
1073 "all"
1074 if OptionsInfo["OutfilesFilteredByRulesAllMode"]
1075 else "top %s" % OptionsInfo["OutfilesFilteredByRulesMaxCount"]
1076 )
1077 MiscUtil.PrintInfo(
1078 "\nGenerating output files %s*.%s for %s torsion rules triggering alerts..."
1079 % (OutfilesRoot, OutfilesExt, MsgTxt)
1080 )
1081
1082 # Delete any existing output files...
1083 Outfiles = glob.glob("%s*.%s" % (OutfilesRoot, OutfilesExt))
1084 if len(Outfiles):
1085 MiscUtil.PrintInfo("Deleting existing output files %s*.%s..." % (OutfilesRoot, OutfilesExt))
1086 for Outfile in Outfiles:
1087 try:
1088 os.remove(Outfile)
1089 except Exception as ErrMsg:
1090 MiscUtil.PrintWarning("Failed to delete file: %s" % ErrMsg)
1091
1092 RuleIndex = 0
1093 for RuleID in RuleIDs:
1094 RuleIndex += 1
1095 Outfile = "%s%s.%s" % (OutfilesRoot, RuleIndex, OutfilesExt)
1096 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
1097 if Writer is None:
1098 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
1099
1100 OutfilesWriters[RuleID] = Writer
1101
1102 return OutfilesWriters
1103
1104
1105 def CloseByRuleOutfilesWriters(OutfilesWriters):
1106 """Close by rule outfile writers."""
1107
1108 for RuleID, Writer in OutfilesWriters.items():
1109 if Writer is not None:
1110 Writer.close()
1111
1112
1113 def InstantiateTorsionStrainEnergyAlertsClass(Quiet=False):
1114 """Initialize torsion strain energy alerts class."""
1115
1116 try:
1117 TorsionStrainEnergyAlertsHandle = TorsionStrainEnergyAlerts(
1118 AlertsMode=OptionsInfo["AlertsMode"],
1119 TotalEnergyCutoff=OptionsInfo["TotalEnergyCutoff"],
1120 MaxSingleEnergyCutoff=OptionsInfo["MaxSingleEnergyCutoff"],
1121 RotBondsSMARTSMode=OptionsInfo["RotBondsSMARTSMode"],
1122 RotBondsSMARTSPattern=OptionsInfo["RotBondsSMARTSPattern"],
1123 TorsionLibraryFilePath=OptionsInfo["TorsionEnergyLibraryFile"],
1124 AlertTorsionsNotObserved=OptionsInfo["FilterTorsionsNotObserved"],
1125 )
1126 except Exception as ErrMsg:
1127 MiscUtil.PrintError("Failed to instantiate TorsionStrainEnergyAlerts:\n%s\n" % (ErrMsg))
1128
1129 if not Quiet:
1130 MiscUtil.PrintInfo(
1131 "\nRetrieving data from torsion strain energy library file %s..."
1132 % TorsionStrainEnergyAlertsHandle.GetTorsionLibraryFilePath()
1133 )
1134 TorsionStrainEnergyAlertsHandle.ListTorsionLibraryInfo()
1135
1136 return TorsionStrainEnergyAlertsHandle
1137
1138
1139 def ProcessRotatableBondsSMARTSMode():
1140 """ "Process SMARTS pattern for rotatable bonds."""
1141
1142 RotBondsMode = OptionsInfo["RotBondsSMARTSMode"]
1143
1144 RotBondsSMARTSPattern = None
1145 RotBondsSMARTSPatternSpecified = OptionsInfo["RotBondsSMARTSPatternSpecified"]
1146
1147 if re.match("^(NonStrict|SemiStrict|Strict)$", RotBondsMode, re.I):
1148 RotBondsSMARTSPattern = None
1149 elif re.match("Specify", RotBondsMode, re.I):
1150 RotBondsSMARTSPatternSpecified = RotBondsSMARTSPatternSpecified.strip()
1151 if not len(RotBondsSMARTSPatternSpecified):
1152 MiscUtil.PrintError(
1153 'Empty value specified for SMILES/SMARTS pattern in "--rotBondsSMARTSPattern" option, %s.'
1154 % RotBondsMode
1155 )
1156
1157 RotBondsPatternMol = Chem.MolFromSmarts(RotBondsSMARTSPatternSpecified)
1158 if RotBondsPatternMol is None:
1159 MiscUtil.PrintError(
1160 'Failed to create rotatable bonds pattern molecule. The rotatable bonds SMARTS pattern, "%s", specified using "--rotBondsSMARTSPattern" option is not valid.'
1161 % (RotBondsSMARTSPatternSpecified)
1162 )
1163 else:
1164 MiscUtil.PrintError(
1165 'The value, %s, specified for option "-r, --rotBondsSMARTSMode" is not valid. ' % RotBondsMode
1166 )
1167
1168 OptionsInfo["RotBondsSMARTSPattern"] = RotBondsSMARTSPattern
1169
1170
1171 def ProcessSDFieldLabelsOption():
1172 """Process SD data field label option."""
1173
1174 ParamsOptionName = "--outfileSDFieldLabels"
1175 ParamsOptionValue = Options["--outfileSDFieldLabels"]
1176
1177 ParamsIDsToLabels = {
1178 "RotBondsCountLabel": "RotBondsCount",
1179 "TotalEnergyLabel": "TotalEnergy",
1180 "TotalEnergyLowerBoundCILabel": "TotalEnergyLowerBoundCI",
1181 "TotalEnergyUpperBoundCILabel": "TotalEnergyUpperBoundCI",
1182 "MaxSingleEnergyLabel": "MaxSingleEnergy",
1183 "MaxSingleEnergyAlertsCountLabel": "MaxSingleEnergyAlertsCount",
1184 "AnglesNotObservedCountLabel": "AnglesNotObservedCount",
1185 "TorsionAlertsLabel": "TorsionAlerts(RotBondIndices TorsionIndices TorsionAngle Energy EnergyLowerBoundCI EnergyUpperBoundCI HierarchyClass HierarchySubClass TorsionRule EnergyMethod AngleNotObserved MaxSingleEnergyAlert)",
1186 "TorsionRuleLabel": "TorsionRule (HierarchyClass HierarchySubClass TorsionRule EnergyMethod)",
1187 "TorsionRuleMaxSingleEnergyAlertsCountLabel": "TorsionRuleMaxSingleEnergyAlertsCount",
1188 "TorsionRuleAnglesNotObservedCountLabel": "TorsionRuleAnglesNotObservedCount",
1189 "TorsionRuleAlertsLabel": "TorsionRuleAlerts (RotBondIndices TorsionIndices TorsionAngle Energy EnergyLowerBoundCI EnergyUpperBoundCI AngleNotObserved MaxSingleEnergyAlert)",
1190 }
1191
1192 if re.match("^auto$", ParamsOptionValue, re.I):
1193 OptionsInfo["SDFieldIDsToLabels"] = ParamsIDsToLabels
1194 return
1195
1196 # Setup a canonical paramater names...
1197 ValidParamNames = []
1198 CanonicalParamNamesMap = {}
1199 for ParamName in sorted(ParamsIDsToLabels):
1200 ValidParamNames.append(ParamName)
1201 CanonicalParamNamesMap[ParamName.lower()] = ParamName
1202
1203 ParamsOptionValue = ParamsOptionValue.strip()
1204 if not ParamsOptionValue:
1205 MiscUtil.PrintError('No valid parameter name and value pairs specified using "%s" option' % ParamsOptionName)
1206
1207 ParamsOptionValueWords = ParamsOptionValue.split(",")
1208 if len(ParamsOptionValueWords) % 2:
1209 MiscUtil.PrintError(
1210 'The number of comma delimited paramater names and values, %d, specified using "%s" option must be an even number.'
1211 % (len(ParamsOptionValueWords), ParamsOptionName)
1212 )
1213
1214 # Validate paramater name and value pairs...
1215 for Index in range(0, len(ParamsOptionValueWords), 2):
1216 Name = ParamsOptionValueWords[Index].strip()
1217 Value = ParamsOptionValueWords[Index + 1].strip()
1218
1219 CanonicalName = Name.lower()
1220 if CanonicalName not in CanonicalParamNamesMap:
1221 MiscUtil.PrintError(
1222 'The parameter name, %s, specified using "%s" is not a valid name. Supported parameter names: %s'
1223 % (Name, ParamsOptionName, " ".join(ValidParamNames))
1224 )
1225
1226 ParamName = CanonicalParamNamesMap[CanonicalName]
1227 ParamValue = Value
1228
1229 # Set value...
1230 ParamsIDsToLabels[ParamName] = ParamValue
1231
1232 OptionsInfo["SDFieldIDsToLabels"] = ParamsIDsToLabels
1233
1234
1235 def ProcessOptions():
1236 """Process and validate command line arguments and options."""
1237
1238 MiscUtil.PrintInfo("Processing options...")
1239
1240 # Validate options...
1241 ValidateOptions()
1242
1243 TotalEnergyMode, MaxSingleEnergyMode, TotalOrMaxSingleEnergyMode = [False] * 3
1244 AlertsMode = Options["--alertsMode"]
1245 if re.match("^TotalEnergy$", AlertsMode, re.I):
1246 TotalEnergyMode = True
1247 elif re.match("^MaxSingleEnergy$", AlertsMode, re.I):
1248 MaxSingleEnergyMode = True
1249 elif re.match("^TotalOrMaxSingleEnergy$", AlertsMode, re.I):
1250 TotalOrMaxSingleEnergyMode = True
1251 OptionsInfo["AlertsMode"] = AlertsMode
1252 OptionsInfo["TotalEnergyMode"] = TotalEnergyMode
1253 OptionsInfo["MaxSingleEnergyMode"] = MaxSingleEnergyMode
1254 OptionsInfo["TotalOrMaxSingleEnergyMode"] = TotalOrMaxSingleEnergyMode
1255
1256 OptionsInfo["FilterTorsionsNotObserved"] = (
1257 True if re.match("^yes$", Options["--filterTorsionsNotObserved"], re.I) else False
1258 )
1259
1260 OptionsInfo["MaxSingleEnergyCutoff"] = float(Options["--alertsMaxSingleEnergyCutoff"])
1261 OptionsInfo["TotalEnergyCutoff"] = float(Options["--alertsTotalEnergyCutoff"])
1262
1263 OptionsInfo["Infile"] = Options["--infile"]
1264 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
1265 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
1266 "--infileParams",
1267 Options["--infileParams"],
1268 InfileName=Options["--infile"],
1269 ParamsDefaultInfo=ParamsDefaultInfoOverride,
1270 )
1271
1272 OptionsInfo["Outfile"] = Options["--outfile"]
1273 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
1274 "--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"]
1275 )
1276
1277 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1278 OutfileFiltered = "%s_Filtered.%s" % (FileName, FileExt)
1279 OptionsInfo["OutfileFiltered"] = OutfileFiltered
1280 OptionsInfo["OutfileFilteredMode"] = True if re.match("^yes$", Options["--outfileFiltered"], re.I) else False
1281
1282 OptionsInfo["OutfileSummary"] = "%s_AlertsSummary.csv" % (FileName)
1283
1284 OutfileSummaryMode = Options["--outfileSummary"]
1285 if re.match("^auto$", OutfileSummaryMode, re.I):
1286 OutfileSummaryMode = "yes" if re.match("^MaxSingleEnergy$", Options["--alertsMode"], re.I) else "no"
1287 OptionsInfo["OutfileSummaryMode"] = True if re.match("^yes$", OutfileSummaryMode, re.I) else False
1288
1289 if re.match("^yes$", Options["--outfileSummary"], re.I):
1290 if not re.match("^MaxSingleEnergy$", Options["--alertsMode"], re.I):
1291 MiscUtil.PrintError(
1292 'The value "%s" specified for "--outfileSummary" option is not valid. The specified value is only allowed during "MaxSingleEnergy" value of "-a, --alertsMode" option.'
1293 % (Options["--outfileSummary"])
1294 )
1295
1296 OutfilesFilteredByRulesMode = Options["--outfilesFilteredByRules"]
1297 if re.match("^auto$", OutfilesFilteredByRulesMode, re.I):
1298 OutfilesFilteredByRulesMode = "yes" if re.match("^MaxSingleEnergy$", Options["--alertsMode"], re.I) else "no"
1299 OptionsInfo["OutfilesFilteredByRulesMode"] = True if re.match("^yes$", OutfilesFilteredByRulesMode, re.I) else False
1300
1301 if re.match("^yes$", Options["--outfilesFilteredByRules"], re.I):
1302 if not re.match("^MaxSingleEnergy$", Options["--alertsMode"], re.I):
1303 MiscUtil.PrintError(
1304 'The value "%s" specified for "--outfilesFilteredByRules" option is not valid. The specified value is only allowed during "MaxSingleEnergy" value of "-a, --alertsMode" option.'
1305 % (Options["--outfileSummary"])
1306 )
1307
1308 OptionsInfo["TrackAlertsSummaryInfo"] = (
1309 True if (OptionsInfo["OutfileSummaryMode"] or OptionsInfo["OutfilesFilteredByRulesMode"]) else False
1310 )
1311
1312 OutfilesFilteredByRulesMaxCount = Options["--outfilesFilteredByRulesMaxCount"]
1313 if not re.match("^All$", OutfilesFilteredByRulesMaxCount, re.I):
1314 OutfilesFilteredByRulesMaxCount = int(OutfilesFilteredByRulesMaxCount)
1315 OptionsInfo["OutfilesFilteredByRulesMaxCount"] = OutfilesFilteredByRulesMaxCount
1316 OptionsInfo["OutfilesFilteredByRulesAllMode"] = (
1317 True if re.match("^All$", Options["--outfilesFilteredByRulesMaxCount"], re.I) else False
1318 )
1319
1320 OptionsInfo["OutfileAlerts"] = True if re.match("^yes$", Options["--outfileAlerts"], re.I) else False
1321
1322 if re.match("^yes$", Options["--outfilesFilteredByRules"], re.I):
1323 if not re.match("^yes$", Options["--outfileAlerts"], re.I):
1324 MiscUtil.PrintError(
1325 'The value "%s" specified for "--outfilesFilteredByRules" option is not valid. The specified value is only allowed during "yes" value of "--outfileAlerts" option.'
1326 % (Options["--outfilesFilteredByRules"])
1327 )
1328
1329 OptionsInfo["OutfileAlertsMode"] = Options["--outfileAlertsMode"]
1330 OptionsInfo["OutfileAlertsOnly"] = True if re.match("^AlertsOnly$", Options["--outfileAlertsMode"], re.I) else False
1331
1332 ProcessSDFieldLabelsOption()
1333
1334 OptionsInfo["Overwrite"] = Options["--overwrite"]
1335 OptionsInfo["CountMode"] = True if re.match("^count$", Options["--mode"], re.I) else False
1336
1337 OptionsInfo["Precision"] = int(Options["--precision"])
1338
1339 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1340 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1341
1342 OptionsInfo["RotBondsSMARTSMode"] = Options["--rotBondsSMARTSMode"]
1343 OptionsInfo["RotBondsSMARTSPatternSpecified"] = Options["--rotBondsSMARTSPattern"]
1344 ProcessRotatableBondsSMARTSMode()
1345
1346 OptionsInfo["TorsionEnergyLibraryFile"] = Options["--torsionEnergyLibraryFile"]
1347
1348 # Setup delimiter for writing out torsion alert information to output files...
1349 OptionsInfo["IntraSetValuesDelim"] = ","
1350 OptionsInfo["InterSetValuesDelim"] = " "
1351
1352
1353 def RetrieveOptions():
1354 """Retrieve command line arguments and options."""
1355
1356 # Get options...
1357 global Options
1358 Options = docopt(_docoptUsage_)
1359
1360 # Set current working directory to the specified directory...
1361 WorkingDir = Options["--workingdir"]
1362 if WorkingDir:
1363 os.chdir(WorkingDir)
1364
1365 # Handle examples option...
1366 if "--examples" in Options and Options["--examples"]:
1367 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1368 sys.exit(0)
1369
1370
1371 def ProcessListTorsionLibraryOption():
1372 """Process list torsion library information."""
1373
1374 # Validate and process dataFile option for listing torsion library information...
1375 OptionsInfo["TorsionEnergyLibraryFile"] = Options["--torsionEnergyLibraryFile"]
1376 if not re.match("^auto$", Options["--torsionEnergyLibraryFile"], re.I):
1377 MiscUtil.ValidateOptionFilePath("-t, --torsionEnergyLibraryFile", Options["--torsionEnergyLibraryFile"])
1378
1379 # Instantiate TorsionStrainEnergyAlerts using defaults...
1380 TorsionStrainEnergyAlertsHandle = TorsionStrainEnergyAlerts(
1381 TorsionLibraryFilePath=OptionsInfo["TorsionEnergyLibraryFile"]
1382 )
1383 MiscUtil.PrintInfo(
1384 "\nRetrieving data from torsion strain energy library file %s..."
1385 % TorsionStrainEnergyAlertsHandle.GetTorsionLibraryFilePath()
1386 )
1387 TorsionStrainEnergyAlertsHandle.ListTorsionLibraryInfo()
1388
1389
1390 def ValidateOptions():
1391 """Validate option values."""
1392
1393 MiscUtil.ValidateOptionTextValue(
1394 "-a, --alertsMode", Options["--alertsMode"], "TotalEnergy MaxSingleEnergy TotalOrMaxSingleEnergy"
1395 )
1396
1397 MiscUtil.ValidateOptionFloatValue(
1398 "--alertsMaxSingleEnergyCutoff", Options["--alertsMaxSingleEnergyCutoff"], {">": 0.0}
1399 )
1400 MiscUtil.ValidateOptionFloatValue("--alertsTotalEnergyCutoff", Options["--alertsTotalEnergyCutoff"], {">": 0.0})
1401
1402 MiscUtil.ValidateOptionTextValue("--filterTorsionsNotObserved", Options["--filterTorsionsNotObserved"], "yes no")
1403
1404 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1405 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
1406
1407 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
1408 if re.match("^filter$", Options["--mode"], re.I):
1409 MiscUtil.ValidateOptionsOutputFileOverwrite(
1410 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
1411 )
1412 MiscUtil.ValidateOptionsDistinctFileNames(
1413 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
1414 )
1415
1416 MiscUtil.ValidateOptionTextValue("--outfileFiltered", Options["--outfileFiltered"], "yes no")
1417
1418 MiscUtil.ValidateOptionTextValue("--outfilesFilteredByRules", Options["--outfilesFilteredByRules"], "yes no auto")
1419 if not re.match("^All$", Options["--outfilesFilteredByRulesMaxCount"], re.I):
1420 MiscUtil.ValidateOptionIntegerValue(
1421 "--outfilesFilteredByRulesMaxCount", Options["--outfilesFilteredByRulesMaxCount"], {">": 0}
1422 )
1423
1424 MiscUtil.ValidateOptionTextValue("--outfileSummary", Options["--outfileSummary"], "yes no auto")
1425 MiscUtil.ValidateOptionTextValue("--outfileAlerts", Options["--outfileAlerts"], "yes no")
1426 MiscUtil.ValidateOptionTextValue("--outfileAlertsMode", Options["--outfileAlertsMode"], "All AlertsOnly")
1427
1428 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "filter count")
1429 if re.match("^filter$", Options["--mode"], re.I):
1430 if not Options["--outfile"]:
1431 MiscUtil.PrintError(
1432 'The outfile must be specified using "-o, --outfile" during "filter" value of "-m, --mode" option'
1433 )
1434
1435 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1436
1437 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
1438
1439 MiscUtil.ValidateOptionTextValue(
1440 "-r, --rotBondsSMARTSMode", Options["--rotBondsSMARTSMode"], "NonStrict SemiStrict Strict Specify"
1441 )
1442 if re.match("^Specify$", Options["--rotBondsSMARTSMode"], re.I):
1443 if not Options["--rotBondsSMARTSPattern"]:
1444 MiscUtil.PrintError(
1445 'The SMARTS pattern must be specified using "--rotBondsSMARTSPattern" during "Specify" value of "-r, --rotBondsSMARTS" option'
1446 )
1447 else:
1448 if Options["--rotBondsSMARTSPattern"]:
1449 MiscUtil.PrintError(
1450 'The SMARTS pattern must not be specified using "--rotBondsSMARTSPattern" during "%s" value of "-r, --rotBondsSMARTS" option'
1451 % (Options["--rotBondsSMARTSMode"])
1452 )
1453
1454 if not re.match("^auto$", Options["--torsionEnergyLibraryFile"], re.I):
1455 MiscUtil.ValidateOptionFilePath("-t, --torsionEnergyLibraryFile", Options["--torsionEnergyLibraryFile"])
1456
1457
1458 # Setup a usage string for docopt...
1459 _docoptUsage_ = """
1460 RDKitFilterTorsionStrainEnergyAlerts.py - Filter torsion strain energy library alerts
1461
1462 Usage:
1463 RDKitFilterTorsionStrainEnergyAlerts.py [--alertsMode <TotalEnergy, MaxSingleEnergy, or TotalOrMaxSingleEnergy>]
1464 [--alertsMaxSingleEnergyCutoff <Number>] [--alertsTotalEnergyCutoff <Number>]
1465 [--filterTorsionsNotObserved <yes or no>] [--infileParams <Name,Value,...>] [--mode <filter or count>]
1466 [--mp <yes or no>] [--mpParams <Name,Value,...>] [--outfileAlerts <yes or no>]
1467 [--outfileAlertsMode <All or AlertsOnly>] [--outfileFiltered <yes or no>]
1468 [--outfilesFilteredByRules <yes or no>] [--outfilesFilteredByRulesMaxCount <All or number>]
1469 [--outfileSummary <yes or no>] [--outfileSDFieldLabels <Type,Label,...>] [--outfileParams <Name,Value,...>]
1470 [--overwrite] [--precision <number>] [ --rotBondsSMARTSMode <NonStrict, SemiStrict,...>]
1471 [--rotBondsSMARTSPattern <SMARTS>] [--torsionEnergyLibraryFile <FileName or auto>]
1472 [-w <dir>] -i <infile> -o <outfile>
1473 RDKitFilterTorsionStrainEnergyAlerts.py [--torsionEnergyLibraryFile <FileName or auto>] -l | --list
1474 RDKitFilterTorsionStrainEnergyAlerts.py -h | --help | -e | --examples
1475
1476 Description:
1477 Filter strained molecules from an input file for torsion strain energy library
1478 [ Ref 153 ] alerts by matching rotatable bonds against SMARTS patterns specified
1479 for torsion rules in a torsion energy library file and write out appropriate
1480 molecules to output files. The molecules must have 3D coordinates in input file.
1481 The default torsion strain energy library file, TorsionStrainEnergyLibrary.xml,
1482 is available under MAYACHEMTOOLS/lib/python/TorsionAlerts directory.
1483
1484 The data in torsion strain energy library file is organized in a hierarchical
1485 manner. It consists of one generic class and six specific classes at the highest
1486 level. Each class contains multiple subclasses corresponding to named functional
1487 groups or substructure patterns. The subclasses consist of torsion rules sorted
1488 from specific to generic torsion patterns. The torsion rule, in turn, contains a
1489 list of peak values for torsion angles and two tolerance values. A pair of tolerance
1490 values define torsion bins around a torsion peak value.
1491
1492 A strain energy calculation method, 'exact' or 'approximate' [ Ref 153 ], is
1493 associated with each torsion rule for calculating torsion strain energy. The 'exact'
1494 stain energy calculation relies on the energy bins available under the energy histogram
1495 consisting of 36 bins covering angles from -180 to 180. The width of each bin is 10
1496 degree. The energy bins are are defined at the right end points. The first and the
1497 last energy bins correspond to -170 and 180 respectively. The torsion angle is mapped
1498 to a energy bin. An angle offset is calculated for the torsion angle from the the right
1499 end point angle of the bin. The strain energy is estimated for the angle offset based
1500 on the energy difference between the current and previous bins. The torsion strain
1501 energy, in terms of torsion energy units (TEUs), corresponds to the sum of bin strain
1502 energy and the angle offset strain energy.
1503
1504 Energy = BinEnergyDiff/10.0 * BinAngleOffset + BinEnergy[BinNum]
1505
1506 Where:
1507
1508 BinEnergyDiff = BinEnergy[BinNum] - BinEnergy[PreviousBinNum]
1509 BinAngleOffset = TorsionAngle - BinAngleRightSide
1510
1511 The 'approximate' strain energy calculation relies on the angle difference between a
1512 torsion angle and the torsion peaks observed for the torsion rules in the torsion
1513 energy library. The torsion angle is matched to a torsion peak based on the value of
1514 torsion angle difference. It must be less than or equal to the value for the second
1515 tolerance 'tolerance2'. Otherwise, the torsion angle is not observed in the torsion
1516 energy library and a value of 'NA' is assigned for torsion energy along with the lower
1517 and upper bounds on energy at 95% confidence interval. The 'approximate' torsion
1518 energy (TEUs) for observed torsion angle is calculated using the following formula:
1519
1520 Energy = beta_1 * (AngleDiff ** 2) + beta_2 * (AngleDiff ** 4)
1521
1522 The coefficients 'beta_1' and 'beta_2' are available for the observed angles in
1523 the torsion strain energy library. The 'AngleDiff' is the difference between the
1524 torsion angle and the matched torsion peak.
1525
1526 For example:
1527
1528 <library>
1529 <hierarchyClass id1="G" id2="G" name="GG">
1530 ...
1531 </hierarchyClass>
1532 <hierarchyClass id1="C" id2="O" name="CO">
1533 <hierarchySubClass name="Ester bond I" smarts="O=[C:2][O:3]">
1534 <torsionRule method="exact" smarts=
1535 "[O:1]=[C:2]!@[O:3]~[CH0:4]">
1536 <angleList>
1537 <angle score="56.52" tolerance1="20.00"
1538 tolerance2="25.00" value="0.0"/>
1539 </angleList>
1540 <histogram>
1541 <bin count="1"/>
1542 ...
1543 </histogram>
1544 <histogram_shifted>
1545 <bin count="0"/>
1546 ...
1547 </histogram_shifted>
1548 <histogram_converted>
1549 <bin energy="4.67... lower="2.14..." upper="Inf"/>
1550 ...
1551 <bin energy="1.86..." lower="1.58..." upper="2.40..."/>
1552 ...
1553 </histogram_converted>
1554 </torsionRule>
1555 <torsionRule method="approximate" smarts=
1556 "[cH0:1][c:2]([cH0])!@[O:3][p:4]">
1557 <angleList>
1558 <angle beta_1="0.002..." beta_2="-7.843...e-07"
1559 score="27.14" theta_0="-90.0" tolerance1="30.00"
1560 tolerance2="45.00" value="-90.0"/>
1561 ...
1562 </angleList>
1563 <histogram>
1564 <bin count="0"/>
1565 ...
1566 </histogram>
1567 <histogram_shifted>
1568 <bin count="0"/>
1569 ...
1570 </histogram_shifted>
1571 </torsionRule>
1572 ...
1573 ...
1574 </hierarchyClass>
1575 <hierarchyClass id1="N" id2="C" name="NC">
1576 ...
1577 </hierarchyClass>
1578 <hierarchyClass id1="S" id2="N" name="SN">
1579 ...
1580 </hierarchyClass>
1581 <hierarchyClass id1="C" id2="S" name="CS">
1582 ...
1583 </hierarchyClass>
1584 <hierarchyClass id1="C" id2="C" name="CC">
1585 ...
1586 </hierarchyClass>
1587 <hierarchyClass id1="S" id2="S" name="SS">
1588 ...
1589 </hierarchyClass>
1590 </library>
1591
1592 The rotatable bonds in a 3D molecule are identified using a default SMARTS pattern.
1593 A custom SMARTS pattern may be optionally specified to detect rotatable bonds.
1594 Each rotatable bond is matched to a torsion rule in the torsion strain energy library.
1595 The strain energy is calculated for each rotatable bond using the calculation
1596 method, 'exact' or 'approximate', associated with the matched torsion rule.
1597
1598 The total strain energy (TEUs) of a molecule corresponds to the sum of 'exact' and
1599 'approximate' strain energies calculated for all matched rotatable bonds in the
1600 molecule. The total strain energy is set to 'NA' for molecules containing a 'approximate'
1601 energy estimate for a torsion angle not observed in the torsion energy library. In
1602 addition, the lower and upper bounds on energy at 95% confidence interval are
1603 set to 'NA'.
1604
1605 The following output files are generated after the filtering:
1606
1607 <OutfileRoot>.sdf
1608 <OutfileRoot>_Filtered.sdf
1609 <OutfileRoot>_AlertsSummary.csv
1610 <OutfileRoot>_Filtered_TopRule*.sdf
1611
1612 The last two set of outfile files, <OutfileRoot>_AlertsSummary.csv and
1613 <OutfileRoot>_<OutfileRoot>_AlertsSummary.csv, are only generated during filtering
1614 by 'MaxSingleEnergy'.
1615
1616 The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
1617
1618 The supported output file formats are: SD (.sdf, .sd)
1619
1620 Options:
1621 -a, --alertsMode <TotalEnergy,...> [default: TotalEnergy]
1622 Torsion strain energy library alert types to use for filtering molecules
1623 containing rotatable bonds based on the calculated values for the total
1624 torsion strain energy of a molecule and the maximum single strain
1625 energy of a rotatable bond in a molecule.
1626
1627 Possible values: TotalEnergy, MaxSingleEnergy, or TotalOrMaxSingleEnergy
1628
1629 The strain energy cutoff values in terms of torsion energy units (TEUs) are
1630 used to filter molecules as shown below:
1631
1632 AlertsMode AlertsEnergyCutoffs (TEUs)
1633
1634 TotalEnergy >= TotalEnergyCutoff
1635
1636 MaxSingleEnergy >= MaxSingleEnergyCutoff
1637
1638 TotalOrMaxSingleEnergy >= TotalEnergyCutoff
1639 or >= MaxSingleEnergyCutoff
1640
1641 --alertsMaxSingleEnergyCutoff <Number> [default: 1.8]
1642 Maximum single strain energy (TEUs) cutoff [ Ref 153 ] for filtering molecules
1643 based on the maximum value of a single strain energy of a rotatable bond
1644 in a molecule. This option is used during 'MaxSingleEnergy' or
1645 'TotalOrMaxSingleEnergy' values of '-a, --alertsMode' option.
1646
1647 The maximum single strain energy must be greater than or equal to the
1648 specified cutoff value for filtering molecules.
1649 --alertsTotalEnergyCutoff <Number> [default: 6.0]
1650 Total strain strain energy (TEUs) cutoff [ Ref 153 ] for filtering molecules
1651 based on total strain energy for all rotatable bonds in a molecule. This
1652 option is used during 'TotalEnergy' or 'TotalOrMaxSingleEnergy'
1653 values of '-a, --alertsMode' option.
1654
1655 The total strain energy must be greater than or equal to the specified
1656 cutoff value for filtering molecules.
1657 --filterTorsionsNotObserved <yes or no> [default: no]
1658 Filter molecules containing torsion angles not observed in torsion strain
1659 energy library. It's not possible to calculate torsion strain energies for
1660 these torsions during 'approximate' match to a specified torsion in the
1661 library.
1662
1663 The 'approximate' strain energy calculation relies on the angle difference
1664 between a torsion angle and the torsion peaks observed for the torsion
1665 rules in the torsion energy library. The torsion angle is matched to a
1666 torsion peak based on the value of torsion angle difference. It must be
1667 less than or equal to the value for the second tolerance 'tolerance2'.
1668 Otherwise, the torsion angle is not observed in the torsion energy library
1669 and a value of 'NA' is assigned for torsion energy along with the lower and
1670 upper bounds on energy at 95% confidence interval.
1671 -e, --examples
1672 Print examples.
1673 -h, --help
1674 Print this help message.
1675 -i, --infile <infile>
1676 Input file name.
1677 --infileParams <Name,Value,...> [default: auto]
1678 A comma delimited list of parameter name and value pairs for reading
1679 molecules from files. The supported parameter names for different file
1680 formats, along with their default values, are shown below:
1681
1682 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
1683
1684 -l, --list
1685 List torsion library information without performing any filtering.
1686 -m, --mode <filter or count> [default: filter]
1687 Specify whether to filter molecules for torsion strain energy library [ Ref 153 ]
1688 alerts by matching rotatable bonds against SMARTS patterns specified for
1689 torsion rules to calculate torsion strain energies and write out the rest
1690 of the molecules to an outfile or simply count the number of matched
1691 molecules marked for filtering.
1692 --mp <yes or no> [default: no]
1693 Use multiprocessing.
1694
1695 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1696 function employing lazy RDKit data iterable. This allows processing of
1697 arbitrary large data sets without any additional requirements memory.
1698
1699 All input data may be optionally loaded into memory by mp.Pool.map()
1700 before starting worker processes in a process pool by setting the value
1701 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1702
1703 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1704 data mode may adversely impact the performance. The '--mpParams' section
1705 provides additional information to tune the value of 'chunkSize'.
1706 --mpParams <Name,Value,...> [default: auto]
1707 A comma delimited list of parameter name and value pairs to configure
1708 multiprocessing.
1709
1710 The supported parameter names along with their default and possible
1711 values are shown below:
1712
1713 chunkSize, auto
1714 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1715 numProcesses, auto [ Default: mp.cpu_count() ]
1716
1717 These parameters are used by the following functions to configure and
1718 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1719 mp.Pool.imap().
1720
1721 The chunkSize determines chunks of input data passed to each worker
1722 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1723 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1724
1725 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1726 automatically converts RDKit data iterable into a list, loads all data into
1727 memory, and calculates the default chunkSize using the following method
1728 as shown in its code:
1729
1730 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1731 if extra: chunkSize += 1
1732
1733 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1734 and 100 data items.
1735
1736 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1737 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1738 data into memory. Consequently, the size of input data is not known a priori.
1739 It's not possible to estimate an optimal value for the chunkSize. The default
1740 chunkSize is set to 1.
1741
1742 The default value for the chunkSize during 'Lazy' data mode may adversely
1743 impact the performance due to the overhead associated with exchanging
1744 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1745 a larger value during 'Lazy' input data mode, based on the size of your input
1746 data and number of processes in the process pool.
1747
1748 The mp.Pool.map() function waits for all worker processes to process all
1749 the data and return the results. The mp.Pool.imap() function, however,
1750 returns the the results obtained from worker processes as soon as the
1751 results become available for specified chunks of data.
1752
1753 The order of data in the results returned by both mp.Pool.map() and
1754 mp.Pool.imap() functions always corresponds to the input data.
1755 -o, --outfile <outfile>
1756 Output file name.
1757 --outfileAlerts <yes or no> [default: yes]
1758 Write out alerts information to SD output files.
1759 --outfileAlertsMode <All or AlertsOnly> [default: AlertsOnly]
1760 Write alerts information to SD output files for all alerts or only for alerts
1761 specified by '--AlertsMode' option. Possible values: All or AlertsOnly
1762 This option is only valid for 'Yes' value of '--outfileAlerts' option.
1763
1764 The following alerts information is added to SD output files using
1765 'TorsionAlerts' data field:
1766
1767 RotBondIndices TorsionIndices TorsionAngle
1768 Energy EnergyLowerBoundCI EnergyUpperBoundCI
1769 HierarchyClass HierarchySubClass TorsionRule
1770 EnergyMethod AngleNotObserved MaxSingleEnergyAlert
1771
1772 The following data filelds are added to SD output files based on the value of
1773 '--AlertsMode' option:
1774
1775 TotalEnergy
1776 TotalEnergyLowerBoundCI
1777 TotalEnergyUpperBoundCI
1778
1779 MaxSingleEnergy
1780 MaxSingleEnergyAlertsCount
1781
1782 AnglesNotObservedCount
1783
1784 The 'RotBondsCount' is always added to SD output files containing both
1785 remaining and filtered molecules.
1786
1787 Format:
1788
1789 > <RotBondsCount>
1790 Number
1791
1792 > <TotalEnergy>
1793 Number
1794
1795 > <TotalEnergyLowerBoundCI>
1796 Number
1797
1798 > <TotalEnergyUpperBoundCI>
1799 Number
1800
1801 > <MaxSingleEnergy>
1802 Number
1803
1804 > <MaxSingleEnergyAlertsCount>
1805 Number
1806
1807 > <AnglesNotObservedCount>
1808 Number
1809
1810 > <TorsionAlerts (RotBondIndices TorsionIndices TorsionAngle
1811 Energy EnergyLowerBoundCI EnergyUpperBoundCI
1812 HierarchyClass HierarchySubClass TorsionRule
1813 EnergyMethod AngleNotObserved MaxSingleEnergyAlert)>
1814 AtomIndex2,AtomIndex3 AtomIndex1,AtomIndex2,AtomIndex3,AtomIndex4
1815 Angle Energy EnergyLowerBoundCI EnergyUpperBoundCI
1816 ClassName SubClassName SMARTS EnergyMethod Yes|No|NA Yes|No|NA
1817 ... ... ...
1818 ... ... ...
1819
1820 A set of 12 values is written out as value of 'TorsionAlerts' data field for
1821 each torsion in a molecule. The space character is used as a delimiter
1822 to separate values with in a set and across set. The comma character
1823 is used to delimit multiple values for each value in a set.
1824
1825 The 'RotBondIndices' and 'TorsionIndices' contain 2 and 4 comma delimited
1826 values representing atom indices for a rotatable bond and the matched
1827 torsion.
1828
1829 The 'Energy' value is the estimated strain energy for the matched torsion.
1830 The 'EnergyLowerBoundCI' and 'EnergyUpperBoundCI' represent lower and
1831 bound energy estimates at 95% confidence interval. The 'EnergyMethod',
1832 exact or approximate, corresponds to the method employed to estimate
1833 torsion strain energy.
1834
1835 The 'AngleNotObserved' is only valid for 'approximate' value of 'EnergyMethod'.
1836 It has three possible values: Yes, No, or NA. The 'Yes' value indicates that
1837 the 'TorsionAngle' is outside the 'tolerance2' of all peaks for the matched
1838 torsion rule in the torsion library.
1839
1840 The 'MaxSingleEnergyAlert' is valid for the following values of '-a, --alertsMode'
1841 option: 'MaxSingleEnergy' or 'TotalOrMaxSingleEnergy'. It has three possible
1842 values: Yes, No, or NA. It's set to 'NA' for 'Yes' or 'NA' values of
1843 'AngleNotObserved'. The 'Yes' value indicates that the estimated torsion
1844 energy is greater than the specified value for '--alertsMaxSingleEnergyCutoff'
1845 option.
1846
1847 For example:
1848
1849 > <RotBondsCount> (1)
1850 14
1851
1852 > <TotalEnergy> (1)
1853 6.8065
1854
1855 > <TotalEnergyLowerBoundCI> (1)
1856 5.9340
1857
1858 > <TotalEnergyUpperBoundCI> (1)
1859 NA
1860
1861 > <MaxSingleEnergy> (1)
1862 1.7108
1863
1864 > <MaxSingleEnergyAlertsCount> (1)
1865 0
1866
1867 > <AnglesNotObservedCount> (1)
1868 0
1869
1870 > <TorsionAlerts(RotBondIndices TorsionIndices TorsionAngle Energy
1871 EnergyLowerBoundCI EnergyUpperBoundCI HierarchyClass
1872 HierarchySubClass TorsionRule EnergyMethod AngleNotObserved
1873 MaxSingleEnergyAlert)> (1)
1874 2,1 48,2,1,0 61.90 0.0159 -0.0320 0.0674 CO Ether [O:1][CX4:2]!
1875 @[O:3][CX4:4] Exact NA No 2,3 1,2,3,4 109.12 1.5640 1.1175 NA CC
1876 None/[CX4][CX3] [O:1][CX4:2]!@[CX3:3]=[O:4] Exact NA No
1877 ... ... ...
1878
1879 --outfileFiltered <yes or no> [default: yes]
1880 Write out a file containing filtered molecules. Its name is automatically
1881 generated from the specified output file. Default: <OutfileRoot>_
1882 Filtered.<OutfileExt>.
1883 --outfilesFilteredByRules <yes or no> [default: auto]
1884 Write out SD files containing filtered molecules for individual torsion
1885 rules triggering alerts in molecules. The name of SD files are automatically
1886 generated from the specified output file. Default file names: <OutfileRoot>_
1887 Filtered_TopRule*.sdf.
1888
1889 Default value: 'yes' for 'MaxSingleEnergy' of '-a, --alertsMode' option';
1890 otherwise, 'no'.
1891
1892 The output files are only generated for 'MaxSingleEnergy' of
1893 '-a, --alertsMode' option.
1894
1895 The following alerts information is added to SD output files:
1896
1897 > <RotBondsCount>
1898 Number
1899
1900 > <TotalEnergy>
1901 Number
1902
1903 > <TotalEnergyLowerBoundCI>
1904 Number
1905
1906 > <TotalEnergyUpperBoundCI>
1907 Number
1908
1909 > <MaxSingleEnergy>
1910 Number
1911
1912 > <MaxSingleEnergyAlertsCount>
1913 Number
1914
1915 > <AnglesNotObservedCount>
1916 Number
1917
1918 > <TorsionRule (HierarchyClass HierarchySubClass TorsionRule
1919 EnergyMethod)>
1920 ClassName SubClassName EnergyMethod SMARTS
1921 ... ... ...
1922
1923 > <TorsionRuleMaxSingleEnergyAlertsCount>
1924 Number
1925
1926 > <TorsionRuleAnglesNotObservedCount>
1927 Number
1928
1929 > <TorsionRuleAlerts (RotBondIndices TorsionIndices TorsionAngle
1930 Energy EnergyLowerBoundCI EnergyUpperBoundCI
1931 AngleNotObserved MaxSingleEnergyAlert)>
1932 AtomIndex2,AtomIndex3 AtomIndex1,AtomIndex2,AtomIndex3,AtomIndex4
1933 Angle Energy EnergyLowerBoundCI EnergyUpperBoundCI EnergyMethod
1934 Yes|No|NA Yes|No|NA
1935 ... ... ...
1936
1937 For example:
1938
1939 > <RotBondsCount> (1)
1940 8
1941
1942 > <TotalEnergy> (1)
1943 6.1889
1944
1945 > <TotalEnergyLowerBoundCI> (1)
1946 5.1940
1947
1948 > <TotalEnergyUpperBoundCI> (1)
1949 NA
1950
1951 > <MaxSingleEnergy> (1)
1952 1.9576
1953
1954 > <MaxSingleEnergyAlertsCount> (1)
1955 1
1956
1957 > <AnglesNotObservedCount> (1)
1958 0
1959
1960 > <TorsionRule (HierarchyClass HierarchySubClass TorsionRule
1961 EnergyMethod)> (1)
1962 CC None/[CX4:2][CX4:3] [!#1:1][CX4:2]!@[CX4:3][!#1:4] Exact
1963
1964 > <TorsionRuleMaxSingleEnergyAlertsCount> (1)
1965 0
1966
1967 > <TorsionRuleAnglesNotObservedCount> (1)
1968 0
1969
1970 > <TorsionRuleAlerts (RotBondIndices TorsionIndices TorsionAngle
1971 Energy EnergyLowerBoundCI EnergyUpperBoundCI AngleNotObserved
1972 MaxSingleEnergyAlert)> (1)
1973 1,3 0,1,3,4 72.63 0.8946 0.8756 0.9145 NA No
1974
1975 --outfilesFilteredByRulesMaxCount <All or number> [default: 10]
1976 Write out SD files containing filtered molecules for specified number of
1977 top N torsion rules triggering alerts for the largest number of molecules
1978 or for all torsion rules triggering alerts across all molecules.
1979
1980 These output files are only generated for 'MaxSingleEnergy' value of
1981 '-a, --alertsMode' option.
1982 --outfileSummary <yes or no> [default: auto]
1983 Write out a CVS text file containing summary of torsions rules responsible
1984 for triggering torsion alerts. Its name is automatically generated from the
1985 specified output file. Default: <OutfileRoot>_AlertsSummary.csv.
1986
1987 Default value: 'yes' for 'MaxSingleEnergy' of '-a, --alertsMode' option';
1988 otherwise, 'no'.
1989
1990 The summary output file is only generated for 'MaxSingleEnergy' of
1991 '-a, --alertsMode' option.
1992
1993 The following alerts information is written to summary text file:
1994
1995 TorsionRule, HierarchyClass, HierarchySubClass, EnergyMethod,
1996 MaxSingleEnergyTorsionAlertTypes, MaxSingleEnergyTorsionAlertCount,
1997 MaxSingleEnergyTorsionAlertMolCount
1998
1999 The double quotes characters are removed from SMART patterns before
2000 before writing them to a CSV file. In addition, the torsion rules are sorted by
2001 TorsionAlertMolCount.
2002 --outfileSDFieldLabels <Type,Label,...> [default: auto]
2003 A comma delimited list of SD data field type and label value pairs for writing
2004 torsion alerts information along with molecules to SD files.
2005
2006 The supported SD data field label type along with their default values are
2007 shown below:
2008
2009 For all SD files:
2010
2011 RotBondsCountLabel, RotBondsCount,
2012
2013 TotalEnergyLabel, TotalEnergy,
2014 TotalEnergyLowerBoundCILabel, TotalEnergyLowerBoundCI,
2015 TotalEnergyUpperBoundCILabel, TotalEnergyUpperBoundCI,
2016
2017 MaxSingleEnergyLabel, MaxSingleEnergy,
2018 MaxSingleEnergyAlertsCountLabel,
2019 MaxSingleEnergyAlertsCount
2020
2021 AnglesNotObservedCountLabel,
2022 AnglesNotObservedCount
2023
2024 TorsionAlertsLabel, TorsionAlerts(RotBondIndices TorsionIndices
2025 TorsionAngle Energy EnergyLowerBoundCI EnergyUpperBoundCI
2026 HierarchyClass HierarchySubClass TorsionRule
2027 EnergyMethod AngleNotObserved)
2028
2029 For individual SD files filtered by torsion rules:
2030
2031 TorsionRuleLabel, TorsionRule (HierarchyClass HierarchySubClass
2032 EnergyMethod TorsionRule)
2033 TorsionRuleMaxSingleEnergyAlertsCountLabel,
2034 TorsionRuleMaxSingleEnergyAlertsCount,
2035 TorsionRuleAnglesNotObservedCountLabel,
2036 TorsionRuleAnglesNotObservedCount,
2037 TorsionRuleAlertsLabel, TorsionRuleAlerts (RotBondIndices
2038 TorsionIndices TorsionAngle Energy EnergyLowerBoundCI
2039 EnergyUpperBoundCI EnergyMethod AngleObserved)
2040
2041 --outfileParams <Name,Value,...> [default: auto]
2042 A comma delimited list of parameter name and value pairs for writing
2043 molecules to files. The supported parameter names for different file
2044 formats, along with their default values, are shown below:
2045
2046 SD: kekulize,yes,forceV3000,no
2047
2048 --overwrite
2049 Overwrite existing files.
2050 --precision <number> [default: 4]
2051 Floating point precision for writing torsion strain energy values.
2052 -r, --rotBondsSMARTSMode <NonStrict, SemiStrict,...> [default: SemiStrict]
2053 SMARTS pattern to use for identifying rotatable bonds in a molecule
2054 for matching against torsion rules in the torsion library. Possible values:
2055 NonStrict, SemiStrict, Strict or Specify. The rotatable bond SMARTS matches
2056 are filtered to ensure that each atom in the rotatable bond is attached to
2057 at least two heavy atoms.
2058
2059 The following SMARTS patterns are used to identify rotatable bonds for
2060 different modes:
2061
2062 NonStrict: [!$(*#*)&!D1]-&!@[!$(*#*)&!D1]
2063
2064 SemiStrict:
2065 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
2066 &!$(C([CH3])([CH3])[CH3])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
2067 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
2068
2069 Strict:
2070 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
2071 &!$(C([CH3])([CH3])[CH3])&!$([CD3](=[N,O,S])-!@[#7,O,S!D1])
2072 &!$([#7,O,S!D1]-!@[CD3]=[N,O,S])&!$([CD3](=[N+])-!@[#7!D1])
2073 &!$([#7!D1]-!@[CD3]=[N+])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
2074 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
2075
2076 The 'NonStrict' and 'Strict' SMARTS patterns are available in RDKit. The
2077 'NonStrict' SMARTS pattern corresponds to original Daylight SMARTS
2078 specification for rotatable bonds. The 'SemiStrict' SMARTS pattern is
2079 derived from 'Strict' SMARTS pattern for its usage in this script.
2080
2081 You may use any arbitrary SMARTS pattern to identify rotatable bonds by
2082 choosing 'Specify' value for '-r, --rotBondsSMARTSMode' option and providing its
2083 value via '--rotBondsSMARTSPattern' option.
2084 --rotBondsSMARTSPattern <SMARTS>
2085 SMARTS pattern for identifying rotatable bonds. This option is only valid
2086 for 'Specify' value of '-r, --rotBondsSMARTSMode' option.
2087 -t, --torsionEnergyLibraryFile <FileName or auto> [default: auto]
2088 Specify a XML file name containing data for torsion starin energy library
2089 hierarchy or use default file, TorsionEnergyLibrary.xml, available in
2090 MAYACHEMTOOLS/lib/Python/TorsionAlerts directory.
2091
2092 The format of data in local XML file must match format of the data in Torsion
2093 Library [ Ref 153 ] file available in MAYACHEMTOOLS directory.
2094 -w, --workingdir <dir>
2095 Location of working directory which defaults to the current directory.
2096
2097 Examples:
2098 To filter molecules containing rotatable bonds with total strain energy value
2099 of >= 6.0 (TEUs) based on torsion rules in the torsion energy library and write
2100 write out SD files containing remaining and filtered molecules, type:
2101
2102 % RDKitFilterTorsionStrainEnergyAlerts.py -i Sample3D.sdf
2103 -o Sample3DOut.sdf
2104
2105 To filter molecules containing any rotatable bonds with strain energy value of
2106 >= 1.8 (TEUs) based on torsion rules in the torsion energy library and write out
2107 SD files containing remaining and filtered molecules, and individual SD files for
2108 torsion rules triggering alerts along with appropriate torsion information for
2109 red alerts, type:
2110
2111 % RDKitFilterTorsionStrainEnergyAlerts.py -a MaxSingleEnergy
2112 -i Sample3D.sdf -o Sample3DOut.sdf
2113
2114 To filter molecules containing rotatable bonds with total strain energy value
2115 of >= 6.0 (TEUs) or any single strain energy value of >= 1.8 (TEUs) and write out
2116 SD files containing remaining and filtered molecules, type:
2117
2118 % RDKitFilterTorsionStrainEnergyAlerts.py -a TotalOrMaxSingleEnergy
2119 -i Sample3D.sdf -o Sample3DOut.sdf
2120
2121 To filter molecules containing rotatable bonds with specific cutoff values for
2122 total or single torsion strain energy and write out SD files containing
2123 remaining and filtered molecules, type:
2124
2125 % RDKitFilterTorsionStrainEnergyAlerts.py -a TotalOrMaxSingleEnergy
2126 -i Sample3D.sdf -o Sample3DOut.sdf --alertsTotalEnergyCutoff 6.0
2127 --alertsMaxSingleEnergyCutoff 1.8
2128
2129 To run the first example for filtering molecules and writing out torsion
2130 information for all alert types to SD files, type:
2131
2132 % RDKitFilterTorsionStrainEnergyAlerts.py -i Sample3D.sdf
2133 -o Sample3DOut.sdf --outfileAlertsMode All
2134
2135 To run the first example for filtering molecules in multiprocessing mode on
2136 all available CPUs without loading all data into memory and write out SD files,
2137 type:
2138
2139 % RDKitFilterTorsionStrainEnergyAlerts.py --mp yes -i Sample3D.sdf
2140 -o Sample3DOut.sdf
2141
2142 To run the first example for filtering molecules in multiprocessing mode on
2143 all available CPUs by loading all data into memory and write out a SD files,
2144 type:
2145
2146 % RDKitFilterTorsionStrainEnergyAlerts.py --mp yes --mpParams
2147 "inputDataMode, InMemory" -i Sample3D.sdf -o Sample3DOut.sdf
2148
2149 To run the first example for filtering molecules in multiprocessing mode on
2150 specific number of CPUs and chunksize without loading all data into memory
2151 and write out SD files, type:
2152
2153 % RDKitFilterTorsionStrainEnergyAlerts.py --mp yes --mpParams
2154 "inputDataMode,lazy,numProcesses,4,chunkSize,8" -i Sample3D.sdf
2155 -o Sample3DOut.sdf
2156
2157 To list information about default torsion library file without performing any
2158 filtering, type:
2159
2160 % RDKitFilterTorsionStrainEnergyAlerts.py -l
2161
2162 To list information about a local torsion library XML file without performing
2163 any, filtering, type:
2164
2165 % RDKitFilterTorsionStrainEnergyAlerts.py --torsionEnergyLibraryFile
2166 TorsionStrainEnergyLibrary.xml -l
2167
2168 Author:
2169 Manish Sud (msud@san.rr.com)
2170
2171 Collaborator:
2172 Pat Walters
2173
2174 See also:
2175 RDKitFilterChEMBLAlerts.py, RDKitFilterPAINS.py, RDKitFilterTorsionLibraryAlerts.py,
2176 RDKitConvertFileFormat.py, RDKitSearchSMARTS.py
2177
2178 Copyright:
2179 Copyright (C) 2026 Manish Sud. All rights reserved.
2180
2181 This script uses the torsion strain energy library developed by Gu, S.;
2182 Smith, M. S.; Yang, Y.; Irwin, J. J.; Shoichet, B. K. [ Ref 153 ].
2183
2184 The torsion strain enegy library is based on the Torsion Library jointly
2185 developed by the University of Hamburg, Center for Bioinformatics,
2186 Hamburg, Germany and F. Hoffmann-La-Roche Ltd., Basel, Switzerland.
2187
2188 The functionality available in this script is implemented using RDKit, an
2189 open source toolkit for cheminformatics developed by Greg Landrum.
2190
2191 This file is part of MayaChemTools.
2192
2193 MayaChemTools is free software; you can redistribute it and/or modify it under
2194 the terms of the GNU Lesser General Public License as published by the Free
2195 Software Foundation; either version 3 of the License, or (at your option) any
2196 later version.
2197
2198 """
2199
2200 if __name__ == "__main__":
2201 main()
