+++ /dev/null
-
-/* copyright (c) 1998, 1999 William R. Pearson and the U. of Virginia */
-
-/* $Name: fa_34_26_5 $ - $Id: dropfx.c,v 1.68 2007/04/26 18:37:18 wrp Exp $ */
-
-/* implements the fastx algorithm, see:
-
- W. R. Pearson, T. Wood, Z. Zhang, A W. Miller (1997) "Comparison of
- DNA sequences with protein sequences" Genomics 46:24-36
-
- see dropnfa.c for better variable descriptions and comments
-*/
-
-/* 18-Sept-2006 - remove global variables used for alignment */
-
-/* 22-June-2006 - correct incorrect alignment coordinates generated
- after pro_dna() on projected DNA region.
-*/
-
-/* 9-May-2003 -> 3.46 changed lx_band to use projected protein
- boundary end. this fixes some addressing issues on MacOSX, and
- speeds up alignment on very long proteins
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-#include <math.h>
-
-#include "defs.h"
-#include "param.h"
-#define XTERNAL
-#include "upam.h"
-
-/* this must be consistent with upam.h */
-#define MAXHASH 32
-#define NMAP MAXHASH+1
-
-/* globals for fasta */
-#define MAXWINDOW 64
-
-#ifndef MAXSAV
-#define MAXSAV 10
-#endif
-
-#ifndef ALLOCN0
-static char *verstr="3.5 Sept 2006";
-#else
-static char *verstr="3.5an0 May 2006";
-#endif
-
-struct dstruct /* diagonal structure for saving current run */
-{
- int score; /* hash score of current match */
- int start; /* start of current match */
- int stop; /* end of current match */
- struct savestr *dmax; /* location in vmax[] where best score data saved */
-};
-
-struct savestr
-{
- int score; /* pam score with segment optimization */
- int score0; /* pam score of best single segment */
- int gscore; /* score from global match */
- int dp; /* diagonal of match */
- int start; /* start of match in lib seq */
- int stop; /* end of match in lib seq */
-};
-
-struct swstr { int H, E;};
-
-struct bdstr { int CC, DD, CP, DP;};
-
-void savemax();
-void kpsort();
-
-struct sx_s {int C1, C2, C3, I1, I2, I3, flag; };
-
-struct f_struct {
- struct dstruct *diag;
- struct savestr vmax[MAXSAV]; /* best matches saved for one sequence */
- struct savestr *vptr[MAXSAV];
- struct savestr *lowmax;
- int ndo;
- int noff;
- int hmask; /* hash constants */
- int *pamh1; /* pam based array */
- int *pamh2; /* pam based kfact array */
- int *link, *harr; /* hash arrays */
- int kshft; /* shift width */
- int nsav, lowscor; /* number of saved runs, worst saved run */
-#ifndef TFAST
- unsigned char *aa0x; /* contains translated codons 111222333*/
- unsigned char *aa0y; /* contains translated codons 123123123*/
-#else
- unsigned char *aa1x; /* contains translated codons 111222333 */
- unsigned char *aa1y; /* contains translated codons 123123123 */
-#endif
- struct sx_s *cur;
- int *waa0;
- int *waa1;
- int *res;
- int max_res;
-};
-
-#define DROP_INTERN
-#include "drop_func.h"
-
-static int dmatchx(const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int hoff, int window,
- int **pam2, int gdelval, int ggapval, int gshift,
- struct f_struct *f_str);
-
-int shscore(unsigned char *aa0, int n0, int **pam2);
-int saatran(const unsigned char *ntseq, unsigned char *aaseq, int maxs, int frame);
-int spam (const unsigned char *aa0, const unsigned char *aa1,
- struct savestr *dmax, int **pam2,
- struct f_struct *f_str);
-int sconn (struct savestr **v, int n,int cgap, int pgap, struct f_struct *f_str);
-int lx_band(const unsigned char *prot_seq, int len_prot,
- const unsigned char *dna_prot_seq, int len_dna_prot,
- int **pam_matrix, int gopen, int gext,
- int gshift, int start_diag, int width, struct f_struct *f_str);
-
-static void
-update_code(char *al_str, int al_str_max, int op, int op_cnt, char *op_char);
-
-extern void w_abort (char *p, char *p1);
-
-/* initialize for fasta */
-
-void
-init_work (unsigned char *aa0, int n0,
- struct pstruct *ppst,
- struct f_struct **f_arg)
-{
- int mhv, phv;
- int hmax;
- int i0, hv;
- int pamfact;
- int btemp;
- struct f_struct *f_str;
- int ktup; /* word size examined */
- int fact; /* factor used to scale ktup match value */
- int kt1; /* ktup-1 */
- int lkt; /* last ktup - initiall kt1, but can be increased
- for hsq >= NMAP */
-
- int maxn0;
- int *pwaa;
- int i, j, q;
- struct swstr *ss, *r_ss;
- int *waa;
- int *res;
- int nsq, ip, *hsq;
-#ifndef TFAST
- int last_n0, itemp;
- unsigned char *fd, *fs, *aa0x, *aa0y, *aa0s;
- int n0x, n0x3;
-#endif
-
- if (ppst->ext_sq_set) {
- nsq = ppst->nsqx; ip = 1;
- hsq = ppst->hsqx;
- }
- else {
- nsq = ppst->nsq; ip = 0;
- hsq = ppst->hsq;
- }
-
- f_str = (struct f_struct *)calloc(1,sizeof(struct f_struct));
-
- btemp = 2 * ppst->param_u.fa.bestoff / 3 +
- n0 / ppst->param_u.fa.bestscale +
- ppst->param_u.fa.bkfact *
- (ppst->param_u.fa.bktup - ppst->param_u.fa.ktup);
- btemp = min (btemp, ppst->param_u.fa.bestmax);
- if (btemp > 3 * n0) btemp = 3 * shscore(aa0,n0,ppst->pam2[0]) / 5;
-
- ppst->param_u.fa.cgap = btemp + ppst->param_u.fa.bestoff / 3;
- if (ppst->param_u.fa.optcut_set != 1)
-#ifndef TFAST
- ppst->param_u.fa.optcut = (btemp*5)/4;
-#else
- ppst->param_u.fa.optcut = (btemp*4)/3;
-#endif
-
-#ifdef OLD_FASTA_GAP
- ppst->param_u.fa.pgap = ppst->gdelval + ppst->ggapval;
-#else
- ppst->param_u.fa.pgap = ppst->gdelval + 2*ppst->ggapval;
-#endif
- pamfact = ppst->param_u.fa.pamfact;
- ktup = ppst->param_u.fa.ktup;
- fact = ppst->param_u.fa.scfact * ktup;
-
- if (pamfact == -1)
- pamfact = 0;
- else if (pamfact == -2)
- pamfact = 1;
-
- for (i0 = 1, mhv = -1; i0 <=nsq; i0++)
- if (hsq[i0] < NMAP && hsq[i0] > mhv) mhv = hsq[i0];
-
- if (mhv <= 0) {
- fprintf (stderr, " maximum hsq <=0 %d\n", mhv);
- exit (1);
- }
-
- for (f_str->kshft = 0; mhv > 0; mhv /= 2)
- f_str->kshft++;
-
-/* kshft = 2; */
- kt1 = ktup - 1;
- hv = 1;
- for (i0 = 0; i0 < ktup; i0++) {
- hv = hv << f_str->kshft;
- }
- hmax = hv;
- f_str->hmask = (hmax >> f_str->kshft) - 1;
-
-
- if ((f_str->harr = (int *) calloc (hmax, sizeof (int))) == NULL) {
- fprintf (stderr, " cannot allocate hash array\n");
- exit (1);
- }
- if ((f_str->pamh1 = (int *) calloc (nsq+1, sizeof (int))) == NULL) {
- fprintf (stderr, " cannot allocate pamh1 array\n");
- exit (1);
- }
- if ((f_str->pamh2 = (int *) calloc (hmax, sizeof (int))) == NULL) {
- fprintf (stderr, " cannot allocate pamh2 array\n");
- exit (1);
- }
- if ((f_str->link = (int *) calloc (n0, sizeof (int))) == NULL) {
- fprintf (stderr, " cannot allocate hash link array");
- exit (1);
- }
-
-#ifdef TFAST
- if ((f_str->aa1x =(unsigned char *)calloc((size_t)ppst->maxlen+2,
- sizeof(unsigned char)))
- == NULL) {
- fprintf (stderr, "cannot allocate aa1x array %d\n", ppst->maxlen+2);
- exit (1);
- }
- f_str->aa1x++;
-
- if ((f_str->aa1y =(unsigned char *)calloc((size_t)ppst->maxlen+2,
- sizeof(unsigned char)))
- == NULL) {
- fprintf (stderr, "cannot allocate aa1y array %d\n", ppst->maxlen+2);
- exit (1);
- }
- f_str->aa1y++;
-#else /* FASTX */
- maxn0 = n0 + 2;
- if ((aa0x =(unsigned char *)calloc((size_t)maxn0,sizeof(unsigned char)))
- == NULL) {
- fprintf (stderr, "cannot allocate aa0x array %d\n", maxn0);
- exit (1);
- }
- aa0x++;
- f_str->aa0x = aa0x;
-
- if ((aa0y =(unsigned char *)calloc((size_t)maxn0,sizeof(unsigned char)))
- == NULL) {
- fprintf (stderr, "cannot allocate aa0y array %d\n", maxn0);
- exit (1);
- }
- aa0y++;
- f_str->aa0y = aa0y;
-
- last_n0 = 0;
- for (itemp=0; itemp<3; itemp++) {
- n0x = saatran(aa0,&aa0x[last_n0],n0,itemp);
- /*
- for (i=0; i<n0x; i++) {
- fprintf(stderr,"%c",aa[aa0x[last_n0+i]]);
- if ((i%60)==59) fprintf(stderr,"\n");
- }
- fprintf(stderr,"\n");
- */
- last_n0 += n0x+1;
- }
-
- /* fprintf(stderr,"\n"); */
-
- for (itemp=0, fs=aa0x; itemp <3; itemp++,fs++) {
- for (fd = &aa0y[itemp]; *fs!=EOSEQ; fd += 3, fs++) *fd = *fs;
- *fd=EOSEQ;
- }
-
- /* now switch aa0 and aa0x for hashing functions */
- /* this seems dangerous in threaded code, but only the pointer is changed,
- not the data itself */
-
- fs = aa0;
- aa0 = aa0x;
- aa0x = fs;
-
-#endif
-
- for (i0 = 0; i0 < hmax; i0++)
- f_str->harr[i0] = -1;
- for (i0 = 0; i0 < n0; i0++)
- f_str->link[i0] = -1;
-
- /* encode the aa0 array */
-
- phv = hv = 0;
- lkt = kt1;
- for (i0 = 0; i0 < min(lkt,n0); i0++) {
- if (hsq[aa0[i0]] >= NMAP) {hv=phv=0; lkt=i0+ktup; continue;}
- hv = (hv << f_str->kshft) + hsq[aa0[i0]];
- phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup;
- }
-
- for (; i0 < n0; i0++) {
- if (hsq[aa0[i0]] >= NMAP) {
- hv=phv=0;
- lkt = i0+ktup;
- /* restart hv, phv calculation */
- for (; (i0 < lkt || hsq[aa0[i0]]>=NMAP) && i0<n0; i0++) {
- if (hsq[aa0[i0]] >= NMAP) {hv=phv=0; lkt = i0+ktup; continue;}
- hv = (hv << f_str->kshft) + hsq[aa0[i0]];
- phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup;
- }
- }
- if (i0 >= n0) break;
- hv = ((hv & f_str->hmask) << f_str->kshft) + hsq[aa0[i0]];
- f_str->link[i0] = f_str->harr[hv];
- f_str->harr[hv] = i0;
- if (pamfact) {
- f_str->pamh2[hv] = (phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup);
- /* this check should always be true, but just in case */
- if (hsq[aa0[i0-kt1]]<NMAP)
- phv -= ppst->pam2[ip][aa0[i0 - kt1]][aa0[i0 - kt1]] * ktup;
- }
- else f_str->pamh2[hv] = fact * ktup;
- }
-
-#ifndef TFAST
- /* done hashing, now switch aa0, aa0x back */
- fs = aa0;
- aa0 = aa0x;
- aa0x = fs;
-#endif
-
-/* this has been modified from 0..<nsq to 1..<=nsq because the
- pam2[0][0] is now undefined for consistency with blast
-*/
-
- if (pamfact)
- for (i0 = 1; i0 <= nsq; i0++)
- f_str->pamh1[i0] = ppst->pam2[ip][i0][i0] * ktup;
- else
- for (i0 = 1; i0 <= nsq; i0++)
- f_str->pamh1[i0] = fact;
-
- f_str->ndo = 0; /* used to save time on diagonals with long queries */
-
-#ifndef ALLOCN0
- if ((f_str->diag = (struct dstruct *) calloc ((size_t)MAXDIAG,
- sizeof (struct dstruct)))==NULL) {
- fprintf (stderr," cannot allocate diagonal arrays: %ld\n",
- (long) MAXDIAG *sizeof (struct dstruct));
- exit (1);
- };
-#else
- if ((f_str->diag = (struct dstruct *) calloc ((size_t)n0,
- sizeof (struct dstruct)))==NULL) {
- fprintf (stderr," cannot allocate diagonal arrays: %ld\n",
- (long)n0*sizeof (struct dstruct));
- exit (1);
- };
-#endif
-
-
- if ((waa= (int *)malloc (sizeof(int)*(nsq+1)*n0)) == NULL) {
- fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
- exit(1);
- }
-
- pwaa = waa;
- for (i=0; i<=nsq; i++) {
- for (j=0;j<n0; j++) {
- *pwaa = ppst->pam2[ip][i][aa0[j]];
- pwaa++;
- }
- }
- f_str->waa0 = waa;
-
- if ((waa= (int *)malloc (sizeof(int)*(nsq+1)*n0)) == NULL) {
- fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
- exit(1);
- }
-
- pwaa = waa;
- for (i=0; i<=nsq; i++) {
- for (j=0;j<n0; j++) {
- *pwaa = ppst->pam2[0][i][aa0[j]];
- pwaa++;
- }
- }
- f_str->waa1 = waa;
-
-#ifndef TFAST
- maxn0 = max(2*n0,MIN_RES);
-#else
- /* maxn0 needs to be large enough to accomodate introns
- for TFASTX. For all other functions, it will be
- more reasonable. */
- maxn0 = max(4*n0,MIN_RES);
-#endif
- if ((res = (int *)calloc((size_t)maxn0,sizeof(int)))==NULL) {
- fprintf(stderr,"cannot allocate alignment results array %d\n",maxn0);
- exit(1);
- }
- f_str->res = res;
- f_str->max_res = maxn0;
-
- *f_arg = f_str;
-}
-
-
-/* pstring1 is a message to the manager, currently 512 */
-/* pstring2 is the same information, but in a markx==10 format */
-void
-get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
-{
-#ifndef TFAST
- char *pg_str="FASTX";
-#else
- char *pg_str="TFASTX";
-#endif
-
- if (!pstr->param_u.fa.optflag)
-#ifdef OLD_FASTA_GAP
- sprintf (pstring1, "%s (%s) function [%s matrix (%d:%d:%d)%s] ktup: %d\n join: %d, gap-pen: %d/%d, shift: %d width: %3d",pg_str,verstr,
-#else
- sprintf (pstring1, "%s (%s) function [%s matrix (o=%d:%d:%d:%d)%s] ktup: %d\n join: %d, open/ext: %d/%d, shift: %d width: %3d",pg_str,verstr,
-#endif
- pstr->pamfile, pstr->pam_h,pstr->pam_l,pstr->pam_xx,pstr->pam_xm,
- (pstr->ext_sq_set) ? "xS":"\0",
- pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
- pstr->gdelval, pstr->ggapval, pstr->gshift,
- pstr->param_u.fa.optwid);
- else
-#ifdef OLD_FASTA_GAP
- sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d:%d)%s] ktup: %d\n join: %d, opt: %d, gap-pen: %d/%d shift: %3d, width: %3d",pg_str,verstr,
-#else
- sprintf (pstring1, "%s (%s) function [optimized, %s matrix (o=%d:%d:%d:%d)%s] ktup: %d\n join: %d, opt: %d, open/ext: %d/%d shift: %3d, width: %3d",pg_str,verstr,
-#endif
- pstr->pamfile, pstr->pam_h,pstr->pam_l,pstr->pam_xx, pstr->pam_xm,
- (pstr->ext_sq_set) ? "xS":"\0",
- pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
- pstr->param_u.fa.optcut, pstr->gdelval, pstr->ggapval,
- pstr->gshift,pstr->param_u.fa.optwid);
-
- if (pstr->param_u.fa.iniflag) strcat(pstring1," init1");
- /*
- if (pstr->zsflag==0) strcat(pstring1," not-scaled");
- else if (pstr->zsflag==1) strcat(pstring1," reg.-scaled");
- */
-
- if (pstring2 != NULL) {
-#ifdef OLD_FASTA_GAP
- sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)%s\n\
-; pg_gap-pen: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
-#else
- sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)%s\n\
-; pg_open_ext: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
-#endif
- pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
- (pstr->ext_sq_set) ? "xS":"\0", pstr->gdelval,
- pstr->ggapval,pstr->param_u.fa.ktup,pstr->param_u.fa.optcut,
- pstr->param_u.fa.cgap);
- }
-}
-
-void
-close_work (const unsigned char *aa0, int n0,
- struct pstruct *ppst,
- struct f_struct **f_arg)
-{
- struct f_struct *f_str;
-
- f_str = *f_arg;
-
- if (f_str != NULL) {
- free(f_str->cur);
-#ifndef TFAST
- f_str->aa0y--;
- free(f_str->aa0y);
- f_str->aa0x--;
- free(f_str->aa0x);
-#else
- f_str->aa1y--;
- free(f_str->aa1y);
- f_str->aa1x--;
- free(f_str->aa1x);
-#endif
- free(f_str->res);
- free(f_str->waa1);
- free(f_str->waa0);
- free(f_str->diag);
- free(f_str->link);
- free(f_str->pamh2);
- free(f_str->pamh1);
- free(f_str->harr);
- free(f_str);
- *f_arg = NULL;
- }
-}
-
-void do_fastx (const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- struct pstruct *ppst, struct f_struct *f_str,
- struct rstruct *rst, int *hoff)
-{
- int nd; /* diagonal array size */
- int lhval;
- int kfact;
- int i;
- int my_hoff;
- register struct dstruct *dptr;
- register int tscor;
-
-#ifndef ALLOCN0
- register struct dstruct *diagp;
-#else
- register int dpos;
- int lposn0;
-#endif
- struct dstruct *dpmax;
- register int lpos;
- int tpos;
- struct savestr *vmptr;
- int scor, tmp;
- int im, ib, nsave;
- int ktup, kt1, *hsq, ip, lkt;
-#ifndef TFAST
- int n0x31, n0x32;
- n0x31 = (n0-2)/3;
- n0x32 = n0x31+1+(n0-n0x31-1)/2;
-#else
- const unsigned char *fs;
- unsigned char *fd;
- int n1x31, n1x32, last_n1, itemp;
- n1x31 = (n1-2)/3;
- n1x32 = n1x31+1+(n1-n1x31-1)/2;
-#endif
-
- if (ppst->ext_sq_set) {
- ip = 1;
- hsq = ppst->hsqx;
- }
- else {
- ip = 0;
- hsq = ppst->hsq;
- }
-
- ktup = ppst->param_u.fa.ktup;
- kt1 = ktup-1;
-
- if (n1 < ktup) {
- rst->score[0] = rst->score[1] = rst->score[2] = 0;
- return;
- }
-
- if (n0+n1+1 >= MAXDIAG) {
- fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
- rst->score[0] = rst->score[1] = rst->score[2] = -1;
- return;
- }
-
- f_str->noff = n0 - 1;
-
-#ifdef ALLOCN0
- nd = n0;
-#endif
-
-#ifndef ALLOCN0
- nd = n0 + n1;
-#endif
-
- dpmax = &f_str->diag[nd];
- for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;)
- {
- dptr->stop = -1;
- dptr->dmax = NULL;
- dptr++->score = 0;
- }
-
- for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
- vmptr->score = 0;
- f_str->lowmax = f_str->vmax;
- f_str->lowscor = 0;
-
- /* start hashing */
- lhval = 0;
- lkt = kt1;
- for (lpos = 0; (lpos < lkt || hsq[aa1[lpos]]>=NMAP) && lpos<n1; lpos++) {
- if (hsq[aa1[lpos]]>=NMAP) {
- lhval = 0; lkt=lpos+ktup; continue;
-#ifdef ALLOCN0 /* reinitialize dptr */
- dptr = &f_str->diag[lpos % nd];
- dptr->stop = -1;
- dptr->dmax = NULL;
- dptr->score = 0;
-#endif
- }
- lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
- }
-
-#ifndef ALLOCN0
- diagp = &f_str->diag[f_str->noff + lkt];
- for (; lpos < n1; lpos++, diagp++) {
- /* if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; continue;} */
- if (hsq[aa1[lpos]]>=NMAP) {
- lpos++ ; diagp++;
- while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
- if (lpos >= n1) break;
- lhval = 0;
- }
- lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
- for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
- if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
-#else
- lposn0 = f_str->noff + lpos;
- for (; lpos < n1; lpos++, lposn0++) {
- if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; goto loopl;}
- lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
- for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
- dpos = lposn0 - tpos;
- if ((tscor = (dptr = &f_str->diag[dpos % nd])->stop) >= 0) {
-#endif
- tscor += ktup;
- if ((tscor -= lpos) <= 0) { /* better to start over */
- scor = dptr->score;
- if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 && f_str->lowscor < scor)
-#ifdef ALLOCN0
- savemax (dptr, dpos, f_str);
-#else
- savemax (dptr, f_str);
-#endif
- if ((tscor += scor) >= kfact) {
- dptr->score = tscor;
- dptr->stop = lpos;
- }
- else {
- dptr->score = kfact;
- dptr->start = (dptr->stop = lpos) - kt1;
- }
- } /* continue current run in diagonal */
- else {
- dptr->score += f_str->pamh1[aa0[tpos]];
- dptr->stop = lpos;
- }
- }
- else {
- dptr->score = f_str->pamh2[lhval];
- dptr->start = (dptr->stop = lpos) - kt1;
- }
- } /* end tpos */
-
-#ifdef ALLOCN0
- /* reinitialize diag structure */
- loopl:
- if ((dptr = &f_str->diag[lpos % nd])->score > f_str->lowscor) {
- savemax (dptr, lpos, f_str);
- }
- dptr->stop = -1;
- dptr->dmax = NULL;
- dptr->score = 0;
-#endif
- } /* end lpos */
-
-#ifdef ALLOCN0
- for (tpos = 0, dpos = f_str->noff + n1 - 1; tpos < n0; tpos++, dpos--) {
- if ((dptr = &f_str->diag[dpos % nd])->score > f_str->lowscor)
- savemax (dptr, dpos, f_str);
- }
-#else
- for (dptr = f_str->diag; dptr < dpmax;) {
- if (dptr->score > f_str->lowscor) savemax (dptr, f_str);
- dptr->stop = -1;
- dptr->dmax = NULL;
- dptr++->score = 0;
- }
- f_str->ndo = nd;
-#endif
-
-/*
- at this point all of the elements of aa1[lpos]
- have been searched for elements of aa0[tpos]
- with the results in diag[dpos]
-*/
-
- for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
- {
- /*
- fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
- f_str->noff+vmptr->start-vmptr->dp,
- f_str->noff+vmptr->stop-vmptr->dp,
- vmptr->start,vmptr->stop,
- vmptr->dp,vmptr->score);
- */
- if (vmptr->score > 0) {
- vmptr->score = spam (aa0, aa1, vmptr, ppst->pam2[ip], f_str);
- f_str->vptr[nsave++] = vmptr;
- }
- }
-
- if (nsave <= 0) {
- rst->score[0] = rst->score[1] = rst->score[2] = 0;
- return;
- }
-
-#ifndef TFAST
- /* FASTX code here to modify the start, stop points for
- the three phases of the translated protein sequence
- */
- /*
- fprintf(stderr,"n0x: %d; n0x31:%d; n0x32: %d\n",n0,n0x31,n0x32);
- for (ib=0; ib<nsave; ib++) {
- fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
- f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
- f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
- f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
- f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
- }
-
- fprintf(stderr,"---\n");
- */
- for (ib=0; ib<nsave; ib++) {
- if (f_str->noff-f_str->vptr[ib]->dp+f_str->vptr[ib]->start >= n0x32)
- f_str->vptr[ib]->dp += n0x32;
- if (f_str->noff-f_str->vptr[ib]->dp +f_str->vptr[ib]->start >= n0x31)
- f_str->vptr[ib]->dp += n0x31;
- }
-
- /*
- for (ib=0; ib<nsave; ib++) {
- fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
- f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
- f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
- f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
- f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
- }
- */
-#else
-
- /* TFASTX code here to modify the start, stop points for
- the three phases of the translated protein sequence
- TFASTX modifies library start points, rather than
- query start points
- */
-
- /*
- fprintf(stderr,"n0: %d; noff: %d; n1: %d; n1x31: %d n1x32 %d\n",n0, f_str->noff,n1,n1x31,n1x32);
- for (ib=0; ib<nsave; ib++) {
- fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
- f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
- f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
- f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
- f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
- }
-
- fprintf(stderr,"---\n");
- */
-
- for (ib=0; ib<nsave; ib++) {
- if (f_str->vptr[ib]->start >= n1x32) {
- f_str->vptr[ib]->start -= n1x32;
- f_str->vptr[ib]->stop -= n1x32;
- f_str->vptr[ib]->dp -= n1x32;
- }
- if (f_str->vptr[ib]->start >= n1x31) {
- f_str->vptr[ib]->start -= n1x31;
- f_str->vptr[ib]->stop -= n1x31;
- f_str->vptr[ib]->dp -= n1x31;
- }
- }
-
- /*
- for (ib=0; ib<nsave; ib++) {
- fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
- f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
- f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
- f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
- f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
- }
- */
-
-#endif /* TFASTX */
-
- scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap,
- ppst->param_u.fa.pgap, f_str);
-
- for (vmptr=f_str->vptr[0],ib=1; ib<nsave; ib++)
- if (f_str->vptr[ib]->score > vmptr->score) vmptr=f_str->vptr[ib];
-
-/* kssort (f_str->vptr, nsave); */
-
- rst->score[1] = vmptr->score; /* best single score - init1*/
- rst->score[0] = max (scor, vmptr->score); /* initn */
- rst->score[2] = rst->score[0]; /* initn */
-
- my_hoff=f_str->noff - vmptr->dp;
-
- /*
- if (n1 > 5000) {
- fprintf(stderr," Long n1: %d\n",n1);
- }
- */
-
- if (ppst->param_u.fa.optflag) {
- if (rst->score[0] > ppst->param_u.fa.optcut) {
-#ifndef TFAST
- rst->score[2] = dmatchx(aa0, n0,aa1,n1,my_hoff,
- ppst->param_u.fa.optwid, ppst->pam2[ip],
- ppst->gdelval,ppst->ggapval,ppst->gshift,f_str);
-#else /* TFASTX */
- /* generate f_str->aa1y */
-/*
- for (i=0; i<n1; i++) {
- fputc(ppst->sq[aa1[i]],stderr);
- if (i%60==59) fputc('\n',stderr);
- }
- fprintf(stderr,"\n-----\n");
-*/
- for (fs=aa1,itemp=0; itemp <3; itemp++,fs++) {
- for (fd= &f_str->aa1y[itemp]; *fs!=EOSEQ; fd += 3, fs++) *fd = *fs;
- *fd=EOSEQ;
- }
-
-/*
- for (i=0; i<n1; i++) {
- fputc(ppst->sq[f_str->aa1y[i]],stderr);
- if (i%60==59) fputc('\n',stderr);
- }
-*/
- rst->score[2] = dmatchx(aa0, n0, aa1, n1, my_hoff=vmptr->dp-f_str->noff,
- ppst->param_u.fa.optwid, ppst->pam2[ip],
- ppst->gdelval,ppst->ggapval,ppst->gshift,f_str);
-#endif /* TFASTX */
- }
- }
- *hoff = my_hoff;
-}
-
-/* returns rst.score[0] - initn
- rst.score[1] - init1
- rst.score[2] - opt
-*/
-
-void do_work (const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int frame,
- struct pstruct *ppst, struct f_struct *f_str,
- int qr_flg, struct rstruct *rst)
-{
- int hoff;
- int last_n1, itx, itt, n10, i;
-
-#ifdef TFAST
- unsigned char *aa1x;
- /* aa0 has a protein sequence */
- /* aa1 has a raw DNA sequence */
-
- itt = frame;
- last_n1 = 0;
- aa1x = f_str->aa1x;
- for (itx= itt*3; itx< itt*3+3; itx++) {
- n10 = saatran(aa1,&aa1x[last_n1],n1,itx);
- /*
- fprintf(stderr," itt %d itx: %d\n",itt,itx);
- for (i=0; i<n10; i++) {
- fprintf(stderr,"%c",aa[f_str->aa1x[last_n1+i]]);
- if ((i%60)==59) fprintf(stderr,"\n");
- }
- fprintf(stderr,"\n");
- */
- last_n1 += n10+1;
- }
- n10 = last_n1-1;
-#endif
-
- rst->score[0] = rst->score[1] = rst->score[2] = 0;
- rst->escore = 1.0;
- rst->segnum = rst->seglen = 1;
-
-#ifndef TFAST
- do_fastx (f_str->aa0x, n0, aa1, n1, ppst, f_str, rst, &hoff);
-#else /* tfastx */
- do_fastx (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
-#endif
-}
-
-void do_opt (const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int frame,
- struct pstruct *ppst,
- struct f_struct *f_str,
- struct rstruct *rst)
-{
- int optflag, tscore, hoff;
-
- optflag = ppst->param_u.fa.optflag;
- ppst->param_u.fa.optflag = 1;
-
-#ifndef TFAST
- do_fastx (f_str->aa0x, n0, aa1, n1, ppst, f_str, rst, &hoff);
-#else
- do_fastx (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff);
-#endif
-
- ppst->param_u.fa.optflag = optflag;
-}
-
-#ifdef ALLOCN0
-void
-savemax (dptr, dpos, f_str)
- register struct dstruct *dptr;
- int dpos;
- struct f_struct *f_str;
-{
- register struct savestr *vmptr;
- register int i;
-
-#else
-void
-savemax (dptr, f_str)
- register struct dstruct *dptr;
- struct f_struct *f_str;
-{
- register int dpos;
- register struct savestr *vmptr;
- register int i;
-
- dpos = (int) (dptr - f_str->diag);
-
-#endif
-
-/* check to see if this is the continuation of a run that is already saved */
-
- if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
- vmptr->start == dptr->start)
- {
- vmptr->stop = dptr->stop;
- if ((i = dptr->score) <= vmptr->score)
- return;
- vmptr->score = i;
- if (vmptr != f_str->lowmax)
- return;
- }
- else
- {
- i = f_str->lowmax->score = dptr->score;
- f_str->lowmax->dp = dpos;
- f_str->lowmax->start = dptr->start;
- f_str->lowmax->stop = dptr->stop;
- dptr->dmax = f_str->lowmax;
- }
-
- for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
- if (vmptr->score < i)
- {
- i = vmptr->score;
- f_str->lowmax = vmptr;
- }
- f_str->lowscor = i;
-}
-
-int spam (const unsigned char *aa0, const unsigned char *aa1,
- struct savestr *dmax, int **pam2,
- struct f_struct *f_str)
-{
- int lpos;
- int tot, mtot;
- struct {
- int start, stop, score;
- } curv, maxv;
- const unsigned char *aa0p, *aa1p;
-
- aa1p = &aa1[lpos = dmax->start];
- aa0p = &aa0[lpos - dmax->dp + f_str->noff];
- curv.start = lpos;
-
- tot = curv.score = maxv.score = 0;
- for (; lpos <= dmax->stop; lpos++) {
- tot += pam2[*aa0p++][*aa1p++];
- if (tot > curv.score) {
- curv.stop = lpos;
- curv.score = tot;
- }
- else if (tot < 0) {
- if (curv.score > maxv.score) {
- maxv.start = curv.start;
- maxv.stop = curv.stop;
- maxv.score = curv.score;
- }
- tot = curv.score = 0;
- curv.start = lpos+1;
- }
- }
-
- if (curv.score > maxv.score) {
- maxv.start = curv.start;
- maxv.stop = curv.stop;
- maxv.score = curv.score;
- }
-
-/* if (maxv.start != dmax->start || maxv.stop != dmax->stop)
- printf(" new region: %3d %3d %3d %3d\n",maxv.start,
- dmax->start,maxv.stop,dmax->stop);
-*/
- dmax->start = maxv.start;
- dmax->stop = maxv.stop;
-
- return maxv.score;
-}
-
-#define XFACT 10
-
-int sconn (struct savestr **v, int n,
- int cgap, int pgap, struct f_struct *f_str)
-{
- int i, si;
- struct slink {
- int score;
- struct savestr *vp;
- struct slink *next;
- } *start, *sl, *sj, *so, sarr[MAXSAV];
- int lstart, tstart, plstop, ptstop;
-
-/* sort the score left to right in lib pos */
-
- kpsort (v, n);
-
- start = NULL;
-
-/* for the remaining runs, see if they fit */
-
- for (i = 0, si = 0; i < n; i++)
- {
-
-/* if the score is less than the gap penalty, it never helps */
- if (v[i]->score < cgap)
- continue;
- lstart = v[i]->start;
- tstart = lstart - v[i]->dp + f_str->noff;
-
-/* put the run in the group */
- sarr[si].vp = v[i];
- sarr[si].score = v[i]->score;
- sarr[si].next = NULL;
-
-/* if it fits, then increase the score */
- for (sl = start; sl != NULL; sl = sl->next)
- {
- plstop = sl->vp->stop;
- ptstop = plstop - sl->vp->dp + f_str->noff;
- if (plstop < lstart+XFACT && ptstop < tstart+XFACT) {
- sarr[si].score = sl->score + v[i]->score + pgap;
- break;
- }
- }
-
-/* now recalculate where the score fits */
- if (start == NULL)
- start = &sarr[si];
- else
- for (sj = start, so = NULL; sj != NULL; sj = sj->next)
- {
- if (sarr[si].score > sj->score)
- {
- sarr[si].next = sj;
- if (so != NULL)
- so->next = &sarr[si];
- else
- start = &sarr[si];
- break;
- }
- so = sj;
- }
- si++;
- }
-
- if (start != NULL)
- return (start->score);
- else
- return (0);
-}
-
-void
-kssort (v, n)
-struct savestr *v[];
-int n;
-{
- int gap, i, j;
- struct savestr *tmp;
-
- for (gap = n / 2; gap > 0; gap /= 2)
- for (i = gap; i < n; i++)
- for (j = i - gap; j >= 0; j -= gap)
- {
- if (v[j]->score >= v[j + gap]->score)
- break;
- tmp = v[j];
- v[j] = v[j + gap];
- v[j + gap] = tmp;
- }
-}
-
-void
-kpsort (v, n)
-struct savestr *v[];
-int n;
-{
- int gap, i, j;
- struct savestr *tmp;
-
- for (gap = n / 2; gap > 0; gap /= 2)
- for (i = gap; i < n; i++)
- for (j = i - gap; j >= 0; j -= gap)
- {
- if (v[j]->start <= v[j + gap]->start)
- break;
- tmp = v[j];
- v[j] = v[j + gap];
- v[j + gap] = tmp;
- }
-}
-
-static int
-dmatchx(const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int hoff, int window,
- int **pam2, int gdelval, int ggapval, int gshift,
- struct f_struct *f_str)
-{
-
- hoff -= window/2;
-
-#ifndef TFAST
- return lx_band(aa1,n1,f_str->aa0y,n0,
- pam2,
-#ifdef OLD_FASTA_GAP
- -(gdelval-ggapval),
-#else
- -gdelval,
-#endif
- -ggapval,-gshift,
- hoff,window,f_str);
-#else
- return lx_band(aa0,n0,f_str->aa1y,n1,
- pam2,
-#ifdef OLD_FASTA_GAP
- -(gdelval-ggapval),
-#else
- -gdelval,
-#endif
- -ggapval,-gshift,
- hoff,window,f_str);
-#endif
-}
-
-static void
-init_row(struct sx_s *row, int sp) {
- int i;
- for (i = 0; i < sp; i++) {
- row[i].C1 = row[i].I1 = 0;
- row[i].C2 = row[i].I2 = 0;
- row[i].C3 = row[i].I3 = 0;
- row[i].flag = 0;
- }
-}
-
-int
-lx_band(const unsigned char *prot_seq, /* array with protein sequence numbers*/
- int len_prot, /* length of prot. seq */
- const unsigned char *dna_prot_seq, /* translated DNA sequence numbers*/
- int len_dna_prot, /* length trans. seq. */
- int **pam_matrix, /* scoring matrix */
- int gopen, int gext, /* gap open, gap extend penalties */
- int gshift, /* frame-shift penalty */
- int start_diag, /* start diagonal of band */
- int width, /* width for band alignment */
- struct f_struct *f_str)
-{
- void *ckalloc();
- int i, j, bd, bd1, x1, sp, p1=0, p2=0, end_prot;
- int sc, del, best = 0, cd,ci, e1, e2, e3, cd1, cd2, cd3, f, gg;
- register int *wt;
- const unsigned char *dp;
- register struct sx_s *ap, *aq;
-
- sp = width+7;
- gg = gopen+gext;
- /* sp = sp/3; */
- if (f_str->cur == NULL)
- f_str->cur = (struct sx_s *) ckalloc(sizeof(struct sx_s)*sp);
-
- init_row(f_str->cur, sp);
-
- /*
- if (start_diag %3 !=0) start_diag = start_diag/3-1;
- else start_diag = start_diag/3;
- */
-
- /*
- if (width % 3 != 0) width = width/3+1;
- else width = width /3;
- */
-
- /* currently, this code assumes that the DNA sequence is longer than the
- protein sequence. This is not always true. len_prot in the loop below
- should be decreased to the projection of the DNA on the protein */
-
- x1 = start_diag; /* x1 = lower bound of DNA */
-
-
- end_prot = max(0,-width-start_diag) + (len_dna_prot+5)/3 + width;
- end_prot = min(end_prot,len_prot);
-
- /* i counts through protein sequence, x1 through DNAp */
-
- for (i = max(0, -width-start_diag), x1+=i; i < end_prot; i++, x1++) {
- bd = min(x1+width, len_dna_prot/3); /* upper bound of band */
- bd1 = max(0,x1); /* lower bound of band */
- wt = pam_matrix[prot_seq[i]];
- del = 1-x1; /*adjustment*/
- bd += del;
- bd1 +=del;
-
- ap = &f_str->cur[bd1];
- aq = ap+1;
- e1 = f_str->cur[bd1-1].C3;
- e2 = ap->C1;
- cd1 = cd2= cd3= 0;
-
- for (dp = &dna_prot_seq[(bd1-del)*3]; ap < &f_str->cur[bd]; ap++) {
- sc = max(max(e1, (e3=ap->C2))-gshift, e2)+wt[*dp++];
- if (cd1 > sc) sc = cd1;
- cd1 -= gext;
- if ((ci = aq->I1) > 0) {
- if (sc < ci) { ap->C1 = ci; ap->I1 = ci-gext;}
- else {
- ap->C1 = sc;
- sc -= gg;
- if (sc > 0) {
- if (sc > best) best =sc;
- if (cd1 < sc) cd1 = sc;
- ap->I1 = max(ci-gext, sc);
- } else ap->I1 = ci-gext;
- }
- } else {
- if (sc <= 0) {
- ap->I1 = ap->C1 = 0;
- } else {
- ap->C1 = sc; sc-=gg;
- if (sc >0) {
- if (sc > best) best =sc;
- if (cd1 < sc) cd1 = sc;
- ap->I1 = sc;
- } else ap->I1 = 0;
- }
- }
- sc = max(max(e2, (e1=ap->C3))-gshift, e3)+wt[*dp++];
- if (cd2 > sc) sc = cd2;
- cd2 -= gext;
- if ((ci = aq->I2) > 0) {
- if (sc < ci) { ap->C2 = ci; ap->I2 = ci-gext;}
- else {
- ap->C2 = sc;
- sc -= gg;
- if (sc > 0) {
- if (sc > best) best =sc;
- if (cd2 < sc) cd2 = sc;
- ap->I2 = max(ci-gext, sc);
- }
- }
- } else {
- if (sc <= 0) {
- ap->I2 = ap->C2 = 0;
- } else {
- ap->C2 = sc; sc-=gg;
- if (sc >0) {
- if (sc > best) best =sc;
- if (cd2 < sc) cd2 = sc;
- ap->I2 = sc;
- } else ap->I2 = 0;
- }
- }
- sc = max(max(e3, (e2=aq->C1))-gshift, e1)+wt[*dp++];
- if (cd3 > sc) sc = cd3;
- cd3 -= gext;
- if ((ci = aq++->I3) > 0) {
- if (sc < ci) { ap->C3 = ci; ap->I3 = ci-gext;}
- else {
- ap->C3 = sc;
- sc -= gg;
- if (sc > 0) {
- if (sc > best) best =sc;
- if (cd3 < sc) cd3 = sc;
- ap->I3 = max(ci-gext, sc);
- }
- }
- } else {
- if (sc <= 0) {
- ap->I3 = ap->C3 = 0;
- } else {
- ap->C3 = sc; sc-=gg;
- if (sc >0) {
- if (sc > best) best =sc;
- if (cd3 < sc) cd3 = sc;
- ap->I3 = sc;
- } else ap->I3 = 0;
- }
- }
- }
- }
- /* printf("The best score is %d\n", best); */
- return best+gopen+gext;
-}
-
-/* ckalloc - allocate space; check for success */
-void *ckalloc(size_t amount)
-{
- void *p;
-
- if ((p = (void *)malloc( (size_t)amount)) == NULL)
- w_abort("Ran out of memory.","");
- return(p);
-}
-
-/* calculate the 100% identical score */
-int
-shscore(unsigned char *aa0, int n0, int **pam2)
-{
- int i, sum;
- for (i=0,sum=0; i<n0; i++)
- sum += pam2[aa0[i]][aa0[i]];
- return sum;
-}
-
-#define SGW1 100
-#define SGW2 300
-#define WIDTH 60
-
-/* code above is to convert sequence into numbers */
-
-typedef struct mat *match_ptr;
-
-typedef struct mat {
- int i, j, l;
- match_ptr next;
-} match_node;
-
-typedef struct {
- int i,j;
-} state;
-
-typedef state *state_ptr;
-
-typedef struct st_s { int C, I, D;} *st_ptr;
-
-/* static st_ptr up=NULL, down, tp; */
-/* static int *st_up; */
-/* static int gop, gext, shift; */
-
-void *ckalloc(size_t);
-static match_ptr small_global(), global();
-static int local_align(), find_best();
-static void init_row2(), init_ROW();
-
-int
-pro_dna(const unsigned char *prot_seq, /* array with prot. seq. numbers*/
- int len_prot, /* length of prot. seq */
- const unsigned char *dna_prot_seq, /* trans. DNA seq. numbers*/
- int len_dna_prot, /* length trans. seq. */
- int **pam_matrix, /* scoring matrix */
- int gopen, int gex, /* gap open, gap extend penalties */
- int gshift, /* frame-shift penalty */
- int max_res,
- struct a_res_str *a_res) /* alignment info */
-{
- match_ptr align, ap, aq;
- int x, y, ex, ey, i, score;
- int *alignment;
- st_ptr up, down, tp;
-
- /* these globals removed */
- /* gext = gex; gop = gopen; shift = gshift; */
-
- /* for fastx (but not tfastx), these could be moved into init_work(),
- and done only once */
-
- up = (st_ptr) ckalloc(sizeof(struct st_s)*(len_dna_prot+10));
- down = (st_ptr) ckalloc(sizeof(struct st_s)*(len_dna_prot+10));
- tp = (st_ptr) ckalloc(sizeof(struct st_s)*(len_dna_prot+10));
-
- /*local alignment find the best local alignment x (prot) and y (DNA)
- is the starting position of the best local alignment
- and ex (prot) ey (DNA) is the ending position */
- score= local_align(&x, &y, &ex, &ey, pam_matrix,
- gopen, gex, gshift,
- dna_prot_seq, len_dna_prot,
- prot_seq, len_prot, up, down);
-
- /* this is very strange, since local_align initialized up, down */
- up += 3; down += 3; tp += 3;
-
- /* x, y - start in prot, dna_prot */
- a_res->min0 = x; /* prot */
- a_res->max0 = ex; /* prot */
-
- a_res->min1 = y; /* DNA-prot */
- a_res->max1 = ey; /* DNA-prot */
-
- align = global(x, y, ex, ey, pam_matrix, gopen, gex, gshift,
- dna_prot_seq, prot_seq, 0, 0, &up, &down, &tp);
-
- alignment = a_res->res;
-
- /* from earlier version */
- /* alignment[0] = x; */ /* start of alignment in prot */
- /* alignment[1] = y; */ /* start of alignment in DNA */
-
- for (ap = align, i= 0; ap; i++) {
- if (i < max_res) {alignment[i] = ap->l;}
- aq = ap->next; free(ap); ap = aq;
- }
-
- if (i >= max_res) {
- fprintf(stderr," alignment truncated: %d/%d\n", max_res,i);
- }
-
- up = &up[-3]; down = &down[-3]; tp = &tp[-3];
- free(up); free(tp); free(down);
- /* free(st_up); */ /* moved into local align */
-
- a_res->nres = i; /* i has the length of the alignment */
- return score;
-}
-
-static void
-swap(void **a, void **b) {
- void *t;
-
- t = *a;
- *a = *b;
- *b = t;
-}
-
-/*
- local alignment find the best local alignment x and y
- is the starting position of the best local alignment
- and ex ey is the ending position
-*/
-static int
-local_align(int *x, int *y, int *ex, int *ey,
- int **wgts, int gop, int gext, int shift,
- unsigned char *dnap, int ld,
- unsigned char *pro, int lp,
- st_ptr up, st_ptr down) {
-
- int i, j, score, x1,x2,x3,x4, e1, e2 = 0, e3,
- sc, del, e, best = 0, *wt, cd, ci;
- state_ptr cur_st, last_st, cur_i_st;
- st_ptr cur, last;
- unsigned char *dp;
- int *st_up, *cur_d_st;
-
-/*
- Array rowiC store the best scores of alignment ending at a position
- Arrays rowiD, and rowiI store the best scores of alignment ending
- at a position with a deletion or insrtion
- Arrays sti stores the starting position of the best alignment whose
- score stored in the corresponding row array.
- The program stores two rows to complete the computation, same is
- for the global alignment routine.
-*/
-
- /* for fastx (but not tfastx), this could be moved into init_work(),
- and done only once */
- st_up = (int *) ckalloc(sizeof(int)*(ld+10));
- init_row2(st_up, ld+5);
-
- ld += 2;
- init_ROW(up, ld+1); /* set to zero */
- init_ROW(down, ld+1); /* set to zero */
-
-
- cur = up+1;
- last = down+1;
-
- /* for fastx (but not tfastx), these could be moved into init_work(),
- and done only once */
- cur_st = (state_ptr) ckalloc(sizeof(state)*(ld+1));
- last_st = (state_ptr) ckalloc(sizeof(state)*(ld+1));
- cur_i_st = (state_ptr) ckalloc(sizeof(state)*(ld+1));
-
- cur_d_st = st_up;
-
- dp = dnap-2;
- for (i = 0; i < lp; i++) {
- wt = &wgts[pro[i]][0];
- for (j = 0; j < 2; j++) {
- cur_st[j].i = i+1;
- cur_st[j].j = j+1;
- }
- for (j = 2; j < ld; j++) {
- score = wt[dp[j]];
- del = -1;
- if (j >= 3) {
- sc = -score;
- e3 = e2-shift; e2 = last[j-3].C;
- e1 = last[j-2].C-shift;
- if (e1 > sc) {sc = e1; del = 2;}
- if (e2 > sc) {sc = e2; del = 3;}
- if (e3 > sc) {sc = e3; del = 4;}
- } else {
- sc = e2 = 0;
- if (sc < -score) sc=-score;
- else del = 3;
- }
- sc += score;
- if (sc < (ci=last[j].I)) {
- sc = ci; del = 0;
- }
- if (sc < (cd=cur[j].D)) {
- sc = cd; del = 5;
- }
- cur[j].C = sc;
- e = sc - gop;
- if (e > cd) {
- cur[j+3].D = e-gext;
- cur_d_st[j+3] = 3;
- } else {
- cur[j+3].D = cd-gext;
- cur_d_st[j+3] = cur_d_st[j]+3;
- }
- switch(del) {
- case 5:
- e1 = cur_d_st[j];
- cur_st[j].i = cur_st[j-e1].i;
- cur_st[j].j = cur_st[j-e1].j;
- break;
- case 0:
- cur_st[j].i = cur_i_st[j].i;
- cur_st[j].j = cur_i_st[j].j;
- break;
- case 2:
- case 3:
- case 4:
- if (i) {
- if (j-del >= 0) {
- cur_st[j].i = last_st[j-del].i;
- cur_st[j].j = last_st[j-del].j;
- } else {
- cur_st[j].i = i;
- cur_st[j].j = 0;
- }
- } else {
- cur_st[j].i = 0;
- cur_st[j].j = max(0, j-del+1);
- }
- break;
- case -1:
- cur_st[j].i = i+1;
- cur_st[j].j = j+1;
- break;
- }
- if (e > ci) {
- cur[j].I = e -gext;
- cur_i_st[j].i = cur_st[j].i;
- cur_i_st[j].j = cur_st[j].j;
- } else {
- cur[j].I = ci- gext;
- }
- if (sc > best) {
- x1 = cur_st[j].i;
- x2 = cur_st[j].j;
- best =sc;
- x3 = i;
- x4 = j;
- }
- }
- swap((void **)&last, (void **)&cur);
- swap((void **)&cur_st, (void **)&last_st);
- }
- /* printf("The best score is %d\n", best); */
- *x = x1; *y = x2; *ex = x3; *ey = x4;
- free(cur_st); free(last_st); free(cur_i_st);
- free(st_up);
- return best;
-}
-
-/*
- Both global_up and global_down do linear space score only global
- alignments on subsequence pro[x]...pro[ex], and dna[y]...dna[ey].
- global_up do the algorithm upwards, from row x towards row y.
- global_down do the algorithm downwards, from row y towards x.
-*/
-
-static void
-global_up(st_ptr *row1, st_ptr *row2,
- int x, int y, int ex, int ey,
- int **wgts, int gop, int gext, int shift,
- unsigned char *dnap,
- unsigned char *pro,
- int N) {
- int i, j, k, sc, e, e1, e2, e3, t, ci, cd, score, *wt;
- st_ptr cur, last;
-
- cur = *row1; last = *row2;
- sc = -gop-gext;
- for (j = 1; j <= ey-y+1; j++) {
- if (j % 3 == 0) {last[j].C = sc; sc -= gext; last[j].I = sc-gop;}
- else { last[j].I = last[j].C = -10000;}
- cur[j].I = -10000;
- }
- last[0].C = 0; cur[0].D = cur[1].D = cur[2].D = -10000;
- last[0].D = last[1].D = last[2].D = -10000;
- if (N) last[0].I = -gext; else last[0].I = -gop-gext;
- for (i = 1; i <= ex-x+1; i++) {
- wt = &wgts[pro[i+x-1]][0]; e2 = last[0].C; e1 = -10000;
- for (j = 0; j <= ey-y+1; j++) {
- t = j+y;
- sc = -10000;
- if (t < 3) score = -10000;
- else score = wt[dnap[t-3]];
- if (j < 4) {
- if (j == 3) sc = e2;
- else if (j == 2) sc = e2-shift;
- } else {
- e3 = e2; e2 = e1;
- e1 = last[j-2].C;
- sc = max(max(e1, e3)-shift, e2);
- }
- sc += score;
- sc = max(sc, max(ci=last[j].I, cd = cur[j].D));
- cur[j].C = sc;
- cur[j+3].D = max(cd, sc-gop)-gext;
- cur[j].I = max(ci, sc-gop)-gext;
- }
- swap((void **)&last, (void **)&cur);
- }
- for (i = 0; i <= ey-y+1; i++) last[i].I = cur[i].I;
- if (*row1 != last) swap((void **)row1, (void **)row2);
-}
-
-static void
-global_down(st_ptr *row1, st_ptr *row2,
- int x, int y, int ex, int ey,
- int **wgts, int gop, int gext, int shift,
- unsigned char *dnap, unsigned char *pro,
- int N) {
- int i, j, k, sc, del, *tmp, e, t, e1,e2,e3, ci,cd, s1, s2, s3, *wt;
- st_ptr cur, last;
-
- cur = (*row1); last = *row2;
- sc = -gop-gext;
- for (j = ey-y; j >= 0; j--) {
- if ((ey-y+1-j) % 3) {last[j].C = sc; sc-=gext; last[j].I = sc-gop;}
- else last[j].I = last[j].C = -10000;
- }
- last[ey-y+1].C = 0;
- cur[ey-y+1].D = cur[ey-y].D = cur[ey-y-1].D = -10000;
- last[ey-y+1].D = last[ey-y].D = last[ey-y-1].D = -10000;
- if (N) last[ey-y+1].I = -gext; else last[ey-y+1].I = -gop-gext;
- for (i = ex-x; i >= 0; i--) {
- wt = &wgts[pro[i+x]][0]; e2 = last[ey-y+1].C;
- e1 = s2 = s3 = -10000;
- for (j = ey-y+1; j >= 0; j--) {
- t = j+y;
- s1 = wt[dnap[t-1]];
- sc = -10000;
- if (t+3 > ey) {
- if (t+2==ey) sc = e2+s2;
- else if (t+1==ey) sc = e2-shift+s1;
- } else {
- e3 = e2; e2 = e1;
- e1 = last[j+2].C;
- sc = max(max(e1+s1, e3+s3)-shift, e2+s2);
- }
- if (sc < (cd= cur[j].D)) {
- sc = cd;
- cur[j-3].D = cd-gext;
- } else cur[j-3].D =max(cd, sc-gop)-gext;
- if (sc < (ci= last[j].I)) {
- sc = ci; del = 0;
- cur[j].I = ci - gext;
- } else cur[j].I = max(sc-gop,ci)-gext;
- cur[j].C = sc;
- s3 = s2; s2 = s1;
- }
- swap((void **)&last, (void **)&cur);
- }
- for (i = 0; i <= ey-y+1; i++) last[i].I = cur[i].I;
- if (*row1 != last) swap((void **)row1, (void **)row2);
-}
-
-static void
-init_row2(int *row, int ld) {
- int i;
- for (i = 0; i < ld; i++) row[i] = 0;
-}
-
-static void
-init_ROW(st_ptr row, int ld) {
- int i;
- for (i = 0; i < ld; i++) row[i].I = row[i].D = row[i].C = 0;
-}
-
-static match_ptr
-combine(match_ptr x1, match_ptr x2, int st) {
- match_ptr x;
-
- if (x1 == NULL) return x2;
- for (x = x1; x->next; x = x->next);
- x->next = x2;
- if (st) {
- for (x = x2; x; x = x->next) {
- x->j++;
- if (x->l == 3 || x->l == 4) break;
- }
- x->l--;
- }
- return x1;
-}
-
-/*
- global use the two upwards and downwards score only linear
- space global alignment subroutine to recursively build the
- alignment.
-*/
-
-match_ptr
-global(int x, int y, int ex, int ey,
- int **wgts, int gop, int gext, int shift,
- unsigned char *dnap,
- unsigned char *pro,
- int N1, int N2,
- st_ptr *up_stp, st_ptr *dn_stp, st_ptr *tp_stp
- )
-{
- int m;
- int m1, m2;
- match_ptr x1, x2, mm1, mm2;
- /*printf("%d %d %d %d\n", x,y, ex, ey);*/
- /*
- if the space required is limited, we can do a quadratic space
- algorithm to find the alignment.
- */
- if (ex <= x) {
- mm1 = NULL; mm2= NULL;
- for (m = y+3; m <= ey; m+=3) {
- x1 = (match_ptr) ckalloc(sizeof(match_node));
- x1->l = 5; x1->next = mm1;
- if (mm1== NULL) mm2 = x1;
- mm1 = x1;
- }
- if (ex == x) {
- if ((ey-y) % 3 != 0) {
- x1 = (match_ptr) ckalloc(sizeof(match_node));
- x1->l = ((ey-y) % 3) +1; x1->next = NULL;
- if (mm2) mm2->next = x1;
- else mm1 = x1;
- } else {
- if (mm2) mm2->l = 4;
- }
- }
- return mm1;
- }
- if (ey <= y) {
- mm1 = NULL;
- for (m = x; m <= ex; m++) {
- x1 = (match_ptr) ckalloc(sizeof(match_node));
- x1->l = 0; x1->next = mm1; mm1 = x1;
- }
- return mm1;
- }
- if (ex -x < SGW1-1 && ey-y < SGW2-1)
- return small_global(x,y,ex,ey,
- wgts, gop, gext, shift,
- dnap, pro, N1, N2);
- m = (x+ex)/2;
- /*
- Do the score only global alignment from row x to row m, m is
- the middle row of x and ex. Store the information of row m in
- upC, upD, and upI.
- */
- global_up(up_stp, tp_stp, x, y, m, ey,
- wgts, gop, gext, shift,
- dnap, pro, N1);
-
- /*
- Do the score only global alignment downwards from row ex
- to row m+1, store information of row m+1 in downC downI and downD
- */
- global_down(dn_stp, tp_stp, m+1, y, ex, ey,
- wgts, gop, gext, shift,
- dnap, pro, N2);
-
- /*
- Use these information of row m and m+1, to find the crossing
- point of the best alignment with the middle row. The crossing
- point is given by m1 and m2. Then we recursively call global
- itself to compute alignments in two smaller regions found by
- the crossing point and combine the two alignments to form a
- whole alignment. Return that alignment.
- */
- if (find_best(*up_stp, *dn_stp, &m1, &m2, ey-y+1, y, gop)) {
- x1 = global(x, y, m, m1, wgts, gop, gext, shift, dnap, pro, N1, 0,
- up_stp, dn_stp, tp_stp);
- x2 = global(m+1, m2, ex, ey, wgts, gop, gext, shift, dnap, pro, 0, N2,
- up_stp, dn_stp, tp_stp);
- if (m1 == m2) x1 = combine(x1,x2,1);
- else x1 = combine(x1, x2,0);
- } else {
- x1 = global(x, y, m-1, m1, wgts, gop, gext, shift, dnap, pro, N1, 1,
- up_stp, dn_stp, tp_stp);
- x2 = global(m+2, m2, ex, ey, wgts, gop, gext, shift, dnap, pro, 1, N2,
- up_stp, dn_stp, tp_stp);
- mm1 = (match_ptr) ckalloc(sizeof(match_node));
- mm1->i = m; mm1->l = 0; mm1->j = m1;
- mm2 = (match_ptr) ckalloc(sizeof(match_node));
- mm2->i = m+1; mm2->l = 0; mm2->j = m1;
- mm1->next = mm2; mm2->next = x2;
- x1 = combine(x1, mm1, 0);
- }
- return x1;
-}
-
-static int
-find_best(st_ptr up, st_ptr down,
- int *m1, int *m2,
- int ld, int y, int gop) {
- int i, best = -100000, j = 0, s1, s2, s3, s4, st;
- up++;
- for (i = 1; i < ld; i++) {
- s2 = up[i-1].C + down[i].C;
- s4 = up[i-1].I + down[i].I + gop;
- if (best < s2) {
- best = s2; j = i; st = 1;
- }
- if (best < s4) {
- best = s4; j = i; st = 0;
- }
- }
- *m1 = j-1+y;
- *m2 = j+y;
- /*printf("find best score =%d\n", best);*/
- return st;
-}
-
-/*
- An alignment is represented as a linked list whose element
- is of type match_node. Each element represent an edge in the
- path of the alignment graph. The fields of match_node are
- l --- gives the type of the edge.
- i, j --- give the end position.
-*/
-
-static match_ptr
-small_global(int x, int y, int ex, int ey,
- int **wgts, int gop, int gext, int shift,
- unsigned char *dnap, unsigned char *pro,
- int N1, int N2) {
- static int C[SGW1+1][SGW2+1], st[SGW1+1][SGW2+1], D[SGW2+7], I[SGW2+1];
- int i, j, e, sc, score, del, k, t, *wt, ci, cd;
- int *cI, *cD, *cC, *lC, *cst, e2, e3, e4;
- match_ptr mp, first;
-
- /*printf("small_global %d %d %d %d\n", x, y, ex, ey);*/
- sc = -gop-gext; C[0][0] = 0;
- if (N1) I[0] = -gext; else I[0] = sc;
- for (j = 1; j <= ey-y+1; j++) {
- if (j % 3== 0) {
- C[0][j] = sc; sc -= gext; I[j] = sc-gop;
- } else I[j] = C[0][j] = -10000;
- st[0][j] = 5;
- }
- lC = &C[0][0]; cD = D; D[0] = D[1] = D[2] = -10000;
- cI = I;
- for (i = 1; i <= ex-x+1; i++) {
- cC = &C[i][0];
- wt = &wgts[pro[i+x-1]][0]; cst = &st[i][0];
- for (j = 0; j <=ey-y+1; j++) {
- sc = -10000; del = 0;
- ci = cI[j];
- cd= cD[j];
- t = j+y;
- if (t < 3) score = -10000;
- else score = wt[dnap[t-3]];
- if (j >= 4) {
- e2 = lC[j-2]-shift; sc = lC[j-3]; e4 = lC[j-4]-shift;
- del = 3;
- if (e2 > sc) { sc = e2; del = 2;}
- if (e4 >= sc) { sc = e4; del = 4;}
- } else {
- if (j ==3) {sc= lC[0]; del = 3;}
- else if (j == 2) {sc = lC[0]-shift; del = 2;}
- }
- sc = sc+score;
- if (sc < ci) {
- sc = ci; del = 0;
- }
- if (sc <= cd) {
- sc = cd;
- del = 5;
- }
- cC[j] = sc;
- sc -= gop;
- if (sc < cd) {
- del += 10;
- cD[j+3] = cd - gext;
- } else cD[j+3] = sc -gext;
- if (sc < ci) {
- del += 20;
- cI[j] = ci-gext;
- } else cI[j] = sc-gext;
- *(cst++) = del;
- }
- lC = cC;
- }
- if (N2 && ci +gop > cC[ey-y+1]) {
- st[ex-x+1][ey-y+1] = 0;
- /*printf("small score = %d\n", ci+gop);*/
- } /*else printf("small score =%d\n", cC[ey-y+1]);*/
- first = NULL; e = 1;
- for (i = ex+1, j = ey+1; i > x || j > y; i--) {
- mp = (match_ptr) ckalloc(sizeof(match_node));
- mp->i = i-1;
- k = (t=st[i-x][j-y])%10;
- mp->j = j-1;
- if (e == 5 && (t/10)%2 == 1) k = 5;
- if (e == 0 && (t/20)== 1) k = 0;
- if (k == 5) { j -= 3; i++; e=5;}
- else {j -= k;if (k==0) e= 0; else e = 1;}
- mp->l = k;
- mp->next = first;
- first = mp;
- }
-
- /* for (i = 0; i <= ex-x; i++) {
- for (j = 0; j <= ey-y; j++)
- printf("%d ", C[i][j]);
- printf("\n");
- }
- */
- return first;
-}
-
-
-#define XTERNAL
-#include "upam.h"
-
-extern void display_alig(a, dna, pro,length, ld)
-int *a;
-unsigned char *dna, *pro;
-int length, ld;
-{
- int len = 0, i, j, x, y, lines, k;
- static char line1[100], line2[100], line3[100],
- tmp[10] = " ";
- unsigned char *dna1, c1, c2, c3, *st;
-
- dna1 = ckalloc((size_t)ld);
- for (st = dna, i = 0; i < ld; i++, st++) dna1[i] = aa[*st];
- line1[0] = line2[0] = line3[0] = '\0'; x= a[0]; y = a[1]-1;
-
- for (len = 0, j = 2, lines = 0; j < length; j++) {
- i = a[j];
- /*printf("%d %d %d\n", i, len, b->j);*/
- if (i > 0 && i < 5) tmp[i-2] = aa[pro[x++]];
- if (i == 5) {
- i = 3; tmp[0] = tmp[1] = tmp[2] = '-';
- if (a[j+1] == 2) tmp[2] = ' ';
- }
- if (i > 0) {
- strncpy(&line1[len], (const char *)&dna1[y], i); y+=i;
- } else {line1[len] = '-'; i = 1; tmp[0] = aa[pro[x++]];}
- strncpy(&line2[len], tmp, i);
- for (k = 0; k < i; k++) {
- if (tmp[k] != ' ' && tmp[k] != '-') {
- if (k == 2) tmp[k] = '\\';
- else if (k == 1) tmp[k] = '|';
- else tmp[k] = '/';
- } else tmp[k] = ' ';
- }
- if (i == 1) tmp[0] = ' ';
- strncpy(&line3[len], tmp, i);
- tmp[0] = tmp[1] = tmp[2] = ' ';
- len += i;
- line1[len] = line2[len] =line3[len] = '\0';
- if (len >= WIDTH) {
- printf("\n%5d", WIDTH*lines++);
- for (k = 10; k <= WIDTH; k+=10)
- printf(" . :");
- if (k-5 < WIDTH) printf(" .");
- c1 = line1[WIDTH]; c2 = line2[WIDTH]; c3 = line3[WIDTH];
- line1[WIDTH] = line2[WIDTH] = line3[WIDTH] = '\0';
- printf("\n %s\n %s\n %s\n", line1, line3, line2);
- line1[WIDTH] = c1; line2[WIDTH] = c2; line3[WIDTH] = c3;
- strncpy(line1, &line1[WIDTH], sizeof(line1)-1);
- strncpy(line2, &line2[WIDTH], sizeof(line2)-1);
- strncpy(line3, &line3[WIDTH], sizeof(line3)-1);
- len = len - WIDTH;
- }
- }
- printf("\n%5d", WIDTH*lines);
- for (k = 10; k < len; k+=10)
- printf(" . :");
- if (k-5 < len) printf(" .");
- printf("\n %s\n %s\n %s\n", line1, line3, line2);
-}
-
-
-/* alignment store the operation that align the protein and dna sequence.
- The code of the number in the array is as follows:
- 0: delete of an amino acid.
- 2: frame shift, 2 nucleotides match with an amino acid
- 3: match an amino acid with a codon
- 4: the other type of frame shift
- 5: delete of a codon
-
-
- Also the first two element of the array stores the starting point
- in the protein and dna sequences in the local alignment.
-
- Display looks like where WIDTH is assumed to be divisible by 10.
-
- 0 . : . : . : . : . : . :
- CCTATGATACTGGGATACTGGAACGTCCGCGGACTGACACACCCGATCCGCATGCTCCTG
- P M I L G Y W N V R G L T H P I R M L L
-
- 60 . : . : . : . : . : . :
- GAATACACAGACTCAAGCTATGATGAGAAGAGATACACCATGGGTGACGCTCCCGACTTT
- E Y T D S S Y D E K R Y T M G D A P D F
-*/
-
-
-/* fatal - print message and die */
-void fatal(msg)
-char *msg;
-{
- fprintf(stderr, "%s\n", msg);
- exit(1);
-}
-
-int do_walign (const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int frame,
- struct pstruct *ppst,
- struct f_struct *f_str,
- struct a_res_str *a_res,
- int *have_ares)
-{
- int score;
- int i, last_n1, itemp, n10;
- int n_aa, n_nt, hoff, nt_min, nt_max, w_fact;
- unsigned char *fs, *fd;
- struct rstruct rst;
- int itx;
-
-#ifndef TFAST /* FASTX */
- n_aa = n1;
- n_nt = n0;
-
- /* check for large differences in sequence length */
- nt_min = 0; nt_max = n_nt;
- if (n_nt > 6 * n_aa) {
- /* find out where the diagonal is - get hoff
- hoff < 0 => seq0 is in the middle of seq1
- */
- do_fastx(f_str->aa0x, n0, aa1, n1, ppst, f_str, &rst, &hoff);
- if (rst.score[0] > 2 * rst.score[2]) {w_fact = 4;}
- else w_fact = 2;
-
- if (hoff > n_aa) { /* hoff > 0 => seq1 is in the middle of seq0 */
- nt_min = max(0,(hoff-w_fact*n_aa)*3);
- nt_max = min((hoff+w_fact*n_aa)*3,n_nt);
- }
- else {
- nt_max = min(3*w_fact*n_aa,n_nt);
- }
- }
-
- a_res->res = f_str->res;
-
- score = pro_dna(aa1, n1, f_str->aa0y+nt_min, nt_max-nt_min, ppst->pam2[0],
-#ifdef OLD_FASTA_GAP
- -(ppst->gdelval - ppst->ggapval),
-#else
- -ppst->gdelval,
-#endif
- -ppst->ggapval,
- -ppst->gshift,
- f_str->max_res, a_res);
-
- /* correct for nt_min missing residues in alignment */
-
-#else /* TFASTX */
-
- /*
- for (i=0; i<n1; i++) {
- fputc(ppst->sq[f_str->aa1x[i]],stderr);
- if (i%60==59) fputc('\n',stderr);
- }
- fprintf(stderr,"\n-----\n");
- */
-
- last_n1 = 0;
- for (itx=3*frame; itx<3+3*frame; itx++) {
- n10 = saatran(aa1,&f_str->aa1x[last_n1],n1,itx);
-/*
- for (i=0; i<n10; i++) {
- fprintf(stderr,"%c",pst.sq[aa10[last_n1+i]]);
- if ((i%60)==59) fprintf(stderr,"\n");
- }
- fprintf(stderr,"\n");
-*/
- last_n1 += n10+1;
- }
- n10 = last_n1-1;
-
- /* create aa1y from aa1x */
- for (fs=f_str->aa1x,itemp=0; itemp <3; itemp++,fs++) {
- for (fd= &f_str->aa1y[itemp]; *fs!=EOSEQ; fd += 3, fs++) *fd = *fs;
- *fd=EOSEQ;
- }
- /*
- for (i=0; i<n1; i++) {
- fputc(ppst->sq[f_str->aa1y[i]],stderr);
- if (i%60==59) fputc('\n',stderr);
- }
- fprintf(stderr,"\n-----\n");
- */
-
- n_aa = n0;
- n_nt = n1;
-
- /* check for large differences in sequence length */
- nt_min = 0; nt_max = n_nt;
- if (n_nt > 6 * n_aa) {
- /* find out where the diagonal is - get hoff
- hoff < 0 => seq0 is in the middle of seq1
- */
- do_fastx(aa0, n0, f_str->aa1x, n10, ppst, f_str, &rst, &hoff);
- if (rst.score[0] > 2 * rst.score[2]) {w_fact = 4;}
- else w_fact = 2;
-
- if ( hoff > n_aa) { /* hoff > 0 => seq1 is in the middle of seq0 */
- nt_min = max(0,(hoff-w_fact*n_aa)*3);
- nt_max = min((hoff+w_fact*n_aa)*3,n_nt);
- }
- else {
- nt_max = min(3*w_fact*n_aa,n_nt);
- }
- }
-
- a_res->res = f_str->res;
-
- score = pro_dna(aa0, n0, f_str->aa1y+nt_min, nt_max-nt_min, ppst->pam2[0],
-#ifdef OLD_FASTA_GAP
- -(ppst->gdelval - ppst->ggapval),
-#else
- -ppst->gdelval,
-#endif
- -ppst->ggapval,
- -ppst->gshift,
- f_str->max_res, a_res);
-
-#endif /* TFASTX */
-
- /* pro_dna always compares protein to DNA, and returns protein
- coordinates in a_res->min0,max0 */
-
- a_res->min1 += nt_min;
- a_res->max1 += nt_min;
-
- /* display_alig(f_str->res,f_str->aa0y,aa1,*nres,n0); */
-
- *have_ares = 1;
- return score;
-}
-
-/* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
-/* call from calcons, calc_id, calc_code */
-void
-aln_func_vals(int frame, struct a_struct *aln) {
-
-#ifndef TFAST
- aln->llrev = 0;
- aln->llfact = 1;
- aln->llmult = 1;
- aln->qlfact = 3;
- aln->frame = 0;
- if (frame > 0) aln->qlrev = 1;
- else aln->qlrev = 0;
-#else /* TFASTX */
- aln->qlfact = 1;
- aln->qlrev = 0;
- aln->llfact = 3;
- aln->llmult = 1;
- aln->frame = 0;
- if (frame > 0) aln->llrev = 1;
- else aln->llrev = 0;
-#endif /* TFASTX */
-}
-
-/* this function is required for programs like tfastx/y/s that do
- translations on DNA sequences and save them in f_str->aa1??
-*/
-
-void
-pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {
-#ifdef TFAST
- int i, last_n1, itemp, n10;
- unsigned char *fs, *fd;
- int itx;
-
- last_n1 = 0;
- for (itx=3*frame; itx<3+3*frame; itx++) {
- n10 = saatran(aa1,&f_str->aa1x[last_n1],n1,itx);
-/*
- for (i=0; i<n10; i++) {
- fprintf(stderr,"%c",pst.sq[aa10[last_n1+i]]);
- if ((i%60)==59) fprintf(stderr,"\n");
- }
- fprintf(stderr,"\n");
-*/
- last_n1 += n10+1;
- }
- n10 = last_n1-1;
-
- /* create aa1y from aa1x */
- for (fs=f_str->aa1x,itemp=0; itemp <3; itemp++,fs++) {
- for (fd= &f_str->aa1y[itemp]; *fs!=EOSEQ; fd += 3, fs++) *fd = *fs;
- *fd=EOSEQ;
- }
-#endif
-}
-
-
-/*
- Alignment: store the operation that align the protein and dna sequence.
- The code of the number in the array is as follows:
- 0: delete of an amino acid.
- 2: frame shift, 2 nucleotides match with an amino acid
- 3: match an amino acid with a codon
- 4: the other type of frame shift
- 5: delete of a codon
-
- The first two elements of the array stores the starting point
- in the protein and dna sequences in the local alignment.
-*/
-
-#include "a_mark.h"
-
-int calcons(const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- int *nc,
- struct a_struct *aln,
- struct a_res_str a_res,
- struct pstruct pst,
- char *seqc0, char *seqc1, char *seqca,
- struct f_struct *f_str)
-{
- int i0, i1, i, j;
- int lenc, not_c, itmp, ngap_p, ngap_d, nfs;
- char *sp0, *sp1, *spa, *sq;
- const unsigned char *ap0, *ap1;
- int *rp, *rpmax;
-
- if (pst.ext_sq_set) {sq = pst.sqx;}
- else {sq = pst.sq;}
-
-
-
-#ifndef TFAST /* FASTX */
- aln->amin1 = aln->smin1 = a_res.min0; /* prot */
- aln->amin0 = aln->smin0 = a_res.min1; /* DNA */
-
- ap0 = f_str->aa0y; /* translated DNA */
- ap1 = aa1; /* protein */
-
- sp0 = seqc0;
- sp1 = seqc1;
-#else /* TFASTX */
- aln->amin0 = aln->smin0 = a_res.min0; /* DNA */
- aln->amin1 = aln->smin1 = a_res.min1; /* prot */
-
- ap1 = aa0; /* protein */
- ap0 = f_str->aa1y; /* translated DNA */
-
- sp1 = seqc0;
- sp0 = seqc1;
-#endif
-
- rp = a_res.res;
- rpmax = rp+a_res.nres;
-
- spa = seqca;
-
- lenc = not_c = aln->nident = aln->nsim = ngap_p = ngap_d = nfs= 0;
- i0 = a_res.min1;
- i1 = a_res.min0;
-
- while (rp < rpmax) {
- /* fprintf(stderr,"%d %d %d (%c) %d (%c)\n"
- ,(int)(rp-res),*rp,i0,sq[ap0[i0]],i1,sq[ap1[i1]]);
- */
- switch (*rp++) {
- case 0: /* aa insertion */
- *sp0++ = '-';
- *sp1++ = sq[ap1[i1++]];
- *spa++ = M_DEL;
- lenc++;
- ngap_d++;
- break;
- case 2: /* -1 frameshift */
- nfs++;
- *sp0++ = '/';
- i0 -= 1;
- *sp1++ = '-';
- *spa++ = M_DEL;
- not_c++;
-
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 3: /* codon/aa match */
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 4: /* +1 frameshift */
- nfs++;
- *sp0++ = '\\';
- i0 += 1;
- *sp1++ = '-';
- *spa++ = M_DEL;
- not_c++;
-
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 5: /* codon insertion */
- *sp0++ = sq[ap0[i0]];
- i0 += 3;
- *sp1++ = '-';
- *spa++ = M_DEL;
- lenc++;
- ngap_p++;
- break;
- }
- }
- *spa = '\0';
-
-#ifndef TFAST /* FASTX */
- aln->amax0 = i0;
- aln->amax1 = i1;
- aln->ngap_q = ngap_d;
- aln->ngap_l = ngap_p;
-#else
- aln->amax1 = i0;
- aln->amax0 = i1;
- aln->amin1 = aln->smin1;
- aln->amin0 = aln->smin0;
- aln->ngap_q = ngap_p;
- aln->ngap_l = ngap_d;
-#endif
- aln->nfs = nfs;
-
- if (lenc < 0) lenc = 1;
- *nc = lenc;
-/* now we have the middle, get the right end */
- return lenc+not_c;
-}
-
-int calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
- const unsigned char *aa1, int n1,
- int *nc,
- struct a_struct *aln,
- struct a_res_str a_res,
- struct pstruct pst,
- char *seqc0, char *seqc0a, char *seqc1, char *seqca,
- char *ann_arr, struct f_struct *f_str)
-{
- int i0, i1, i, j;
- int lenc, not_c, itmp, ngap_p, ngap_d, nfs;
- char *sp0, *sp0a, *sp1, *spa, *sq;
- const unsigned char *ap0, *ap1;
- int *rp, *rpmax;
-
- if (pst.ext_sq_set) {sq = pst.sqx;}
- else {sq = pst.sq;}
-
-#ifndef TFAST /* FASTX */
- aln->amin1 = aln->smin1 = a_res.min0; /* prot */
- aln->amin0 = aln->smin0 = a_res.min1; /* DNA */
-
- ap0 = f_str->aa0y; /* translated DNA */
- ap1 = aa1; /* protein */
-#else /* TFASTX */
- aln->amin0 = aln->smin0 = a_res.min0; /* DNA */
- aln->amin1 = aln->smin1 = a_res.min1; /* prot */
-
- ap1 = aa0;
- ap0 = f_str->aa1y;
-#endif
-
- rp = a_res.res;
- rpmax = &a_res.res[a_res.nres];
-
-#ifndef TFAST
- sp0 = seqc0;
- sp1 = seqc1;
-#else
- sp1 = seqc0;
- sp0 = seqc1;
-#endif
- spa = seqca;
- sp0a = seqc0a;
-
- lenc = not_c = aln->nident = aln->nsim = ngap_p = ngap_d = nfs= 0;
- i0 = a_res.min1;
- i1 = a_res.min0;
-
- while (rp < rpmax) {
- /* fprintf(stderr,"%d %d %d (%c) %d (%c)\n"
- ,(int)(rp-res),*rp,i0,sq[ap0[i0]],i1,sq[ap1[i1]]);
- */
- switch (*rp++) {
- case 0: /* aa insertion */
- *sp0++ = '-';
- *sp1++ = sq[ap1[i1++]];
- *spa++ = M_DEL;
- *sp0a++ = ' ';
- lenc++;
- ngap_d++;
- break;
- case 2: /* -1 frameshift */
- nfs++;
- *sp0++ = '/';
- i0 -= 1;
- *sp1++ = '-';
- *spa++ = M_DEL;
- *sp0a++ = ' ';
- not_c++;
-
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
-#ifndef TFAST
- *sp0a++ = ' ';
-#else
- *sp0a++ = ann_arr[aa0a[i1]];
-#endif
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 3: /* codon/aa match */
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
-#ifndef TFAST
- *sp0a++ = ' ';
-#else
- *sp0a++ = ann_arr[aa0a[i1]];
-#endif
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 4: /* +1 frameshift */
- nfs++;
- *sp0a++ = ' ';
- *sp0++ = '\\';
- i0 += 1;
- *sp1++ = '-';
- *spa++ = M_DEL;
- not_c++;
-
- if ((itmp=pst.pam2[0][ap0[i0]][ap1[i1]])<0) { *spa = M_NEG; }
- else if (itmp == 0) { *spa = M_ZERO;}
- else {*spa = M_POS;}
- if (*spa == M_POS || *spa == M_ZERO) { aln->nsim++;}
-
- *sp0 = sq[ap0[i0]];
- i0 += 3;
- *sp1 = sq[ap1[i1++]];
- if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
- sp0++; sp1++; spa++;
- lenc++;
- break;
- case 5: /* codon insertion */
- *sp0a++ = ' ';
- *sp0++ = sq[ap0[i0]];
- i0 += 3;
- *sp1++ = '-';
- *spa++ = M_DEL;
- lenc++;
- ngap_p++;
- break;
- }
- }
- *sp0a = *spa = '\0';
-
-#ifndef TFAST
- aln->amax0 = i0;
- aln->amax1 = i1;
- aln->ngap_q = ngap_d;
- aln->ngap_l = ngap_p;
-#else
- aln->amax1 = i0;
- aln->amax0 = i1;
- aln->ngap_q = ngap_p;
- aln->ngap_l = ngap_d;
-#endif
- aln->nfs = nfs;
-
- if (lenc < 0) lenc = 1;
- *nc = lenc;
-/* now we have the middle, get the right end */
- return lenc+not_c;
-}
-
-/* build an array of match/ins/del - length strings */
-int calc_code(const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- struct a_struct *aln,
- struct a_res_str a_res,
- struct pstruct pst,
- char *al_str, int al_str_n, struct f_struct *f_str)
-{
- int i0, i1, i, j;
- int lenc, not_c, itmp, ngap_p, ngap_d, nfs;
- char op_char[10];
- int op, op_cnt;
- char sp0, sp1, *sq;
- const unsigned char *ap0, *ap1;
- int *rp, *rpmax;
-
- if (pst.ext_sq_set) {sq = pst.sqx;}
- else {sq = pst.sq;}
-
-
-#ifndef TFAST /* FASTX */
- strncpy(op_char,"- /=\\+*",sizeof(op_char));
- aln->amin1 = aln->smin1 = a_res.min0; /* prot */
- aln->amin0 = aln->smin0 = a_res.min1; /* DNA */
-
- ap0 = f_str->aa0y;
- ap1 = aa1;
-#else /* TFASTX */
- strncpy(op_char,"+ /=\\-*",sizeof(op_char));
- aln->amin0 = aln->smin0 = a_res.min0; /* DNA */
- aln->amin1 = aln->smin1 = a_res.min1; /* prot */
-
- ap1 = aa0;
- ap0 = f_str->aa1y;
-#endif
-
- rp = a_res.res;
- rpmax = &a_res.res[a_res.nres];
-
- op_cnt = lenc = not_c = aln->nident = aln->nsim = ngap_p = ngap_d = nfs = 0;
- op = 3; /* code for a match - all alignments start with a match */
-
- i0 = a_res.min1;
- i1 = a_res.min0;
-
- while (rp < rpmax) {
- switch (*rp++) {
- case 0: /* aa insertion */
- if (op == 0) op_cnt++;
- else {
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 0; op_cnt = 1;
- }
- i1++;
- lenc++;
- ngap_d++;
- break;
- case 2: /* -1 frameshift */
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
-
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 2; op_cnt = 1;
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 3; op_cnt = 1;
- nfs++;
- i0 -= 1;
- not_c++;
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
- lenc++;
- break;
- case 3: /* codon/aa match */
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
-
- if (op == 3 || op == 6) {
- if (sp0 != '*' && sp1 != '*') {
- if (op == 6 ) {
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op_cnt = 1; op = 3;
- }
- else {op_cnt++;}
- }
- else {
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op_cnt = 1; op = 6;
- }
- }
- else {
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- if (op == 2 || op == 4) op_cnt = 2;
- else op_cnt = 1;
- op = 3;
- }
- lenc++;
- break;
- case 4: /* +1 frameshift */
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 4; op_cnt = 1;
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 3; op_cnt = 1;
-
- nfs++;
- i0 += 1;
- not_c++;
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
- lenc++;
- break;
- case 5: /* codon insertion */
- if (op == 5) op_cnt++;
- else {
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
- op = 5; op_cnt = 1;
- }
- i0 += 3;
- lenc++;
- ngap_p++;
- break;
- }
- }
- update_code(al_str, al_str_n-strlen(al_str),op, op_cnt,op_char);
-
-#ifndef TFAST
- aln->amax0 = i0;
- aln->amax1 = i1;
- aln->ngap_q = ngap_d;
- aln->ngap_l = ngap_p;
-#else
- aln->amax1 = i0;
- aln->amax0 = i1;
- aln->ngap_q = ngap_p;
- aln->ngap_l = ngap_d;
-#endif
- aln->nfs = nfs;
-
- if (lenc < 0) lenc = 1;
-
- return lenc;
-}
-
-static void
-update_code(char *al_str, int al_str_max, int op, int op_cnt, char *op_char) {
-
- char tmp_cnt[20];
-
- sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
- strncat(al_str,tmp_cnt,al_str_max-1);
- al_str[al_str_max-1]='\0';
-}
-
-int calc_id(const unsigned char *aa0, int n0,
- const unsigned char *aa1, int n1,
- struct a_struct *aln,
- struct a_res_str a_res,
- struct pstruct pst,
- struct f_struct *f_str)
-{
- int i0, i1, i, j;
- int lenc, not_c, itmp, ngap_p, ngap_d, nfs;
- char sp0, sp1, *sq;
- const unsigned char *ap0, *ap1;
- int *rp, *rpmax;
-
- if (pst.ext_sq_set) {sq = pst.sqx;}
- else {sq = pst.sq;}
-
-
-#ifndef TFAST /* FASTX */
- aln->amin1 = aln->smin1 = a_res.min0; /* prot */
- aln->amin0 = aln->smin0 = a_res.min1; /* DNA */
-
- ap0 = f_str->aa0y;
- ap1 = aa1;
-#else /* TFASTX */
- aln->amin0 = aln->smin0 = a_res.min0; /* DNA */
- aln->amin1 = aln->smin1 = a_res.min1; /* prot */
-
- ap1 = aa0;
- ap0 = f_str->aa1y;
-#endif
-
- rp = a_res.res;
- rpmax = &a_res.res[a_res.nres];
-
- lenc = not_c = aln->nident = aln->nsim = ngap_p = ngap_d = nfs = 0;
- i0 = a_res.min1;
- i1 = a_res.min0;
-
- while (rp < rpmax) {
- /* fprintf(stderr,"%d %d %d (%c) %d (%c)\n"
- ,(int)(rp-res),*rp,i0,sq[ap0[i0]],i1,sq[ap1[i1]]);
- */
- switch (*rp++) {
- case 0: /* aa insertion */
- i1++;
- lenc++;
- ngap_d++;
- break;
- case 2: /* -1 frameshift */
- nfs++;
- i0 -= 1;
- not_c++;
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
- lenc++;
- break;
- case 3: /* codon/aa match */
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
- lenc++;
- break;
- case 4: /* +1 frameshift */
- nfs++;
- i0 += 1;
- not_c++;
- if (pst.pam2[0][ap0[i0]][ap1[i1]]>=0) { aln->nsim++;}
- sp0 = sq[ap0[i0]];
- i0 += 3;
- sp1 = sq[ap1[i1++]];
- if (toupper(sp0) == toupper(sp1)) aln->nident++;
- lenc++;
- break;
- case 5: /* codon insertion */
- i0 += 3;
- lenc++;
- ngap_p++;
- break;
- }
- }
-
-#ifndef TFAST
- aln->amax0 = i0;
- aln->amax1 = i1;
- aln->ngap_q = ngap_d;
- aln->ngap_l = ngap_p;
-#else
- aln->amax1 = i0;
- aln->amax0 = i1;
- aln->ngap_q = ngap_p;
- aln->ngap_l = ngap_d;
-#endif
- aln->nfs = nfs;
-
- if (lenc < 0) lenc = 1;
-/* now we have the middle, get the right end */
- return lenc;
-}
-
-#ifdef PCOMPLIB
-#include "p_mw.h"
-void
-update_params(struct qmng_str *qm_msg, struct pstruct *ppst)
-{
- ppst->n0 = qm_msg->n0;
-}
-#endif
git://source.jalview.org / jabaws.git / blobdiff