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/czmasek/SOFTWARE/";
 170
 171
 172 # Java virtual machine:
 173 # ---------------------
 174 our $JAVA                      = $SOFTWARE_DIR."JAVA/jdk1.6.0_03/bin/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."PHYLIP/phylip-3.68/src/seqboot";
 185 our $NEIGHBOR                  = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/neighbor";
 186 our $PROTPARS                  = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/protpars";
 187 our $PROML                     = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/proml";
 188 our $FITCH                     = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/fitch";
 189 our $CONSENSE                  = $SOFTWARE_DIR."PHYLIP/phylip-3.68/src/consense";
 190 our $PHYLIP_VERSION            = "3.68";
 191
 192 # TREE-PUZZLE:
 193 # ------------
 194 our $PUZZLE                    = $SOFTWARE_DIR."TREE_PUZZLE/tree-puzzle-5.2/src/puzzle";
 195 our $PUZZLE_VERSION            = "5.2";
 196
 197 # FASTME:
 198 # -----------------------------------------------------
 199 our $FASTME                    = $SOFTWARE_DIR."FASTME/fastme2.0/fastme";
 200 our $FASTME_VERSION            = "2.0";
 201
 202 # BIONJ:
 203 # -----------------------------------------------------
 204 our $BIONJ                    = $SOFTWARE_DIR."BIONJ/bionj";
 205 our $BIONJ_VERSION            = "[1997]";
 206
 207 # WEIGHBOR:
 208 # -----------------------------------------------------
 209 our $WEIGHBOR                 = $SOFTWARE_DIR."WEIGHBOR/Weighbor/weighbor";
 210 our $WEIGHBOR_VERSION         = "1.2.1";
 211
 212 # PHYML:
 213 # -----------------------------------------------------
 214 our $PHYML                    = $SOFTWARE_DIR."PHYML/phyml_v2.4.4/exe/phyml_linux";
 215 our $PHYML_VERSION            = "2.4.4";
 216
 217 # RAXML:
 218 # -----------------------------------------------------
 219 our $RAXML                    = $SOFTWARE_DIR."RAXML/RAxML-7.0.4/raxmlHPC";
 220 our $RAXML_VERSION            = "7.0.4";
 221
 222
 223 # forester.jar. This jar file is currently available at: http://www.phylosoft.org
 224 # -------------------------------------------------------------------------------
 225
 226 our $FORESTER_JAR             = $SOFTWARE_DIR."FORESTER/DEV/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 > 0 ) {
 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
 999 # "Model of substitution" order for DQO TREE-PUZZLE 5.0:
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 # One argument:
1009 # matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1010 # 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1011 # Last modified: 07/07/01
1012 sub setModelForPuzzle {
1013     my $matrix_option = $_[ 0 ];
1014     my $matr          = "";
1015
1016     if ( $matrix_option == 0 ) { # JTT
1017         $matr = "
1018 m
1019 m";
1020     }
1021     elsif ( $matrix_option == 2 ) { # BLOSUM 62
1022         $matr = "
1023 m
1024 m
1025 m
1026 m";
1027     }
1028     elsif ( $matrix_option == 3 ) { # mtREV24
1029         $matr = "
1030 m
1031 m
1032 m";
1033     }
1034     elsif ( $matrix_option == 5 ) { # VT
1035         $matr = "
1036 m
1037 m
1038 m
1039 m
1040 m";
1041     }
1042     elsif ( $matrix_option == 6 ) { # WAG
1043         $matr = "
1044 m
1045 m
1046 m
1047 m
1048 m
1049 m";
1050     }
1051     elsif ( $matrix_option == 7 ) { # auto
1052         $matr = "";
1053     }
1054     else { # PAM
1055         $matr = "
1056 m"
1057     }
1058
1059     return $matr;
1060
1061 } ## setModelForPuzzle
1062
1063 # One argument:
1064 # Model of rate heterogeneity:
1065 #    1 for "8 Gamma distributed rates"
1066 #    2 for "Two rates (1 invariable + 1 variable)"
1067 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1068 #    otherwise: Uniform rate
1069 # Last modified: 09/08/03
1070 sub setRateHeterogeneityOptionForPuzzle {
1071     my $rate_heterogeneity_option  = $_[ 0 ];
1072     my $opt                        = "";
1073
1074     if ( $rate_heterogeneity_option == 1 ) {
1075         $opt = "
1076 w";
1077     }
1078     elsif ( $rate_heterogeneity_option == 2 ) {
1079         $opt = "
1080 w
1081 w";
1082     }
1083     elsif ( $rate_heterogeneity_option == 3 ) {
1084         $opt = "
1085 w
1086 w
1087 w";
1088     }
1089     else {
1090         $opt = "";
1091     }
1092
1093     return $opt;
1094 } ## setRateHeterogeneityOptionForPuzzle
1095
1096
1097 # One argument:
1098 # Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1099 # Last modified: 09/08/03
1100 sub setParameterEstimatesOptionForPuzzle {
1101     my $parameter_estimates_option  = $_[ 0 ];
1102     my $opt                         = "";
1103
1104     if ( $parameter_estimates_option == 1 ) {
1105         $opt = "
1106 e";
1107     }
1108     else {
1109         $opt = "";
1110     }
1111
1112     return $opt;
1113 } ## setParameterEstimatesOptionForPuzzle
1114
1115
1116
1117 # three/four/five arguments:
1118 # 1. Name of inputfile
1119 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1120 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1121 # 3. Number of sequences in alignment
1122 # 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1123 # 5. Model of rate heterogeneity:
1124 #    1 for "8 Gamma distributed rates"
1125 #    2 for "Two rates (1 invariable + 1 variable)"
1126 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1127 #    otherwise: Uniform rate
1128 sub executePuzzleBootstrapped {
1129     my $in                         = $_[ 0 ];
1130     my $matrix_option              = $_[ 1 ];
1131     my $number_of_seqs             = $_[ 2 ];
1132     my $parameter_estimates_option = $_[ 3 ];
1133     my $rate_heterogeneity_option  = $_[ 4 ];
1134
1135     my $l             = 0;
1136     my $slen          = 0;
1137     my $counter       = 0;
1138     my $mat           = "";
1139     my $est           = "";
1140     my $rate          = "";
1141     my $a             = "";
1142     my @a             = ();
1143
1144     &testForTextFilePresence( $in );
1145
1146     open( GRP, "<$in" ) || die "\n\n$0: Unexpected error: Cannot open file <<$in>>: $!";
1147     while( <GRP> ) {
1148         if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
1149             $counter++;
1150         }
1151     }
1152     close( GRP );
1153
1154     $l   = `cat $in | wc -l`;
1155     $slen   = $l / $counter;
1156
1157     system( "split --suffix-length=4 -$slen $in $in.splt." )
1158     && die "\n\n$0: executePuzzleDQObootstrapped: Could not execute \"split --suffix-length=4 -$slen $in $in.splt.\": $!";
1159
1160     @a = <$in.splt.*>;
1161
1162     $mat = setModelForPuzzle( $matrix_option );
1163      if ( $parameter_estimates_option ) {
1164         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
1165     }
1166     if ( $rate_heterogeneity_option ) {
1167         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
1168     }
1169
1170     my $k="";
1171     if (  $number_of_seqs <= 257 ) {
1172         $k = "k";
1173     }
1174
1175     foreach $a ( @a ) {
1176         print "-".$a."\n";
1177         system( "$PUZZLE $a << !
1178 $k
1179 k
1180 k$mat$est$rate
1181 y
1182 !" )
1183         && die "$0: Could not execute \"$PUZZLE $a\"";
1184
1185         system( "cat $a.dist >> $in.dist" )
1186         && die "$0: Could not execute \"cat outdist >> $in.dist\"";
1187
1188         unlink( $a, $a.".dist", $a.".puzzle" );
1189     }
1190
1191     return;
1192
1193 } ## executePuzzleBootstrapped
1194
1195
1196
1197
1198
1199 # three/four/five arguments:
1200 # 1. Name of inputfile
1201 # 2. Matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
1202 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
1203 # 3. Number of sequences in alignment
1204 # 4. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
1205 # 5. Model of rate heterogeneity:
1206 #    1 for "8 Gamma distributed rates"
1207 #    2 for "Two rates (1 invariable + 1 variable)"
1208 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
1209 #    otherwise: Uniform rate
1210 sub executePuzzle {
1211     my $in                         = $_[ 0 ];
1212     my $matrix_option              = $_[ 1 ];
1213     my $number_of_seqs             = $_[ 2 ];
1214     my $parameter_estimates_option = $_[ 3 ];
1215     my $rate_heterogeneity_option  = $_[ 4 ];
1216     my $mat                        = "";
1217     my $est                        = "";
1218     my $rate                       = "";
1219
1220     &testForTextFilePresence( $in );
1221
1222     $mat = &setModelForPuzzle( $matrix_option );
1223     if ( $parameter_estimates_option ) {
1224         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
1225     }
1226     if ( $rate_heterogeneity_option ) {
1227         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
1228     }
1229
1230     my $k="";
1231     if (  $number_of_seqs <= 257 ) {
1232         $k = "k";
1233     }
1234
1235
1236     system( "$PUZZLE $in << !
1237 $k
1238 k
1239 k$mat$est$rate
1240 y
1241 !" )
1242     && die "$0: Could not execute \"$PUZZLE\"";
1243
1244     return;
1245
1246 } ## executePuzzle
1247
1248
1249
1250
1251 # Preparation of the pwd file
1252 sub addDistsToQueryToPWDfile {
1253     my $pwd_file          = $_[ 0 ];
1254     my $disttoquery_file  = $_[ 1 ];
1255     my $outfile           = $_[ 2 ];
1256     my $name_of_query     = $_[ 3 ];
1257     my $name_of_query_    = "";
1258     my $return_line_pwd   = "";
1259     my $return_line_dq    = "";
1260     my $num_of_sqs        = 0;
1261     my $block             = 0;
1262     my $name_from_pwd     = "X";
1263     my $name_from_dq      = "Y";
1264     my @dists_to_query    = ();
1265     my $i                 = 0;
1266
1267     &testForTextFilePresence( $pwd_file );
1268     &testForTextFilePresence( $disttoquery_file );
1269
1270     $name_of_query_ = $name_of_query;
1271     for ( my $j = 0; $j <= ( $LENGTH_OF_NAME - length( $name_of_query ) - 1 ); ++$j ) {
1272         $name_of_query_ .= " ";
1273     }
1274
1275     open( OUT_AD, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
1276     open( IN_PWD, "$pwd_file" ) || &dieWithUnexpectedError( "Cannot open file \"$pwd_file\"" );
1277     open( IN_DQ, "$disttoquery_file" ) || &dieWithUnexpectedError( "Cannot open file \"$disttoquery_file\"" );
1278
1279     W: while ( $return_line_pwd = <IN_PWD> ) {
1280
1281
1282         if ( $return_line_pwd =~ /^\s*(\d+)\s*$/ ) {
1283             $num_of_sqs = $1;
1284             $num_of_sqs++;
1285             if ( $block > 0 ) {
1286                 print OUT_AD "$name_of_query_  ";
1287                 for ( my $j = 0; $j < $i; ++$j ) {
1288                     print OUT_AD "$dists_to_query[ $j ]  ";
1289                 }
1290                 print OUT_AD "0.0\n";
1291             }
1292             print OUT_AD "  $num_of_sqs\n";
1293             $block++;
1294             @dists_to_query = ();
1295             $i = 0;
1296         }
1297
1298         if ( $block == 1
1299         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
1300             $name_from_pwd = $1;
1301
1302             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1303                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1304             }
1305
1306             if ( $return_line_dq !~ /\S/ ) {
1307                 if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1308                     &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1309                 }
1310             }
1311             $return_line_dq =~ /^\s*(\S+)\s+(\S+)/;
1312             $name_from_dq = $1;
1313             $dists_to_query[ $i++ ] = $2;
1314
1315
1316             if ( $name_from_pwd ne $name_from_dq ) {
1317                 &dieWithUnexpectedError( "Order of sequence names in \"$pwd_file\" and \"$disttoquery_file\" is not the same" );
1318             }
1319             print OUT_AD $return_line_pwd;
1320
1321         }
1322         elsif ( $block > 1
1323         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
1324             $name_from_pwd = $1;
1325             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1326                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
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             $dists_to_query[ $i++ ] = $1;
1335             print OUT_AD $return_line_pwd;
1336         }
1337     }
1338     print OUT_AD "$name_of_query_  ";
1339     for ( my $j = 0; $j < $i; ++$j ) {
1340         print OUT_AD "$dists_to_query[ $j ]  ";
1341     }
1342     print OUT_AD "0.0\n";
1343
1344     close( OUT_AD );
1345     close( IN_PWD );
1346     close( IN_DQ );
1347     return $block;
1348
1349 } ## addDistsToQueryToPWDfile
1350
1351
1352
1353
1354 # Three arguments:
1355 # 1. HMMER model db
1356 # 2. name of HMM
1357 # 3. outputfile name
1358 # Last modified: 02/27/01
1359 sub executeHmmfetch {
1360
1361     my $db      = $_[ 0 ];
1362     my $name    = $_[ 1 ];
1363     my $outfile = $_[ 2 ];
1364
1365     system( "$HMMFETCH $db $name > $outfile" )
1366     && &dieWithUnexpectedError( "Could not execute \"$HMMFETCH $db $name > $outfile\"" );
1367     return;
1368
1369 } ## executeHmmfetch
1370
1371
1372
1373 # Checks wether a file is present, not empty and a plain textfile.
1374 # One argument: name of file.
1375 # Last modified: 07/07/01
1376 sub testForTextFilePresence {
1377     my $file = $_[ 0 ];
1378     unless ( ( -s $file ) && ( -f $file ) && ( -T $file ) ) {
1379         dieWithUnexpectedError( "File \"$file\" does not exist, is empty, or is not a plain textfile" );
1380     }
1381 } ## testForTextFilePresence
1382
1383
1384 # Last modified: 02/21/03
1385 sub addSlashAtEndIfNotPresent {
1386     my $filename = $_[ 0 ];
1387     $filename =~ s/\s+//g;
1388     unless ( $filename =~ /\/$/ ) {
1389        $filename = $filename."/";
1390     }
1391     return $filename;
1392 } ## addSlashAtEndIfNotPresent
1393
1394
1395
1396 # Last modified: 02/15/02
1397 sub exitWithWarning {
1398
1399     my $text = $_[ 0 ];
1400     if ( defined( $_[ 1 ] ) && $_[ 1 ] == 1 ) {
1401         print( "<H4 class=\"error\">user error</H4>\n" );
1402         print( "<P>\n" );
1403         print( "<B>$text</B>\n" );
1404         print( "</P>\n" );
1405         print( "<P> &nbsp</P>\n" );
1406     }
1407     else {
1408         print( "\n\n$text\n\n" );
1409     }
1410
1411     exit( 0 );
1412
1413 } ## exit_with_warning
1414
1415
1416
1417 # Last modified: 02/15/02
1418 sub dieWithUnexpectedError {
1419
1420     my $text = $_[ 0 ];
1421
1422     die( "\n\n$0:\nUnexpected error (should not have happened):\n$text\n$!\n\n" );
1423
1424 } ## dieWithUnexpectedError
1425
1426
1427
1428 1;