1 /*****************************************************************
2 * SQUID - a library of functions for biological sequence analysis
3 * Copyright (C) 1992-2002 Washington University School of Medicine
5 * This source code is freely distributed under the terms of the
6 * GNU General Public License. See the files COPYRIGHT and LICENSE
8 *****************************************************************/
11 * SRE, Sun Jun 6 17:50:45 1999 [bus from Madison, 1999 worm mtg]
13 * Import/export of ClustalV/W multiple sequence alignment
14 * formatted files. Derivative of msf.c; MSF is a pretty
15 * generic interleaved format.
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)
28 /* needed for PACKAGE_VERSION */
29 #include "../config.h"
32 #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
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
43 char *amino_strong[] = {"STA", "NEQK", "NHQK", "NDEQ", "QHRK", "MILV",
44 "MILF", "HY", "FYW", NULL};
46 char *amino_weak[] = {"CSA", "ATV", "SAG", "STNK", "STPA", "SGND",
47 "SNDEQK", "NDEQHK", "NEQHRK", "FVLIM", "HFY", NULL};
51 #ifdef TESTDRIVE_CLUSTAL
52 /*****************************************************************
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
58 main(int argc, char **argv)
66 if ((afp = MSAFileOpen(file, MSAFILE_CLUSTAL, NULL)) == NULL)
67 Die("Couldn't open %s\n", file);
69 while ((msa = ReadClustal(afp)) != NULL)
71 WriteClustal(stdout, msa);
78 /******************************************************************/
79 #endif /* testdrive_clustal */
82 /* Function: ReadClustal()
83 * Date: SRE, Sun Jun 6 17:53:49 1999 [bus from Madison, 1999 worm mtg]
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.
89 * Args: afp - open alignment file
91 * Returns: MSA * - an alignment object
92 * caller responsible for an MSAFree()
93 * NULL if no more alignments
96 * Will Die() here with a (potentially) useful message
97 * if a parsing error occurs.
100 ReadClustal(MSAFILE *afp)
110 if (feof(afp->f)) return NULL;
112 /* Skip until we see the CLUSTAL header
114 while ((s = MSAFileGetLine(afp)) != NULL)
116 if (strncmp(s, "CLUSTAL", 7) == 0 &&
117 strstr(s, "multiple sequence alignment") != NULL)
120 if (s == NULL) return NULL;
122 msa = MSAAlloc(10, 0);
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
129 while ((s = MSAFileGetLine(afp)) != NULL)
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);
135 /* The test for a conservation markup line
137 if (strpbrk(name, ".*:") != NULL && strpbrk(seq, ".*:") != NULL)
140 Die("Parse failed at line %d, file %s: possibly using spaces as gaps",
141 afp->linenumber, afp->fname);
143 /* It's not blank, and it's not a coord line: must be sequence
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);
150 MSAVerifyParse(msa); /* verifies, and also sets alen and wgt. */
155 /* Function: WriteClustal()
156 * Date: SRE, Sun Jun 6 18:12:47 1999 [bus from Madison, worm mtg 1999]
158 * Purpose: Write an alignment in Clustal format to an open file.
160 * Args: fp - file that's open for writing.
161 * msa - alignment to write.
166 WriteClustal(FILE *fp, MSA *msa)
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) */
175 /* consensus line stuff */
184 /* calculate max namelen used */
186 for (idx = 0; idx < msa->nseq; idx++)
187 if ((len = strlen(msa->sqname[idx])) > namelen)
191 fprintf(fp, "CLUSTAL O(%s) multiple sequence alignment\n", PACKAGE_VERSION);
193 fprintf(fp, "CLUSTAL W(1.5) multiple sequence alignment\n");
196 /*****************************************************
197 * Write the sequences
198 *****************************************************/
201 fprintf(fp, "\n"); /* original had two blank lines */
204 for (pos = 0; pos < msa->alen; pos += cpl)
206 fprintf(fp, "\n"); /* Blank line between sequence blocks */
207 for (idx = 0; idx < msa->nseq; idx++)
209 strncpy(buf, msa->aseq[idx] + pos, cpl);
212 fprintf(fp, "%-*s %s\n", namelen+5, msa->sqname[idx], buf);
214 fprintf(fp, "%*s %s\n", namelen, msa->sqname[idx], buf);
218 /* do consensus dots */
220 /* print namelen+5 spaces */
221 for(subpos = 0; subpos <= namelen+5; subpos++)
224 for(subpos = pos; subpos < min(pos + cpl, msa->alen); subpos++)
226 /* see if 100% conservation */
227 first = msa->aseq[0][subpos];
229 for (idx = 1; idx < msa->nseq; idx++)
231 if(msa->aseq[idx][subpos] != first)
241 /* if not then check strong */
242 for(bin = 0; bin < 9; bin++)
243 strong_bins[bin] = 0; /* clear the bins */
245 for(idx = 0; idx < msa->nseq; idx++)
247 switch(msa->aseq[idx][subpos])
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;
269 for(bin = 0; bin < 9; bin++)
270 if(strong_bins[bin] == msa->nseq)
280 for(bin = 0; bin < 11; bin++)
281 weak_bins[bin] = 0; /* clear the bins */
283 for(idx = 0; idx < msa->nseq; idx++)
285 switch(msa->aseq[idx][subpos])
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;
308 for(bin = 0; bin < 11; bin++)
309 if(weak_bins[bin] == msa->nseq)