+++ /dev/null
-/*****************************************************************
- * SQUID - a library of functions for biological sequence analysis
- * Copyright (C) 1992-2002 Washington University School of Medicine
- *
- * This source code is freely distributed under the terms of the
- * GNU General Public License. See the files COPYRIGHT and LICENSE
- * for details.
- *****************************************************************/
-
-/* clustal.c
- * SRE, Sun Jun 6 17:50:45 1999 [bus from Madison, 1999 worm mtg]
- *
- * Import/export of ClustalV/W multiple sequence alignment
- * formatted files. Derivative of msf.c; MSF is a pretty
- * generic interleaved format.
- *
- * 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)
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include "squid.h"
-#include "msa.h"
-
-#ifdef CLUSTALO
-/* needed for PACKAGE_VERSION */
-#include "../config.h"
-
-#ifndef min
- #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
-#endif
-
-/*These are all the positively scoring groups that occur in the Gonnet Pam250
-matrix. There are strong and weak groups, defined as strong score >0.5 and
-weak score =<0.5. Strong matching columns to be assigned ':' and weak matches
-assigned '.' in the clustal output format.
-amino_strong = res_cat1
-amino_weak = res_cat2
-*/
-
-char *amino_strong[] = {"STA", "NEQK", "NHQK", "NDEQ", "QHRK", "MILV",
- "MILF", "HY", "FYW", NULL};
-
-char *amino_weak[] = {"CSA", "ATV", "SAG", "STNK", "STPA", "SGND",
- "SNDEQK", "NDEQHK", "NEQHRK", "FVLIM", "HFY", NULL};
-
-#endif
-
-#ifdef TESTDRIVE_CLUSTAL
-/*****************************************************************
- * msf.c test driver:
- * 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
- *
- */
-int
-main(int argc, char **argv)
-{
- MSAFILE *afp;
- MSA *msa;
- char *file;
-
- file = argv[1];
-
- if ((afp = MSAFileOpen(file, MSAFILE_CLUSTAL, NULL)) == NULL)
- Die("Couldn't open %s\n", file);
-
- while ((msa = ReadClustal(afp)) != NULL)
- {
- WriteClustal(stdout, msa);
- MSAFree(msa);
- }
-
- MSAFileClose(afp);
- exit(0);
-}
-/******************************************************************/
-#endif /* testdrive_clustal */
-
-
-/* Function: ReadClustal()
- * Date: SRE, Sun Jun 6 17:53:49 1999 [bus from Madison, 1999 worm mtg]
- *
- * Purpose: Parse an alignment read from an open Clustal format
- * alignment file. Clustal is a single-alignment format.
- * Return the alignment, or NULL if we have no data.
- *
- * Args: afp - open alignment file
- *
- * Returns: MSA * - an alignment object
- * caller responsible for an MSAFree()
- * NULL if no more alignments
- *
- * Diagnostics:
- * Will Die() here with a (potentially) useful message
- * if a parsing error occurs.
- */
-MSA *
-ReadClustal(MSAFILE *afp)
-{
- MSA *msa;
- char *s;
- int slen;
- int sqidx;
- char *name;
- char *seq;
- char *s2;
-
- if (feof(afp->f)) return NULL;
-
- /* Skip until we see the CLUSTAL header
- */
- while ((s = MSAFileGetLine(afp)) != NULL)
- {
- if (strncmp(s, "CLUSTAL", 7) == 0 &&
- strstr(s, "multiple sequence alignment") != NULL)
- break;
- }
- if (s == NULL) return NULL;
-
- msa = MSAAlloc(10, 0);
-
- /* Now we're in the sequence section.
- * As discussed above, if we haven't seen a sequence name, then we
- * don't include the sequence in the alignment.
- * Watch out for conservation markup lines that contain *.: chars
- */
- while ((s = MSAFileGetLine(afp)) != NULL)
- {
- if ((name = sre_strtok(&s, WHITESPACE, NULL)) == NULL) continue;
- if ((seq = sre_strtok(&s, WHITESPACE, &slen)) == NULL) continue;
- s2 = sre_strtok(&s, "\n", NULL);
-
- /* The test for a conservation markup line
- */
- if (strpbrk(name, ".*:") != NULL && strpbrk(seq, ".*:") != NULL)
- continue;
- if (s2 != NULL)
- Die("Parse failed at line %d, file %s: possibly using spaces as gaps",
- afp->linenumber, afp->fname);
-
- /* It's not blank, and it's not a coord line: must be sequence
- */
- sqidx = MSAGetSeqidx(msa, name, msa->lastidx+1);
- msa->lastidx = sqidx;
- msa->sqlen[sqidx] = sre_strcat(&(msa->aseq[sqidx]), msa->sqlen[sqidx], seq, slen);
- }
-
- MSAVerifyParse(msa); /* verifies, and also sets alen and wgt. */
- return msa;
-}
-
-
-/* Function: WriteClustal()
- * Date: SRE, Sun Jun 6 18:12:47 1999 [bus from Madison, worm mtg 1999]
- *
- * Purpose: Write an alignment in Clustal format to an open file.
- *
- * Args: fp - file that's open for writing.
- * msa - alignment to write.
- *
- * Returns: (void)
- */
-void
-WriteClustal(FILE *fp, MSA *msa)
-{
- int idx; /* counter for sequences */
- int len; /* tmp variable for name lengths */
- int namelen; /* maximum name length used */
- int pos; /* position counter */
- char buf[80]; /* buffer for writing seq */
- int cpl = 60; /* char per line (< 64) */
-
- /* consensus line stuff */
- int subpos;
- char first;
- int bail;
- int strong_bins[9];
- int weak_bins[11];
- int cons;
- int bin;
-
- /* calculate max namelen used */
- namelen = 0;
- for (idx = 0; idx < msa->nseq; idx++)
- if ((len = strlen(msa->sqname[idx])) > namelen)
- namelen = len;
-
-#ifdef CLUSTALO
- fprintf(fp, "CLUSTAL O(%s) multiple sequence alignment\n", PACKAGE_VERSION);
-#else
- fprintf(fp, "CLUSTAL W(1.5) multiple sequence alignment\n");
-#endif
-
- /*****************************************************
- * Write the sequences
- *****************************************************/
-
-#ifdef CLUSTALO
- fprintf(fp, "\n"); /* original had two blank lines */
-#endif
-
- for (pos = 0; pos < msa->alen; pos += cpl)
- {
- fprintf(fp, "\n"); /* Blank line between sequence blocks */
- for (idx = 0; idx < msa->nseq; idx++)
- {
- strncpy(buf, msa->aseq[idx] + pos, cpl);
- buf[cpl] = '\0';
-#ifdef CLUSTALO
- fprintf(fp, "%-*s %s\n", namelen+5, msa->sqname[idx], buf);
-#else
- fprintf(fp, "%*s %s\n", namelen, msa->sqname[idx], buf);
-#endif
- }
-#ifdef CLUSTALO
- /* do consensus dots */
-
- /* print namelen+5 spaces */
- for(subpos = 0; subpos <= namelen+5; subpos++)
- fprintf(fp, " ");
-
- for(subpos = pos; subpos < min(pos + cpl, msa->alen); subpos++)
- {
- /* see if 100% conservation */
- first = msa->aseq[0][subpos];
- bail = 0;
- for (idx = 1; idx < msa->nseq; idx++)
- {
- if(msa->aseq[idx][subpos] != first)
- {
- bail = 1;
- break;
- }
- }
- if(!bail)
- fprintf(fp, "*");
- else
- {
- /* if not then check strong */
- for(bin = 0; bin < 9; bin++)
- strong_bins[bin] = 0; /* clear the bins */
-
- for(idx = 0; idx < msa->nseq; idx++)
- {
- switch(msa->aseq[idx][subpos])
- {
- case 'S': strong_bins[0]++; break;
- case 'T': strong_bins[0]++; break;
- case 'A': strong_bins[0]++; break;
- case 'N': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; break;
- case 'E': strong_bins[1]++; strong_bins[3]++; break;
- case 'Q': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; strong_bins[4]++; break;
- case 'K': strong_bins[1]++; strong_bins[2]++; strong_bins[4]++; break;
- case 'D': strong_bins[3]++; break;
- case 'R': strong_bins[4]++; break;
- case 'H': strong_bins[4]++; strong_bins[7]++; break;
- case 'M': strong_bins[5]++; strong_bins[6]++; break;
- case 'I': strong_bins[5]++; strong_bins[6]++; break;
- case 'L': strong_bins[5]++; strong_bins[6]++; break;
- case 'V': strong_bins[5]++; break;
- case 'F': strong_bins[6]++; strong_bins[8]++; break;
- case 'Y': strong_bins[7]++; strong_bins[8]++; break;
- case 'W': strong_bins[8]++; break;
- }
- }
- bail = 0;
- for(bin = 0; bin < 9; bin++)
- if(strong_bins[bin] == msa->nseq)
- {
- bail = 1;
- break;
- }
- if(bail)
- fprintf(fp, ":");
- else
- {
- /* check weak */
- for(bin = 0; bin < 11; bin++)
- weak_bins[bin] = 0; /* clear the bins */
-
- for(idx = 0; idx < msa->nseq; idx++)
- {
- switch(msa->aseq[idx][subpos])
- {
- case 'C': weak_bins[0]++; break;
- case 'S': weak_bins[0]++; weak_bins[2]++; weak_bins[3]++; weak_bins[4]++; weak_bins[5]++; weak_bins[6]++; break;
- case 'A': weak_bins[0]++; weak_bins[1]++; weak_bins[2]++; weak_bins[4]++; break;
- case 'T': weak_bins[1]++; weak_bins[3]++; weak_bins[4]++; break;
- case 'V': weak_bins[1]++; weak_bins[9]++; break;
- case 'G': weak_bins[2]++; break;
- case 'N': weak_bins[3]++; weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
- case 'K': weak_bins[3]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
- case 'D': weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; break;
- case 'E': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
- case 'Q': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
- case 'H': weak_bins[7]++; weak_bins[8]++; weak_bins[10]++; break;
- case 'R': weak_bins[8]++; break;
- case 'F': weak_bins[9]++; weak_bins[10]++; break;
- case 'L': weak_bins[9]++; break;
- case 'I': weak_bins[9]++; break;
- case 'M': weak_bins[9]++; break;
- case 'Y': weak_bins[10]++; break;
- }
- }
- bail = 0;
- for(bin = 0; bin < 11; bin++)
- if(weak_bins[bin] == msa->nseq)
- {
- bail = 1;
- break;
- }
- if(bail)
- fprintf(fp, ".");
- else
- fprintf(fp, " ");
- }
- }
- }
- fprintf(fp,"\n");
-#endif
- }
-
- return;
-}
-
-
-
git://source.jalview.org / jabaws.git / blobdiff