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