shuffle.c

这是一个基于HMM 模型的生物多序列比对算法的linux实现版本。hmmer

 * Purpose:  Returns a random string s1 with the same
 *           length and first order Markov properties
 *           as s2. 
 *           
 *           s1 and s2 may be identical, to randomize s2
 *           in place.
 *
 * Args:     s1 - allocated space for random string
 *           s2 - string to base s1's properties on.
 *
 * Returns:  1 on success; 0 if s2 doesn't look alphabetical.
 */
int 
StrMarkov1(char *s1, char *s2)
{
  int   len;
  int   pos; 
  int   x,y;
  int   i;			/* initial symbol */
  float p[26][26];		/* symbol probabilities */

  /* First, verify that the string is entirely alphabetic.
   */
  len = strlen(s2);
  for (pos = 0; pos < len; pos++)
    if (! isalpha((int) s2[pos])) return 0;

  /* Collect first order counts and convert to frequencies.
   */
  for (x = 0; x < 26; x++) FSet(p[x], 26, 0.);

  i = x = toupper((int) s2[0]) - 'A';
  for (pos = 1; pos < len; pos++)
    {
      y = toupper((int) s2[pos]) - 'A';
      p[x][y] += 1.0; 
      x = y;
    }
  for (x = 0; x < 26; x++) 
    FNorm(p[x], 26);

  /* Generate a random string using those p's.
   */
  x = i;
  s1[0] = x + 'A';
  for (pos = 1; pos < len; pos++)
    {
      y = FChoose(p[x], 26);
      s1[pos] = y + 'A';
      x = y;
    } 
  s1[pos] = '\0';

  return 1;
}



/* Function: StrReverse()
 * Date:     SRE, Thu Nov 20 10:54:52 1997 [St. Louis]
 * 
 * Purpose:  Returns a reversed version of s2, in s1.
 *           (s1 and s2 can be identical, to reverse in place)
 * 
 * Args:     s1 - allocated space for reversed string.
 *           s2 - string to reverse.
 *           
 * Return:   1.
 */                
int
StrReverse(char *s1, char *s2)
{
  int  len;
  int  pos;
  char c;
  
  len = strlen(s2);
  for (pos = 0; pos < len/2; pos++)
    {				/* swap ends */
      c             = s2[len-pos-1];
      s1[len-pos-1] = s2[pos];
      s1[pos]       = c;
    }
  if (len%2) { s1[pos] = s2[pos]; } /* copy middle residue in odd-len s2 */
  s1[len] = '\0';
  return 1;
}

/* Function: StrRegionalShuffle()
 * Date:     SRE, Thu Nov 20 11:02:34 1997 [St. Louis]
 * 
 * Purpose:  Returns a regionally shuffled version of s2, in s1.
 *           (s1 and s2 can be identical to regionally 
 *           shuffle in place.) See [Pearson88].
 *           
 * Args:     s1 - allocated space for regionally shuffled string.
 *           s2 - string to regionally shuffle
 *           w  - window size (typically 10 or 20)      
 *           
 * Return:   1.
 */
int
StrRegionalShuffle(char *s1, char *s2, int w)
{
  int  len;
  char c;
  int  pos;
  int  i, j;

  if (s1 != s2) strcpy(s1, s2);
  len = strlen(s1);

  for (i = 0; i < len; i += w)
    for (j = MIN(len-1, i+w-1); j > i; j--)
      {
	pos     = i + CHOOSE(j-i);
	c       = s1[pos];
	s1[pos] = s1[j];
	s1[j]   = c;
      }
  return 1;
}


/* Function: AlignmentShuffle()
 * Date:     SRE, Sun Apr 22 18:37:15 2001 [St. Louis]
 *
 * Purpose:  Returns a shuffled version of ali2, in ali1.
 *           (ali1 and ali2 can be identical, to shuffle
 *           in place.) The alignment columns are shuffled,
 *           preserving % identity within the columns.
 *
 * Args:     ali1 - allocated space for shuffled alignment
 *                  [0..nseq-1][0..alen-1]
 *           ali2 - alignment to be shuffled
 *           nseq - number of sequences in the alignment       
 *           alen - length of alignment, in columns.
 *
 * Returns:  int
 */
int
AlignmentShuffle(char **ali1, char **ali2, int nseq, int alen)
{
  int  i;
  int  pos;
  char c;

  if (ali1 != ali2) 
    {
      for (i = 0; i < nseq; i++) strcpy(ali1[i], ali2[i]);
    }

  for (i = 0; i < nseq; i++)
    ali1[i][alen] = '\0';

  for (; alen > 1; alen--) 
    {
      pos = CHOOSE(alen);
      for (i = 0; i < nseq; i++) 
	{
	  c               = ali1[i][pos];
	  ali1[i][pos]    = ali1[i][alen-1];
	  ali1[i][alen-1] = c;
	}
    }

  return 1;
}

/* Function: AlignmentBootstrap()
 * Date:     SRE, Sun Apr 22 18:49:14 2001 [St. Louis]
 *
 * Purpose:  Returns a bootstrapped alignment sample in ali1, 
 *           constructed from ali2 by sampling columns with 
 *           replacement. 
 *           
 *           Unlike the other shuffling routines, ali1 and 
 *           ali2 cannot be the same. ali2 is left unchanged.
 *           ali1 must be a properly allocated space for an
 *           alignment the same size as ali2.
 *
 * Args:     ali1 - allocated space for bootstrapped alignment
 *                  [0..nseq-1][0..alen-1]
 *           ali2 - alignment to be bootstrapped
 *           nseq - number of sequences in the alignment       
 *           alen - length of alignment, in columns. 
 *                  
 * Returns:  1 on success.                 
 */
int
AlignmentBootstrap(char **ali1, char **ali2, int nseq, int alen)
{
  int  pos;
  int  col;
  int  i;

  for (pos = 0; pos < alen; pos++)
    {
      col = CHOOSE(alen);
      for (i = 0; i < nseq; i++) 
	ali1[i][pos] = ali2[i][col];
    }
  for (i = 0; i < nseq; i++)
    ali1[i][alen] = '\0';

  return 1;
}


/* Function: QRNAShuffle()
 * Date:     SRE, Mon Dec 10 10:14:12 2001 [St. Louis]
 *
 * Purpose:  Shuffle a pairwise alignment x,y while preserving the
 *           position of gaps; return the shuffled alignment in xs,
 *           ys.
 *           
 *           Works by doing three separate
 *           shuffles, of (1) columns with residues in both
 *           x and y, (2) columns with residue in x and gap in y,
 *           and (3) columns with gap in x and residue in y.
 *           
 *           xs,x and ys,y may be identical: that is, to shuffle
 *           an alignment "in place", destroying the original
 *           alignment, just call:
 *              QRNAShuffle(x,y,x,y);
 *
 * Args:     xs, ys: allocated space for shuffled pairwise ali of x,y [L+1]
 *           x, y: pairwise alignment to be shuffled [0..L-1]
 *
 * Returns:  1 on success, 0 on failure.
 *           The shuffled alignment is returned in xs, ys.
 */
int
QRNAShuffle(char *xs, char *ys, char *x, char *y)
{
  int  L;
  int *xycol, *xcol, *ycol;
  int  nxy, nx, ny;
  int  i;
  int  pos, c;
  char xsym, ysym;

  if (xs != x) strcpy(xs, x);
  if (ys != y) strcpy(ys, y);

  /* First, construct three arrays containing lists of the column positions
   * of the three types of columns. (If a column contains gaps in both x and y,
   * we've already simply copied it to the shuffled sequence.)
   */
  L = strlen(x);
  xycol = MallocOrDie(sizeof(int) * L);
  xcol  = MallocOrDie(sizeof(int) * L);
  ycol  = MallocOrDie(sizeof(int) * L);
  nxy = nx = ny = 0;

  for (i = 0; i < L; i++)
    {
      if      (isgap(x[i]) && isgap(y[i]))     { continue; }
      else if (! isgap(x[i]) && ! isgap(y[i])) { xycol[nxy] = i; nxy++; }
      else if (isgap(x[i]))                    { ycol[ny] = i;   ny++;  }
      else if (isgap(y[i]))                    { xcol[nx] = i;   nx++;  }
    }

  /* Second, shuffle the sequences indirectly, via shuffling these arrays.
   * Yow, careful with those indices, and with order of the statements...
   */
  for (; nxy > 1; nxy--) {
    pos          = CHOOSE(nxy);
    xsym             = xs[xycol[pos]];   ysym             = ys[xycol[pos]];    c            = xycol[pos];   
    xs[xycol[pos]]   = xs[xycol[nxy-1]]; ys[xycol[pos]]   = ys[xycol[nxy-1]];  xycol[pos]   = xycol[nxy-1];
    xs[xycol[nxy-1]] = xsym;             ys[xycol[nxy-1]] = ysym;              xycol[pos]   = xycol[nxy-1];
  }
  for (; nx > 1; nx--) {
    pos        = CHOOSE(nx); 
    xsym           = xs[xcol[pos]];  ysym           = ys[xcol[pos]];  c          = xcol[pos];  
    xs[xcol[pos]]  = xs[xcol[nx-1]]; ys[xcol[pos]]  = ys[xcol[nx-1]]; xcol[pos]  = xcol[nx-1]; 
    xs[xcol[nx-1]] = xsym;           ys[xcol[nx-1]] = ysym;           xcol[nx-1] = c;          
  }
  for (; ny > 1; ny--) {
    pos        = CHOOSE(ny); 
    xsym           = xs[ycol[pos]];  ysym           = ys[ycol[pos]];  c          = ycol[pos]; 
    xs[ycol[pos]]  = xs[ycol[ny-1]]; ys[ycol[pos]]  = ys[ycol[ny-1]]; ycol[pos]  = ycol[ny-1];
    xs[ycol[ny-1]] = xsym;           ys[ycol[ny-1]] = ysym;           ycol[ny-1] = c;          
  }

  free(xycol); free(xcol); free(ycol);
  return 1;
}


#ifdef TESTDRIVER
/*
 * cc -g -o testdriver -DTESTDRIVER -L. shuffle.c -lsquid -lm
 */
int 
main(int argc, char **argv)
{
  char s1[100];
  char s2[100];

  sre_srandom(42);
  strcpy(s2, "GGGGGGGGGGCCCCCCCCCC");
  /*  strcpy(s2, "AGACATAAAGTTCCGTACTGCCGGGAT");
   */
  StrDPShuffle(s1, s2);
  printf("DPshuffle: %s\n", s1);
  StrMarkov0(s1,s2);
  printf("Markov 0 : %s\n", s1);
  StrMarkov1(s1,s2);
  printf("Markov 1 : %s\n", s1);

  strcpy(s1, "ACGTACGT--------ACGTACGT----ACGTACGT");
  strcpy(s2, "ACGTACGTACGTACGT------------ACGTACGT");
  QRNAShuffle(s1,s2,s1,s2);
  printf("QRNA : %s\n", s1);
  printf("     : %s\n", s2);

  return 0;
}
#endif