+++ /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.
- *****************************************************************/
-
-/* alignio.c
- * SRE, Mon Jul 12 11:57:37 1993
- * RCS $Id: alignio.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: alignio.c,v 1.11 2002/10/09 14:26:09 eddy Exp)
- *
- * Input/output of sequence alignments.
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include "squid.h"
-#include "sre_random.h"
-
-/* Function: AllocAlignment()
- *
- * Purpose: Allocate space for an alignment, given the number
- * of sequences and the alignment length in columns.
- *
- * Args: nseq - number of sequences
- * alen - width of alignment
- * ret_aseq - RETURN: alignment itself
- * ainfo - RETURN: other info associated with alignment
- *
- * Return: (void)
- * aseq, ainfo free'd by caller: FreeAlignment(aseq, &ainfo).
- * note that ainfo itself is alloc'ed in caller, usually
- * just by a "AINFO ainfo" definition.
- */
-void
-AllocAlignment(int nseq, int alen, char ***ret_aseq, AINFO *ainfo)
-{
- char **aseq;
- int idx;
-
- InitAinfo(ainfo);
-
- aseq = (char **) MallocOrDie (sizeof(char *) * nseq);
- for (idx = 0; idx < nseq; idx++)
- aseq[idx] = (char *) MallocOrDie (sizeof(char) * (alen+1));
-
- ainfo->alen = alen;
- ainfo->nseq = nseq;
-
- ainfo->wgt = (float *) MallocOrDie (sizeof(float) * nseq);
- FSet(ainfo->wgt, nseq, 1.0);
-
- ainfo->sqinfo = (SQINFO *) MallocOrDie (sizeof(SQINFO) * nseq);
- for (idx = 0; idx < nseq; idx++)
- ainfo->sqinfo[idx].flags = 0;
-
- *ret_aseq = aseq;
-}
-
-
-/* Function: InitAinfo()
- * Date: SRE, Tue Jan 19 10:16:02 1999 [St. Louis]
- *
- * Purpose: Initialize the fields in ainfo structure to
- * default (null) values. Does nothing with
- * fields that are dependent on nseq or alen.
- *
- * Args: ainfo - optional info structure for an alignment
- *
- * Returns: (void). ainfo is modified.
- */
-void
-InitAinfo(AINFO *ainfo)
-{
- ainfo->name = NULL;
- ainfo->desc = NULL;
- ainfo->cs = NULL;
- ainfo->rf = NULL;
- ainfo->acc = NULL;
- ainfo->au = NULL;
- ainfo->flags = 0;
-
- ainfo->tc1 = ainfo->tc2 = 0.0;
- ainfo->nc1 = ainfo->nc2 = 0.0;
- ainfo->ga1 = ainfo->ga2 = 0.0;
-}
-
-
-/* Function: FreeAlignment()
- *
- * Purpose: Free the space allocated to alignment, names, and optional
- * information.
- *
- * Args: aseqs - sequence alignment
- * ainfo - associated alignment data.
- */
-void
-FreeAlignment(char **aseqs, AINFO *ainfo)
-{
- int i;
-
- for (i = 0; i < ainfo->nseq; i++)
- {
- if (ainfo->sqinfo[i].flags & SQINFO_SS) free(ainfo->sqinfo[i].ss);
- if (ainfo->sqinfo[i].flags & SQINFO_SA) free(ainfo->sqinfo[i].sa);
- }
- if (ainfo->cs != NULL) free(ainfo->cs);
- if (ainfo->rf != NULL) free(ainfo->rf);
- if (ainfo->name != NULL) free(ainfo->name);
- if (ainfo->desc != NULL) free(ainfo->desc);
- if (ainfo->acc != NULL) free(ainfo->acc);
- if (ainfo->au != NULL) free(ainfo->au);
-
- free(ainfo->sqinfo);
- free(ainfo->wgt);
- Free2DArray((void **) aseqs, ainfo->nseq);
-}
-
-
-
-/* Function: SAMizeAlignment()
- * Date: SRE, Tue Jun 30 09:49:40 1998 [St. Louis]
- *
- * Purpose: Make a "best effort" attempt to convert an alignment
- * to SAM gap format: - in delete col, . in insert col.
- * Only works if alignment adheres to SAM's upper/lower
- * case convention, which is true for instance of old
- * HMMER alignments.
- *
- * Args: aseq - alignment to convert
- * nseq - number of seqs in alignment
- * alen - length of alignment
- *
- * Returns: (void)
- */
-void
-SAMizeAlignment(char **aseq, int nseq, int alen)
-{
- int col; /* counter for aligned columns */
- int i; /* counter for seqs */
- int sawlower, sawupper, sawgap;
- char gapchar;
-
- for (col = 0; col < alen; col++)
- {
- sawlower = sawupper = sawgap = 0;
- /* pass 1: do we see only upper or lower? */
- for (i = 0; i < nseq; i++)
- {
- if (isgap(aseq[i][col])) { sawgap = 1; continue; }
- if (isupper((int) aseq[i][col])) { sawupper = 1; continue; }
- if (islower((int) aseq[i][col])) sawlower = 1;
- }
- /* select gap character for column */
- gapchar = '-'; /* default */
- if (sawlower && ! sawupper) gapchar = '.';
-
- /* pass 2: set gap char */
- for (i = 0; i < nseq; i++)
- if (isgap(aseq[i][col])) aseq[i][col] = gapchar;
- }
-}
-
-
-/* Function: SAMizeAlignmentByGapFrac()
- * Date: SRE, Tue Jun 30 10:58:38 1998 [St. Louis]
- *
- * Purpose: Convert an alignment to SAM's gap and case
- * conventions, using gap fraction in a column
- * to choose match versus insert columns. In match columns,
- * residues are upper case and gaps are '-'.
- * In insert columns, residues are lower case and
- * gaps are '.'
- *
- * Args: aseq - aligned sequences
- * nseq - number of sequences
- * alen - length of alignment
- * maxgap - if more gaps than this fraction, column is insert.
- *
- * Returns: (void) Characters in aseq may be altered.
- */
-void
-SAMizeAlignmentByGapFrac(char **aseq, int nseq, int alen, float maxgap)
-{
- int apos; /* counter over columns */
- int idx; /* counter over sequences */
- int ngap; /* number of gaps seen */
-
- for (apos = 0; apos < alen; apos++)
- {
- /* count gaps */
- ngap = 0;
- for (idx = 0; idx < nseq; idx++)
- if (isgap(aseq[idx][apos])) ngap++;
-
- /* convert to SAM conventions */
- if ((float) ngap / (float) nseq > maxgap)
- { /* insert column */
- for (idx = 0; idx < nseq; idx++)
- if (isgap(aseq[idx][apos])) aseq[idx][apos] = '.';
- else aseq[idx][apos] = (char) tolower((int) aseq[idx][apos]);
- }
- else
- { /* match column */
- for (idx = 0; idx < nseq; idx++)
- if (isgap(aseq[idx][apos])) aseq[idx][apos] = '-';
- else aseq[idx][apos] = (char) toupper((int) aseq[idx][apos]);
- }
- }
-}
-
-
-
-
-/* Function: MakeAlignedString()
- *
- * Purpose: Given a raw string of some type (secondary structure, say),
- * align it to a given aseq by putting gaps wherever the
- * aseq has gaps.
- *
- * Args: aseq: template for alignment
- * alen: length of aseq
- * ss: raw string to align to aseq
- * ret_s: RETURN: aligned ss
- *
- * Return: 1 on success, 0 on failure (and squid_errno is set.)
- * ret_ss is malloc'ed here and must be free'd by caller.
- */
-int
-MakeAlignedString(char *aseq, int alen, char *ss, char **ret_s)
-{
- char *new;
- int apos, rpos;
-
- new = (char *) MallocOrDie ((alen+1) * sizeof(char));
- for (apos = rpos = 0; apos < alen; apos++)
- if (! isgap(aseq[apos]))
- {
- new[apos] = ss[rpos];
- rpos++;
- }
- else
- new[apos] = '.';
- new[apos] = '\0';
-
- if (rpos != strlen(ss))
- { squid_errno = SQERR_PARAMETER; free(new); return 0; }
- *ret_s = new;
- return 1;
-}
-
-
-/* Function: MakeDealignedString()
- *
- * Purpose: Given an aligned string of some type (either sequence or
- * secondary structure, for instance), dealign it relative
- * to a given aseq. Return a ptr to the new string.
- *
- * Args: aseq : template alignment
- * alen : length of aseq
- * ss: : string to make dealigned copy of; same length as aseq
- * ret_s : RETURN: dealigned copy of ss
- *
- * Return: 1 on success, 0 on failure (and squid_errno is set)
- * ret_s is alloc'ed here and must be freed by caller
- */
-int
-MakeDealignedString(char *aseq, int alen, char *ss, char **ret_s)
-{
- char *new;
- int apos, rpos;
-
- new = (char *) MallocOrDie ((alen+1) * sizeof(char));
- for (apos = rpos = 0; apos < alen; apos++)
- if (! isgap(aseq[apos]))
- {
- new[rpos] = ss[apos];
- rpos++;
- }
- new[rpos] = '\0';
- if (alen != strlen(ss))
- { squid_errno = SQERR_PARAMETER; free(new); return 0; }
- *ret_s = new;
- return 1;
-}
-
-
-/* Function: DealignedLength()
- *
- * Purpose: Count the number of non-gap symbols in seq.
- * (i.e. find the length of the unaligned sequence)
- *
- * Args: aseq - aligned sequence to count symbols in, \0 terminated
- *
- * Return: raw length of seq.
- */
-int
-DealignedLength(char *aseq)
-{
- int rlen;
- for (rlen = 0; *aseq; aseq++)
- if (! isgap(*aseq)) rlen++;
- return rlen;
-}
-
-
-/* Function: WritePairwiseAlignment()
- *
- * Purpose: Write a nice formatted pairwise alignment out,
- * with a BLAST-style middle line showing identities
- * as themselves (single letter) and conservative
- * changes as '+'.
- *
- * Args: ofp - open fp to write to (stdout, perhaps)
- * aseq1, aseq2 - alignments to write (not necessarily
- * flushed right with gaps)
- * name1, name2 - names of sequences
- * spos1, spos2 - starting position in each (raw) sequence
- * pam - PAM matrix; positive values define
- * conservative changes
- * indent - how many extra spaces to print on left
- *
- * Return: 1 on success, 0 on failure
- */
-int
-WritePairwiseAlignment(FILE *ofp,
- char *aseq1, char *name1, int spos1,
- char *aseq2, char *name2, int spos2,
- int **pam, int indent)
-{
- char sname1[11]; /* shortened name */
- char sname2[11];
- int still_going; /* True if writing another block */
- char buf1[61]; /* buffer for writing seq1; CPL+1*/
- char bufmid[61]; /* buffer for writing consensus */
- char buf2[61];
- char *s1, *s2; /* ptrs into each sequence */
- int count1, count2; /* number of symbols we're writing */
- int rpos1, rpos2; /* position in raw seqs */
- int rawcount1, rawcount2; /* number of nongap symbols written */
- int apos;
-
- strncpy(sname1, name1, 10);
- sname1[10] = '\0';
- strtok(sname1, WHITESPACE);
-
- strncpy(sname2, name2, 10);
- sname2[10] = '\0';
- strtok(sname2, WHITESPACE);
-
- s1 = aseq1;
- s2 = aseq2;
- rpos1 = spos1;
- rpos2 = spos2;
-
- still_going = TRUE;
- while (still_going)
- {
- still_going = FALSE;
-
- /* get next line's worth from both */
- strncpy(buf1, s1, 60); buf1[60] = '\0';
- strncpy(buf2, s2, 60); buf2[60] = '\0';
- count1 = strlen(buf1);
- count2 = strlen(buf2);
-
- /* is there still more to go? */
- if ((count1 == 60 && s1[60] != '\0') ||
- (count2 == 60 && s2[60] != '\0'))
- still_going = TRUE;
-
- /* shift seq ptrs by a line */
- s1 += count1;
- s2 += count2;
-
- /* assemble the consensus line */
- for (apos = 0; apos < count1 && apos < count2; apos++)
- {
- if (!isgap(buf1[apos]) && !isgap(buf2[apos]))
- {
- if (buf1[apos] == buf2[apos])
- bufmid[apos] = buf1[apos];
- else if (pam[buf1[apos] - 'A'][buf2[apos] - 'A'] > 0)
- bufmid[apos] = '+';
- else
- bufmid[apos] = ' ';
- }
- else
- bufmid[apos] = ' ';
- }
- bufmid[apos] = '\0';
-
- rawcount1 = 0;
- for (apos = 0; apos < count1; apos++)
- if (!isgap(buf1[apos])) rawcount1++;
-
- rawcount2 = 0;
- for (apos = 0; apos < count2; apos++)
- if (!isgap(buf2[apos])) rawcount2++;
-
- (void) fprintf(ofp, "%*s%-10.10s %5d %s %5d\n", indent, "",
- sname1, rpos1, buf1, rpos1 + rawcount1 -1);
- (void) fprintf(ofp, "%*s %s\n", indent, "",
- bufmid);
- (void) fprintf(ofp, "%*s%-10.10s %5d %s %5d\n", indent, "",
- sname2, rpos2, buf2, rpos2 + rawcount2 -1);
- (void) fprintf(ofp, "\n");
-
- rpos1 += rawcount1;
- rpos2 += rawcount2;
- }
-
- return 1;
-}
-
-
-/* Function: MingapAlignment()
- *
- * Purpose: Remove all-gap columns from a multiple sequence alignment
- * and its associated data. The alignment is assumed to be
- * flushed (all aseqs the same length).
- */
-int
-MingapAlignment(char **aseqs, AINFO *ainfo)
-{
- int apos; /* position in original alignment */
- int mpos; /* position in new alignment */
- int idx;
-
- /* We overwrite aseqs, using its allocated memory.
- */
- for (apos = 0, mpos = 0; aseqs[0][apos] != '\0'; apos++)
- {
- /* check for all-gap in column */
- for (idx = 0; idx < ainfo->nseq; idx++)
- if (! isgap(aseqs[idx][apos]))
- break;
- if (idx == ainfo->nseq) continue;
-
- /* shift alignment and ainfo */
- if (mpos != apos)
- {
- for (idx = 0; idx < ainfo->nseq; idx++)
- aseqs[idx][mpos] = aseqs[idx][apos];
-
- if (ainfo->cs != NULL) ainfo->cs[mpos] = ainfo->cs[apos];
- if (ainfo->rf != NULL) ainfo->rf[mpos] = ainfo->rf[apos];
- }
- mpos++;
- }
- /* null terminate everything */
- for (idx = 0; idx < ainfo->nseq; idx++)
- aseqs[idx][mpos] = '\0';
- ainfo->alen = mpos; /* set new length */
- if (ainfo->cs != NULL) ainfo->cs[mpos] = '\0';
- if (ainfo->rf != NULL) ainfo->rf[mpos] = '\0';
- return 1;
-}
-
-
-
-/* Function: RandomAlignment()
- *
- * Purpose: Create a random alignment from raw sequences.
- *
- * Ideally, we would like to sample an alignment from the
- * space of possible alignments according to its probability,
- * given a prior probability distribution for alignments.
- * I don't see how to describe such a distribution, let alone
- * sample it.
- *
- * This is a rough approximation that tries to capture some
- * desired properties. We assume the alignment is generated
- * by a simple HMM composed of match and insert states.
- * Given parameters (pop, pex) for the probability of opening
- * and extending an insertion, we can find the expected number
- * of match states, M, in the underlying model for each sequence.
- * We use an average M taken over all the sequences (this is
- * an approximation. The expectation of M given all the sequence
- * lengths is a nasty-looking summation.)
- *
- * M = len / ( 1 + pop ( 1 + 1/ (1-pex) ) )
- *
- * Then, we assign positions in each raw sequence onto the M match
- * states and M+1 insert states of this "HMM", by rolling random
- * numbers and inserting the (rlen-M) inserted positions randomly
- * into the insert slots, taking into account the relative probability
- * of open vs. extend.
- *
- * The resulting alignment has two desired properties: insertions
- * tend to follow the HMM-like exponential distribution, and
- * the "sparseness" of the alignment is controllable through
- * pop and pex.
- *
- * Args: rseqs - raw sequences to "align", 0..nseq-1
- * sqinfo - array of 0..nseq-1 info structures for the sequences
- * nseq - number of sequences
- * pop - probability to open insertion (0<pop<1)
- * pex - probability to extend insertion (0<pex<1)
- * ret_aseqs - RETURN: alignment (flushed)
- * ainfo - fill in: alignment info
- *
- * Return: 1 on success, 0 on failure. Sets squid_errno to indicate cause
- * of failure.
- */
-int
-RandomAlignment(char **rseqs, SQINFO *sqinfo, int nseq, float pop, float pex,
- char ***ret_aseqs, AINFO *ainfo)
-{
- char **aseqs; /* RETURN: alignment */
- int alen; /* length of alignment */
- int *rlen; /* lengths of each raw sequence */
- int M; /* length of "model" */
- int **ins; /* insertion counts, 0..nseq-1 by 0..M */
- int *master_ins; /* max insertion counts, 0..M */
- int apos, rpos, idx;
- int statepos;
- int count;
- int minlen;
-
- /* calculate expected length of model, M
- */
- rlen = (int *) MallocOrDie (sizeof(int) * nseq);
- M = 0;
- minlen = 9999999;
- for (idx = 0; idx < nseq; idx++)
- {
- rlen[idx] = strlen(rseqs[idx]);
- M += rlen[idx];
- minlen = (rlen[idx] < minlen) ? rlen[idx] : minlen;
- }
- M = (int) ((float) M / (1.0 + pop * (1.0 + 1.0 / (1.0 - pex))));
- M /= nseq;
- if (M > minlen) M = minlen;
-
- /* make arrays that count insertions in M+1 possible insert states
- */
- ins = (int **) MallocOrDie (sizeof(int *) * nseq);
- master_ins = (int *) MallocOrDie (sizeof(int) * (M+1));
- for (idx = 0; idx < nseq; idx++)
- {
- ins[idx] = (int *) MallocOrDie (sizeof(int) * (M+1));
- for (rpos = 0; rpos <= M; rpos++)
- ins[idx][rpos] = 0;
- }
- /* normalize */
- pop = pop / (pop+pex);
- pex = 1.0 - pop;
- /* make insertions for individual sequences */
- for (idx = 0; idx < nseq; idx++)
- {
- apos = -1;
- for (rpos = 0; rpos < rlen[idx]-M; rpos++)
- {
- if (sre_random() < pop || apos == -1) /* open insertion */
- apos = CHOOSE(M+1); /* choose 0..M */
- ins[idx][apos]++;
- }
- }
- /* calculate master_ins, max inserts */
- alen = M;
- for (apos = 0; apos <= M; apos++)
- {
- master_ins[apos] = 0;
- for (idx = 0; idx < nseq; idx++)
- if (ins[idx][apos] > master_ins[apos])
- master_ins[apos] = ins[idx][apos];
- alen += master_ins[apos];
- }
-
-
- /* Now, construct alignment
- */
- aseqs = (char **) MallocOrDie (sizeof (char *) * nseq);
- for (idx = 0; idx < nseq; idx++)
- aseqs[idx] = (char *) MallocOrDie (sizeof(char) * (alen+1));
- for (idx = 0; idx < nseq; idx++)
- {
- apos = rpos = 0;
-
- for (statepos = 0; statepos <= M; statepos++)
- {
- for (count = 0; count < ins[idx][statepos]; count++)
- aseqs[idx][apos++] = rseqs[idx][rpos++];
- for (; count < master_ins[statepos]; count++)
- aseqs[idx][apos++] = ' ';
-
- if (statepos != M)
- aseqs[idx][apos++] = rseqs[idx][rpos++];
- }
- aseqs[idx][alen] = '\0';
- }
- ainfo->flags = 0;
- ainfo->alen = alen;
- ainfo->nseq = nseq;
- ainfo->sqinfo = (SQINFO *) MallocOrDie (sizeof(SQINFO) * nseq);
- for (idx = 0; idx < nseq; idx++)
- SeqinfoCopy(&(ainfo->sqinfo[idx]), &(sqinfo[idx]));
-
- free(rlen);
- free(master_ins);
- Free2DArray((void **) ins, nseq);
- *ret_aseqs = aseqs;
- return 1;
-}
-
-/* Function: AlignmentHomogenousGapsym()
- * Date: SRE, Sun Mar 19 19:37:12 2000 [wren, St. Louis]
- *
- * Purpose: Sometimes we've got to convert alignments to
- * a lowest common denominator, and we need
- * a single specific gap character -- for example,
- * PSI-BLAST blastpgp -B takes a very simplistic
- * alignment input format which appears to only
- * allow '-' as a gap symbol.
- *
- * Anything matching the isgap() macro is
- * converted.
- *
- * Args: aseq - aligned character strings, [0..nseq-1][0..alen-1]
- * nseq - number of aligned strings
- * alen - length of alignment
- * gapsym - character to use for gaps.
- *
- * Returns: void ("never fails")
- */
-void
-AlignmentHomogenousGapsym(char **aseq, int nseq, int alen, char gapsym)
-{
- int i, apos;
-
- for (i = 0; i < nseq; i++)
- for (apos = 0; apos < alen; apos++)
- if (isgap(aseq[i][apos])) aseq[i][apos] = gapsym;
-}
git://source.jalview.org / jabaws.git / blobdiff