       1 #!/usr/bin/perl -w
       2 #
       3 # File: ExtractFromSequenceFiles.pl
       4 # Author: Manish Sud <msud@san.rr.com>
       5 #
       6 # Copyright (C) 2018 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
      36 my($ScriptName, %Options, $StartTime, $EndTime, $TotalTime);
      37
      38 # Autoflush STDOUT
      39 $| = 1;
      40
      41 # Starting message...
      42 $ScriptName = basename($0);
      43 print "\n$ScriptName: Starting...\n\n";
      44 $StartTime = new Benchmark;
      45
      46 # Setup script usage message...
      47 SetupScriptUsage();
      48 if ($Options{help} || @ARGV < 1) {
      49   die GetUsageFromPod("$FindBin::Bin/$ScriptName");
      50 }
      51
      52 # Expand wild card file names...
      53 my(@SequenceFilesList);
      54 @SequenceFilesList = ExpandFileNames(\@ARGV, "aln msf fasta fta pir");
      55
      56 # Process options...
      57 print "Processing options...\n";
      58 my(%OptionsInfo);
      59 ProcessOptions();
      60
      61 # Set up information about input files...
      62 print "Checking input sequence file(s)...\n";
      63 my(%SequenceFilesInfo);
      64 RetrieveSequenceFilesInfo();
      65
      66 # Process input files..
      67 my($FileIndex);
      68 if (@SequenceFilesList > 1) {
      69   print "\nProcessing sequence files...\n";
      70 }
      71 for $FileIndex (0 .. $#SequenceFilesList) {
      72   if ($SequenceFilesInfo{FilesOkay}[$FileIndex]) {
      73     print "\nProcessing file $SequenceFilesList[$FileIndex]...\n";
      74     ExtractFromSequenceFiles($FileIndex);
      75   }
      76 }
      77 print "\n$ScriptName:Done...\n\n";
      78
      79 $EndTime = new Benchmark;
      80 $TotalTime = timediff ($EndTime, $StartTime);
      81 print "Total time: ", timestr($TotalTime), "\n";
      82
      83 ###############################################################################
      84
      85 # Extract from sequence files...
      86 sub ExtractFromSequenceFiles {
      87   my($FileIndex) = @_;
      88   my($OutSequenceFile, $SequenceFile, $SequenceDataRef, $SpecifiedSequenceDataRef);
      89
      90   # Read sequence file...
      91   $SequenceFile = $SequenceFilesList[$FileIndex];
      92   open SEQUENCEFILE, "$SequenceFile" or die "Error: Can't open $SequenceFile: $! \n";
      93   $SequenceDataRef = ReadSequenceFile($SequenceFile);
      94   close SEQUENCEFILE;
      95
      96   $OutSequenceFile = $SequenceFilesInfo{OutFile}[$FileIndex];
      97   print "Generating sequence file $OutSequenceFile...\n";
      98
      99   # Retrieve sequence data for specified sequences...
     100   $SpecifiedSequenceDataRef = GetSpecifiedSequenceData($SequenceDataRef);
     101
     102   # Handle gaps...
     103   if ($OptionsInfo{IgnoreGaps}) {
     104     if (@{$SpecifiedSequenceDataRef->{IDs}} > 1) {
     105       if (AreSequenceLengthsIdentical($SpecifiedSequenceDataRef)) {
     106         $SpecifiedSequenceDataRef = RemoveSequenceAlignmentGapColumns($SpecifiedSequenceDataRef);
     107       }
     108     }
     109     else {
     110       # Remove the gaps from the sequence...
     111       my($ID, $Sequence);
     112       $ID = $SpecifiedSequenceDataRef->{IDs}[0];
     113       $Sequence = $SpecifiedSequenceDataRef->{Sequence}{$ID};
     114       $SpecifiedSequenceDataRef->{Sequence}{$ID} = RemoveSequenceGaps($Sequence);
     115     }
     116   }
     117
     118   # Write out the file...
     119   WritePearsonFastaSequenceFile($OutSequenceFile, $SpecifiedSequenceDataRef, $OptionsInfo{MaxSequenceLength});
     120 }
     121
     122 # Get specified sequence data...
     123 sub GetSpecifiedSequenceData {
     124   my($SequenceDataRef) = @_;
     125
     126   if ($OptionsInfo{Mode} =~ /^SequenceID$/i) {
     127     return GetDataBySequenceIDs($SequenceDataRef);
     128   }
     129   elsif ($Options{mode} =~ /^SequenceNum$/i) {
     130     return GetDataBySequenceNums($SequenceDataRef);
     131   }
     132   elsif ($Options{mode} =~ /^SequenceNumRange$/i) {
     133     return GetDataBySequenceNumRange($SequenceDataRef);
     134   }
     135   else {
     136     return undef;
     137   }
     138 }
     139
     140 # Get specified sequence data...
     141 sub GetDataBySequenceIDs {
     142   my($SequenceDataRef) = @_;
     143   my($ID, $SequenceCount, $IDMatched, $SpecifiedID, %SpecifiedSequenceDataMap);
     144
     145   # Go over sequences and collect sequences for writing out a new sequence file...
     146   %SpecifiedSequenceDataMap = ();
     147   @{$SpecifiedSequenceDataMap{IDs}} = ();
     148   %{$SpecifiedSequenceDataMap{Description}} = ();
     149   %{$SpecifiedSequenceDataMap{Sequence}} = ();
     150
     151   $SequenceCount = 0;
     152   ID: for $ID (@{$SequenceDataRef->{IDs}}) {
     153     if ($OptionsInfo{MatchExactSequenceIDs}) {
     154       if (!exists $OptionsInfo{SpecifiedSequenceIDsMap}{lc($ID)}) {
     155         next ID;
     156       }
     157       if ($SequenceCount >= scalar @{$OptionsInfo{SpecifiedSequenceIDs}}) {
     158         last ID;
     159       }
     160       $SequenceCount++;
     161     }
     162     else {
     163       # Does this ID contains specified ID as substring...
     164       $IDMatched = 0;
     165       SPECIFIEDID: for $SpecifiedID (@{$OptionsInfo{SpecifiedSequenceIDs}}) {
     166         if ($ID =~ /$SpecifiedID/i) {
     167           $IDMatched = 1;
     168           last SPECIFIEDID;
     169         }
     170       }
     171       if (!$IDMatched) {
     172         next ID;
     173       }
     174       $SequenceCount++;
     175     }
     176     # Collect sequence data...
     177     push @{$SpecifiedSequenceDataMap{IDs}}, $ID;
     178     $SpecifiedSequenceDataMap{Description}{$ID} = $SequenceDataRef->{Description}{$ID};
     179     $SpecifiedSequenceDataMap{Sequence}{$ID} = $SequenceDataRef->{Sequence}{$ID};
     180   }
     181
     182   return \%SpecifiedSequenceDataMap;
     183 }
     184
     185 # Get specified sequence data...
     186 sub GetDataBySequenceNums {
     187   my($SequenceDataRef) = @_;
     188   my($ID, $SequenceNum, $SequenceCount, %SpecifiedSequenceDataMap);
     189
     190   # Go over sequences and collect sequences for writing out a new sequence file...
     191   %SpecifiedSequenceDataMap = ();
     192   @{$SpecifiedSequenceDataMap{IDs}} = ();
     193   %{$SpecifiedSequenceDataMap{Description}} = ();
     194   %{$SpecifiedSequenceDataMap{Sequence}} = ();
     195
     196   $SequenceNum = 0;
     197   $SequenceCount = 0;
     198   ID: for $ID (@{$SequenceDataRef->{IDs}}) {
     199     $SequenceNum++;
     200     if (!exists $OptionsInfo{SpecifiedSequenceIDsMap}{$SequenceNum}) {
     201       next ID;
     202     }
     203     if ($SequenceCount >= scalar @{$OptionsInfo{SpecifiedSequenceIDs}}) {
     204       last ID;
     205     }
     206     $SequenceCount++;
     207
     208     # Collect sequence data...
     209     push @{$SpecifiedSequenceDataMap{IDs}}, $ID;
     210     $SpecifiedSequenceDataMap{Description}{$ID} = $SequenceDataRef->{Description}{$ID};
     211     $SpecifiedSequenceDataMap{Sequence}{$ID} = $SequenceDataRef->{Sequence}{$ID};
     212   }
     213
     214   return \%SpecifiedSequenceDataMap;
     215 }
     216
     217 # Get specified sequence data...
     218 sub GetDataBySequenceNumRange {
     219   my($SequenceDataRef) = @_;
     220   my($ID, $SequenceNum, $SequenceCount, %SpecifiedSequenceDataMap);
     221
     222   # Go over sequences and collect sequences for writing out a new sequence file...
     223   %SpecifiedSequenceDataMap = ();
     224   @{$SpecifiedSequenceDataMap{IDs}} = ();
     225   %{$SpecifiedSequenceDataMap{Description}} = ();
     226   %{$SpecifiedSequenceDataMap{Sequence}} = ();
     227
     228   $SequenceNum = 0;
     229   $SequenceCount = 0;
     230   ID: for $ID (@{$SequenceDataRef->{IDs}}) {
     231     $SequenceNum++;
     232
     233     if (!($SequenceNum >= $OptionsInfo{SpecifiedSequenceIDs}[0] && $SequenceNum <= $OptionsInfo{SpecifiedSequenceIDs}[1])) {
     234       next ID;
     235     }
     236     if ($SequenceNum > $OptionsInfo{SpecifiedSequenceIDs}[1]) {
     237       last ID;
     238     }
     239     $SequenceCount++;
     240     # Collect sequence data...
     241     push @{$SpecifiedSequenceDataMap{IDs}}, $ID;
     242     $SpecifiedSequenceDataMap{Description}{$ID} = $SequenceDataRef->{Description}{$ID};
     243     $SpecifiedSequenceDataMap{Sequence}{$ID} = $SequenceDataRef->{Sequence}{$ID};
     244   }
     245
     246   return \%SpecifiedSequenceDataMap;
     247 }
     248
     249
     250 # Process option values...
     251 sub ProcessOptions {
     252   %OptionsInfo = ();
     253
     254   # Miscellaneous options...
     255   $OptionsInfo{IgnoreGaps} = ($Options{ignoregaps} =~ /Yes/i) ? 1 : 0;
     256
     257   $OptionsInfo{Mode} = $Options{mode};
     258   $OptionsInfo{MatchExactSequenceIDs} = $Options{sequenceidmatch} =~ /Exact/i ? 1 :0;
     259
     260   # Check specified sequences value...
     261   $OptionsInfo{SpecifiedSequences} = $Options{sequences};
     262   @{$OptionsInfo{SpecifiedSequenceIDs}} = ();
     263   %{$OptionsInfo{SpecifiedSequenceIDsMap}} = ();
     264
     265   my(@SpecifiedSequenceIDs) = ();
     266   if ($Options{mode} =~ /^SequenceID$/i) {
     267     if (!$Options{sequences}) {
     268       die "Error: No value specified for option \"-s, --Sequences\" during \"SequenceID\" of \"-m, --mode\" option\n";
     269     }
     270     @SpecifiedSequenceIDs = split /\,/, $Options{sequences};
     271   }
     272   elsif ($Options{mode} =~ /^SequenceNum$/i) {
     273     if ($Options{sequences}) {
     274       @SpecifiedSequenceIDs = split /\,/, $Options{sequences};
     275       my($SequenceNum);
     276       for $SequenceNum (@SpecifiedSequenceIDs) {
     277         if (!IsPositiveInteger($SequenceNum)) {
     278           die "Error: The value specified, $SequenceNum, in \"$Options{sequences}\" for option \"-s, --Sequences\" is not valid: Valid values: > 0\n";
     279         }
     280       }
     281     }
     282     else {
     283       push @SpecifiedSequenceIDs, "1";
     284     }
     285   }
     286   elsif ($Options{mode} =~ /^SequenceNumRange$/i) {
     287     if (!$Options{sequences}) {
     288       die "Error: No value specified for option \"-s, --Sequences\" during \"SequenceNumRange\" of \"-m, --mode\" option\n";
     289     }
     290     @SpecifiedSequenceIDs = split /\,/, $Options{sequences};
     291     if (@SpecifiedSequenceIDs != 2) {
     292       die "Error: The number of values", scalar @SpecifiedSequenceIDs, " specified, $Options{sequences}, for option \"-s, --Sequences\" are not valid. Number of values must be 2 to indicate starting and ending sequence number.\n";
     293     }
     294     my($SequenceNum);
     295     for $SequenceNum (@SpecifiedSequenceIDs) {
     296       if (!IsPositiveInteger($SequenceNum)) {
     297         die "Error: The value specified, $SequenceNum, in \"$Options{sequences}\" for option \"-s, --Sequences\" is not valid: Valid values: > 0\n";
     298       }
     299     }
     300     if ($SpecifiedSequenceIDs[0] > $SpecifiedSequenceIDs[1]) {
     301       die "Error: The value specified \"$Options{sequences}\" for option \"-s, --Sequences\" are not valid: Starting sequence number $SpecifiedSequenceIDs[0] must be smaller than ending sequence number $SpecifiedSequenceIDs[1]\n";
     302     }
     303   }
     304   push @{$OptionsInfo{SpecifiedSequenceIDs}}, @SpecifiedSequenceIDs;
     305   my($SequenceID);
     306   for $SequenceID (@SpecifiedSequenceIDs) {
     307     if ($Options{mode} =~ /^SequenceID$/i) {
     308       $OptionsInfo{SpecifiedSequenceIDsMap}{lc($SequenceID)} = $SequenceID;
     309     }
     310     else {
     311       $OptionsInfo{SpecifiedSequenceIDsMap}{$SequenceID} = $SequenceID;
     312     }
     313   }
     314
     315   $OptionsInfo{MaxSequenceLength} = $Options{sequencelength};
     316   $OptionsInfo{OverwriteFiles} = $Options{overwrite} ? 1 : 0;
     317   $OptionsInfo{OutFileRoot} = $Options{root} ? $Options{root} : 0;
     318 }
     319
     320 # Retrieve information about sequence files...
     321 sub RetrieveSequenceFilesInfo {
     322   my($Index, $SequenceFile, $FileSupported, $FileFormat, $SequenceCount, $FileDir, $FileName, $FileExt, $OutFileRoot, $OutFileExt, $OutFileMode, $SequenceDataRef);
     323
     324   %SequenceFilesInfo = ();
     325   @{$SequenceFilesInfo{FilesOkay}} = ();
     326   @{$SequenceFilesInfo{OutFileRoot}} = ();
     327   @{$SequenceFilesInfo{OutFileExt}} = ();
     328   @{$SequenceFilesInfo{OutFile}} = ();
     329   @{$SequenceFilesInfo{Format}} = ();
     330   @{$SequenceFilesInfo{SequenceCount}} = ();
     331
     332   FILELIST: for $Index (0 .. $#SequenceFilesList) {
     333     $SequenceFile = $SequenceFilesList[$Index];
     334     $SequenceFilesInfo{FilesOkay}[$Index] = 0;
     335     $SequenceFilesInfo{OutFileRoot}[$Index] = '';
     336     $SequenceFilesInfo{OutFileExt}[$Index] = '';
     337     $SequenceFilesInfo{OutFile}[$Index] = '';
     338     $SequenceFilesInfo{Format}[$Index] = 'NotSupported';
     339     $SequenceFilesInfo{SequenceCount}[$Index] = 0;
     340
     341     if (! open SEQUENCEFILE, "$SequenceFile") {
     342       warn "Warning: Ignoring file $SequenceFile: Couldn't open it: $! \n";
     343       next FILELIST;
     344     }
     345     close SEQUENCEFILE;
     346
     347     ($FileSupported, $FileFormat) = IsSupportedSequenceFile($SequenceFile);
     348     if (!$FileSupported) {
     349       warn "Warning: Ignoring file $SequenceFile: Sequence file format is not supported.\n";
     350       next FILELIST;
     351     }
     352     $SequenceDataRef = ReadSequenceFile($SequenceFile);
     353
     354     $SequenceCount = $SequenceDataRef->{Count};
     355     if (!$SequenceCount) {
     356       warn "Warning: Ignoring file $SequenceFile: Sequence data is missing.\n";
     357       next FILELIST;
     358     }
     359
     360     # Setup output file names...
     361     $FileDir = ""; $FileName = ""; $FileExt = "";
     362     ($FileDir, $FileName, $FileExt) = ParseFileName($SequenceFile);
     363     $OutFileExt = 'fasta';
     364     if ($OptionsInfo{OutFileRoot} && (@SequenceFilesList == 1)) {
     365       my ($RootFileDir, $RootFileName, $RootFileExt) = ParseFileName($OptionsInfo{OutFileRoot});
     366       if ($RootFileName && $RootFileExt) {
     367         $FileName = $RootFileName;
     368       }
     369       else {
     370         $FileName = $OptionsInfo{OutFileRoot};
     371       }
     372       $OutFileRoot = $FileName;
     373     }
     374     else {
     375       $OutFileRoot = $FileName;
     376     }
     377     MODE: {
     378         if ($OptionsInfo{Mode} =~ /^SequenceID$/i) { $OutFileMode = 'SequenceID'; last MODE;}
     379         if ($OptionsInfo{Mode} =~ /^SequenceNum$/i) { $OutFileMode = 'SequenceNum'; last MODE;}
     380         if ($OptionsInfo{Mode} =~ /^SequenceNumRange$/i) { $OutFileMode = 'SequenceNumRange'; last MODE;}
     381         $OutFileMode = '';
     382     }
     383     if (!$OptionsInfo{OverwriteFiles}) {
     384       if (-e "${OutFileRoot}${OutFileMode}.${OutFileExt}") {
     385         warn "Warning: Ignoring file $SequenceFile: The file ${OutFileRoot}${OutFileMode}.${OutFileExt} already exists\n";
     386         next FILELIST;
     387       }
     388     }
     389
     390     $SequenceFilesInfo{FilesOkay}[$Index] = 1;
     391     $SequenceFilesInfo{OutFileRoot}[$Index] = $OutFileRoot;
     392     $SequenceFilesInfo{OutFileExt}[$Index] = $OutFileExt;
     393     $SequenceFilesInfo{OutFile}[$Index] = "${OutFileRoot}${OutFileMode}.${OutFileExt}";
     394
     395     $SequenceFilesInfo{Format}[$Index] = $FileFormat;
     396     $SequenceFilesInfo{SequenceCount}[$Index] = $SequenceCount;
     397   }
     398 }
     399
     400 # Setup script usage  and retrieve command line arguments specified using various options...
     401 sub SetupScriptUsage {
     402
     403   # Retrieve all the options...
     404   %Options = ();
     405   $Options{ignoregaps} = 'Yes';
     406   $Options{mode} = 'SequenceNum';
     407   $Options{sequenceidmatch} = 'Relaxed';
     408   $Options{sequencelength} = 80;
     409
     410   if (!GetOptions(\%Options, "help|h", "ignoregaps|i=s", "mode|m=s", "overwrite|o", "root|r=s", "sequences|s=s", "sequenceidmatch=s", "sequencelength=i", "workingdir|w=s")) {
     411     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";
     412   }
     413   if ($Options{workingdir}) {
     414     if (! -d $Options{workingdir}) {
     415       die "Error: The value specified, $Options{workingdir}, for option \"-w --workingdir\" is not a directory name.\n";
     416     }
     417     chdir $Options{workingdir} or die "Error: Couldn't chdir $Options{workingdir}: $! \n";
     418   }
     419   if ($Options{ignoregaps} !~ /^(yes|no)$/i) {
     420     die "Error: The value specified, $Options{ignoregaps}, for option \"-q --quote\" is not valid. Allowed values: yes or no\n";
     421   }
     422   if ($Options{mode} !~ /^(SequenceID|SequenceNum|SequenceNumRange)$/i) {
     423     die "Error: The value specified, $Options{mode}, for option \"-m --mode\" is not valid. Allowed values: SequenceID, SequenceNum, or SequenceNumRange\n";
     424   }
     425   if ($Options{sequenceidmatch} !~ /^(Exact|Relaxed)$/i) {
     426     die "Error: The value specified, $Options{sequenceidmatch}, for option \"--SequenceIDMatch\" is not valid. Allowed values: Exact or Relaxed\n";
     427   }
     428   if (!IsPositiveInteger($Options{sequencelength})) {
     429     die "Error: The value specified, $Options{sequencelength}, for option \"--SequenceLength\" is not valid. Allowed values: >0\n";
     430   }
     431 }
     432