+++ /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.
- *****************************************************************/
-
-/* weight.c
- * SRE, Thu Mar 3 07:56:01 1994
- *
- * Calculate weights for sequences in an alignment.
- * RCS $Id: weight.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: weight.c,v 1.9 2002/10/09 14:26:09 eddy Exp)
- */
-
-#include <ctype.h>
-#include <string.h>
-#include "squid.h"
-#include "sre_random.h"
-
-static void upweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, int node);
-static void downweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt,
- float *fwt, int node);
-static float simple_distance(char *s1, char *s2);
-static int simple_diffmx(char **aseqs,int num, float ***ret_dmx);
-
-/* Function: GSCWeights()
- *
- * Purpose: Use Erik's tree-based algorithm to set weights for
- * sequences in an alignment. upweight() and downweight()
- * are derived from Graeme Mitchison's code.
- *
- * Args: aseq - array of (0..nseq-1) aligned sequences
- * nseq - number of seqs in alignment
- * alen - length of alignment
- * wgt - allocated [0..nseq-1] array of weights to be returned
- *
- * Return: (void)
- * wgt is filled in.
- */
-void
-GSCWeights(char **aseq, int nseq, int alen, float *wgt)
-{
- float **dmx; /* distance (difference) matrix */
- struct phylo_s *tree;
- float *lwt, *rwt; /* weight on left, right of this tree node */
- float *fwt; /* final weight assigned to this node */
- int i;
-
- /* Sanity check first
- */
- if (nseq == 1) { wgt[0] = 1.0; return; }
-
- /* I use a simple fractional difference matrix derived by
- * pairwise identity. Perhaps I should include a Poisson
- * distance correction.
- */
- MakeDiffMx(aseq, nseq, &dmx);
- if (! Cluster(dmx, nseq, CLUSTER_MIN, &tree)) Die("Cluster() failed");
-
- /* Allocations
- */
- lwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
- rwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
- fwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
-
- /* lwt and rwt are the total branch weight to the left and
- * right of a node or sequence. They are 0..2N-2. 0..N-1 are
- * the sequences; these have weight 0. N..2N-2 are the actual
- * tree nodes.
- */
- for (i = 0; i < nseq; i++)
- lwt[i] = rwt[i] = 0.0;
- /* recursively calculate rwt, lwt, starting
- at node nseq (the root) */
- upweight(tree, nseq, lwt, rwt, nseq);
-
- /* recursively distribute weight across the
- tree */
- fwt[nseq] = nseq;
- downweight(tree, nseq, lwt, rwt, fwt, nseq);
- /* collect the weights */
- for (i = 0; i < nseq; i++)
- wgt[i] = fwt[i];
-
- FMX2Free(dmx);
- FreePhylo(tree, nseq);
- free(lwt); free(rwt); free(fwt);
-}
-
-static void
-upweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, int node)
-{
- int ld,rd;
-
- ld = tree[node-nseq].left;
- if (ld >= nseq) upweight(tree, nseq, lwt, rwt, ld);
- rd = tree[node-nseq].right;
- if (rd >= nseq) upweight(tree, nseq, lwt, rwt, rd);
- lwt[node] = lwt[ld] + rwt[ld] + tree[node-nseq].lblen;
- rwt[node] = lwt[rd] + rwt[rd] + tree[node-nseq].rblen;
-}
-
-
-static void
-downweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, float *fwt, int node)
-{
- int ld,rd;
- float lnum, rnum;
-
- ld = tree[node-nseq].left;
- rd = tree[node-nseq].right;
- if (lwt[node] + rwt[node] > 0.0)
- {
- fwt[ld] = fwt[node] * (lwt[node] / (lwt[node] + rwt[node]));
- fwt[rd] = fwt[node] * (rwt[node] / (lwt[node] + rwt[node]));
- }
- else
- {
- lnum = (ld >= nseq) ? tree[ld-nseq].incnum : 1.0;
- rnum = (rd >= nseq) ? tree[rd-nseq].incnum : 1.0;
- fwt[ld] = fwt[node] * lnum / (lnum + rnum);
- fwt[rd] = fwt[node] * rnum / (lnum + rnum);
- }
-
- if (ld >= nseq) downweight(tree, nseq, lwt, rwt, fwt, ld);
- if (rd >= nseq) downweight(tree, nseq, lwt, rwt, fwt, rd);
-}
-
-
-
-
-/* Function: VoronoiWeights()
- *
- * Purpose: Calculate weights using the scheme of Sibbald &
- * Argos (JMB 216:813-818 1990). The scheme is
- * slightly modified because the original algorithm
- * actually doesn't work on gapped alignments.
- * The sequences are assumed to be protein.
- *
- * Args: aseq - array of (0..nseq-1) aligned sequences
- * nseq - number of sequences
- * alen - length of alignment
- * wgt - allocated [0..nseq-1] array of weights to be returned
- *
- * Return: void
- * wgt is filled in.
- */
-void
-VoronoiWeights(char **aseq, int nseq, int alen, float *wgt)
-{
- float **dmx; /* distance (difference) matrix */
- float *halfmin; /* 1/2 minimum distance to other seqs */
- char **psym; /* symbols seen in each column */
- int *nsym; /* # syms seen in each column */
- int symseen[27]; /* flags for observed syms */
- char *randseq; /* randomly generated sequence */
- int acol; /* pos in aligned columns */
- int idx; /* index in sequences */
- int symidx; /* 0..25 index for symbol */
- int i; /* generic counter */
- float min; /* minimum distance */
- float dist; /* distance between random and real */
- float challenge, champion; /* for resolving ties */
- int itscale; /* how many iterations per seq */
- int iteration;
- int best; /* index of nearest real sequence */
-
- /* Sanity check first
- */
- if (nseq == 1) { wgt[0] = 1.0; return; }
-
- itscale = 50;
-
- /* Precalculate 1/2 minimum distance to other
- * sequences for each sequence
- */
- if (! simple_diffmx(aseq, nseq, &dmx))
- Die("simple_diffmx() failed");
- halfmin = MallocOrDie (sizeof(float) * nseq);
- for (idx = 0; idx < nseq; idx++)
- {
- for (min = 1.0, i = 0; i < nseq; i++)
- {
- if (i == idx) continue;
- if (dmx[idx][i] < min) min = dmx[idx][i];
- }
- halfmin[idx] = min / 2.0;
- }
- Free2DArray((void **) dmx, nseq);
-
- /* Set up the random sequence generating model.
- */
- psym = MallocOrDie (alen * sizeof(char *));
- nsym = MallocOrDie (alen * sizeof(int));
- for (acol = 0; acol < alen; acol++)
- psym[acol] = MallocOrDie (27 * sizeof(char));
-
-/* #ifdef ORIGINAL_SIBBALD_ALGORITHM_IS_BROKEN */
- for (acol = 0; acol < alen; acol++)
- {
- memset(symseen, 0, sizeof(int) * 27);
- for (idx = 0; idx < nseq; idx++)
- if (! isgap(aseq[idx][acol]))
- {
- if (isupper((int) aseq[idx][acol]))
- symidx = aseq[idx][acol] - 'A';
- else
- symidx = aseq[idx][acol] - 'a';
- if (symidx >= 0 && symidx < 26)
- symseen[symidx] = 1;
- }
- else
- symseen[26] = 1; /* a gap */
-
- for (nsym[acol] = 0, i = 0; i < 26; i++)
- if (symseen[i])
- {
- psym[acol][nsym[acol]] = 'A'+i;
- nsym[acol]++;
- }
- if (symseen[26]) { psym[acol][nsym[acol]] = ' '; nsym[acol]++; }
- }
-/* #endif ORIGINAL_SIBBALD_ALGORITHM_IS_BROKEN */
-
- /* Note: the original Sibbald&Argos algorithm calls for
- * bounding the sampled space using a template-like random
- * sequence generator. However, this leads to one minor
- * and one major problem. The minor problem is that
- * exceptional amino acids in a column can have a
- * significant effect by altering the amount of sampled
- * sequence space; the larger the data set, the worse
- * this problem becomes. The major problem is that
- * there is no reasonable way to deal with gaps.
- * Gapped sequences simply inhabit a different dimensionality
- * and it's pretty painful to imagine calculating Voronoi
- * volumes when the N in your N-space is varying.
- * Note that all the examples shown by Sibbald and Argos
- * are *ungapped* examples.
- *
- * The best way I've found to circumvent this problem is
- * just not to bound the sampled space; count gaps as
- * symbols and generate completely random sequences.
- */
-#ifdef ALL_SEQUENCE_SPACE
- for (acol = 0; acol < alen; acol++)
- {
- strcpy(psym[acol], "ACDEFGHIKLMNPQRSTVWY ");
- nsym[acol] = 21;
- }
-#endif
-
- /* Sibbald and Argos algorithm:
- * 1) assign all seqs weight 0.
- * 2) generate a "random" sequence
- * 3) calculate distance to every other sequence
- * (if we get a distance < 1/2 minimum distance
- * to other real seqs, we can stop)
- * 4) if unique closest sequence, increment its weight 1.
- * if multiple closest seq, choose one randomly
- * 5) repeat 2-4 for lots of iterations
- * 6) normalize all weights to sum to nseq.
- */
- randseq = MallocOrDie ((alen+1) * sizeof(char));
-
- best = 42.; /* solely to silence GCC uninit warnings. */
- FSet(wgt, nseq, 0.0);
- for (iteration = 0; iteration < itscale * nseq; iteration++)
- {
- for (acol = 0; acol < alen; acol++)
- randseq[acol] = (nsym[acol] == 0) ? ' ' : psym[acol][CHOOSE(nsym[acol])];
- randseq[acol] = '\0';
-
- champion = sre_random();
- for (min = 1.0, idx = 0; idx < nseq; idx++)
- {
- dist = simple_distance(aseq[idx], randseq);
- if (dist < halfmin[idx])
- {
- best = idx;
- break;
- }
- if (dist < min)
- { champion = sre_random(); best = idx; min = dist; }
- else if (dist == min)
- {
- challenge = sre_random();
- if (challenge > champion)
- { champion = challenge; best = idx; min = dist; }
- }
- }
- wgt[best] += 1.0;
- }
-
- for (idx = 0; idx < nseq; idx++)
- wgt[idx] = wgt[idx] / (float) itscale;
-
- free(randseq);
- free(nsym);
- free(halfmin);
- Free2DArray((void **) psym, alen);
-}
-
-
-/* Function: simple_distance()
- *
- * Purpose: For two identical-length null-terminated strings, return
- * the fractional difference between them. (0..1)
- * (Gaps don't count toward anything.)
- */
-static float
-simple_distance(char *s1, char *s2)
-{
- int diff = 0;
- int valid = 0;
-
- for (; *s1 != '\0'; s1++, s2++)
- {
- if (isgap(*s1) || isgap(*s2)) continue;
- if (*s1 != *s2) diff++;
- valid++;
- }
- return (valid > 0 ? ((float) diff / (float) valid) : 0.0);
-}
-
-/* Function: simple_diffmx()
- *
- * Purpose: Given a set of flushed, aligned sequences, construct
- * an NxN fractional difference matrix using the
- * simple_distance rule.
- *
- * Args: aseqs - flushed, aligned sequences
- * num - number of aseqs
- * ret_dmx - RETURN: difference matrix (caller must free)
- *
- * Return: 1 on success, 0 on failure.
- */
-static int
-simple_diffmx(char **aseqs,
- int num,
- float ***ret_dmx)
-{
- float **dmx; /* RETURN: distance matrix */
- int i,j; /* counters over sequences */
-
- /* Allocate
- */
- if ((dmx = (float **) malloc (sizeof(float *) * num)) == NULL)
- Die("malloc failed");
- for (i = 0; i < num; i++)
- if ((dmx[i] = (float *) malloc (sizeof(float) * num)) == NULL)
- Die("malloc failed");
-
- /* Calculate distances, symmetric matrix
- */
- for (i = 0; i < num; i++)
- for (j = i; j < num; j++)
- dmx[i][j] = dmx[j][i] = simple_distance(aseqs[i], aseqs[j]);
-
- /* Return
- */
- *ret_dmx = dmx;
- return 1;
-}
-
-
-
-/* Function: BlosumWeights()
- * Date: SRE, Fri Jul 16 17:33:59 1999 (St. Louis)
- *
- * Purpose: Assign weights to a set of aligned sequences
- * using the BLOSUM rule:
- * - do single linkage clustering at some pairwise identity
- * - in each cluster, give each sequence 1/clustsize
- * total weight.
- *
- * The clusters have no pairwise link >= maxid.
- *
- * O(N) in memory. Probably ~O(NlogN) in time; O(N^2)
- * in worst case, which is no links between sequences
- * (e.g., values of maxid near 1.0).
- *
- * Args: aseqs - alignment
- * nseq - number of seqs in alignment
- * alen - # of columns in alignment
- * maxid - fractional identity (e.g. 0.62 for BLOSUM62)
- * wgt - [0..nseq-1] array of weights to be returned
- */
-void
-BlosumWeights(char **aseqs, int nseq, int alen, float maxid, float *wgt)
-{
- int *c, nc;
- int *nmem; /* number of seqs in each cluster */
- int i; /* loop counter */
-
- SingleLinkCluster(aseqs, nseq, alen, maxid, &c, &nc);
-
- FSet(wgt, nseq, 1.0);
- nmem = MallocOrDie(sizeof(int) * nc);
-
- for (i = 0; i < nc; i++) nmem[i] = 0;
- for (i = 0; i < nseq; i++) nmem[c[i]]++;
- for (i = 0; i < nseq; i++) wgt[i] = 1. / (float) nmem[c[i]];
-
- free(nmem);
- free(c);
- return;
-}
-
-
-/* Function: PositionBasedWeights()
- * Date: SRE, Fri Jul 16 17:47:22 1999 [St. Louis]
- *
- * Purpose: Implementation of Henikoff and Henikoff position-based
- * weights (JMB 243:574-578, 1994) [Henikoff94b].
- *
- * A significant advantage of this approach that Steve and Jorja
- * don't point out is that it is O(N) in memory, unlike
- * many other approaches like GSC weights or Voronoi.
- *
- * A potential disadvantage that they don't point out
- * is that in the theoretical limit of infinite sequences
- * in the alignment, weights go flat: eventually every
- * column has at least one representative of each of 20 aa (or 4 nt)
- * in it.
- *
- * They also don't give a rule for how to handle gaps.
- * The rule used here seems the obvious and sensible one
- * (ignore them). This means that longer sequences
- * initially get more weight; hence a "double
- * normalization" in which the weights are first divided
- * by sequence length (to compensate for that effect),
- * then normalized to sum to nseq.
- *
- * Limitations:
- * Implemented in a way that's alphabet-independent:
- * it uses the 26 upper case letters as "residues".
- * Any alphabetic character in aseq is interpreted as
- * a unique "residue" (case insensitively; lower case
- * mapped to upper case). All other characters are
- * interpreted as gaps.
- *
- * This way, we don't have to pass around any alphabet
- * type info (DNA vs. RNA vs. protein) and don't have
- * to deal with remapping IUPAC degenerate codes
- * probabilistically. However, on the down side,
- * a sequence with a lot of degenerate IUPAC characters
- * will get an artifactually high PB weight.
- *
- * Args: aseq - sequence alignment to weight
- * nseq - number of sequences in alignment
- * alen - length of alignment
- * wgt - RETURN: weights filled in (pre-allocated 0..nseq-1)
- *
- * Returns: (void)
- * wgt is allocated (0..nseq-1) by caller, and filled in here.
- */
-void
-PositionBasedWeights(char **aseq, int nseq, int alen, float *wgt)
-{
- int rescount[26]; /* count of A-Z residues in a column */
- int nres; /* number of different residues in col */
- int idx, pos; /* indices into aseq */
- int x;
- float norm;
-
- FSet(wgt, nseq, 0.0);
- for (pos = 0; pos < alen; pos++)
- {
- for (x = 0; x < 26; x++) rescount[x] = 0;
- for (idx = 0; idx < nseq; idx++)
- if (isalpha(aseq[idx][pos]))
- rescount[toupper(aseq[idx][pos]) - 'A'] ++;
-
- nres = 0;
- for (x = 0; x < 26; x++)
- if (rescount[x] > 0) nres++;
-
- for (idx = 0; idx < nseq; idx++)
- if (isalpha(aseq[idx][pos]))
- wgt[idx] += 1. / (float) (nres * rescount[toupper(aseq[idx][pos]) - 'A']);
- }
-
- for (idx = 0; idx < nseq; idx++)
- wgt[idx] /= (float) DealignedLength(aseq[idx]);
- norm = (float) nseq / FSum(wgt, nseq);
- FScale(wgt, nseq, norm);
- return;
-}
-
-
-
-
-/* Function: FilterAlignment()
- * Date: SRE, Wed Jun 30 09:19:30 1999 [St. Louis]
- *
- * Purpose: Constructs a new alignment by removing near-identical
- * sequences from a given alignment (where identity is
- * calculated *based on the alignment*).
- * Does not affect the given alignment.
- * Keeps earlier sequence, discards later one.
- *
- * Usually called as an ad hoc sequence "weighting" mechanism.
- *
- * Limitations:
- * Unparsed Stockholm markup is not propagated into the
- * new alignment.
- *
- * Args: msa -- original alignment
- * cutoff -- fraction identity cutoff. 0.8 removes sequences > 80% id.
- * ret_new -- RETURN: new MSA, usually w/ fewer sequences
- *
- * Return: (void)
- * ret_new must be free'd by caller: MSAFree().
- */
-void
-FilterAlignment(MSA *msa, float cutoff, MSA **ret_new)
-{
- int nnew; /* number of seqs in new alignment */
- int *list;
- int *useme;
- float ident;
- int i,j;
- int remove;
-
- /* find which seqs to keep (list) */
- /* diff matrix; allow ragged ends */
- list = MallocOrDie (sizeof(int) * msa->nseq);
- useme = MallocOrDie (sizeof(int) * msa->nseq);
- for (i = 0; i < msa->nseq; i++) useme[i] = FALSE;
-
- nnew = 0;
- for (i = 0; i < msa->nseq; i++)
- {
- remove = FALSE;
- for (j = 0; j < nnew; j++)
- {
- ident = PairwiseIdentity(msa->aseq[i], msa->aseq[list[j]]);
- if (ident > cutoff)
- {
- remove = TRUE;
- printf("removing %12s -- fractional identity %.2f to %s\n",
- msa->sqname[i], ident,
- msa->sqname[list[j]]);
- break;
- }
- }
- if (remove == FALSE) {
- list[nnew++] = i;
- useme[i] = TRUE;
- }
- }
-
- MSASmallerAlignment(msa, useme, ret_new);
- free(list);
- free(useme);
- return;
-}
-
-
-/* Function: SampleAlignment()
- * Date: SRE, Wed Jun 30 10:13:56 1999 [St. Louis]
- *
- * Purpose: Constructs a new, smaller alignment by sampling a given
- * number of sequences at random. Does not change the
- * alignment nor the order of the sequences.
- *
- * If you ask for a sample that is larger than nseqs,
- * it silently returns the original alignment.
- *
- * Not really a weighting method, but this is as good
- * a place as any to keep it, since it's similar in
- * construction to FilterAlignment().
- *
- * Args: msa -- original alignment
- * sample -- number of sequences in new alignment (0 < sample <= nseq)
- * ret_new -- RETURN: new MSA
- *
- * Return: (void)
- * ret_new must be free'd by caller: MSAFree().
- */
-void
-SampleAlignment(MSA *msa, int sample, MSA **ret_new)
-{
- int *list; /* array for random selection w/o replace */
- int *useme; /* array of flags 0..nseq-1: TRUE to use */
- int i, idx;
- int len;
-
- /* Allocations
- */
- list = (int *) MallocOrDie (sizeof(int) * msa->nseq);
- useme = (int *) MallocOrDie (sizeof(int) * msa->nseq);
- for (i = 0; i < msa->nseq; i++)
- {
- list[i] = i;
- useme[i] = FALSE;
- }
-
- /* Sanity check.
- */
- if (sample >= msa->nseq) sample = msa->nseq;
-
- /* random selection w/o replacement */
- for (len = msa->nseq, i = 0; i < sample; i++)
- {
- idx = CHOOSE(len);
- printf("chose %d: %s\n", list[idx], msa->sqname[list[idx]]);
- useme[list[idx]] = TRUE;
- list[idx] = list[--len];
- }
-
- MSASmallerAlignment(msa, useme, ret_new);
- free(list);
- free(useme);
- return;
-}
-
-
-/* Function: SingleLinkCluster()
- * Date: SRE, Fri Jul 16 15:02:57 1999 [St. Louis]
- *
- * Purpose: Perform simple single link clustering of seqs in a
- * sequence alignment. A pairwise identity threshold
- * defines whether two sequences are linked or not.
- *
- * Important: runs in O(N) memory, unlike standard
- * graph decomposition algorithms that use O(N^2)
- * adjacency matrices or adjacency lists. Requires
- * O(N^2) time in worst case (which is when you have
- * no links at all), O(NlogN) in "average"
- * case, and O(N) in best case (when there is just
- * one cluster in a completely connected graph.
- *
- * (Developed because hmmbuild could no longer deal
- * with GP120, a 16,013 sequence alignment.)
- *
- * Limitations:
- * CASE-SENSITIVE. Assumes aseq have been put into
- * either all lower or all upper case; or at least,
- * within a column, there's no mixed case.
- *
- * Algorithm:
- * I don't know if this algorithm is published. I
- * haven't seen it in graph theory books, but that might
- * be because it's so obvious that nobody's bothered.
- *
- * In brief, we're going to do a breadth-first search
- * of the graph, and we're going to calculate links
- * on the fly rather than precalculating them into
- * some sort of standard adjacency structure.
- *
- * While working, we keep two stacks of maximum length N:
- * a : list of vertices that are still unconnected.
- * b : list of vertices that we've connected to
- * in our current breadth level, but we haven't
- * yet tested for other connections to a.
- * The current length (number of elements in) a and b are
- * kept in na, nb.
- *
- * We store our results in an array of length N:
- * c : assigns each vertex to a component. for example
- * c[4] = 1 means that vertex 4 is in component 1.
- * nc is the number of components. Components
- * are numbered from 0 to nc-1. We return c and nc
- * to our caller.
- *
- * The algorithm is:
- *
- * Initialisation:
- * a <-- all the vertices
- * na <-- N
- * b <-- empty set
- * nb <-- 0
- * nc <-- 0
- *
- * Then:
- * while (a is not empty)
- * pop a vertex off a, push onto b
- * while (b is not empty)
- * pop vertex v off b
- * assign c[v] = nc
- * for each vertex w in a:
- * compare v,w. If w is linked to v, remove w
- * from a, push onto b.
- * nc++
- * q.e.d. :)
- *
- * Args: aseq - aligned sequences
- * nseq - number of sequences in aseq
- * alen - alignment length
- * maxid - fractional identity threshold 0..1. if id >= maxid, seqs linked
- * ret_c - RETURN: 0..nseq-1 assignments of seqs to components (clusters)
- * ret_nc - RETURN: number of components
- *
- * Returns: void.
- * ret_c is allocated here. Caller free's with free(*ret_c)
- */
-void
-SingleLinkCluster(char **aseq, int nseq, int alen, float maxid,
- int **ret_c, int *ret_nc)
-{
- int *a, na; /* stack of available vertices */
- int *b, nb; /* stack of working vertices */
- int *c; /* array of results */
- int nc; /* total number of components */
- int v,w; /* index of a working vertices */
- int i; /* loop counter */
-
- /* allocations and initializations
- */
- a = MallocOrDie (sizeof(int) * nseq);
- b = MallocOrDie (sizeof(int) * nseq);
- c = MallocOrDie (sizeof(int) * nseq);
- for (i = 0; i < nseq; i++) a[i] = i;
- na = nseq;
- nb = 0;
- nc = 0;
-
- /* Main algorithm
- */
- while (na > 0)
- {
- v = a[na-1]; na--; /* pop a vertex off a, */
- b[nb] = v; nb++; /* and push onto b */
- while (nb > 0)
- {
- v = b[nb-1]; nb--; /* pop vertex off b */
- c[v] = nc; /* assign it to component nc */
- for (i = na-1; i >= 0; i--)/* backwards, becase of deletion/swapping we do*/
- if (simple_distance(aseq[v], aseq[a[i]]) <= 1. - maxid) /* linked? */
- {
- w = a[i]; a[i] = a[na-1]; na--; /* delete w from a (note swap) */
- b[nb] = w; nb++; /* push w onto b */
- }
- }
- nc++;
- }
-
- /* Cleanup and return
- */
- free(a);
- free(b);
- *ret_c = c;
- *ret_nc = nc;
- return;
-}
git://source.jalview.org / jabaws.git / blobdiff