forester/archive/RIO/others/hmmer/squid/a2m.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 /* a2m.c
  12  *
  13  * reading/writing A2M (aligned FASTA) files.
  14  *
  15  * RCS $Id: a2m.c,v 1.1.1.1 2005/03/22 08:34:17 cmzmasek Exp $
  16  */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include "squid.h"
  22 #include "msa.h"
  23
  24 /* Function: ReadA2M()
  25  * Date:     SRE, Sun Jun  6 17:11:29 1999 [bus from Madison 1999 worm mtg]
  26  *
  27  * Purpose:  Parse an alignment read from an open A2M format
  28  *           alignment file. A2M is a single alignment format.
  29  *           Return the alignment, or NULL if we've already
  30  *           read the alignment.
  31  *
  32  * Args:     afp - open alignment file
  33  *
  34  * Returns:  MSA *  - an alignment object.
  35  *                    Caller responsible for an MSAFree()
  36  */
  37 MSA *
  38 ReadA2M(MSAFILE *afp)
  39 {
  40   MSA  *msa;
  41   char *buf;
  42   char *name;
  43   char *desc;
  44   char *seq;
  45   int   idx;
  46   int   len1, len2;
  47
  48   if (feof(afp->f)) return NULL;
  49
  50   name = NULL;
  51   msa  = MSAAlloc(10, 0);
  52   idx  = 0;
  53   while ((buf = MSAFileGetLine(afp)) != NULL)
  54     {
  55       if (*buf == '>')
  56         {
  57           buf++;                /* skip the '>' */
  58           if ((name = sre_strtok(&buf, WHITESPACE, &len1)) == NULL)
  59             Die("Blank name in A2M file %s (line %d)\n", afp->fname, afp->linenumber);
  60           desc = sre_strtok(&buf, "\n", &len2);
  61
  62           idx = GKIStoreKey(msa->index, name);
  63           if (idx >= msa->nseqalloc) MSAExpand(msa);
  64
  65           msa->sqname[idx] = sre_strdup(name, len1);
  66           if (desc != NULL) MSASetSeqDescription(msa, idx, desc);
  67           msa->nseq++;
  68         }
  69       else if (name != NULL)
  70         {
  71           if ((seq = sre_strtok(&buf, WHITESPACE, &len1)) == NULL) continue;
  72           msa->sqlen[idx] = sre_strcat(&(msa->aseq[idx]), msa->sqlen[idx], seq, len1);
  73         }
  74     }
  75   if (name == NULL) { MSAFree(msa); return NULL; }
  76
  77   MSAVerifyParse(msa);
  78   return msa;
  79 }
  80
  81
  82 /* Function: WriteA2M()
  83  * Date:     SRE, Sun Jun  6 17:40:35 1999 [bus from Madison, 1999 worm mtg]
  84  *
  85  * Purpose:  Write an "aligned FASTA" (aka a2m, to UCSC) formatted
  86  *           alignment.
  87  *
  88  * Args:     fp    - open FILE to write to.
  89  *           msa   - alignment to write
  90  *
  91  * Returns:  void
  92  */
  93 void
  94 WriteA2M(FILE *fp, MSA *msa)
  95 {
  96   int  idx;                     /* sequence index */
  97   int  pos;                     /* position in sequence */
  98   char buf[64];                 /* buffer for individual lines */
  99   int  cpl = 60;                /* char per line; must be < 64 unless buf is bigger */
 100
 101   buf[cpl] = '\0';
 102   for (idx = 0; idx < msa->nseq; idx++)
 103     {
 104       fprintf(fp, ">%s %s\n",
 105               msa->sqname[idx],
 106               (msa->sqdesc != NULL && msa->sqdesc[idx] != NULL) ? msa->sqdesc[idx] : "");
 107       for (pos = 0; pos < msa->alen; pos+=cpl)
 108         {
 109           strncpy(buf, &(msa->aseq[idx][pos]), cpl);
 110           fprintf(fp, "%s\n", buf);
 111         }
 112     }
 113 }