+++ /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.
- *****************************************************************/
-
-/* selex.c
- *
- * SRE, Mon Jun 14 11:08:38 1999
- * SELEX obsolete as the preferred HMMER/SQUID format
- * replaced by Stockholm format
- * selex support retained for backwards compatibility
- * kludged to use the MSA interface
- *
- * SRE, Mon Jan 30 14:41:49 1995:
- * #=SA side chain % surface accessibility annotation supported
- *
- * SRE, Tue Nov 9 17:40:50 1993:
- * major revision. #= special comments and aliinfo_s optional
- * alignment info support added. Support for #=CS (consensus
- * secondary structure), #=SS (individual secondary structure),
- * #=RF (reference coordinate system), #=SQ (per-sequence header info),
- * and #=AU ("author") added.
- *
- * Fri Dec 4 17:43:24 1992, SRE:
- * Reading and writing aligned sequences to/from disk files.
- * Implements a new, broader specification of SELEX format
- * and supercedes alignio.c.
- *
- * SELEX format is documented in Docs/formats.tex.
- ****************************************************************************
- * RCS $Id: selex.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: selex.c,v 1.11 2002/10/12 04:40:35 eddy Exp)
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include <memory.h>
-#include "squid.h"
-#include "msa.h"
-
-static int copy_alignment_line(char *aseq, int apos, int name_rcol,
- char *buffer, int lcol, int rcol, char gapsym);
-static void actually_write_selex(FILE *fp, MSA *msa, int cpl);
-
-static char commentsyms[] = "%#";
-
-/* Function: ReadSELEX()
- * Date: SRE, Sun Jun 6 18:24:09 1999 [St. Louis]
- *
- * Purpose: Parse an alignment read from an open SELEX format
- * alignment file. (SELEX is a single alignment format).
- * Return the alignment, or NULL if we've already read the
- * alignment or there's no alignment data in the file.
- *
- * Limitations: SELEX is the only remaining multipass parser for
- * alignment files. It cannot read from gzip or from stdin.
- * It Die()'s here if you try. The reason for this
- * that SELEX allows space characters as gaps, so we don't
- * know the borders of an alignment block until we've seen
- * the whole block. I could rewrite to allow single-pass
- * parsing (by storing the whole block in memory) but
- * since SELEX is now legacy, why bother.
- *
- * Note that the interface is totally kludged: fastest
- * possible adaptation of old ReadSELEX() to the new
- * MSA interface.
- *
- * Args: afp - open alignment file
- *
- * Returns: MSA * - an alignment object
- * caller responsible for an MSAFree()
- * NULL if no alignment data.
- */
-MSA *
-ReadSELEX(MSAFILE *afp)
-{
- MSA *msa; /* RETURN: mult seq alignment */
- FILE *fp; /* ptr to opened seqfile */
- char **aseqs; /* aligned seqs */
- int num = 0; /* number of seqs read */
- char buffer[LINEBUFLEN]; /* input buffer for lines */
- char bufcpy[LINEBUFLEN]; /* strtok'able copy of buffer */
- struct block_struc { /** alignment data for a block: */
- int lcol; /* furthest left aligned sym */
- int rcol; /* furthest right aligned sym */
- } *blocks = NULL;
- int blocknum; /* number of blocks in file */
- char *nptr; /* ptr to start of name on line */
- char *sptr; /* ptr into sequence on line */
- int currnum; /* num. seqs in given block */
- int currblock; /* index for blocks */
- int i; /* loop counter */
- int seqidx; /* counter for seqs */
- int alen; /* length of alignment */
- int warn_names; /* becomes TRUE if names don't match between blocks */
- int headnum; /* seqidx in per-sequence header info */
- int currlen;
- int count;
- int have_cs = 0;
- int have_rf = 0;
- AINFO base_ainfo, *ainfo; /* hack: used to be passed ptr to AINFO */
-
-
- /* Convert from MSA interface to what old ReadSELEX() did:
- * - copy our open fp, rather than opening file
- * - verify that we're not reading a gzip or stdin
- */
- if (feof(afp->f)) return NULL;
- if (afp->do_gzip || afp->do_stdin)
- Die("Can't read a SELEX format alignment from a pipe, stdin, or gzip'ed file");
- fp = afp->f;
- ainfo = &base_ainfo;
-
- /***************************************************
- * First pass across file.
- * Count seqs, get names, determine column info
- * Determine what sorts of info are active in this file.
- ***************************************************/
-
- InitAinfo(ainfo);
- /* get first line of the block
- * (non-comment, non-blank) */
- do
- {
- if (fgets(buffer, LINEBUFLEN, fp) == NULL)
- { squid_errno = SQERR_NODATA; return 0; }
- strcpy(bufcpy, buffer);
- if (*buffer == '#')
- {
- if (strncmp(buffer, "#=CS", 4) == 0) have_cs = 1;
- else if (strncmp(buffer, "#=RF", 4) == 0) have_rf = 1;
- }
- }
- while ((nptr = strtok(bufcpy, WHITESPACE)) == NULL ||
- (strchr(commentsyms, *nptr) != NULL));
-
- blocknum = 0;
- warn_names = FALSE;
- while (!feof(fp))
- {
- /* allocate for info about this block. */
- if (blocknum == 0)
- blocks = (struct block_struc *) MallocOrDie (sizeof(struct block_struc));
- else
- blocks = (struct block_struc *) ReallocOrDie (blocks, (blocknum+1) * sizeof(struct block_struc));
- blocks[blocknum].lcol = LINEBUFLEN+1;
- blocks[blocknum].rcol = -1;
-
- currnum = 0;
- while (nptr != NULL) /* becomes NULL when this block ends. */
- {
- /* First block only: save names */
- if (blocknum == 0)
- {
- if (currnum == 0)
- ainfo->sqinfo = (SQINFO *) MallocOrDie (sizeof(SQINFO));
- else
- ainfo->sqinfo = (SQINFO *) ReallocOrDie (ainfo->sqinfo, (currnum + 1) * sizeof(SQINFO));
-
- ainfo->sqinfo[currnum].flags = 0;
- SetSeqinfoString(&(ainfo->sqinfo[currnum]), nptr, SQINFO_NAME);
- }
- else /* in each additional block: check names */
- {
- if (strcmp(ainfo->sqinfo[currnum].name, nptr) != 0)
- warn_names = TRUE;
- }
- currnum++;
-
- /* check rcol, lcol */
- if ((sptr = strtok(NULL, WHITESPACE)) != NULL)
- {
- /* is this the furthest left we've
- seen word 2 in this block? */
- if (sptr - bufcpy < blocks[blocknum].lcol)
- blocks[blocknum].lcol = sptr - bufcpy;
- /* look for right side in buffer */
- for (sptr = buffer + strlen(buffer) - 1;
- strchr(WHITESPACE, *sptr) != NULL;
- sptr --)
- /* do nothing */ ;
- if (sptr - buffer > blocks[blocknum].rcol)
- blocks[blocknum].rcol = sptr - buffer;
- }
-
- /* get the next line; blank line means end of block */
- do
- {
- if (fgets(buffer, LINEBUFLEN, fp) == NULL)
- { nptr = NULL; break; }
- strcpy(bufcpy, buffer);
-
- if (strncmp(buffer, "#=SS", 4) == 0) ainfo->sqinfo[currnum-1].flags |= SQINFO_SS;
- else if (strncmp(buffer, "#=SA", 4) == 0) ainfo->sqinfo[currnum-1].flags |= SQINFO_SA;
- else if (strncmp(buffer, "#=CS", 4) == 0) have_cs = 1;
- else if (strncmp(buffer, "#=RF", 4) == 0) have_rf = 1;
-
- if ((nptr = strtok(bufcpy, WHITESPACE)) == NULL)
- break;
- } while (strchr(commentsyms, *nptr) != NULL);
- }
-
-
- /* check that number of sequences matches expected */
- if (blocknum == 0)
- num = currnum;
- else if (currnum != num)
- Die("Parse error in ReadSELEX()");
- blocknum++;
-
- /* get first line of next block
- * (non-comment, non-blank) */
- do
- {
- if (fgets(buffer, LINEBUFLEN, fp) == NULL) { nptr = NULL; break; }
- strcpy(bufcpy, buffer);
- }
- while ((nptr = strtok(bufcpy, WHITESPACE)) == NULL ||
- (strchr(commentsyms, *nptr) != NULL));
- }
-
-
- /***************************************************
- * Get ready for second pass:
- * figure out the length of the alignment
- * malloc space
- * rewind the file
- ***************************************************/
-
- alen = 0;
- for (currblock = 0; currblock < blocknum; currblock++)
- alen += blocks[currblock].rcol - blocks[currblock].lcol + 1;
-
- rewind(fp);
-
- /* allocations. we can't use AllocateAlignment because of
- * the way we already used ainfo->sqinfo.
- */
- aseqs = (char **) MallocOrDie (num * sizeof(char *));
- if (have_cs)
- ainfo->cs = (char *) MallocOrDie ((alen+1) * sizeof(char));
- if (have_rf)
- ainfo->rf = (char *) MallocOrDie ((alen+1) * sizeof(char));
-
-
-
- for (i = 0; i < num; i++)
- {
- aseqs[i] = (char *) MallocOrDie ((alen+1) * sizeof(char));
- if (ainfo->sqinfo[i].flags & SQINFO_SS)
- ainfo->sqinfo[i].ss = (char *) MallocOrDie ((alen+1) * sizeof(char));
- if (ainfo->sqinfo[i].flags & SQINFO_SA)
- ainfo->sqinfo[i].sa = (char *) MallocOrDie ((alen+1) * sizeof(char));
- }
-
- ainfo->alen = alen;
- ainfo->nseq = num;
- ainfo->wgt = (float *) MallocOrDie (sizeof(float) * num);
- FSet(ainfo->wgt, num, 1.0);
-
- /***************************************************
- * Second pass across file. Parse header; assemble sequences
- ***************************************************/
- /* We've now made a complete first pass over the file. We know how
- * many blocks it contains, we know the number of seqs in the first
- * block, and we know every block has the same number of blocks;
- * so we can be a bit more cavalier about error-checking as we
- * make the second pass.
- */
-
- /* Look for header
- */
- headnum = 0;
- for (;;)
- {
- if (fgets(buffer, LINEBUFLEN, fp) == NULL)
- Die("Parse error in ReadSELEX()");
- strcpy(bufcpy, buffer);
- if ((nptr = strtok(bufcpy, WHITESPACE)) == NULL) continue; /* skip blank lines */
-
- if (strcmp(nptr, "#=AU") == 0 && (sptr = strtok(NULL, "\n")) != NULL)
- ainfo->au = Strdup(sptr);
- else if (strcmp(nptr, "#=ID") == 0 && (sptr = strtok(NULL, "\n")) != NULL)
- ainfo->name = Strdup(sptr);
- else if (strcmp(nptr, "#=AC") == 0 && (sptr = strtok(NULL, "\n")) != NULL)
- ainfo->acc = Strdup(sptr);
- else if (strcmp(nptr, "#=DE") == 0 && (sptr = strtok(NULL, "\n")) != NULL)
- ainfo->desc = Strdup(sptr);
- else if (strcmp(nptr, "#=GA") == 0)
- {
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=GA line in ReadSELEX()");
- ainfo->ga1 = atof(sptr);
-
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=GA line in ReadSELEX()");
- ainfo->ga2 = atof(sptr);
-
- ainfo->flags |= AINFO_GA;
- }
- else if (strcmp(nptr, "#=TC") == 0)
- {
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=TC line in ReadSELEX()");
- ainfo->tc1 = atof(sptr);
-
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=TC line in ReadSELEX()");
- ainfo->tc2 = atof(sptr);
-
- ainfo->flags |= AINFO_TC;
- }
- else if (strcmp(nptr, "#=NC") == 0)
- {
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=NC line in ReadSELEX()");
- ainfo->nc1 = atof(sptr);
-
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=NC line in ReadSELEX()");
- ainfo->nc2 = atof(sptr);
-
- ainfo->flags |= AINFO_NC;
- }
- else if (strcmp(nptr, "#=SQ") == 0) /* per-sequence header info */
- {
- /* first field is the name */
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX()");
- if (strcmp(sptr, ainfo->sqinfo[headnum].name) != 0) warn_names = TRUE;
-
- /* second field is the weight */
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX()");
- if (!IsReal(sptr))
- Die("Parse error in #=SQ line in ReadSELEX(): weight is not a number");
- ainfo->wgt[headnum] = atof(sptr);
-
- /* third field is database source id */
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX(): incomplete line");
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_ID);
-
- /* fourth field is database accession number */
- if ((sptr = strtok(NULL, WHITESPACE)) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX(): incomplete line");
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_ACC);
-
- /* fifth field is start..stop::olen */
- if ((sptr = strtok(NULL, ".:")) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX(): incomplete line");
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_START);
-
- if ((sptr = strtok(NULL, ".:")) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX(): incomplete line");
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_STOP);
-
- if ((sptr = strtok(NULL, ":\t ")) == NULL)
- Die("Parse error in #=SQ line in ReadSELEX(): incomplete line");
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_OLEN);
-
- /* rest of line is optional description */
- if ((sptr = strtok(NULL, "\n")) != NULL)
- SetSeqinfoString(&(ainfo->sqinfo[headnum]), sptr, SQINFO_DESC);
-
- headnum++;
- }
- else if (strcmp(nptr, "#=CS") == 0) break;
- else if (strcmp(nptr, "#=RF") == 0) break;
- else if (strchr(commentsyms, *nptr) == NULL) break; /* non-comment, non-header */
- }
-
-
- currlen = 0;
- for (currblock = 0 ; currblock < blocknum; currblock++)
- {
- /* parse the block */
- seqidx = 0;
- while (nptr != NULL)
- {
- /* Consensus structure */
- if (strcmp(nptr, "#=CS") == 0)
- {
- if (! copy_alignment_line(ainfo->cs, currlen, strlen(nptr)-1,
- buffer, blocks[currblock].lcol, blocks[currblock].rcol, (char) '.'))
- Die("Parse error in #=CS line in ReadSELEX()");
- }
-
- /* Reference coordinates */
- else if (strcmp(nptr, "#=RF") == 0)
- {
- if (! copy_alignment_line(ainfo->rf, currlen, strlen(nptr)-1,
- buffer, blocks[currblock].lcol, blocks[currblock].rcol, (char) '.'))
- Die("Parse error in #=RF line in ReadSELEX()");
- }
- /* Individual secondary structure */
- else if (strcmp(nptr, "#=SS") == 0)
- {
- if (! copy_alignment_line(ainfo->sqinfo[seqidx-1].ss, currlen, strlen(nptr)-1,
- buffer, blocks[currblock].lcol,
- blocks[currblock].rcol, (char) '.'))
- Die("Parse error in #=SS line in ReadSELEX()");
- }
-
- /* Side chain % surface accessibility code */
- else if (strcmp(nptr, "#=SA") == 0)
- {
- if (! copy_alignment_line(ainfo->sqinfo[seqidx-1].sa, currlen, strlen(nptr)-1,
- buffer, blocks[currblock].lcol,
- blocks[currblock].rcol, (char) '.'))
- Die("Parse error in #=SA line in ReadSELEX()");
- }
- /* Aligned sequence; avoid unparsed machine comments */
- else if (strncmp(nptr, "#=", 2) != 0)
- {
- if (! copy_alignment_line(aseqs[seqidx], currlen, strlen(nptr)-1,
- buffer, blocks[currblock].lcol, blocks[currblock].rcol, (char) '.'))
- Die("Parse error in alignment line in ReadSELEX()");
- seqidx++;
- }
-
- /* get next line */
- for (;;)
- {
- nptr = NULL;
- if (fgets(buffer, LINEBUFLEN, fp) == NULL) break; /* EOF */
- strcpy(bufcpy, buffer);
- if ((nptr = strtok(bufcpy, WHITESPACE)) == NULL) break; /* blank */
- if (strncmp(buffer, "#=", 2) == 0) break; /* machine comment */
- if (strchr(commentsyms, *nptr) == NULL) break; /* data */
- }
- } /* end of a block */
-
- currlen += blocks[currblock].rcol - blocks[currblock].lcol + 1;
-
- /* get line 1 of next block */
- for (;;)
- {
- if (fgets(buffer, LINEBUFLEN, fp) == NULL) break; /* no data */
- strcpy(bufcpy, buffer);
- if ((nptr = strtok(bufcpy, WHITESPACE)) == NULL) continue; /* blank */
- if (strncmp(buffer, "#=", 2) == 0) break; /* machine comment */
- if (strchr(commentsyms, *nptr) == NULL) break; /* non-comment */
- }
- } /* end of the file */
-
- /* Lengths in sqinfo are for raw sequence (ungapped),
- * and SS, SA are 0..rlen-1 not 0..alen-1.
- * Only the seqs with structures come out of here with lengths set.
- */
- for (seqidx = 0; seqidx < num; seqidx++)
- {
- int apos, rpos;
- /* secondary structures */
- if (ainfo->sqinfo[seqidx].flags & SQINFO_SS)
- {
- for (apos = rpos = 0; apos < alen; apos++)
- if (! isgap(aseqs[seqidx][apos]))
- {
- ainfo->sqinfo[seqidx].ss[rpos] = ainfo->sqinfo[seqidx].ss[apos];
- rpos++;
- }
- ainfo->sqinfo[seqidx].ss[rpos] = '\0';
- }
- /* Surface accessibility */
- if (ainfo->sqinfo[seqidx].flags & SQINFO_SA)
- {
- for (apos = rpos = 0; apos < alen; apos++)
- if (! isgap(aseqs[seqidx][apos]))
- {
- ainfo->sqinfo[seqidx].sa[rpos] = ainfo->sqinfo[seqidx].sa[apos];
- rpos++;
- }
- ainfo->sqinfo[seqidx].sa[rpos] = '\0';
- }
- }
-
- /* NULL-terminate all the strings */
- if (ainfo->rf != NULL) ainfo->rf[alen] = '\0';
- if (ainfo->cs != NULL) ainfo->cs[alen] = '\0';
- for (seqidx = 0; seqidx < num; seqidx++)
- aseqs[seqidx][alen] = '\0';
-
- /* find raw sequence lengths for sqinfo */
- for (seqidx = 0; seqidx < num; seqidx++)
- {
- count = 0;
- for (sptr = aseqs[seqidx]; *sptr != '\0'; sptr++)
- if (!isgap(*sptr)) count++;
- ainfo->sqinfo[seqidx].len = count;
- ainfo->sqinfo[seqidx].flags |= SQINFO_LEN;
- }
-
-
- /***************************************************
- * Garbage collection and return
- ***************************************************/
- free(blocks);
- if (warn_names)
-#ifdef CLUSTALO
- Warning("sequences may be in different orders in blocks of %s?", afp->fname);
-#else
- Warn("sequences may be in different orders in blocks of %s?", afp->fname);
-#endif
-
- /* Convert back to MSA structure. (Wasteful kludge.)
- */
- msa = MSAFromAINFO(aseqs, ainfo);
- MSAVerifyParse(msa);
- FreeAlignment(aseqs, ainfo);
- return msa;
-}
-
-
-/* Function: WriteSELEX()
- * Date: SRE, Mon Jun 14 13:13:14 1999 [St. Louis]
- *
- * Purpose: Write a SELEX file in multiblock format.
- *
- * Args: fp - file that's open for writing
- * msa - multiple sequence alignment object
- *
- * Returns: (void)
- */
-void
-WriteSELEX(FILE *fp, MSA *msa)
-{
- actually_write_selex(fp, msa, 50); /* 50 char per block */
-}
-
-/* Function: WriteSELEXOneBlock()
- * Date: SRE, Mon Jun 14 13:14:56 1999 [St. Louis]
- *
- * Purpose: Write a SELEX alignment file in Pfam's single-block
- * format style. A wrapper for actually_write_selex().
- *
- * Args: fp - file that's open for writing
- * msa- alignment to write
- *
- * Returns: (void)
- */
-void
-WriteSELEXOneBlock(FILE *fp, MSA *msa)
-{
- actually_write_selex(fp, msa, msa->alen); /* one big block */
-}
-
-
-/* Function: actually_write_selex()
- * Date: SRE, Mon Jun 14 12:54:46 1999 [St. Louis]
- *
- * Purpose: Write an alignment in SELEX format to an open
- * file. This is the function that actually does
- * the work. The API's WriteSELEX() and
- * WriteSELEXOneBlock() are wrappers.
- *
- * Args: fp - file that's open for writing
- * msa - alignment to write
- * cpl - characters to write per line in alignment block
- *
- * Returns: (void)
- */
-static void
-actually_write_selex(FILE *fp, MSA *msa, int cpl)
-{
- int i;
- int len = 0;
- int namewidth;
- char *buf;
- int currpos;
-
- buf = malloc(sizeof(char) * (cpl+101)); /* 100 chars allowed for name, etc. */
-
- /* Figure out how much space we need for name + markup
- * to keep the alignment in register, for easier human viewing --
- * even though Stockholm format doesn't care about visual
- * alignment.
- */
- namewidth = 0;
- for (i = 0; i < msa->nseq; i++)
- if ((len = strlen(msa->sqname[i])) > namewidth)
- namewidth = len;
- if (namewidth < 6) namewidth = 6; /* minimum space for markup tags */
-
- /* Free text comments
- */
- for (i = 0; i < msa->ncomment; i++)
- fprintf(fp, "# %s\n", msa->comment[i]);
- if (msa->ncomment > 0) fprintf(fp, "\n");
-
- /* Per-file annotation
- */
- if (msa->name != NULL) fprintf(fp, "#=ID %s\n", msa->name);
- if (msa->acc != NULL) fprintf(fp, "#=AC %s\n", msa->acc);
- if (msa->desc != NULL) fprintf(fp, "#=DE %s\n", msa->desc);
- if (msa->au != NULL) fprintf(fp, "#=AU %s\n", msa->au);
-
- /* Thresholds are hacky. Pfam has two. Rfam has one.
- */
- if (msa->cutoff_is_set[MSA_CUTOFF_GA1] && msa->cutoff_is_set[MSA_CUTOFF_GA2])
- fprintf(fp, "#=GA %.1f %.1f\n", msa->cutoff[MSA_CUTOFF_GA1], msa->cutoff[MSA_CUTOFF_GA2]);
- else if (msa->cutoff_is_set[MSA_CUTOFF_GA1])
- fprintf(fp, "#=GA %.1f\n", msa->cutoff[MSA_CUTOFF_GA1]);
- if (msa->cutoff_is_set[MSA_CUTOFF_NC1] && msa->cutoff_is_set[MSA_CUTOFF_NC2])
- fprintf(fp, "#=NC %.1f %.1f\n", msa->cutoff[MSA_CUTOFF_NC1], msa->cutoff[MSA_CUTOFF_NC2]);
- else if (msa->cutoff_is_set[MSA_CUTOFF_NC1])
- fprintf(fp, "#=NC %.1f\n", msa->cutoff[MSA_CUTOFF_NC1]);
- if (msa->cutoff_is_set[MSA_CUTOFF_TC1] && msa->cutoff_is_set[MSA_CUTOFF_TC2])
- fprintf(fp, "#=TC %.1f %.1f\n", msa->cutoff[MSA_CUTOFF_TC1], msa->cutoff[MSA_CUTOFF_TC2]);
- else if (msa->cutoff_is_set[MSA_CUTOFF_TC1])
- fprintf(fp, "#=TC %.1f\n", msa->cutoff[MSA_CUTOFF_TC1]);
-
- /* Per-sequence annotation
- */
- for (i = 0; i < msa->nseq; i++)
- fprintf(fp, "#=SQ %-*.*s %6.4f %s %s %d..%d::%d %s\n",
- namewidth, namewidth, msa->sqname[i],
- msa->wgt[i],
- "-", /* MSA has no ID field */
- (msa->sqacc != NULL && msa->sqacc[i] != NULL) ? msa->sqacc[i] : "-",
- 0, 0, 0, /* MSA has no start, stop, olen field */
- (msa->sqdesc != NULL && msa->sqdesc[i] != NULL) ? msa->sqdesc[i] : "-");
- fprintf(fp, "\n");
-
- /* Alignment section:
- */
- for (currpos = 0; currpos < msa->alen; currpos += cpl)
- {
- if (currpos > 0) fprintf(fp, "\n");
-
- if (msa->ss_cons != NULL) {
- strncpy(buf, msa->ss_cons + currpos, cpl);
- buf[cpl] = '\0';
- fprintf(fp, "%-*.*s %s\n", namewidth, namewidth, "#=CS", buf);
- }
- if (msa->rf != NULL) {
- strncpy(buf, msa->rf + currpos, cpl);
- buf[cpl] = '\0';
- fprintf(fp, "%-*.*s %s\n", namewidth, namewidth, "#=RF", buf);
- }
- for (i = 0; i < msa->nseq; i++)
- {
- strncpy(buf, msa->aseq[i] + currpos, cpl);
- buf[cpl] = '\0';
- fprintf(fp, "%-*.*s %s\n", namewidth, namewidth, msa->sqname[i], buf);
-
- if (msa->ss != NULL && msa->ss[i] != NULL) {
- strncpy(buf, msa->ss[i] + currpos, cpl);
- buf[cpl] = '\0';
- fprintf(fp, "%-*.*s %s\n", namewidth, namewidth, "#=SS", buf);
- }
- if (msa->sa != NULL && msa->sa[i] != NULL) {
- strncpy(buf, msa->sa[i] + currpos, cpl);
- buf[cpl] = '\0';
- fprintf(fp, "%-*.*s %s\n", namewidth, namewidth, "#=SA", buf);
- }
- }
- }
- free(buf);
-}
-
-
-/* Function: copy_alignment_line()
- *
- * Purpose: Given a line from an alignment file, and bounds lcol,rcol
- * on what part of it may be sequence, save the alignment into
- * aseq starting at position apos.
- *
- * name_rcol is set to the rightmost column this aseqs's name
- * occupies; if name_rcol >= lcol, we have a special case in
- * which the name intrudes into the sequence zone.
- */
-static int
-copy_alignment_line(char *aseq, int apos, int name_rcol,
- char *buffer, int lcol, int rcol, char gapsym)
-{
- char *s1, *s2;
- int i;
-
- s1 = aseq + apos;
- s2 = buffer; /* be careful that buffer doesn't end before lcol! */
- for (i = 0; i < lcol; i++)
- if (*s2) s2++;
-
- for (i = lcol; i <= rcol; i++)
- {
- if (*s2 == '\t') {
-#ifdef CLUSTALO
- Warning("TAB characters will corrupt a SELEX alignment! Please remove them first.");
-#else
- Warn("TAB characters will corrupt a SELEX alignment! Please remove them first.");
-#endif
- return 0;
- }
- if (name_rcol >= i) /* name intrusion special case: pad left w/ gaps */
- *s1 = gapsym;
- /* short buffer special case: pad right w/ gaps */
- else if (*s2 == '\0' || *s2 == '\n')
- *s1 = gapsym;
-
- else if (*s2 == ' ') /* new: disallow spaces as gap symbols */
- *s1 = gapsym;
-
- else /* normal case: copy buffer into aseq */
- *s1 = *s2;
-
- s1++;
- if (*s2) s2++;
- }
- return 1;
-}
-
-
-
-
-
-/* Function: DealignAseqs()
- *
- * Given an array of (num) aligned sequences aseqs,
- * strip the gaps. Store the raw sequences in a new allocated array.
- *
- * Caller is responsible for free'ing the memory allocated to
- * rseqs.
- *
- * Returns 1 on success. Returns 0 and sets squid_errno on
- * failure.
- */
-int
-DealignAseqs(char **aseqs, int num, char ***ret_rseqs)
-{
- char **rseqs; /* de-aligned sequence array */
- int idx; /* counter for sequences */
- int depos; /* position counter for dealigned seq*/
- int apos; /* position counter for aligned seq */
- int seqlen; /* length of aligned seq */
-
- /* alloc space */
- rseqs = (char **) MallocOrDie (num * sizeof(char *));
- /* main loop */
- for (idx = 0; idx < num; idx++)
- {
- seqlen = strlen(aseqs[idx]);
- /* alloc space */
- rseqs[idx] = (char *) MallocOrDie ((seqlen + 1) * sizeof(char));
-
- /* strip gaps */
- depos = 0;
- for (apos = 0; aseqs[idx][apos] != '\0'; apos++)
- if (!isgap(aseqs[idx][apos]))
- {
- rseqs[idx][depos] = aseqs[idx][apos];
- depos++;
- }
- rseqs[idx][depos] = '\0';
- }
- *ret_rseqs = rseqs;
- return 1;
-}
-
-
-/* Function: IsSELEXFormat()
- *
- * Return TRUE if filename may be in SELEX format.
- *
- * Accuracy is sacrificed for speed; a TRUE return does
- * *not* guarantee that the file will pass the stricter
- * error-checking of ReadSELEX(). All it checks is that
- * the first 500 non-comment lines of a file are
- * blank, or if there's a second "word" on the line
- * it looks like sequence (i.e., it's not kOtherSeq).
- *
- * Returns TRUE or FALSE.
- */
-int
-IsSELEXFormat(char *filename)
-{
- FILE *fp; /* ptr to open sequence file */
- char buffer[LINEBUFLEN];
- char *sptr; /* ptr to first word */
- int linenum;
-
-
- if ((fp = fopen(filename, "r")) == NULL)
- { squid_errno = SQERR_NOFILE; return 0; }
-
- linenum = 0;
- while (linenum < 500 &&
- fgets(buffer, LINEBUFLEN, fp) != NULL)
- {
- linenum++;
- /* dead giveaways for extended SELEX */
- if (strncmp(buffer, "#=AU", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=ID", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=AC", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=DE", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=GA", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=TC", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=NC", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=SQ", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=SS", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=CS", 4) == 0) goto DONE;
- else if (strncmp(buffer, "#=RF", 4) == 0) goto DONE;
-
- /* a comment? */
- if (strchr(commentsyms, *buffer) != NULL) continue;
-
- /* a blank line? */
- if ((sptr = strtok(buffer, WHITESPACE)) == NULL) continue;
-
- /* a one-word line (name only)
- is possible, though rare */
- if ((sptr = strtok(NULL, "\n")) == NULL) continue;
-
- if (Seqtype(sptr) == kOtherSeq) {fclose(fp); return 0;}
- }
-
- DONE:
- fclose(fp);
- return 1;
-}
-
-
-
-
-
-
-
-
git://source.jalview.org / jabaws.git / blobdiff