--- /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