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