1 #!/usr/bin/perl -w
   2 #
   3 # File: AnalyzeSequenceFilesData.pl
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2024 Manish Sud. All rights reserved.
   7 #
   8 # This file is part of MayaChemTools.
   9 #
  10 # MayaChemTools is free software; you can redistribute it and/or modify it under
  11 # the terms of the GNU Lesser General Public License as published by the Free
  12 # Software Foundation; either version 3 of the License, or (at your option) any
  13 # later version.
  14 #
  15 # MayaChemTools is distributed in the hope that it will be useful, but without
  16 # any warranty; without even the implied warranty of merchantability of fitness
  17 # for a particular purpose.  See the GNU Lesser General Public License for more
  18 # details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public License
  21 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  22 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  23 # Boston, MA, 02111-1307, USA.
  24 #
  25 
  26 use strict;
  27 use FindBin; use lib "$FindBin::Bin/../lib";
  28 use Getopt::Long;
  29 use File::Basename;
  30 use Text::ParseWords;
  31 use Benchmark;
  32 use FileUtil;
  33 use TextUtil;
  34 use SequenceFileUtil;
  35 use AminoAcids;
  36 use NucleicAcids;
  37 
  38 my($ScriptName, %Options, $StartTime, $EndTime, $TotalTime);
  39 
  40 # Autoflush STDOUT
  41 $| = 1;
  42 
  43 # Starting message...
  44 $ScriptName = basename($0);
  45 print "\n$ScriptName: Starting...\n\n";
  46 $StartTime = new Benchmark;
  47 
  48 # Setup script usage message...
  49 SetupScriptUsage();
  50 if ($Options{help} || @ARGV < 1) {
  51   die GetUsageFromPod("$FindBin::Bin/$ScriptName");
  52 }
  53 
  54 # Expand wild card file names...
  55 my(@SequenceFilesList);
  56 @SequenceFilesList = ExpandFileNames(\@ARGV, "aln msf fasta fta pir");
  57 
  58 print "Processing options...\n";
  59 my(%OptionsInfo);
  60 ProcessOptions();
  61 
  62 # Set up information about input files...
  63 print "Checking input sequence file(s)...\n";
  64 my(%SequenceFilesInfo);
  65 RetrieveSequenceFilesInfo();
  66 SetupSequenceRegionsData();
  67 
  68 # Process input files..
  69 my($FileIndex);
  70 if (@SequenceFilesList > 1) {
  71   print "\nProcessing sequence files...\n";
  72 }
  73 for $FileIndex (0 .. $#SequenceFilesList) {
  74   if ($SequenceFilesInfo{FilesOkay}[$FileIndex]) {
  75     print "\nProcessing file $SequenceFilesList[$FileIndex]...\n";
  76     AnalyzeSequenceFileData($FileIndex);
  77   }
  78 }
  79 print "\n$ScriptName:Done...\n\n";
  80 
  81 $EndTime = new Benchmark;
  82 $TotalTime = timediff ($EndTime, $StartTime);
  83 print "Total time: ", timestr($TotalTime), "\n";
  84 
  85 ###############################################################################
  86 
  87 # Analyze sequence file...
  88 sub AnalyzeSequenceFileData {
  89   my($FileIndex) = @_;
  90   my($SequenceFile, $SequenceDataRef);
  91 
  92   $SequenceFile = $SequenceFilesList[$FileIndex];
  93 
  94   open SEQUENCEFILE, "$SequenceFile" or die "Error: Can't open $SequenceFile: $! \n";
  95   $SequenceDataRef = ReadSequenceFile($SequenceFile);
  96   close SEQUENCEFILE;
  97 
  98   if ($OptionsInfo{CalculatePercentIdentityMatrix}) {
  99     CalculatePercentIdentityMatrix($FileIndex, $SequenceDataRef);
 100   }
 101   if ($OptionsInfo{PerformResidueFrequencyAnalysis}) {
 102     PerformResidueFrequencyAnalysis($FileIndex, $SequenceDataRef);
 103   }
 104 }
 105 
 106 # Calculate percent identity matrix...
 107 sub CalculatePercentIdentityMatrix {
 108   my($FileIndex, $SequenceDataRef) = @_;
 109   my($PercentIdentity, $PercentIdentityMatrixFile, $PercentIdentityMatrixRef, $RowID, $ColID, $Line, @LineWords);
 110 
 111   $PercentIdentityMatrixFile = $SequenceFilesInfo{OutFileRoot}[$FileIndex] . 'PercentIdentityMatrix.' . $SequenceFilesInfo{OutFileExt}[$FileIndex];
 112   $PercentIdentityMatrixRef = CalculatePercentSequenceIdentityMatrix($SequenceDataRef, $OptionsInfo{IgnoreGaps}, $OptionsInfo{Precision});
 113 
 114   print "Generating percent identity matrix file $PercentIdentityMatrixFile...\n";
 115   open OUTFILE, ">$PercentIdentityMatrixFile" or die "Can't open $PercentIdentityMatrixFile: $!\n";
 116 
 117   # Write out column labels...
 118   @LineWords = ();
 119   push @LineWords, '';
 120   for $ColID (@{$PercentIdentityMatrixRef->{IDs}}) {
 121     push @LineWords, $ColID;
 122   }
 123   $Line = JoinWords(\@LineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 124   print OUTFILE "$Line\n";
 125 
 126   # Write out rows...
 127   for $RowID (@{$PercentIdentityMatrixRef->{IDs}}) {
 128     @LineWords = ();
 129     push @LineWords, $RowID;
 130     for $ColID (@{$PercentIdentityMatrixRef->{IDs}}) {
 131       $PercentIdentity = $PercentIdentityMatrixRef->{PercentIdentity}{$RowID}{$ColID};
 132       push @LineWords, $PercentIdentity;
 133     }
 134     $Line = JoinWords(\@LineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 135     print OUTFILE "$Line\n";
 136   }
 137   close OUTFILE;
 138 }
 139 
 140 # Perform frequency analysis...
 141 sub PerformResidueFrequencyAnalysis {
 142   my($FileIndex, $SequenceDataRef) = @_;
 143 
 144   CountResiduesInRegions($FileIndex, $SequenceDataRef);
 145   CalculatePercentResidueFrequencyInRegions($FileIndex, $SequenceDataRef);
 146   GeneratePercentResidueFrequencyOutFilesForRegions($FileIndex, $SequenceDataRef);
 147 }
 148 
 149 # Count residues...
 150 sub CountResiduesInRegions {
 151   my($FileIndex, $SequenceDataRef) = @_;
 152 
 153   # Setup rerfernce sequence data...
 154   my($RefereceSequenceID, $RefereceSequence);
 155   $RefereceSequenceID = $SequenceFilesInfo{RefereceSequenceID}[$FileIndex];
 156   $RefereceSequence = $SequenceFilesInfo{RefereceSequence}[$FileIndex];
 157 
 158   # Count residues...
 159   my($RegionID, $StartResNum, $EndResNum, $ResNum, $ResIndex, $ID, $Sequence, $Residue);
 160   for $RegionID (@{$SequenceFilesInfo{RegionsData}[$FileIndex]{RegionIDs}}) {
 161     $StartResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{StartResNum};
 162     $EndResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{EndResNum};
 163     RESNUM: for $ResNum ($StartResNum .. $EndResNum) {
 164       $ResIndex = $ResNum - 1;
 165       if ($OptionsInfo{IgnoreGaps} && $SequenceFilesInfo{RefereceSequenceResNums}[$FileIndex]{IsGap}{$ResNum}) {
 166         next RESNUM;
 167       }
 168       # Go over residues in column $ResNum in all the sequences...
 169       ID: for $ID (@{$SequenceDataRef->{IDs}}) {
 170         $Sequence = $SequenceDataRef->{Sequence}{$ID};
 171         $Residue = substr($Sequence, $ResIndex, 1);
 172         if (IsGapResidue($Residue)) {
 173           $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{Gap} += 1;
 174         }
 175         else {
 176           if (exists $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{$Residue}) {
 177             $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{$Residue} += 1;
 178           }
 179           else {
 180             # Internal error...
 181             print "Warning: Ignoring residue $Residue in sequence $ID during ResidueFrequencyAnalysis calculation: Unknown residue...\n";
 182           }
 183         }
 184       }
 185     }
 186   }
 187 }
 188 
 189 # Calculate percent frequency for various residues in the sequence regions...
 190 sub CalculatePercentResidueFrequencyInRegions {
 191   my($FileIndex, $SequenceDataRef) = @_;
 192   my($RegionID, $StartResNum, $EndResNum, $ResNum, $Residue, $Count, $PercentCount, $MaxResiduesCount, $Precision);
 193 
 194   $MaxResiduesCount = $SequenceDataRef->{Count};
 195   $Precision = $OptionsInfo{Precision};
 196   for $RegionID (@{$SequenceFilesInfo{RegionsData}[$FileIndex]{RegionIDs}}) {
 197     $StartResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{StartResNum};
 198     $EndResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{EndResNum};
 199     RESNUM: for $ResNum ($StartResNum .. $EndResNum) {
 200       if ($OptionsInfo{IgnoreGaps} && $SequenceFilesInfo{RefereceSequenceResNums}[$FileIndex]{IsGap}{$ResNum}) {
 201         next RESNUM;
 202       }
 203       for $Residue (keys %{$SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}}) {
 204         $Count = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{$Residue};
 205         $PercentCount = ($Count / $MaxResiduesCount) * 100;
 206         $PercentCount = sprintf("%.${Precision}f", $PercentCount) + 0;
 207         $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{PercentCount}{$ResNum}{$Residue} = $PercentCount;
 208       }
 209     }
 210   }
 211 }
 212 
 213 # Generate output files...
 214 sub GeneratePercentResidueFrequencyOutFilesForRegions {
 215   my($FileIndex, $SequenceDataRef) = @_;
 216 
 217   # Setup rerfernce sequence data...
 218   my($RefereceSequenceID, $RefereceSequence);
 219   $RefereceSequenceID = $SequenceFilesInfo{RefereceSequenceID}[$FileIndex];
 220   $RefereceSequence = $SequenceFilesInfo{RefereceSequence}[$FileIndex];
 221 
 222   my($RegionID, $StartResNum, $EndResNum, $ResNum, $Count, $PercentCount, $Residue, $RegionNum, $RegionOutFile, $PercentRegionOutFile, $OutFileRoot, $OutFileExt, $Line, @LineWords, @PercentLineWords);
 223 
 224   $OutFileRoot = $SequenceFilesInfo{OutFileRoot}[$FileIndex];
 225   $OutFileExt = $SequenceFilesInfo{OutFileExt}[$FileIndex];
 226   $RegionNum = 0;
 227   for $RegionID (@{$SequenceFilesInfo{RegionsData}[$FileIndex]{RegionIDs}}) {
 228     $RegionNum++;
 229     $StartResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{StartResNum};
 230     $EndResNum = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{EndResNum};
 231 
 232     $RegionOutFile = "${OutFileRoot}ResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt}";
 233     $PercentRegionOutFile = "${OutFileRoot}PercentResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt}";
 234 
 235     print "Generating $RegionOutFile and $PercentRegionOutFile...\n";
 236     open REGIONOUTFILE, ">$RegionOutFile" or die "Error: Can't open $RegionOutFile: $! \n";
 237     open PERCENTREGIONOUTFILE, ">$PercentRegionOutFile" or die "Error: Can't open $PercentRegionOutFile: $! \n";
 238 
 239     # Write out reference residue positions as column values....
 240     @LineWords = ();
 241     push @LineWords, '';
 242     RESNUM: for $ResNum ($StartResNum .. $EndResNum) {
 243       if ($OptionsInfo{IgnoreGaps} && $SequenceFilesInfo{RefereceSequenceResNums}[$FileIndex]{IsGap}{$ResNum}) {
 244         next RESNUM;
 245       }
 246       push @LineWords, $ResNum;
 247     }
 248     $Line = JoinWords(\@LineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 249     print REGIONOUTFILE "$Line\n";
 250     print PERCENTREGIONOUTFILE "$Line\n";
 251 
 252 
 253     # Write out row data for each residue; Gap residue is written last...
 254     RESIDUE: for $Residue (sort @{$SequenceFilesInfo{ResidueCodes}[$FileIndex]}) {
 255       if ($Residue =~ /^Gap$/i) {
 256         next RESIDUE;
 257       }
 258       @LineWords = ();
 259       push @LineWords, $Residue;
 260       @PercentLineWords = ();
 261       push @PercentLineWords, $Residue;
 262 
 263       RESNUM: for $ResNum ($StartResNum .. $EndResNum) {
 264         if ($OptionsInfo{IgnoreGaps} && $SequenceFilesInfo{RefereceSequenceResNums}[$FileIndex]{IsGap}{$ResNum}) {
 265           next RESNUM;
 266         }
 267         $Count = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{$Residue};
 268         push @LineWords, $Count;
 269         $PercentCount = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{PercentCount}{$ResNum}{$Residue};
 270         push @PercentLineWords, $PercentCount;
 271       }
 272       $Line = JoinWords(\@LineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 273       print REGIONOUTFILE "$Line\n";
 274 
 275       $Line = JoinWords(\@PercentLineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 276       print PERCENTREGIONOUTFILE "$Line\n";
 277     }
 278 
 279     # Write out data for gap...
 280     $Residue = 'Gap';
 281     @LineWords = ();
 282     push @LineWords, $Residue;
 283     @PercentLineWords = ();
 284     push @PercentLineWords, $Residue;
 285 
 286     RESNUM: for $ResNum ($StartResNum .. $EndResNum) {
 287       if ($OptionsInfo{IgnoreGaps} && $SequenceFilesInfo{RefereceSequenceResNums}[$FileIndex]{IsGap}{$ResNum}) {
 288         next RESNUM;
 289       }
 290       $Count = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{Count}{$ResNum}{$Residue};
 291       push @LineWords, $Count;
 292 
 293       $PercentCount = $SequenceFilesInfo{RegionsData}[$FileIndex]{$RegionID}{PercentCount}{$ResNum}{$Residue};
 294       push @PercentLineWords, $PercentCount;
 295     }
 296     $Line = JoinWords(\@LineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 297     print REGIONOUTFILE "$Line\n";
 298 
 299     $Line = JoinWords(\@PercentLineWords, $OptionsInfo{OutDelim}, $OptionsInfo{OutQuote});
 300     print PERCENTREGIONOUTFILE "$Line\n";
 301 
 302     close REGIONOUTFILE;
 303     close PERCENTREGIONOUTFILE;
 304   }
 305 }
 306 
 307 # Process option values...
 308 sub ProcessOptions {
 309   %OptionsInfo = ();
 310 
 311   # Setup analysis mode...
 312   $OptionsInfo{CalculatePercentIdentityMatrix} = ($Options{mode} =~ /^(PercentIdentityMatrix|All)$/i) ? 1 : 0;
 313   $OptionsInfo{PerformResidueFrequencyAnalysis} = ($Options{mode} =~ /^(ResidueFrequencyAnalysis|All)$/i) ? 1 : 0;
 314 
 315   # Setup delimiter and quotes...
 316   $OptionsInfo{OutDelim} = ($Options{outdelim} =~ /tab/i ) ? "\t" : (($Options{outdelim} =~ /semicolon/i) ? "\;" : "\,");
 317   $OptionsInfo{OutQuote} = ($Options{quote} =~ /yes/i ) ? 1 : 0;
 318 
 319   # Setup reference sequence and regions for residue frequence analysis...
 320   $OptionsInfo{SpecifiedRefereceSequence} = $Options{referencesequence};
 321   $OptionsInfo{SpecifiedRegion} = $Options{region};
 322   @{$OptionsInfo{SpecifiedRegions}} = ();
 323 
 324   my(@SpecifiedRegions);
 325   @SpecifiedRegions = ();
 326   if ($Options{region} =~ /\,/) {
 327     @SpecifiedRegions = split /\,/, $OptionsInfo{SpecifiedRegion};
 328     if (@SpecifiedRegions % 2) {
 329       die "Error: The value specified, $Options{region}, for option \"--region\" is not valid. Allowed values: \"StartResNum,EndResNum,[StartResNum,EndResNum...]\" or UseCompleteSequence\n";
 330     }
 331     # Make sure EndResNum > StartResNum...
 332     my($StartResNum, $EndResNum, $Index, $RegionNum);
 333     $RegionNum = 0;
 334     for ($Index = 0; $Index <= $#SpecifiedRegions; $Index += 2) {
 335       $StartResNum = $SpecifiedRegions[$Index];
 336       $EndResNum = $SpecifiedRegions[$Index + 1];
 337       $RegionNum++;
 338       if (!IsPositiveInteger($StartResNum)) {
 339         die "Error: The value specified, $Options{region}, for option \"--region\" is not valid: The start residue number, $StartResNum, must be a positive integer for region $RegionNum.\n";
 340       }
 341       if (!IsPositiveInteger($EndResNum)) {
 342         die "Error: The value specified, $Options{region}, for option \"--region\" is not valid: The start residue number, $EndResNum, must be a positive integer for region $RegionNum.\n";
 343       }
 344       if ($StartResNum >= $EndResNum) {
 345         die "Error: The value specified, $Options{region}, for option \"--region\" is not valid: The start residue number, $StartResNum, must be smaller than end residue number, $EndResNum, for region $RegionNum.\n";
 346       }
 347     }
 348   }
 349   else {
 350     if ($Options{region} !~ /^UseCompleteSequence$/i) {
 351       die "Error: The value specified, $Options{region}, for option \"--region\" is not valid. Allowed values: \"StartResNum,EndResNum,[StartResNum,EndResNum...]\" or UseCompleteSequence\n";
 352     }
 353   }
 354   push @{$OptionsInfo{SpecifiedRegions}}, @SpecifiedRegions;
 355 
 356   # Miscellaneous options...
 357   $OptionsInfo{Precision} = $Options{precision};
 358   $OptionsInfo{IgnoreGaps} = ($Options{ignoregaps} =~ /Yes/i) ? 1 : 0;
 359   $OptionsInfo{RegionResiduesMode} = $Options{regionresiduesmode};
 360 
 361   $OptionsInfo{OverwriteFiles} = $Options{overwrite} ? 1 : 0;
 362   $OptionsInfo{OutFileRoot} = $Options{root} ? $Options{root} : 0;
 363 }
 364 
 365 # Retrieve information about sequence files...
 366 sub RetrieveSequenceFilesInfo {
 367   my($Index, $SequenceFile, $FileSupported, $FileFormat, $SequenceCount, $RefereceSequence, $RefereceSequenceID, $RefereceSequenceLen, $RefereceSequenceWithNoGaps, $RefereceSequenceWithNoGapsLen, $RefereceSequenceRegionCount, $FileDir, $FileName, $FileExt, $OutFileRoot, $OutFileExt, $SequenceDataRef, $SpecifiedRefereceSequence, @SpecifiedRegions, @RefereceSequenceRegions);
 368 
 369   %SequenceFilesInfo = ();
 370   @{$SequenceFilesInfo{FilesOkay}} = ();
 371   @{$SequenceFilesInfo{OutFileRoot}} = ();
 372   @{$SequenceFilesInfo{OutFileExt}} = ();
 373   @{$SequenceFilesInfo{Format}} = ();
 374   @{$SequenceFilesInfo{SequenceCount}} = ();
 375   @{$SequenceFilesInfo{RefereceSequenceID}} = ();
 376   @{$SequenceFilesInfo{RefereceSequence}} = ();
 377   @{$SequenceFilesInfo{RefereceSequenceLen}} = ();
 378   @{$SequenceFilesInfo{RefereceSequenceWithNoGaps}} = ();
 379   @{$SequenceFilesInfo{RefereceSequenceWithNoGapsLen}} = ();
 380   @{$SequenceFilesInfo{RefereceSequenceRegions}} = ();
 381   @{$SequenceFilesInfo{RefereceSequenceRegionCount}} = ();
 382   @{$SequenceFilesInfo{ResidueCodes}} = ();
 383 
 384   FILELIST: for $Index (0 .. $#SequenceFilesList) {
 385     $SequenceFile = $SequenceFilesList[$Index];
 386     $SequenceFilesInfo{FilesOkay}[$Index] = 0;
 387     $SequenceFilesInfo{OutFileRoot}[$Index] = '';
 388     $SequenceFilesInfo{OutFileExt}[$Index] = '';
 389     $SequenceFilesInfo{Format}[$Index] = 'NotSupported';
 390     $SequenceFilesInfo{SequenceCount}[$Index] = 0;
 391     $SequenceFilesInfo{RefereceSequenceID}[$Index] = '';
 392     $SequenceFilesInfo{RefereceSequence}[$Index] = '';
 393     $SequenceFilesInfo{RefereceSequenceLen}[$Index] = '';
 394     $SequenceFilesInfo{RefereceSequenceWithNoGaps}[$Index] = '';
 395     $SequenceFilesInfo{RefereceSequenceWithNoGapsLen}[$Index] = '';
 396     @{$SequenceFilesInfo{RefereceSequenceRegions}[$Index]} = ();
 397     $SequenceFilesInfo{RefereceSequenceRegionCount}[$Index] = 0;
 398     @{$SequenceFilesInfo{ResidueCodes}[$Index]} = ();
 399 
 400     if (! open SEQUENCEFILE, "$SequenceFile") {
 401       warn "Warning: Ignoring file $SequenceFile: Couldn't open it: $! \n";
 402       next FILELIST;
 403     }
 404     close SEQUENCEFILE;
 405 
 406     ($FileSupported, $FileFormat) = IsSupportedSequenceFile($SequenceFile);
 407     if (!$FileSupported) {
 408       warn "Warning: Ignoring file $SequenceFile: Sequence file format is not supported.\n";
 409       next FILELIST;
 410     }
 411 
 412     $SequenceDataRef = ReadSequenceFile($SequenceFile);
 413 
 414     $SequenceCount = $SequenceDataRef->{Count};
 415     if (!$SequenceCount) {
 416       warn "Warning: Ignoring file $SequenceFile: Sequence data is missing.\n";
 417       next FILELIST;
 418     }
 419 
 420     # Make sure all sequence lengths are identical...
 421     if (!AreSequenceLengthsIdentical($SequenceDataRef)) {
 422       warn "Warning: Ignoring file $SequenceFile: Sequence legths are not identical.\n";
 423       next FILELIST;
 424     }
 425     $SpecifiedRefereceSequence = $OptionsInfo{SpecifiedRefereceSequence};
 426     # Make sure reference sequence ID is valid...
 427     if ($SpecifiedRefereceSequence =~ /^UseFirstSequenceID$/i) {
 428       $RefereceSequenceID = $SequenceDataRef->{IDs}[0];
 429     }
 430     else {
 431       if (!exists($SequenceDataRef->{Sequence}{$SpecifiedRefereceSequence})) {
 432         warn "Warning: Ignoring file $SequenceFile: Rreference sequence ID, $SpecifiedRefereceSequence, specified using option \"--referencesequence\" is missing.\n";
 433         next FILELIST;
 434       }
 435       $RefereceSequenceID = $SpecifiedRefereceSequence;
 436     }
 437 
 438     # Make sure sequence regions corresponding to reference sequence are valid...
 439     @RefereceSequenceRegions = ();
 440     $RefereceSequenceRegionCount = 0;
 441     $RefereceSequence = $SequenceDataRef->{Sequence}{$RefereceSequenceID};
 442     $RefereceSequenceLen = length $RefereceSequence;
 443 
 444     $RefereceSequenceWithNoGaps = RemoveSequenceGaps($RefereceSequence);
 445     $RefereceSequenceWithNoGapsLen = length $RefereceSequenceWithNoGaps;
 446 
 447     @SpecifiedRegions = ();
 448     push @SpecifiedRegions, @{$OptionsInfo{SpecifiedRegions}};
 449     if (@SpecifiedRegions) {
 450       # Make sure specified regions are valid...
 451       my($StartResNum, $EndResNum, $RegionIndex, $RegionNum);
 452       $RegionNum = 0;
 453       for ($RegionIndex = 0; $RegionIndex <= $#SpecifiedRegions; $RegionIndex += 2) {
 454         $StartResNum = $SpecifiedRegions[$RegionIndex];
 455         $EndResNum = $SpecifiedRegions[$RegionIndex + 1];
 456         $RegionNum++;
 457         if ($OptionsInfo{IgnoreGaps}) {
 458           if ($StartResNum > $RefereceSequenceWithNoGapsLen) {
 459             warn "Warning: Ignoring file $SequenceFile: The value specified, $Options{region}, for option \"--region\" is not valid: The start residue number, $StartResNum, must be smaller the sequence length, $RefereceSequenceWithNoGapsLen, of reference sequence ID,  $RefereceSequenceID, in region $RegionNum. The reference sequence residue numbers correspond to the sequence with no gaps. Specify \"No\" value for \"-i, --ignoregaps\" option to use residue numbers corresponding to reference sequence including gaps.\n";
 460             next FILELIST;
 461           }
 462           if ($EndResNum > $RefereceSequenceWithNoGapsLen) {
 463             warn "Warning: Ignoring file $SequenceFile: The value specified, $Options{region}, for option \"--region\" is not valid: The end residue number, $EndResNum, must be smaller the sequence length, $RefereceSequenceWithNoGapsLen, of reference sequence ID,  $RefereceSequenceID, in region $RegionNum. The reference sequence residue numbers correspond to the sequence with no gaps. Specify \"No\" value for \"-i, --ignoregaps\" option to use residue numbers corresponding to reference sequence including gaps.\n";
 464             next FILELIST;
 465           }
 466         }
 467         else {
 468           if ($StartResNum > $RefereceSequenceLen) {
 469             warn "Warning: Ignoring file $SequenceFile: The value specified, $Options{region}, for option \"--region\" is not valid: The start residue number, $StartResNum, must be smaller the sequence length, $RefereceSequenceLen, of reference sequence ID,  $RefereceSequenceID, in region $RegionNum.\n";
 470             next FILELIST;
 471           }
 472           if ($EndResNum > $RefereceSequenceLen) {
 473             warn "Warning: Ignoring file $SequenceFile: The value specified, $Options{region}, for option \"--region\" is not valid: The end residue number, $EndResNum, must be smaller the sequence length, $RefereceSequenceLen, of reference sequence ID,  $RefereceSequenceID, in region $RegionNum.\n";
 474             next FILELIST;
 475           }
 476         }
 477         push @RefereceSequenceRegions, ($StartResNum, $EndResNum);
 478       }
 479       $RefereceSequenceRegionCount = $RegionNum;
 480     }
 481     else {
 482       # Use complete sequence corresponding to reference sequence ID...
 483       if ($OptionsInfo{IgnoreGaps}) {
 484         push @RefereceSequenceRegions, (1, $RefereceSequenceWithNoGapsLen);
 485       }
 486       else {
 487         push @RefereceSequenceRegions, (1, $RefereceSequenceLen);
 488       }
 489       $RefereceSequenceRegionCount = 1;
 490     }
 491     # Setup output file names...
 492     $FileDir = ""; $FileName = ""; $FileExt = "";
 493     ($FileDir, $FileName, $FileExt) = ParseFileName($SequenceFile);
 494     $FileExt = "csv";
 495     if ($Options{outdelim} =~ /^tab$/i) {
 496       $FileExt = "tsv";
 497     }
 498     $OutFileExt = $FileExt;
 499     if ($OptionsInfo{OutFileRoot} && (@SequenceFilesList == 1)) {
 500       my ($RootFileDir, $RootFileName, $RootFileExt) = ParseFileName($OptionsInfo{OutFileRoot});
 501       if ($RootFileName && $RootFileExt) {
 502         $FileName = $RootFileName;
 503       }
 504       else {
 505         $FileName = $OptionsInfo{OutFileRoot};
 506       }
 507       $OutFileRoot = $FileName;
 508     }
 509     else {
 510       $OutFileRoot = $FileName;
 511     }
 512     if (!$OptionsInfo{OverwriteFiles}) {
 513       if ($OptionsInfo{CalculatePercentIdentityMatrix}) {
 514         if (-e "${OutFileRoot}PercentIdentityMatrix.${OutFileExt}") {
 515           warn "Warning: Ignoring file $SequenceFile: The file ${OutFileRoot}PercentIdentityMatrix.${OutFileExt} already exists\n";
 516           next FILELIST;
 517         }
 518       }
 519       if ($OptionsInfo{PerformResidueFrequencyAnalysis}) {
 520         my($RegionNum);
 521         for $RegionNum (1 .. $RefereceSequenceRegionCount) {
 522           if (-e "${OutFileRoot}ResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt}") {
 523             warn "Warning: Ignoring file $SequenceFile: The file ${OutFileRoot}ResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt} already exists\n";
 524             next FILELIST;
 525           }
 526           if (-e "${OutFileRoot}PercentResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt}") {
 527             warn "Warning: Ignoring file $SequenceFile: The file ${OutFileRoot}PercentResidueFrequencyAnalysisRegion${RegionNum}.${OutFileExt} already exists\n";
 528             next FILELIST;
 529           }
 530         }
 531       }
 532     }
 533 
 534     $SequenceFilesInfo{FilesOkay}[$Index] = 1;
 535     $SequenceFilesInfo{OutFileRoot}[$Index] = $OutFileRoot;
 536     $SequenceFilesInfo{OutFileExt}[$Index] = $OutFileExt;
 537 
 538     $SequenceFilesInfo{Format}[$Index] = $FileFormat;
 539     $SequenceFilesInfo{SequenceCount}[$Index] = $SequenceCount;
 540     $SequenceFilesInfo{RefereceSequenceID}[$Index] = $RefereceSequenceID;
 541     $SequenceFilesInfo{RefereceSequence}[$Index] = $RefereceSequence;
 542     $SequenceFilesInfo{RefereceSequenceLen}[$Index] = $RefereceSequenceLen;
 543     $SequenceFilesInfo{RefereceSequenceWithNoGaps}[$Index] = $RefereceSequenceWithNoGaps;
 544     $SequenceFilesInfo{RefereceSequenceWithNoGapsLen}[$Index] = $RefereceSequenceWithNoGapsLen;
 545     push @{$SequenceFilesInfo{RefereceSequenceRegions}[$Index]}, @RefereceSequenceRegions;
 546     $SequenceFilesInfo{RefereceSequenceRegionCount}[$Index] = $RefereceSequenceRegionCount;
 547 
 548     # Setup residue codes...
 549     SetupSequenceFileResidueCodes($SequenceDataRef, $Index);
 550   }
 551 }
 552 
 553 sub SetupSequenceFileResidueCodes {
 554   my($SequenceDataRef, $FileIndex) = @_;
 555   my($Residue, @ResidueCodesList, %ResidueCodesMap);
 556 
 557   # Initialize
 558   @{$SequenceFilesInfo{ResidueCodes}[$FileIndex]} = ();
 559 
 560   %ResidueCodesMap = ();
 561   @ResidueCodesList = ();
 562 
 563   # Setup standard amino acids and nucleic acids one letter codes...
 564   if ($OptionsInfo{RegionResiduesMode} =~ /^AminoAcids$/i) {
 565     @ResidueCodesList = AminoAcids::GetAminoAcids('OneLetterCode');
 566   }
 567   elsif ($OptionsInfo{RegionResiduesMode} =~ /^NucleicAcids$/i) {
 568     push @ResidueCodesList, ('A', 'G', 'T', 'U', 'C');
 569   }
 570   push @ResidueCodesList, 'Gap';
 571   for $Residue (@ResidueCodesList) {
 572     $ResidueCodesMap{$Residue} = $Residue;
 573   }
 574 
 575   # Go over all the residues in all the sequences and add missing ones to the list...
 576   my($ID, $Sequence, $ResIndex);
 577   for $ID (@{$SequenceDataRef->{IDs}}) {
 578     $Sequence = $SequenceDataRef->{Sequence}{$ID};
 579     RES: for $ResIndex (0 .. (length($Sequence) - 1)) {
 580       $Residue = substr($Sequence, $ResIndex, 1);
 581       if (IsGapResidue($Residue)) {
 582         next RES;
 583       }
 584       if (exists $ResidueCodesMap{$Residue}) {
 585         next RES;
 586       }
 587       push @ResidueCodesList, $Residue;
 588       $ResidueCodesMap{$Residue} = $Residue;
 589     }
 590   }
 591   push @{$SequenceFilesInfo{ResidueCodes}[$FileIndex]}, @ResidueCodesList;
 592 }
 593 
 594 # Setup regions data for performing residue frequency analysis...
 595 sub SetupSequenceRegionsData {
 596   my($Index, $RefereceSequence, $RefereceSequenceLen, $RegionID, $StartResNum, $EndResNum, $RegionIndex, $RegionNum, $NoGapResNum, $ResNum, $ResIndex, $Residue, $ResidueCode, @RefereceSequenceRegions);
 597 
 598 
 599   @{$SequenceFilesInfo{RefereceSequenceResNums}} = ();
 600   @{$SequenceFilesInfo{RegionsData}} = ();
 601 
 602   FILELIST: for $Index (0 .. $#SequenceFilesList) {
 603     %{$SequenceFilesInfo{RefereceSequenceResNums}[$Index]{IsGap}} = ();
 604     %{$SequenceFilesInfo{RefereceSequenceResNums}[$Index]{NoGapToGap}} = ();
 605     %{$SequenceFilesInfo{RefereceSequenceResNums}[$Index]{GapToNoGap}} = ();
 606     %{$SequenceFilesInfo{RegionsData}[$Index]} = ();
 607     @{$SequenceFilesInfo{RegionsData}[$Index]{RegionIDs}} = ();
 608 
 609     if (!$SequenceFilesInfo{FilesOkay}[$Index]) {
 610       next FILELIST;
 611     }
 612     if (!$OptionsInfo{PerformResidueFrequencyAnalysis}) {
 613       next FILELIST;
 614     }
 615 
 616     $RefereceSequence = $SequenceFilesInfo{RefereceSequence}[$Index];
 617     $RefereceSequenceLen = $SequenceFilesInfo{RefereceSequenceLen}[$Index];
 618 
 619     # Setup residue number mapping and gap status for refernece sequence...
 620     $NoGapResNum = 0;
 621     $ResNum = 0;
 622     for $ResIndex (0 .. ($RefereceSequenceLen - 1)) {
 623       $ResNum++;
 624       $Residue = substr($RefereceSequence, $ResIndex, 1);
 625       $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{IsGap}{$ResNum} = 1;
 626       if (!IsGapResidue($Residue)) {
 627         $NoGapResNum++;
 628         $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{IsGap}{$ResNum} = 0;
 629         $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{NoGapToGap}{$NoGapResNum} = $ResNum;
 630         $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{GapToNoGap}{$ResNum} = $NoGapResNum;
 631       }
 632     }
 633     # Map residue numbers for specified regions to the reference residue in input sequence/alignment files
 634     $RegionNum = 0;
 635     @RefereceSequenceRegions = ();
 636     push @RefereceSequenceRegions, @{$SequenceFilesInfo{RefereceSequenceRegions}[$Index]};
 637     for ($RegionIndex = 0; $RegionIndex <= $#RefereceSequenceRegions; $RegionIndex += 2) {
 638       $StartResNum = $RefereceSequenceRegions[$RegionIndex];
 639       $EndResNum = $RefereceSequenceRegions[$RegionIndex + 1];
 640       $RegionNum++;
 641       $RegionID = "Region${RegionNum}";
 642       if ($OptionsInfo{IgnoreGaps}) {
 643         # Map residue numbers to the reference sequence residue numbers to account for any ignored gaps...
 644         $StartResNum = $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{NoGapToGap}{$StartResNum};
 645         $EndResNum = $SequenceFilesInfo{RefereceSequenceResNums}[$Index]{NoGapToGap}{$EndResNum};
 646       }
 647       push @{$SequenceFilesInfo{RegionsData}[$Index]{RegionIDs}}, $RegionID;
 648       $SequenceFilesInfo{RegionsData}[$Index]{$RegionID}{StartResNum} = $StartResNum;
 649       $SequenceFilesInfo{RegionsData}[$Index]{$RegionID}{EndResNum} = $EndResNum;
 650 
 651       # Initialize data for residue codes...
 652       for $ResNum ($StartResNum .. $EndResNum) {
 653         for $ResidueCode (@{$SequenceFilesInfo{ResidueCodes}[$Index]}) {
 654           $SequenceFilesInfo{RegionsData}[$Index]{$RegionID}{Count}{$ResNum}{$ResidueCode} = 0;
 655         }
 656       }
 657     }
 658   }
 659 }
 660 
 661 # Setup script usage  and retrieve command line arguments specified using various options...
 662 sub SetupScriptUsage {
 663 
 664   # Retrieve all the options...
 665   %Options = ();
 666   $Options{ignoregaps} = 'yes';
 667   $Options{regionresiduesmode} = 'None';
 668   $Options{mode} = 'PercentIdentityMatrix';
 669   $Options{outdelim} = 'comma';
 670   $Options{precision} = 2;
 671   $Options{quote} = 'yes';
 672   $Options{referencesequence} = 'UseFirstSequenceID';
 673   $Options{region} = 'UseCompleteSequence';
 674 
 675   if (!GetOptions(\%Options, "help|h", "ignoregaps|i=s", "mode|m=s", "outdelim=s", "overwrite|o", "precision|p=i", "quote|q=s", "referencesequence=s", "region=s", "regionresiduesmode=s", "root|r=s", "workingdir|w=s")) {
 676     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";
 677   }
 678   if ($Options{workingdir}) {
 679     if (! -d $Options{workingdir}) {
 680       die "Error: The value specified, $Options{workingdir}, for option \"-w --workingdir\" is not a directory name.\n";
 681     }
 682     chdir $Options{workingdir} or die "Error: Couldn't chdir $Options{workingdir}: $! \n";
 683   }
 684   if ($Options{ignoregaps} !~ /^(yes|no)$/i) {
 685     die "Error: The value specified, $Options{ignoregaps}, for option \"-q --quote\" is not valid. Allowed values: yes or no\n";
 686   }
 687   if ($Options{regionresiduesmode} !~ /^(AminoAcids|NucleicAcids|None)$/i) {
 688     die "Error: The value specified, $Options{regionresiduesmode}, for option \"--regionresiduesmode\" is not valid. Allowed values: AminoAcids, NucleicAcids or None\n";
 689   }
 690   if ($Options{mode} !~ /^(PercentIdentityMatrix|ResidueFrequencyAnalysis|All)$/i) {
 691     die "Error: The value specified, $Options{mode}, for option \"-m --mode\" is not valid. Allowed values: PercentIdentityMatrix, ResidueFrequencyAnalysis  or All\n";
 692   }
 693   if ($Options{outdelim} !~ /^(comma|semicolon|tab)$/i) {
 694     die "Error: The value specified, $Options{outdelim}, for option \"--outdelim\" is not valid. Allowed values: comma, tab, or semicolon\n";
 695   }
 696   if ($Options{quote} !~ /^(yes|no)$/i) {
 697     die "Error: The value specified, $Options{quote}, for option \"-q --quote\" is not valid. Allowed values: yes or no\n";
 698   }
 699   if (!IsPositiveInteger($Options{precision})) {
 700     die "Error: The value specified, $Options{precision}, for option \"-p --precision\" is not valid. Allowed values: > 0 \n";
 701   }
 702 }
 703