binaries/src/clustalo/src/squid/clustal.c

   1 /*****************************************************************
   2  * SQUID - a library of functions for biological sequence analysis
   3  * Copyright (C) 1992-2002 Washington University School of Medicine
   4  *
   5  *     This source code is freely distributed under the terms of the
   6  *     GNU General Public License. See the files COPYRIGHT and LICENSE
   7  *     for details.
   8  *****************************************************************/
   9
  10 /* clustal.c
  11  * SRE, Sun Jun  6 17:50:45 1999 [bus from Madison, 1999 worm mtg]
  12  *
  13  * Import/export of ClustalV/W multiple sequence alignment
  14  * formatted files. Derivative of msf.c; MSF is a pretty
  15  * generic interleaved format.
  16  *
  17  * RCS $Id: clustal.c 228 2011-03-29 14:05:27Z dave $ (Original squid RCS Id: clustal.c,v 1.1 1999/07/15 22:26:53 eddy 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 CLUSTALO
  28 /* needed for PACKAGE_VERSION */
  29 #include "../config.h"
  30
  31 #ifndef min
  32         #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
  33 #endif
  34
  35 /*These are all the positively scoring groups that occur in the Gonnet Pam250
  36 matrix. There are strong and weak groups, defined as strong score >0.5 and
  37 weak score =<0.5. Strong matching columns to be assigned ':' and weak matches
  38 assigned '.' in the clustal output format.
  39 amino_strong = res_cat1
  40 amino_weak = res_cat2
  41 */
  42
  43 char *amino_strong[] = {"STA", "NEQK", "NHQK", "NDEQ", "QHRK", "MILV",
  44     "MILF", "HY", "FYW", NULL};
  45
  46 char *amino_weak[] = {"CSA", "ATV", "SAG", "STNK", "STPA", "SGND",
  47     "SNDEQK", "NDEQHK", "NEQHRK", "FVLIM", "HFY", NULL};
  48
  49 #endif
  50
  51 #ifdef TESTDRIVE_CLUSTAL
  52 /*****************************************************************
  53  * msf.c test driver:
  54  * 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
  55  *
  56  */
  57 int
  58 main(int argc, char **argv)
  59 {
  60   MSAFILE *afp;
  61   MSA     *msa;
  62   char    *file;
  63
  64   file = argv[1];
  65
  66   if ((afp = MSAFileOpen(file, MSAFILE_CLUSTAL, NULL)) == NULL)
  67     Die("Couldn't open %s\n", file);
  68
  69   while ((msa = ReadClustal(afp)) != NULL)
  70     {
  71       WriteClustal(stdout, msa);
  72       MSAFree(msa);
  73     }
  74
  75   MSAFileClose(afp);
  76   exit(0);
  77 }
  78 /******************************************************************/
  79 #endif /* testdrive_clustal */
  80
  81
  82 /* Function: ReadClustal()
  83  * Date:     SRE, Sun Jun  6 17:53:49 1999 [bus from Madison, 1999 worm mtg]
  84  *
  85  * Purpose:  Parse an alignment read from an open Clustal format
  86  *           alignment file. Clustal is a single-alignment format.
  87  *           Return the alignment, or NULL if we have no data.
  88  *
  89  * Args:     afp  - open alignment file
  90  *
  91  * Returns:  MSA * - an alignment object
  92  *                   caller responsible for an MSAFree()
  93  *           NULL if no more alignments
  94  *
  95  * Diagnostics:
  96  *           Will Die() here with a (potentially) useful message
  97  *           if a parsing error occurs.
  98  */
  99 MSA *
 100 ReadClustal(MSAFILE *afp)
 101 {
 102   MSA    *msa;
 103   char   *s;
 104   int     slen;
 105   int     sqidx;
 106   char   *name;
 107   char   *seq;
 108   char   *s2;
 109
 110   if (feof(afp->f)) return NULL;
 111
 112   /* Skip until we see the CLUSTAL header
 113    */
 114   while ((s = MSAFileGetLine(afp)) != NULL)
 115     {
 116       if (strncmp(s, "CLUSTAL", 7) == 0 &&
 117           strstr(s, "multiple sequence alignment") != NULL)
 118         break;
 119     }
 120   if (s == NULL) return NULL;
 121
 122   msa = MSAAlloc(10, 0);
 123
 124   /* Now we're in the sequence section.
 125    * As discussed above, if we haven't seen a sequence name, then we
 126    * don't include the sequence in the alignment.
 127    * Watch out for conservation markup lines that contain *.: chars
 128    */
 129   while ((s = MSAFileGetLine(afp)) != NULL)
 130     {
 131       if ((name = sre_strtok(&s, WHITESPACE, NULL))  == NULL) continue;
 132       if ((seq  = sre_strtok(&s, WHITESPACE, &slen)) == NULL) continue;
 133       s2 = sre_strtok(&s, "\n", NULL);
 134
 135       /* The test for a conservation markup line
 136        */
 137       if (strpbrk(name, ".*:") != NULL && strpbrk(seq, ".*:") != NULL)
 138         continue;
 139       if (s2 != NULL)
 140         Die("Parse failed at line %d, file %s: possibly using spaces as gaps",
 141             afp->linenumber, afp->fname);
 142
 143       /* It's not blank, and it's not a coord line: must be sequence
 144        */
 145       sqidx = MSAGetSeqidx(msa, name, msa->lastidx+1);
 146       msa->lastidx = sqidx;
 147       msa->sqlen[sqidx] = sre_strcat(&(msa->aseq[sqidx]), msa->sqlen[sqidx], seq, slen);
 148     }
 149
 150   MSAVerifyParse(msa);          /* verifies, and also sets alen and wgt. */
 151   return msa;
 152 }
 153
 154
 155 /* Function: WriteClustal()
 156  * Date:     SRE, Sun Jun  6 18:12:47 1999 [bus from Madison, worm mtg 1999]
 157  *
 158  * Purpose:  Write an alignment in Clustal format to an open file.
 159  *
 160  * Args:     fp    - file that's open for writing.
 161  *           msa   - alignment to write.
 162  *
 163  * Returns:  (void)
 164  */
 165 void
 166 WriteClustal(FILE *fp, MSA *msa)
 167 {
 168   int    idx;                   /* counter for sequences         */
 169   int    len;                   /* tmp variable for name lengths */
 170   int    namelen;               /* maximum name length used      */
 171   int    pos;                   /* position counter              */
 172   char   buf[80];               /* buffer for writing seq        */
 173   int    cpl = 60;              /* char per line (< 64)          */
 174
 175   /* consensus line stuff */
 176   int subpos;
 177   char first;
 178   int bail;
 179   int strong_bins[9];
 180   int weak_bins[11];
 181   int cons;
 182   int bin;
 183
 184   /* calculate max namelen used */
 185   namelen = 0;
 186   for (idx = 0; idx < msa->nseq; idx++)
 187     if ((len = strlen(msa->sqname[idx])) > namelen)
 188       namelen = len;
 189
 190 #ifdef CLUSTALO
 191   fprintf(fp, "CLUSTAL O(%s) multiple sequence alignment\n", PACKAGE_VERSION);
 192 #else
 193   fprintf(fp, "CLUSTAL W(1.5) multiple sequence alignment\n");
 194 #endif
 195
 196   /*****************************************************
 197    * Write the sequences
 198    *****************************************************/
 199
 200 #ifdef CLUSTALO
 201     fprintf(fp, "\n");  /* original had two blank lines */
 202 #endif
 203
 204   for (pos = 0; pos < msa->alen; pos += cpl)
 205     {
 206       fprintf(fp, "\n");        /* Blank line between sequence blocks */
 207       for (idx = 0; idx < msa->nseq; idx++)
 208       {
 209         strncpy(buf, msa->aseq[idx] + pos, cpl);
 210             buf[cpl] = '\0';
 211 #ifdef CLUSTALO
 212             fprintf(fp, "%-*s %s\n", namelen+5, msa->sqname[idx], buf);
 213 #else
 214             fprintf(fp, "%*s %s\n", namelen, msa->sqname[idx], buf);
 215 #endif
 216           }
 217 #ifdef CLUSTALO
 218       /* do consensus dots */
 219
 220       /* print namelen+5 spaces */
 221       for(subpos = 0; subpos <= namelen+5; subpos++)
 222         fprintf(fp, " ");
 223
 224       for(subpos = pos; subpos < min(pos + cpl, msa->alen); subpos++)
 225       {
 226           /* see if 100% conservation */
 227           first = msa->aseq[0][subpos];
 228           bail = 0;
 229           for (idx = 1; idx < msa->nseq; idx++)
 230           {
 231             if(msa->aseq[idx][subpos] != first)
 232             {
 233               bail = 1;
 234               break;
 235             }
 236           }
 237           if(!bail)
 238             fprintf(fp, "*");
 239           else
 240           {
 241             /* if not then check strong */
 242             for(bin = 0; bin < 9; bin++)
 243               strong_bins[bin] = 0; /* clear the bins */
 244
 245             for(idx = 0; idx < msa->nseq; idx++)
 246             {
 247               switch(msa->aseq[idx][subpos])
 248               {
 249                 case 'S': strong_bins[0]++; break;
 250                 case 'T': strong_bins[0]++; break;
 251                 case 'A': strong_bins[0]++; break;
 252                 case 'N': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; break;
 253                 case 'E': strong_bins[1]++; strong_bins[3]++; break;
 254                 case 'Q': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; strong_bins[4]++; break;
 255                 case 'K': strong_bins[1]++; strong_bins[2]++; strong_bins[4]++; break;
 256                 case 'D': strong_bins[3]++; break;
 257                 case 'R': strong_bins[4]++; break;
 258                 case 'H': strong_bins[4]++; strong_bins[7]++; break;
 259                 case 'M': strong_bins[5]++; strong_bins[6]++; break;
 260                 case 'I': strong_bins[5]++; strong_bins[6]++; break;
 261                 case 'L': strong_bins[5]++; strong_bins[6]++; break;
 262                 case 'V': strong_bins[5]++; break;
 263                 case 'F': strong_bins[6]++; strong_bins[8]++; break;
 264                 case 'Y': strong_bins[7]++; strong_bins[8]++; break;
 265                 case 'W': strong_bins[8]++; break;
 266               }
 267             }
 268             bail = 0;
 269             for(bin = 0; bin < 9; bin++)
 270               if(strong_bins[bin] == msa->nseq)
 271               {
 272                   bail = 1;
 273                   break;
 274               }
 275             if(bail)
 276               fprintf(fp, ":");
 277             else
 278             {
 279               /* check weak */
 280               for(bin = 0; bin < 11; bin++)
 281                 weak_bins[bin] = 0; /* clear the bins */
 282
 283               for(idx = 0; idx < msa->nseq; idx++)
 284               {
 285                 switch(msa->aseq[idx][subpos])
 286                 {
 287                   case 'C': weak_bins[0]++; break;
 288                   case 'S': weak_bins[0]++; weak_bins[2]++; weak_bins[3]++; weak_bins[4]++; weak_bins[5]++; weak_bins[6]++; break;
 289                   case 'A': weak_bins[0]++; weak_bins[1]++; weak_bins[2]++; weak_bins[4]++; break;
 290                   case 'T': weak_bins[1]++; weak_bins[3]++; weak_bins[4]++; break;
 291                   case 'V': weak_bins[1]++; weak_bins[9]++; break;
 292                   case 'G': weak_bins[2]++; break;
 293                   case 'N': weak_bins[3]++; weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
 294                   case 'K': weak_bins[3]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
 295                   case 'D': weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; break;
 296                   case 'E': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
 297                   case 'Q': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
 298                   case 'H': weak_bins[7]++; weak_bins[8]++; weak_bins[10]++; break;
 299                   case 'R': weak_bins[8]++; break;
 300                   case 'F': weak_bins[9]++; weak_bins[10]++; break;
 301                   case 'L': weak_bins[9]++; break;
 302                   case 'I': weak_bins[9]++; break;
 303                   case 'M': weak_bins[9]++; break;
 304                   case 'Y': weak_bins[10]++; break;
 305                 }
 306               }
 307               bail = 0;
 308               for(bin = 0; bin < 11; bin++)
 309                 if(weak_bins[bin] == msa->nseq)
 310                 {
 311                     bail = 1;
 312                     break;
 313                 }
 314               if(bail)
 315                 fprintf(fp, ".");
 316               else
 317                 fprintf(fp, " ");
 318             }
 319           }
 320       }
 321       fprintf(fp,"\n");
 322 #endif
 323     }
 324
 325   return;
 326 }
 327
 328
 329