forester/archive/RIO/others/hmmer/squid/alignio.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 /* alignio.c
  12  * SRE, Mon Jul 12 11:57:37 1993
  13  * RCS $Id: alignio.c,v 1.1.1.1 2005/03/22 08:34:27 cmzmasek Exp $
  14  *
  15  * Input/output of sequence alignments.
  16  */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include <ctype.h>
  22 #include "squid.h"
  23
  24 #ifdef MEMDEBUG
  25 #include "dbmalloc.h"
  26 #endif
  27
  28 /* Function: AllocAlignment()
  29  *
  30  * Purpose:  Allocate space for an alignment, given the number
  31  *           of sequences and the alignment length in columns.
  32  *
  33  * Args:     nseq     - number of sequences
  34  *           alen     - width of alignment
  35  *           ret_aseq - RETURN: alignment itself
  36  *           ainfo    - RETURN: other info associated with alignment
  37  *
  38  * Return:   (void)
  39  *           aseq, ainfo free'd by caller: FreeAlignment(aseq, &ainfo).
  40  *           note that ainfo itself is alloc'ed in caller, usually
  41  *           just by a "AINFO ainfo" definition.
  42  */
  43 void
  44 AllocAlignment(int nseq, int alen, char ***ret_aseq, AINFO *ainfo)
  45 {
  46   char **aseq;
  47   int idx;
  48
  49   InitAinfo(ainfo);
  50
  51   aseq = (char **) MallocOrDie (sizeof(char *) * nseq);
  52   for (idx = 0; idx < nseq; idx++)
  53     aseq[idx] = (char *) MallocOrDie (sizeof(char) * (alen+1));
  54
  55   ainfo->alen  = alen;
  56   ainfo->nseq  = nseq;
  57
  58   ainfo->wgt   = (float *) MallocOrDie (sizeof(float) * nseq);
  59   FSet(ainfo->wgt, nseq, 1.0);
  60
  61   ainfo->sqinfo = (SQINFO *) MallocOrDie (sizeof(SQINFO) * nseq);
  62   for (idx = 0; idx < nseq; idx++)
  63     ainfo->sqinfo[idx].flags = 0;
  64
  65   *ret_aseq = aseq;
  66 }
  67
  68
  69 /* Function: InitAinfo()
  70  * Date:     SRE, Tue Jan 19 10:16:02 1999 [St. Louis]
  71  *
  72  * Purpose:  Initialize the fields in ainfo structure to
  73  *           default (null) values. Does nothing with
  74  *           fields that are dependent on nseq or alen.
  75  *
  76  * Args:     ainfo  - optional info structure for an alignment
  77  *
  78  * Returns:  (void). ainfo is modified.
  79  */
  80 void
  81 InitAinfo(AINFO *ainfo)
  82 {
  83   ainfo->name  = NULL;
  84   ainfo->desc  = NULL;
  85   ainfo->cs    = NULL;
  86   ainfo->rf    = NULL;
  87   ainfo->acc   = NULL;
  88   ainfo->au    = NULL;
  89   ainfo->flags = 0;
  90
  91   ainfo->tc1  = ainfo->tc2 = 0.0;
  92   ainfo->nc1  = ainfo->nc2 = 0.0;
  93   ainfo->ga1  = ainfo->ga2 = 0.0;
  94 }
  95
  96
  97 /* Function: FreeAlignment()
  98  *
  99  * Purpose:  Free the space allocated to alignment, names, and optional
 100  *           information.
 101  *
 102  * Args:     aseqs - sequence alignment
 103  *           ainfo - associated alignment data.
 104  */
 105 void
 106 FreeAlignment(char **aseqs, AINFO *ainfo)
 107 {
 108   int i;
 109
 110   for (i = 0; i < ainfo->nseq; i++)
 111     {
 112       if (ainfo->sqinfo[i].flags & SQINFO_SS) free(ainfo->sqinfo[i].ss);
 113       if (ainfo->sqinfo[i].flags & SQINFO_SA) free(ainfo->sqinfo[i].sa);
 114     }
 115   if (ainfo->cs   != NULL) free(ainfo->cs);
 116   if (ainfo->rf   != NULL) free(ainfo->rf);
 117   if (ainfo->name != NULL) free(ainfo->name);
 118   if (ainfo->desc != NULL) free(ainfo->desc);
 119   if (ainfo->acc  != NULL) free(ainfo->acc);
 120   if (ainfo->au   != NULL) free(ainfo->au);
 121
 122   free(ainfo->sqinfo);
 123   free(ainfo->wgt);
 124   Free2DArray((void **) aseqs, ainfo->nseq);
 125 }
 126
 127
 128
 129 /* Function: SAMizeAlignment()
 130  * Date:     SRE, Tue Jun 30 09:49:40 1998 [St. Louis]
 131  *
 132  * Purpose:  Make a "best effort" attempt to convert an alignment
 133  *           to SAM gap format: - in delete col, . in insert col.
 134  *           Only works if alignment adheres to SAM's upper/lower
 135  *           case convention, which is true for instance of old
 136  *           HMMER alignments.
 137  *
 138  * Args:     aseq  - alignment to convert
 139  *           nseq  - number of seqs in alignment
 140  *           alen  - length of alignment
 141  *
 142  * Returns:  (void)
 143  */
 144 void
 145 SAMizeAlignment(char **aseq, int nseq, int alen)
 146 {
 147   int col;                      /* counter for aligned columns */
 148   int i;                        /* counter for seqs */
 149   int sawlower, sawupper, sawgap;
 150   char gapchar;
 151
 152   for (col = 0; col < alen; col++)
 153     {
 154       sawlower = sawupper = sawgap = 0;
 155                                 /* pass 1: do we see only upper or lower? */
 156       for (i = 0; i < nseq; i++)
 157         {
 158           if (isgap(aseq[i][col]))         { sawgap   = 1; continue; }
 159           if (isupper((int) aseq[i][col])) { sawupper = 1; continue; }
 160           if (islower((int) aseq[i][col]))   sawlower = 1;
 161         }
 162                                 /* select gap character for column */
 163       gapchar = '-';            /* default */
 164       if (sawlower && ! sawupper) gapchar = '.';
 165
 166                                 /* pass 2: set gap char */
 167       for (i = 0; i < nseq; i++)
 168         if (isgap(aseq[i][col])) aseq[i][col] = gapchar;
 169     }
 170 }
 171
 172
 173 /* Function: SAMizeAlignmentByGapFrac()
 174  * Date:     SRE, Tue Jun 30 10:58:38 1998 [St. Louis]
 175  *
 176  * Purpose:  Convert an alignment to SAM's gap and case
 177  *           conventions, using gap fraction in a column
 178  *           to choose match versus insert columns. In match columns,
 179  *           residues are upper case and gaps are '-'.
 180  *           In insert columns, residues are lower case and
 181  *           gaps are '.'
 182  *
 183  * Args:     aseq   - aligned sequences
 184  *           nseq   - number of sequences
 185  *           alen   - length of alignment
 186  *           maxgap - if more gaps than this fraction, column is insert.
 187  *
 188  * Returns:  (void) Characters in aseq may be altered.
 189  */
 190 void
 191 SAMizeAlignmentByGapFrac(char **aseq, int nseq, int alen, float maxgap)
 192 {
 193   int apos;                     /* counter over columns */
 194   int idx;                      /* counter over sequences */
 195   int ngap;                     /* number of gaps seen */
 196
 197   for (apos = 0; apos < alen; apos++)
 198     {
 199                                 /* count gaps */
 200       ngap = 0;
 201       for (idx = 0; idx < nseq; idx++)
 202         if (isgap(aseq[idx][apos])) ngap++;
 203
 204                                 /* convert to SAM conventions */
 205       if ((float) ngap / (float) nseq > maxgap)
 206         {                       /* insert column */
 207           for (idx = 0; idx < nseq; idx++)
 208             if (isgap(aseq[idx][apos])) aseq[idx][apos] = '.';
 209             else aseq[idx][apos] = (char) tolower((int) aseq[idx][apos]);
 210         }
 211       else
 212         {                       /* match column */
 213           for (idx = 0; idx < nseq; idx++)
 214             if (isgap(aseq[idx][apos])) aseq[idx][apos] = '-';
 215             else aseq[idx][apos] = (char) toupper((int) aseq[idx][apos]);
 216         }
 217     }
 218 }
 219
 220
 221
 222
 223 /* Function: MakeAlignedString()
 224  *
 225  * Purpose:  Given a raw string of some type (secondary structure, say),
 226  *           align it to a given aseq by putting gaps wherever the
 227  *           aseq has gaps.
 228  *
 229  * Args:     aseq:  template for alignment
 230  *           alen:  length of aseq
 231  *           ss:    raw string to align to aseq
 232  *           ret_s: RETURN: aligned ss
 233  *
 234  * Return:   1 on success, 0 on failure (and squid_errno is set.)
 235  *           ret_ss is malloc'ed here and must be free'd by caller.
 236  */
 237 int
 238 MakeAlignedString(char *aseq, int alen, char *ss, char **ret_s)
 239 {
 240   char *new;
 241   int   apos, rpos;
 242
 243   new = (char *) MallocOrDie ((alen+1) * sizeof(char));
 244   for (apos = rpos = 0; apos < alen; apos++)
 245     if (! isgap(aseq[apos]))
 246       {
 247         new[apos] = ss[rpos];
 248         rpos++;
 249       }
 250     else
 251       new[apos] = '.';
 252   new[apos] = '\0';
 253
 254   if (rpos != strlen(ss))
 255     { squid_errno = SQERR_PARAMETER; free(new); return 0; }
 256   *ret_s = new;
 257   return 1;
 258 }
 259
 260
 261 /* Function: MakeDealignedString()
 262  *
 263  * Purpose:  Given an aligned string of some type (either sequence or
 264  *           secondary structure, for instance), dealign it relative
 265  *           to a given aseq. Return a ptr to the new string.
 266  *
 267  * Args:     aseq  : template alignment
 268  *           alen  : length of aseq
 269  *           ss:   : string to make dealigned copy of; same length as aseq
 270  *           ret_s : RETURN: dealigned copy of ss
 271  *
 272  * Return:   1 on success, 0 on failure (and squid_errno is set)
 273  *           ret_s is alloc'ed here and must be freed by caller
 274  */
 275 int
 276 MakeDealignedString(char *aseq, int alen, char *ss, char **ret_s)
 277 {
 278   char *new;
 279   int   apos, rpos;
 280
 281   new = (char *) MallocOrDie ((alen+1) * sizeof(char));
 282   for (apos = rpos = 0; apos < alen; apos++)
 283     if (! isgap(aseq[apos]))
 284       {
 285         new[rpos] = ss[apos];
 286         rpos++;
 287       }
 288   new[rpos] = '\0';
 289   if (alen != strlen(ss))
 290     { squid_errno = SQERR_PARAMETER; free(new); return 0; }
 291   *ret_s = new;
 292   return 1;
 293 }
 294
 295
 296 /* Function: DealignedLength()
 297  *
 298  * Purpose:  Count the number of non-gap symbols in seq.
 299  *           (i.e. find the length of the unaligned sequence)
 300  *
 301  * Args:     aseq - aligned sequence to count symbols in, \0 terminated
 302  *
 303  * Return:   raw length of seq.
 304  */
 305 int
 306 DealignedLength(char *aseq)
 307 {
 308   int rlen;
 309   for (rlen = 0; *aseq; aseq++)
 310     if (! isgap(*aseq)) rlen++;
 311   return rlen;
 312 }
 313
 314
 315 /* Function: WritePairwiseAlignment()
 316  *
 317  * Purpose:  Write a nice formatted pairwise alignment out,
 318  *           with a BLAST-style middle line showing identities
 319  *           as themselves (single letter) and conservative
 320  *           changes as '+'.
 321  *
 322  * Args:     ofp          - open fp to write to (stdout, perhaps)
 323  *           aseq1, aseq2 - alignments to write (not necessarily
 324  *                          flushed right with gaps)
 325  *           name1, name2 - names of sequences
 326  *           spos1, spos2 - starting position in each (raw) sequence
 327  *           pam          - PAM matrix; positive values define
 328  *                          conservative changes
 329  *           indent       - how many extra spaces to print on left
 330  *
 331  * Return:  1 on success, 0 on failure
 332  */
 333 int
 334 WritePairwiseAlignment(FILE *ofp,
 335                        char *aseq1, char *name1, int spos1,
 336                        char *aseq2, char *name2, int spos2,
 337                        int **pam, int indent)
 338 {
 339   char sname1[11];              /* shortened name               */
 340   char sname2[11];
 341   int  still_going;             /* True if writing another block */
 342   char buf1[61];                /* buffer for writing seq1; CPL+1*/
 343   char bufmid[61];              /* buffer for writing consensus  */
 344   char buf2[61];
 345   char *s1, *s2;                /* ptrs into each sequence          */
 346   int  count1, count2;          /* number of symbols we're writing  */
 347   int  rpos1, rpos2;            /* position in raw seqs             */
 348   int  rawcount1, rawcount2;    /* number of nongap symbols written */
 349   int  apos;
 350
 351   strncpy(sname1, name1, 10);
 352   sname1[10] = '\0';
 353   strtok(sname1, WHITESPACE);
 354
 355   strncpy(sname2, name2, 10);
 356   sname2[10] = '\0';
 357   strtok(sname2, WHITESPACE);
 358
 359   s1 = aseq1;
 360   s2 = aseq2;
 361   rpos1 = spos1;
 362   rpos2 = spos2;
 363
 364   still_going = TRUE;
 365   while (still_going)
 366     {
 367       still_going = FALSE;
 368
 369                                 /* get next line's worth from both */
 370       strncpy(buf1, s1, 60); buf1[60] = '\0';
 371       strncpy(buf2, s2, 60); buf2[60] = '\0';
 372       count1 = strlen(buf1);
 373       count2 = strlen(buf2);
 374
 375                                 /* is there still more to go? */
 376       if ((count1 == 60 && s1[60] != '\0') ||
 377           (count2 == 60 && s2[60] != '\0'))
 378         still_going = TRUE;
 379
 380                                 /* shift seq ptrs by a line */
 381       s1 += count1;
 382       s2 += count2;
 383
 384                                 /* assemble the consensus line */
 385       for (apos = 0; apos < count1 && apos < count2; apos++)
 386         {
 387           if (!isgap(buf1[apos]) && !isgap(buf2[apos]))
 388             {
 389               if (buf1[apos] == buf2[apos])
 390                 bufmid[apos] = buf1[apos];
 391               else if (pam[buf1[apos] - 'A'][buf2[apos] - 'A'] > 0)
 392                 bufmid[apos] = '+';
 393               else
 394                 bufmid[apos] = ' ';
 395             }
 396           else
 397             bufmid[apos] = ' ';
 398         }
 399       bufmid[apos] = '\0';
 400
 401       rawcount1 = 0;
 402       for (apos = 0; apos < count1; apos++)
 403         if (!isgap(buf1[apos])) rawcount1++;
 404
 405       rawcount2 = 0;
 406       for (apos = 0; apos < count2; apos++)
 407         if (!isgap(buf2[apos])) rawcount2++;
 408
 409       (void) fprintf(ofp, "%*s%-10.10s %5d %s %5d\n", indent, "",
 410                      sname1, rpos1, buf1, rpos1 + rawcount1 -1);
 411       (void) fprintf(ofp, "%*s                 %s\n", indent, "",
 412                      bufmid);
 413       (void) fprintf(ofp, "%*s%-10.10s %5d %s %5d\n", indent, "",
 414                      sname2, rpos2, buf2, rpos2 + rawcount2 -1);
 415       (void) fprintf(ofp, "\n");
 416
 417       rpos1 += rawcount1;
 418       rpos2 += rawcount2;
 419     }
 420
 421   return 1;
 422 }
 423
 424
 425 /* Function: MingapAlignment()
 426  *
 427  * Purpose:  Remove all-gap columns from a multiple sequence alignment
 428  *           and its associated data. The alignment is assumed to be
 429  *           flushed (all aseqs the same length).
 430  */
 431 int
 432 MingapAlignment(char **aseqs, AINFO *ainfo)
 433 {
 434   int apos;                     /* position in original alignment */
 435   int mpos;                     /* position in new alignment      */
 436   int idx;
 437
 438   /* We overwrite aseqs, using its allocated memory.
 439    */
 440   for (apos = 0, mpos = 0; aseqs[0][apos] != '\0'; apos++)
 441     {
 442                                 /* check for all-gap in column */
 443       for (idx = 0; idx < ainfo->nseq; idx++)
 444         if (! isgap(aseqs[idx][apos]))
 445           break;
 446       if (idx == ainfo->nseq) continue;
 447
 448                                 /* shift alignment and ainfo */
 449       if (mpos != apos)
 450         {
 451           for (idx = 0; idx < ainfo->nseq; idx++)
 452             aseqs[idx][mpos] = aseqs[idx][apos];
 453
 454           if (ainfo->cs != NULL) ainfo->cs[mpos] = ainfo->cs[apos];
 455           if (ainfo->rf != NULL) ainfo->rf[mpos] = ainfo->rf[apos];
 456         }
 457       mpos++;
 458     }
 459                                 /* null terminate everything */
 460   for (idx = 0; idx < ainfo->nseq; idx++)
 461     aseqs[idx][mpos] = '\0';
 462   ainfo->alen = mpos;   /* set new length */
 463   if (ainfo->cs != NULL) ainfo->cs[mpos] = '\0';
 464   if (ainfo->rf != NULL) ainfo->rf[mpos] = '\0';
 465   return 1;
 466 }
 467
 468
 469
 470 /* Function: RandomAlignment()
 471  *
 472  * Purpose:  Create a random alignment from raw sequences.
 473  *
 474  *           Ideally, we would like to sample an alignment from the
 475  *           space of possible alignments according to its probability,
 476  *           given a prior probability distribution for alignments.
 477  *           I don't see how to describe such a distribution, let alone
 478  *           sample it.
 479  *
 480  *           This is a rough approximation that tries to capture some
 481  *           desired properties. We assume the alignment is generated
 482  *           by a simple HMM composed of match and insert states.
 483  *           Given parameters (pop, pex) for the probability of opening
 484  *           and extending an insertion, we can find the expected number
 485  *           of match states, M, in the underlying model for each sequence.
 486  *           We use an average M taken over all the sequences (this is
 487  *           an approximation. The expectation of M given all the sequence
 488  *           lengths is a nasty-looking summation.)
 489  *
 490  *           M = len / ( 1 + pop ( 1 + 1/ (1-pex) ) )
 491  *
 492  *           Then, we assign positions in each raw sequence onto the M match
 493  *           states and M+1 insert states of this "HMM", by rolling random
 494  *           numbers and inserting the (rlen-M) inserted positions randomly
 495  *           into the insert slots, taking into account the relative probability
 496  *           of open vs. extend.
 497  *
 498  *           The resulting alignment has two desired properties: insertions
 499  *           tend to follow the HMM-like exponential distribution, and
 500  *           the "sparseness" of the alignment is controllable through
 501  *           pop and pex.
 502  *
 503  * Args:     rseqs     - raw sequences to "align", 0..nseq-1
 504  *           sqinfo    - array of 0..nseq-1 info structures for the sequences
 505  *           nseq      - number of sequences
 506  *           pop       - probability to open insertion (0<pop<1)
 507  *           pex       - probability to extend insertion (0<pex<1)
 508  *           ret_aseqs - RETURN: alignment (flushed)
 509  *           ainfo     - fill in: alignment info
 510  *
 511  * Return:   1 on success, 0 on failure. Sets squid_errno to indicate cause
 512  *           of failure.
 513  */
 514 int
 515 RandomAlignment(char **rseqs, SQINFO *sqinfo, int nseq, float pop, float pex,
 516                 char ***ret_aseqs, AINFO *ainfo)
 517 {
 518   char **aseqs;                 /* RETURN: alignment   */
 519   int    alen;                  /* length of alignment */
 520   int   *rlen;                  /* lengths of each raw sequence */
 521   int    M;                     /* length of "model"   */
 522   int  **ins;                   /* insertion counts, 0..nseq-1 by 0..M */
 523   int   *master_ins;            /* max insertion counts, 0..M */
 524   int    apos, rpos, idx;
 525   int    statepos;
 526   int    count;
 527   int    minlen;
 528
 529   /* calculate expected length of model, M
 530    */
 531   rlen = (int *) MallocOrDie (sizeof(int) * nseq);
 532   M = 0;
 533   minlen = 9999999;
 534   for (idx = 0; idx < nseq; idx++)
 535     {
 536       rlen[idx] = strlen(rseqs[idx]);
 537       M += rlen[idx];
 538       minlen = (rlen[idx] < minlen) ? rlen[idx] : minlen;
 539     }
 540   M = (int) ((float) M / (1.0 + pop * (1.0 + 1.0 / (1.0 - pex))));
 541   M /= nseq;
 542   if (M > minlen) M = minlen;
 543
 544   /* make arrays that count insertions in M+1 possible insert states
 545    */
 546   ins = (int **) MallocOrDie (sizeof(int *) * nseq);
 547   master_ins = (int *) MallocOrDie (sizeof(int) * (M+1));
 548   for (idx = 0; idx < nseq; idx++)
 549     {
 550       ins[idx] = (int *) MallocOrDie (sizeof(int) * (M+1));
 551       for (rpos = 0; rpos <= M; rpos++)
 552         ins[idx][rpos] = 0;
 553     }
 554                                 /* normalize */
 555   pop = pop / (pop+pex);
 556   pex = 1.0 - pop;
 557                                 /* make insertions for individual sequences */
 558   for (idx = 0; idx < nseq; idx++)
 559     {
 560       apos = -1;
 561       for (rpos = 0; rpos < rlen[idx]-M; rpos++)
 562         {
 563           if (sre_random() < pop || apos == -1) /* open insertion */
 564             apos = CHOOSE(M+1);        /* choose 0..M */
 565           ins[idx][apos]++;
 566         }
 567     }
 568                                 /* calculate master_ins, max inserts */
 569   alen = M;
 570   for (apos = 0; apos <= M; apos++)
 571     {
 572       master_ins[apos] = 0;
 573       for (idx = 0; idx < nseq; idx++)
 574         if (ins[idx][apos] > master_ins[apos])
 575           master_ins[apos] = ins[idx][apos];
 576       alen += master_ins[apos];
 577     }
 578
 579
 580   /* Now, construct alignment
 581    */
 582   aseqs = (char **) MallocOrDie (sizeof (char *) * nseq);
 583   for (idx = 0; idx < nseq; idx++)
 584     aseqs[idx] = (char *) MallocOrDie (sizeof(char) * (alen+1));
 585   for (idx = 0; idx < nseq; idx++)
 586     {
 587       apos = rpos = 0;
 588
 589       for (statepos = 0; statepos <= M; statepos++)
 590         {
 591           for (count = 0; count < ins[idx][statepos]; count++)
 592             aseqs[idx][apos++] = rseqs[idx][rpos++];
 593           for (; count < master_ins[statepos]; count++)
 594             aseqs[idx][apos++] = ' ';
 595
 596           if (statepos != M)
 597             aseqs[idx][apos++] = rseqs[idx][rpos++];
 598         }
 599       aseqs[idx][alen] = '\0';
 600     }
 601   ainfo->flags = 0;
 602   ainfo->alen  = alen;
 603   ainfo->nseq  = nseq;
 604   ainfo->sqinfo = (SQINFO *) MallocOrDie (sizeof(SQINFO) * nseq);
 605   for (idx = 0; idx < nseq; idx++)
 606     SeqinfoCopy(&(ainfo->sqinfo[idx]), &(sqinfo[idx]));
 607
 608   free(rlen);
 609   free(master_ins);
 610   Free2DArray((void **) ins, nseq);
 611   *ret_aseqs = aseqs;
 612   return 1;
 613 }
 614
 615 /* Function: AlignmentHomogenousGapsym()
 616  * Date:     SRE, Sun Mar 19 19:37:12 2000 [wren, St. Louis]
 617  *
 618  * Purpose:  Sometimes we've got to convert alignments to
 619  *           a lowest common denominator, and we need
 620  *           a single specific gap character -- for example,
 621  *           PSI-BLAST blastpgp -B takes a very simplistic
 622  *           alignment input format which appears to only
 623  *           allow '-' as a gap symbol.
 624  *
 625  *           Anything matching the isgap() macro is
 626  *           converted.
 627  *
 628  * Args:     aseq   - aligned character strings, [0..nseq-1][0..alen-1]
 629  *           nseq   - number of aligned strings
 630  *           alen   - length of alignment
 631  *           gapsym - character to use for gaps.
 632  *
 633  * Returns:  void ("never fails")
 634  */
 635 void
 636 AlignmentHomogenousGapsym(char **aseq, int nseq, int alen, char gapsym)
 637 {
 638   int i, apos;
 639
 640   for (i = 0; i < nseq; i++)
 641     for (apos = 0; apos < alen; apos++)
 642       if (isgap(aseq[i][apos])) aseq[i][apos] = gapsym;
 643 }