1 #!/usr/bin/perl -w 2 # 3 # $RCSfile: ElementalAnalysisSDFiles.pl,v $ 4 # $Date: 2010/01/03 00:59:50 $ 5 # $Revision: 1.13 $ 6 # 7 # Author: Manish Sud <msud@san.rr.com> 8 # 9 # Copyright (C) 2004-2010 Manish Sud. All rights reserved. 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 use strict; 30 use FindBin; use lib "$FindBin::Bin/../lib"; 31 use Getopt::Long; 32 use File::Basename; 33 use Text::ParseWords; 34 use Benchmark; 35 use FileUtil; 36 use SDFileUtil; 37 use TextUtil; 38 use MolecularFormula; 39 40 my($ScriptName, %Options, $StartTime, $EndTime, $TotalTime); 41 42 # Autoflush STDOUT 43 $| = 1; 44 45 # Starting message... 46 $ScriptName = basename($0); 47 print "\n$ScriptName: Starting...\n\n"; 48 $StartTime = new Benchmark; 49 50 # Get the options and setup script... 51 SetupScriptUsage(); 52 if ($Options{help} || @ARGV < 1) { 53 die GetUsageFromPod("$FindBin::Bin/$ScriptName"); 54 } 55 56 my(@SDFilesList); 57 @SDFilesList = ExpandFileNames(\@ARGV, "sdf sd"); 58 59 my($SpecifiedFormulaFieldName, $CheckFormula, $Precision, $DetailLevel, @SpecifiedCalculations, @SpecifiedValueFieldNames); 60 ProcessOptions(); 61 62 print "Checking input SD file(s)...\n"; 63 my(@SDFilesOkay, @SDFilesOutFile, @SDFilesFormulaFieldName, @SDFilesValueFieldNamesMap); 64 RetrieveSDFilesInfo(); 65 66 # Generate output files... 67 my($Index, $SDFile); 68 if (@SDFilesList > 1) { 69 print "Processing SD file(s)...\n"; 70 } 71 for $Index (0 .. $#SDFilesList) { 72 if ($SDFilesOkay[$Index]) { 73 $SDFile = $SDFilesList[$Index]; 74 if (@SDFilesList > 1) { 75 print "\nProcessing file $SDFile...\n"; 76 } 77 else { 78 print "Processing file $SDFile...\n" 79 } 80 PerformElementalAnalysis($Index); 81 } 82 } 83 84 print "$ScriptName:Done...\n\n"; 85 86 $EndTime = new Benchmark; 87 $TotalTime = timediff ($EndTime, $StartTime); 88 print "Total time: ", timestr($TotalTime), "\n"; 89 90 ############################################################################### 91 92 # Perform analysis... 93 sub PerformElementalAnalysis { 94 my($Index) = @_; 95 my($SDFile, $NewSDFile, $KeyDataFieldName); 96 97 $SDFile = $SDFilesList[$Index]; 98 $NewSDFile = $SDFilesOutFile[$Index]; 99 100 print "Generating new SD file $NewSDFile...\n"; 101 open NEWSDFILE, ">$NewSDFile" or die "Error: Couldn't open $NewSDFile: $! \n"; 102 open SDFILE, "$SDFile" or die "Error: Can't open $SDFile: $! \n"; 103 104 # Go over all compound records and store 'em using key value as hash... 105 my($CmpdCount, $FormulaFieldName, $FormulaFieldValue, $CmpdString, $Value, $CalculationType, $CalculatedValue, $ErrorMsg, $Status, $ElementsRef, $ElementCompositionRef, @CalculatedValues, @CmpdLines, %DataFieldValuesMap); 106 $CmpdCount = 0; 107 $FormulaFieldName = $SDFilesFormulaFieldName[$Index]; 108 109 COMPOUND: while ($CmpdString = ReadCmpdString(\*SDFILE)) { 110 $CmpdCount++; 111 @CmpdLines = split "\n", $CmpdString; 112 %DataFieldValuesMap = GetCmpdDataHeaderLabelsAndValues(\@CmpdLines); 113 114 @CalculatedValues = (); 115 for $Value (@SpecifiedCalculations) { 116 push @CalculatedValues, ''; 117 } 118 119 if (!exists $DataFieldValuesMap{$FormulaFieldName}) { 120 $ErrorMsg = "Ignoring compound record $CmpdCount: Formula field $FormulaFieldName not found"; 121 PrintErrorMsg($CmpdString, $ErrorMsg); 122 WriteNewCompoundRecord($Index, \*NEWSDFILE, \@CmpdLines, \@CalculatedValues); 123 next COMPOUND; 124 } 125 126 # Make sure it's a valid molecular formula... 127 $FormulaFieldValue = $DataFieldValuesMap{$FormulaFieldName}; 128 if ($CheckFormula) { 129 ($Status, $ErrorMsg) = MolecularFormula::IsMolecularFormula($FormulaFieldValue); 130 if (!$Status) { 131 $ErrorMsg = "Ignoring compound record $CmpdCount: Formula field value $FormulaFieldValue is not valid: $ErrorMsg"; 132 PrintErrorMsg($CmpdString, $ErrorMsg); 133 WriteNewCompoundRecord($Index, \*NEWSDFILE, \@CmpdLines, \@CalculatedValues); 134 next COMPOUND; 135 } 136 } 137 # Calculate appropriate values and write 'em out... 138 @CalculatedValues = (); 139 for $CalculationType (@SpecifiedCalculations) { 140 if ($CalculationType =~ /^ElementalAnalysis$/i) { 141 ($ElementsRef, $ElementCompositionRef) = MolecularFormula::CalculateElementalComposition($FormulaFieldValue); 142 $CalculatedValue = (defined($ElementsRef) && defined($ElementCompositionRef)) ? MolecularFormula::FormatCompositionInfomation($ElementsRef, $ElementCompositionRef, $Precision) : ''; 143 } 144 elsif ($CalculationType =~ /^MolecularWeight$/i) { 145 $CalculatedValue = MolecularFormula::CalculateMolecularWeight($FormulaFieldValue); 146 $CalculatedValue = (defined($CalculatedValue) && length($CalculatedValue)) ? (sprintf("%.${Precision}f", $CalculatedValue)) : ""; 147 } 148 elsif ($CalculationType =~ /^ExactMass$/i) { 149 $CalculatedValue = MolecularFormula::CalculateExactMass($FormulaFieldValue); 150 $CalculatedValue = (defined($CalculatedValue) && length($CalculatedValue)) ? (sprintf("%.${Precision}f", $CalculatedValue)) : ""; 151 } 152 else { 153 $CalculatedValue = ''; 154 } 155 push @CalculatedValues, $CalculatedValue; 156 } 157 WriteNewCompoundRecord($Index, \*NEWSDFILE, \@CmpdLines, \@CalculatedValues); 158 } 159 close NEWSDFILE; 160 close SDFILE; 161 } 162 163 # Write out compound record with calculated values... 164 sub WriteNewCompoundRecord { 165 my($Index, $SDFileRef, $CmpdLinesRef, $CalculatedValuesRef) = @_; 166 167 # Write out compound lines except the last line which contains $$$$... 168 my($LineIndex); 169 for $LineIndex (0 .. ($#{$CmpdLinesRef} - 1)) { 170 print $SDFileRef "$CmpdLinesRef->[$LineIndex]\n"; 171 } 172 173 # Write out calculated values... 174 my($CalcIndex, $FieldName, $FieldValue); 175 for $CalcIndex (0 .. $#SpecifiedCalculations) { 176 $FieldName = $SDFilesValueFieldNamesMap[$Index]{$SpecifiedCalculations[$CalcIndex]}; 177 $FieldValue = $CalculatedValuesRef->[$CalcIndex]; 178 print $SDFileRef "> <$FieldName>\n$FieldValue\n\n"; 179 } 180 print $SDFileRef "\$\$\$\$\n"; 181 } 182 183 # Print out error message... 184 sub PrintErrorMsg { 185 my($CmpdString, $ErrorMsg) = @_; 186 187 if ($DetailLevel >= 2 ) { 188 print "$ErrorMsg:\n$CmpdString\n"; 189 } 190 elsif ($DetailLevel >= 1) { 191 print "$ErrorMsg\n"; 192 } 193 } 194 195 # Process option values... 196 sub ProcessOptions { 197 $DetailLevel = $Options{detail}; 198 $CheckFormula = $Options{fast} ? 0 : 1; 199 $Precision = $Options{precision}; 200 201 $SpecifiedFormulaFieldName = ""; 202 if (defined $Options{formulafield}) { 203 $SpecifiedFormulaFieldName = $Options{formulafield}; 204 } 205 206 # Setup what to calculate... 207 @SpecifiedCalculations = (); 208 if ($Options{mode} =~ /^All$/i) { 209 @SpecifiedCalculations = qw(ElementalAnalysis MolecularWeight ExactMass); 210 } 211 else { 212 my($Mode, $ModeValue, @SpecifiedModeValues); 213 $Mode = $Options{mode}; 214 $Mode =~ s/ //g; 215 @SpecifiedModeValues = split /\,/, $Mode; 216 for $ModeValue (@SpecifiedModeValues) { 217 if ($ModeValue !~ /^(ElementalAnalysis|MolecularWeight|ExactMass)$/i) { 218 if ($ModeValue =~ /^All$/i) { 219 die "Error: All value for option \"-m --mode\" is not allowed with other valid values.\n"; 220 } 221 else { 222 die "Error: The value specified, $ModeValue, for option \"-m --mode\" is not valid. Allowed values: ElementalAnalysis, MolecularWeight, or ExactMass\n"; 223 } 224 } 225 push @SpecifiedCalculations, $ModeValue; 226 } 227 } 228 229 @SpecifiedValueFieldNames = (); 230 if ($Options{valuefieldnames}) { 231 my($Value, $Label, @ValueLabels); 232 @ValueLabels = split /\,/, $Options{valuefieldnames}; 233 if (@ValueLabels % 2) { 234 die "Error: The value specified, $Options{valuefieldnames}, for option \"-v --valuefieldnames\" is not valid: It must contain even number of comma delimited values\n"; 235 } 236 for ($Index = 0; $Index < @ValueLabels; $Index +=2) { 237 $Value = $ValueLabels[$Index]; 238 $Value =~ s/ //g; 239 $Label = $ValueLabels[$Index + 1]; 240 if ($Value !~ /^(ElementalAnalysis|MolecularWeight|ExactMass)$/i) { 241 die "Error: The value specified, $Value, using option \"-v --valuefieldnames\" is not valid. Allowed values: ElementalAnalysis, MolecularWeight, or ExactMass\n"; 242 } 243 push @SpecifiedValueFieldNames, ($Value, $Label); 244 } 245 } 246 } 247 248 # Retrieve information about input SD files... 249 sub RetrieveSDFilesInfo { 250 my($Index, $SDFile, $FileDir, $FileName, $FileExt, $OutFileRoot, $OutFile, $FormulaFieldName, $Value, $FieldName, $NewFieldName, $Count); 251 252 my(%NewValueFieldNames) = (ElementalAnalysis => 'ElementalAnalysis', MolecularWeight => 'MolecularWeight', ExactMass => 'ExactMass'); 253 if (@SpecifiedValueFieldNames) { 254 for ($Index = 0; $Index < @SpecifiedValueFieldNames; $Index +=2) { 255 $Value = $SpecifiedValueFieldNames[$Index]; 256 $FieldName = $SpecifiedValueFieldNames[$Index + 1]; 257 if (exists $NewValueFieldNames{$Value}) { 258 $NewValueFieldNames{$Value} = $FieldName; 259 } 260 } 261 } 262 263 @SDFilesOkay = (); 264 @SDFilesOutFile = (); 265 @SDFilesFormulaFieldName = (); 266 @SDFilesValueFieldNamesMap = (); 267 268 FILELIST: for $Index (0 .. $#SDFilesList) { 269 $SDFile = $SDFilesList[$Index]; 270 $SDFilesOkay[$Index] = 0; 271 $SDFilesOutFile[$Index] = ''; 272 $SDFilesFormulaFieldName[$Index] = ''; 273 %{$SDFilesValueFieldNamesMap[$Index]} = (); 274 275 if (!(-e $SDFile)) { 276 warn "Warning: Ignoring file $SDFile: It doesn't exist\n"; 277 next FILELIST; 278 } 279 if (!CheckFileType($SDFile, "sd sdf")) { 280 warn "Warning: Ignoring file $SDFile: It's not a SD file\n"; 281 next FILELIST; 282 } 283 $FileDir = ""; $FileName = ""; $FileExt = ""; 284 ($FileDir, $FileName, $FileExt) = ParseFileName($SDFile); 285 if ($Options{root} && (@SDFilesList == 1)) { 286 my ($RootFileDir, $RootFileName, $RootFileExt) = ParseFileName($Options{root}); 287 if ($RootFileName && $RootFileExt) { 288 $FileName = $RootFileName; 289 } 290 else { 291 $FileName = $Options{root}; 292 } 293 $OutFileRoot = $FileName; 294 } 295 else { 296 $OutFileRoot = $FileName . "ElementalAnalysis"; 297 } 298 299 $OutFile = $OutFileRoot . ".$FileExt"; 300 if (lc($OutFile) eq lc($SDFile)) { 301 warn "Warning: Ignoring file $SDFile:Output file name, $OutFile, is same as input SD file name, $SDFile\n"; 302 next FILELIST; 303 } 304 if (!$Options{overwrite}) { 305 if (-e $OutFile) { 306 warn "Warning: Ignoring file $SDFile: The file $OutFile already exists\n"; 307 next FILELIST; 308 } 309 } 310 # Get data field names and values... 311 my($CmpdString, $FieldName, @CmpdLines, @DataFieldNames, %DataFieldNamesMap); 312 @DataFieldNames = (); 313 if (!open(SDFILE, "$SDFile")) { 314 warn "Warning: Ignoring file $SDFile: Couldn't open it: $! \n"; 315 next FILELIST; 316 } 317 $CmpdString = ReadCmpdString(\*SDFILE); 318 close SDFILE; 319 320 @CmpdLines = split "\n", $CmpdString; 321 @DataFieldNames = GetCmpdDataHeaderLabels(\@CmpdLines); 322 %DataFieldNamesMap = GetCmpdDataHeaderLabelsAndValues(\@CmpdLines); 323 324 # Setup formula field name... 325 $FormulaFieldName = ''; 326 if ($SpecifiedFormulaFieldName) { 327 $FormulaFieldName = $SpecifiedFormulaFieldName; 328 } 329 else { 330 FIELDNAME: for $FieldName (@DataFieldNames) { 331 if ($FieldName =~ /Formula/i) { 332 $FormulaFieldName = $FieldName; 333 last FIELDNAME; 334 } 335 } 336 if (!$FormulaFieldName) { 337 warn "Warning: Ignoring file $SDFile: Data field label containing the word Formula doesn't exist\n"; 338 next FILELIST; 339 } 340 } 341 $SDFilesOkay[$Index] = 1; 342 $SDFilesOutFile[$Index] = $OutFile; 343 $SDFilesFormulaFieldName[$Index] = $FormulaFieldName; 344 345 # Setup value data field names for calculated values... 346 for $Value (keys %NewValueFieldNames) { 347 $FieldName = $NewValueFieldNames{$Value}; 348 349 # Make sure it doesn't already exists... 350 $Count = 1; 351 $NewFieldName = $FieldName; 352 while (exists $DataFieldNamesMap{$NewFieldName}) { 353 $Count++; 354 $NewFieldName = $FieldName . $Count; 355 } 356 $SDFilesValueFieldNamesMap[$Index]{$Value} = $NewFieldName; 357 } 358 } 359 } 360 361 # Setup script usage and retrieve command line arguments specified using various options... 362 sub SetupScriptUsage { 363 364 # Retrieve all the options... 365 %Options = (); 366 $Options{detail} = 1; 367 $Options{mode} = "All"; 368 $Options{precision} = 2; 369 370 if (!GetOptions(\%Options, "detail|d=i", "fast", "formulafield|f=s", "mode|m=s", "help|h", "overwrite|o", "precision|p=i", "root|r=s", "valuefieldnames|v=s", "workingdir|w=s")) { 371 die "\nTo get a list of valid options and their values, use \"$ScriptName -h\" or\n\"perl -S $ScriptName -h\" command and try again...\n"; 372 } 373 if ($Options{workingdir}) { 374 if (! -d $Options{workingdir}) { 375 die "Error: The value specified, $Options{workingdir}, for option \"-w --workingdir\" is not a directory name.\n"; 376 } 377 chdir $Options{workingdir} or die "Error: Couldn't chdir $Options{workingdir}: $! \n"; 378 } 379 if (!IsPositiveInteger($Options{detail})) { 380 die "Error: The value specified, $Options{detail}, for option \"-d --detail\" is not valid. Allowed values: > 0\n"; 381 } 382 if (!IsPositiveInteger($Options{precision})) { 383 die "Error: The value specified, $Options{precision}, for option \"-p --precision\" is not valid. Allowed values: > 0 \n"; 384 } 385 } 386
