1 #!/bin/env python 2 # 3 # File: RDKitCalculateRMSD.py 4 # Author: Manish Sud <msud@san.rr.com> 5 # 6 # Copyright (C) 2026 Manish Sud. All rights reserved. 7 # 8 # The functionality available in this script is implemented using RDKit, an 9 # open source toolkit for cheminformatics developed by Greg Landrum. 10 # 11 # This file is part of MayaChemTools. 12 # 13 # MayaChemTools is free software; you can redistribute it and/or modify it under 14 # the terms of the GNU Lesser General Public License as published by the Free 15 # Software Foundation; either version 3 of the License, or (at your option) any 16 # later version. 17 # 18 # MayaChemTools is distributed in the hope that it will be useful, but without 19 # any warranty; without even the implied warranty of merchantability of fitness 20 # for a particular purpose. See the GNU Lesser General Public License for more 21 # details. 22 # 23 # You should have received a copy of the GNU Lesser General Public License 24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or 25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330, 26 # Boston, MA, 02111-1307, USA. 27 # 28 29 from __future__ import print_function 30 31 import os 32 import sys 33 import time 34 import re 35 36 # RDKit imports... 37 try: 38 from rdkit import rdBase 39 from rdkit.Chem import AllChem, rdMolAlign 40 except ImportError as ErrMsg: 41 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg) 42 sys.stderr.write("Check/update your RDKit environment and try again.\n\n") 43 sys.exit(1) 44 45 # MayaChemTools imports... 46 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python")) 47 try: 48 from docopt import docopt 49 import MiscUtil 50 import RDKitUtil 51 except ImportError as ErrMsg: 52 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg) 53 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n") 54 sys.exit(1) 55 56 ScriptName = os.path.basename(sys.argv[0]) 57 Options = {} 58 OptionsInfo = {} 59 60 61 def main(): 62 """Start execution of the script.""" 63 64 MiscUtil.PrintInfo( 65 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 66 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()) 67 ) 68 69 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 70 71 # Retrieve command line arguments and options... 72 RetrieveOptions() 73 74 # Process and validate command line arguments and options... 75 ProcessOptions() 76 77 # Perform actions required by the script... 78 CalculateRMSD() 79 80 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 81 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 82 83 84 def CalculateRMSD(): 85 """Calculate RMSD values.""" 86 87 if not re.match("^(OneToOne|AllToAll|FirstToAll)$", OptionsInfo["Mode"], re.I): 88 MiscUtil.PrintError("RMSD couldn't be calculated: Specified mode, %s, is not supported" % OptionsInfo["Mode"]) 89 90 RefFile = OptionsInfo["RefFile"] 91 ProbeFile = OptionsInfo["ProbeFile"] 92 93 Outfile = OptionsInfo["Outfile"] 94 OutDelim = OptionsInfo["OutDelim"] 95 96 # Read reference and probe molecules... 97 OptionsInfo["InfileParams"]["AllowEmptyMols"] = False 98 99 MiscUtil.PrintInfo("\nProcessing file %s..." % (RefFile)) 100 ValidRefMols, RefMolCount, ValidRefMolCount = RDKitUtil.ReadAndValidateMolecules( 101 RefFile, **OptionsInfo["InfileParams"] 102 ) 103 104 MiscUtil.PrintInfo("Processing file %s..." % (ProbeFile)) 105 ValidProbeMols, ProbeMolCount, ValidProbeMolCount = RDKitUtil.ReadAndValidateMolecules( 106 ProbeFile, **OptionsInfo["InfileParams"] 107 ) 108 109 # Set up output file... 110 MiscUtil.PrintInfo("\nGenerating file %s...\n" % Outfile) 111 OutFH = open(Outfile, "w") 112 if OutFH is None: 113 MiscUtil.PrintError("Couldn't open output file: %s.\n" % (Outfile)) 114 115 Line = "RefMolID%sProbeMolID%sRMSD\n" % (OutDelim, OutDelim) 116 OutFH.write(Line) 117 118 if re.match("^OneToOne$", OptionsInfo["Mode"], re.I): 119 CalculateOneToOneRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim) 120 elif re.match("^AllToAll$", OptionsInfo["Mode"], re.I): 121 CalculateAllToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim) 122 elif re.match("^FirstToAll$", OptionsInfo["Mode"], re.I): 123 CalculateFirstToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim) 124 else: 125 MiscUtil.PrintError("RMSD couldn't be calculated: Specified mode, %s, is not supported" % OptionsInfo["Mode"]) 126 127 OutFH.close() 128 129 MiscUtil.PrintInfo("\nTotal number of molecules: Reference - %d; Probe - %d" % (RefMolCount, ProbeMolCount)) 130 MiscUtil.PrintInfo("Number of valid molecules: Reference - %d; Probe - %d" % (ValidRefMolCount, ValidProbeMolCount)) 131 MiscUtil.PrintInfo( 132 "Number of ignored molecules: Reference - %d; Probe - %d" 133 % ((RefMolCount - ValidRefMolCount), (ProbeMolCount - ValidProbeMolCount)) 134 ) 135 136 137 def CalculateOneToOneRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim): 138 """Calculate pairwise RMSD values.""" 139 140 ValidRefMolCount = len(ValidRefMols) 141 ValidProbeMolCount = len(ValidProbeMols) 142 143 MolCount = ValidRefMolCount 144 if ValidRefMolCount > ValidProbeMolCount: 145 MolCount = ValidProbeMolCount 146 147 if ValidRefMolCount != ValidProbeMolCount: 148 MiscUtil.PrintWarning( 149 "Number of valid reference molecules, %d, is not equal to number of valid probe molecules, %d .\n" 150 % (ValidRefMolCount, ValidProbeMolCount) 151 ) 152 MiscUtil.PrintWarning("Pairwise RMSD will be calculated only for first %s molecules.\n" % (MolCount)) 153 154 # Process molecules... 155 for MolIndex in range(0, MolCount): 156 RefMol = ValidRefMols[MolIndex] 157 ProbeMol = ValidProbeMols[MolIndex] 158 159 RefMolName = RDKitUtil.GetMolName(RefMol, (MolIndex + 1)) 160 ProbeMolName = RDKitUtil.GetMolName(ProbeMol, (MolIndex + 1)) 161 162 RMSD = CalculateRMSDValue(RefMol, ProbeMol) 163 164 Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD) 165 OutFH.write(Line) 166 167 168 def CalculateAllToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim): 169 """Calculate RMSD values between all pairs of molecules.""" 170 171 # Process molecules... 172 RefMolCount = 0 173 for RefMol in ValidRefMols: 174 RefMolCount += 1 175 RefMolName = RDKitUtil.GetMolName(RefMol, RefMolCount) 176 177 ProbeMolCount = 0 178 for ProbeMol in ValidProbeMols: 179 ProbeMolCount += 1 180 ProbeMolName = RDKitUtil.GetMolName(ProbeMol, ProbeMolCount) 181 182 RMSD = CalculateRMSDValue(RefMol, ProbeMol) 183 184 Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD) 185 OutFH.write(Line) 186 187 188 def CalculateFirstToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim): 189 """Calculate RMSD values between first reference molecues and all probe molecules.""" 190 191 # Process molecules... 192 RefMol = ValidRefMols[0] 193 RefMolCount = 1 194 RefMolName = RDKitUtil.GetMolName(RefMol, RefMolCount) 195 196 ProbeMolCount = 0 197 for ProbeMol in ValidProbeMols: 198 ProbeMolCount += 1 199 ProbeMolName = RDKitUtil.GetMolName(ProbeMol, ProbeMolCount) 200 201 RMSD = CalculateRMSDValue(RefMol, ProbeMol) 202 203 Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD) 204 OutFH.write(Line) 205 206 207 def CalculateRMSDValue(RefMol, ProbeMol): 208 """Calculate RMSD value for a pair of molecules and return it as a string.""" 209 210 try: 211 if OptionsInfo["UseBestRMSD"]: 212 RMSD = AllChem.GetBestRMS(ProbeMol, RefMol) 213 else: 214 RMSD = rdMolAlign.AlignMol(ProbeMol, RefMol, maxIters=OptionsInfo["MaxIters"]) 215 RMSD = "%.2f" % RMSD 216 except (RuntimeError, ValueError): 217 RMSD = "None" 218 219 return RMSD 220 221 222 def ProcessOptions(): 223 """Process and validate command line arguments and options.""" 224 225 MiscUtil.PrintInfo("Processing options...") 226 227 # Validate options... 228 ValidateOptions() 229 230 OptionsInfo["CalcRMSD"] = Options["--calcRMSD"] 231 OptionsInfo["UseBestRMSD"] = False 232 if re.match("^BestRMSD$", OptionsInfo["CalcRMSD"], re.I): 233 OptionsInfo["UseBestRMSD"] = True 234 235 OptionsInfo["MaxIters"] = int(Options["--maxIters"]) 236 237 OptionsInfo["Mode"] = Options["--mode"] 238 239 OptionsInfo["RefFile"] = Options["--reffile"] 240 OptionsInfo["ProbeFile"] = Options["--probefile"] 241 242 # No need for any RDKit specific --outfileParams.... 243 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"]) 244 245 OptionsInfo["Outfile"] = Options["--outfile"] 246 247 OptionsInfo["Overwrite"] = Options["--overwrite"] 248 249 OptionsInfo["OutDelim"] = " " 250 if MiscUtil.CheckFileExt(OptionsInfo["Outfile"], "csv"): 251 OptionsInfo["OutDelim"] = "," 252 elif MiscUtil.CheckFileExt(OptionsInfo["Outfile"], "tsv txt"): 253 OptionsInfo["OutDelim"] = "\t" 254 else: 255 MiscUtil.PrintError( 256 'The file name specified , %s, for option "--outfile" is not valid. Supported file formats: csv tsv txt\n' 257 % (OptionsInfo["Outfile"]) 258 ) 259 260 261 def RetrieveOptions(): 262 """Retrieve command line arguments and options.""" 263 264 # Get options... 265 global Options 266 Options = docopt(_docoptUsage_) 267 268 # Set current working directory to the specified directory... 269 WorkingDir = Options["--workingdir"] 270 if WorkingDir: 271 os.chdir(WorkingDir) 272 273 # Handle examples option... 274 if "--examples" in Options and Options["--examples"]: 275 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 276 sys.exit(0) 277 278 279 def ValidateOptions(): 280 """Validate option values.""" 281 282 MiscUtil.ValidateOptionTextValue("--calcRMSD", Options["--calcRMSD"], "RMSD BestRMSD") 283 284 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0}) 285 MiscUtil.ValidateOptionTextValue("--mode", Options["--mode"], "OneToOne AllToAll FirstToAll") 286 287 MiscUtil.ValidateOptionFilePath("-r, --reffile", Options["--reffile"]) 288 MiscUtil.ValidateOptionFileExt("-r, --reffile", Options["--reffile"], "sdf sd mol") 289 290 MiscUtil.ValidateOptionFilePath("-p, --probefile", Options["--probefile"]) 291 MiscUtil.ValidateOptionFileExt("-p, --probefile", Options["--probefile"], "sdf sd mol") 292 293 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "csv tsv txt") 294 MiscUtil.ValidateOptionsOutputFileOverwrite( 295 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 296 ) 297 MiscUtil.ValidateOptionsDistinctFileNames( 298 "-r, --reffile", Options["--reffile"], "-o, --outfile", Options["--outfile"] 299 ) 300 MiscUtil.ValidateOptionsDistinctFileNames( 301 "-p, --probefile", Options["--probefile"], "-o, --outfile", Options["--outfile"] 302 ) 303 304 305 # Setup a usage string for docopt... 306 _docoptUsage_ = """ 307 RDKitCalculateRMSD.py - Calculate RMSD between molecules 308 309 Usage: 310 RDKitCalculateRMSD.py [--calcRMSD <RMSD, BestRMSD>] [--infileParams <Name,Value,...>] 311 [--maxIters <number>] [--mode <OneToOne, AllToAll, FirstToAll>] 312 [--overwrite] [-w <dir>] -r <reffile> -p <probefile> -o <outfile> 313 RDKitCalculateRMSD.py -h | --help | -e | --examples 314 315 Description: 316 Calculate Root Mean Square Distance (RMSD) between a set of similar molecules in 317 reference and probe input files. The RDKit function fails to calculate RMSD values for 318 dissimilar molecules. Consequently, a text string 'None' is written out as a RMSD value 319 for dissimilar molecule pairs. 320 321 The supported input file formats are: Mol (.mol), SD (.sdf, .sd) 322 323 The supported output file formats are: CSV (.csv), TSV (.tsv, .txt) 324 325 Options: 326 -c, --calcRMSD <RMSD, BestRMSD> [default: RMSD] 327 Methodology for calculating RMSD values. Possible values: RMSD, BestRMSD. 328 During BestRMSMode mode, the RDKit 'function AllChem.GetBestRMS' is used to 329 align and calculate RMSD. This function calculates optimal RMSD for aligning two 330 molecules, taking symmetry into account. Otherwise, the RMSD value is calculated 331 using 'AllChem.AlignMol function' without changing the atom order. A word to the 332 wise from RDKit documentation: The AllChem.GetBestRMS function will attempt to 333 align all permutations of matching atom orders in both molecules, for some molecules 334 it will lead to 'combinatorial explosion'. 335 --infileParams <Name,Value,...> [default: auto] 336 A comma delimited list of parameter name and value pairs for reading 337 molecules from files. The supported parameter names for different file 338 formats, along with their default values, are shown below: 339 340 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes 341 342 --maxIters <number> [default: 50] 343 Maximum number of iterations to perform for each molecule pair during minimization 344 of RMSD values. This option is ignored during BestRMSD mode. 345 -m, --mode <OneToOne, AllToAll, FirstToAll> [default: OneToOne] 346 Specify how molecules are handled in reference and probe input files during 347 calculation of RMSD between reference and probe molecules. Possible values: 348 OneToOne, AllToAll and AllToFirst. For OneToOne mode, the number of molecules 349 in reference file must be equal to the number of molecules in probe file. The RMSD 350 is calculated for each pair of molecules in the reference and probe file and written 351 to the output file. For AllToAll mode, the RMSD is calculated for each reference 352 molecule against all probe molecules. For FirstToAll mode, however, the RMSD 353 is only calculated for the first reference molecule against all probe molecules. 354 -e, --examples 355 Print examples. 356 -h, --help 357 Print this help message. 358 -p, --probefile <probefile> 359 Probe input file name. 360 -r, --reffile <reffile> 361 Reference input file name. 362 -o, --outfile <outfile> 363 Output file name for writing out RMSD values. Supported text file extensions: csv or tsv. 364 --overwrite 365 Overwrite existing files. 366 -w, --workingdir <dir> 367 Location of working directory which defaults to the current directory. 368 369 Examples: 370 To calculate RMSD between pair of molecules in reference and probe input 371 3D SD files and write out a CSV file containing calculated RMSD values along with 372 appropriate molecule IDs, type: 373 374 % RDKitCalculateRMSD.py -r Sample3DRef.sdf -p Sample3DProb.sdf 375 -o SampleOut.csv 376 377 To calculate RMSD between all molecules in reference and probe input 378 3D SD files and write out a CSV file containing calculated RMSD values along with 379 appropriate molecule IDs, type: 380 381 % RDKitCalculateRMSD.py -m AllToAll -r Sample3DRef.sdf -p 382 Sample3DProb.sdf -o SampleOut.csv 383 384 To calculate best RMSD between first molecule in reference all probe molecules 385 in 3D SD files and write out a TSV file containing calculated RMSD values along with 386 appropriate molecule IDs, type: 387 388 % RDKitCalculateRMSD.py -m FirstToAll --calcRMSD BestRMSD -r 389 Sample3DRef.sdf -p Sample3DProb.sdf -o SampleOut.tsv 390 391 To calculate RMSD between all molecules in reference and probe molecules input 392 3D SD files without removing hydrogens and write out a TSV file containing 393 calculated RMSD values along with appropriate molecule IDs, type: 394 395 % RDKitCalculateRMSD.py -m AllToAll --infileParams 396 "removeHydrogens,no" -r Sample3DRef.sdf -p Sample3DProb.sdf 397 -o SampleOut.tsv 398 399 Author: 400 Manish Sud(msud@san.rr.com) 401 402 See also: 403 RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py, RDKitConvertFileFormat.py, 404 RDKitGenerateConformers.py, RDKitPerformMinimization.py 405 406 Copyright: 407 Copyright (C) 2026 Manish Sud. All rights reserved. 408 409 The functionality available in this script is implemented using RDKit, an 410 open source toolkit for cheminformatics developed by Greg Landrum. 411 412 This file is part of MayaChemTools. 413 414 MayaChemTools is free software; you can redistribute it and/or modify it under 415 the terms of the GNU Lesser General Public License as published by the Free 416 Software Foundation; either version 3 of the License, or (at your option) any 417 later version. 418 419 """ 420 421 if __name__ == "__main__": 422 main()
