1 #!/usr/bin/perl -w 2 # 3 # File: ExtractFromSequenceFiles.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 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
