forester/archive/perl/rio_module.pm

   1 # Copyright (C) 2002-2003 Washington University School of Medicine
   2 # and Howard Hughes Medical Institute
   3 # All rights reserved
   4 #
   5 # Author: Christian M. Zmasek
   6 # zmasek@genetics.wustl.edu
   7 # http://www.genetics.wustl.edu/eddy/people/zmasek/
   8 #
   9 # Last modified 03/13/03
  10
  11
  12 package rio_module2;
  13 use strict;
  14 require Exporter;
  15
  16 our $VERSION = 3.20;
  17
  18 our @ISA    = qw( Exporter );
  19
  20 our @EXPORT = qw( executeConsense
  21                   executeMakeTree
  22                   executePuzzleDQO
  23                   executePuzzleDQObootstrapped
  24                   pfam2phylipMatchOnly
  25                   startsWithSWISS_PROTname
  26                   isPfamSequenceLine
  27                   isPfamCommentLine
  28                   containsPfamNamedSequence
  29                   isRFline
  30                   executeNeighbor
  31                   executeProtpars
  32                   setModelForPuzzle
  33                   setRateHeterogeneityOptionForPuzzle
  34                   setParameterEstimatesOptionForPuzzle
  35                   executePuzzleBootstrapped
  36                   executePuzzle
  37                   executeHmmfetch
  38                   addDistsToQueryToPWDfile
  39                   testForTextFilePresence
  40                   exitWithWarning
  41                   dieWithUnexpectedError
  42                   addSlashAtEndIfNotPresent
  43                   $LENGTH_OF_NAME
  44                   $MIN_NUMBER_OF_AA
  45                   $TREMBL_ACDEOS_FILE
  46                   $SWISSPROT_ACDEOS_FILE
  47                   $SPECIES_NAMES_FILE
  48                   $SPECIES_TREE_FILE_DEFAULT
  49                   $MULTIPLE_TREES_FILE_SUFFIX
  50                   $LOG_FILE_SUFFIX
  51                   $ALIGN_FILE_SUFFIX
  52                   $TREE_FILE_SUFFIX
  53                   $ADDITION_FOR_RIO_ANNOT_TREE
  54                   $SUFFIX_PWD
  55                   $SUFFIX_BOOT_STRP_POS
  56                   $SUFFIX_PWD_NOT_BOOTS
  57                   $SUFFIX_HMM
  58                   $MATRIX_FOR_PWD
  59                   $RIO_PWD_DIRECTORY
  60                   $RIO_BSP_DIRECTORY
  61                   $RIO_NBD_DIRECTORY
  62                   $RIO_ALN_DIRECTORY
  63                   $RIO_HMM_DIRECTORY
  64                   $PFAM_FULL_DIRECTORY
  65                   $PFAM_SEED_DIRECTORY
  66                   $PRIOR_FILE_DIR
  67                   $PFAM_HMM_DB
  68                   $SEQBOOT
  69                   $NEIGHBOR
  70                   $PROTPARS
  71                   $CONSENSE
  72                   $PUZZLE
  73                   $HMMALIGN
  74                   $HMMSEARCH
  75                   $HMMBUILD
  76                   $HMMFETCH
  77                   $SFE
  78                   $HMMCALIBRATE
  79                   $P7EXTRACT
  80                   $MULTIFETCH
  81                   $BOOTSTRAP_CZ
  82                   $BOOTSTRAP_CZ_PL
  83                   $TRANSFERSBRANCHLENGHTS
  84                   $MAKETREE
  85                   $RIO_PL
  86                   $DORIO
  87                   $PUZZLE_DQO
  88                   $BOOTSTRAPS
  89                   $PATH_TO_FORESTER
  90                   $JAVA
  91                   $NODE_LIST
  92                   $RIO_SLAVE_DRIVER
  93                   $RIO_SLAVE
  94                   $TEMP_DIR_DEFAULT
  95                   $EXPASY_SPROT_SEARCH_DE
  96                   $EXPASY_SPROT_SEARCH_AC
  97  );
  98
  99
 100
 101
 102 # =============================================================================
 103 # =============================================================================
 104 #
 105 # THESE VARIABLES ARE ENVIRONMENT DEPENDENT, AND NEED TO BE SET ACCORDINGLY
 106 # BY THE USER
 107 # -------------------------------------------------------------------------
 108 #
 109
 110
 111
 112 # RIO itself:
 113 # -----------
 114 our $PATH_TO_FORESTER          = "/nfs/dm3/homedir1/czmasek/RIO1.24/";
 115
 116
 117 # Java virtual machine:
 118 # ---------------------
 119 our $JAVA                      = "/usr/local/java/jdk/bin/java";
 120
 121
 122
 123 # Where all the temporary files can be created:
 124 # ---------------------------------------------
 125 our $TEMP_DIR_DEFAULT          = "/tmp/";
 126
 127
 128
 129 # Pfam data:
 130 # ----------
 131 our $PFAM_FULL_DIRECTORY       = "/path/to/Pfam/Full/";
 132 our $PFAM_SEED_DIRECTORY       = "/path/to/Pfam/Seed/";
 133 our $PFAM_HMM_DB               = "/path/to/Pfam/Pfam_ls"; # Need to run "hmmindex" on this
 134                                                           # to produce .ssi file.
 135                                                           # Then, for example
 136                                                           # "setenv HMMERDB /home/rio/pfam-6.6/"
 137
 138
 139 $PATH_TO_FORESTER = &addSlashAtEndIfNotPresent( $PATH_TO_FORESTER );
 140
 141
 142 # Description lines and species from SWISS-PROT and TrEMBL:
 143 # ---------------------------------------------------------
 144 our $TREMBL_ACDEOS_FILE        = $PATH_TO_FORESTER."data/trembl22_ACDEOS_1-6";
 145
 146 our $SWISSPROT_ACDEOS_FILE     = $PATH_TO_FORESTER."data/sp40_ACDEOS_1-6";
 147
 148
 149
 150 # Names of species which can be analyzed and analyzed
 151 # against (must also be in tree $SPECIES_TREE_FILE_DEFAULT).
 152 # By using a list with less species, RIO analyses become faster
 153 # but lose phylogenetic resolution.
 154 # For many purposes, list "tree_of_life_bin_1-6_species_list"
 155 # in "data/species/" might be sufficient:
 156 # --------------------------------------------------------------
 157 our $SPECIES_NAMES_FILE        = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6_species_list";
 158
 159
 160
 161 # A default species tree in NHX format.
 162 # For many purposes, tree "tree_of_life_bin_1-6.nhx"
 163 # in "data/species/" might be fine:
 164 # --------------------------------------------------
 165 our $SPECIES_TREE_FILE_DEFAULT = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6.nhx";
 166
 167
 168
 169 # Data for using precalculated distances:
 170 # ---------------------------------------
 171 our $MATRIX_FOR_PWD            = 2;  # The matrix which has been used for the pwd in $RIO_PWD_DIRECTORY.
 172                                      # 0=JTT, 1=PAM, 2=BLOSUM 62, 3=mtREV24, 5=VT, 6=WAG.
 173
 174 our $RIO_PWD_DIRECTORY         = $PATH_TO_FORESTER."example_data/";  # all must end with "/"
 175 our $RIO_BSP_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 176 our $RIO_NBD_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 177 our $RIO_ALN_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 178 our $RIO_HMM_DIRECTORY         = $PATH_TO_FORESTER."example_data/";
 179
 180
 181
 182 #
 183 # End of variables which need to be set by the user.
 184 #
 185 # =============================================================================
 186 # =============================================================================
 187
 188
 189
 190
 191
 192 $TEMP_DIR_DEFAULT    = &addSlashAtEndIfNotPresent( $TEMP_DIR_DEFAULT );
 193 $PFAM_FULL_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_FULL_DIRECTORY );
 194 $PFAM_SEED_DIRECTORY = &addSlashAtEndIfNotPresent( $PFAM_SEED_DIRECTORY );
 195
 196
 197
 198 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 199 # These variables should normally not be changed:
 200 #
 201
 202 our $PRIOR_FILE_DIR            = $PATH_TO_FORESTER."data/priors_for_hmmbuild/";
 203                                  # Directory containing dirichlet prior
 204                                  # files needed for certain aligments
 205                                  # by hmmbuild (e.g. Collagen).
 206
 207
 208 # PHYLIP:
 209 our $SEQBOOT                   = $PATH_TO_FORESTER."phylip_mod/exe/seqboot";
 210 our $NEIGHBOR                  = $PATH_TO_FORESTER."phylip_mod/exe/neighbor";
 211 our $PROTPARS                  = $PATH_TO_FORESTER."phylip_mod/exe/protpars";
 212 our $CONSENSE                  = $PATH_TO_FORESTER."phylip_mod/exe/consense";
 213
 214 # TREE-PUZZLE:
 215 our $PUZZLE                    = $PATH_TO_FORESTER."puzzle_mod/src/puzzle";
 216 our $PUZZLE_DQO                = $PATH_TO_FORESTER."puzzle_dqo/src/puzzle";
 217
 218 # HMMER:
 219 our $HMMALIGN                  = $PATH_TO_FORESTER."hmmer/binaries/hmmalign";
 220 our $HMMSEARCH                 = $PATH_TO_FORESTER."hmmer/binaries/hmmsearch";
 221 our $HMMBUILD                  = $PATH_TO_FORESTER."hmmer/binaries/hmmbuild";
 222 our $HMMFETCH                  = $PATH_TO_FORESTER."hmmer/binaries/hmmfetch";
 223 our $SFE                       = $PATH_TO_FORESTER."hmmer/binaries/sfetch";
 224 our $HMMCALIBRATE              = $PATH_TO_FORESTER."hmmer/binaries/hmmcalibrate";
 225
 226 our $P7EXTRACT                 = $PATH_TO_FORESTER."perl/p7extract.pl";
 227 our $MULTIFETCH                = $PATH_TO_FORESTER."perl/multifetch.pl";
 228
 229
 230 # RIO/FORESTER:
 231 our $BOOTSTRAP_CZ              = $PATH_TO_FORESTER."C/bootstrap_cz";
 232 our $BOOTSTRAP_CZ_PL           = $PATH_TO_FORESTER."perl/bootstrap_cz.pl";
 233 our $TRANSFERSBRANCHLENGHTS    = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.transfersBranchLenghts";
 234 our $MAKETREE                  = $PATH_TO_FORESTER."perl/makeTree.pl";
 235 our $RIO_PL                    = $PATH_TO_FORESTER."perl/rio.pl";
 236 our $DORIO                     = $JAVA." -cp $PATH_TO_FORESTER"."java forester.tools.DoRIO";
 237 # parallel RIO:
 238 our $RIO_SLAVE_DRIVER          = $PATH_TO_FORESTER."perl/rio_slave_driver.pl";
 239 our $RIO_SLAVE                 = $PATH_TO_FORESTER."perl/rio_slave.pl";
 240 our $NODE_LIST                 = $PATH_TO_FORESTER."data/node_list.dat";
 241
 242 #
 243 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 244
 245
 246 our $BOOTSTRAPS         = 100;
 247 our $MIN_NUMBER_OF_AA   = 20;  # After removal of gaps, if less, gaps are not removed.
 248 our $LENGTH_OF_NAME     = 26;
 249
 250
 251
 252
 253 our $MULTIPLE_TREES_FILE_SUFFIX  = ".mlt";
 254 our $LOG_FILE_SUFFIX             = ".log";
 255 our $ALIGN_FILE_SUFFIX           = ".aln";
 256 our $TREE_FILE_SUFFIX            = ".nhx";
 257 our $ADDITION_FOR_RIO_ANNOT_TREE = ".rio";
 258 our $SUFFIX_PWD                  = ".pwd";
 259 our $SUFFIX_BOOT_STRP_POS        = ".bsp";
 260 our $SUFFIX_PWD_NOT_BOOTS        = ".nbd";
 261 our $SUFFIX_HMM                  = ".hmm";
 262
 263 our $EXPASY_SPROT_SEARCH_DE      = "http://www.expasy.org/cgi-bin/sprot-search-de?";
 264 our $EXPASY_SPROT_SEARCH_AC      = "http://www.expasy.org/cgi-bin/sprot-search-ac?";
 265
 266
 267
 268 # One argument: input multiple trees file
 269 # Last modified: 07/05/01
 270 sub executeConsense {
 271     my $in  = $_[ 0 ];
 272
 273     &testForTextFilePresence( "$in" );
 274
 275     system( "$CONSENSE >/dev/null 2>&1 << !
 276 $in
 277 Y
 278 !" )
 279     && &dieWithUnexpectedError( "Could not execute \"$CONSENSE \"" );
 280
 281     return;
 282 }
 283
 284
 285
 286 # Four arguments:
 287 # 1. options ("-" is not necessary)
 288 # 2. alignment or pwd file
 289 # 3. outfile
 290 # 4. temp dir
 291 # Last modified: 07/05/01
 292 sub executeMakeTree {
 293
 294     my $opts         = $_[ 0 ];
 295     my $B            = $_[ 1 ];
 296     my $C            = $_[ 2 ];
 297     my $D            = $_[ 3 ];
 298
 299     &testForTextFilePresence( $B );
 300
 301     $opts = "-".$opts;
 302
 303     system( "$MAKETREE $opts $B $C $D" )
 304     && &dieWithUnexpectedError( "Could not execute \"$MAKETREE $opts $B $C $D\"" );
 305
 306 } ## executeMakeTree
 307
 308
 309
 310
 311 # Two arguments:
 312 # 1. Name of inputfile
 313 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 314 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 315 sub executePuzzleDQO {
 316     my $in            = $_[ 0 ];
 317     my $matrix_option = $_[ 1 ];
 318     my $mat           = "";
 319
 320     &testForTextFilePresence( $in );
 321
 322     $mat = setModelForPuzzle( $matrix_option );
 323
 324     system( "$PUZZLE_DQO $in >/dev/null 2>&1 << !$mat
 325 y
 326 !" )
 327     && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO\"" );
 328
 329     return;
 330
 331 } ## executePuzzleDQO
 332
 333
 334
 335
 336 # Two arguments:
 337 # 1. Name of inputfile
 338 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 339 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 340 # Last modified: 01/28/02
 341 sub executePuzzleDQObootstrapped {
 342     my $in            = $_[ 0 ];
 343     my $matrix_option = $_[ 1 ];
 344
 345
 346     my $l             = 0;
 347     my $slen          = 0;
 348     my $counter       = 0;
 349     my $mat           = "";
 350     my $a             = "";
 351     my @a             = ();
 352
 353     &testForTextFilePresence( $in );
 354
 355     open( GRP, "<$in" ) || &dieWithUnexpectedError( "Cannot open file \"$in\"" );
 356     while( <GRP> ) {
 357         if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
 358             $counter++;
 359         }
 360     }
 361     close( GRP );
 362
 363     $l   = `cat $in | wc -l`;
 364     $slen   = $l / $counter;
 365
 366     system( "split -$slen $in $in.splt." )
 367     && &dieWithUnexpectedError( "Could not execute \"split -$slen $in $in.splt.\"" );
 368
 369     @a = <$in.splt.*>;
 370
 371     $mat = setModelForPuzzle( $matrix_option );
 372
 373     foreach $a ( @a ) {
 374
 375         system( "$PUZZLE_DQO $a >/dev/null 2>&1 << !$mat
 376 y
 377 !" )
 378         && &dieWithUnexpectedError( "Could not execute \"$PUZZLE_DQO $a\"" );
 379
 380         system( "cat $a.dist >> $in.dist" )
 381         && &dieWithUnexpectedError( "Could not execute \"cat outdist >> $in.dist\"" );
 382
 383         unlink( $a, $a.".dist" );
 384     }
 385
 386     return;
 387
 388 } ## executePuzzleDQObootstrapped
 389
 390
 391
 392 # Transfers a Pfam (SELEX) alignment to a
 393 # PHYLIP sequential style alignment.
 394 # It only writes "match columns" as indicated by the
 395 # "# RF" line ('x' means match).
 396 #
 397 # Three arguments:
 398 # 1. infile name
 399 # 2. outfile name
 400 # 3. 1 to NOT ensure that match states contain only 'A'-'Z' or '-'
 401 #
 402 # Returns the number of match states (=length of output alignment),
 403 #         the length of the input alignment,
 404 #         the number of seqs in the input alignment
 405 #
 406 # Last modified: 07/07/01
 407 #
 408 sub pfam2phylipMatchOnly {
 409
 410     my $infile          = $_[ 0 ];
 411     my $outfile         = $_[ 1 ];
 412     my $ne              = $_[ 2 ];
 413     my @seq_name        = ();
 414     my @seq_array       = ();
 415     my $return_line     = "";
 416     my $seq             = "";
 417     my $x               = 0;
 418     my $y               = 0;
 419     my $i               = 0;
 420     my $x_offset        = 0;
 421     my $max_x           = 0;
 422     my $rf_y            = 0;
 423     my $number_colum    = 0;
 424     my $not_ensure      = 0;
 425     my $saw_rf_line     = 0;
 426
 427     if ( $ne && $ne == 1 ) {
 428         $not_ensure = 1;
 429     }
 430
 431     &testForTextFilePresence( $infile );
 432
 433     open( INPP, "$infile" ) || &dieWithUnexpectedError( "Cannot open file \"$infile\"" );
 434
 435     # This reads in the first block. It reads in the seq names.
 436     while ( 1 ) {
 437         if ( &isPfamSequenceLine( $return_line ) ) {
 438             $return_line =~ /^(\S+)\s+(\S+)/;
 439             $seq_name[ $y ] = substr( $1, 0, $LENGTH_OF_NAME );
 440             $seq = $2;
 441             for ( $x = 0; $x < length( $seq ); $x++ ) {
 442                 $seq_array[ $x ][ $y ] = substr( $seq, $x, 1 );
 443             }
 444             $y++;
 445         }
 446         elsif ( &isRFline( $return_line ) ) {
 447             $saw_rf_line = 1;
 448             $return_line =~ /\s+(\S+)\s*$/;
 449             $seq = $1;
 450             $x_offset = length( $seq );
 451             $rf_y = $y;
 452             for ( $x = 0; $x < $x_offset; $x++ ) {
 453                 $seq_array[ $x ][ $rf_y ] = substr( $seq, $x, 1 );
 454             }
 455             last;
 456         }
 457
 458         $return_line = <INPP>;
 459
 460         if ( !$return_line ) {
 461              &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
 462         }
 463     }
 464
 465     if ( $saw_rf_line != 1 ) {
 466          &dieWithUnexpectedError( "Alignment not in expected format (no RF line)" );
 467     }
 468
 469     $y = 0;
 470     $max_x = 0;
 471
 472     # This reads all blocks after the 1st one.
 473     while ( $return_line = <INPP> ) {
 474         if ( &isPfamSequenceLine( $return_line ) ) {
 475             $return_line =~ /^\S+\s+(\S+)/;
 476             $seq = $1;
 477             for ( $x = 0; $x < length( $seq ); $x++ ) {
 478                 $seq_array[ $x + $x_offset ][ $y % $rf_y ] = substr( $seq, $x, 1 );
 479             }
 480             $y++;
 481         }
 482         elsif ( &isRFline( $return_line ) ) {
 483             if ( $y != $rf_y ) {
 484                 &dieWithUnexpectedError( "Alignment not in expected format" );
 485             }
 486
 487             $return_line =~ /\s+(\S+)\s*$/;
 488             $seq = $1;
 489             $max_x = length( $seq );
 490
 491             for ( $x = 0; $x < length( $seq ); $x++ ) {
 492                 $seq_array[ $x + $x_offset ][ $rf_y ] = substr( $seq, $x, 1 );
 493             }
 494
 495             $y = 0;
 496             $x_offset = $x_offset + $max_x;
 497             $max_x = 0;
 498         }
 499     }
 500
 501     close( INPP );
 502
 503     # Counts the match states, and hence the number of aa in the alignment:
 504     for ( $x = 0; $x < $x_offset; $x++ ) {
 505         if ( !$seq_array[ $x ][ $rf_y ] ) {
 506             &dieWithUnexpectedError( "Alignment not in expected format" );
 507         }
 508         if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
 509             $number_colum++;
 510         }
 511     }
 512
 513     # Writes the file:
 514
 515     open( OUTPP, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
 516     print OUTPP "$rf_y $number_colum\n";
 517     for ( $y = 0; $y < $rf_y; $y++ ) {
 518         print OUTPP "$seq_name[ $y ]";
 519         for ( $i = 0; $i < ( $LENGTH_OF_NAME - length( $seq_name[ $y ] ) ); $i++ ) {
 520                 print OUTPP " ";
 521         }
 522         for ( $x = 0; $x < $x_offset; $x++ ) {
 523             if ( $seq_array[ $x ][ $rf_y ] eq 'x' ) {
 524                 if ( !$seq_array[ $x ][ $y ] ) {
 525                     &dieWithUnexpectedError( "Alignment not in expected format" );
 526                 }
 527                 if ( $not_ensure != 1 && $seq_array[ $x ][ $y ] !~ /[A-Z]|-/ ) {
 528                     &dieWithUnexpectedError( "Alignment not in expected format (match states must only contain 'A'-'Z' or '-')" );
 529                 }
 530                 print OUTPP "$seq_array[ $x ][ $y ]";
 531             }
 532         }
 533         print OUTPP "\n";
 534     }
 535     close( OUTPP );
 536
 537     return $number_colum, $x_offset, $rf_y;
 538
 539 } ## pfam2phylipMatchOnly
 540
 541
 542
 543 # Returns whether the argument (a String)
 544 # starts with a SWISS-PROT name (SEQN_SPECI).
 545 # Last modified: 06/21/01
 546 sub startsWithSWISS_PROTname {
 547     return ( $_[ 0 ] =~ /^[A-Z0-9]{1,4}_[A-Z0-9]{1,5}/ );
 548 }
 549
 550
 551
 552 # Returns whether the argument starts with XXX.. XXXXX.. and the first
 553 # character is not a "#".
 554 # Last modified: 06/21/01
 555 sub isPfamSequenceLine {
 556     return( !&isPfamCommentLine( $_[ 0 ] )
 557     && &containsPfamNamedSequence( $_[ 0 ] ) );
 558 }
 559
 560
 561
 562 # Returns whether the argument does start with a "#".
 563 # Last modified: 06/21/01
 564 sub isPfamCommentLine {
 565     return ( $_[ 0 ] =~ /^#/ );
 566 }
 567
 568
 569
 570 # Returns whether the argument starts with XXX XXXXX.
 571 # Last modified: 06/21/01
 572 sub containsPfamNamedSequence {
 573     return ( $_[ 0 ] =~ /^\S+\s+\S+/ );
 574 }
 575
 576
 577 # Returns whether the argument starts with XXX XXXXX.
 578 # Last modified: 06/21/01
 579 sub isRFline {
 580     return ( $_[ 0 ] =~ /^#.*RF/ );
 581 }
 582
 583
 584
 585
 586 # Five arguments:
 587 # 1. pairwise distance file
 588 # 2. number of bootstraps
 589 # 3. randomize_input_order: 0: do not randomize input order; >=1 jumble
 590 # 4. seed for random number generator
 591 # 5. lower-triangular data matrix? 1: yes; no, otherwise
 592 # Last modified: 06/08/01
 593 sub executeNeighbor {
 594     my $inpwd  = $_[ 0 ];
 595     my $bs     = $_[ 1 ];
 596     my $rand   = $_[ 2 ];
 597     my $s      = $_[ 3 ];
 598     my $l      = $_[ 4 ];
 599     my $jumble = "";
 600     my $multi  = "";
 601     my $lower  =  "";
 602
 603
 604     &testForTextFilePresence( $inpwd );
 605
 606     if ( $rand >= 1 ) {
 607         $jumble = "
 608 J
 609 $s";
 610     }
 611
 612     if (  $bs >= 2 ) {
 613         $multi = "
 614 M
 615 $bs
 616 $s";
 617     }
 618     if ( $l == 1 ) {
 619         $lower = "
 620 L";
 621     }
 622
 623
 624     system( "$NEIGHBOR >/dev/null 2>&1 << !
 625 $inpwd$jumble$multi$lower
 626 2
 627 3
 628 Y
 629 !" )
 630     && &dieWithUnexpectedError( "Could not execute \"$NEIGHBOR $inpwd$jumble$multi$lower\"" );
 631     # 3: Do NOT print out tree
 632
 633
 634     return;
 635
 636 } ## executeNeighbor
 637
 638
 639
 640 # Four arguments:
 641 # 1. name of alignment file (in correct format!)
 642 # 2. number of bootstraps
 643 # 3. jumbles: 0: do not jumble; >=1 number of jumbles
 644 # 4. seed for random number generator
 645 # Last modified: 03/13/04
 646 sub executeProtpars {
 647     my $alin   = $_[ 0 ];
 648     my $bs     = $_[ 1 ];
 649     my $rand   = $_[ 2 ];
 650     my $s      = $_[ 3 ];
 651     my $jumble = "";
 652     my $multi  = "";
 653
 654
 655     &testForTextFilePresence( $alin );
 656
 657     if ( $bs >= 2 && $rand < 1 ) {
 658          $rand = 1;
 659     }
 660
 661     if ( $rand >= 1 ) {
 662         $jumble = "
 663 J
 664 $s
 665 $rand";
 666     }
 667
 668     if (  $bs >= 2 ) {
 669         $multi = "
 670 M
 671 D
 672 $bs";
 673     }
 674
 675
 676
 677     system( "$PROTPARS  2>&1 << !
 678 $alin$jumble$multi
 679 I
 680 3
 681 Y
 682 !" )
 683     && &dieWithUnexpectedError( "Could not execute \"$PROTPARS $alin$jumble$multi\"" );
 684     # 3: Do NOT print out tree
 685     # I: Interleaved
 686
 687     return;
 688
 689 } ## executeProtpars
 690
 691
 692
 693 # "Model of substitution" order for DQO TREE-PUZZLE 5.0:
 694 # Auto
 695 # m -> Dayhoff (Dayhoff et al. 1978)
 696 # m -> JTT (Jones et al. 1992)
 697 # m -> mtREV24 (Adachi-Hasegawa 1996)
 698 # m -> BLOSUM62 (Henikoff-Henikoff 92)
 699 # m -> VT (Mueller-Vingron 2000)
 700 # m -> WAG (Whelan-Goldman 2000)
 701 # m -> Auto
 702 # One argument:
 703 # matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 704 # 5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 705 # Last modified: 07/07/01
 706 sub setModelForPuzzle {
 707     my $matrix_option = $_[ 0 ];
 708     my $matr          = "";
 709
 710     if ( $matrix_option == 0 ) { # JTT
 711         $matr = "
 712 m
 713 m";
 714     }
 715     elsif ( $matrix_option == 2 ) { # BLOSUM 62
 716         $matr = "
 717 m
 718 m
 719 m
 720 m";
 721     }
 722     elsif ( $matrix_option == 3 ) { # mtREV24
 723         $matr = "
 724 m
 725 m
 726 m";
 727     }
 728     elsif ( $matrix_option == 5 ) { # VT
 729         $matr = "
 730 m
 731 m
 732 m
 733 m
 734 m";
 735     }
 736     elsif ( $matrix_option == 6 ) { # WAG
 737         $matr = "
 738 m
 739 m
 740 m
 741 m
 742 m
 743 m";
 744     }
 745     elsif ( $matrix_option == 7 ) { # auto
 746         $matr = "";
 747     }
 748     else { # PAM
 749         $matr = "
 750 m"
 751     }
 752
 753     return $matr;
 754
 755 } ## setModelForPuzzle
 756
 757 # One argument:
 758 # Model of rate heterogeneity:
 759 #    1 for "8 Gamma distributed rates"
 760 #    2 for "Two rates (1 invariable + 1 variable)"
 761 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
 762 #    otherwise: Uniform rate
 763 # Last modified: 09/08/03
 764 sub setRateHeterogeneityOptionForPuzzle {
 765     my $rate_heterogeneity_option  = $_[ 0 ];
 766     my $opt                        = "";
 767
 768     if ( $rate_heterogeneity_option == 1 ) {
 769         $opt = "
 770 w";
 771     }
 772     elsif ( $rate_heterogeneity_option == 2 ) {
 773         $opt = "
 774 w
 775 w";
 776     }
 777     elsif ( $rate_heterogeneity_option == 3 ) {
 778         $opt = "
 779 w
 780 w
 781 w";
 782     }
 783     else {
 784         $opt = "";
 785     }
 786
 787     return $opt;
 788 } ## setRateHeterogeneityOptionForPuzzle
 789
 790
 791 # One argument:
 792 # Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
 793 # Last modified: 09/08/03
 794 sub setParameterEstimatesOptionForPuzzle {
 795     my $parameter_estimates_option  = $_[ 0 ];
 796     my $opt                         = "";
 797
 798     if ( $parameter_estimates_option == 1 ) {
 799         $opt = "
 800 e";
 801     }
 802     else {
 803         $opt = "";
 804     }
 805
 806     return $opt;
 807 } ## setParameterEstimatesOptionForPuzzle
 808
 809
 810
 811 # Two/three/four arguments:
 812 # 1. Name of inputfile
 813 # 2. matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 814 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 815 # 3. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
 816 # 4. Model of rate heterogeneity:
 817 #    1 for "8 Gamma distributed rates"
 818 #    2 for "Two rates (1 invariable + 1 variable)"
 819 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
 820 #    otherwise: Uniform rate
 821 # Last modified: 09/08/03 (added 3rd and 4th parameter)
 822 sub executePuzzleBootstrapped {
 823     my $in                         = $_[ 0 ];
 824     my $matrix_option              = $_[ 1 ];
 825     my $parameter_estimates_option = $_[ 2 ];
 826     my $rate_heterogeneity_option  = $_[ 3 ];
 827
 828     my $l             = 0;
 829     my $slen          = 0;
 830     my $counter       = 0;
 831     my $mat           = "";
 832     my $est           = "";
 833     my $rate          = "";
 834     my $a             = "";
 835     my @a             = ();
 836
 837     &testForTextFilePresence( $in );
 838
 839     open( GRP, "<$in" ) || die "\n\n$0: Unexpected error: Cannot open file <<$in>>: $!";
 840     while( <GRP> ) {
 841         if ( $_ =~ /^\s*\d+\s+\d+\s*$/ ) {
 842             $counter++;
 843         }
 844     }
 845     close( GRP );
 846
 847         $l   = `cat $in | wc -l`;
 848         $slen   = $l / $counter;
 849
 850         system( "split -$slen $in $in.splt." )
 851     && die "\n\n$0: executePuzzleDQObootstrapped: Could not execute \"split -$slen $in $in.splt.\": $!";
 852
 853     @a = <$in.splt.*>;
 854
 855     $mat = setModelForPuzzle( $matrix_option );
 856      if ( $parameter_estimates_option ) {
 857         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
 858     }
 859     if ( $rate_heterogeneity_option ) {
 860         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
 861     }
 862
 863     foreach $a ( @a ) {
 864         print "-".$a."\n";
 865         system( "$PUZZLE $a << !
 866 k
 867 k$mat$est$rate
 868 y
 869 !" )
 870         && die "$0: Could not execute \"$PUZZLE $a\"";
 871
 872         system( "cat $a.dist >> $in.dist" )
 873         && die "$0: Could not execute \"cat outdist >> $in.dist\"";
 874
 875         unlink( $a, $a.".dist", $a.".tree" );
 876     }
 877
 878     return;
 879
 880 } ## executePuzzleBootstrapped
 881
 882
 883
 884
 885
 886 # Two/three/four arguments:
 887 # 1. Name of inputfile
 888 # 2. Matrix option: 0 = JTT; 2 = BLOSUM 62; 3 = mtREV24;
 889 #    5 = VT; 6 = WAG; 7 = auto; PAM otherwise
 890 # 3. Parameter estimates: 1 for "Exact (slow)"; "Approximate (faster)" otherwise
 891 # 4. Model of rate heterogeneity:
 892 #    1 for "8 Gamma distributed rates"
 893 #    2 for "Two rates (1 invariable + 1 variable)"
 894 #    3 for "Mixed (1 invariable + 8 Gamma rates)"
 895 #    otherwise: Uniform rate
 896 # Last modified: 09/08/03 (added 3rd and 4th parameter)
 897 sub executePuzzle {
 898     my $in                         = $_[ 0 ];
 899     my $matrix_option              = $_[ 1 ];
 900     my $parameter_estimates_option = $_[ 2 ];
 901     my $rate_heterogeneity_option  = $_[ 3 ];
 902     my $mat                        = "";
 903     my $est                        = "";
 904     my $rate                       = "";
 905
 906     &testForTextFilePresence( $in );
 907
 908     $mat = &setModelForPuzzle( $matrix_option );
 909     if ( $parameter_estimates_option ) {
 910         $est = &setParameterEstimatesOptionForPuzzle( $parameter_estimates_option );
 911     }
 912     if ( $rate_heterogeneity_option ) {
 913         $rate = &setRateHeterogeneityOptionForPuzzle( $rate_heterogeneity_option );
 914     }
 915
 916
 917     system( "$PUZZLE $in << !
 918 k
 919 k$mat$est$rate
 920 y
 921 !" )
 922     && die "$0: Could not execute \"$PUZZLE\"";
 923
 924     return;
 925
 926 } ## executePuzzle
 927
 928
 929
 930
 931 # Preparation of the pwd file
 932 sub addDistsToQueryToPWDfile {
 933     my $pwd_file          = $_[ 0 ];
 934     my $disttoquery_file  = $_[ 1 ];
 935     my $outfile           = $_[ 2 ];
 936     my $name_of_query     = $_[ 3 ];
 937     my $name_of_query_    = "";
 938     my $return_line_pwd   = "";
 939     my $return_line_dq    = "";
 940     my $num_of_sqs        = 0;
 941     my $block             = 0;
 942     my $name_from_pwd     = "X";
 943     my $name_from_dq      = "Y";
 944     my @dists_to_query    = ();
 945     my $i                 = 0;
 946
 947     &testForTextFilePresence( $pwd_file );
 948     &testForTextFilePresence( $disttoquery_file );
 949
 950     $name_of_query_ = $name_of_query;
 951     for ( my $j = 0; $j <= ( $LENGTH_OF_NAME - length( $name_of_query ) - 1 ); ++$j ) {
 952             $name_of_query_ .= " ";
 953     }
 954
 955     open( OUT_AD, ">$outfile" ) || &dieWithUnexpectedError( "Cannot create file \"$outfile\"" );
 956     open( IN_PWD, "$pwd_file" ) || &dieWithUnexpectedError( "Cannot open file \"$pwd_file\"" );
 957     open( IN_DQ, "$disttoquery_file" ) || &dieWithUnexpectedError( "Cannot open file \"$disttoquery_file\"" );
 958
 959     W: while ( $return_line_pwd = <IN_PWD> ) {
 960
 961
 962         if ( $return_line_pwd =~ /^\s*(\d+)\s*$/ ) {
 963             $num_of_sqs = $1;
 964             $num_of_sqs++;
 965             if ( $block > 0 ) {
 966                 print OUT_AD "$name_of_query_  ";
 967                 for ( my $j = 0; $j < $i; ++$j ) {
 968                     print OUT_AD "$dists_to_query[ $j ]  ";
 969                 }
 970                 print OUT_AD "0.0\n";
 971             }
 972             print OUT_AD "  $num_of_sqs\n";
 973             $block++;
 974             @dists_to_query = ();
 975             $i = 0;
 976         }
 977
 978         if ( $block == 1
 979         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
 980             $name_from_pwd = $1;
 981
 982             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
 983                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
 984             }
 985
 986             if ( $return_line_dq !~ /\S/ ) {
 987                 if ( !defined( $return_line_dq = <IN_DQ> ) ) {
 988                     &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
 989                 }
 990             }
 991             $return_line_dq =~ /^\s*(\S+)\s+(\S+)/;
 992             $name_from_dq = $1;
 993             $dists_to_query[ $i++ ] = $2;
 994
 995
 996             if ( $name_from_pwd ne $name_from_dq ) {
 997                 &dieWithUnexpectedError( "Order of sequence names in \"$pwd_file\" and \"$disttoquery_file\" is not the same" );
 998             }
 999             print OUT_AD $return_line_pwd;
1000
1001         }
1002         elsif ( $block > 1
1003         && $return_line_pwd =~ /^\s*(\S+)\s+\S+/ ) {
1004             $name_from_pwd = $1;
1005             if ( !defined( $return_line_dq = <IN_DQ> ) ) {
1006                 &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1007             }
1008             if ( $return_line_dq !~ /\S/ ) {
1009                 if ( !defined( $return_line_dq = <IN_DQ>) ) {
1010                     &dieWithUnexpectedError( "\"$disttoquery_file\" seems too short" );
1011                 }
1012             }
1013             $return_line_dq =~ /^\s*\S+\s+(\S+)/;
1014             $dists_to_query[ $i++ ] = $1;
1015             print OUT_AD $return_line_pwd;
1016         }
1017     }
1018     print OUT_AD "$name_of_query_  ";
1019     for ( my $j = 0; $j < $i; ++$j ) {
1020         print OUT_AD "$dists_to_query[ $j ]  ";
1021     }
1022     print OUT_AD "0.0\n";
1023
1024     close( OUT_AD );
1025     close( IN_PWD );
1026     close( IN_DQ );
1027     return $block;
1028
1029 } ## addDistsToQueryToPWDfile
1030
1031
1032
1033
1034 # Three arguments:
1035 # 1. HMMER model db
1036 # 2. name of HMM
1037 # 3. outputfile name
1038 # Last modified: 02/27/01
1039 sub executeHmmfetch {
1040
1041     my $db      = $_[ 0 ];
1042     my $name    = $_[ 1 ];
1043     my $outfile = $_[ 2 ];
1044
1045     system( "$HMMFETCH $db $name > $outfile" )
1046     && &dieWithUnexpectedError( "Could not execute \"$HMMFETCH $db $name > $outfile\"" );
1047     return;
1048
1049 } ## executeHmmfetch
1050
1051
1052
1053 # Checks wether a file is present, not empty and a plain textfile.
1054 # One argument: name of file.
1055 # Last modified: 07/07/01
1056 sub testForTextFilePresence {
1057     my $file = $_[ 0 ];
1058     unless ( ( -s $file ) && ( -f $file ) && ( -T $file ) ) {
1059         dieWithUnexpectedError( "File \"$file\" does not exist, is empty, or is not a plain textfile" );
1060     }
1061 } ## testForTextFilePresence
1062
1063
1064 # Last modified: 02/21/03
1065 sub addSlashAtEndIfNotPresent {
1066     my $filename = $_[ 0 ];
1067     $filename =~ s/\s+//g;
1068     unless ( $filename =~ /\/$/ ) {
1069        $filename = $filename."/";
1070     }
1071     return $filename;
1072 } ## addSlashAtEndIfNotPresent
1073
1074
1075
1076 # Last modified: 02/15/02
1077 sub exitWithWarning {
1078
1079     my $text = $_[ 0 ];
1080     if ( defined( $_[ 1 ] ) && $_[ 1 ] == 1 ) {
1081         print( "<H4 class=\"error\">user error</H4>\n" );
1082         print( "<P>\n" );
1083         print( "<B>$text</B>\n" );
1084         print( "</P>\n" );
1085         print( "<P> &nbsp</P>\n" );
1086     }
1087     else {
1088         print( "\n\n$text\n\n" );
1089     }
1090
1091     exit( 0 );
1092
1093 } ## exit_with_warning
1094
1095
1096
1097 # Last modified: 02/15/02
1098 sub dieWithUnexpectedError {
1099
1100     my $text = $_[ 0 ];
1101
1102     die( "\n\n$0:\nUnexpected error (should not have happened):\n$text\n$!\n\n" );
1103
1104 } ## dieWithUnexpectedError
1105
1106
1107
1108 1;