forester/perl/forester.pm

   1 # $Id: forester.pm,v 1.26 2010/12/13 19:00:22 cmzmasek Exp $
   2 #
   3 # FORESTER -- software libraries and applications
   4 # for evolutionary biology research and applications.
   5 #
   6 # Copyright (C) 2007-2009 Christian M. Zmasek
   7 # Copyright (C) 2007-2009 Burnham Institute for Medical Research
   8 # All rights reserved
   9 #
  10 # This library is free software; you can redistribute it and/or
  11 # modify it under the terms of the GNU Lesser General Public
  12 # License as published by the Free Software Foundation; either
  13 # version 2.1 of the License, or (at your option) any later version.
  14 #
  15 # This library is distributed in the hope that it will be useful,
  16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18 # Lesser General Public License for more details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public
  21 # License along with this library; if not, write to the Free Software
  22 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  23 #
  24 # Contact: phylosoft @ gmail . com
  25 #     WWW: www.phylosoft.org/forester
  26 #
  27 #
  28 #
  29
  30
  31 package forester;
  32 use strict;
  33 require Exporter;
  34
  35 our $VERSION = 1.000;
  36
  37 our @ISA    = qw( Exporter );
  38
  39 our @EXPORT = qw( executeConsense
  40                   executePhyloPl
  41                   executePuzzleDQO
  42                   executePuzzleDQObootstrapped
  43                   pfam2phylipMatchOnly
  44                   startsWithSWISS_PROTname
  45                   isPfamSequenceLine
  46                   isPfamCommentLine
  47                   containsPfamNamedSequence
  48                   isRFline
  49                   executeProtpars
  50                   setModelForPuzzle
  51                   setRateHeterogeneityOptionForPuzzle
  52                   setParameterEstimatesOptionForPuzzle
  53                   executePuzzleBootstrapped
  54                   executePuzzle
  55                   executeFastme
  56                   executeNeighbor
  57                   executeFitch
  58                   executeBionj
  59                   executeWeighbor
  60                   executePhyml
  61                   executeHmmfetch
  62                   addDistsToQueryToPWDfile
  63                   testForTextFilePresence
  64                   exitWithWarning
  65                   dieWithUnexpectedError
  66                   addSlashAtEndIfNotPresent
  67                   $LENGTH_OF_NAME
  68                   $MIN_NUMBER_OF_AA
  69                   $TREMBL_ACDEOS_FILE
  70                   $SWISSPROT_ACDEOS_FILE
  71                   $SPECIES_NAMES_FILE
  72                   $SPECIES_TREE_FILE_DEFAULT
  73                   $MULTIPLE_TREES_FILE_SUFFIX
  74                   $LOG_FILE_SUFFIX
  75                   $ALIGN_FILE_SUFFIX
  76                   $TREE_FILE_SUFFIX
  77                   $ADDITION_FOR_RIO_ANNOT_TREE
  78                   $SUFFIX_PWD
  79                   $SUFFIX_BOOT_STRP_POS
  80                   $MULTIPLE_PWD_FILE_SUFFIX
  81                   $SUFFIX_PWD_NOT_BOOTS
  82                   $SUFFIX_HMM
  83                   $MATRIX_FOR_PWD
  84                   $RIO_PWD_DIRECTORY
  85                   $RIO_BSP_DIRECTORY
  86                   $RIO_NBD_DIRECTORY
  87                   $RIO_ALN_DIRECTORY
  88                   $RIO_HMM_DIRECTORY
  89                   $PFAM_FULL_DIRECTORY
  90                   $PFAM_SEED_DIRECTORY
  91                   $PRIOR_FILE_DIR
  92                   $PFAM_HMM_DB
  93                   $FORESTER_JAR
  94                   $SEQBOOT
  95                   $NEIGHBOR
  96                   $PROTPARS
  97                   $CONSENSE
  98                   $PROML
  99                   $PHYLIP_VERSION
 100                   $PUZZLE
 101                   $PUZZLE_VERSION
 102                   $FASTME
 103                   $FASTME_VERSION
 104                   $BIONJ
 105                   $BIONJ_VERSION
 106                   $WEIGHBOR
 107                   $WEIGHBOR_VERSION
 108                   $RAXML
 109                   $RAXML_VERSION
 110                   $PHYML
 111                   $PHYML_VERSION
 112                   $HMMALIGN
 113                   $HMMSEARCH
 114                   $HMMBUILD
 115                   $HMMFETCH
 116                   $SFE
 117                   $HMMCALIBRATE
 118                   $P7EXTRACT
 119                   $MULTIFETCH
 120                   $BOOTSTRAP_CZ
 121                   $BOOTSTRAP_CZ_PL
 122                   $SUPPORT_TRANSFER
 123                   $SUPPORT_STATISTICS
 124                   $NEWICK_TO_PHYLOXML
 125                   $PHYLO_PL
 126                   $RIO_PL
 127                   $DORIO
 128                   $PUZZLE_DQO
 129                   $BOOTSTRAPS
 130                   $PATH_TO_FORESTER
 131                   $JAVA
 132                   $NODE_LIST
 133                   $RIO_SLAVE_DRIVER
 134                   $RIO_SLAVE
 135                   $TEMP_DIR_DEFAULT
 136                   $EXPASY_SPROT_SEARCH_DE
 137                   $EXPASY_SPROT_SEARCH_AC
 138  );
 139
 140
 141
 142
 143 # =============================================================================
 144 # =============================================================================
 145 #
 146 # THESE VARIABLES ARE ENVIRONMENT DEPENDENT, AND NEED TO BE SET ACCORDINGLY
 147 # BY THE USER
 148 # -------------------------------------------------------------------------
 149 #
 150
 151 # For using just "phylo_pl.pl", only the following variables need to be set
 152 # $JAVA
 153 # $FORESTER_JAR
 154 # $TEMP_DIR_DEFAULT
 155 # $SEQBOOT
 156 # $CONSENSE
 157 # $PUZZLE
 158 # $FASTME
 159 # $NEIGHBOR
 160 # $FITCH
 161 # $BIONJ
 162 # $WEIGHBOR
 163 # $PHYML
 164 # $PROTPARS
 165
 166 # Software directory:
 167 # ---------------------
 168
 169 our $SOFTWARE_DIR              = "/home/zma/SOFTWARE/";
 170
 171
 172 # Java virtual machine:
 173 # ---------------------
 174 our $JAVA                      = "java";
 175
 176
 177 # Where all the temporary files can be created:
 178 # ---------------------------------------------
 179 our $TEMP_DIR_DEFAULT          = "/tmp/";
 180
 181
 182 # Programs from Joe Felsenstein's PHYLIP package:
 183 # -----------------------------------------------
 184 our $SEQBOOT                   = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/seqboot";
 185 our $NEIGHBOR                  = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/neighbor";
 186 our $PROTPARS                  = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/protpars";
 187 our $PROML                     = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/proml";
 188 our $FITCH                     = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/fitch";
 189 our $CONSENSE                  = $SOFTWARE_DIR."PHYLO/Phylip/Phylip3.695/phylip-3.696/exe/consense";
 190 our $PHYLIP_VERSION            = "3.695";
 191
 192 # TREE-PUZZLE:
 193 # ------------
 194 our $PUZZLE                    = $SOFTWARE_DIR."PHYLO/TREE-PUZZLE/tree-puzzle-5.2/src/puzzle";
 195 our $PUZZLE_VERSION            = "5.2";
 196
 197 # FASTME:
 198 # -----------------------------------------------------
 199 our $FASTME                    = $SOFTWARE_DIR."PHYLO/FastME/fastme2.0/fastme";
 200 our $FASTME_VERSION            = "2.0";
 201
 202 # BIONJ:
 203 # -----------------------------------------------------
 204 our $BIONJ                     = "";
 205 our $BIONJ_VERSION             = "";
 206
 207 # WEIGHBOR:
 208 # -----------------------------------------------------
 209 our $WEIGHBOR                  = "";
 210 our $WEIGHBOR_VERSION          = "";
 211
 212 # PHYML:
 213 # -----------------------------------------------------
 214 our $PHYML                     = $SOFTWARE_DIR."PHYLO/PhyML/PhyML-3.1/PhyML-3.1/PhyML-3.1_linux64";
 215 our $PHYML_VERSION             = "3.1";
 216
 217 # RAXML:
 218 # -----------------------------------------------------
 219 our $RAXML                     = $SOFTWARE_DIR."PHYLO/RAxML/20161215/standard-RAxML-master/raxmlHPC-AVX";
 220 our $RAXML_VERSION             = "8.2.9";
 221
 222
 223 # forester.jar. This jar file is currently available at: https://sites.google.com/site/cmzmasek/home/software/forester
 224 # --------------------------------------------------------------------------------------------------------------------
 225
 226 our $FORESTER_JAR              = "/home/zma/git/forester/forester/java/forester.jar";
 227
 228
 229
 230 # End of variables which need to be set by the user for using "phylo_pl.pl".
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245 # Tool from forester.jar to transfer support values:
 246 # -------------------------------------------------
 247 our $SUPPORT_TRANSFER            = $JAVA." -cp $FORESTER_JAR org.forester.application.support_transfer";
 248
 249
 250
 251 # Tool from forester.jar for simple statistics for support values:
 252 # ----------------------------------------------------------------
 253 our $SUPPORT_STATISTICS          = $JAVA." -cp $FORESTER_JAR org.forester.application.support_statistics";
 254
 255
 256 # Tool from forester.jar to transfer nh to phyloXML:
 257 # -------------------------------------------------
 258 our $NEWICK_TO_PHYLOXML          = $JAVA." -cp $FORESTER_JAR org.forester.application.phyloxml_converter";
 259
 260
 261
 262 # FORESTER itself (currently not needed for "phylo_pl.pl"):
 263 # ---------------------------------------------------------
 264 our $PATH_TO_FORESTER          = "";
 265
 266
 267 # Pfam data (not needed for phylo_pl.pl):
 268 # --------------------------------------
 269 our $PFAM_FULL_DIRECTORY       = "/path/to/Pfam/Full/";
 270 our $PFAM_SEED_DIRECTORY       = "/path/to/Pfam/Seed/";
 271 our $PFAM_HMM_DB               = "/path/to/Pfam/Pfam_ls"; # Need to run "hmmindex" on this
 272                                                           # to produce .ssi file.
 273                                                           # Then, for example
 274                                                           # "setenv HMMERDB /home/rio/pfam-6.6/"
 275
 276
 277 $PATH_TO_FORESTER = &addSlashAtEndIfNotPresent( $PATH_TO_FORESTER );
 278
 279
 280 # Description lines and species from SWISS-PROT and TrEMBL (not needed for phylo_pl.pl):
 281 # -------------------------------------------------------------------------------------
 282 our $TREMBL_ACDEOS_FILE        = $PATH_TO_FORESTER."data/trembl22_ACDEOS_1-6";
 283
 284 our $SWISSPROT_ACDEOS_FILE     = $PATH_TO_FORESTER."data/sp40_ACDEOS_1-6";
 285
 286
 287
 288 # Names of species which can be analyzed and analyzed
 289 # against (must also be in tree $SPECIES_TREE_FILE_DEFAULT).
 290 # By using a list with less species, RIO analyses become faster
 291 # but lose phylogenetic resolution.
 292 # For many purposes, list "tree_of_life_bin_1-6_species_list"
 293 # in "data/species/" might be sufficient:
 294 # (not needed for phylo_pl.pl)
 295 # --------------------------------------------------------------
 296 our $SPECIES_NAMES_FILE        = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6_species_list";
 297
 298
 299
 300 # A default species tree in NHX format.
 301 # For many purposes, tree "tree_of_life_bin_1-6.nhx"
 302 # in "data/species/" might be fine:
 303 # (not needed for phylo_pl.pl)
 304 # --------------------------------------------------
 305 our $SPECIES_TREE_FILE_DEFAULT = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6.nhx";
 306
 307
 308
 309 # Data for using precalculated distances:
 310 # (not needed for phylo_pl.pl)
 311 # ---------------------------------------
 312 our $MATRIX_FOR_PWD            = 2;  # The matrix which has been used for the pwd in $RIO_PWD_DIRECTORY.
 313                                      # 0=JTT, 1=PAM, 2=BLOSUM 62, 3=mtREV24, 5=VT, 6=WAG.
 314
 315 our $RIO_PWD_DIRECTORY         = $PATH_TO_FORESTER."example_data/";  # all must end with "/"
 316 our $RIO_BSP_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 317 our $RIO_NBD_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 318 our $RIO_ALN_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 319 our $RIO_HMM_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 320
 321
 322
 323 #
 324 # End of variables which need to be set by the user.
 325 #
 326 # =============================================================================
 327 # =============================================================================
 328
 329
 330
 331
 332
 333 $TEMP_DIR_DEFAULT    = &addSlashAtEndIfNotPresent( $TEMP_DIR_DEFAULT );
 334 $PFAM_FULL_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_FULL_DIRECTORY );
 335 $PFAM_SEED_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_SEED_DIRECTORY );
 336
 337
 338
 339 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 340 # These variables should normally not be changed:
 341 #
 342
 343 our $PRIOR_FILE_DIR            = $PATH_TO_FORESTER."data/priors_for_hmmbuild/";
 344                                  # Directory containing dirichlet prior
 345                                  # files needed for certain aligments
 346                                  # by hmmbuild (e.g. Collagen).
 347
 348
 349
 350
 351
 352 # TREE-PUZZLE:
 353 our $PUZZLE_DQO                = $PATH_TO_FORESTER."puzzle_dqo/src/puzzle";
 354
 355 # HMMER:
 356 our $HMMALIGN                  = $PATH_TO_FORESTER."hmmer/binaries/hmmalign";
 357 our $HMMSEARCH                 = $PATH_TO_FORESTER."hmmer/binaries/hmmsearch";
 358 our $HMMBUILD                  = $PATH_TO_FORESTER."hmmer/binaries/hmmbuild";
 359 our $HMMFETCH                  = $PATH_TO_FORESTER."hmmer/binaries/hmmfetch";
 360 our $SFE                       = $PATH_TO_FORESTER."hmmer/binaries/sfetch";
 361 our $HMMCALIBRATE              = $PATH_TO_FORESTER."hmmer/binaries/hmmcalibrate";
 362
 363 our $P7EXTRACT                 = $PATH_TO_FORESTER."perl/p7extract.pl";
 364 our $MULTIFETCH                = $PATH_TO_FORESTER."perl/multifetch.pl";
 365
 366
 367 # RIO/FORESTER:
 368 our $BOOTSTRAP_CZ              = $PATH_TO_FORESTER."C/bootstrap_cz";
 369 our $BOOTSTRAP_CZ_PL           = $PATH_TO_FORESTER."perl/bootstrap_cz.pl";
 370 #our $SUPPORT_TRANSFER         = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.transfersBranchLenghts";
 371 #our $SUPPORT_TRANSFER         = $JAVA." -cp /home/czmasek/SOFTWARE/FORESTER/forester3/forester.jar org.forester.tools.SupportTransfer";
 372
 373 our $PHYLO_PL                  = $PATH_TO_FORESTER."perl/phylo_pl.pl";
 374 our $RIO_PL                    = $PATH_TO_FORESTER."perl/rio.pl";
 375 our $DORIO                     = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.DoRIO";
 376 # parallel RIO:
 377 our $RIO_SLAVE_DRIVER          = $PATH_TO_FORESTER."perl/rio_slave_driver.pl";
 378 our $RIO_SLAVE                 = $PATH_TO_FORESTER."perl/rio_slave.pl";
 379 our $NODE_LIST                 = $PATH_TO_FORESTER."data/node_list.dat";
 380
 381 #
 382 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 383
 384
 385 our $BOOTSTRAPS         = 100;
 386 our $MIN_NUMBER_OF_AA   = 20;  # After removal of gaps, if less, gaps are not removed.
 387 our $LENGTH_OF_NAME     = 10;
 388
 389
 390
 391
 392 our $MULTIPLE_TREES_FILE_SUFFIX  = ".mlt";
 393 our $LOG_FILE_SUFFIX             = ".log";
 394 our $ALIGN_FILE_SUFFIX           = ".aln";
 395 our $TREE_FILE_SUFFIX            = ".nhx";
 396 our $ADDITION_FOR_RIO_ANNOT_TREE = ".rio";
 397 our $SUFFIX_PWD                  = ".pwd";
 398 our $MULTIPLE_PWD_FILE_SUFFIX    = ".mpwd";
 399 our $SUFFIX_BOOT_STRP_POS        = ".bsp";
 400 our $SUFFIX_PWD_NOT_BOOTS        = ".nbd";
 401 our $SUFFIX_HMM                  = ".hmm";
 402
 403 our $EXPASY_SPROT_SEARCH_DE      = "http://www.expasy.org/cgi-bin/sprot-search-de?";
 404 our $EXPASY_SPROT_SEARCH_AC      = "http://www.expasy.org/cgi-bin/sprot-search-ac?";
 405
 406
 407
 408 # One argument: input multiple trees file
 409 # Last modified: 07/05/01
 410 sub executeConsense {
 411     my $in  = $_[ 0 ];
 412
 413     &testForTextFilePresence( $in );
 414
 415     system( "$CONSENSE >/dev/null 2>&1 << !
 416 $in
 417 Y
 418 !" )
 419     && &dieWithUnexpectedError( "Could not execute \"$CONSENSE $in\"" );
 420
 421     return;
 422 }
 423
 424
 425
 426 # Four arguments:
 427 # 1. options ("-" is not necessary)
 428 # 2. alignment or pwd file
 429 # 3. outfile
 430 # 4. temp dir
 431 # Last modified: 07/05/01
 432 sub executePhyloPl {
 433
 434     my $opts         = $_[ 0 ];
 435     my $B            = $_[ 1 ];
 436     my $C            = $_[ 2 ];
 437     my $D            = $_[ 3 ];
 438
 439     &testForTextFilePresence( $B );
 440
 441     $opts = "-".$opts;
 442
 443     system( "$PHYLO_PL $opts $B $C $D" )
 444     && &dieWithUnexpectedError( "Could not execute \"$PHYLO_PL $opts $B $C $D\"" );
 445
 446 } ## executePhyloPl
 447
 448
 449
 450
 451 # Two arguments:
 452 # 1. Name of inputfile
 453 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 454 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 455 sub executePuzzleDQO {
 456     my $in            = $_[ 0 ];
 457     my $matrix_option = $_[ 1 ];
 458     my $mat           = "";
 459
 460     &testForTextFilePresence( $in );
 461
 462     $mat = setModelForPuzzle( $matrix_option );
 463
 464     system( "$PUZZLE_DQO $in >/dev/null 2>&1 << !$mat
 465 y
 466 !" )
 467     && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO\"" );
 468
 469     return;
 470
 471 } ## executePuzzleDQO
 472
 473
 474
 475
 476 # Two arguments:
 477 # 1. Name of inputfile
 478 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 479 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 480 # Last modified: 01/28/02
 481 sub executePuzzleDQObootstrapped {
 482     my $in            = $_[ 0 ];
 483     my $matrix_option = $_[ 1 ];
 484
 485
 486     my $l             = 0;
 487     my $slen          = 0;
 488     my $counter       = 0;
 489     my $mat           = "";
 490     my $a             = "";
 491     my @a             = ();
 492
 493     &testForTextFilePresence( $in );
 494
 495     open( GRP, "<$in" ) || &dieWithUnexpectedError( "Cannot open file \"$in\"" );
 496     while( <GRP> ) {
 497         if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
 498             $counter++;
 499         }
 500     }
 501     close( GRP );
 502
 503     $l   = `cat $in | wc -l`;
 504     $slen   = $l / $counter;
 505
 506     system( "split -$slen $in $in.splt." )
 507     && &dieWithUnexpectedError( "Could not execute \"split -$slen $in $in.splt.\"" );
 508
 509     @a = <$in.splt.*>;
 510
 511     $mat = setModelForPuzzle( $matrix_option );
 512
 513     foreach $a ( @a ) {
 514
 515         system( "$PUZZLE_DQO $a >/dev/null 2>&1 << !$mat
 516 y
 517 !" )
 518         && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO $a\"" );
 519
 520         system( "cat $a.dist >> $in.dist" )
 521         && &dieWithUnexpectedError( "Could not execute \"cat outdist >> $in.dist\"" );
 522
 523         unlink( $a, $a.".dist" );
 524     }
 525
 526     return;
 527
 528 } ## executePuzzleDQObootstrapped
 529
 530
 531
 532 # Transfers a Pfam (SELEX) alignment to a
 533 # PHYLIP sequential style alignment.
 534 # It only writes "match columns" as indicated by the
 535 # "# RF" line ('x' means match).
 536 #
 537 # Three arguments:
 538 # 1. infile name
 539 # 2. outfile name
 540 # 3. 1 to NOT ensure that match states contain only 'A'-'Z' or '-'
 541 #
 542 # Returns the number of match states (=length of output alignment),
 543 #         the length of the input alignment,
 544 #         the number of seqs in the input alignment
 545 #
 546 # Last modified: 07/07/01
 547 #
 548 sub pfam2phylipMatchOnly {
 549
 550     my $infile          = $_[ 0 ];
 551     my $outfile         = $_[ 1 ];
 552     my $ne              = $_[ 2 ];
 553     my @seq_name        = ();
 554     my @seq_array       = ();
 555     my $return_line     = "";
 556     my $seq             = "";
 557     my $x               = 0;
 558     my $y               = 0;
 559     my $i               = 0;
 560     my $x_offset        = 0;
 561     my $max_x           = 0;
 562     my $rf_y            = 0;
 563     my $number_colum    = 0;
 564     my $not_ensure      = 0;
 565     my $saw_rf_line     = 0;
 566
 567     if ( $ne && $ne == 1 ) {
 568         $not_ensure = 1;
 569     }
 570
 571     &testForTextFilePresence( $infile );
 572
 573     open( INPP, "$infile" ) || &dieWithUnexpectedError( "Cannot open file \"$infile\"" );
 574
 575     # This reads in the first block. It reads in the seq names.
 576     while ( 1 ) {
 577         if ( &isPfamSequenceLine( $return_line ) ) {
 578             $return_line =~ /^(\S+)\s+(\S+)/;
 579             $seq_name[ $y ] = substr( $1, 0, $LENGTH_OF_NAME );
 580             $seq = $2;
 581             for ( $x = 0; $x < length( $seq ); $x++ ) {
 582                 $seq_array[ $x ][ $y ] = substr( $seq, $x, 1 );
 583             }
 584             $y++;
 585         }
 586         elsif ( &isRFline( $return_line ) ) {
 587             $saw_rf_line = 1;
 588             $return_line =~ /\s+(\S+)\s*$/;
 589             $seq = $1;
 590             $x_offset = length( $seq );
 591             $rf_y = $y;
 592             for ( $x = 0; $x < $x_offset; $x++ ) {
 593                 $seq_array[ $x ][ $rf_y ] = substr( $seq, $x, 1 );
 594             }
 595             last;
 596         }
 597
 598         $return_line = <INPP>;
 599
 600         if ( !$return_line ) {
 601              &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
 602         }
 603     }
 604
 605     if ( $saw_rf_line != 1 ) {
 606          &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
 607     }
 608
 609     $y = 0;
 610     $max_x = 0;
 611
 612     # This reads all blocks after the 1st one.
 613     while ( $return_line = <INPP> ) {
 614         if ( &isPfamSequenceLine( $return_line ) ) {
 615             $return_line =~ /^\S+\s+(\S+)/;
 616             $seq = $1;
 617             for ( $x = 0; $x < length( $seq ); $x++ ) {
 618                 $seq_array[ $x + $x_offset ][ $y % $rf_y ] = substr( $seq, $x, 1 );
 619             }
 620             $y++;
 621         }
 622         elsif ( &isRFline( $return_line ) ) {
 623             if ( $y != $rf_y ) {
 624                 &dieWithUnexpectedError( "Alignment not in expected format" );
 625             }
 626
 627             $return_line =~ /\s+(\S+)\s*$/;
 628             $seq = $1;
 629             $max_x = length( $seq );
 630
 631             for ( $x = 0; $x < length( $seq ); $x++ ) {
 632                 $seq_array[ $x + $x_offset ][ $rf_y ] = substr( $seq, $x, 1 );
 633             }
 634
 635             $y = 0;
 636             $x_offset = $x_offset + $max_x;
 637             $max_x = 0;
 638         }
 639     }
 640
 641     close( INPP );
 642
 643     # Counts the match states, and hence the number of aa in the alignment:
 644     for ( $x = 0; $x < $x_offset; $x++ ) {
 645         if ( !$seq_array[ $x ][ $rf_y ] ) {
 646             &dieWithUnexpectedError( "Alignment not in expected format" );
 647         }
 648         if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
 649             $number_colum++;
 650         }
 651     }
 652
 653     # Writes the file:
 654
 655     open( OUTPP, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
 656     print OUTPP "$rf_y $number_colum\n";
 657     for ( $y = 0; $y < $rf_y; $y++ ) {
 658         print OUTPP "$seq_name[ $y ]";
 659         for ( $i = 0; $i < ( $LENGTH_OF_NAME - length( $seq_name[ $y ] ) ); $i++ ) {
 660             print OUTPP " ";
 661         }
 662         for ( $x = 0; $x < $x_offset; $x++ ) {
 663             if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
 664                 if ( !$seq_array[ $x ][ $y ] ) {
 665                     &dieWithUnexpectedError( "Alignment not in expected format" );
 666                 }
 667                 if ( $not_ensure != 1 && $seq_array[ $x ][ $y ] !~ /[A-Z]|-/ ) {
 668                     &dieWithUnexpectedError( "Alignment not in expected format (match states must only contain 'A'-'Z' or '-')" );
 669                 }
 670                 print OUTPP "$seq_array[ $x ][ $y ]";
 671             }
 672         }
 673         print OUTPP "\n";
 674     }
 675     close( OUTPP );
 676
 677     return $number_colum, $x_offset, $rf_y;
 678
 679 } ## pfam2phylipMatchOnly
 680
 681
 682
 683 # Returns whether the argument (a String)
 684 # starts with a SWISS-PROT name (SEQN_SPECI).
 685 # Last modified: 06/21/01
 686 sub startsWithSWISS_PROTname {
 687     return ( $_[ 0 ] =~ /^[A-Z0-9]{1,4}_[A-Z0-9]{1,5}/ );
 688 }
 689
 690
 691
 692 # Returns whether the argument starts with XXX.. XXXXX.. and the first
 693 # character is not a "#".
 694 # Last modified: 06/21/01
 695 sub isPfamSequenceLine {
 696     return( !&isPfamCommentLine( $_[ 0 ] )
 697     && &containsPfamNamedSequence( $_[ 0 ] ) );
 698 }
 699
 700
 701
 702 # Returns whether the argument does start with a "#".
 703 # Last modified: 06/21/01
 704 sub isPfamCommentLine {
 705     return ( $_[ 0 ] =~ /^#/ );
 706 }
 707
 708
 709
 710 # Returns whether the argument starts with XXX XXXXX.
 711 # Last modified: 06/21/01
 712 sub containsPfamNamedSequence {
 713     return ( $_[ 0 ] =~ /^\S+\s+\S+/ );
 714 }
 715
 716
 717 # Returns whether the argument starts with XXX XXXXX.
 718 # Last modified: 06/21/01
 719 sub isRFline {
 720     return ( $_[ 0 ] =~ /^#.*RF/ );
 721 }
 722
 723
 724
 725 # Three arguments:
 726 # 1. pairwise distance file
 727 # 2. number of bootstraps
 728 # 3. initial tree: BME, GME or NJ
 729 # Last modified: 2008/12/31
 730 sub executeFastme {
 731     my $inpwd    = $_[ 0 ];
 732     my $bs       = $_[ 1 ];
 733     my $init_opt = $_[ 2 ];
 734
 735     &testForTextFilePresence( $inpwd );
 736     my $command = "";
 737     if ( $bs > 1 ) {
 738         $command = "$FASTME -b $init_opt -i $inpwd -n $bs -s b";
 739     }
 740     else {
 741         $command = "$FASTME -b $init_opt -i $inpwd -s b";
 742     }
 743     print $command;
 744
 745     system( $command );
 746
 747
 748 } ## executeFastme
 749
 750
 751 # Four arguments:
 752 # 1. pairwise distance file
 753 # 2. number of bootstraps
 754 # 3. seed for random number generator
 755 # 4. lower-triangular data matrix? 1: yes; no, otherwise
 756 sub executeNeighbor {
 757     my $inpwd  = $_[ 0 ];
 758     my $bs     = $_[ 1 ];
 759     my $s      = $_[ 2 ];
 760     my $l      = $_[ 3 ];
 761     my $multi  = "";
 762     my $lower  = "";
 763
 764     &testForTextFilePresence( $inpwd );
 765
 766     if (  $bs >= 2 ) {
 767         $multi = "
 768 M
 769 $bs
 770 $s";
 771     }
 772     if ( $l == 1 ) {
 773         $lower = "
 774 L";
 775     }
 776
 777     system( "$NEIGHBOR >/dev/null 2>&1 << !
 778 $inpwd$multi$lower
 779 2
 780 3
 781 Y
 782 !" )
 783     && &dieWithUnexpectedError( "Could not execute \"$NEIGHBOR $inpwd$multi$lower\"" );
 784
 785 } ## executeNeighbor
 786
 787
 788 # Seven arguments:
 789 # 1. pairwise distance file
 790 # 2. number of bootstraps
 791 # 3. seed for random number generator
 792 # 4. number of jumbles for input order
 793 # 5. lower-triangular data matrix? 1: yes; no, otherwise
 794 # 6. FM for Fitch-Margoliash, ME for ME# 6.
 795 # 7. 1 to use globale rearragements
 796 sub executeFitch {
 797     my $inpwd            = $_[ 0 ];
 798     my $bs               = $_[ 1 ];
 799     my $s                = $_[ 2 ];
 800     my $j                = $_[ 3 ];
 801     my $l                = $_[ 4 ];
 802     my $m                = $_[ 5 ];
 803     my $use_global_rearr = $_[ 6 ];
 804     my $jumble = "";
 805     my $multi  = "";
 806     my $lower  = "";
 807     my $method = "";
 808
 809     my $global = "";
 810     if ( $use_global_rearr == 1 ) {
 811         $global = "
 812 G";
 813     }
 814
 815     &testForTextFilePresence( $inpwd );
 816
 817     if ( $m eq "FM" ) {
 818         $method = "";
 819     }
 820     elsif ( $m eq "ME" ) {
 821         $method = "
 822 D";
 823     }
 824     else {
 825         &dieWithUnexpectedError( "method for FITCH must be either FM or ME" );
 826     }
 827
 828     if ( $j >= 1 ) {
 829         $jumble = "
 830 J
 831 $s
 832 $j";
 833     }
 834
 835     if (  $bs >= 2 ) {
 836         $multi = "
 837 M
 838 $bs
 839 $s";
 840     }
 841     if ( $l == 1 ) {
 842         $lower = "
 843 L";
 844     }
 845
 846     # jumble must be set BEFORE multi!
 847     system( "$FITCH 2>&1 << !
 848 $inpwd$method$global$jumble$multi$lower
 849 3
 850 Y
 851 !" )
 852     && &dieWithUnexpectedError( "Could not execute \"$FITCH $inpwd$method$global$jumble$multi$lower\"" );
 853     # 3: Do NOT print out tree
 854
 855 } ## executeFitch
 856
 857
 858
 859 # Two arguments:
 860 # 1. pairwise distance file
 861 # 2. outfile
 862 sub executeBionj {
 863     my $inpwd    = $_[ 0 ];
 864     my $out      = $_[ 1 ];
 865
 866     &testForTextFilePresence( $inpwd );
 867     my $command = "$BIONJ $inpwd $out";
 868
 869     system( $command )
 870     && &dieWithUnexpectedError( $command );
 871
 872 }
 873
 874 # Four arguments:
 875 # 1. (effective) sequence length
 876 # 2. (effective) number of bases
 877 # 3. pairwise distance file
 878 # 4. outfile
 879 sub executeWeighbor {
 880     my $L = $_[ 0 ];
 881     my $b = $_[ 1 ];
 882     my $i = $_[ 2 ];
 883     my $o = $_[ 3 ];
 884
 885     &testForTextFilePresence( $i );
 886     my $command = "$WEIGHBOR -L $L -b $b -i $i -o $o";
 887
 888     system( $command )
 889     && &dieWithUnexpectedError( $command );
 890
 891 }
 892
 893 # Six arguments:
 894 # 1. DNA or Amino-Acids sequence filename (PHYLIP format)
 895 # 2. number of data sets to analyse (ex:3)
 896 # 3. Model: JTT | MtREV | Dayhoff | WAG | VT | DCMut | Blosum62 (Amino-Acids)
 897 # 4. number of relative substitution rate categories (ex:4), positive integer
 898 # 5. starting tree filename (Newick format), your tree filename | BIONJ for a distance-based tree
 899 # 6. 1 to estimate proportion of invariable sites, otherwise, fixed proportion "0.0" is used
 900 # PHYML produces several results files :
 901 # <sequence file name>_phyml_lk.txt : likelihood value(s)
 902 # <sequence file name>_phyml_tree.txt : inferred tree(s)
 903 # <sequence file name>_phyml_stat.txt : detailed execution stats
 904 sub executePhyml {
 905     my $sequences = $_[ 0 ];
 906     my $data_sets = $_[ 1 ];
 907     my $model     = $_[ 2 ];
 908     my $nb_categ  = $_[ 3 ];
 909     my $tree      = $_[ 4 ];
 910     my $estimate_invar_sites = $_[ 5 ];
 911
 912     if ( $data_sets < 1 ) {
 913         $data_sets = 1
 914     }
 915
 916     my $invar          = "0.0";   # proportion of invariable sites,
 917                                   # a fixed value (ex:0.0) | e to get the maximum likelihood estimate
 918     if ( $estimate_invar_sites == 1 ) {
 919         $invar = "e";
 920     }
 921
 922     my $data_type      = "1";     # 0 = DNA | 1 = Amino-Acids
 923     my $format         = "i";     # i = interleaved sequence format | s = sequential
 924     my $bootstrap_sets = "0";     # number of bootstrap data sets to generate (ex:2)
 925                                   # only works with one data set to analyse
 926
 927     my $alpha          = "e";     # gamma distribution parameter,
 928                                   # a fixed value (ex:1.0) | e to get the maximum likelihood estimate
 929
 930     my $opt_topology   = "y";     # optimise tree topology ? y | n
 931     my $opt_lengths    = "y";     # optimise branch lengths and rate parameters ? y | n
 932
 933     if ( $data_sets > 1 ) {
 934         # No need to calc branch lengths for bootstrapped analysis
 935         $opt_lengths = "n";
 936     }
 937
 938     &testForTextFilePresence( $sequences );
 939     my $command = "$PHYML $sequences $data_type $format $data_sets $bootstrap_sets $model $invar $nb_categ $alpha $tree $opt_topology $opt_lengths";
 940
 941     print( "\n$command\n");
 942
 943     system( $command )
 944     && &dieWithUnexpectedError( $command );
 945
 946 }
 947
 948
 949
 950
 951 # Four arguments:
 952 # 1. name of alignment file (in correct format!)
 953 # 2. number of bootstraps
 954 # 3. jumbles: 0: do not jumble; >=1 number of jumbles
 955 # 4. seed for random number generator
 956 sub executeProtpars {
 957     my $align  = $_[ 0 ];
 958     my $bs     = $_[ 1 ];
 959     my $rand   = $_[ 2 ];
 960     my $s      = $_[ 3 ];
 961     my $jumble = "";
 962     my $multi  = "";
 963
 964
 965     &testForTextFilePresence( $align );
 966
 967     if ( $bs > 1 && $rand < 1 ) {
 968         $rand = 1;
 969     }
 970
 971     if ( $rand >= 1 ) {
 972         $jumble = "
 973 J
 974 $s
 975 $rand";
 976     }
 977
 978     if (  $bs > 1 ) {
 979         $multi = "
 980 M
 981 D
 982 $bs";
 983     }
 984
 985     system( "$PROTPARS  2>&1 << !
 986 $align$jumble$multi
 987 Y
 988 !" )
 989     && &dieWithUnexpectedError( "Could not execute \"$PROTPARS $align$jumble$multi\"" );
 990     # 3: Do NOT print out tree
 991
 992
 993     return;
 994
 995 } ## executeProtpars
 996
 997
 998 # "Model of substitution" order for TREE-PUZZLE 5.2:
 999 # For amino acids:
1000 # Auto
1001 # m -> Dayhoff (Dayhoff et al. 1978)
1002 # m -> JTT (Jones et al. 1992)
1003 # m -> mtREV24 (Adachi-Hasegawa 1996)
1004 # m -> BLOSUM62 (Henikoff-Henikoff 92)
1005 # m -> VT (Mueller-Vingron 2000)
1006 # m -> WAG (Whelan-Goldman 2000)
1007 # m -> Auto
1008 #
1009 # For nucleotides:
1010 # HKY (Hasegawa et al. 1985)
1011 # m -> TN (Tamura-Nei 1993)
1012 # m -> GTR (e.g. Lanave et al. 1980)
1013 # m -> SH (Schoeniger-von Haeseler 1994)
1014 # m -> HKY (Hasegawa et al. 1985)
1015 #
1016 # One argument:
1017 # matrix option:
1018 # 0 = JTT
1019 # 2 = BLOSUM 62
1020 # 3 = mtREV24
1021 # 5 = VT
1022 # 6 = WAG
1023 # 7 = auto
1024 # 9 = HKY [na]
1025 # 10 = TN [na]
1026 # 11 = GTR [na]
1027 # 12 = SH [na]
1028 #     PAM otherwise
1029 # Last modified: 17/04/26
1030 sub setModelForPuzzle {
1031     my $matrix_option = $_[ 0 ];
1032     my $matr          = "";
1033
1034     if ( $matrix_option == 0 || $matrix_option == 11 ) { # JTT or GTR
1035         $matr = "
1036 m
1037 m";
1038     }
1039     elsif ( $matrix_option == 2 ) { # BLOSUM 62
1040         $matr = "
1041 m
1042 m
1043 m
1044 m";
1045     }
1046     elsif ( $matrix_option == 3 || $matrix_option == 12) { # mtREV24 or SH
1047         $matr = "
1048 m
1049 m
1050 m";
1051     }
1052     elsif ( $matrix_option == 5 ) { # VT
1053         $matr = "
1054 m
1055 m
1056 m
1057 m
1058 m";
1059     }
1060     elsif ( $matrix_option == 6 ) { # WAG
1061         $matr = "
1062 m
1063 m
1064 m
1065 m
1066 m
1067 m";
1068     }
1069     elsif ( $matrix_option == 7 || $matrix_option == 9 ) { # auto or HKY
1070         $matr = "";
1071     }
1072     elsif ( $matrix_option == 10 ) { # TN
1073         $matr = "
1074 m"
1075     }
1076     else { # PAM
1077         $matr = "
1078 m"
1079     }
1080
1081     return $matr;
1082
1083 } ## setModelForPuzzle
1084
1085 # One argument:
1086 # Model of rate heterogeneity:
1087 #    1 for "8 Gamma distributed rates"
1088 #    2 for "Two rates (1 invariable + 1 variable)"
1089 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1090 #    otherwise: Uniform rate
1091 # Last modified: 09/08/03
1092 sub setRateHeterogeneityOptionForPuzzle {
1093     my $rate_heterogeneity_option  = $_[ 0 ];
1094     my $opt                        = "";
1095
1096     if ( $rate_heterogeneity_option == 1 ) {
1097         $opt = "
1098 w";
1099     }
1100     elsif ( $rate_heterogeneity_option == 2 ) {
1101         $opt = "
1102 w
1103 w";
1104     }
1105     elsif ( $rate_heterogeneity_option == 3 ) {
1106         $opt = "
1107 w
1108 w
1109 w";
1110     }
1111     else {
1112         $opt = "";
1113     }
1114
1115     return $opt;
1116 } ## setRateHeterogeneityOptionForPuzzle
1117
1118
1119 # One argument:
1120 # Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1121 # Last modified: 09/08/03
1122 sub setParameterEstimatesOptionForPuzzle {
1123     my $parameter_estimates_option  = $_[ 0 ];
1124     my $opt                         = "";
1125
1126     if ( $parameter_estimates_option == 1 ) {
1127         $opt = "
1128 e";
1129     }
1130     else {
1131         $opt = "";
1132     }
1133
1134     return $opt;
1135 } ## setParameterEstimatesOptionForPuzzle
1136
1137
1138
1139 # three/four/five arguments:
1140 # 1. Name of inputfile
1141 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1142 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1143 # 3. Number of sequences in alignment
1144 # 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1145 # 5. Model of rate heterogeneity:
1146 #    1 for "8 Gamma distributed rates"
1147 #    2 for "Two rates (1 invariable + 1 variable)"
1148 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1149 #    otherwise: Uniform rate
1150 sub executePuzzleBootstrapped {
1151     my $in                         = $_[ 0 ];
1152     my $matrix_option              = $_[ 1 ];
1153     my $number_of_seqs             = $_[ 2 ];
1154     my $parameter_estimates_option = $_[ 3 ];
1155     my $rate_heterogeneity_option  = $_[ 4 ];
1156
1157     my $l             = 0;
1158     my $slen          = 0;
1159     my $counter       = 0;
1160     my $mat           = "";
1161     my $est           = "";
1162     my $rate          = "";
1163     my $a             = "";
1164     my @a             = ();
1165
1166     &testForTextFilePresence( $in );
1167
1168     open( GRP, "<$in" ) || die "\n\n$0: Unexpected error: Cannot open file <<$in>>: $!";
1169     while( <GRP> ) {
1170         if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
1171             $counter++;
1172         }
1173     }
1174     close( GRP );
1175
1176     $l   = `cat $in | wc -l`;
1177     $slen   = $l / $counter;
1178
1179     system( "split --suffix-length=4 -$slen $in $in.splt." )
1180     && die "\n\n$0: executePuzzleDQObootstrapped: Could not execute \"split --suffix-length=4 -$slen $in $in.splt.\": $!";
1181
1182     @a = <$in.splt.*>;
1183
1184     $mat = setModelForPuzzle( $matrix_option );
1185      if ( $parameter_estimates_option ) {
1186         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
1187     }
1188     if ( $rate_heterogeneity_option ) {
1189         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
1190     }
1191
1192     my $k="";
1193     if (  $number_of_seqs <= 257 ) {
1194         $k = "k";
1195     }
1196
1197     foreach $a ( @a ) {
1198         print "-".$a."\n";
1199         system( "$PUZZLE $a << !
1200 $k
1201 k
1202 k$mat$est$rate
1203 y
1204 !" )
1205         && die "$0: Could not execute \"$PUZZLE $a\"";
1206
1207         system( "cat $a.dist >> $in.dist" )
1208         && die "$0: Could not execute \"cat outdist >> $in.dist\"";
1209
1210         unlink( $a, $a.".dist", $a.".puzzle" );
1211     }
1212
1213     return;
1214
1215 } ## executePuzzleBootstrapped
1216
1217
1218
1219
1220
1221 # three/four/five arguments:
1222 # 1. Name of inputfile
1223 # 2. Matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1224 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1225 # 3. Number of sequences in alignment
1226 # 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1227 # 5. Model of rate heterogeneity:
1228 #    1 for "8 Gamma distributed rates"
1229 #    2 for "Two rates (1 invariable + 1 variable)"
1230 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1231 #    otherwise: Uniform rate
1232 sub executePuzzle {
1233     my $in                         = $_[ 0 ];
1234     my $matrix_option              = $_[ 1 ];
1235     my $number_of_seqs             = $_[ 2 ];
1236     my $parameter_estimates_option = $_[ 3 ];
1237     my $rate_heterogeneity_option  = $_[ 4 ];
1238     my $mat                        = "";
1239     my $est                        = "";
1240     my $rate                       = "";
1241
1242     &testForTextFilePresence( $in );
1243
1244     $mat = &setModelForPuzzle( $matrix_option );
1245     if ( $parameter_estimates_option ) {
1246         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
1247     }
1248     if ( $rate_heterogeneity_option ) {
1249         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
1250     }
1251
1252     my $k="";
1253     if (  $number_of_seqs <= 257 ) {
1254         $k = "k";
1255     }
1256
1257
1258     system( "$PUZZLE $in << !
1259 $k
1260 k
1261 k$mat$est$rate
1262 y
1263 !" )
1264     && die "$0: Could not execute \"$PUZZLE\"";
1265
1266     return;
1267
1268 } ## executePuzzle
1269
1270
1271
1272
1273 # Preparation of the pwd file
1274 sub addDistsToQueryToPWDfile {
1275     my $pwd_file          = $_[ 0 ];
1276     my $disttoquery_file  = $_[ 1 ];
1277     my $outfile           = $_[ 2 ];
1278     my $name_of_query     = $_[ 3 ];
1279     my $name_of_query_    = "";
1280     my $return_line_pwd   = "";
1281     my $return_line_dq    = "";
1282     my $num_of_sqs        = 0;
1283     my $block             = 0;
1284     my $name_from_pwd     = "X";
1285     my $name_from_dq      = "Y";
1286     my @dists_to_query    = ();
1287     my $i                 = 0;
1288
1289     &testForTextFilePresence( $pwd_file );
1290     &testForTextFilePresence( $disttoquery_file );
1291
1292     $name_of_query_ = $name_of_query;
1293     for ( my $j = 0; $j <= ( $LENGTH_OF_NAME - length( $name_of_query ) - 1 ); ++$j ) {
1294         $name_of_query_ .= " ";
1295     }
1296
1297     open( OUT_AD, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
1298     open( IN_PWD, "$pwd_file" ) || &dieWithUnexpectedError( "Cannot open file \"$pwd_file\"" );
1299     open( IN_DQ, "$disttoquery_file" ) || &dieWithUnexpectedError( "Cannot open file \"$disttoquery_file\"" );
1300
1301     W: while ( $return_line_pwd = <IN_PWD> ) {
1302
1303
1304         if ( $return_line_pwd =~ /^\s*(\d+)\s*$/ ) {
1305             $num_of_sqs = $1;
1306             $num_of_sqs++;
1307             if ( $block > 0 ) {
1308                 print OUT_AD "$name_of_query_  ";
1309                 for ( my $j = 0; $j < $i; ++$j ) {
1310                     print OUT_AD "$dists_to_query[ $j ]  ";
1311                 }
1312                 print OUT_AD "0.0\n";
1313             }
1314             print OUT_AD "  $num_of_sqs\n";
1315             $block++;
1316             @dists_to_query = ();
1317             $i = 0;
1318         }
1319
1320         if ( $block == 1
1321         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
1322             $name_from_pwd = $1;
1323
1324             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1325                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1326             }
1327
1328             if ( $return_line_dq !~ /\S/ ) {
1329                 if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1330                     &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1331                 }
1332             }
1333             $return_line_dq =~ /^\s*(\S+)\s+(\S+)/;
1334             $name_from_dq = $1;
1335             $dists_to_query[ $i++ ] = $2;
1336
1337
1338             if ( $name_from_pwd ne $name_from_dq ) {
1339                 &dieWithUnexpectedError( "Order of sequence names in \"$pwd_file\" and \"$disttoquery_file\" is not the same" );
1340             }
1341             print OUT_AD $return_line_pwd;
1342
1343         }
1344         elsif ( $block > 1
1345         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
1346             $name_from_pwd = $1;
1347             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1348                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1349             }
1350             if ( $return_line_dq !~ /\S/ ) {
1351                 if ( !defined( $return_line_dq = <IN_DQ>) ) {
1352                     &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1353                 }
1354             }
1355             $return_line_dq =~ /^\s*\S+\s+(\S+)/;
1356             $dists_to_query[ $i++ ] = $1;
1357             print OUT_AD $return_line_pwd;
1358         }
1359     }
1360     print OUT_AD "$name_of_query_  ";
1361     for ( my $j = 0; $j < $i; ++$j ) {
1362         print OUT_AD "$dists_to_query[ $j ]  ";
1363     }
1364     print OUT_AD "0.0\n";
1365
1366     close( OUT_AD );
1367     close( IN_PWD );
1368     close( IN_DQ );
1369     return $block;
1370
1371 } ## addDistsToQueryToPWDfile
1372
1373
1374
1375
1376 # Three arguments:
1377 # 1. HMMER model db
1378 # 2. name of HMM
1379 # 3. outputfile name
1380 # Last modified: 02/27/01
1381 sub executeHmmfetch {
1382
1383     my $db      = $_[ 0 ];
1384     my $name    = $_[ 1 ];
1385     my $outfile = $_[ 2 ];
1386
1387     system( "$HMMFETCH $db $name > $outfile" )
1388     && &dieWithUnexpectedError( "Could not execute \"$HMMFETCH $db $name > $outfile\"" );
1389     return;
1390
1391 } ## executeHmmfetch
1392
1393
1394
1395 # Checks wether a file is present, not empty and a plain textfile.
1396 # One argument: name of file.
1397 # Last modified: 07/07/01
1398 sub testForTextFilePresence {
1399     my $file = $_[ 0 ];
1400     unless ( ( -s $file ) && ( -f $file ) && ( -T $file ) ) {
1401         dieWithUnexpectedError( "File \"$file\" does not exist, is empty, or is not a plain textfile" );
1402     }
1403 } ## testForTextFilePresence
1404
1405
1406 # Last modified: 02/21/03
1407 sub addSlashAtEndIfNotPresent {
1408     my $filename = $_[ 0 ];
1409     $filename =~ s/\s+//g;
1410     unless ( $filename =~ /\/$/ ) {
1411        $filename = $filename."/";
1412     }
1413     return $filename;
1414 } ## addSlashAtEndIfNotPresent
1415
1416
1417
1418 # Last modified: 02/15/02
1419 sub exitWithWarning {
1420
1421     my $text = $_[ 0 ];
1422     if ( defined( $_[ 1 ] ) && $_[ 1 ] == 1 ) {
1423         print( "<H4 class=\"error\">user error</H4>\n" );
1424         print( "<P>\n" );
1425         print( "<B>$text</B>\n" );
1426         print( "</P>\n" );
1427         print( "<P> &nbsp</P>\n" );
1428     }
1429     else {
1430         print( "\n\n$text\n\n" );
1431     }
1432
1433     exit( 0 );
1434
1435 } ## exit_with_warning
1436
1437
1438
1439 # Last modified: 02/15/02
1440 sub dieWithUnexpectedError {
1441
1442     my $text = $_[ 0 ];
1443
1444     die( "\n\n$0:\nUnexpected error (should not have happened):\n$text\n$!\n\n" );
1445
1446 } ## dieWithUnexpectedError
1447
1448
1449
1450 1;