+++ /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.
- *****************************************************************/
-
-/* msf.c
- * SRE, Sun Jul 11 16:17:32 1993
- *
- * Import/export of GCG MSF multiple sequence alignment
- * formatted files. Designed using format specifications
- * kindly provided by Steve Smith of Genetics Computer Group.
- *
- * RCS $Id: msf.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: msf.c,v 1.4 2001/04/23 00:35:33 eddy Exp)
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include <time.h>
-#include "squid.h"
-#include "msa.h"
-
-#ifdef TESTDRIVE_MSF
-/*****************************************************************
- * msf.c test driver:
- * cc -DTESTDRIVE_MSF -g -O2 -Wall -o test msf.c msa.c gki.c sqerror.c sre_string.c file.c hsregex.c sre_math.c sre_ctype.c sqio.c alignio.c selex.c interleaved.c types.c -lm
- *
- */
-int
-main(int argc, char **argv)
-{
- MSAFILE *afp;
- MSA *msa;
- char *file;
-
- file = argv[1];
-
- if ((afp = MSAFileOpen(file, MSAFILE_STOCKHOLM, NULL)) == NULL)
- Die("Couldn't open %s\n", file);
-
- while ((msa = ReadMSF(afp)) != NULL)
- {
- WriteMSF(stdout, msa);
- MSAFree(msa);
- }
-
- MSAFileClose(afp);
- exit(0);
-}
-/******************************************************************/
-#endif /* testdrive_msf */
-
-
-
-/* Function: ReadMSF()
- * Date: SRE, Tue Jun 1 08:07:22 1999 [St. Louis]
- *
- * Purpose: Parse an alignment read from an open MSF format
- * alignment file. (MSF is a single-alignment format.)
- * Return the alignment, or NULL if we've already
- * read the alignment.
- *
- * 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 *
-ReadMSF(MSAFILE *afp)
-{
- MSA *msa;
- char *s;
- int alleged_alen;
- int alleged_type;
- int alleged_checksum;
- char *tok;
- char *sp;
- int slen;
- int sqidx;
- char *name;
- char *seq;
-
- if (feof(afp->f)) return NULL;
- if ((s = MSAFileGetLine(afp)) == NULL) return NULL;
-
- /* The first line is the header.
- * This is a new-ish GCG feature. Don't count on it, so
- * we can be a bit more tolerant towards non-GCG software
- * generating "MSF" files.
- */
- msa = MSAAlloc(10, 0);
- if (strncmp(s, "!!AA_MULTIPLE_ALIGNMENT", 23) == 0) {
- msa->type = kAmino;
- if ((s = MSAFileGetLine(afp)) == NULL) return NULL;
- } else if (strncmp(s, "!!NA_MULTIPLE_ALIGNMENT", 23) == 0) {
- msa->type = kRNA;
- if ((s = MSAFileGetLine(afp)) == NULL) return NULL;
- }
-
- /* Now we're in the free text comment section of the MSF file.
- * It ends when we see the "MSF: Type: Check: .." line.
- * This line must be present.
- */
- do
- {
- if ((strstr(s, "..") != NULL && strstr(s, "MSF:") != NULL) &&
- Strparse("^.+MSF: +([0-9]+) +Type: +([PNX]).+Check: +([0-9]+) +\\.\\.", s, 3))
- {
- alleged_alen = atoi(sqd_parse[0]);
- switch (*(sqd_parse[1])) {
- case 'N' : alleged_type = kRNA; break;
- case 'P' : alleged_type = kAmino; break;
- case 'X' : alleged_type = kOtherSeq; break;
- default : alleged_type = kOtherSeq;
- }
- alleged_checksum = atoi(sqd_parse[3]);
- if (msa->type == kOtherSeq) msa->type = alleged_type;
- break; /* we're done with comment section. */
- }
- if (! IsBlankline(s))
- MSAAddComment(msa, s);
- } while ((s = MSAFileGetLine(afp)) != NULL);
-
- /* Now we're in the name section.
- * GCG has a relatively poorly documented feature: only sequences that
- * appear in this list will be read from the alignment section. Commenting
- * out sequences in the name list (by preceding them with "!") is
- * allowed as a means of manually defining subsets of sequences in
- * the alignment section. We can support this feature reasonably
- * easily because of the hash table for names in the MSA: we
- * only add names to the hash table when we see 'em in the name section.
- */
- while ((s = MSAFileGetLine(afp)) != NULL)
- {
- while ((*s == ' ' || *s == '\t') && *s) s++; /* skip leading whitespace */
-
- if (*s == '\n') continue; /* skip blank lines */
- else if (*s == '!') MSAAddComment(msa, s);
- else if ((sp = strstr(s, "Name:")) != NULL)
- {
- /* We take the name and the weigh, and that's it */
- sp += 5;
- tok = sre_strtok(&sp, " \t", &slen); /* <sequence name> */
- sqidx = GKIStoreKey(msa->index, tok);
- if (sqidx >= msa->nseqalloc) MSAExpand(msa);
- msa->sqname[sqidx] = sre_strdup(tok, slen);
- msa->nseq++;
-
- if ((sp = strstr(sp, "Weight:")) == NULL)
- Die("No Weight: on line %d for %s in name section of MSF file %s\n",
- afp->linenumber, msa->sqname[sqidx], afp->fname);
- sp += 7;
- tok = sre_strtok(&sp, " \t", &slen);
- msa->wgt[sqidx] = atof(tok);
- msa->flags |= MSA_SET_WGT;
- }
- else if (strncmp(s, "//", 2) == 0)
- break;
- else
- {
- Die("Invalid line (probably %d) in name section of MSF file %s:\n%s\n",
- afp->linenumber, afp->fname, s);
- squid_errno = SQERR_FORMAT; /* NOT THREADSAFE */
- return NULL;
- }
-
- }
-
- /* And 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.
- * Also, watch out for coordinate-only lines.
- */
- while ((s = MSAFileGetLine(afp)) != NULL)
- {
- sp = s;
- if ((name = sre_strtok(&sp, " \t", NULL)) == NULL) continue;
- if ((seq = sre_strtok(&sp, "\n", &slen)) == NULL) continue;
-
- /* The test for a coord line: digits starting both fields
- */
- if (isdigit(*name) && isdigit(*seq))
- continue;
-
- /* It's not blank, and it's not a coord line: must be sequence
- */
- sqidx = GKIKeyIndex(msa->index, name);
- if (sqidx < 0) continue; /* not a sequence we recognize */
-
- msa->sqlen[sqidx] = sre_strcat(&(msa->aseq[sqidx]), msa->sqlen[sqidx], seq, slen);
- }
-
- /* We've left blanks in the aseqs; take them back out.
- */
- for (sqidx = 0; sqidx < msa->nseq; sqidx++)
- {
- if (msa->aseq[sqidx] == NULL)
- Die("Didn't find a sequence for %s in MSF file %s\n", msa->sqname[sqidx], afp->fname);
-
- for (s = sp = msa->aseq[sqidx]; *s != '\0'; s++)
- {
- if (*s == ' ' || *s == '\t') {
- msa->sqlen[sqidx]--;
- } else {
- *sp = *s;
- sp++;
- }
- }
- *sp = '\0';
- }
-
- MSAVerifyParse(msa); /* verifies, and also sets alen and wgt. */
- return msa;
-}
-
-
-/* Function: WriteMSF()
- * Date: SRE, Mon May 31 11:25:18 1999 [St. Louis]
- *
- * Purpose: Write an alignment in MSF format to an open file.
- *
- * Args: fp - file that's open for writing.
- * msa - alignment to write.
- *
- * Note that msa->type, usually optional, must be
- * set for WriteMSF to work. If it isn't, a fatal
- * error is generated.
- *
- * Returns: (void)
- */
-void
-WriteMSF(FILE *fp, MSA *msa)
-{
- time_t now; /* current time as a time_t */
- char date[64]; /* today's date in GCG's format "October 3, 1996 15:57" */
- char **gcg_aseq; /* aligned sequences with gaps converted to GCG format */
- char **gcg_sqname; /* sequence names with GCG-valid character sets */
- int idx; /* counter for sequences */
- char *s; /* pointer into sqname or seq */
- int len; /* tmp variable for name lengths */
- int namelen; /* maximum name length used */
- int pos; /* position counter */
- char buffer[51]; /* buffer for writing seq */
- int i; /* another position counter */
-
- /*****************************************************************
- * Make copies of sequence names and sequences.
- * GCG recommends that name characters should only contain
- * alphanumeric characters, -, or _
- * Some GCG and GCG-compatible software is sensitive to this.
- * We silently convert all other characters to '_'.
- *
- * For sequences, GCG allows only ~ and . for gaps.
- * Otherwise, everthing is interpreted as a residue;
- * so squid's IUPAC-restricted chars are fine. ~ means
- * an external gap. . means an internal gap.
- *****************************************************************/
-
- /* make copies that we can edit */
- gcg_aseq = MallocOrDie(sizeof(char *) * msa->nseq);
- gcg_sqname = MallocOrDie(sizeof(char *) * msa->nseq);
- for (idx = 0; idx < msa->nseq; idx++)
- {
- gcg_aseq[idx] = sre_strdup(msa->aseq[idx], msa->alen);
- gcg_sqname[idx] = sre_strdup(msa->sqname[idx], -1);
- }
- /* alter names as needed */
- for (idx = 0; idx < msa->nseq; idx++)
- for (s = gcg_sqname[idx]; *s != '\0'; s++)
- if (! isalnum((int) *s) && *s != '-' && *s != '_')
- *s = '_';
- /* alter gap chars in seq */
- for (idx = 0; idx < msa->nseq; idx++)
- {
- for (s = gcg_aseq[idx]; *s != '\0' && isgap(*s); s++)
- *s = '~';
- for (; *s != '\0'; s++)
- if (isgap(*s)) *s = '.';
- for (pos = msa->alen-1; pos > 0 && isgap(gcg_aseq[idx][pos]); pos--)
- gcg_aseq[idx][pos] = '~';
- }
- /* calculate max namelen used */
- namelen = 0;
- for (idx = 0; idx < msa->nseq; idx++)
- if ((len = strlen(msa->sqname[idx])) > namelen)
- namelen = len;
-
- /*****************************************************
- * Write the MSF header
- *****************************************************/
- /* required file type line */
- if (msa->type == kOtherSeq)
- msa->type = GuessAlignmentSeqtype(msa->aseq, msa->nseq);
-
- if (msa->type == kRNA) fprintf(fp, "!!NA_MULTIPLE_ALIGNMENT 1.0\n");
- else if (msa->type == kDNA) fprintf(fp, "!!NA_MULTIPLE_ALIGNMENT 1.0\n");
- else if (msa->type == kAmino) fprintf(fp, "!!AA_MULTIPLE_ALIGNMENT 1.0\n");
- else if (msa->type == kOtherSeq)
- Die("WriteMSF(): couldn't guess whether that alignment is RNA or protein.\n");
- else
- Die("Invalid sequence type %d in WriteMSF()\n", msa->type);
-
- /* free text comments */
- if (msa->ncomment > 0)
- {
- for (idx = 0; idx < msa->ncomment; idx++)
- fprintf(fp, "%s\n", msa->comment[idx]);
- fprintf(fp, "\n");
- }
- /* required checksum line */
- now = time(NULL);
- if (strftime(date, 64, "%B %d, %Y %H:%M", localtime(&now)) == 0)
- Die("What time is it on earth? strftime() failed in WriteMSF().\n");
- fprintf(fp, " %s MSF: %d Type: %c %s Check: %d ..\n",
- msa->name != NULL ? msa->name : "squid.msf",
- msa->alen,
- msa->type == kRNA ? 'N' : 'P',
- date,
- GCGMultchecksum(gcg_aseq, msa->nseq));
- fprintf(fp, "\n");
-
- /*****************************************************
- * Names/weights section
- *****************************************************/
-
- for (idx = 0; idx < msa->nseq; idx++)
- {
- fprintf(fp, " Name: %-*.*s Len: %5d Check: %4d Weight: %.2f\n",
- namelen, namelen,
- gcg_sqname[idx],
- msa->alen,
- GCGchecksum(gcg_aseq[idx], msa->alen),
- msa->wgt[idx]);
- }
- fprintf(fp, "\n");
- fprintf(fp, "//\n");
-
- /*****************************************************
- * Write the sequences
- *****************************************************/
-
- for (pos = 0; pos < msa->alen; pos += 50)
- {
- fprintf(fp, "\n"); /* Blank line between sequence blocks */
-
- /* Coordinate line */
- len = (pos + 50) > msa->alen ? msa->alen - pos : 50;
- if (len > 10)
- fprintf(fp, "%*s %-6d%*s%6d\n", namelen, "",
- pos+1,
- len + ((len-1)/10) - 12, "",
- pos + len);
- else
- fprintf(fp, "%*s %-6d\n", namelen, "", pos+1);
-
- for (idx = 0; idx < msa->nseq; idx++)
- {
- fprintf(fp, "%-*s ", namelen, gcg_sqname[idx]);
- /* get next line's worth of 50 from seq */
- strncpy(buffer, gcg_aseq[idx] + pos, 50);
- buffer[50] = '\0';
- /* draw the sequence line */
- for (i = 0; i < len; i++)
- {
- if (! (i % 10)) fputc(' ', fp);
- fputc(buffer[i], fp);
- }
- fputc('\n', fp);
- }
- }
-
- Free2DArray((void **) gcg_aseq, msa->nseq);
- Free2DArray((void **) gcg_sqname, msa->nseq);
- return;
-}
-
-
-
git://source.jalview.org / jabaws.git / blobdiff