forester/archive/RIO/C/bootstrap_cz.c

   1 /*
   2 # bootstrap_cz
   3 # ------------
   4 # Copyright (C) 1999-2002 Washington University School of Medicine
   5 # and Howard Hughes Medical Institute
   6 # All rights reserved
   7 #
   8 # Author: Christian M. Zmasek
   9 #         zmasek@genetics.wustl.edu
  10 #         http://www.genetics.wustl.edu/eddy/people/zmasek/
  11 #
  12 # Created: 06/06/01
  13 #
  14 # Last modified: 01/27/02
  15 #
  16 # Purpose:
  17 # Bootstrap resamples an alignment in PHYLIP sequential format <bootstraps> times.
  18 # Bootstrapping is not done randomly but according to a BSP (bootstrap
  19 # positions) file.
  20 # The BSP file can be created with the Perl program "bootstrap_cz.pl"
  21 # in mode 0.
  22 # This prgram has the same functionality as "bootstrap_cz.pl" in mode 1.
  23 # Sequence names are normalized to LENGTH_OF_NAME characters.
  24 # The output alignment is in PHYLIP's sequential or interleaved format.
  25 # (These two are the same in this case, since all the seqs will be one
  26 # line in length (no returns in seq).)
  27 #
  28 # Usage: bootstrap_cz <bootstraps> <alignment infile> <BSP file> <outfile>
  29 #        [number of processors]
  30 */
  31
  32
  33
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36 #include <ctype.h>
  37 #include <string.h>
  38
  39
  40 #define LENGTH_OF_NAME 26
  41
  42
  43 static char **names,       /* This stores the sequence names */
  44             **sequences;   /* This stores the sequences */
  45 static int  number_of_seqs,
  46             number_of_colm;
  47
  48
  49 void readInAlignmnet( const char * );
  50 void bootstrapAccordingToBSPfile( int, const char *, const char * );
  51 void checkForMemAllocFailure( void * );
  52 int  fileExists( const char *);
  53 void errorInCommandLine();
  54
  55
  56
  57
  58
  59 /* Reads the seqs and seq-names from inalignment     */
  60 /* into **sequences and **sequences.                 */
  61 /* Inalignment must be in PHYLIP sequential format.  */
  62 /* Last modified: 06/25/01                           */
  63 void readInAlignment( const char *inalignment ) {
  64
  65     FILE *inalignment_fp              = NULL;
  66     char *str                         = NULL;
  67     int  max_length                   = 0;
  68     register char c                   = ' ';
  69     register int  i                   = 0,
  70                   ii                  = 0,
  71                   z                   = 0,
  72                   seq                 = 0;
  73
  74     number_of_seqs = 0;
  75     number_of_colm = 0;
  76
  77     inalignment_fp = fopen( inalignment, "r" );
  78     if ( inalignment_fp == NULL ) {
  79         printf( "\nbootstrap_cz: Error: Could not open alignment file for reading.\n" );
  80         exit( -1 );
  81     }
  82
  83     if ( fscanf( inalignment_fp, "%d", &number_of_seqs ) != 1 ) {
  84         printf( "\nbootstrap_cz: Error: Could not read in number of seqs.\n" );
  85         exit( -1 );
  86     }
  87     if ( fscanf( inalignment_fp, "%d", &number_of_colm ) != 1 ) {
  88         printf( "\nbootstrap_cz: Error: Could not read in number of columns.\n" );
  89         exit( -1 );
  90     }
  91
  92     names = malloc( number_of_seqs * sizeof( char *) );
  93     checkForMemAllocFailure( names );
  94     for ( i = 0; i < number_of_seqs; ++i ) {
  95         names[ i ] = malloc( LENGTH_OF_NAME * sizeof( char ) );
  96         checkForMemAllocFailure( names[ i ] );
  97     }
  98
  99     sequences = malloc( number_of_seqs * sizeof( char * ) );
 100     checkForMemAllocFailure( sequences );
 101     for ( i = 0; i < number_of_seqs; ++i ) {
 102         sequences[ i ] = malloc( number_of_colm * sizeof( char ) );
 103         checkForMemAllocFailure( sequences[ i ] );
 104     }
 105
 106     max_length = ( 30 * LENGTH_OF_NAME ) + number_of_colm;
 107
 108     str = malloc( max_length * sizeof( char * ) );
 109     checkForMemAllocFailure( str );
 110
 111     while ( fgets( str, max_length, inalignment_fp ) != NULL ) {
 112
 113         if ( !isspace( str[ 0 ] ) != 0 ) {
 114
 115             i = 0;
 116             while ( str[ i ] != ' ' ) {
 117                 names[ seq ][ i ] = str[ i ];
 118                 i++;
 119             }
 120
 121             ii = i;
 122             while ( ii < LENGTH_OF_NAME ) {
 123                 names[ seq ][ ii ] = ' ';
 124                 ii++;
 125             }
 126
 127             z = 0;
 128
 129             while ( str[ i ] != '\n' && str[ i ] != '\r' && str[ i ] != '\0' ) {
 130                 c = str[ i ];
 131                 if ( c != ' ' ) {
 132                     if ( isupper( c ) != 0 || c == '-' ) {
 133                         sequences[ seq ][ z++ ] = c;
 134                     }
 135                     else {
 136                         printf( "\nbootstrap_cz: Error: Sequence must be represented by uppercase letters A-Z and \"-\" only.\n" );
 137                         exit( -1 );
 138                     }
 139                 }
 140                 i++;
 141                 if ( z > number_of_colm ) {
 142                     printf( "\nbootstrap_cz: Error: line in \"%s\" contains more than %d columns.\n",
 143                             inalignment, number_of_colm );
 144                     exit( -1 );
 145                 }
 146             }
 147             if ( z != number_of_colm ) {
 148                 printf( "\nbootstrap_cz: Error: line in \"%s\" contains a incorrect number of columns.\n",
 149                         inalignment );
 150                 exit( -1 );
 151             }
 152
 153             seq++;
 154
 155             if ( seq > number_of_seqs ) {
 156                 printf( "\nbootstrap_cz: Error: \"%s\" contains more than %d seqs.\n",
 157                          inalignment, number_of_seqs );
 158                 exit( -1 );
 159             }
 160         }
 161
 162
 163     } /* while ( fgets ) */
 164
 165     if ( seq != number_of_seqs ) {
 166         printf( "\nbootstrap_cz: Error: \"%s\" contains a incorrect number of seqs.\n",
 167                 inalignment );
 168         exit( -1 );
 169     }
 170
 171     fclose( inalignment_fp );
 172
 173     return;
 174
 175 } /* readInAlignment */
 176
 177
 178
 179 /* Rearrenges the aa in sequences according to    */
 180 /* the bsp (bootstrap positions) file bsp_file.   */
 181 /* Writes the results to outfile                  */
 182 /* Last modified: 06/07/01                        */
 183 void bootstrapAccordingToBSPfile( int bootstraps,
 184                                   const char *bsp_file,
 185                                   const char *outfile ) {
 186
 187     FILE          *bsp_file_fp = NULL,
 188                   *outfile_fp  = NULL;
 189     int           *positions   = NULL,
 190                   p            = 0;
 191     register int  boot         = 0,
 192                   seq          = 0,
 193                   i            = 0;
 194
 195     positions = malloc( number_of_colm * sizeof( int ) );
 196     checkForMemAllocFailure( positions );
 197
 198
 199     bsp_file_fp = fopen( bsp_file, "r" );
 200     if ( bsp_file_fp == NULL ) {
 201         printf( "\nbootstrap_cz: Error: could not open file \"%s\" for reading.\n",
 202                 bsp_file  );
 203         exit( -1 );
 204     }
 205
 206     outfile_fp = fopen( outfile, "w" );
 207     if ( outfile_fp == NULL ) {
 208         printf( "\nbootstrap_cz: Error: could not open file \"%s\" for writing.\n",
 209                 outfile );
 210         exit( -1 );
 211     }
 212
 213     for ( boot = 0; boot < bootstraps; ++boot ) {
 214
 215         for ( i = 0; i < number_of_colm; ++i ) {
 216             if ( fscanf( bsp_file_fp, "%d", &p ) != 1 ) {
 217                 printf( "\nbootstrap_cz: Error: file \"%s\" does not correspond to alignment.\n",
 218                         bsp_file );
 219                 exit( -1 );
 220             }
 221             positions[ i ] = p;
 222         }
 223
 224         fprintf( outfile_fp, " %d  %d\n", number_of_seqs, number_of_colm );
 225         for ( seq = 0; seq < number_of_seqs; ++seq ) {
 226             for ( i = 0; i < LENGTH_OF_NAME; ++i ) {
 227                 fprintf( outfile_fp, "%c", names[ seq ][ i ] );
 228             }
 229             for ( i = 0; i < number_of_colm; ++i ) {
 230                 fprintf( outfile_fp, "%c", sequences[ seq ][ positions[ i ] ] );
 231             }
 232             fprintf( outfile_fp, "\n" );
 233         }
 234     }
 235
 236     /* Now, the bsp file must not contain any more numbers */
 237     if ( fscanf( bsp_file_fp, "%d", &p ) == 1 ) {
 238         printf( "\nbootstrap_cz: Error: file \"%s\" does not correspond to alignment (too long).\n",
 239                 bsp_file );
 240         printf( ">%d<\n", p );
 241         printf( "number of seqs=%d\n", number_of_seqs );
 242         exit( -1 );
 243     }
 244
 245     fclose( bsp_file_fp );
 246     fclose( outfile_fp );
 247     free( positions );
 248     return;
 249
 250 } /* bootstrapAccordingToBSPfile */
 251
 252
 253
 254 /* Rearrenges the aa in sequences according to    */
 255 /* the bsp (bootstrap positions) file bsp_file.   */
 256 /* Writes the results to outfile                  */
 257 /* Last modified: 01/25/02                        */
 258 void bootstrapAccordingToBSPfileP( int bootstraps,
 259                                    int processors,
 260                                    const char *bsp_file,
 261                                    const char *outfile ) {
 262
 263     FILE          *bsp_file_fp  = NULL,
 264                   *outfile_fp   = NULL;
 265     int           *positions    = NULL,
 266                   p             = 0;
 267     char          *outfile_     = NULL;
 268     register int  boot          = 0,
 269                   seq           = 0,
 270                   i             = 0,
 271                   j             = 0,
 272                   z             = 0,
 273                   flag          = 0;
 274     int           block_size    = 0,
 275                   larger_blocks = 0;
 276
 277     block_size    = ( int ) bootstraps / processors;
 278     larger_blocks = bootstraps - ( block_size * processors ); /* number of blocks which have a size of
 279                                                                  block_size + 1 */
 280
 281     positions = malloc( number_of_colm * sizeof( int ) );
 282     checkForMemAllocFailure( positions );
 283
 284     outfile_ = malloc( ( strlen( outfile ) + 20 ) * sizeof( char ) );
 285     checkForMemAllocFailure( outfile_ );
 286
 287     bsp_file_fp = fopen( bsp_file, "r" );
 288     if ( bsp_file_fp == NULL ) {
 289         printf( "\nbootstrap_cz: Error: could not open file \"%s\" for reading.\n",
 290                 bsp_file  );
 291         exit( -1 );
 292     }
 293
 294     j    = -1;
 295     flag = 1;
 296     z    = 0;
 297
 298     for ( boot = 0; boot < bootstraps; ++boot ) {
 299
 300         for ( i = 0; i < number_of_colm; ++i ) {
 301             if ( fscanf( bsp_file_fp, "%d", &p ) != 1 ) {
 302                 printf( "\nbootstrap_cz: Error: file \"%s\" does not correspond to alignment.\n",
 303                         bsp_file );
 304                 exit( -1 );
 305             }
 306             positions[ i ] = p;
 307         }
 308
 309         j++;
 310
 311         if ( larger_blocks > 0 ) {
 312             if ( j >= block_size + 1 ) {
 313                 flag = 1;
 314                 j    = 0;
 315                 larger_blocks--;
 316             }
 317         }
 318         else if ( j >= block_size ) {
 319             flag = 1;
 320             j    = 0;
 321         }
 322
 323         if ( flag == 1 ) {
 324             if ( boot > 0 ) {
 325                 fclose( outfile_fp );
 326             }
 327             sprintf( outfile_, "%s%d", outfile, z++ );
 328             if ( fileExists( outfile_ ) == 1 ) {
 329                 printf( "\nbootstrap_cz: Error: outfile \"%s\" already exists.\n",
 330                         outfile_ );
 331                 exit( -1 );
 332             }
 333             outfile_fp = fopen( outfile_, "w" );
 334             if ( outfile_fp == NULL ) {
 335                 printf( "\nbootstrap_cz: Error: could not open file \"%s\" for writing.\n",
 336                         outfile_ );
 337                 exit( -1 );
 338             }
 339             flag = 0;
 340         }
 341
 342         fprintf( outfile_fp, " %d  %d\n", number_of_seqs, number_of_colm );
 343         for ( seq = 0; seq < number_of_seqs; ++seq ) {
 344             for ( i = 0; i < LENGTH_OF_NAME; ++i ) {
 345                 fprintf( outfile_fp, "%c", names[ seq ][ i ] );
 346             }
 347             for ( i = 0; i < number_of_colm; ++i ) {
 348                 fprintf( outfile_fp, "%c", sequences[ seq ][ positions[ i ] ] );
 349             }
 350             fprintf( outfile_fp, "\n" );
 351         }
 352     }
 353
 354     /* Now, the bsp file must not contain any more numbers */
 355     if ( fscanf( bsp_file_fp, "%d", &p ) == 1 ) {
 356         printf( "\nbootstrap_cz: Error: file \"%s\" does not correspond to alignment (too long).\n",
 357                 bsp_file );
 358         printf( ">%d<\n", p );
 359         printf( "number of seqs=%d\n", number_of_seqs );
 360         exit( -1 );
 361     }
 362
 363     fclose( bsp_file_fp );
 364     fclose( outfile_fp );
 365
 366     free( positions );
 367     free( outfile_ );
 368
 369     return;
 370
 371 } /* bootstrapAccordingToBSPfileP */
 372
 373
 374
 375
 376 /* Exits if *p is NULL.        */
 377 /* Last modified: 06/06/01     */
 378 void checkForMemAllocFailure( void *p ) {
 379     if ( p == NULL ) {
 380         printf( "\nbootstrap_cz: Memory allocation failed.\n" );
 381         exit( -1 );
 382     }
 383     else {
 384         return;
 385     }
 386 } /* checkForMemAllocFailure */
 387
 388
 389
 390 /* Returns 1 if filename can be opened. */
 391 /* Returns 0 otherwise.                 */
 392 /* Last modified: 06/07/01              */
 393 int fileExists( const char *filename ) {
 394     FILE *fp = NULL;
 395     if ( ( fp = fopen( filename, "r" ) ) != NULL ) {
 396         fclose( fp );
 397         return 1;
 398     }
 399     else {
 400         return 0;
 401     }
 402 } /* fileExists */
 403
 404
 405
 406 void errorInCommandLine() {
 407     printf( "\n" );
 408     printf( " bootstrap_cz  version 3.000\n" );
 409     printf( " ---------------------------\n\n" );
 410     printf( " Purpose:\n" );
 411     printf( " Bootstrap resamples an alignment in PHYLIP sequential format <bootstraps> times.\n" );
 412     printf( " Bootstrapping is not done randomly but according to a BSP (bootstrap\n" );
 413     printf( " positions) file.\n" );
 414     printf( " The BSP file can be created with the Perl program \"bootstrap_cz.pl\"\n" );
 415     printf( " in mode 0.\n" );
 416     printf( " This prgram has the same functionality as \"bootstrap_cz.pl\" in mode 1.\n" );
 417     printf( " Sequence names are normalized to LENGTH_OF_NAME characters.\n" );
 418     printf( " The output alignment is in PHYLIP's sequential or interleaved format.\n" );
 419     printf( " (These two are the same in this case, since all the seqs will be one\n" );
 420     printf( " line in length (no returns in seq).)\n\n" );
 421     printf( " Usage: bootstrap_cz <bootstraps> <alignment infile> <BSP file> <outfile>\n" );
 422     printf( "        [number of processors]\n\n" );
 423 } /* errorInCommandLine */
 424
 425
 426
 427 int main( int argc, char *argv[] ) {
 428
 429     char *inalign   = NULL,
 430          *bsp_file  = NULL,
 431          *outfile   = NULL;
 432     int  bootstraps = 0,
 433          processors = 0;
 434
 435
 436     if ( argc != 5 && argc != 6 ) {
 437         errorInCommandLine();
 438         exit( -1 );
 439     }
 440
 441     bootstraps = atoi( argv[ 1 ] );
 442     inalign    = argv[ 2 ];
 443     bsp_file   = argv[ 3 ];
 444     outfile    = argv[ 4 ];
 445
 446     if ( bootstraps < 1 ) {
 447         errorInCommandLine();
 448         exit( -1 );
 449     }
 450
 451     if ( argc == 6 ) {
 452         processors = atoi( argv[ 5 ] );
 453         if ( processors < 1 ) {
 454             errorInCommandLine();
 455             exit( -1 );
 456         }
 457         if ( processors > bootstraps ) {
 458             processors = bootstraps;
 459         }
 460     }
 461
 462     if ( argc == 5 && fileExists( outfile ) == 1 ) {
 463         printf( "\nbootstrap_cz: Error: outfile \"%s\" already exists.\n",
 464                 outfile );
 465         exit( -1 );
 466     }
 467
 468     readInAlignment( inalign );
 469
 470     if ( argc == 5 ) {
 471         bootstrapAccordingToBSPfile( bootstraps,
 472                                      bsp_file,
 473                                      outfile );
 474     }
 475     else {
 476         bootstrapAccordingToBSPfileP( bootstraps,
 477                                       processors,
 478                                       bsp_file,
 479                                       outfile );
 480     }
 481
 482     return 0;
 483
 484 } /* main */