forester/archive/RIO/others/hmmer/squid/phylip.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 /* phylip.c
  12  * SRE, Mon Jun 14 14:08:33 1999 [St. Louis]
  13  *
  14  * Import/export of PHYLIP interleaved multiple sequence alignment
  15  * format files.
  16  *
  17  * RCS $Id: phylip.c,v 1.1.1.1 2005/03/22 08:34:25 cmzmasek Exp $
  18  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <ctype.h>
  24 #include "squid.h"
  25 #include "msa.h"
  26
  27 #ifdef TESTDRIVE_PHYLIP
  28 /*****************************************************************
  29  * phylip.c test driver:
  30  *
  31  */
  32 int
  33 main(int argc, char **argv)
  34 {
  35   MSAFILE *afp;
  36   MSA     *msa;
  37   char    *file;
  38
  39   file = argv[1];
  40
  41   if ((afp = MSAFileOpen(file, MSAFILE_UNKNOWN, NULL)) == NULL)
  42     Die("Couldn't open %s\n", file);
  43
  44   printf("format %d\n", afp->format);
  45
  46   while ((msa = ReadPhylip(afp)) != NULL)
  47     {
  48       WritePhylip(stdout, msa);
  49       MSAFree(msa);
  50     }
  51
  52   MSAFileClose(afp);
  53   exit(0);
  54 }
  55 /******************************************************************/
  56 #endif /* testdrive_phylip */
  57
  58
  59
  60 /* Function: ReadPhylip()
  61  * Date:     SRE, Fri Jun 18 12:59:37 1999 [Sanger Centre]
  62  *
  63  * Purpose:  Parse an alignment from an open Phylip format
  64  *           alignment file. Phylip is a single-alignment format.
  65  *           Return the alignment, or NULL if we have no data.
  66  *
  67  * Args:     afp - open alignment file
  68  *
  69  * Returns:  MSA * - an alignment object
  70  *                   Caller responsible for an MSAFree()
  71  *           NULL if no more alignments
  72  */
  73 MSA *
  74 ReadPhylip(MSAFILE *afp)
  75 {
  76   MSA  *msa;
  77   char *s, *s1, *s2;
  78   char  name[11];               /* seq name max len = 10 char */
  79   int   nseq, alen;
  80   int   idx;                    /* index of current sequence */
  81   int   slen;
  82   int   nblock;
  83
  84   if (feof(afp->f)) return NULL;
  85
  86   /* Skip until we see a nonblank line; it's the header,
  87    * containing nseq/alen
  88    */
  89   nseq = 0; alen = 0;
  90   while ((s = MSAFileGetLine(afp)) != NULL)
  91     {
  92       if ((s1 = sre_strtok(&s, WHITESPACE, NULL)) == NULL) continue;
  93       if ((s2 = sre_strtok(&s, WHITESPACE, NULL)) == NULL)
  94         Die("Failed to parse nseq/alen from first line of PHYLIP file %s\n", afp->fname);
  95       if (! IsInt(s1) || ! IsInt(s2))
  96         Die("nseq and/or alen not an integer in first line of PHYLIP file %s\n", afp->fname);
  97       nseq = atoi(s1);
  98       alen = atoi(s2);
  99       break;
 100     }
 101
 102   msa = MSAAlloc(nseq, 0);
 103   idx    = 0;
 104   nblock = 0;
 105   while ((s = MSAFileGetLine(afp)) != NULL)
 106     {
 107       /* ignore blank lines. nonblank lines start w/ nonblank char */
 108       if (isspace(*s)) continue;
 109                                 /* First block has seq names */
 110       if (nblock == 0) {
 111         strncpy(name, s, 10);
 112         name[10] = '\0';
 113         GKIStoreKey(msa->index, name);
 114         msa->sqname[idx] = sre_strdup(name, -1);
 115         s += 10;
 116       }
 117                                 /* be careful of trailing whitespace on lines */
 118       if ((s1 = sre_strtok(&s, WHITESPACE, &slen)) == NULL)
 119         Die("Failed to parse sequence at line %d of PHYLIP file %s\n",
 120             afp->linenumber, afp->fname);
 121       msa->sqlen[idx] = sre_strcat(&(msa->aseq[idx]), msa->sqlen[idx], s1, slen);
 122
 123       idx++;
 124       if (idx == nseq) { idx = 0; nblock++; }
 125     }
 126   msa->nseq = nseq;
 127   MSAVerifyParse(msa);          /* verifies; sets alen, wgt; frees sqlen[] */
 128   return msa;
 129 }
 130
 131
 132
 133 /* Function: WritePhylip()
 134  * Date:     SRE, Fri Jun 18 12:07:41 1999 [Sanger Centre]
 135  *
 136  * Purpose:  Write an alignment in Phylip format to an open file.
 137  *
 138  * Args:     fp    - file that's open for writing.
 139  *           msa   - alignment to write.
 140  *
 141  * Returns:  (void)
 142  */
 143 void
 144 WritePhylip(FILE *fp, MSA *msa)
 145 {
 146   int    idx;                   /* counter for sequences         */
 147   int    cpl = 50;              /* 50 seq char per line          */
 148   char   buf[51];               /* buffer for writing seq        */
 149   int    pos;
 150
 151   /* First line has nseq, alen
 152    */
 153   fprintf(fp, " %d  %d\n", msa->nseq, msa->alen);
 154
 155   /* Alignment section.
 156    * PHYLIP is a multiblock format, blocks (optionally) separated
 157    * by blanks; names only attached to first block. Names are
 158    * restricted to ten char; we achieve this by simple truncation (!).
 159    * (Do we need to convert gap characters from our ./- convention?)
 160    */
 161   for (pos = 0; pos < msa->alen; pos += cpl)
 162     {
 163       if (pos > 0) fprintf(fp, "\n");
 164
 165       for (idx = 0; idx < msa->nseq; idx++)
 166         {
 167           strncpy(buf, msa->aseq[idx] + pos, cpl);
 168           buf[cpl] = '\0';
 169           if (pos > 0) fprintf(fp, "%s\n", buf);
 170           else         fprintf(fp, "%-10.10s%s\n", msa->sqname[idx], buf);
 171         }
 172     }
 173   return;
 174 }