forester/archive/RIO/others/hmmer/squid/clustal.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 /* clustal.c
  12  * SRE, Sun Jun  6 17:50:45 1999 [bus from Madison, 1999 worm mtg]
  13  *
  14  * Import/export of ClustalV/W multiple sequence alignment
  15  * formatted files. Derivative of msf.c; MSF is a pretty
  16  * generic interleaved format.
  17  *
  18  * RCS $Id: clustal.c,v 1.1.1.1 2005/03/22 08:34:27 cmzmasek Exp $
  19  */
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <ctype.h>
  25 #include "squid.h"
  26 #include "msa.h"
  27
  28 #ifdef TESTDRIVE_CLUSTAL
  29 /*****************************************************************
  30  * msf.c test driver:
  31  * cc -DTESTDRIVE_CLUSTAL -g -O2 -Wall -o test clustal.c msa.c gki.c sqerror.c sre_string.c file.c hsregex.c sre_math.c sre_ctype.c -lm
  32  *
  33  */
  34 int
  35 main(int argc, char **argv)
  36 {
  37   MSAFILE *afp;
  38   MSA     *msa;
  39   char    *file;
  40
  41   file = argv[1];
  42
  43   if ((afp = MSAFileOpen(file, MSAFILE_CLUSTAL, NULL)) == NULL)
  44     Die("Couldn't open %s\n", file);
  45
  46   while ((msa = ReadClustal(afp)) != NULL)
  47     {
  48       WriteClustal(stdout, msa);
  49       MSAFree(msa);
  50     }
  51
  52   MSAFileClose(afp);
  53   exit(0);
  54 }
  55 /******************************************************************/
  56 #endif /* testdrive_clustal */
  57
  58
  59 /* Function: ReadClustal()
  60  * Date:     SRE, Sun Jun  6 17:53:49 1999 [bus from Madison, 1999 worm mtg]
  61  *
  62  * Purpose:  Parse an alignment read from an open Clustal format
  63  *           alignment file. Clustal is a single-alignment format.
  64  *           Return the alignment, or NULL if we have no data.
  65  *
  66  * Args:     afp  - open alignment file
  67  *
  68  * Returns:  MSA * - an alignment object
  69  *                   caller responsible for an MSAFree()
  70  *           NULL if no more alignments
  71  *
  72  * Diagnostics:
  73  *           Will Die() here with a (potentially) useful message
  74  *           if a parsing error occurs.
  75  */
  76 MSA *
  77 ReadClustal(MSAFILE *afp)
  78 {
  79   MSA    *msa;
  80   char   *s;
  81   int     slen;
  82   int     sqidx;
  83   char   *name;
  84   char   *seq;
  85   char   *s2;
  86
  87   if (feof(afp->f)) return NULL;
  88
  89   /* Skip until we see the CLUSTAL header
  90    */
  91   while ((s = MSAFileGetLine(afp)) != NULL)
  92     {
  93       if (strncmp(s, "CLUSTAL", 7) == 0 &&
  94           strstr(s, "multiple sequence alignment") != NULL)
  95         break;
  96     }
  97   if (s == NULL) return NULL;
  98
  99   msa = MSAAlloc(10, 0);
 100
 101   /* Now we're in the sequence section.
 102    * As discussed above, if we haven't seen a sequence name, then we
 103    * don't include the sequence in the alignment.
 104    * Watch out for conservation markup lines that contain *.: chars
 105    */
 106   while ((s = MSAFileGetLine(afp)) != NULL)
 107     {
 108       if ((name = sre_strtok(&s, WHITESPACE, NULL))  == NULL) continue;
 109       if ((seq  = sre_strtok(&s, WHITESPACE, &slen)) == NULL) continue;
 110       s2 = sre_strtok(&s, "\n", NULL);
 111
 112       /* The test for a conservation markup line
 113        */
 114       if (strpbrk(name, ".*:") != NULL && strpbrk(seq, ".*:") != NULL)
 115         continue;
 116       if (s2 != NULL)
 117         Die("Parse failed at line %d, file %s: possibly using spaces as gaps",
 118             afp->linenumber, afp->fname);
 119
 120       /* It's not blank, and it's not a coord line: must be sequence
 121        */
 122       sqidx = MSAGetSeqidx(msa, name, msa->lastidx+1);
 123       msa->lastidx = sqidx;
 124       msa->sqlen[sqidx] = sre_strcat(&(msa->aseq[sqidx]), msa->sqlen[sqidx], seq, slen);
 125     }
 126
 127   MSAVerifyParse(msa);          /* verifies, and also sets alen and wgt. */
 128   return msa;
 129 }
 130
 131
 132 /* Function: WriteClustal()
 133  * Date:     SRE, Sun Jun  6 18:12:47 1999 [bus from Madison, worm mtg 1999]
 134  *
 135  * Purpose:  Write an alignment in Clustal format to an open file.
 136  *
 137  * Args:     fp    - file that's open for writing.
 138  *           msa   - alignment to write.
 139  *
 140  * Returns:  (void)
 141  */
 142 void
 143 WriteClustal(FILE *fp, MSA *msa)
 144 {
 145   int    idx;                   /* counter for sequences         */
 146   int    len;                   /* tmp variable for name lengths */
 147   int    namelen;               /* maximum name length used      */
 148   int    pos;                   /* position counter              */
 149   char   buf[64];               /* buffer for writing seq        */
 150   int    cpl = 50;              /* char per line (< 64)          */
 151
 152                                 /* calculate max namelen used */
 153   namelen = 0;
 154   for (idx = 0; idx < msa->nseq; idx++)
 155     if ((len = strlen(msa->sqname[idx])) > namelen)
 156       namelen = len;
 157
 158   fprintf(fp, "CLUSTAL W(1.5) multiple sequence alignment\n");
 159
 160   /*****************************************************
 161    * Write the sequences
 162    *****************************************************/
 163
 164   for (pos = 0; pos < msa->alen; pos += cpl)
 165     {
 166       fprintf(fp, "\n");        /* Blank line between sequence blocks */
 167       for (idx = 0; idx < msa->nseq; idx++)
 168         {
 169           strncpy(buf, msa->aseq[idx] + pos, cpl);
 170           buf[cpl] = '\0';
 171           fprintf(fp, "%*s %s\n", namelen, msa->sqname[idx], buf);
 172         }
 173     }
 174
 175   return;
 176 }
 177
 178
 179