1 package SequenceFileUtil;
   2 #
   3 # File: SequenceFileUtil.pm
   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 Exporter;
  28 use Text::ParseWords;
  29 use TextUtil;
  30 use FileUtil;
  31 
  32 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  33 
  34 @ISA = qw(Exporter);
  35 @EXPORT = qw(AreSequenceLengthsIdentical CalcuatePercentSequenceIdentity CalculatePercentSequenceIdentityMatrix GetLongestSequence GetShortestSequence GetSequenceLength IsGapResidue IsSupportedSequenceFile IsClustalWSequenceFile IsPearsonFastaSequenceFile IsMSFSequenceFile ReadSequenceFile RemoveSequenceGaps RemoveSequenceAlignmentGapColumns WritePearsonFastaSequenceFile);
  36 @EXPORT_OK = qw();
  37 
  38 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
  39 
  40 # Compare lengths of all sequences...
  41 sub AreSequenceLengthsIdentical {
  42   my($SequencesDataRef) = @_;
  43   my($Status, $ID, $FirstID, $FirstSeqLen, $FirstDifferentLenID, $SeqLen);
  44 
  45   $Status = 1;
  46   $FirstID = '';
  47   $FirstDifferentLenID = '';
  48 
  49   ID: for $ID (@{$SequencesDataRef->{IDs}}) {
  50     if (!$FirstID) {
  51       $FirstID = $ID;
  52       $FirstSeqLen = length($SequencesDataRef->{Sequence}{$ID});
  53       next ID;
  54     }
  55     $SeqLen = length($SequencesDataRef->{Sequence}{$ID});
  56     if ($SeqLen != $FirstSeqLen) {
  57       $Status = 0;
  58       $FirstDifferentLenID = $ID;
  59       last ID;
  60     }
  61   }
  62   return ($Status);
  63 }
  64 
  65 # Calculate percent identity between two sequences. By default, gaps are ignored.
  66 sub CalcuatePercentSequenceIdentity {
  67   my($Sequence1, $Sequence2, $PercentIdentity, $IgnoreGaps, $Precision);
  68 
  69   $PercentIdentity = '';
  70   $Precision = 1;
  71   $IgnoreGaps = 1;
  72   if (@_ == 4) {
  73     ($Sequence1, $Sequence2, $IgnoreGaps, $Precision) = @_;
  74   }
  75   elsif (@_ == 3) {
  76     ($Sequence1, $Sequence2, $IgnoreGaps) = @_;
  77   }
  78   elsif (@_ == 2) {
  79     ($Sequence1, $Sequence2) = @_;
  80   }
  81   else {
  82     return $PercentIdentity;
  83   }
  84   if (!(IsNotEmpty($Sequence1) && IsNotEmpty($Sequence2))) {
  85     return $PercentIdentity;
  86   }
  87   my($Index, $Identity, $Sequence1Len, $Sequence2Len, $Residue1, $Residue2, $ResMatchCount, $ResCount);
  88 
  89   $Sequence1Len = length($Sequence1);
  90   $Sequence2Len = length($Sequence2);
  91 
  92   $ResMatchCount = 0;
  93   $ResCount = 0;
  94   RESIDUE: for $Index (0 .. ($Sequence1Len - 1)) {
  95     $Residue1 = substr($Sequence1, $Index, 1);
  96     $Residue2 = ($Index < $Sequence2Len) ? substr($Sequence2, $Index, 1) : '';
  97     if ($IgnoreGaps) {
  98       if ($Residue1 !~ /[A-Z]/i || $Residue2 !~ /[A-Z]/i) {
  99         next RESIDUE;
 100       }
 101     }
 102     if ($Residue1 eq $Residue2) {
 103       $ResMatchCount++;
 104     }
 105     $ResCount++;
 106   }
 107   $Identity = $ResCount ? ($ResMatchCount/$ResCount) : 0.0;
 108   $PercentIdentity = sprintf("%.${Precision}f", ($Identity * 100));
 109 
 110   return $PercentIdentity;
 111 }
 112 
 113 # Calculate pairwise identify matrix for all the sequences and return a reference
 114 # to a hash with the following keys:
 115 #
 116 # {IDs} - Sequence IDs
 117 # {Count} - Number of IDs
 118 # {PercentIdentity}{$RowID}{$ColID} - Percent identify for a pair of sequences
 119 #
 120 sub CalculatePercentSequenceIdentityMatrix {
 121   my($SequencesDataRef, $IgnoreGaps, , $Precision, $ID, $RowID, $ColID, $RowIDSeq, $ColIDSeq, $PercentIdentity, %IdentityMatrixData);
 122 
 123   $IgnoreGaps = 1;
 124   $Precision = 1;
 125   if (@_ == 3) {
 126     ($SequencesDataRef, $IgnoreGaps, $Precision) = @_;
 127   }
 128   elsif (@_ == 2) {
 129     ($SequencesDataRef, $IgnoreGaps) = @_;
 130   }
 131   else {
 132     ($SequencesDataRef) = @_;
 133   }
 134 
 135   %IdentityMatrixData = ();
 136   @{$IdentityMatrixData{IDs}} = ();
 137   %{$IdentityMatrixData{PercentIdentity}} = ();
 138   $IdentityMatrixData{Count} = 0;
 139 
 140   for $ID (@{$SequencesDataRef->{IDs}}) {
 141     push @{$IdentityMatrixData{IDs}}, $ID;
 142     $IdentityMatrixData{Count} += 1;
 143   }
 144   # Initialize and calculate percent identity data values...
 145   for $RowID (@{$SequencesDataRef->{IDs}}) {
 146     %{$IdentityMatrixData{PercentIdentity}{$RowID}} = ();
 147     $RowIDSeq = $SequencesDataRef->{Sequence}{$RowID};
 148     for $ColID (@{$SequencesDataRef->{IDs}}) {
 149       $IdentityMatrixData{$RowID}{$ColID} = '';
 150       $ColIDSeq = $SequencesDataRef->{Sequence}{$ColID};
 151       $PercentIdentity = CalcuatePercentSequenceIdentity($RowIDSeq, $ColIDSeq, $IgnoreGaps, $Precision);
 152       $IdentityMatrixData{PercentIdentity}{$RowID}{$ColID} = $PercentIdentity;
 153     }
 154   }
 155   return \%IdentityMatrixData;
 156 }
 157 
 158 # Retrieve information about shortest sequence...
 159 sub GetShortestSequence {
 160   my($SequencesDataRef, $IgnoreGaps, $ID, $Sequence, $SeqLen, $Description);
 161 
 162   $IgnoreGaps = 1;
 163   if (@_ == 2) {
 164     ($SequencesDataRef, $IgnoreGaps) = @_;
 165   }
 166   else {
 167     ($SequencesDataRef) = @_;
 168   }
 169 
 170   ($ID, $Sequence, $SeqLen, $Description) =  _GetShortestOrLongestSequence($SequencesDataRef, 'Shortest', $IgnoreGaps);
 171   return ($ID, $Sequence, $SeqLen, $Description);
 172 }
 173 
 174 # Retrieve information about longest sequence..
 175 sub GetLongestSequence {
 176   my($SequencesDataRef, $IgnoreGaps, $ID, $Sequence, $SeqLen, $Description);
 177 
 178   $IgnoreGaps = 1;
 179   if (@_ == 2) {
 180     ($SequencesDataRef, $IgnoreGaps) = @_;
 181   }
 182   else {
 183     ($SequencesDataRef) = @_;
 184   }
 185 
 186   ($ID, $Sequence, $SeqLen, $Description) =  _GetShortestOrLongestSequence($SequencesDataRef, 'Longest', $IgnoreGaps);
 187   return ($ID, $Sequence, $SeqLen, $Description);
 188 }
 189 
 190 # Get sequence length...
 191 sub GetSequenceLength {
 192   my($Seq, $SeqLen, $IgnoreGaps);
 193 
 194   $SeqLen = ''; $IgnoreGaps = 1;
 195   if (@_ == 2) {
 196     ($Seq, $IgnoreGaps) = @_;
 197   }
 198   else {
 199     ($Seq) = @_;
 200   }
 201   if ($IgnoreGaps) {
 202     my($Index, $Residue);
 203     $SeqLen = 0;
 204     for $Index (0 .. (length($Seq) - 1)) {
 205       $Residue = substr($Seq, $Index, 1);
 206       if ($Residue =~ /[A-Z]/i) {
 207         $SeqLen++;
 208       }
 209     }
 210   }
 211   else {
 212     $SeqLen = length($Seq);
 213   }
 214 
 215   return $SeqLen;
 216 }
 217 
 218 # Is it a gap residue...
 219 sub IsGapResidue {
 220   my($Residue) = @_;
 221   my($Status);
 222 
 223   $Status = ($Residue !~ /[A-Z]/i ) ? 1 : 0;
 224 
 225   return $Status;
 226 }
 227 
 228 # Is it a supported sequence file?
 229 #
 230 # Supported seqence formats are:
 231 #
 232 # ALN/ClustalW   .aln
 233 # GCG/MSF         .msf
 234 # PILEUP/MSF     .msf
 235 # Fasts(Pearson) .fasta, .fta
 236 # NBRF/PIR         .pir
 237 #
 238 sub IsSupportedSequenceFile {
 239   my($SequenceFile) = @_;
 240   my($Status, $SequenceFormat);
 241   $Status = 0; $SequenceFormat = 'NotSupported';
 242 
 243   SEQFORMAT: {
 244       if (IsClustalWSequenceFile($SequenceFile)) {$Status = 1; $SequenceFormat = 'ClustalW'; last SEQFORMAT}
 245       if (IsPearsonFastaSequenceFile($SequenceFile)) {$Status = 1; $SequenceFormat = 'Pearson'; last SEQFORMAT}
 246       if (IsPIRFastaSequenceFile($SequenceFile)) {$Status = 1; $SequenceFormat = 'PIR'; last SEQFORMAT}
 247       if (IsMSFSequenceFile($SequenceFile)) {$Status = 1; $SequenceFormat = 'MSF'; last SEQFORMAT}
 248       $Status = 0; $SequenceFormat = 'NotSupported';
 249   }
 250   return ($Status, $SequenceFormat);
 251 }
 252 
 253 # Is it a ClustalW multiple sequence sequence file...
 254 sub IsClustalWSequenceFile {
 255   my($SequenceFile) = @_;
 256   my($Status, $Line);
 257 
 258   $Status = 0;
 259 
 260   open SEQUENCEFILE, "$SequenceFile" or die "Couldn't open $SequenceFile: $!\n";
 261   $Line = GetTextLine(\*SEQUENCEFILE);
 262   $Status = ($Line =~ /(ClustalW|Clustal W|Clustal)/i ) ? 1 : 0;
 263   close SEQUENCEFILE;
 264 
 265   return $Status;
 266 }
 267 
 268 # Is it a valid Pearson fasta sequence or alignment file?
 269 #
 270 sub IsPearsonFastaSequenceFile {
 271   my($FastaFile, $Line, $Status);
 272 
 273   ($FastaFile) = @_;
 274   $Status = 0;
 275 
 276   open FASTAFILE, "$FastaFile" or die "Couldn't open $FastaFile: $!\n";
 277   $Line = GetTextLine(\*FASTAFILE);
 278 
 279   # First line starts with > and the fourth character is not ';'; otherwise, it's
 280   # PIR FASTA format.
 281   if ($Line =~ /^>/) {
 282     my($FourthChar);
 283     $FourthChar = substr($Line, 3, 1);
 284     $Status = ($FourthChar !~ /\;/) ? 1 : 0;
 285   }
 286   close FASTAFILE;
 287 
 288   return $Status;
 289 }
 290 
 291 # Is it a valid NBRF/PIR fasta sequence or alignment file?
 292 #
 293 sub IsPIRFastaSequenceFile {
 294   my($FastaFile, $Line, $Status);
 295 
 296   ($FastaFile) = @_;
 297   $Status = 0;
 298 
 299   open FASTAFILE, "$FastaFile" or die "Couldn't open $FastaFile: $!\n";
 300   $Line = GetTextLine(\*FASTAFILE);
 301 
 302   # First line starts with > and the fourth character is ';'; otherwise, it's
 303   # a Pearson FASTA format.
 304   if ($Line =~ /^>/) {
 305     my($FourthChar);
 306     $FourthChar = substr($Line, 3, 1);
 307     $Status = ($FourthChar =~ /\;/) ? 1 : 0;
 308   }
 309   close FASTAFILE;
 310 
 311   return $Status;
 312 }
 313 
 314 # Is it a valid MSF sequence or alignment file?
 315 #
 316 sub IsMSFSequenceFile {
 317   my($MSFFile) = @_;
 318 
 319   open MSFFILE, "$MSFFile" or die "Couldn't open $MSFFile: $!\n";
 320 
 321   my($Line, $Status);
 322 
 323   $Status = 0;
 324   # Find a line that contains MSF: keyword and ends with '..'
 325   LINE: while ($Line = GetTextLine(\*MSFFILE)) {
 326     $Line = RemoveLeadingWhiteSpaces($Line);
 327     if ($Line =~ /MSF:/i && $Line =~ /\.\.[ ]*$/) {
 328       $Status = 1;
 329       last LINE;
 330     }
 331     elsif ($Line =~ /(!!AA_MULTIPLE_ALIGNMENT|!!NA_MULTIPLE_ALIGNMENT|PILEUP)/i) {
 332       # Pileup MSF...
 333       $Status = 1;
 334       last LINE;
 335     }
 336   }
 337   close MSFFILE;
 338 
 339   return $Status;
 340 }
 341 
 342 # Read sequence or sequence alignment file...
 343 sub ReadSequenceFile {
 344   my($SequenceFile) = @_;
 345 
 346   if (IsPearsonFastaSequenceFile($SequenceFile)) {
 347     return ReadPearsonFastaSequenceFile($SequenceFile);
 348   }
 349   elsif (IsPIRFastaSequenceFile($SequenceFile)) {
 350     return ReadPIRFastaSequenceFile($SequenceFile);
 351   }
 352   elsif (IsMSFSequenceFile($SequenceFile)) {
 353     return ReadMSFSequenceFile($SequenceFile);
 354   }
 355   elsif (IsClustalWSequenceFile($SequenceFile)) {
 356     return ReadClustalWSequenceFile($SequenceFile);
 357   }
 358   else {
 359     return undef;
 360   }
 361 }
 362 
 363 # Read file and setup alignment data...
 364 sub ReadClustalWSequenceFile {
 365   my($SequenceFile) = @_;
 366 
 367   return _ReadFileAndSetupSequencesData($SequenceFile, 'ClustalW');
 368 }
 369 
 370 # Read file and setup alignment data...
 371 sub ReadPearsonFastaSequenceFile {
 372   my($SequenceFile) = @_;
 373 
 374   return _ReadFileAndSetupSequencesData($SequenceFile, 'Pearson');
 375 }
 376 
 377 # Read file and setup alignment data...
 378 sub ReadPIRFastaSequenceFile {
 379   my($SequenceFile) = @_;
 380 
 381   return _ReadFileAndSetupSequencesData($SequenceFile, 'PIR');
 382 }
 383 
 384 
 385 # Read file and setup sequence data...
 386 sub ReadMSFSequenceFile {
 387   my($SequenceFile) = @_;
 388 
 389   return _ReadFileAndSetupSequencesData($SequenceFile, 'MSF');
 390 }
 391 
 392 # Write out a Pearson FASTA file...
 393 sub WritePearsonFastaSequenceFile {
 394   my($SequenceFileName, $SequenceDataRef, $MaxLength, $ID, $Description, $Sequence, $WrappedSequence);
 395 
 396   $MaxLength = 80;
 397   if (@_ == 3) {
 398     ($SequenceFileName, $SequenceDataRef, $MaxLength) = @_;
 399   }
 400   elsif (@_ == 2) {
 401     ($SequenceFileName, $SequenceDataRef) = @_;
 402   }
 403   open SEQUENCEFILE, ">$SequenceFileName" or die "Can't open $SequenceFileName: $!\n";
 404   for $ID (@{$SequenceDataRef->{IDs}}) {
 405     $Description = $SequenceDataRef->{Description}{$ID};
 406     $Sequence = $SequenceDataRef->{Sequence}{$ID};
 407     $WrappedSequence = WrapText($Sequence, $MaxLength, "\n");
 408 
 409     # Description also contains ID...
 410     print SEQUENCEFILE ">$Description\n";
 411     print SEQUENCEFILE "$WrappedSequence\n";
 412   }
 413   close SEQUENCEFILE;
 414 }
 415 
 416 # Get ID, Sequence and Length for smallest or longest sequence
 417 sub _GetShortestOrLongestSequence {
 418   my($SequencesDataRef, $SequenceType, $IgnoreGaps) = @_;
 419   my($ID, $Seq, $SeqLen, $Description, $FirstID, $FirstSeqLen, $CurrentID, $CurrentSeq, $CurrentSeqLen, $CurrentDescription);
 420 
 421   ($ID, $Seq, $SeqLen) = ('', '', '');
 422   $FirstID = '';
 423 
 424   ID: for $CurrentID (@{$SequencesDataRef->{IDs}}) {
 425     $CurrentSeq = $IgnoreGaps ? RemoveSequenceGaps($SequencesDataRef->{Sequence}{$CurrentID}) : $SequencesDataRef->{Sequence}{$CurrentID};
 426     $CurrentSeqLen = GetSequenceLength($CurrentSeq, $IgnoreGaps);
 427     $CurrentDescription = $SequencesDataRef->{Description}{$CurrentID};
 428     if (!$FirstID) {
 429       $FirstID = $ID; $FirstSeqLen = $CurrentSeqLen;
 430       ($ID, $Seq, $SeqLen, $Description) = ($CurrentID, $CurrentSeq, $CurrentSeqLen, $CurrentDescription);
 431       next ID;
 432     }
 433     if ($CurrentSeqLen != $SeqLen) {
 434       if (($SequenceType =~ /Shortest/i) && ($CurrentSeqLen < $SeqLen)) {
 435         ($ID, $Seq, $SeqLen, $Description) = ($CurrentID, $CurrentSeq, $CurrentSeqLen, $CurrentDescription);
 436       }
 437       elsif (($SequenceType =~ /Longest/i) && ($CurrentSeqLen > $SeqLen) ) {
 438         ($ID, $Seq, $SeqLen, $Description) = ($CurrentID, $CurrentSeq, $CurrentSeqLen, $CurrentDescription);
 439       }
 440     }
 441   }
 442   return ($ID, $Seq, $SeqLen, $Description);
 443 }
 444 
 445 # Remove gaps in the sequence and return new sequence...
 446 sub RemoveSequenceGaps {
 447   my($Seq) = @_;
 448   my($SeqWithoutGaps, $SeqLen, $Index, $Residue);
 449 
 450   $SeqWithoutGaps = '';
 451   $SeqLen = length($Seq);
 452   for $Index (0 .. ($SeqLen - 1)) {
 453     $Residue = substr($Seq, $Index, 1);
 454     if ($Residue =~ /[A-Z]/i) {
 455       $SeqWithoutGaps .= $Residue;
 456     }
 457   }
 458 
 459   return $SeqWithoutGaps;
 460 }
 461 
 462 # Using input alignment data map ref containing following keys, generate
 463 # a new hash with same set of keys after residue columns containg only
 464 # gaps have been removed:
 465 #
 466 # {IDs} : Array of IDs in order as they appear in file
 467 # {Count}: ID count...
 468 # {Description}{$ID} : Description data...
 469 # {Sequence}{$ID} : Sequence data...
 470 #
 471 sub RemoveSequenceAlignmentGapColumns {
 472   my($ID, $AlignmentDataMapRef, %NewAlignmentDataMap);
 473 
 474   ($AlignmentDataMapRef) = @_;
 475 
 476   %NewAlignmentDataMap = ();
 477   @{$NewAlignmentDataMap{IDs}} =();
 478   %{$NewAlignmentDataMap{Description}} =();
 479   %{$NewAlignmentDataMap{Sequence}} =();
 480   $NewAlignmentDataMap{Count} = 0;
 481 
 482   # Transfer ID and count information...
 483   for $ID (@{$AlignmentDataMapRef->{IDs}}) {
 484     push @{$NewAlignmentDataMap{IDs}}, $ID;
 485     $NewAlignmentDataMap{Description}{$ID} = $AlignmentDataMapRef->{Description}{$ID};
 486     $NewAlignmentDataMap{Sequence}{$ID} = '';
 487     $NewAlignmentDataMap{Count} += 1;
 488   }
 489 
 490   # Go over residue columns and transfer the data...
 491   my($FirstID, $FirstSeq, $FirstSeqLen, $Index, $Res, $GapColumn);
 492 
 493   $FirstID = $AlignmentDataMapRef->{IDs}[0];
 494   $FirstSeq = $AlignmentDataMapRef->{Sequence}{$FirstID};
 495   $FirstSeqLen = length($FirstSeq);
 496 
 497   RES: for $Index (0 .. ($FirstSeqLen - 1)) {
 498     # Is this a gap column?
 499     $GapColumn = 1;
 500     ID: for $ID (@{$AlignmentDataMapRef->{IDs}}) {
 501       $Res = substr($AlignmentDataMapRef->{Sequence}{$ID}, $Index, 1);
 502       if ($Res =~ /[A-Z]/i) {
 503         $GapColumn = 0;
 504         last ID;
 505       }
 506     }
 507     if ($GapColumn) {
 508       next RES;
 509     }
 510     # Transfer this residue...
 511     for $ID (@{$AlignmentDataMapRef->{IDs}}) {
 512       $Res = substr($AlignmentDataMapRef->{Sequence}{$ID}, $Index, 1);
 513       $NewAlignmentDataMap{Sequence}{$ID} .= $Res;
 514     }
 515   }
 516 
 517   return (\%NewAlignmentDataMap);
 518 }
 519 
 520 #
 521 # Read sequences file and return a reference to hash with the following keys:
 522 #
 523 # {IDs} - Array of sequence IDs
 524 # {Count} - Number of sequences
 525 # {Description}{$ID} - Sequence description
 526 # {Sequence}{$ID} - Sequence for a specific ID
 527 # {InputFileType} - Sequence file format
 528 # {ConservedAnnotation} - Conserved residue annonation
 529 #
 530 # Note:
 531 #   . Conserved residue annotation either exist in the input sequence alignment file or set
 532 #     for a file containing same number of residues for all the sequence using the following
 533 #     notation: * - Residue conserved; ' ' - Residue not conserved.
 534 #
 535 sub _ReadFileAndSetupSequencesData {
 536   my($SequenceFile, $SequenceType) = @_;
 537   my($SequenceDataMapRef);
 538 
 539   $SequenceDataMapRef = undef;
 540 
 541   # Read sequence file...
 542   $SequenceDataMapRef = '';
 543   if ($SequenceType =~ /^ClustalW$/i) {
 544     $SequenceDataMapRef = _ReadClustalWFile($SequenceFile);
 545   }
 546   elsif ($SequenceType =~ /^Pearson$/i) {
 547     $SequenceDataMapRef = _ReadPearsonFastaFile($SequenceFile);
 548   }
 549   elsif ($SequenceType =~ /^PIR$/i) {
 550     $SequenceDataMapRef = _ReadPIRFastaFile($SequenceFile);
 551   }
 552   elsif ($SequenceType =~ /^MSF$/i) {
 553     $SequenceDataMapRef = _ReadMSFFile($SequenceFile);
 554   }
 555   else {
 556     return $SequenceDataMapRef;
 557   }
 558 
 559   if (exists $SequenceDataMapRef->{ConservedAnnotation}) {
 560     return ($SequenceDataMapRef);
 561   }
 562   if (!(($SequenceDataMapRef->{Count} > 1) && (AreSequenceLengthsIdentical($SequenceDataMapRef)))) {
 563     return ($SequenceDataMapRef);
 564   }
 565 
 566   # Use the first sequence to setup an empty ConservedAnnotation key...
 567   # And mark fully conserved residues...
 568   #
 569   my($ID, $Sequence, $FirstSequence, $FirstSeqLen, $Res, $FirstRes, $ResConserved, $Index);
 570   $ID = $SequenceDataMapRef->{IDs}[0];
 571   $FirstSequence = $SequenceDataMapRef->{Sequence}{$ID};
 572   $FirstSeqLen = length($FirstSequence);
 573   $SequenceDataMapRef->{ConservedAnnotation} = '';
 574   for $Index (0 .. ($FirstSeqLen - 1)) {
 575     $FirstRes = '';
 576     $ResConserved = 1;
 577     ID: for $ID (@{$SequenceDataMapRef->{IDs}}) {
 578       $Sequence = $SequenceDataMapRef->{Sequence}{$ID};
 579       $Res = substr($Sequence, $Index, 1);
 580       if (!$FirstRes) {
 581         $FirstRes = $Res;
 582         next ID;
 583       }
 584       if (($Res !~ /[A-Z]/i) || ($Res ne $FirstRes)) {
 585         $ResConserved = 0;
 586         last ID;
 587       }
 588     }
 589     if ($ResConserved) {
 590       $SequenceDataMapRef->{ConservedAnnotation} .= '*';
 591     }
 592     else {
 593       $SequenceDataMapRef->{ConservedAnnotation} .= ' ';
 594     }
 595   }
 596 
 597   return ($SequenceDataMapRef);
 598 }
 599 
 600 # Read sequence data in ClustalW multiple sequence alignment file and
 601 # return a reference to hash with these keys and values:
 602 #
 603 # {IDs} - Array of sequence IDs
 604 # {Count} - Number of sequences
 605 # {Description}{$ID} - Sequence description
 606 # {Sequence}{$ID} - Sequence for a specific ID
 607 # {InputFileType} - Sequence file format
 608 # {ConservedAnnotation} - Conserved residue annonations: space, *, : , .
 609 #
 610 #
 611 #
 612 # And based on ClustalW/X manual, here is what the ConservedAnnonations mean:
 613 #
 614 # '*' indicates positions which have a single, fully conserved residue
 615 #
 616 # ':' indicates that one of the following 'strong' groups is fully conserved: STA
 617 #    NEQK NHQK NDEQ QHRK MILV MILF HY FYW
 618 
 619 # '.' indicates that one of the following 'weaker' groups is fully conserved:
 620 #     CSA ATV SAG STNK STPA SGND SNDEQK NDEQHK NEQHRK FVLIM HFY
 621 #
 622 # These are all the positively scoring groups that occur in the Gonnet Pam250
 623 # matrix. The strong and weak groups are defined as strong score >0.5 and weak
 624 # score =<0.5 respectively.
 625 #
 626 sub _ReadClustalWFile {
 627   my($SequenceFile) = @_;
 628   my(%SequencesDataMap);
 629 
 630   # Initialize data...
 631   %SequencesDataMap = ();
 632   @{$SequencesDataMap{IDs}} = ();
 633   %{$SequencesDataMap{Description}} = ();
 634   %{$SequencesDataMap{Sequence}} = ();
 635   $SequencesDataMap{Count} = 0;
 636   $SequencesDataMap{ConservedAnnotation} = '';
 637   $SequencesDataMap{InputFileType} = 'ClustalW';
 638 
 639   open SEQUENCEFILE, "$SequenceFile" or die "Couldn't open $SequenceFile: $!\n";
 640 
 641   my($Line, $LineLength, $AnnotationStart, $AnnotationLength, $Annotation, $Sequence, $SequenceLength, $ID, $IDIndex);
 642 
 643   # Ignore the header line...
 644   $Line = <SEQUENCEFILE>;
 645 
 646   LINE: while ($Line = GetTextLine(\*SEQUENCEFILE)) {
 647     if (($Line =~ /^[ \*\:\.]/) && ($Line !~ /[A-Z]/i)) {
 648       # Annotation for sequences: fully conserverd, weaker or stronger group conserverd.
 649       # Extract it and save...
 650       $LineLength = length($Line);
 651       $AnnotationStart = $LineLength - $SequenceLength;
 652       $AnnotationLength = $SequenceLength;
 653       $Annotation = substr($Line, $AnnotationStart, $AnnotationLength);
 654       $SequencesDataMap{ConservedAnnotation} .= $Annotation;
 655     }
 656     else {
 657       # Extract ID and sequences...
 658       ($ID, $Sequence)= $Line =~ /^[ ]*(.*?)[ ]+(.*?)[ 01-9]*$/;
 659       $Sequence =~ s/ //g;
 660       if (!($ID && $Sequence)) {
 661         next LINE;
 662       }
 663 
 664       if (exists $SequencesDataMap{Sequence}{$ID}) {
 665         # Append to existing alignment value...
 666         $SequenceLength = length($Sequence);
 667         $SequencesDataMap{Sequence}{$ID} .= $Sequence;
 668       }
 669       else {
 670         # New alignment data...
 671         $SequencesDataMap{Count} += 1;
 672         push @{$SequencesDataMap{IDs}}, $ID;
 673         $SequencesDataMap{Description}{$ID} = $ID;
 674         $SequencesDataMap{Sequence}{$ID} = $Sequence;
 675         $SequenceLength = length($Sequence);
 676       }
 677     }
 678   }
 679   close SEQUENCEFILE;
 680   return (\%SequencesDataMap);
 681 }
 682 
 683 # Read Pearson fasta file and return a reference to hash with these keys:
 684 #
 685 # {IDs} - Array of sequence IDs
 686 # {Count} - Number of sequences
 687 # {Description}{$ID} - Sequence description
 688 # {Sequence}{$ID} - Sequence for a specific ID
 689 # {InputFileType} - Sequence file format
 690 # {ConservedAnnotation} - Conserved residue annonation
 691 #
 692 sub _ReadPearsonFastaFile {
 693   my($FastaFileName, $ID, $Description, $Line, $IgnoreID, @LineWords, %FastaDataMap);
 694 
 695   ($FastaFileName) = @_;
 696 
 697   %FastaDataMap = ();
 698   @{$FastaDataMap{IDs}} =();
 699   %{$FastaDataMap{Description}} =();
 700   %{$FastaDataMap{Sequence}} =();
 701   $FastaDataMap{Count} = 0;
 702   $FastaDataMap{InputFileType} = 'Pearson';
 703 
 704   open FASTAFILE, "$FastaFileName" or die "Couldn't open $FastaFileName: $!\n";
 705   $ID = '';
 706   $IgnoreID = 0;
 707   LINE: while ($Line = GetTextLine(\*FASTAFILE)) {
 708     if ($Line =~ /^\>/) {
 709       # Start of a new ID...
 710       $Line =~ s/^\>//;
 711       $Line = RemoveLeadingWhiteSpaces($Line);
 712       @LineWords = ();
 713       @LineWords = split / /, $Line;
 714 
 715       $ID = $LineWords[0];
 716       $ID =~ s/ //g;
 717       $Description = $Line;
 718 
 719       $IgnoreID = 0;
 720       if (exists $FastaDataMap{Sequence}{$ID}) {
 721         $IgnoreID = 1;
 722         warn "Warning: ID, $ID, in Fasta file already exists. Ignoring ID and sequence data...\n";
 723         next LINE;
 724       }
 725       push @{$FastaDataMap{IDs}}, $ID;
 726       $FastaDataMap{Description}{$ID} = $Description;
 727       $FastaDataMap{Count} += 1;
 728       next LINE;
 729     }
 730     if ($IgnoreID) { next LINE; }
 731 
 732     # Remove any spaces in the sequence...
 733     $Line =~ s/ //g;
 734     # Sequence data for active ID...
 735     if (exists $FastaDataMap{Sequence}{$ID}) {
 736       $FastaDataMap{Sequence}{$ID} .= $Line;
 737     }
 738     else {
 739       $FastaDataMap{Sequence}{$ID} = $Line;
 740     }
 741   }
 742   close FASTAFILE;
 743   return \%FastaDataMap;
 744 }
 745 
 746 # Read PIR fasta file and return a reference to hash with these keys:
 747 #
 748 # {IDs} - Array of sequence IDs
 749 # {Count} - Number of sequences
 750 # {Description}{$ID} - Sequence description
 751 # {Sequence}{$ID} - Sequence for a specific ID
 752 # {InputFileType} - Sequence file format
 753 # {ConservedAnnotation} - Conserved residue annonation
 754 #
 755 # Format:
 756 # A sequence in PIR format consists of:
 757 # One line starting with
 758 #   a ">" (greater-than) sign, followed by
 759 #   a two-letter code describing the sequence type code (P1, F1, DL, DC, RL, RC, N3, N1 or XX), followed by
 760 #   a semicolon, followed by
 761 #   the sequence identification code (the database ID-code).
 762 # One line containing a textual description of the sequence.
 763 # One or more lines containing the sequence itself. The end of the
 764 # sequence is marked by a "*" (asterisk) character.
 765 #
 766 # A file in PIR format may comprise more than one sequence.
 767 #
 768 # The PIR format is also often referred to as the NBRF format.
 769 #
 770 # Code SequenceType
 771 # P1    Protein (complete)
 772 # F1    Protein (fragment)
 773 # DL    DNA (linear)
 774 # DC    DNA (circular)
 775 # RL    RNA (linear)
 776 # RC   RNA (circular)
 777 # N3    tRNA
 778 # N1    Other functional RNA
 779 #
 780 
 781 sub _ReadPIRFastaFile {
 782   my($FastaFileName, $ID, $Description, $Line, $SequenceTypeCode, $ReadingSequenceData, %FastaDataMap);
 783 
 784   ($FastaFileName) = @_;
 785 
 786   %FastaDataMap = ();
 787   @{$FastaDataMap{IDs}} =();
 788   %{$FastaDataMap{Description}} =();
 789   %{$FastaDataMap{Sequence}} =();
 790   %{$FastaDataMap{SequenceTypeCode}} =();
 791   $FastaDataMap{Count} = 0;
 792   $FastaDataMap{InputFileType} = 'PIR';
 793 
 794   open FASTAFILE, "$FastaFileName" or die "Couldn't open $FastaFileName: $!\n";
 795   $ID = '';
 796   $ReadingSequenceData = 0;
 797   LINE: while ($Line = GetTextLine(\*FASTAFILE)) {
 798     if ($Line =~ /^\>/) {
 799       # Start of a new ID...
 800       $Line =~ s/^\>//;
 801       $Line = RemoveLeadingWhiteSpaces($Line);
 802       ($SequenceTypeCode, $ID) = /^\>(.*?)\;(.*?)$/;
 803 
 804       # Use next line to retrieve sequence description...
 805       $Line = GetTextLine(\*FASTAFILE);
 806       $Line = RemoveLeadingWhiteSpaces($Line);
 807       $Description = $Line;
 808 
 809       if (exists $FastaDataMap{Sequence}{$ID}) {
 810         warn "Warning: ID, $ID, in Fasta file already exists. Ignoring ID and sequence data...\n";
 811         next LINE;
 812       }
 813       $ReadingSequenceData = 1;
 814       push @{$FastaDataMap{IDs}}, $ID;
 815       $FastaDataMap{SequenceTypeCode}{$ID} = $SequenceTypeCode;
 816       $FastaDataMap{Description}{$ID} = $Description;
 817       $FastaDataMap{Count} += 1;
 818       next LINE;
 819     }
 820     if (!$ReadingSequenceData) { next LINE; }
 821 
 822     # Remove any spaces in the sequence...
 823     $Line =~ s/ //g;
 824     if ($Line =~ /[\*]$/) {
 825       # End of sequence...
 826       $ReadingSequenceData = 0;
 827       $Line =~ s/[\*]$//;
 828     }
 829     # Sequence data for active ID...
 830     if (exists $FastaDataMap{Sequence}{$ID}) {
 831       $FastaDataMap{Sequence}{$ID} .= $Line;
 832     }
 833     else {
 834       $FastaDataMap{Sequence}{$ID} = $Line;
 835     }
 836   }
 837   close FASTAFILE;
 838   return \%FastaDataMap;
 839 }
 840 
 841 # Read MSF file and return a reference to hash with these keys:
 842 #
 843 # {IDs} : Array of IDs in order as they appear in file
 844 # {Count}: ID count...
 845 # {Description}{$ID} : Description data...
 846 # {Sequence}{$ID} : Sequence data...
 847 #
 848 sub _ReadMSFFile {
 849   my($MSFFileName, $Line, @LineWords, %MSFDataMap);
 850 
 851   ($MSFFileName) = @_;
 852 
 853   %MSFDataMap = ();
 854   @{$MSFDataMap{IDs}} =();
 855   %{$MSFDataMap{Description}} =();
 856   %{$MSFDataMap{Sequence}} =();
 857   $MSFDataMap{Count} = 0;
 858   $MSFDataMap{InputFileType} = 'MSF';
 859 
 860   open MSFFILE, "$MSFFileName" or die "Couldn't open $MSFFileName: $!\n";
 861 
 862   # Collect sequences and IDs...
 863   #
 864   # '//' after the name fields indicates end of header list and start of sequence data.
 865   #
 866   my($ID, $Len, $Check, $Weight, $Sequence, $NameFieldsFound, %MSFIDsMap);
 867   %MSFIDsMap = ();
 868   $NameFieldsFound = 0;
 869   LINE: while ($Line = GetTextLine(\*MSFFILE)) {
 870     if ($Line =~ /Name:/) {
 871       $NameFieldsFound++;
 872       ($ID, $Len, $Check, $Weight) = $Line =~ /^[ ]*Name:[ ]+(.*?)[ ]+Len:[ ]+(.*?)[ ]+Check:[ ]+(.*?)[ ]+Weight:[ ]+(.*?)[ ]*$/;
 873       if ($ID =~ / /) {
 874         ($ID) = $ID =~ /^(.*?)[ ]+/
 875       }
 876       if (exists $MSFIDsMap{$ID}) {
 877         warn "Warning: ID, $ID, in MSF file already exists. Ignoring ID and sequence data...\n";
 878         next LINE;
 879       }
 880       $MSFIDsMap{$ID} = $ID;
 881       push @{$MSFDataMap{IDs}}, $ID;
 882       $MSFDataMap{Description}{$ID} = $ID;
 883       $MSFDataMap{Count} += 1;
 884     }
 885     elsif ( /\/\// && $NameFieldsFound) {
 886       # End of header list...
 887       last LINE;
 888     }
 889   }
 890   # Collect all sequences...
 891   #
 892   my($FirstField, $SecondField);
 893   while ($Line = GetTextLine(\*MSFFILE)) {
 894     ($FirstField, $SecondField) = $Line =~ /^[ ]*(.*?)[ ]+(.*?)$/;
 895     if (exists $MSFIDsMap{$FirstField}) {
 896       # It's ID and sequence data...
 897       $ID = $FirstField;
 898       $Sequence = $SecondField;
 899       # Take out spaces and leave the gap characters...
 900       $Sequence =~ s/ //g;
 901       if ($MSFDataMap{Sequence}{$ID}) {
 902         $MSFDataMap{Sequence}{$ID} .= $Sequence;
 903       }
 904       else {
 905         $MSFDataMap{Sequence}{$ID} = $Sequence;
 906       }
 907     }
 908   }
 909 
 910   close MSFFILE;
 911   return \%MSFDataMap;
 912 }
 913 
 914