--- /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.
+ *****************************************************************/
+
+/* msa.c
+ * SRE, Mon May 17 10:48:47 1999
+ *
+ * SQUID's interface for multiple sequence alignment
+ * manipulation: access to the MSA object.
+ *
+ * RCS $Id: msa.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: msa.c,v 1.18 2002/10/12 04:40:35 eddy Exp)
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "squid.h"
+#include "msa.h" /* multiple sequence alignment object support */
+#include "gki.h" /* string indexing hashtable code */
+#include "ssi.h" /* SSI sequence file indexing code */
+
+/* Function: MSAAlloc()
+ * Date: SRE, Tue May 18 10:45:47 1999 [St. Louis]
+ *
+ * Purpose: Allocate an MSA structure, return a pointer
+ * to it.
+ *
+ * Designed to be used in three ways:
+ * 1) We know exactly the dimensions of the alignment:
+ * both nseq and alen.
+ * msa = MSAAlloc(nseq, alen);
+ *
+ * 2) We know the number of sequences but not alen.
+ * (We add sequences later.)
+ * msa = MSAAlloc(nseq, 0);
+ *
+ * 3) We even don't know the number of sequences, so
+ * we'll have to dynamically expand allocations.
+ * We provide a blocksize for the allocation expansion,
+ * and expand when needed.
+ * msa = MSAAlloc(10, 0);
+ * if (msa->nseq == msa->nseqalloc) MSAExpand(msa);
+ *
+ * Args: nseq - number of sequences, or nseq allocation blocksize
+ * alen - length of alignment in columns, or 0
+ *
+ * Returns: pointer to new MSA object, w/ all values initialized.
+ * Note that msa->nseq is initialized to 0, though space
+ * is allocated.
+ *
+ * Diagnostics: "always works". Die()'s on memory allocation failure.
+ *
+ */
+MSA *
+MSAAlloc(int nseq, int alen)
+{
+ MSA *msa;
+ int i;
+
+ msa = MallocOrDie(sizeof(MSA));
+ msa->aseq = MallocOrDie(sizeof(char *) * nseq);
+ msa->sqname = MallocOrDie(sizeof(char *) * nseq);
+ msa->sqlen = MallocOrDie(sizeof(int) * nseq);
+ msa->wgt = MallocOrDie(sizeof(float) * nseq);
+
+ for (i = 0; i < nseq; i++)
+ {
+ msa->sqname[i] = NULL;
+ msa->sqlen[i] = 0;
+ msa->wgt[i] = -1.0;
+
+ if (alen != 0) msa->aseq[i] = MallocOrDie(sizeof(char) * (alen+1));
+ else msa->aseq[i] = NULL;
+ }
+
+ msa->alen = alen;
+ msa->nseq = 0;
+ msa->nseqalloc = nseq;
+ msa->nseqlump = nseq;
+
+ msa->flags = 0;
+ msa->type = kOtherSeq;
+ msa->name = NULL;
+ msa->desc = NULL;
+ msa->acc = NULL;
+ msa->au = NULL;
+ msa->ss_cons = NULL;
+ msa->sa_cons = NULL;
+ msa->rf = NULL;
+ msa->sqacc = NULL;
+ msa->sqdesc = NULL;
+ msa->ss = NULL;
+ msa->sslen = NULL;
+ msa->sa = NULL;
+ msa->salen = NULL;
+ msa->index = GKIInit();
+ msa->lastidx = 0;
+
+ for (i = 0; i < MSA_MAXCUTOFFS; i++) {
+ msa->cutoff[i] = 0.;
+ msa->cutoff_is_set[i] = FALSE;
+ }
+
+ /* Initialize unparsed optional markup
+ */
+ msa->comment = NULL;
+ msa->ncomment = 0;
+ msa->alloc_ncomment = 0;
+
+ msa->gf_tag = NULL;
+ msa->gf = NULL;
+ msa->ngf = 0;
+
+ msa->gs_tag = NULL;
+ msa->gs = NULL;
+ msa->gs_idx = NULL;
+ msa->ngs = 0;
+
+ msa->gc_tag = NULL;
+ msa->gc = NULL;
+ msa->gc_idx = NULL;
+ msa->ngc = 0;
+
+ msa->gr_tag = NULL;
+ msa->gr = NULL;
+ msa->gr_idx = NULL;
+ msa->ngr = 0;
+
+ /* Done. Return the alloced, initialized structure
+ */
+ return msa;
+}
+
+/* Function: MSAExpand()
+ * Date: SRE, Tue May 18 11:06:53 1999 [St. Louis]
+ *
+ * Purpose: Increase the sequence allocation in an MSA
+ * by msa->nseqlump. (Typically used when we're reading
+ * in an alignment sequentially from a file,
+ * so we don't know nseq until we're done.)
+ *
+ * Args: msa - the MSA object
+ *
+ * Returns: (void)
+ *
+ */
+void
+MSAExpand(MSA *msa)
+{
+ int i,j;
+
+ msa->nseqalloc += msa->nseqlump;
+
+ msa->aseq = ReallocOrDie(msa->aseq, sizeof(char *) * msa->nseqalloc);
+ msa->sqname = ReallocOrDie(msa->sqname, sizeof(char *) * msa->nseqalloc);
+ msa->sqlen = ReallocOrDie(msa->sqlen, sizeof(char *) * msa->nseqalloc);
+ msa->wgt = ReallocOrDie(msa->wgt, sizeof(float) * msa->nseqalloc);
+
+ if (msa->ss != NULL) {
+ msa->ss = ReallocOrDie(msa->ss, sizeof(char *) * msa->nseqalloc);
+ msa->sslen = ReallocOrDie(msa->sslen, sizeof(int) * msa->nseqalloc);
+ }
+ if (msa->sa != NULL) {
+ msa->sa = ReallocOrDie(msa->sa, sizeof(char *) * msa->nseqalloc);
+ msa->salen = ReallocOrDie(msa->salen, sizeof(int) * msa->nseqalloc);
+ }
+ if (msa->sqacc != NULL)
+ msa->sqacc = ReallocOrDie(msa->sqacc, sizeof(char *) * msa->nseqalloc);
+ if (msa->sqdesc != NULL)
+ msa->sqdesc =ReallocOrDie(msa->sqdesc,sizeof(char *) * msa->nseqalloc);
+
+ for (i = msa->nseqalloc-msa->nseqlump; i < msa->nseqalloc; i++)
+ {
+ msa->sqname[i] = NULL;
+ msa->wgt[i] = -1.0;
+
+ if (msa->sqacc != NULL) msa->sqacc[i] = NULL;
+ if (msa->sqdesc != NULL) msa->sqdesc[i] = NULL;
+
+ if (msa->alen != 0)
+ msa->aseq[i] = ReallocOrDie(msa->aseq[i], sizeof(char) * (msa->alen+1));
+ else msa->aseq[i] = NULL;
+ msa->sqlen[i] = 0;
+
+ if (msa->ss != NULL) {
+ if (msa->alen != 0)
+ msa->ss[i] = ReallocOrDie(msa->ss[i], sizeof(char) * (msa->alen+1));
+ else msa->ss[i] = NULL;
+ msa->sslen[i] = 0;
+ }
+ if (msa->sa != NULL) {
+ if (msa->alen != 0)
+ msa->sa[i] = ReallocOrDie(msa->ss[i], sizeof(char) * (msa->alen+1));
+ else
+ msa->sa[i] = NULL;
+ msa->salen[i] = 0;
+ }
+ }
+
+ /* Reallocate and re-init for unparsed #=GS tags, if we have some.
+ * gs is [0..ngs-1][0..nseq-1][], so we're reallocing the middle
+ * set of pointers.
+ */
+ if (msa->gs != NULL)
+ for (i = 0; i < msa->ngs; i++)
+ {
+ if (msa->gs[i] != NULL)
+ {
+ msa->gs[i] = ReallocOrDie(msa->gs[i], sizeof(char *) * msa->nseqalloc);
+ for (j = msa->nseqalloc-msa->nseqlump; j < msa->nseqalloc; j++)
+ msa->gs[i][j] = NULL;
+ }
+ }
+
+ /* Reallocate and re-init for unparsed #=GR tags, if we have some.
+ * gr is [0..ngs-1][0..nseq-1][], so we're reallocing the middle
+ * set of pointers.
+ */
+ if (msa->gr != NULL)
+ for (i = 0; i < msa->ngr; i++)
+ {
+ if (msa->gr[i] != NULL)
+ {
+ msa->gr[i] = ReallocOrDie(msa->gr[i], sizeof(char *) * msa->nseqalloc);
+ for (j = msa->nseqalloc-msa->nseqlump; j < msa->nseqalloc; j++)
+ msa->gr[i][j] = NULL;
+ }
+ }
+
+ return;
+}
+
+/* Function: MSAFree()
+ * Date: SRE, Tue May 18 11:20:16 1999 [St. Louis]
+ *
+ * Purpose: Free a multiple sequence alignment structure.
+ *
+ * Args: msa - the alignment
+ *
+ * Returns: (void)
+ */
+void
+MSAFree(MSA *msa)
+{
+ Free2DArray((void **) msa->aseq, msa->nseq);
+ Free2DArray((void **) msa->sqname, msa->nseq);
+ Free2DArray((void **) msa->sqacc, msa->nseq);
+ Free2DArray((void **) msa->sqdesc, msa->nseq);
+ Free2DArray((void **) msa->ss, msa->nseq);
+ Free2DArray((void **) msa->sa, msa->nseq);
+
+ if (msa->sqlen != NULL) free(msa->sqlen);
+ if (msa->wgt != NULL) free(msa->wgt);
+
+ if (msa->name != NULL) free(msa->name);
+ if (msa->desc != NULL) free(msa->desc);
+ if (msa->acc != NULL) free(msa->acc);
+ if (msa->au != NULL) free(msa->au);
+ if (msa->ss_cons != NULL) free(msa->ss_cons);
+ if (msa->sa_cons != NULL) free(msa->sa_cons);
+ if (msa->rf != NULL) free(msa->rf);
+ if (msa->sslen != NULL) free(msa->sslen);
+ if (msa->salen != NULL) free(msa->salen);
+
+ Free2DArray((void **) msa->comment, msa->ncomment);
+ Free2DArray((void **) msa->gf_tag, msa->ngf);
+ Free2DArray((void **) msa->gf, msa->ngf);
+ Free2DArray((void **) msa->gs_tag, msa->ngs);
+ Free3DArray((void ***)msa->gs, msa->ngs, msa->nseq);
+ Free2DArray((void **) msa->gc_tag, msa->ngc);
+ Free2DArray((void **) msa->gc, msa->ngc);
+ Free2DArray((void **) msa->gr_tag, msa->ngr);
+ Free3DArray((void ***)msa->gr, msa->ngr, msa->nseq);
+
+ GKIFree(msa->index);
+ GKIFree(msa->gs_idx);
+ GKIFree(msa->gc_idx);
+ GKIFree(msa->gr_idx);
+
+ free(msa);
+}
+
+
+/* Function: MSASetSeqAccession()
+ * Date: SRE, Mon Jun 21 04:13:33 1999 [Sanger Centre]
+ *
+ * Purpose: Set a sequence accession in an MSA structure.
+ * Handles some necessary allocation/initialization.
+ *
+ * Args: msa - multiple alignment to add accession to
+ * seqidx - index of sequence to attach accession to
+ * acc - accession
+ *
+ * Returns: void
+ */
+void
+MSASetSeqAccession(MSA *msa, int seqidx, char *acc)
+{
+ int x;
+
+ if (msa->sqacc == NULL) {
+ msa->sqacc = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (x = 0; x < msa->nseqalloc; x++)
+ msa->sqacc[x] = NULL;
+ }
+ msa->sqacc[seqidx] = sre_strdup(acc, -1);
+}
+
+/* Function: MSASetSeqDescription()
+ * Date: SRE, Mon Jun 21 04:21:09 1999 [Sanger Centre]
+ *
+ * Purpose: Set a sequence description in an MSA structure.
+ * Handles some necessary allocation/initialization.
+ *
+ * Args: msa - multiple alignment to add accession to
+ * seqidx - index of sequence to attach accession to
+ * desc - description
+ *
+ * Returns: void
+ */
+void
+MSASetSeqDescription(MSA *msa, int seqidx, char *desc)
+{
+ int x;
+
+ if (msa->sqdesc == NULL) {
+ msa->sqdesc = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (x = 0; x < msa->nseqalloc; x++)
+ msa->sqdesc[x] = NULL;
+ }
+ msa->sqdesc[seqidx] = sre_strdup(desc, -1);
+}
+
+
+/* Function: MSAAddComment()
+ * Date: SRE, Tue Jun 1 17:37:21 1999 [St. Louis]
+ *
+ * Purpose: Add an (unparsed) comment line to the MSA structure,
+ * allocating as necessary.
+ *
+ * Args: msa - a multiple alignment
+ * s - comment line to add
+ *
+ * Returns: (void)
+ */
+void
+MSAAddComment(MSA *msa, char *s)
+{
+ /* If this is our first recorded comment, we need to malloc();
+ * and if we've filled available space, we need to realloc().
+ * Note the arbitrary lumpsize of 10 lines per allocation...
+ */
+ if (msa->comment == NULL) {
+ msa->comment = MallocOrDie (sizeof(char *) * 10);
+ msa->alloc_ncomment = 10;
+ }
+ if (msa->ncomment == msa->alloc_ncomment) {
+ msa->alloc_ncomment += 10;
+ msa->comment = ReallocOrDie(msa->comment, sizeof(char *) * msa->alloc_ncomment);
+ }
+
+ msa->comment[msa->ncomment] = sre_strdup(s, -1);
+ msa->ncomment++;
+ return;
+}
+
+/* Function: MSAAddGF()
+ * Date: SRE, Wed Jun 2 06:53:54 1999 [bus to Madison]
+ *
+ * Purpose: Add an unparsed #=GF markup line to the MSA
+ * structure, allocating as necessary.
+ *
+ * Args: msa - a multiple alignment
+ * tag - markup tag (e.g. "AU")
+ * value - free text markup (e.g. "Alex Bateman")
+ *
+ * Returns: (void)
+ */
+void
+MSAAddGF(MSA *msa, char *tag, char *value)
+{
+ /* If this is our first recorded unparsed #=GF line, we need to malloc();
+ * if we've filled availabl space If we already have a hash index, and the GF
+ * Note the arbitrary lumpsize of 10 lines per allocation...
+ */
+ if (msa->gf_tag == NULL) {
+ msa->gf_tag = MallocOrDie (sizeof(char *) * 10);
+ msa->gf = MallocOrDie (sizeof(char *) * 10);
+ msa->alloc_ngf = 10;
+ }
+ if (msa->ngf == msa->alloc_ngf) {
+ msa->alloc_ngf += 10;
+ msa->gf_tag = ReallocOrDie(msa->gf_tag, sizeof(char *) * msa->alloc_ngf);
+ msa->gf = ReallocOrDie(msa->gf, sizeof(char *) * msa->alloc_ngf);
+ }
+
+ msa->gf_tag[msa->ngf] = sre_strdup(tag, -1);
+ msa->gf[msa->ngf] = sre_strdup(value, -1);
+ msa->ngf++;
+
+ return;
+}
+
+
+/* Function: MSAAddGS()
+ * Date: SRE, Wed Jun 2 06:57:03 1999 [St. Louis]
+ *
+ * Purpose: Add an unparsed #=GS markup line to the MSA
+ * structure, allocating as necessary.
+ *
+ * It's possible that we could get more than one
+ * of the same type of GS tag per sequence; for
+ * example, "DR PDB;" structure links in Pfam.
+ * Hack: handle these by appending to the string,
+ * in a \n separated fashion.
+ *
+ * Args: msa - multiple alignment structure
+ * tag - markup tag (e.g. "AC")
+ * sqidx - index of sequence to assoc markup with (0..nseq-1)
+ * value - markup (e.g. "P00666")
+ *
+ * Returns: 0 on success
+ */
+void
+MSAAddGS(MSA *msa, char *tag, int sqidx, char *value)
+{
+ int tagidx;
+ int i;
+
+ /* Is this an unparsed tag name that we recognize?
+ * If not, handle adding it to index, and reallocating
+ * as needed.
+ */
+ if (msa->gs_tag == NULL) /* first tag? init w/ malloc */
+ {
+ msa->gs_idx = GKIInit();
+ tagidx = GKIStoreKey(msa->gs_idx, tag);
+ SQD_DASSERT1((tagidx == 0));
+ msa->gs_tag = MallocOrDie(sizeof(char *));
+ msa->gs = MallocOrDie(sizeof(char **));
+ msa->gs[0] = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (i = 0; i < msa->nseqalloc; i++)
+ msa->gs[0][i] = NULL;
+ }
+ else
+ {
+ /* new tag? */
+ tagidx = GKIKeyIndex(msa->gs_idx, tag);
+ if (tagidx < 0) { /* it's a new tag name; realloc */
+ tagidx = GKIStoreKey(msa->gs_idx, tag);
+ /* since we alloc in blocks of 1,
+ we always realloc upon seeing
+ a new tag. */
+ SQD_DASSERT1((tagidx == msa->ngs));
+ msa->gs_tag = ReallocOrDie(msa->gs_tag, (msa->ngs+1) * sizeof(char *));
+ msa->gs = ReallocOrDie(msa->gs, (msa->ngs+1) * sizeof(char **));
+ msa->gs[msa->ngs] = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (i = 0; i < msa->nseqalloc; i++)
+ msa->gs[msa->ngs][i] = NULL;
+ }
+ }
+
+ if (tagidx == msa->ngs) {
+ msa->gs_tag[tagidx] = sre_strdup(tag, -1);
+ msa->ngs++;
+ }
+
+ if (msa->gs[tagidx][sqidx] == NULL) /* first annotation of this seq with this tag? */
+ msa->gs[tagidx][sqidx] = sre_strdup(value, -1);
+ else {
+ /* >1 annotation of this seq with this tag; append */
+ int len;
+ if ((len = sre_strcat(&(msa->gs[tagidx][sqidx]), -1, "\n", 1)) < 0)
+ Die("failed to sre_strcat()");
+ if (sre_strcat(&(msa->gs[tagidx][sqidx]), len, value, -1) < 0)
+ Die("failed to sre_strcat()");
+ }
+ return;
+}
+
+/* Function: MSAAppendGC()
+ * Date: SRE, Thu Jun 3 06:25:14 1999 [Madison]
+ *
+ * Purpose: Add an unparsed #=GC markup line to the MSA
+ * structure, allocating as necessary.
+ *
+ * When called multiple times for the same tag,
+ * appends value strings together -- used when
+ * parsing multiblock alignment files, for
+ * example.
+ *
+ * Args: msa - multiple alignment structure
+ * tag - markup tag (e.g. "CS")
+ * value - markup, one char per aligned column
+ *
+ * Returns: (void)
+ */
+void
+MSAAppendGC(MSA *msa, char *tag, char *value)
+{
+ int tagidx;
+
+ /* Is this an unparsed tag name that we recognize?
+ * If not, handle adding it to index, and reallocating
+ * as needed.
+ */
+ if (msa->gc_tag == NULL) /* first tag? init w/ malloc */
+ {
+ msa->gc_tag = MallocOrDie(sizeof(char *));
+ msa->gc = MallocOrDie(sizeof(char *));
+ msa->gc_idx = GKIInit();
+ tagidx = GKIStoreKey(msa->gc_idx, tag);
+ SQD_DASSERT1((tagidx == 0));
+ msa->gc[0] = NULL;
+ }
+ else
+ { /* new tag? */
+ tagidx = GKIKeyIndex(msa->gc_idx, tag);
+ if (tagidx < 0) { /* it's a new tag name; realloc */
+ tagidx = GKIStoreKey(msa->gc_idx, tag);
+ /* since we alloc in blocks of 1,
+ we always realloc upon seeing
+ a new tag. */
+ SQD_DASSERT1((tagidx == msa->ngc));
+ msa->gc_tag = ReallocOrDie(msa->gc_tag, (msa->ngc+1) * sizeof(char **));
+ msa->gc = ReallocOrDie(msa->gc, (msa->ngc+1) * sizeof(char **));
+ msa->gc[tagidx] = NULL;
+ }
+ }
+
+ if (tagidx == msa->ngc) {
+ msa->gc_tag[tagidx] = sre_strdup(tag, -1);
+ msa->ngc++;
+ }
+ sre_strcat(&(msa->gc[tagidx]), -1, value, -1);
+ return;
+}
+
+/* Function: MSAGetGC()
+ * Date: SRE, Fri Aug 13 13:25:57 1999 [St. Louis]
+ *
+ * Purpose: Given a tagname for a miscellaneous #=GC column
+ * annotation, return a pointer to the annotation
+ * string.
+ *
+ * Args: msa - alignment and its annotation
+ * tag - name of the annotation
+ *
+ * Returns: ptr to the annotation string. Caller does *not*
+ * free; is managed by msa object still.
+ */
+char *
+MSAGetGC(MSA *msa, char *tag)
+{
+ int tagidx;
+
+ if (msa->gc_idx == NULL) return NULL;
+ if ((tagidx = GKIKeyIndex(msa->gc_idx, tag)) < 0) return NULL;
+ return msa->gc[tagidx];
+}
+
+
+/* Function: MSAAppendGR()
+ * Date: SRE, Thu Jun 3 06:34:38 1999 [Madison]
+ *
+ * Purpose: Add an unparsed #=GR markup line to the
+ * MSA structure, allocating as necessary.
+ *
+ * When called multiple times for the same tag,
+ * appends value strings together -- used when
+ * parsing multiblock alignment files, for
+ * example.
+ *
+ * Args: msa - multiple alignment structure
+ * tag - markup tag (e.g. "SS")
+ * sqidx - index of seq to assoc markup with (0..nseq-1)
+ * value - markup, one char per aligned column
+ *
+ * Returns: (void)
+ */
+void
+MSAAppendGR(MSA *msa, char *tag, int sqidx, char *value)
+{
+ int tagidx;
+ int i;
+
+ /* Is this an unparsed tag name that we recognize?
+ * If not, handle adding it to index, and reallocating
+ * as needed.
+ */
+ if (msa->gr_tag == NULL) /* first tag? init w/ malloc */
+ {
+ msa->gr_tag = MallocOrDie(sizeof(char *));
+ msa->gr = MallocOrDie(sizeof(char **));
+ msa->gr[0] = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (i = 0; i < msa->nseqalloc; i++)
+ msa->gr[0][i] = NULL;
+ msa->gr_idx = GKIInit();
+ tagidx = GKIStoreKey(msa->gr_idx, tag);
+ SQD_DASSERT1((tagidx == 0));
+ }
+ else
+ {
+ /* new tag? */
+ tagidx = GKIKeyIndex(msa->gr_idx, tag);
+ if (tagidx < 0) { /* it's a new tag name; realloc */
+ tagidx = GKIStoreKey(msa->gr_idx, tag);
+ /* since we alloc in blocks of 1,
+ we always realloc upon seeing
+ a new tag. */
+ SQD_DASSERT1((tagidx == msa->ngr));
+ msa->gr_tag = ReallocOrDie(msa->gr_tag, (msa->ngr+1) * sizeof(char *));
+ msa->gr = ReallocOrDie(msa->gr, (msa->ngr+1) * sizeof(char **));
+ msa->gr[msa->ngr] = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ for (i = 0; i < msa->nseqalloc; i++)
+ msa->gr[msa->ngr][i] = NULL;
+ }
+ }
+
+ if (tagidx == msa->ngr) {
+ msa->gr_tag[tagidx] = sre_strdup(tag, -1);
+ msa->ngr++;
+ }
+ sre_strcat(&(msa->gr[tagidx][sqidx]), -1, value, -1);
+ return;
+}
+
+
+/* Function: MSAVerifyParse()
+ * Date: SRE, Sat Jun 5 14:24:24 1999 [Madison, 1999 worm mtg]
+ *
+ * Purpose: Last function called after a multiple alignment is
+ * parsed. Checks that parse was successful; makes sure
+ * required information is present; makes sure required
+ * information is consistent. Some fields that are
+ * only use during parsing may be freed (sqlen, for
+ * example).
+ *
+ * Some fields in msa may be modified (msa->alen is set,
+ * for example).
+ *
+ * Args: msa - the multiple alignment
+ * sqname, aseq must be set
+ * nseq must be correct
+ * alen need not be set; will be set here.
+ * wgt will be set here if not already set
+ *
+ * Returns: (void)
+ * Will Die() here with diagnostics on error.
+ *
+ * Example:
+ */
+void
+MSAVerifyParse(MSA *msa)
+{
+ int idx;
+
+ if (msa->nseq == 0) Die("Parse error: no sequences were found for alignment %s",
+ msa->name != NULL ? msa->name : "");
+
+ msa->alen = msa->sqlen[0];
+
+ /* We can rely on msa->sqname[] being valid for any index,
+ * because of the way the line parsers always store any name
+ * they add to the index.
+ */
+ for (idx = 0; idx < msa->nseq; idx++)
+ {
+ /* aseq is required. */
+ if (msa->aseq[idx] == NULL)
+ Die("Parse error: No sequence for %s in alignment %s", msa->sqname[idx],
+ msa->name != NULL ? msa->name : "");
+ /* either all weights must be set, or none of them */
+ if ((msa->flags & MSA_SET_WGT) && msa->wgt[idx] == -1.0)
+ Die("Parse error: some weights are set, but %s doesn't have one in alignment %s",
+ msa->sqname[idx],
+ msa->name != NULL ? msa->name : "");
+ /* all aseq must be same length. */
+ if (msa->sqlen[idx] != msa->alen)
+ Die("Parse error: sequence %s: length %d, expected %d in alignment %s",
+ msa->sqname[idx], msa->sqlen[idx], msa->alen,
+ msa->name != NULL ? msa->name : "");
+ /* if SS is present, must have length right */
+ if (msa->ss != NULL && msa->ss[idx] != NULL && msa->sslen[idx] != msa->alen)
+ Die("Parse error: #=GR SS annotation for %s: length %d, expected %d in alignment %s",
+ msa->sqname[idx], msa->sslen[idx], msa->alen,
+ msa->name != NULL ? msa->name : "");
+ /* if SA is present, must have length right */
+ if (msa->sa != NULL && msa->sa[idx] != NULL && msa->salen[idx] != msa->alen)
+ Die("Parse error: #=GR SA annotation for %s: length %d, expected %d in alignment %s",
+ msa->sqname[idx], msa->salen[idx], msa->alen,
+ msa->name != NULL ? msa->name : "");
+ }
+
+ /* if cons SS is present, must have length right */
+ if (msa->ss_cons != NULL && strlen(msa->ss_cons) != msa->alen)
+ Die("Parse error: #=GC SS_cons annotation: length %d, expected %d in alignment %s",
+ strlen(msa->ss_cons), msa->alen,
+ msa->name != NULL ? msa->name : "");
+
+ /* if cons SA is present, must have length right */
+ if (msa->sa_cons != NULL && strlen(msa->sa_cons) != msa->alen)
+ Die("Parse error: #=GC SA_cons annotation: length %d, expected %d in alignment %s",
+ strlen(msa->sa_cons), msa->alen,
+ msa->name != NULL ? msa->name : "");
+
+ /* if RF is present, must have length right */
+ if (msa->rf != NULL && strlen(msa->rf) != msa->alen)
+ Die("Parse error: #=GC RF annotation: length %d, expected %d in alignment %s",
+ strlen(msa->rf), msa->alen,
+ msa->name != NULL ? msa->name : "");
+
+ /* Check that all or no weights are set */
+ if (!(msa->flags & MSA_SET_WGT))
+ FSet(msa->wgt, msa->nseq, 1.0); /* default weights */
+
+ /* Clean up a little from the parser */
+ if (msa->sqlen != NULL) { free(msa->sqlen); msa->sqlen = NULL; }
+ if (msa->sslen != NULL) { free(msa->sslen); msa->sslen = NULL; }
+ if (msa->salen != NULL) { free(msa->salen); msa->salen = NULL; }
+
+ return;
+}
+
+
+
+
+/* Function: MSAFileOpen()
+ * Date: SRE, Tue May 18 13:22:01 1999 [St. Louis]
+ *
+ * Purpose: Open an alignment database file and prepare
+ * for reading one alignment, or sequentially
+ * in the (rare) case of multiple MSA databases
+ * (e.g. Stockholm format).
+ *
+ * Args: filename - name of file to open
+ * if "-", read stdin
+ * if it ends in ".gz", read from pipe to gunzip -dc
+ * format - format of file (e.g. MSAFILE_STOCKHOLM)
+ * env - environment variable for path (e.g. BLASTDB)
+ *
+ * Returns: opened MSAFILE * on success.
+ * NULL on failure:
+ * usually, because the file doesn't exist;
+ * for gzip'ed files, may also mean that gzip isn't in the path.
+ */
+MSAFILE *
+MSAFileOpen(char *filename, int format, char *env)
+{
+ MSAFILE *afp;
+
+ afp = MallocOrDie(sizeof(MSAFILE));
+ if (strcmp(filename, "-") == 0)
+ {
+ afp->f = stdin;
+ afp->do_stdin = TRUE;
+ afp->do_gzip = FALSE;
+ afp->fname = sre_strdup("[STDIN]", -1);
+ afp->ssi = NULL; /* can't index stdin because we can't seek*/
+ }
+#ifndef SRE_STRICT_ANSI
+ /* popen(), pclose() aren't portable to non-POSIX systems; disable */
+ else if (Strparse("^.*\\.gz$", filename, 0))
+ {
+ char cmd[256];
+
+ /* Note that popen() will return "successfully"
+ * if file doesn't exist, because gzip works fine
+ * and prints an error! So we have to check for
+ * existence of file ourself.
+ */
+ if (! FileExists(filename))
+ Die("%s: file does not exist", filename);
+ if (strlen(filename) + strlen("gzip -dc ") >= 256)
+ Die("filename > 255 char in MSAFileOpen()");
+ sprintf(cmd, "gzip -dc %s", filename);
+ if ((afp->f = popen(cmd, "r")) == NULL)
+ return NULL;
+
+ afp->do_stdin = FALSE;
+ afp->do_gzip = TRUE;
+ afp->fname = sre_strdup(filename, -1);
+ /* we can't index a .gz file, because we can't seek in a pipe afaik */
+ afp->ssi = NULL;
+ }
+#endif /*SRE_STRICT_ANSI*/
+ else
+ {
+ char *ssifile;
+ char *dir;
+
+ /* When we open a file, it may be either in the current
+ * directory, or in the directory indicated by the env
+ * argument - and we have to construct the SSI filename accordingly.
+ */
+ if ((afp->f = fopen(filename, "r")) != NULL)
+ {
+ ssifile = MallocOrDie(sizeof(char) * (strlen(filename) + 5));
+ sprintf(ssifile, "%s.ssi", filename);
+ }
+ else if ((afp->f = EnvFileOpen(filename, env, &dir)) != NULL)
+ {
+ char *full;
+ full = FileConcat(dir, filename);
+ ssifile = MallocOrDie(sizeof(char) * (strlen(full) + strlen(filename) + 5));
+ sprintf(ssifile, "%s.ssi", full);
+ free(dir);
+ }
+ else return NULL;
+
+ afp->do_stdin = FALSE;
+ afp->do_gzip = FALSE;
+ afp->fname = sre_strdup(filename, -1);
+ afp->ssi = NULL;
+
+ /* Open the SSI index file. If it doesn't exist, or
+ * it's corrupt, or some error happens, afp->ssi stays NULL.
+ */
+ SSIOpen(ssifile, &(afp->ssi));
+ free(ssifile);
+ }
+
+ /* Invoke autodetection if we haven't already been told what
+ * to expect.
+ */
+ if (format == MSAFILE_UNKNOWN)
+ {
+ if (afp->do_stdin == TRUE || afp->do_gzip)
+ Die("Can't autodetect alignment file format from a stdin or gzip pipe");
+ format = MSAFileFormat(afp);
+ if (format == MSAFILE_UNKNOWN)
+ Die("Can't determine format of multiple alignment file %s", afp->fname);
+ }
+
+ afp->format = format;
+ afp->linenumber = 0;
+ afp->buf = NULL;
+ afp->buflen = 0;
+
+ return afp;
+}
+
+
+/* Function: MSAFilePositionByKey()
+ * MSAFilePositionByIndex()
+ * MSAFileRewind()
+ *
+ * Date: SRE, Tue Nov 9 19:02:54 1999 [St. Louis]
+ *
+ * Purpose: Family of functions for repositioning in
+ * open MSA files; analogous to a similarly
+ * named function series in HMMER's hmmio.c.
+ *
+ * Args: afp - open alignment file
+ * offset - disk offset in bytes
+ * key - key to look up in SSI indices
+ * idx - index of alignment.
+ *
+ * Returns: 0 on failure.
+ * 1 on success.
+ * If called on a non-fseek()'able file (e.g. a gzip'ed
+ * or pipe'd alignment), returns 0 as a failure flag.
+ */
+int
+MSAFileRewind(MSAFILE *afp)
+{
+ if (afp->do_gzip || afp->do_stdin) return 0;
+ rewind(afp->f);
+ return 1;
+}
+int
+MSAFilePositionByKey(MSAFILE *afp, char *key)
+{
+ int fh; /* filehandle is ignored */
+ SSIOFFSET offset; /* offset of the key alignment */
+
+ if (afp->ssi == NULL) return 0;
+ if (SSIGetOffsetByName(afp->ssi, key, &fh, &offset) != 0) return 0;
+ if (SSISetFilePosition(afp->f, &offset) != 0) return 0;
+ return 1;
+}
+int
+MSAFilePositionByIndex(MSAFILE *afp, int idx)
+{
+ int fh; /* filehandled is passed but ignored */
+ SSIOFFSET offset; /* disk offset of desired alignment */
+
+ if (afp->ssi == NULL) return 0;
+ if (SSIGetOffsetByNumber(afp->ssi, idx, &fh, &offset) != 0) return 0;
+ if (SSISetFilePosition(afp->f, &offset) != 0) return 0;
+ return 1;
+}
+
+
+/* Function: MSAFileRead()
+ * Date: SRE, Fri May 28 16:01:43 1999 [St. Louis]
+ *
+ * Purpose: Read the next msa from an open alignment file.
+ * This is a wrapper around format-specific calls.
+ *
+ * Args: afp - open alignment file
+ *
+ * Returns: next alignment, or NULL if out of alignments
+ */
+MSA *
+MSAFileRead(MSAFILE *afp)
+{
+ MSA *msa = NULL;
+
+ switch (afp->format) {
+ case MSAFILE_STOCKHOLM: msa = ReadStockholm(afp); break;
+ case MSAFILE_MSF: msa = ReadMSF(afp); break;
+ case MSAFILE_A2M: msa = ReadA2M(afp); break;
+ case MSAFILE_CLUSTAL: msa = ReadClustal(afp); break;
+ case MSAFILE_SELEX: msa = ReadSELEX(afp); break;
+ case MSAFILE_PHYLIP: msa = ReadPhylip(afp); break;
+ default:
+ Die("MSAFILE corrupted: bad format index");
+ }
+ return msa;
+}
+
+/* Function: MSAFileClose()
+ * Date: SRE, Tue May 18 14:05:28 1999 [St. Louis]
+ *
+ * Purpose: Close an open MSAFILE.
+ *
+ * Args: afp - ptr to an open MSAFILE.
+ *
+ * Returns: void
+ */
+void
+MSAFileClose(MSAFILE *afp)
+{
+#ifndef SRE_STRICT_ANSI /* gzip functionality only on POSIX systems */
+ if (afp->do_gzip) pclose(afp->f);
+#endif
+ if (! afp->do_stdin) fclose(afp->f);
+ if (afp->buf != NULL) free(afp->buf);
+ if (afp->ssi != NULL) SSIClose(afp->ssi);
+ if (afp->fname != NULL) free(afp->fname);
+ free(afp);
+}
+
+char *
+MSAFileGetLine(MSAFILE *afp)
+{
+ char *s;
+ if ((s = sre_fgets(&(afp->buf), &(afp->buflen), afp->f)) == NULL)
+ return NULL;
+ afp->linenumber++;
+ return afp->buf;
+}
+
+void
+MSAFileWrite(FILE *fp, MSA *msa, int outfmt, int do_oneline)
+{
+ switch (outfmt) {
+#ifdef CLUSTALO
+ case MSAFILE_A2M: WriteA2M(stdout, msa, 0); break;
+ case MSAFILE_VIENNA: WriteA2M(stdout, msa, 1); break;
+#else
+ case MSAFILE_A2M: WriteA2M(stdout, msa); break;
+#endif
+ case MSAFILE_CLUSTAL: WriteClustal(stdout, msa); break;
+ case MSAFILE_MSF: WriteMSF(stdout, msa); break;
+ case MSAFILE_PHYLIP: WritePhylip(stdout, msa); break;
+ case MSAFILE_SELEX: WriteSELEX(stdout, msa); break;
+ case MSAFILE_STOCKHOLM:
+ if (do_oneline) WriteStockholmOneBlock(stdout, msa);
+ else WriteStockholm(stdout, msa);
+ break;
+ default:
+ Die("can't write. no such alignment format %d\n", outfmt);
+ }
+}
+
+/* Function: MSAGetSeqidx()
+ * Date: SRE, Wed May 19 15:08:25 1999 [St. Louis]
+ *
+ * Purpose: From a sequence name, return seqidx appropriate
+ * for an MSA structure.
+ *
+ * 1) try to guess the index. (pass -1 if you can't guess)
+ * 2) Look up name in msa's hashtable.
+ * 3) If it's a new name, store in msa's hashtable;
+ * expand allocs as needed;
+ * save sqname.
+ *
+ * Args: msa - alignment object
+ * name - a sequence name
+ * guess - a guess at the right index, or -1 if no guess.
+ *
+ * Returns: seqidx
+ */
+int
+MSAGetSeqidx(MSA *msa, char *name, int guess)
+{
+ int seqidx;
+ /* can we guess? */
+ if (guess >= 0 && guess < msa->nseq && strcmp(name, msa->sqname[guess]) == 0)
+ return guess;
+ /* else, a lookup in the index */
+ if ((seqidx = GKIKeyIndex(msa->index, name)) >= 0)
+ return seqidx;
+ /* else, it's a new name */
+ seqidx = GKIStoreKey(msa->index, name);
+ if (seqidx >= msa->nseqalloc) MSAExpand(msa);
+
+ msa->sqname[seqidx] = sre_strdup(name, -1);
+ msa->nseq++;
+ return seqidx;
+}
+
+
+/* Function: MSAFromAINFO()
+ * Date: SRE, Mon Jun 14 11:22:24 1999 [St. Louis]
+ *
+ * Purpose: Convert the old aseq/ainfo alignment structure
+ * to new MSA structure. Enables more rapid conversion
+ * of codebase to the new world order.
+ *
+ * Args: aseq - [0..nseq-1][0..alen-1] alignment
+ * ainfo - old-style optional info
+ *
+ * Returns: MSA *
+ */
+MSA *
+MSAFromAINFO(char **aseq, AINFO *ainfo)
+{
+ MSA *msa;
+ int i, j;
+
+ msa = MSAAlloc(ainfo->nseq, ainfo->alen);
+ for (i = 0; i < ainfo->nseq; i++)
+ {
+ strcpy(msa->aseq[i], aseq[i]);
+ msa->wgt[i] = ainfo->wgt[i];
+ msa->sqname[i] = sre_strdup(ainfo->sqinfo[i].name, -1);
+ msa->sqlen[i] = msa->alen;
+ GKIStoreKey(msa->index, msa->sqname[i]);
+
+ if (ainfo->sqinfo[i].flags & SQINFO_ACC)
+ MSASetSeqAccession(msa, i, ainfo->sqinfo[i].acc);
+
+ if (ainfo->sqinfo[i].flags & SQINFO_DESC)
+ MSASetSeqDescription(msa, i, ainfo->sqinfo[i].desc);
+
+ if (ainfo->sqinfo[i].flags & SQINFO_SS) {
+ if (msa->ss == NULL) {
+ msa->ss = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ msa->sslen = MallocOrDie(sizeof(int) * msa->nseqalloc);
+ for (j = 0; j < msa->nseqalloc; j++) {
+ msa->ss[j] = NULL;
+ msa->sslen[j] = 0;
+ }
+ }
+ MakeAlignedString(msa->aseq[i], msa->alen, ainfo->sqinfo[i].ss, &(msa->ss[i]));
+ msa->sslen[i] = msa->alen;
+ }
+
+ if (ainfo->sqinfo[i].flags & SQINFO_SA) {
+ if (msa->sa == NULL) {
+ msa->sa = MallocOrDie(sizeof(char *) * msa->nseqalloc);
+ msa->salen = MallocOrDie(sizeof(int) * msa->nseqalloc);
+ for (j = 0; j < msa->nseqalloc; j++) {
+ msa->sa[j] = NULL;
+ msa->salen[j] = 0;
+ }
+ }
+ MakeAlignedString(msa->aseq[i], msa->alen, ainfo->sqinfo[i].sa, &(msa->sa[i]));
+ msa->salen[i] = msa->alen;
+ }
+ }
+ /* note that sre_strdup() returns NULL when passed NULL */
+ msa->name = sre_strdup(ainfo->name, -1);
+ msa->desc = sre_strdup(ainfo->desc, -1);
+ msa->acc = sre_strdup(ainfo->acc, -1);
+ msa->au = sre_strdup(ainfo->au, -1);
+ msa->ss_cons = sre_strdup(ainfo->cs, -1);
+ msa->rf = sre_strdup(ainfo->rf, -1);
+ if (ainfo->flags & AINFO_TC) {
+ msa->cutoff[MSA_CUTOFF_TC1] = ainfo->tc1; msa->cutoff_is_set[MSA_CUTOFF_TC1] = TRUE;
+ msa->cutoff[MSA_CUTOFF_TC2] = ainfo->tc2; msa->cutoff_is_set[MSA_CUTOFF_TC2] = TRUE;
+ }
+ if (ainfo->flags & AINFO_NC) {
+ msa->cutoff[MSA_CUTOFF_NC1] = ainfo->nc1; msa->cutoff_is_set[MSA_CUTOFF_NC1] = TRUE;
+ msa->cutoff[MSA_CUTOFF_NC2] = ainfo->nc2; msa->cutoff_is_set[MSA_CUTOFF_NC2] = TRUE;
+ }
+ if (ainfo->flags & AINFO_GA) {
+ msa->cutoff[MSA_CUTOFF_GA1] = ainfo->ga1; msa->cutoff_is_set[MSA_CUTOFF_GA1] = TRUE;
+ msa->cutoff[MSA_CUTOFF_GA2] = ainfo->ga2; msa->cutoff_is_set[MSA_CUTOFF_GA2] = TRUE;
+ }
+ msa->nseq = ainfo->nseq;
+ msa->alen = ainfo->alen;
+ return msa;
+}
+
+
+
+
+/* Function: MSAFileFormat()
+ * Date: SRE, Fri Jun 18 14:26:49 1999 [Sanger Centre]
+ *
+ * Purpose: (Attempt to) determine the format of an alignment file.
+ * Since it rewinds the file pointer when it's done,
+ * cannot be used on a pipe or gzip'ed file. Works by
+ * calling SeqfileFormat() from sqio.c, then making sure
+ * that the format is indeed an alignment. If the format
+ * comes back as FASTA, it assumes that the format as A2M
+ * (e.g. aligned FASTA).
+ *
+ * Args: fname - file to evaluate
+ *
+ * Returns: format code; e.g. MSAFILE_STOCKHOLM
+ */
+int
+MSAFileFormat(MSAFILE *afp)
+{
+ int fmt;
+
+ fmt = SeqfileFormat(afp->f);
+
+ if (fmt == SQFILE_FASTA) fmt = MSAFILE_A2M;
+
+ if (fmt != MSAFILE_UNKNOWN && ! IsAlignmentFormat(fmt))
+ Die("File %s does not appear to be an alignment file;\n\
+rather, it appears to be an unaligned file in %s format.\n\
+I'm expecting an alignment file in this context.\n",
+ afp->fname,
+ SeqfileFormat2String(fmt));
+ return fmt;
+}
+
+
+/* Function: MSAMingap()
+ * Date: SRE, Mon Jun 28 18:57:54 1999 [on jury duty, St. Louis Civil Court]
+ *
+ * Purpose: Remove all-gap columns from a multiple sequence alignment
+ * and its associated per-residue data.
+ *
+ * Args: msa - the alignment
+ *
+ * Returns: (void)
+ */
+void
+MSAMingap(MSA *msa)
+{
+ int *useme; /* array of TRUE/FALSE flags for which columns to keep */
+ int apos; /* position in original alignment */
+ int idx; /* sequence index */
+
+ useme = MallocOrDie(sizeof(int) * msa->alen);
+ for (apos = 0; apos < msa->alen; apos++)
+ {
+ for (idx = 0; idx < msa->nseq; idx++)
+ if (! isgap(msa->aseq[idx][apos]))
+ break;
+ if (idx == msa->nseq) useme[apos] = FALSE; else useme[apos] = TRUE;
+ }
+ MSAShorterAlignment(msa, useme);
+ free(useme);
+ return;
+}
+
+/* Function: MSANogap()
+ * Date: SRE, Wed Nov 17 09:59:51 1999 [St. Louis]
+ *
+ * Purpose: Remove all columns from a multiple sequence alignment that
+ * contain any gaps -- used for filtering before phylogenetic
+ * analysis.
+ *
+ * Args: msa - the alignment
+ *
+ * Returns: (void). The alignment is modified, so if you want to keep
+ * the original for something, make a copy.
+ */
+void
+MSANogap(MSA *msa)
+{
+ int *useme; /* array of TRUE/FALSE flags for which columns to keep */
+ int apos; /* position in original alignment */
+ int idx; /* sequence index */
+
+ useme = MallocOrDie(sizeof(int) * msa->alen);
+ for (apos = 0; apos < msa->alen; apos++)
+ {
+ for (idx = 0; idx < msa->nseq; idx++)
+ if (isgap(msa->aseq[idx][apos]))
+ break;
+ if (idx == msa->nseq) useme[apos] = TRUE; else useme[apos] = FALSE;
+ }
+ MSAShorterAlignment(msa, useme);
+ free(useme);
+ return;
+}
+
+
+/* Function: MSAShorterAlignment()
+ * Date: SRE, Wed Nov 17 09:49:32 1999 [St. Louis]
+ *
+ * Purpose: Given an array "useme" (0..alen-1) of TRUE/FALSE flags,
+ * where TRUE means "keep this column in the new alignment":
+ * Remove all columns annotated as "FALSE" in the useme
+ * array.
+ *
+ * Args: msa - the alignment. The alignment is changed, so
+ * if you don't want the original screwed up, make
+ * a copy of it first.
+ * useme - TRUE/FALSE flags for columns to keep: 0..alen-1
+ *
+ * Returns: (void)
+ */
+void
+MSAShorterAlignment(MSA *msa, int *useme)
+{
+ int apos; /* position in original alignment */
+ int mpos; /* position in new alignment */
+ int idx; /* sequence index */
+ int i; /* markup index */
+
+ /* Since we're minimizing, we can overwrite, using already allocated
+ * memory.
+ */
+ for (apos = 0, mpos = 0; apos < msa->alen; apos++)
+ {
+ if (useme[apos] == FALSE) continue;
+
+ /* shift alignment and associated per-column+per-residue markup */
+ if (mpos != apos)
+ {
+ for (idx = 0; idx < msa->nseq; idx++)
+ {
+ msa->aseq[idx][mpos] = msa->aseq[idx][apos];
+ if (msa->ss != NULL && msa->ss[idx] != NULL) msa->ss[idx][mpos] = msa->ss[idx][apos];
+ if (msa->sa != NULL && msa->sa[idx] != NULL) msa->sa[idx][mpos] = msa->sa[idx][apos];
+
+ for (i = 0; i < msa->ngr; i++)
+ if (msa->gr[i][idx] != NULL) msa->gr[i][idx][mpos] = msa->gr[i][idx][apos];
+ }
+
+ if (msa->ss_cons != NULL) msa->ss_cons[mpos] = msa->ss_cons[apos];
+ if (msa->sa_cons != NULL) msa->sa_cons[mpos] = msa->sa_cons[apos];
+ if (msa->rf != NULL) msa->rf[mpos] = msa->rf[apos];
+
+ for (i = 0; i < msa->ngc; i++)
+ msa->gc[i][mpos] = msa->gc[i][apos];
+ }
+ mpos++;
+ }
+
+ msa->alen = mpos; /* set new length */
+ /* null terminate everything */
+ for (idx = 0; idx < msa->nseq; idx++)
+ {
+ msa->aseq[idx][mpos] = '\0';
+ if (msa->ss != NULL && msa->ss[idx] != NULL) msa->ss[idx][mpos] = '\0';
+ if (msa->sa != NULL && msa->sa[idx] != NULL) msa->sa[idx][mpos] = '\0';
+
+ for (i = 0; i < msa->ngr; i++)
+ if (msa->gr[i][idx] != NULL) msa->gr[i][idx][mpos] = '\0';
+ }
+
+ if (msa->ss_cons != NULL) msa->ss_cons[mpos] = '\0';
+ if (msa->sa_cons != NULL) msa->sa_cons[mpos] = '\0';
+ if (msa->rf != NULL) msa->rf[mpos] = '\0';
+
+ for (i = 0; i < msa->ngc; i++)
+ msa->gc[i][mpos] = '\0';
+
+ return;
+}
+
+
+/* Function: MSASmallerAlignment()
+ * Date: SRE, Wed Jun 30 09:56:08 1999 [St. Louis]
+ *
+ * Purpose: Given an array "useme" of TRUE/FALSE flags for
+ * each sequence in an alignment, construct
+ * and return a new alignment containing only
+ * those sequences that are flagged useme=TRUE.
+ *
+ * Used by routines such as MSAFilterAlignment()
+ * and MSASampleAlignment().
+ *
+ * Limitations:
+ * Does not copy unparsed Stockholm markup.
+ *
+ * Does not make assumptions about meaning of wgt;
+ * if you want the new wgt vector renormalized, do
+ * it yourself with FNorm(new->wgt, new->nseq).
+ *
+ * Args: msa -- the original (larger) alignment
+ * useme -- [0..nseq-1] array of TRUE/FALSE flags; TRUE means include
+ * this seq in new alignment
+ * ret_new -- RETURN: new alignment
+ *
+ * Returns: void
+ * ret_new is allocated here; free with MSAFree()
+ */
+void
+MSASmallerAlignment(MSA *msa, int *useme, MSA **ret_new)
+{
+ MSA *new; /* RETURN: new alignment */
+ int nnew; /* number of seqs in new msa (e.g. # of TRUEs) */
+ int oidx, nidx; /* old, new indices */
+ int i;
+
+ nnew = 0;
+ for (oidx = 0; oidx < msa->nseq; oidx++)
+ if (useme[oidx]) nnew++;
+ if (nnew == 0) { *ret_new = NULL; return; }
+
+ new = MSAAlloc(nnew, 0);
+ nidx = 0;
+ for (oidx = 0; oidx < msa->nseq; oidx++)
+ if (useme[oidx])
+ {
+ new->aseq[nidx] = sre_strdup(msa->aseq[oidx], msa->alen);
+ new->sqname[nidx] = sre_strdup(msa->sqname[oidx], msa->alen);
+ GKIStoreKey(new->index, msa->sqname[oidx]);
+ new->wgt[nidx] = msa->wgt[oidx];
+ if (msa->sqacc != NULL)
+ MSASetSeqAccession(new, nidx, msa->sqacc[oidx]);
+ if (msa->sqdesc != NULL)
+ MSASetSeqDescription(new, nidx, msa->sqdesc[oidx]);
+ if (msa->ss != NULL && msa->ss[oidx] != NULL)
+ {
+ if (new->ss == NULL) new->ss = MallocOrDie(sizeof(char *) * new->nseq);
+ new->ss[nidx] = sre_strdup(msa->ss[oidx], -1);
+ }
+ if (msa->sa != NULL && msa->sa[oidx] != NULL)
+ {
+ if (new->sa == NULL) new->sa = MallocOrDie(sizeof(char *) * new->nseq);
+ new->sa[nidx] = sre_strdup(msa->sa[oidx], -1);
+ }
+ nidx++;
+ }
+
+ new->nseq = nnew;
+ new->alen = msa->alen;
+ new->flags = msa->flags;
+ new->type = msa->type;
+ new->name = sre_strdup(msa->name, -1);
+ new->desc = sre_strdup(msa->desc, -1);
+ new->acc = sre_strdup(msa->acc, -1);
+ new->au = sre_strdup(msa->au, -1);
+ new->ss_cons = sre_strdup(msa->ss_cons, -1);
+ new->sa_cons = sre_strdup(msa->sa_cons, -1);
+ new->rf = sre_strdup(msa->rf, -1);
+ for (i = 0; i < MSA_MAXCUTOFFS; i++) {
+ new->cutoff[i] = msa->cutoff[i];
+ new->cutoff_is_set[i] = msa->cutoff_is_set[i];
+ }
+ free(new->sqlen);
+
+ MSAMingap(new);
+ *ret_new = new;
+ return;
+}
+
+
+/*****************************************************************
+ * Retrieval routines
+ *
+ * Access to MSA structure data is possible through these routines.
+ * I'm not doing this because of object oriented design, though
+ * it might work in my favor someday.
+ * I'm doing this because lots of MSA data is optional, and
+ * checking through the chain of possible NULLs is a pain.
+ *****************************************************************/
+
+char *
+MSAGetSeqAccession(MSA *msa, int idx)
+{
+ if (msa->sqacc != NULL && msa->sqacc[idx] != NULL)
+ return msa->sqacc[idx];
+ else
+ return NULL;
+}
+char *
+MSAGetSeqDescription(MSA *msa, int idx)
+{
+ if (msa->sqdesc != NULL && msa->sqdesc[idx] != NULL)
+ return msa->sqdesc[idx];
+ else
+ return NULL;
+}
+char *
+MSAGetSeqSS(MSA *msa, int idx)
+{
+ if (msa->ss != NULL && msa->ss[idx] != NULL)
+ return msa->ss[idx];
+ else
+ return NULL;
+}
+char *
+MSAGetSeqSA(MSA *msa, int idx)
+{
+ if (msa->sa != NULL && msa->sa[idx] != NULL)
+ return msa->sa[idx];
+ else
+ return NULL;
+}
+
+
+/*****************************************************************
+ * Information routines
+ *
+ * Access information about the MSA.
+ *****************************************************************/
+
+/* Function: MSAAverageSequenceLength()
+ * Date: SRE, Sat Apr 6 09:41:34 2002 [St. Louis]
+ *
+ * Purpose: Return the average length of the (unaligned) sequences
+ * in the MSA.
+ *
+ * Args: msa - the alignment
+ *
+ * Returns: average length
+ */
+float
+MSAAverageSequenceLength(MSA *msa)
+{
+ int i;
+ float avg;
+
+ avg = 0.;
+ for (i = 0; i < msa->nseq; i++)
+ avg += (float) DealignedLength(msa->aseq[i]);
+
+ if (msa->nseq == 0) return 0.;
+ else return (avg / msa->nseq);
+}
+
+
git://source.jalview.org / jabaws.git / blobdiff