forester/archive/RIO/others/hmmer/squid/shuffle.c

   1 /*****************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  *****************************************************************/
  10
  11 /* shuffle.c
  12  *
  13  * Routines for randomizing sequences.
  14  *
  15  * All routines are alphabet-independent (DNA, protein, RNA, whatever);
  16  * they assume that input strings are purely alphabetical [a-zA-Z], and
  17  * will return strings in all upper case [A-Z].
  18  *
  19  * All return 1 on success, and 0 on failure; 0 status invariably
  20  * means the input string was not alphabetical.
  21  *
  22  * StrShuffle()   - shuffled string, preserve mono-symbol composition.
  23  * StrDPShuffle() - shuffled string, preserve mono- and di-symbol composition.
  24  *
  25  * StrMarkov0()   - random string, same zeroth order Markov properties.
  26  * StrMarkov1()   - random string, same first order Markov properties.
  27  *
  28  * StrReverse()   - simple reversal of string
  29  * StrRegionalShuffle() -  mono-symbol shuffled string in regional windows
  30  *
  31  * There are also similar routines for shuffling alignments:
  32  *
  33  * AlignmentShuffle()   - alignment version of StrShuffle().
  34  * AlignmentBootstrap() - sample with replacement; a bootstrap dataset.
  35  *
  36  * CVS $Id: shuffle.c,v 1.1.1.1 2005/03/22 08:34:24 cmzmasek Exp $
  37  */
  38
  39 #include <string.h>
  40 #include <ctype.h>
  41
  42 #include "squid.h"
  43
  44 /* Function: StrShuffle()
  45  *
  46  * Purpose:  Returns a shuffled version of s2, in s1.
  47  *           (s1 and s2 can be identical, to shuffle in place.)
  48  *
  49  * Args:     s1 - allocated space for shuffled string.
  50  *           s2 - string to shuffle.
  51  *
  52  * Return:   1 on success.
  53  */
  54 int
  55 StrShuffle(char *s1, char *s2)
  56 {
  57   int  len;
  58   int  pos;
  59   char c;
  60
  61   if (s1 != s2) strcpy(s1, s2);
  62   for (len = strlen(s1); len > 1; len--)
  63     {
  64       pos       = CHOOSE(len);
  65       c         = s1[pos];
  66       s1[pos]   = s1[len-1];
  67       s1[len-1] = c;
  68     }
  69   return 1;
  70 }
  71
  72 /* Function: StrDPShuffle()
  73  * Date:     SRE, Fri Oct 29 09:15:17 1999 [St. Louis]
  74  *
  75  * Purpose:  Returns a shuffled version of s2, in s1.
  76  *           (s1 and s2 may be identical; i.e. a string
  77  *           may be shuffled in place.) The shuffle is a
  78  *           "doublet-preserving" (DP) shuffle. Both
  79  *           mono- and di-symbol composition are preserved.
  80  *
  81  *           Done by searching for a random Eulerian
  82  *           walk on a directed multigraph.
  83  *           Reference: S.F. Altschul and B.W. Erickson, Mol. Biol.
  84  *           Evol. 2:526-538, 1985. Quoted bits in my comments
  85  *           are from Altschul's outline of the algorithm.
  86  *
  87  * Args:     s1   - RETURN: the string after it's been shuffled
  88  *                    (space for s1 allocated by caller)
  89  *           s2   - the string to be shuffled
  90  *
  91  * Returns:  0 if string can't be shuffled (it's not all [a-zA-z]
  92  *             alphabetic.
  93  *           1 on success.
  94  */
  95 int
  96 StrDPShuffle(char *s1, char *s2)
  97 {
  98   int    len;
  99   int    pos;   /* a position in s1 or s2 */
 100   int    x,y;   /* indices of two characters */
 101   char **E;     /* edge lists: E[0] is the edge list from vertex A */
 102   int   *nE;    /* lengths of edge lists */
 103   int   *iE;    /* positions in edge lists */
 104   int    n;     /* tmp: remaining length of an edge list to be shuffled */
 105   char   sf;    /* last character in s2 */
 106   char   Z[26]; /* connectivity in last edge graph Z */
 107   int    keep_connecting; /* flag used in Z connectivity algorithm */
 108   int    is_eulerian;           /* flag used for when we've got a good Z */
 109
 110   /* First, verify that the string is entirely alphabetic.
 111    */
 112   len = strlen(s2);
 113   for (pos = 0; pos < len; pos++)
 114     if (! isalpha(s2[pos])) return 0;
 115
 116   /* "(1) Construct the doublet graph G and edge ordering E
 117    *      corresponding to S."
 118    *
 119    * Note that these also imply the graph G; and note,
 120    * for any list x with nE[x] = 0, vertex x is not part
 121    * of G.
 122    */
 123   E  = MallocOrDie(sizeof(char *) * 26);
 124   nE = MallocOrDie(sizeof(int)    * 26);
 125   for (x = 0; x < 26; x++)
 126     {
 127       E[x]  = MallocOrDie(sizeof(char) * (len-1));
 128       nE[x] = 0;
 129     }
 130
 131   x = toupper(s2[0]) - 'A';
 132   for (pos = 1; pos < len; pos++)
 133     {
 134       y = toupper(s2[pos]) - 'A';
 135       E[x][nE[x]] = y;
 136       nE[x]++;
 137       x = y;
 138     }
 139
 140   /* Now we have to find a random Eulerian edge ordering.
 141    */
 142   sf = toupper(s2[len-1]) - 'A';
 143   is_eulerian = 0;
 144   while (! is_eulerian)
 145     {
 146       /* "(2) For each vertex s in G except s_f, randomly select
 147        *      one edge from the s edge list of E(S) to be the
 148        *      last edge of the s list in a new edge ordering."
 149        *
 150        * select random edges and move them to the end of each
 151        * edge list.
 152        */
 153       for (x = 0; x < 26; x++)
 154         {
 155           if (nE[x] == 0 || x == sf) continue;
 156
 157           pos           = CHOOSE(nE[x]);
 158           y             = E[x][pos];
 159           E[x][pos]     = E[x][nE[x]-1];
 160           E[x][nE[x]-1] = y;
 161         }
 162
 163       /* "(3) From this last set of edges, construct the last-edge
 164        *      graph Z and determine whether or not all of its
 165        *      vertices are connected to s_f."
 166        *
 167        * a probably stupid algorithm for looking at the
 168        * connectivity in Z: iteratively sweep through the
 169        * edges in Z, and build up an array (confusing called Z[x])
 170        * whose elements are 1 if x is connected to sf, else 0.
 171        */
 172       for (x = 0; x < 26; x++) Z[x] = 0;
 173       Z[(int) sf] = keep_connecting = 1;
 174
 175       while (keep_connecting) {
 176         keep_connecting = 0;
 177         for (x = 0; x < 26; x++)
 178           {
 179             y = E[x][nE[x]-1];            /* xy is an edge in Z */
 180             if (Z[x] == 0 && Z[y] == 1)   /* x is connected to sf in Z */
 181               {
 182                 Z[x] = 1;
 183                 keep_connecting = 1;
 184               }
 185           }
 186       }
 187
 188       /* if any vertex in Z is tagged with a 0, it's
 189        * not connected to sf, and we won't have a Eulerian
 190        * walk.
 191        */
 192       is_eulerian = 1;
 193       for (x = 0; x < 26; x++)
 194         {
 195           if (nE[x] == 0 || x == sf) continue;
 196           if (Z[x] == 0) {
 197             is_eulerian = 0;
 198             break;
 199           }
 200         }
 201
 202       /* "(4) If any vertex is not connected in Z to s_f, the
 203        *      new edge ordering will not be Eulerian, so return to
 204        *      (2). If all vertices are connected in Z to s_f,
 205        *      the new edge ordering will be Eulerian, so
 206        *      continue to (5)."
 207        *
 208        * e.g. note infinite loop while is_eulerian is FALSE.
 209        */
 210     }
 211
 212   /* "(5) For each vertex s in G, randomly permute the remaining
 213    *      edges of the s edge list of E(S) to generate the s
 214    *      edge list of the new edge ordering E(S')."
 215    *
 216    * Essentially a StrShuffle() on the remaining nE[x]-1 elements
 217    * of each edge list; unfortunately our edge lists are arrays,
 218    * not strings, so we can't just call out to StrShuffle().
 219    */
 220   for (x = 0; x < 26; x++)
 221     for (n = nE[x] - 1; n > 1; n--)
 222       {
 223         pos       = CHOOSE(n);
 224         y         = E[x][pos];
 225         E[x][pos] = E[x][n-1];
 226         E[x][n-1] = y;
 227       }
 228
 229   /* "(6) Construct sequence S', a random DP permutation of
 230    *      S, from E(S') as follows. Start at the s_1 edge list.
 231    *      At each s_i edge list, add s_i to S', delete the
 232    *      first edge s_i,s_j of the edge list, and move to
 233    *      the s_j edge list. Continue this process until
 234    *      all edge lists are exhausted."
 235    */
 236   iE = MallocOrDie(sizeof(int) * 26);
 237   for (x = 0; x < 26; x++) iE[x] = 0;
 238
 239   pos = 0;
 240   x = toupper(s2[0]) - 'A';
 241   while (1)
 242     {
 243       s1[pos++] = 'A' + x;      /* add s_i to S' */
 244
 245       y = E[x][iE[x]];
 246       iE[x]++;                  /* "delete" s_i,s_j from edge list */
 247
 248       x = y;                    /* move to s_j edge list. */
 249
 250       if (iE[x] == nE[x])
 251         break;                  /* the edge list is exhausted. */
 252     }
 253   s1[pos++] = 'A' + sf;
 254   s1[pos]   = '\0';
 255
 256   /* Reality checks.
 257    */
 258   if (x   != sf)  Die("hey, you didn't end on s_f.");
 259   if (pos != len) Die("hey, pos (%d) != len (%d).", pos, len);
 260
 261   /* Free and return.
 262    */
 263   Free2DArray((void **) E, 26);
 264   free(nE);
 265   free(iE);
 266   return 1;
 267 }
 268
 269
 270 /* Function: StrMarkov0()
 271  * Date:     SRE, Fri Oct 29 11:08:31 1999 [St. Louis]
 272  *
 273  * Purpose:  Returns a random string s1 with the same
 274  *           length and zero-th order Markov properties
 275  *           as s2.
 276  *
 277  *           s1 and s2 may be identical, to randomize s2
 278  *           in place.
 279  *
 280  * Args:     s1 - allocated space for random string
 281  *           s2 - string to base s1's properties on.
 282  *
 283  * Returns:  1 on success; 0 if s2 doesn't look alphabetical.
 284  */
 285 int
 286 StrMarkov0(char *s1, char *s2)
 287 {
 288   int   len;
 289   int   pos;
 290   float p[26];                  /* symbol probabilities */
 291
 292   /* First, verify that the string is entirely alphabetic.
 293    */
 294   len = strlen(s2);
 295   for (pos = 0; pos < len; pos++)
 296     if (! isalpha(s2[pos])) return 0;
 297
 298   /* Collect zeroth order counts and convert to frequencies.
 299    */
 300   FSet(p, 26, 0.);
 301   for (pos = 0; pos < len; pos++)
 302     p[(int)(toupper(s2[pos]) - 'A')] += 1.0;
 303   FNorm(p, 26);
 304
 305   /* Generate a random string using those p's.
 306    */
 307   for (pos = 0; pos < len; pos++)
 308     s1[pos] = FChoose(p, 26) + 'A';
 309   s1[pos] = '\0';
 310
 311   return 1;
 312 }
 313
 314
 315 /* Function: StrMarkov1()
 316  * Date:     SRE, Fri Oct 29 11:22:20 1999 [St. Louis]
 317  *
 318  * Purpose:  Returns a random string s1 with the same
 319  *           length and first order Markov properties
 320  *           as s2.
 321  *
 322  *           s1 and s2 may be identical, to randomize s2
 323  *           in place.
 324  *
 325  * Args:     s1 - allocated space for random string
 326  *           s2 - string to base s1's properties on.
 327  *
 328  * Returns:  1 on success; 0 if s2 doesn't look alphabetical.
 329  */
 330 int
 331 StrMarkov1(char *s1, char *s2)
 332 {
 333   int   len;
 334   int   pos;
 335   int   x,y;
 336   int   i;                      /* initial symbol */
 337   float p[26][26];              /* symbol probabilities */
 338
 339   /* First, verify that the string is entirely alphabetic.
 340    */
 341   len = strlen(s2);
 342   for (pos = 0; pos < len; pos++)
 343     if (! isalpha(s2[pos])) return 0;
 344
 345   /* Collect first order counts and convert to frequencies.
 346    */
 347   for (x = 0; x < 26; x++) FSet(p[x], 26, 0.);
 348
 349   i = x = toupper(s2[0]) - 'A';
 350   for (pos = 1; pos < len; pos++)
 351     {
 352       y = toupper(s2[pos]) - 'A';
 353       p[x][y] += 1.0;
 354       x = y;
 355     }
 356   for (x = 0; x < 26; x++)
 357     FNorm(p[x], 26);
 358
 359   /* Generate a random string using those p's.
 360    */
 361   x = i;
 362   s1[0] = x + 'A';
 363   for (pos = 1; pos < len; pos++)
 364     {
 365       y = FChoose(p[x], 26);
 366       s1[pos] = y + 'A';
 367       x = y;
 368     }
 369   s1[pos] = '\0';
 370
 371   return 1;
 372 }
 373
 374
 375
 376 /* Function: StrReverse()
 377  * Date:     SRE, Thu Nov 20 10:54:52 1997 [St. Louis]
 378  *
 379  * Purpose:  Returns a reversed version of s2, in s1.
 380  *           (s1 and s2 can be identical, to reverse in place)
 381  *
 382  * Args:     s1 - allocated space for reversed string.
 383  *           s2 - string to reverse.
 384  *
 385  * Return:   1.
 386  */
 387 int
 388 StrReverse(char *s1, char *s2)
 389 {
 390   int  len;
 391   int  pos;
 392   char c;
 393
 394   if (s1 != s2) strcpy(s1, s2);
 395   len = strlen(s1);
 396   for (pos = 0; pos < len/2; pos++)
 397     {                           /* swap ends */
 398       c             = s1[len-pos-1];
 399       s1[len-pos-1] = s1[pos];
 400       s1[pos]       = c;
 401     }
 402   return 1;
 403 }
 404
 405 /* Function: StrRegionalShuffle()
 406  * Date:     SRE, Thu Nov 20 11:02:34 1997 [St. Louis]
 407  *
 408  * Purpose:  Returns a regionally shuffled version of s2, in s1.
 409  *           (s1 and s2 can be identical to regionally
 410  *           shuffle in place.) See [Pearson88].
 411  *
 412  * Args:     s1 - allocated space for regionally shuffled string.
 413  *           s2 - string to regionally shuffle
 414  *           w  - window size (typically 10 or 20)
 415  *
 416  * Return:   1.
 417  */
 418 int
 419 StrRegionalShuffle(char *s1, char *s2, int w)
 420 {
 421   int  len;
 422   char c;
 423   int  pos;
 424   int  i, j;
 425
 426   if (s1 != s2) strcpy(s1, s2);
 427   len = strlen(s1);
 428
 429   for (i = 0; i < len; i += w)
 430     for (j = MIN(len-1, i+w-1); j > i; j--)
 431       {
 432         pos     = i + CHOOSE(j-i);
 433         c       = s1[pos];
 434         s1[pos] = s1[j];
 435         s1[j]   = c;
 436       }
 437   return 1;
 438 }
 439
 440
 441 /* Function: AlignmentShuffle()
 442  * Date:     SRE, Sun Apr 22 18:37:15 2001 [St. Louis]
 443  *
 444  * Purpose:  Returns a shuffled version of ali2, in ali1.
 445  *           (ali1 and ali2 can be identical, to shuffle
 446  *           in place.) The alignment columns are shuffled,
 447  *           preserving % identity within the columns.
 448  *
 449  * Args:     ali1 - allocated space for shuffled alignment
 450  *                  [0..nseq-1][0..alen-1]
 451  *           ali2 - alignment to be shuffled
 452  *           nseq - number of sequences in the alignment
 453  *           alen - length of alignment, in columns.
 454  *
 455  * Returns:  int
 456  */
 457 int
 458 AlignmentShuffle(char **ali1, char **ali2, int nseq, int alen)
 459 {
 460   int  i;
 461   int  pos;
 462   char c;
 463
 464   if (ali1 != ali2)
 465     {
 466       for (i = 0; i < nseq; i++) strcpy(ali1[i], ali2[i]);
 467     }
 468
 469   for (i = 0; i < nseq; i++)
 470     ali1[i][alen] = '\0';
 471
 472   for (; alen > 1; alen--)
 473     {
 474       pos = CHOOSE(alen);
 475       for (i = 0; i < nseq; i++)
 476         {
 477           c               = ali1[i][pos];
 478           ali1[i][pos]    = ali1[i][alen-1];
 479           ali1[i][alen-1] = c;
 480         }
 481     }
 482
 483   return 1;
 484 }
 485
 486 /* Function: AlignmentBootstrap()
 487  * Date:     SRE, Sun Apr 22 18:49:14 2001 [St. Louis]
 488  *
 489  * Purpose:  Returns a bootstrapped alignment sample in ali1,
 490  *           constructed from ali2 by sampling columns with
 491  *           replacement.
 492  *
 493  *           Unlike the other shuffling routines, ali1 and
 494  *           ali2 cannot be the same. ali2 is left unchanged.
 495  *           ali1 must be a properly allocated space for an
 496  *           alignment the same size as ali2.
 497  *
 498  * Args:     ali1 - allocated space for bootstrapped alignment
 499  *                  [0..nseq-1][0..alen-1]
 500  *           ali2 - alignment to be bootstrapped
 501  *           nseq - number of sequences in the alignment
 502  *           alen - length of alignment, in columns.
 503  *
 504  * Returns:  1 on success.
 505  */
 506 int
 507 AlignmentBootstrap(char **ali1, char **ali2, int nseq, int alen)
 508 {
 509   int  pos;
 510   int  col;
 511   int  i;
 512
 513   for (pos = 0; pos < alen; pos++)
 514     {
 515       col = CHOOSE(alen);
 516       for (i = 0; i < nseq; i++)
 517         ali1[i][pos] = ali2[i][col];
 518     }
 519   for (i = 0; i < nseq; i++)
 520     ali1[i][alen] = '\0';
 521
 522   return 1;
 523 }
 524
 525
 526
 527
 528 #ifdef TESTDRIVER
 529 /*
 530  * cc -g -o testdriver -DTESTDRIVER -L. shuffle.c -lsquid -lm
 531  */
 532 int
 533 main(int argc, char **argv)
 534 {
 535   char s1[100];
 536   char s2[100];
 537
 538   sre_srandom(42);
 539   strcpy(s2, "GGGGGGGGGGCCCCCCCCCC");
 540   /*  strcpy(s2, "AGACATAAAGTTCCGTACTGCCGGGAT");
 541    */
 542   StrDPShuffle(s1, s2);
 543   printf("DPshuffle: %s\n", s1);
 544   StrMarkov0(s1,s2);
 545   printf("Markov 0 : %s\n", s1);
 546   StrMarkov1(s1,s2);
 547   printf("Markov 1 : %s\n", s1);
 548   return 0;
 549 }
 550 #endif