forester/archive/perl/pfam2pwd.pl

   1 #!/usr/bin/perl -W
   2
   3 # pfam2pwd.pl
   4 # -----------
   5 # Copyright (C) 1999-2002 Washington University School of Medicine
   6 # and Howard Hughes Medical Institute
   7 # All rights reserved
   8 #
   9 # Author: Christian M. Zmasek
  10 #         zmasek@genetics.wustl.edu
  11 #         http://www.genetics.wustl.edu/eddy/people/zmasek/
  12 #
  13 # Created: 05/17/01
  14 #
  15 # Last modified 02/20/03
  16 #
  17 #
  18 #  See RIO_INSTALL on how to use this program.
  19 #  ------------------------------------------
  20 #
  21
  22 use strict;
  23
  24 use FindBin;
  25 use lib $FindBin::Bin;
  26 use rio_module;
  27
  28 my $VERSION = "3.002";
  29
  30
  31
  32 # =============================================================================
  33 # =============================================================================
  34 #
  35 # THESE VARIABLES NEED TO BE SET BY THE USER
  36 # ------------------------------------------
  37 #
  38
  39
  40 # Pfam alignments to calculate pairwise distances from:
  41 # -----------------------------------------------------
  42 my $MY_PFAM_FULL_DIRECTORY = "/path/to/Pfam/Full/"; # must end with "/"
  43
  44
  45
  46 # This file lists all the alignments for which to calculate pairwise distances
  47 # from. If left empty, ALL the alignments in $MY_PFAM_FULL_DIRECTORY
  48 # will be used:
  49 # ----------------------------------------------------------------------------
  50 my $ALGNS_TO_USE_LIST_FILE = "";
  51
  52
  53
  54 # This is _VERY IMPORTANT_. It determines the species whose sequences
  55 # are being used (sequences from species not listed in $MY_SPECIES_NAMES_FILE
  56 # are ignored). Normally, one would use the same list as RIO uses
  57 # ($SPECIES_NAMES_FILE in "rio_module.pm") -- currently "tree_of_life_bin_1-6.nhx".
  58 #
  59 # For certain large families (such  as protein kinases), one might use a
  60 # species file which contains less species in order to be able to finish
  61 # the calculations in reasonable time.
  62 # For example, to exclude most mammals, use:
  63 # my $MY_SPECIES_NAMES_FILE = $PATH_TO_FORESTER."data/species/tree_of_life_bin_1-6_species_list_NO_RAT_MONKEYS_APES_SHEEP_GOAT_HAMSTER"
  64 # (to only use sequences from SWISS-PROT add this line:
  65 # $TREMBL_ACDEOS_FILE = $PATH_TO_FORESTER."data/NO_TREMBL";)
  66 # ----------------------------------------------------------------------------
  67 my $MY_SPECIES_NAMES_FILE  = $SPECIES_NAMES_FILE;
  68
  69
  70
  71 # This is were the output goes (must end with "/")
  72 # ------------------------------------------------
  73 my $MY_RIO_PWD_DIRECTORY   = "/path/to/pfam2pwd_out/pwd/";
  74 my $MY_RIO_BSP_DIRECTORY   = "/path/to/pfam2pwd_out/bsp/";
  75 my $MY_RIO_NBD_DIRECTORY   = "/path/to/pfam2pwd_out/nbd/";
  76 my $MY_RIO_ALN_DIRECTORY   = "/path/to/pfam2pwd_out/aln/";
  77 my $MY_RIO_HMM_DIRECTORY   = "/path/to/pfam2pwd_out/hmm/";
  78
  79
  80
  81 # A directory to create temporary files in:
  82 # -----------------------------------------
  83 my $MY_TEMP_DIR            = "/tmp/"; # must end with "/"
  84
  85
  86
  87 # Alignments in which the number of sequences after pruning (determined
  88 # by "$MY_SPECIES_NAMES_FILE") is lower than this, are ignored
  89 # (no calculation of pwds):
  90 # ------------------------------------------------------------------
  91 my $MIN_SEQS  = 5;
  92
  93
  94
  95 # Alignments in which the number of sequences after pruning (determined
  96 # by "$MY_SPECIES_NAMES_FILE") is greater than this, are ignored
  97 # (no calculation of pwds):
  98 # ------------------------------------------------------------------
  99 my $MAX_SEQS  = 700;
 100
 101
 102
 103 # Seed for the random number generator for bootstrapping (must be 4n+1):
 104 # ---------------------------------------------------------------------
 105 my $MY_SEED   = 85;
 106
 107
 108
 109 # This is used to choose the model to be used for the (ML)
 110 # distance calculation:
 111 # IMPORTANT: "$MY_MATRIX_FOR_PWD" in "rio_module.pm" needs to
 112 # have the same value, when the pwds calculated are going to
 113 # be used for RIO!
 114 # 0 = JTT
 115 # 2 = BLOSUM 62
 116 # 3 = mtREV24
 117 # 5 = VT
 118 # 6 = WAG
 119 # PAM otherwise
 120 # --------------------------------------------------------
 121 my $MY_MATRIX = 2;
 122
 123
 124
 125 #
 126 # End of variables which need to be set by the user.
 127 #
 128 # =============================================================================
 129 # =============================================================================
 130
 131
 132
 133
 134
 135
 136
 137
 138 my $too_small             = 0;
 139 my $too_large             = 0;
 140 my $i                     = 0;
 141 my $seqs                  = 0;
 142 my $filename              = "";
 143 my $tmp_dir               = "";
 144 my $current_dir           = "";
 145 my $return_line           = "";
 146 my @filenames             = ();
 147 my @too_small_names       = ();
 148 my @too_large_names       = ();
 149 my %Species_names_hash    = ();
 150 my %AC_OS                 = (); # AC -> species name
 151 my %AC_DE                 = (); # AC -> description
 152 my %ALGNS_TO_USE          = (); # name of alignment -> ""
 153 my $use_algns_to_use_list = 0;
 154 my $LOGFILE               = "00_pfam2pwd_LOGFILE";
 155    $HMMBUILD              = $HMMBUILD." --amino";
 156
 157
 158 &createTempdir();
 159
 160
 161 &startLogfile();
 162
 163
 164 opendir( DIR, $MY_PFAM_FULL_DIRECTORY ) || die "\n\n$0: Cannot open directory $MY_PFAM_FULL_DIRECTORY: $!\n\n";
 165 $i = 0;
 166 while( defined( $filename = readdir( DIR ) ) ) {
 167     if ( $filename =~ /^\.\.?$/ ) {
 168         next;
 169     }
 170     $filenames[ $i ] = $filename;
 171     $i++;
 172 }
 173 close( DIR );
 174
 175
 176 &readSpeciesNamesFile( $MY_SPECIES_NAMES_FILE );
 177
 178 &readTrEMBL_ACDEOS_FILE();
 179
 180 if ( defined( $ALGNS_TO_USE_LIST_FILE ) && $ALGNS_TO_USE_LIST_FILE =~ /\w/ ) {
 181     $use_algns_to_use_list = 1;
 182     &readListFile();
 183 }
 184
 185
 186 $current_dir = `pwd`;
 187 $current_dir =~ s/\s//;
 188 chdir ( $tmp_dir )
 189 || die "\n\n$0: Unexpected error: Could not chdir to <<$tmp_dir>>: $!";
 190
 191 $i = 0;
 192
 193 FOREACH_ALIGN: foreach $filename ( @filenames ) {
 194
 195     # If the corresponding pwd, positions, and aln files seem to already exists, do next one.
 196     if ( ( -e $MY_RIO_PWD_DIRECTORY.$filename.$SUFFIX_PWD )
 197     &&   ( -e $MY_RIO_BSP_DIRECTORY.$filename.$SUFFIX_BOOT_STRP_POS )
 198     &&   ( -e $MY_RIO_NBD_DIRECTORY.$filename.$SUFFIX_PWD_NOT_BOOTS )
 199     &&   ( -e $MY_RIO_ALN_DIRECTORY.$filename.$ALIGN_FILE_SUFFIX )
 200     &&   ( -e $MY_RIO_HMM_DIRECTORY.$filename.$SUFFIX_HMM ) ) {
 201         next FOREACH_ALIGN;
 202     }
 203
 204     if ( $use_algns_to_use_list == 1 && !exists( $ALGNS_TO_USE{ $filename } ) ) {
 205         next FOREACH_ALIGN;
 206     }
 207
 208
 209     $seqs = &removeSeqsFromPfamAlign( $MY_PFAM_FULL_DIRECTORY.$filename,
 210                                       "REM_SEQ_OUTFILE",
 211                                       1 );
 212     if ( $seqs < $MIN_SEQS ) {
 213         unlink( "REM_SEQ_OUTFILE" );
 214         $too_small_names[ $too_small++ ] = $filename;
 215         next FOREACH_ALIGN;
 216     }
 217     elsif ( $seqs > $MAX_SEQS ) {
 218         unlink( "REM_SEQ_OUTFILE" );
 219         $too_large_names [ $too_large++ ] = $filename;
 220         next FOREACH_ALIGN;
 221     }
 222
 223
 224     print "\n\n\n";
 225     print " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 226     print " $i: $filename ($seqs seqs)\n";
 227     print " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 228     print "\n";
 229
 230     # If one of the two file exists from a previous (interrupted) run.
 231     unlink( $MY_RIO_PWD_DIRECTORY.$filename.$SUFFIX_PWD );
 232     unlink( $MY_RIO_BSP_DIRECTORY.$filename.$SUFFIX_BOOT_STRP_POS );
 233     unlink( $MY_RIO_NBD_DIRECTORY.$filename.$SUFFIX_PWD_NOT_BOOTS );
 234     unlink( $MY_RIO_ALN_DIRECTORY.$filename.$ALIGN_FILE_SUFFIX );
 235     unlink( $MY_RIO_HMM_DIRECTORY.$filename.$SUFFIX_HMM );
 236
 237
 238     &executeHmmbuild( "REM_SEQ_OUTFILE",
 239                       $MY_RIO_ALN_DIRECTORY.$filename.$ALIGN_FILE_SUFFIX,
 240                       "hmm" );
 241
 242     if ( unlink( "hmm" ) != 1 ) {
 243         die "\n\n$0: Unexpected error: Could not delete <<hmm>>: $!";
 244     }
 245
 246     if ( unlink( "REM_SEQ_OUTFILE" ) != 1 ) {
 247         die "\n\n$0: Unexpected error: Could not delete <<REM_SEQ_OUTFILE>>: $!";
 248     }
 249
 250     executeHmmbuildHand( $MY_RIO_ALN_DIRECTORY.$filename.$ALIGN_FILE_SUFFIX,
 251                          $MY_RIO_HMM_DIRECTORY.$filename.$SUFFIX_HMM );
 252
 253     system( $HMMCALIBRATE, $MY_RIO_HMM_DIRECTORY.$filename.$SUFFIX_HMM )
 254     && die "\n\n$0: Could not execute \"$HMMCALIBRATE $MY_RIO_HMM_DIRECTORY.$filename.$SUFFIX_HMM\": $!";
 255
 256     &pfam2phylipMatchOnly( $MY_RIO_ALN_DIRECTORY.$filename.$ALIGN_FILE_SUFFIX, "infile" );
 257
 258     &executePuzzle( "infile", $MY_MATRIX );
 259
 260     system( "mv", "infile.dist", $MY_RIO_NBD_DIRECTORY.$filename.$SUFFIX_PWD_NOT_BOOTS )
 261     && die "\n\n$0: Unexpected error: $!";
 262
 263     &executeBootstrap( "infile",
 264                        $BOOTSTRAPS,
 265                        "BOOTSTRAPPED_ALGN",
 266                        $MY_RIO_BSP_DIRECTORY.$filename.$SUFFIX_BOOT_STRP_POS,
 267                        $MY_SEED );
 268
 269     if ( unlink( "infile" ) != 1 ) {
 270         die "\n\n$0: Unexpected error: Could not delete <<infile>>: $!";
 271     }
 272
 273
 274     &executePuzzleBootstrapped( "BOOTSTRAPPED_ALGN", $MY_MATRIX );
 275
 276     ##if ( unlink( "outfile" ) != 1 ) {
 277     ##    die "\n\n$0: Unexpected error: Could not delete <<outfile>>: $!";
 278     ##}
 279
 280
 281     system( "mv", "BOOTSTRAPPED_ALGN".".dist", $MY_RIO_PWD_DIRECTORY.$filename.$SUFFIX_PWD )
 282     && die "\n\n$0: Unexpected error: $!\n\n";
 283
 284     if ( unlink( "BOOTSTRAPPED_ALGN" ) != 1 ) {
 285         die "\n\n$0: Unexpected error: Could not delete <<BOOTSTRAPPED_ALGN>>: $!";
 286     }
 287
 288     $i++;
 289
 290 } ## End of FOREACH_ALIGN loop.
 291
 292
 293 chdir( $current_dir )
 294 || die "\n\n$0: Unexpected error: Could not chdir to <<$current_dir>>: $!";
 295
 296 rmdir( $tmp_dir );
 297
 298 &finishLogfile();
 299
 300 print "\n\n\n";
 301 print( "pfam2pwd.pl: Done.\n" );
 302 print( "Successfully calculated $i pairwise distance files.\n" );
 303 print( "Too large alignments (>$MAX_SEQS): $too_large\n" );
 304 print( "Too small alignments (<$MIN_SEQS): $too_small\n" );
 305 print( "See the logfile \"$MY_RIO_PWD_DIRECTORY".$LOGFILE."\"\n" );
 306 print "\n\n\n";
 307
 308 exit( 0 );
 309
 310
 311
 312
 313
 314
 315 # Methods
 316 # -------
 317
 318
 319
 320 # Three arguments:
 321 # 1. Stockholm alignment
 322 # 2. Outalignment
 323 # 3. Outhmm
 324 # Returns the options used.
 325 # Last modified: 06/26/01
 326 sub executeHmmbuild {
 327
 328     my $full         = $_[ 0 ];
 329     my $outalignment = $_[ 1 ];
 330     my $outhmm       = $_[ 2 ];
 331     my $options      = "";
 332
 333     unless ( ( -s $full ) && ( -f $full ) && ( -T $full ) ) {
 334         die "\n\n$0: \"$full\" does not exist, is empty, or is not a plain textfile.\n\n";
 335     }
 336
 337     $options = getHmmbuildOptionsFromPfam( $full );
 338
 339     $options =~ s/-f//;
 340     $options =~ s/-g//;
 341     $options =~ s/-s//;
 342     $options =~ s/-F//;
 343     $options =~ s/-A//;
 344     $options =~ s/-o\s+\S+//;
 345     $options =~ s/(\s|^)[^-]\S+/ /g;
 346
 347     if ( $options =~ /--prior/ ) {
 348         my $basename = basename( $full );
 349         $basename .= ".PRIOR";
 350         $options =~ s/--prior/--prior $PRIOR_FILE_DIR$basename/;
 351     }
 352
 353     # Remove for versions of HMMER lower than 2.2.
 354     if ( $options =~ /--informat\s+\S+/ ) {
 355         $options =~ s/--informat\s+\S+/-/;
 356     }
 357
 358     system( "$HMMBUILD $options -o $outalignment $outhmm $full" )
 359     && die "\n\n$0: Could not execute \"$HMMBUILD $options -o $outalignment $outhmm $full\".\n\n";
 360
 361     return $options;
 362
 363 } ## executeHmmbuild.
 364
 365
 366 # Two arguments:
 367 # 1. Stockholm alignment
 368 # 2. Outhmm
 369 # Returns the options used.
 370 # Last modified: 06/26/01
 371 sub executeHmmbuildHand {
 372
 373     my $full         = $_[ 0 ];
 374     my $outhmm       = $_[ 1 ];
 375     my $options      = "";
 376
 377     unless ( ( -s $full ) && ( -f $full ) && ( -T $full ) ) {
 378         die "\n\n$0: \"$full\" does not exist, is empty, or is not a plain textfile.\n\n";
 379     }
 380
 381     $options = getHmmbuildOptionsFromPfam( $full );
 382
 383     $options =~ s/-f//;
 384     $options =~ s/-g//;
 385     $options =~ s/-s//;
 386     $options =~ s/-F//;
 387     $options =~ s/-A//;
 388     $options =~ s/-o\s+\S+//;
 389     $options =~ s/(\s|^)[^-]\S+/ /g;
 390
 391     if ( $options =~ /--prior/ ) {
 392         my $basename = basename( $full );
 393         $basename .= ".PRIOR";
 394         $options =~ s/--prior/--prior $PRIOR_FILE_DIR$basename/;
 395     }
 396
 397     # Remove for versions of HMMER lower than 2.2.
 398     if ( $options =~ /--informat\s+\S+/ ) {
 399         $options =~ s/--informat\s+\S+/-/;
 400     }
 401
 402     system( "$HMMBUILD --hand $options $outhmm $full" )
 403     && die "\n\n$0: Could not execute \"$HMMBUILD -- hand $options $outhmm $full\".\n\n";
 404
 405     return $options;
 406
 407 } ## executeHmmbuildHand.
 408
 409
 410
 411 # One argument:
 412 # Pfam align name.
 413 # Last modified: 02/26/01
 414 sub getHmmbuildOptionsFromPfam {
 415
 416     my $infile      = $_[ 0 ];
 417     my $return_line = "";
 418     my $result      = "";
 419
 420     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 421         die "\n\n$0: \"$infile\" does not exist, is empty, or is not a plain textfile.\n\n";
 422     }
 423
 424     open( GHO, $infile ) || die "\n\n$0: Unexpected error: Cannot open file <<$infile>>: $!";
 425     while ( $return_line = <GHO> ) {
 426         if ( $return_line =~ /^\s*#.*hmmbuild\s+(.+)\s*$/ ) {
 427             $result = $1;
 428             close( GHO );
 429             return $result;
 430         }
 431     }
 432     close( GHO );
 433     return $result;
 434
 435 } ## getHmmbuildOptionsFromPfam
 436
 437
 438
 439 # Similar to the method with the same name in "rio.pl".
 440 # Removes sequences from a Pfam flat file.
 441 # Adds species to TrEMBL seqs.
 442 # It can remove all sequences not from species listed in a species names file.
 443 # It can remove all sequences which do not have a SWISS-PROT name (XXXX_XXXXX)
 444 # Three arguments:
 445 # 1. Pfam flat file name
 446 # 2. outfile name
 447 # 3. 1 to remove TrEMBL seqs with "(FRAGMENT)" in their DE line.
 448 # Returns the number of sequences in the resulting alignment.
 449 # If a query name is given, it returns -1 if query is not found in alignment,
 450 # -10 if the name is not unique.
 451 # Last modified: 05/24/02
 452 sub removeSeqsFromPfamAlign {
 453     my $infile                   = $_[ 0 ];
 454     my $outfile                  = $_[ 1 ];
 455     my $remove_frags             = $_[ 2 ];
 456     my $return_line              = "";
 457     my $saw_sequence_line        = 0;
 458     my $number_of_seqs           = 0;
 459     my $DE                       = "";
 460     my $OS                       = "";
 461     my $AC                       = "";
 462     my $i                        = 0;
 463     my $length                   = 0;
 464     my $seq_name                 = "";
 465     my $seq                      = "";
 466
 467
 468     open( OUT_RNSP, ">$outfile" ) || die "\n\n$0: Unexpected error: Cannot create file \"$outfile\": $!";
 469     open( IN_RNSP, "$infile" ) || die "\n\n$0: Unexpected error: Cannot open file <<$infile>>: $!";
 470     while ( $return_line = <IN_RNSP> ) {
 471
 472         if ( $saw_sequence_line == 1
 473         && !&containsPfamNamedSequence( $return_line )
 474         && !&isPfamCommentLine( $return_line ) ) {
 475             # This is just for counting purposes.
 476             $saw_sequence_line = 2;
 477         }
 478         if ( &isPfamSequenceLine( $return_line ) ) {
 479             if ( $saw_sequence_line == 0 ) {
 480                 $saw_sequence_line = 1;
 481             }
 482             $return_line =~ /^\s*(\S+)\s+(\S+)/;
 483             $seq_name = $1;
 484             $seq      = $2;
 485             if ( !&startsWithSWISS_PROTname( $return_line ) ) {
 486                 $seq_name =~ /^(\S+)\//;
 487                 $AC = $1;
 488                 unless( exists( $AC_OS{ $AC } ) ) {
 489                     #ACs not present in "ACDEOS" file.
 490                     next;
 491                 }
 492                 $OS = $AC_OS{ $AC };
 493                 if ( !$OS || $OS eq "" ) {
 494                     die "\n\n$0: Unexpected error: species for \"$AC\" not found.\n\n";
 495                 }
 496                 unless( exists( $Species_names_hash{ $OS } ) ) {
 497                     next;
 498                 }
 499                 if ( $remove_frags == 1 ) {
 500                     $DE = $AC_DE{ $AC };
 501                     if ( $DE && $DE =~ /\(FRAGMENT\)/ ) {
 502                         next;
 503                     }
 504                 }
 505                 $seq_name =~ s/\//_$OS\//;
 506             }
 507             else {
 508                 if ( $return_line =~ /_([A-Z0-9]{1,5})\// ) {
 509                     unless( exists( $Species_names_hash{ $1 } ) ) {
 510                         next;
 511                     }
 512                 }
 513                 # remove everything whose species cannot be determined.
 514                 else {
 515                     next;
 516                 }
 517             }
 518             $length = length( $seq_name );
 519             for ( $i = 0; $i <= ( $LENGTH_OF_NAME - $length - 1 ); $i++ ) {
 520                     $seq_name .= " ";
 521             }
 522             $return_line = $seq_name.$seq."\n";
 523         }
 524
 525         if ( !&isPfamCommentLine( $return_line ) ) {
 526             print OUT_RNSP $return_line;
 527         }
 528
 529         if ( $saw_sequence_line == 1 ) {
 530             $number_of_seqs++;
 531         }
 532     } ## while ( $return_line = <IN_RNSP> )
 533     close( IN_RNSP );
 534     close( OUT_RNSP );
 535
 536     return $number_of_seqs;
 537
 538 } ## removeSeqsFromPfamAlign
 539
 540
 541
 542
 543
 544
 545
 546 # Reads in (SWISS-PROT) species names from a file.
 547 # Names must be separated by newlines.
 548 # Lines beginning with "#" are ignored.
 549 # A possible "=" and everything after is ignored.
 550 # One argument: species-names-file name
 551 # Last modified: 04/24/01
 552 sub readSpeciesNamesFile {
 553     my $infile      = $_[ 0 ];
 554     my $return_line = "";
 555     my $species     = "";
 556
 557     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 558         die "\n\n$0: Error: \"$infile\" does not exist, is empty, or is not a plain textfile.\n\n";
 559     }
 560
 561     open( IN_RSNF, "$infile" ) || die "\n\n$0: Unexpected error: Cannot open file <<$infile>>: $!\n\n";
 562     while ( $return_line = <IN_RSNF> ) {
 563         if ( $return_line !~ /^\s*#/ && $return_line =~ /(\S+)/ ) {
 564             $species = $1;
 565             $species =~ s/=.+//;
 566             $Species_names_hash{ $species } = "";
 567         }
 568     }
 569     close( IN_RSNF );
 570
 571     return;
 572 } ## readSpeciesNamesFile
 573
 574
 575
 576 # Last modified: 02/21/03
 577 sub readTrEMBL_ACDEOS_FILE {
 578
 579     my $return_line = "";
 580
 581     unless ( ( -s $TREMBL_ACDEOS_FILE  ) && ( -f $TREMBL_ACDEOS_FILE  ) && ( -T $TREMBL_ACDEOS_FILE ) ) {
 582         die "\n\n$0: Error: \"$TREMBL_ACDEOS_FILE\" does not exist, is empty, or is not a plain textfile.\n\n";
 583     }
 584     # Fill up (huge) hashs.
 585     open( HH, "$TREMBL_ACDEOS_FILE" ) || die "\n\n$0: Unexpected error: Cannot open file <<$TREMBL_ACDEOS_FILE>>: $!\n\n";
 586     while ( $return_line = <HH> ) {
 587
 588         if ( $return_line =~ /(\S+);([^;]*);(\S+)/ ) {
 589              $AC_OS{ $1 } = $3;
 590              $AC_DE{ $1 } = $2;
 591         }
 592     }
 593     close( HH );
 594 } ## readTrEMBL_ACDEOS_FILE
 595
 596
 597 # Last modified: 02/21/03
 598 sub readListFile {
 599
 600     my $return_line = "";
 601
 602     unless ( ( -s $ALGNS_TO_USE_LIST_FILE ) && ( -f $ALGNS_TO_USE_LIST_FILE ) && ( -T $ALGNS_TO_USE_LIST_FILE ) ) {
 603         die "\n\n$0: Error: \"$ALGNS_TO_USE_LIST_FILE\" does not exist, is empty, or is not a plain textfile.\n\n";
 604     }
 605     # Fill up hash.
 606     open( LF, "$ALGNS_TO_USE_LIST_FILE" ) || die "\n\n$0: Unexpected error: Cannot open file <<$ALGNS_TO_USE_LIST_FILE>>: $!\n\n";
 607     while ( $return_line = <LF> ) {
 608         if ( $return_line =~ /^\s*(\S+)\s*$/ ) {  # just a list
 609             $ALGNS_TO_USE{ $1 } = "";
 610         }
 611         elsif ( $return_line =~ /^\s*\S+\s+\S+\s+(\S+)/ ) { # "changes" list from Pfam
 612             $ALGNS_TO_USE{ $1 } = "";
 613         }
 614
 615     }
 616     close( LF );
 617
 618 } ## readListFile
 619
 620
 621
 622 # Five arguments:
 623 # 1. Name of inputfile
 624 # 2. Bootstraps
 625 # 2. Name of output alignment file
 626 # 3. Name of output positions file
 627 # 4. Seed for random number generator
 628 #
 629 # Last modified: 06/23/01
 630 sub executeBootstrap {
 631     my $infile     = $_[ 0 ];
 632     my $bootstraps = $_[ 1 ];
 633     my $outalign   = $_[ 2 ];
 634     my $positions  = $_[ 3 ];
 635     my $seed       = $_[ 4 ];
 636
 637     system( "$BOOTSTRAP_CZ_PL 0 $bootstraps $infile $outalign $positions $seed" )
 638     && die "\n\n$0: executeBootstrap:\nCould not execute \"$BOOTSTRAP_CZ_PL 0 $bootstraps $infile $outalign $positions $seed\".\n\n";
 639
 640 } ## executeBootstrap
 641
 642
 643
 644
 645 # Last modified: 05/22/02
 646 sub createTempdir {
 647
 648     my $ii = 0;
 649     my $time = time;
 650
 651     $tmp_dir = $MY_TEMP_DIR.$time.$ii;
 652
 653     while ( -e $tmp_dir ) {
 654         $ii++;
 655         $tmp_dir = $MY_TEMP_DIR.$time.$ii;
 656     }
 657
 658     mkdir(  $tmp_dir, 0777 )
 659     || die "\n\n$0: Unexpected error: Could not create <<$tmp_dir>>: $!\n\n";
 660
 661     unless ( ( -e $tmp_dir ) && ( -d $tmp_dir ) ) {
 662         die "\n\n$0: Unexpected error: failed to create <<$tmp_dir>>.\n\n";
 663     }
 664
 665 } ## createTempdir
 666
 667
 668
 669 # Last modified: 05/17/01
 670 sub startLogfile {
 671     if ( -e $MY_RIO_PWD_DIRECTORY.$LOGFILE ) {
 672         print "\npfam2pwd.pl:\n";
 673         print "logfile $MY_RIO_PWD_DIRECTORY"."$LOGFILE already exists\n";
 674         print "rename it or place it in another directory\n";
 675         exit( -1 );
 676     }
 677
 678     open( L, ">$MY_RIO_PWD_DIRECTORY".$LOGFILE )
 679     || die "\n\n$0: startLogfile: Cannot create logfile: $!\n\n";
 680     print L "Min seqs           : $MIN_SEQS\n";
 681     print L "Max seqs           : $MAX_SEQS\n";
 682     print L "Seed               : $MY_SEED\n";
 683     print L "TrEMBL ACDEOS file : $TREMBL_ACDEOS_FILE\n";
 684     print L "Species names file : $MY_SPECIES_NAMES_FILE\n";
 685     print L "Pfam directory     : $MY_PFAM_FULL_DIRECTORY\n";
 686     print L "PWD outputdirectory: $MY_RIO_PWD_DIRECTORY\n";
 687     print L "BSP outputdirectory: $MY_RIO_BSP_DIRECTORY\n";
 688     print L "NBD outputdirectory: $MY_RIO_NBD_DIRECTORY\n";
 689     print L "ALN outputdirectory: $MY_RIO_ALN_DIRECTORY\n";
 690     print L "HMM outputdirectory: $MY_RIO_HMM_DIRECTORY\n";
 691     print L "Start date         : ".`date`;
 692     if ( $MY_MATRIX == 0 ) {
 693         print L "Matrix             : JTT\n";
 694     }
 695     elsif ( $MY_MATRIX == 2 ) {
 696         print L "Matrix             : BLOSUM 62\n";
 697     }
 698     elsif ( $MY_MATRIX == 3 ) {
 699         print L "Matrix             : mtREV24\n";
 700     }
 701     elsif ( $MY_MATRIX == 5 ) {
 702         print L "Matrix             : VT\n";
 703     }
 704     elsif ( $MY_MATRIX == 6 ) {
 705         print L "Matrix             : WAG\n";
 706     }
 707     elsif ( $MY_MATRIX == 7 ) {
 708         print L "Matrix             : auto\n";
 709     }
 710     else {
 711         print L "Matrix             : PAM\n";
 712     }
 713 } ## startLogfile
 714
 715
 716
 717 # Last modified: 05/17/01
 718 sub finishLogfile {
 719     my $j = 0;
 720     print L "\n\n";
 721     print L "Successfully calculated $i pairwise distance files.\n";
 722     print L "Too large alignments (>$MAX_SEQS): $too_large\n";
 723     print L "Too small alignments (<$MIN_SEQS): $too_small\n";
 724     print L "Finish date        : ".`date`."\n\n";
 725
 726     print L "List of the $too_large alignments which were ignored because they\n";
 727     print L "contained too many sequences (>$MAX_SEQS) after pruning:\n";
 728     for ( $j = 0; $j < $too_large; ++$j ) {
 729         print L "$too_large_names[ $j ]\n";
 730     }
 731     print L "\n\n";
 732     print L "List of the $too_small alignments which were ignored because they\n";
 733     print L "contained not enough sequences (<$MIN_SEQS) after pruning:\n";
 734     for ( $j = 0; $j < $too_small; ++$j ) {
 735         print L "$too_small_names[ $j ]\n";
 736     }
 737     print L "\n";
 738     close( L );
 739 } ## finishLogfile
 740
 741
 742
 743