forester/archive/perl/bootstrap_cz.pl

   1 #!/usr/bin/perl -w
   2
   3 # bootstrap_cz.pl
   4 # ---------------
   5 # Copyright (C) 1999-2003 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 08/26/03
  16 #
  17 # Purpose:
  18 # Bootstrap resamples an alignment in PHYLIP sequential format <bootstraps>
  19 # times.
  20 # Amino acid sequences must only be represented by uppercase letters (A-Z)
  21 # and '-'.
  22 # In mode 0 it saves the positions which it used to create the
  23 # bootstrapped alignment into <positions outfile>.
  24 # Mode 1 allows to recreate exactly the same boostrapped alignment
  25 # by reading in a <positions infile>.
  26 # Sequence names are normalized to $LENGTH_OF_NAME characters.
  27 # The output alignment is in PHYLIP's sequential or interleaved format.
  28 # (These two are the same in this case, since all the seqs will be one
  29 # line in length (no returns in seq).)
  30 #
  31 # Usage:
  32 # bootstrap_cz.pl <mode (0 or 1)> <bootstraps> <alignment infile>
  33 # <alignment outfile> <positions out- (mode 0) or in-file (mode 1)>
  34 # [random number seed (mode 0 only)]
  35 #
  36
  37 use strict;
  38 use FindBin;
  39 use lib $FindBin::Bin;
  40
  41 use rio_module;
  42
  43 my $VERSION        = "2.001";
  44
  45 my $modus          = -1;  # 0 to create pos. file, 1 to use premade pos. file
  46 my $bootstraps     = -1;
  47 my $infile         = "";
  48 my $outalign_file  = "";
  49 my $positions_file = "";
  50 my $seed           = -1;
  51
  52
  53 $modus          = $ARGV[ 0 ];
  54 $bootstraps     = $ARGV[ 1 ];
  55 $infile         = $ARGV[ 2 ];
  56 $outalign_file  = $ARGV[ 3 ];
  57 $positions_file = $ARGV[ 4 ];
  58 $seed           = $ARGV[ 5 ];
  59
  60 if ( @ARGV != 5 && @ARGV != 6 ) {
  61     &printUsage();
  62     exit( -1 );
  63 }
  64
  65 if ( $modus != 0 && $modus != 1 ) {
  66     &printUsage();
  67     exit( -1 );
  68 }
  69
  70 if ( $modus == 0 && @ARGV != 6 ) {
  71     &printUsage();
  72     exit( -1 );
  73 }
  74
  75 if ( $modus == 1 && @ARGV != 5 ) {
  76     &printUsage();
  77     exit( -1 );
  78 }
  79
  80 if ( $bootstraps < 1 ) {
  81     &printUsage();
  82     exit( -1 );
  83 }
  84
  85 if ( $seed && $seed < 0 ) {
  86     &printUsage();
  87     exit( -1 );
  88 }
  89
  90
  91 unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
  92     die "\n\nbootstrap_cz.pl: \"$infile\" does not exist, is empty, or is not a plain textfile.\n\n";
  93 }
  94 if ( -e $outalign_file ) {
  95     die "\n\nbootstrap_cz.pl: \"$outalign_file\" already exists.\n\n";
  96 }
  97
  98 if ( $modus == 0 ) {
  99     if ( -e $positions_file ) {
 100         die "\n\nbootstrap_cz.pl: \"$positions_file\" already exists.\n\n";
 101     }
 102 }
 103 else {
 104     unless ( ( -s $positions_file ) && ( -f $positions_file ) && ( -T $positions_file ) ) {
 105         die "\n\nbootstrap_cz.pl: \"$positions_file\" does not exist, is empty, or is not a plain textfile.\n\n";
 106     }
 107 }
 108
 109 if ( $modus == 0 ) {
 110     &bootstrap( $modus, $bootstraps, $infile, $outalign_file, $positions_file, $seed );
 111 }
 112 else {
 113     &bootstrap( $modus, $bootstraps, $infile, $outalign_file, $positions_file );
 114 }
 115
 116
 117 exit( 0 );
 118
 119
 120
 121 # Methods
 122 # -------
 123
 124
 125 # Five/six arguemnts:
 126 # 1. Mode: 0 to create pos. file, 1 to use premade pos. file
 127 # 2. bootstraps
 128 # 3. Alignment infile name
 129 # 4. Outfile name
 130 # 5. file name for positions file (created if mode is 0, read if mode is 1)
 131 # [6. If modus is 0: seed for random number generator]
 132 #
 133 # This method is very similar to method "pfam2phylip" "in makeTree.pl".
 134 #
 135 # Last modified: 05/17/01
 136 #
 137 sub bootstrap {
 138
 139     my $modus           = $_[ 0 ];
 140     my $bootstraps      = $_[ 1 ];
 141     my $infile          = $_[ 2 ];
 142     my $outalign_file   = $_[ 3 ];
 143     my $positions_file  = $_[ 4 ];
 144
 145     my @seq_name        = ();
 146     my @seq_array       = ();
 147     my @random_numbers  = ();
 148     my $return_line     = "";
 149     my $seq             = "";
 150     my $x               = 0;
 151     my $y               = 0;
 152     my $seq_no          = 0;
 153     my $original_length = 0;
 154     my $max_x           = 0;
 155     my $n               = 0;
 156     my $i               = 0;
 157     my $random          = -1;
 158     my $length          = 0;
 159     my $number_of_seqs  = 0;
 160     my $number_of_colm  = 0;
 161
 162
 163     # Checks the arguments
 164     # --------------------
 165
 166     if ( $modus == 0 ) {
 167         if ( !$_[ 5 ] ) {
 168             die "\n\n$0: bootstrap: Failed to give a seed for random number generator.\n\n";
 169         }
 170         srand( $_[ 5 ] );
 171     }
 172     elsif( $modus == 1 ) {
 173         if ( $_[ 5 ] ) {
 174             die "\n\n$0: bootstrap: Must not give a seed for random number generator.\n\n";
 175         }
 176         unless ( ( -s $positions_file ) && ( -f $positions_file ) && ( -T $positions_file ) ) {
 177             die "\n\n$0: bootstrap: <<$positions_file>> does not exist, is empty, or is not a plain textfile.\n\n";
 178         }
 179     }
 180     else {
 181         die "\n\n$0: bootstrap: modus must be either 0 or 1.\n\n";
 182     }
 183
 184     unless ( ( -s $infile ) && ( -f $infile ) && ( -T $infile ) ) {
 185         die "\n\n$0: bootstrap: <<$infile>> does not exist, is empty, or is not a plain textfile.\n\n";
 186     }
 187
 188
 189
 190     # Reads in the alignment
 191     # ----------------------
 192
 193     open( IN, "$infile" ) || die "\n$0: bootstrap: Cannot open file <<$infile>>: $!";
 194     while ( $return_line = <IN> ) {
 195
 196         if ( $return_line =~ /^\s*(\d+)\s+(\d+)/ ) {
 197             $number_of_seqs = $1;
 198             $number_of_colm = $2;
 199         }
 200         elsif ( $return_line =~ /^(\S+)\s+(\S+)/ ) {
 201             $seq_name[ $seq_no ] = substr( $1, 0, $LENGTH_OF_NAME );
 202             $seq = $2;
 203             if ( $original_length == 0 ) {
 204                 $original_length = length( $seq );
 205             }
 206             elsif ( $original_length != length( $seq ) ) {
 207                 die "\n\n$0: Sequences do not have the same length.\n\n";
 208             }
 209             for ( $x = 0; $x < $original_length; $x++ ) {
 210                 $seq_array[ $x ][ $seq_no ] = substr( $seq, $x, 1 );
 211             }
 212             $seq_no++;
 213         }
 214     }
 215     close( IN );
 216
 217     if ( ( $number_of_seqs != $seq_no )
 218     || ( $number_of_colm != $original_length ) ) {
 219         die "\n\n$0: Number of sequences or number of columns are inconsisten with the values given in the alignment.\n\n";
 220     }
 221
 222     # Adusts the length of the names to $LENGTH_OF_NAME
 223     # -------------------------------------------------
 224
 225     for ( $y = 0; $y < $seq_no; $y++ ) {
 226         $length = length( $seq_name[ $y ] );
 227         for ( $i = 0; $i <= ( $LENGTH_OF_NAME - $length - 1 ); $i++ ) {
 228                 $seq_name[ $y ] .= " ";
 229         }
 230     }
 231
 232
 233
 234     # Bootstraps $bootstraps times and writes the outputfiles
 235     # -------------------------------------------------------
 236
 237     open( OUT, ">$outalign_file" ) || die "\n\n$0: bootstrap: Cannot create file <<$outalign_file>>: $!";
 238     if ( $modus == 0 ) {
 239         open( OUT_P, ">$positions_file" ) || die "\n\n$0: bootstrap: Cannot create file <<$positions_file>>: $!";
 240     }
 241     else {
 242         open( IN_P, "$positions_file" ) || die "\n\n$0: bootstrap: Cannot open file <<$positions_file>>: $!";
 243     }
 244
 245     for ( $n = 0; $n < $bootstraps; $n++ ) {
 246
 247         if ( $modus == 0 ) {
 248             for ( $x = 0; $x < $original_length; $x++ ) {
 249                 $random = int( rand( $original_length ) );
 250                 print OUT_P "$random ";
 251                 $random_numbers[ $x ] = $random;
 252             }
 253             print OUT_P "\n";
 254         }
 255         else {
 256             $return_line = <IN_P>;
 257             if ( !$return_line || $return_line !~ /\d/ ) {
 258                 die "\n\n$0: bootstrap: <<$positions_file>> seems too short or otherwise unsuitable.\n\n";
 259             }
 260             $return_line =~ s/^\s+//;
 261             $return_line =~ s/\s+$//;
 262             @random_numbers = split( /\s+/, $return_line );
 263             if ( scalar( @random_numbers ) != $original_length ) {
 264                die "\n\n$0: bootstrap: <<$positions_file>> seems not to correspond to <<$infile>>.\n\n";
 265             }
 266         }
 267
 268         print OUT " $seq_no  $original_length\n";
 269
 270         for ( $y = 0; $y < $seq_no; $y++ ) {
 271             print OUT "$seq_name[ $y ]";
 272
 273             for ( $x = 0; $x < $original_length; $x++ ) {
 274                 $random = $random_numbers[ $x ];
 275                 if ( !$seq_array[ $random ][ $y ] || $seq_array[ $random ][ $y ] !~ /[A-Z]|-/ ) {
 276                     die "\n\n$0: Sequence must be represented by uppercase letters A-Z and \"-\" only.\n\n";
 277                 }
 278                 print OUT $seq_array[ $random ][ $y ];
 279             }
 280             print OUT "\n";
 281         }
 282     }
 283
 284     close( OUT );
 285
 286     if ( $modus == 0 ) {
 287         print OUT_P "\n";
 288         close( OUT_P );
 289     }
 290     else {
 291         close( IN_P );
 292     }
 293
 294     return;
 295
 296 } ## bootstrap
 297
 298
 299
 300 sub printUsage {
 301     print "\n";
 302     print " bootstrap_cz.pl  $VERSION\n";
 303     print " ---------------\n";
 304     print "\n";
 305     print " Christian Zmasek (zmasek\@genetics.wustl.edu)\n";
 306     print "\n";
 307     print " Purpose:\n";
 308     print " Bootstrap resamples an alignment in PHYLIP sequential format\n";
 309     print " <bootstraps> times.\n";
 310     print " In mode 0 it saves the positions which it used to create the\n";
 311     print " bootstrapped alignment into <positions outfile>.\n";
 312     print " Mode 1 allows to recreate exactly the same boostrapped alignment\n";
 313     print " by reading in a <positions infile>.\n";
 314     print " Sequence names are normalized to $LENGTH_OF_NAME characters.\n";
 315     print " The output alignment is in PHYLIP's sequential or interleaved format.\n";
 316     print " (These two are the same in this case, since all the seqs will be one\n";
 317     print " line in length (no returns in seq).)\n";
 318     print "\n";
 319     print " Usage:\n";
 320     print " bootstrap_cz.pl <mode (0 or 1)> <bootstraps> <alignment infile>\n";
 321     print " <alignment outfile> <positions out (mode 0) or infile (mode 1)>\n";
 322     print " [random number seed (mode 0 only)]\n";
 323     print "\n";
 324 } ## printUsage
 325