2 /* copyright (c) 1998, 1999 William R. Pearson and the U. of Virginia */
4 /* $Name: fa_34_26_5 $ - $Id: dropnfa.c,v 1.81 2007/04/26 18:37:19 wrp Exp $ */
6 /* 18-Sep-2006 - removed global variables for alignment from nw_align
9 /* 18-Oct-2005 - converted to use a_res and aln for alignment coordinates */
11 /* 14-May-2003 - modified to return alignment start at 0, rather than
12 1, for begin:end alignments
16 implements the fasta algorithm, see:
18 W. R. Pearson, D. J. Lipman (1988) "Improved tools for biological
19 sequence comparison" Proc. Natl. Acad. Sci. USA 85:2444-2448
21 This version uses Smith-Waterman for final protein alignments
23 W. R. Pearson (1996) "Effective protein sequence comparison"
24 Methods Enzymol. 266:227-258
27 26-April-2001 - -DGAP_OPEN redefines -f, as gap open penalty
29 4-Nov-2001 - modify spam() while(1).
41 /* this must be consistent with upam.h */
43 #define NMAP MAXHASH+1
45 /* globals for fasta */
53 static char *verstr="3.5 Sept 2006";
55 static char *verstr="3.5an0 Sept 2006";
58 extern void w_abort(char *, char *);
59 int shscore(const unsigned char *aa0, int n0, int **pam2);
60 extern void init_karlin(const unsigned char *aa0, int n0, struct pstruct *ppst,
61 double *aa0_f, double **kp);
62 extern void init_karlin_a(struct pstruct *, double *, double **);
63 extern int do_karlin(const unsigned char *, int n1, int **,
64 struct pstruct *, double *, double *,
66 extern void aancpy(char *to, char *from, int count, struct pstruct pst);
67 char *ckalloc(size_t);
70 extern int aatran(const unsigned char *ntseq, unsigned char *aaseq, int maxs, int frame);
76 #include "drop_func.h"
78 struct swstr { int H, E;};
81 dmatch (const unsigned char *aa0, int n0,
82 const unsigned char *aa1, int n1,
84 int **pam2, int gdelval, int ggapval,
85 struct f_struct *f_str);
87 /* initialize for fasta */
90 init_work (unsigned char *aa0, int n0,
92 struct f_struct **f_arg)
99 struct f_struct *f_str;
100 /* these used to be globals, but do not need to be */
101 int ktup; /* word size examined */
102 int fact; /* factor used to scale ktup match value */
103 int kt1; /* ktup-1 */
104 int lkt; /* last ktup - initiall kt1, but can be increased
107 int maxn0; /* used in band alignment */
108 int *pwaa; /* pam[aa0[]] profile */
114 if (ppst->ext_sq_set) {
115 nsq = ppst->nsqx; ip = 1;
119 nsq = ppst->nsq; ip = 0;
123 f_str = (struct f_struct *)calloc(1,sizeof(struct f_struct));
126 if((ppst->zsflag%10) == 6) {
128 init_karlin(aa0, n0, ppst, &f_str->aa0_f[0], &f_str->kar_p);
132 btemp = 2 * ppst->param_u.fa.bestoff / 3 +
133 n0 / ppst->param_u.fa.bestscale +
134 ppst->param_u.fa.bkfact *
135 (ppst->param_u.fa.bktup - ppst->param_u.fa.ktup);
138 btemp = (btemp*ppst->pam_h)/5; /* normalize to standard +5/-4 */
140 btemp = min (btemp, ppst->param_u.fa.bestmax);
141 if (btemp > 3 * n0) btemp = 3 * shscore(aa0,n0,ppst->pam2[0]) / 5;
143 ppst->param_u.fa.cgap = btemp + ppst->param_u.fa.bestoff / 3;
145 if (ppst->param_u.fa.optcut_set != 1)
147 ppst->param_u.fa.optcut = btemp;
149 ppst->param_u.fa.optcut = (btemp*3)/2;
152 #ifndef OLD_FASTA_GAP
153 ppst->param_u.fa.pgap = ppst->gdelval + 2*ppst->ggapval;
155 ppst->param_u.fa.pgap = ppst->gdelval + ppst->ggapval;
157 pamfact = ppst->param_u.fa.pamfact;
158 ktup = ppst->param_u.fa.ktup;
159 fact = ppst->param_u.fa.scfact * ktup;
161 if (pamfact == -1) pamfact = 0;
162 else if (pamfact == -2) pamfact = 1;
164 for (i0 = 1, mhv = -1; i0 <= ppst->nsq; i0++)
165 if (hsq[i0] < NMAP && hsq[i0] > mhv) mhv = hsq[i0];
168 fprintf (stderr, " maximum hsq <=0 %d\n", mhv);
172 for (f_str->kshft = 0; mhv > 0; mhv /= 2) f_str->kshft++;
177 for (i0 = 0; i0 < ktup; i0++) hv = hv << f_str->kshft;
179 f_str->hmask = (hmax >> f_str->kshft) - 1;
181 if ((f_str->harr = (int *) calloc (hmax, sizeof (int))) == NULL) {
182 fprintf (stderr, " cannot allocate hash array: hmax: %d hmask: %d\n",
187 if ((f_str->pamh1 = (int *) calloc (nsq+1, sizeof (int))) == NULL) {
188 fprintf (stderr, " cannot allocate pamh1 array nsq=%d\n",nsq);
192 if ((f_str->pamh2 = (int *) calloc (hmax, sizeof (int))) == NULL) {
193 fprintf (stderr, " cannot allocate pamh2 array hmax=%d\n",hmax);
197 if ((f_str->link = (int *) calloc (n0, sizeof (int))) == NULL) {
198 fprintf (stderr, " cannot allocate hash link array n0=%d",n0);
202 for (i0 = 0; i0 < hmax; i0++) f_str->harr[i0] = -1;
203 for (i0 = 0; i0 < n0; i0++) f_str->link[i0] = -1;
205 /* encode the aa0 array */
208 /* restart hv, phv calculation */
209 for (i0 = 0; i0 < min(lkt,n0); i0++) {
210 if (hsq[aa0[i0]] >= NMAP) {hv=phv=0; lkt = i0+ ktup; continue;}
211 hv = (hv << f_str->kshft) + hsq[aa0[i0]];
212 phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
215 for (; i0 < n0; i0++) {
216 if (hsq[aa0[i0]] >= NMAP) {
218 /* restart hv, phv calculation */
219 for (lkt = i0+kt1; (i0 < lkt || hsq[aa0[i0]]>=NMAP) && i0<n0; i0++) {
220 if (hsq[aa0[i0]] >= NMAP) {
225 hv = (hv << f_str->kshft) + hsq[aa0[i0]];
226 phv += ppst->pam2[ip][aa0[i0]][aa0[i0]]*ktup;
230 hv = ((hv & f_str->hmask) << f_str->kshft) + hsq[aa0[i0]];
231 f_str->link[i0] = f_str->harr[hv];
232 f_str->harr[hv] = i0;
234 f_str->pamh2[hv] = (phv += ppst->pam2[ip][aa0[i0]][aa0[i0]] * ktup);
235 /* this check should always be true, but just in case */
236 if (hsq[aa0[i0-kt1]]<NMAP)
237 phv -= ppst->pam2[ip][aa0[i0 - kt1]][aa0[i0 - kt1]] * ktup;
239 else f_str->pamh2[hv] = fact * ktup;
242 /* this has been modified from 0..<ppst->nsq to 1..<=ppst->nsq because the
243 pam2[0][0] is now undefined for consistency with blast
247 for (i0 = 1; i0 <= nsq; i0++)
248 f_str->pamh1[i0] = ppst->pam2[ip][i0][i0] * ktup;
250 for (i0 = 1; i0 <= nsq; i0++)
251 f_str->pamh1[i0] = fact;
256 if ((f_str->diag = (struct dstruct *) calloc ((size_t)MAXDIAG,
257 sizeof (struct dstruct)))==NULL) {
258 fprintf (stderr," cannot allocate diagonal arrays: %lu\n",
259 MAXDIAG *sizeof (struct dstruct));
263 if ((f_str->diag = (struct dstruct *) calloc ((size_t)n0,
264 sizeof (struct dstruct)))==NULL) {
265 fprintf (stderr," cannot allocate diagonal arrays: %ld\n",
266 (long)n0*sizeof (struct dstruct));
273 if ((f_str->aa1x =(unsigned char *)calloc((size_t)ppst->maxlen+2,
274 sizeof(unsigned char)))
276 fprintf (stderr, "cannot allocate aa1x array %d\n", ppst->maxlen+2);
282 f_str->bss = (struct bdstr *) calloc((size_t)ppst->param_u.fa.optwid*2+4,
283 sizeof(struct bdstr));
286 /* allocate space for the scoring arrays */
288 if ((ss = (struct swstr *) calloc (maxn0, sizeof (struct swstr)))
290 fprintf (stderr, "cannot allocate ss array %3d\n", n0);
296 /* initialize the "variable" pam array */
298 if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
299 fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
304 for (i=0; i<=nsq; i++) {
305 for (j=0;j<n0; j++) {
306 *pwaa = ppst->pam2[ip][aa0[j]][i];
312 /* initialize the "conventional" pam array used for alignments */
314 if ((waa= (int *)calloc ((size_t)(nsq+1)*n0,sizeof(int))) == NULL) {
315 fprintf(stderr,"cannot allocate waa struct %3d\n",nsq*n0);
320 for (i=0; i<=nsq; i++) {
321 for (j=0;j<n0; j++) {
322 *pwaa = ppst->pam2[0][aa0[j]][i];
328 f_str->max_res = max(3*n0/2,MIN_RES);
330 /* now we need alignment storage - get it */
331 if ((f_str->res = (int *)calloc((size_t)f_str->max_res,sizeof(int)))==NULL) {
332 fprintf(stderr,"cannot allocate alignment results array %d\n",f_str->max_res);
340 /* pstring1 is a message to the manager, currently 512 */
341 /* pstring2 is the same information, but in a markx==10 format */
343 get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
346 char *pg_str="FASTA";
348 char *pg_str="TFASTA";
351 if (!pstr->param_u.fa.optflag)
353 sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, gap-pen: %d/%d, width: %3d",
355 sprintf (pstring1, "%s (%s) function [%s matrix, (%d:%d)%s] ktup: %d\n join: %d, open/ext: %d/%d, width: %3d",
357 pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
358 (pstr->ext_sq_set) ? "xS":"\0",
359 pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
360 pstr->gdelval, pstr->ggapval, pstr->param_u.fa.optwid);
363 sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, gap-pen: %d/%d, width: %3d",
365 sprintf (pstring1, "%s (%s) function [optimized, %s matrix (%d:%d)%s] ktup: %d\n join: %d, opt: %d, open/ext: %d/%d, width: %3d",
367 pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l,
368 (pstr->ext_sq_set) ? "xS":"\0",
369 pstr->param_u.fa.ktup, pstr->param_u.fa.cgap,
370 pstr->param_u.fa.optcut, pstr->gdelval, pstr->ggapval,
371 pstr->param_u.fa.optwid);
372 if (pstr->param_u.fa.iniflag) strcat(pstring1," init1");
374 if (pstr->zsflag==0) strcat(pstring1," not-scaled");
375 else if (pstr->zsflag==1) strcat(pstring1," reg.-scaled");
378 if (pstring2 != NULL) {
380 sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
381 ; pg_gap-pen: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
383 sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
384 ; pg_open-ext: %d %d\n; pg_ktup: %d\n; pg_optcut: %d\n; pg_cgap: %d\n",
386 pg_str,verstr,pstr->pamfile, pstr->pam_h,pstr->pam_l, pstr->gdelval,
387 pstr->ggapval,pstr->param_u.fa.ktup,pstr->param_u.fa.optcut,
388 pstr->param_u.fa.cgap);
393 close_work (const unsigned char *aa0, int n0,
394 struct pstruct *ppst,
395 struct f_struct **f_arg)
397 struct f_struct *f_str;
404 if (f_str->kar_p!=NULL) free(f_str->kar_p);
425 void savemax (struct dstruct *, int, struct f_struct *);
427 void savemax (struct dstruct *, struct f_struct *);
430 int spam (const unsigned char *, const unsigned char *, struct savestr *,
431 int **, int, int, int);
432 int sconn(struct savestr **, int nsave, int cgap, int pgap, int noff);
433 void kpsort(struct savestr **, int);
436 ALIGN(const unsigned char *, const unsigned char *, int, int,
437 int **, int, int, int *, int *, struct f_struct *);
440 LOCAL_ALIGN(const unsigned char *, const unsigned char *,
442 int **, int, int, int *, int *, int *, int *, int,
446 B_ALIGN(const unsigned char *A, const unsigned char *B, int M,
447 int N, int low, int up, int **W, int G, int H, int *S,
448 int *nS, int MW, int MX, struct bdstr *bss);
451 do_fasta (const unsigned char *aa0, int n0,
452 const unsigned char *aa1, int n1,
453 struct pstruct *ppst, struct f_struct *f_str,
454 struct rstruct *rst, int *hoff)
456 int nd; /* diagonal array size */
459 register struct dstruct *dptr;
463 register struct dstruct *diagp;
469 struct dstruct *dpmax;
472 struct savestr *vmptr;
474 int xdrop, do_extend;
475 int ktup, kt1, lkt, *hsq, ip;
477 if (ppst->ext_sq_set) {
486 xdrop = -ppst->pam_l;
487 /* do extended alignment in spam iff protein or short sequences */
488 do_extend = !ppst->nt_align || (n0 < 50) || (n1 < 50);
490 ktup = ppst->param_u.fa.ktup;
494 rst->score[0] = rst->score[1] = rst->score[2] = 0;
498 if (n0+n1+1 >= MAXDIAG) {
499 fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
500 rst->score[0] = rst->score[1] = rst->score[2] = -1;
510 dpmax = &f_str->diag[nd];
511 for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;)
518 for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
520 f_str->lowmax = f_str->vmax;
526 for (lpos = 0; (lpos < lkt || hsq[aa1[lpos]]>=NMAP) && lpos <n1; lpos++) {
527 /* restart lhval calculation */
528 if (hsq[aa1[lpos]]>=NMAP) {
529 lhval = 0; lkt = lpos + ktup;
532 lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
537 diagp = &f_str->diag[noff + lkt];
538 for (; lpos < n1; lpos++, diagp++) {
539 if (hsq[aa1[lpos]]>=NMAP) {
541 while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
542 if (lpos >= n1) break;
545 lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
546 for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
547 if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
549 lposn0 = noff + lpos;
550 for (; lpos < n1; lpos++, lposn0++) {
551 if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; goto loopl;}
553 if (hsq[aa1[lpos]]>=NMAP) {
555 while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; lposn0++;}
558 lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
559 for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
560 dpos = lposn0 - tpos;
561 if ((tscor = (dptr = &f_str->diag[dpos % nd])->stop) >= 0) {
564 if ((tscor -= lpos) <= 0) {
566 if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 && f_str->lowscor < scor)
568 savemax (dptr, dpos, f_str);
570 savemax (dptr, f_str);
572 if ((tscor += scor) >= kfact) {
578 dptr->start = (dptr->stop = lpos) - kt1;
582 dptr->score += f_str->pamh1[aa0[tpos]];
587 dptr->score = f_str->pamh2[lhval];
588 dptr->start = (dptr->stop = lpos) - kt1;
593 /* reinitialize diag structure */
595 if ((dptr = &f_str->diag[lpos % nd])->score > f_str->lowscor)
596 savemax (dptr, lpos, f_str);
604 for (tpos = 0, dpos = noff + n1 - 1; tpos < n0; tpos++, dpos--) {
605 if ((dptr = &f_str->diag[dpos % nd])->score > f_str->lowscor)
606 savemax (dptr, dpos, f_str);
609 for (dptr = f_str->diag; dptr < dpmax;) {
610 if (dptr->score > f_str->lowscor) savemax (dptr, f_str);
619 at this point all of the elements of aa1[lpos]
620 have been searched for elements of aa0[tpos]
621 with the results in diag[dpos]
623 for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
625 fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
626 noff+vmptr->start-vmptr->dp,
627 noff+vmptr->stop-vmptr->dp,
628 vmptr->start,vmptr->stop,
629 vmptr->dp,vmptr->score);
632 if (vmptr->score > 0) {
633 vmptr->score = spam (aa0, aa1, vmptr, ppst->pam2[ip], xdrop,
635 f_str->vptr[nsave++] = vmptr;
640 rst->score[0] = rst->score[1] = rst->score[2] = 0;
645 fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,noff);
646 for (ib=0; ib<nsave; ib++) {
647 fprintf(stderr,"0: %4d-%4d 1: %4d-%4d dp: %d score: %d\n",
648 noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
649 noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
650 f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
651 f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
653 fprintf(stderr,"---\n");
656 scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap,
657 ppst->param_u.fa.pgap, noff);
659 for (vmptr=f_str->vptr[0],ib=1; ib<nsave; ib++)
660 if (f_str->vptr[ib]->score > vmptr->score) vmptr=f_str->vptr[ib];
662 /* kssort (f_str->vptr, nsave); */
664 rst->score[1] = vmptr->score;
665 rst->score[0] = max (scor, vmptr->score);
666 rst->score[2] = rst->score[0]; /* initn */
668 if (ppst->param_u.fa.optflag) {
669 if (rst->score[0] > ppst->param_u.fa.optcut)
670 rst->score[2] = dmatch (aa0, n0, aa1, n1, *hoff=noff - vmptr->dp,
671 ppst->param_u.fa.optwid, ppst->pam2[ip],
672 ppst->gdelval,ppst->ggapval,f_str);
676 void do_work (const unsigned char *aa0, int n0,
677 const unsigned char *aa1, int n1,
679 struct pstruct *ppst, struct f_struct *f_str,
680 int qr_flg, struct rstruct *rst)
686 rst->score[0] = rst->score[1] = rst->score[2] = 0;
688 rst->segnum = rst->seglen = 1;
690 if (n1 < ppst->param_u.fa.ktup) return;
693 n10=aatran(aa1,f_str->aa1x,n1,frame);
694 do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
696 do_fasta (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff);
700 if((ppst->zsflag%10) == 6 &&
701 (do_karlin(aa1, n1, ppst->pam2[0], ppst,f_str->aa0_f,
702 f_str->kar_p, &lambda, &H)>0)) {
703 rst->comp = 1.0/lambda;
706 else {rst->comp = rst->H = -1.0;}
708 rst->comp = rst->H = -1.0;
712 void do_opt (const unsigned char *aa0, int n0,
713 const unsigned char *aa1, int n1,
715 struct pstruct *ppst,
716 struct f_struct *f_str,
719 int optflag, tscore, hoff, n10;
721 optflag = ppst->param_u.fa.optflag;
722 ppst->param_u.fa.optflag = 1;
725 n10=aatran(aa1,f_str->aa1x,n1,frame);
726 do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
728 do_fasta(aa0,n0,aa1,n1,ppst,f_str,rst, &hoff);
730 ppst->param_u.fa.optflag = optflag;
735 savemax (dptr, dpos, f_str)
736 register struct dstruct *dptr;
738 struct f_struct *f_str;
740 register struct savestr *vmptr;
745 savemax (dptr, f_str)
746 register struct dstruct *dptr;
747 struct f_struct *f_str;
750 register struct savestr *vmptr;
753 dpos = (int) (dptr - f_str->diag);
757 /* check to see if this is the continuation of a run that is already saved */
759 if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
760 vmptr->start == dptr->start)
762 vmptr->stop = dptr->stop;
763 if ((i = dptr->score) <= vmptr->score)
766 if (vmptr != f_str->lowmax)
771 i = f_str->lowmax->score = dptr->score;
772 f_str->lowmax->dp = dpos;
773 f_str->lowmax->start = dptr->start;
774 f_str->lowmax->stop = dptr->stop;
775 dptr->dmax = f_str->lowmax;
778 for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
779 if (vmptr->score < i)
782 f_str->lowmax = vmptr;
787 int spam (const unsigned char *aa0, const unsigned char *aa1,
788 struct savestr *dmax, int **pam2, int xdrop,
789 int noff, int do_extend)
791 register int lpos, tot;
792 register const unsigned char *aa0p, *aa1p;
797 int start, stop, score;
800 aa1p = &aa1[lpos= dmax->start]; /* get the start of lib seq */
801 aa0p = &aa0[lpos - dmax->dp + noff]; /* start of query */
803 /* also add check in calling routine */
804 if (aa0p < aa0) { return -99; }
806 curv.start = lpos; /* start index in lib seq */
808 tot = curv.score = maxv.score = 0;
810 for (; lpos <= dmax->stop; lpos++) {
811 tot += pam2[*aa0p++][*aa1p++];
812 if (tot > curv.score) { /* update current score */
817 if (curv.score > maxv.score) { /* save score, start, stop */
818 maxv.start = curv.start;
819 maxv.stop = curv.stop;
820 maxv.score = curv.score;
822 tot = curv.score = 0; /* reset running score */
823 curv.start = lpos+1; /* reset start */
824 if(lpos >= dmax->stop) break; /* if the zero is beyond stop, quit */
828 if (curv.score > maxv.score) {
829 maxv.start = curv.start;
830 maxv.stop = curv.stop;
831 maxv.score = curv.score;
836 /* now check to see if the score gets better by extending */
837 if (do_extend && maxv.score > xdrop) {
839 if (maxv.stop == dmax->stop) {
841 drop_thresh = maxv.score - xdrop;
842 aa1p = &aa1[lpos= dmax->stop];
843 aa0p = &aa0[lpos - dmax->dp + noff];
844 while (tot > drop_thresh ) {
846 tot += pam2[*(++aa0p)][*(++aa1p)];
847 if (tot > maxv.score) {
850 drop_thresh = tot - xdrop;
855 /* scan backwards now */
857 if (maxv.start == dmax->start) {
859 drop_thresh = maxv.score - xdrop;
860 aa1p = &aa1[lpos= dmax->start];
861 aa0p = &aa0[lpos - dmax->dp + noff];
862 while (tot > drop_thresh) {
864 tot += pam2[*(--aa0p)][*(--aa1p)];
865 if (tot > maxv.score) {
868 drop_thresh = tot - xdrop;
875 /* if (maxv.start != dmax->start || maxv.stop != dmax->stop)
876 printf(" new region: %3d %3d %3d %3d\n",maxv.start,
877 dmax->start,maxv.stop,dmax->stop);
879 dmax->start = maxv.start;
880 dmax->stop = maxv.stop;
885 int sconn (struct savestr **v, int n, int cgap, int pgap, int noff)
893 } *start, *sl, *sj, *so, sarr[MAXSAV];
894 int lstart, tstart, plstop, ptstop;
896 /* sort the score left to right in lib pos */
902 /* for the remaining runs, see if they fit */
904 for (i = 0, si = 0; i < n; i++)
907 /* if the score is less than the gap penalty, it never helps */
908 if (v[i]->score < cgap)
910 lstart = v[i]->start;
911 tstart = lstart - v[i]->dp + noff;
913 /* put the run in the group */
915 sarr[si].score = v[i]->score;
916 sarr[si].next = NULL;
918 /* if it fits, then increase the score */
919 for (sl = start; sl != NULL; sl = sl->next)
921 plstop = sl->vp->stop;
922 ptstop = plstop - sl->vp->dp + noff;
923 if (plstop < lstart && ptstop < tstart)
925 sarr[si].score = sl->score + v[i]->score + pgap;
930 /* now recalculate where the score fits */
934 for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
935 if (sarr[si].score > sj->score) {
937 if (so != NULL) so->next = &sarr[si];
938 else start = &sarr[si];
947 return (start->score);
960 for (gap = n / 2; gap > 0; gap /= 2)
961 for (i = gap; i < n; i++)
962 for (j = i - gap; j >= 0; j -= gap)
964 if (v[j]->score >= v[j + gap]->score)
980 for (gap = n / 2; gap > 0; gap /= 2)
981 for (i = gap; i < n; i++)
982 for (j = i - gap; j >= 0; j -= gap)
984 if (v[j]->start <= v[j + gap]->start)
992 static int dmatch (const unsigned char *aa0, int n0,
993 const unsigned char *aa1, int n1,
994 int hoff, int window,
995 int **pam2, int gdelval, int ggapval,
996 struct f_struct *f_str)
1000 window = min (n1, window);
1001 /* hoff is the offset found from aa1 to seq 2 by hmatch */
1003 low = -window/2-hoff;
1006 return FLOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
1008 #ifdef OLD_FASTA_GAP
1013 -ggapval,window,f_str);
1017 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1019 To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,MW).
1020 The parameters are explained as follows:
1021 A, B : two sequences to be aligned
1022 M : the length of sequence A
1023 N : the length of sequence B
1024 L : lower bound of the band
1025 U : upper bound of the band
1026 W : scoring table for matches and mismatches
1027 G : gap-opening penalty
1028 H : gap-extension penalty
1029 MW : maximum window size
1034 #define MININT -9999999
1037 FLOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
1038 int M, int N, int low, int up,
1039 int **W, int G,int H, int MW,
1040 struct f_struct *f_str)
1043 register struct bdstr *bssp;
1057 if (N <= 0) return 0;
1059 if (M <= 0) return 0;
1063 fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
1068 if (low > 0) leftd = 1;
1069 else if (up < 0) leftd = band;
1072 si = max(0,-up); /* start index -1 */
1073 ei = min(M,N-low); /* end index */
1075 for (j = leftd+1; j <= rightd; j++) {
1080 bssp[rightd+1].CC = MININT;
1081 bssp[rightd+1].DD = MININT;
1084 bssp[leftd-1].CC = MININT;
1085 bssp[leftd].DD = -G;
1087 for (i = si+1; i <= ei; i++) {
1088 if (i > N-up) rightd--;
1089 if (leftd > 1) leftd--;
1091 if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
1092 if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
1099 if (c > best_score) best_score = c;
1101 for (curd=leftd+1; curd <= rightd; curd++) {
1102 if ((c = c-m) > (e = e-H)) e = c;
1103 if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
1104 c = bssp[curd].CC + wa[B[curd+low-1+i]];
1110 if (c > best_score) best_score = c;
1117 /* ckalloc - allocate space; check for success */
1118 char *ckalloc(size_t amount)
1122 if ((p = malloc( (unsigned) amount)) == NULL)
1123 w_abort("Ran out of memory.","");
1127 /* calculate the 100% identical score */
1129 shscore(const unsigned char *aa0, int n0, int **pam2)
1132 for (i=0,sum=0; i<n0; i++)
1133 sum += pam2[aa0[i]][aa0[i]];
1137 int sw_walign (const unsigned char *aa0, int n0,
1138 const unsigned char *aa1, int n1,
1139 struct pstruct *ppst,
1140 struct f_struct *f_str,
1141 struct a_res_str *a_res
1144 register const unsigned char *aa0p, *aa1p;
1147 register struct swstr *ssj;
1153 int cost, I, J, K, L;
1158 #ifdef OLD_FASTA_GAP
1159 q = -(ppst->gdelval - ppst->ggapval);
1166 /* initialize 0th row */
1167 for (ssj=ss; ssj<ss+n0; ssj++) {
1178 pwaa = waa + (*aa1p++ * n0);
1179 for (ssj = ss, aa0p = aa0; ssj < ss+n0; ssj++) {
1180 if ((h = h - m) > (f = f - r)) f = h;
1181 if ((h = ssj->H - m) > (e = ssj->E - r)) e = h;
1196 } /* done with forward pass */
1197 if (score <= 0) return 0;
1199 /* to get the start point, go backwards */
1202 for (ssj=ss+J; ssj>=ss; ssj--) ssj->H= ssj->E= -1;
1204 for (i=I; i>=0; i--) {
1206 p = (i == I) ? 0 : -1;
1207 ssj = ss+J; /* bug in compiler */
1208 for (aa0p = &aa0[J]; ssj>=ss; ssj--,aa0p--) {
1210 ssj->E=max(ssj->E,ssj->H-q)-r;
1211 h = max(max(ssj->E,f),p+ppst->pam2[0][*aa0p][aa1[i]]);
1218 if (cost >= score) goto found;
1225 /* printf(" %d: L: %3d-%3d/%3d; K: %3d-%3d/%3d\n",score,L,J,n0,K,I,n1); */
1227 /* in the f_str version, the *res array is already allocated at 4*n0/3 */
1229 a_res->max0 = J+1; a_res->min0 = L; a_res->max1 = I+1; a_res->min1 = K;
1231 /* the seq array arguments in this call have been reversed to allow
1232 assymetric scoring matrices - this affects the score decoding,
1233 and allocation of the score row matrix */
1234 ALIGN(&aa0[L-1],&aa1[K-1],J-L+1,I-K+1,ppst->pam2[0],q,r,a_res->res,&a_res->nres,f_str);
1236 /* DISPLAY(&aa1[K-1],&aa0[L-1],I-K+1,J-L+1,res,L,K,ppst->sq); */
1241 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
1243 int *S, int **W, int G, int H, int *nres);
1245 #define gap(k) ((k) <= 0 ? 0 : g+h*(k)) /* k-symbol indel cost */
1247 /* static int *sapp; */ /* Current script append ptr */
1248 /* static int last; */ /* Last script op appended */
1250 /* Append "Delete k" op */
1253 *last = (*sapp)[-1] -= (k); \
1255 *last = (*sapp)[0] = -(k); \
1259 /* Append "Insert k" op */
1262 *last = (*sapp)[-1] += (k); \
1264 *last = (*sapp)[0] = (k); \
1269 #define REP { *last = (*sapp)[0] = 0; (*sapp)++;} /* Append "Replace" op */
1271 /* align(A,B,M,N,tb,te) returns the cost of an optimum conversion between
1272 A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
1273 and appends such a conversion to the current script. */
1276 nw_align(const unsigned char *A, const unsigned char *B,
1278 int tb, int te, int **w, int g, int h,
1279 struct f_struct *f_str,
1280 int **sapp, int *last)
1282 int midi, midj, type; /* Midpoint, type, and cost */
1286 register int c, e, d, s;
1288 struct swstr *f_ss, *r_ss;
1295 /* Boundary cases: M <= 1 or N == 0 */
1307 if (tb < te) tb = te;
1308 midc = (tb-h) - gap(N);
1310 wa = w[A[1]]; /* in the original version of this code, A[]
1311 is the second sequence */
1312 for (j = 1; j <= N; j++) {
1313 c = -gap(j-1) + wa[B[j]] - gap(N-j);
1319 if (midj == 0) { DEL(1) INS(N) }
1321 if (midj > 1) { INS(midj-1) }
1323 if (midj < N) { INS(N-midj) }
1328 /* Divide: Find optimum midpoint (midi,midj) of cost midc */
1330 midi = M/2; /* Forward phase: */
1331 f_ss[0].H = 0; /* Compute H(M/2,k) & E(M/2,k) for all k */
1333 for (j = 1; j <= N; j++) {
1334 f_ss[j].H = t = t-h;
1338 for (i = 1; i <= midi; i++) {
1340 f_ss[0].H = c = t = t-h;
1343 for (j = 1; j <= N; j++) {
1344 if ((c = c - m) > (e = e - h)) e = c;
1345 if ((c = f_ss[j].H - m) > (d = f_ss[j].E - h)) d = c;
1354 f_ss[0].E = f_ss[0].H;
1356 r_ss[N].H = 0; /* Reverse phase: */
1357 t = -g; /* Compute R(M/2,k) & S(M/2,k) for all k */
1358 for (j = N-1; j >= 0; j--)
1359 { r_ss[j].H = t = t-h;
1363 for (i = M-1; i >= midi; i--)
1365 r_ss[N].H = c = t = t-h;
1368 for (j = N-1; j >= 0; j--)
1369 { if ((c = c - m) > (e = e - h)) e = c;
1370 if ((c = r_ss[j].H - m) > (d = r_ss[j].E - h)) d = c;
1379 r_ss[N].E = r_ss[N].H;
1381 midc = f_ss[0].H+r_ss[0].H; /* Find optimal midpoint */
1384 for (j = 0; j <= N; j++)
1385 if ((c = f_ss[j].H + r_ss[j].H) >= midc)
1386 if (c > midc || (f_ss[j].H != f_ss[j].E && r_ss[j].H == r_ss[j].E))
1390 for (j = N; j >= 0; j--)
1391 if ((c = f_ss[j].E + r_ss[j].E + g) > midc)
1398 /* Conquer: recursively around midpoint */
1401 nw_align(A,B,midi,midj,tb,-g,w,g,h,f_str, sapp, last);
1402 nw_align(A+midi,B+midj,M-midi,N-midj,-g,te,w,g,h,f_str,sapp, last);
1405 nw_align(A,B,midi-1,midj,tb,0,w,g,h,f_str, sapp, last);
1407 nw_align(A+midi+1,B+midj,M-midi-1,N-midj,0,te,w,g,h,f_str, sapp, last);
1412 /* Interface and top level of comparator */
1415 ALIGN(const unsigned char *A, const unsigned char *B,
1417 int **W, int G, int H, int *S, int *nS,
1418 struct f_struct *f_str)
1421 struct swstr *f_ss, *r_ss;
1427 if ((f_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
1429 fprintf (stderr, "cannot allocate f_ss array %3d\n", N+2);
1435 if ((r_ss = (struct swstr *) calloc (N+2, sizeof (struct swstr)))
1437 fprintf (stderr, "cannot allocate r_ss array %3d\n", N+2);
1443 c = nw_align(A,B,M,N,-G,-G,W,G,H,f_str,&sapp, &last); /* OK, do it */
1445 ck = CHECK_SCORE(A,B,M,N,S,W,G,H,nS);
1446 if (c != ck) fprintf(stderr,"Check_score error %d != %d\n",c,ck);
1449 free(r_ss); free(f_ss);
1454 /* Alignment display routine */
1456 static char ALINE[51], BLINE[51], CLINE[51];
1458 void DISPLAY(unsigned char *A, unsigned char *B, int M, int N,
1459 int *S, int AP, int BP, char *sq)
1460 { register char *a, *b, *c;
1461 register int i, j, op;
1464 i = j = op = lines = 0;
1470 while (i < M || j < N)
1471 { if (op == 0 && *S == 0)
1475 *c++ = (*a++ == *b++) ? '|' : ' ';
1486 { *a++ = sq[A[++i]];
1492 if (a >= ALINE+50 || (i >= M && j >= N))
1493 { *a = *b = *c = '\0';
1494 printf("\n%5d ",50*lines++);
1495 for (b = ALINE+10; b <= a; b += 10)
1499 printf("\n%5d %s\n %s\n%5d %s\n",ap,ALINE,CLINE,bp,BLINE);
1509 /* CHECK_SCORE - return the score of the alignment stored in S */
1511 static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
1512 int M, int N, int *S, int **w, int g, int h,
1515 register int i, j, op, nc;
1518 score = i = j = op = nc = 0;
1519 while (i < M || j < N) {
1522 score = w[A[++i]][B[++j]] + score;
1524 /* fprintf(stderr,"=%4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
1527 score = score - (g+op*h);
1528 /* fprintf(stderr,">%4d %4d %4d %4d\n",i,j,-(g+op*h),score); */
1532 score = score - (g-op*h);
1533 /* fprintf(stderr,"<%4d %4d %4d %4d\n",i,j,-(g-op*h),score); */
1544 BCHECK_SCORE(const unsigned char *A, const unsigned char *B,
1545 int M, int N, int *S, int **w, int g, int h,
1548 register int i, j, op, nc;
1552 score = i = j = op = nc = 0;
1554 while (i < M || j < N) {
1557 score = w[A[++i]][B[++j]] + score;
1559 /* fprintf(stderr,"op0 %4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
1562 score = score - (g+op*h);
1563 /* fprintf(stderr,"op> %4d %4d %4d %4d %4d\n",i,j,op,-(g+op*h),score); */
1567 score = score - (g-op*h);
1568 /* fprintf(stderr,"op< %4d %4d %4d %4d %4d\n",i,j,op,-(g-op*h),score); */
1578 /* A PACKAGE FOR LOCALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1580 To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,dflag,&SI,&SJ,&EI,&EJ,MW).
1581 The parameters are explained as follows:
1582 A, B : two sequences to be aligned
1583 M : the length of sequence A
1584 N : the length of sequence B
1585 L : lower bound of the band
1586 U : upper bound of the band
1587 W : scoring table for matches and mismatches
1588 G : gap-opening penalty
1589 H : gap-extension penalty
1590 dflag : 0 - no display or backward pass
1591 *SI : starting position of sequence A in the optimal local alignment
1592 *SJ : starting position of sequence B in the optimal local alignment
1593 *EI : ending position of sequence A in the optimal local alignment
1594 *EJ : ending position of sequence B in the optimal local alignment
1595 MW : maximum window size
1598 int bd_walign (const unsigned char *aa0, int n0,
1599 const unsigned char *aa1, int n1,
1600 struct pstruct *ppst,
1601 struct f_struct *f_str, int hoff,
1602 struct a_res_str *a_res)
1605 int min0, min1, max0, max1;
1608 window = min (n1, ppst->param_u.fa.optwid);
1609 /* hoff is the offset found from aa1 to seq 2 by hmatch */
1611 low = -window/2-hoff;
1614 score=LOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
1616 #ifdef OLD_FASTA_GAP
1617 -(ppst->gdelval-ppst->ggapval),
1622 &min0,&min1,&max0,&max1,ppst->param_u.fa.optwid,f_str);
1625 fprintf(stderr,"n0/n1: %d/%d hoff: %d window: %d\n",
1626 n0, n1, hoff, window);
1631 fprintf(stderr," ALIGN: start0: %d start1: %d stop0: %d stop1: %d, bot: %d top: %d, win: %d MX %d\n",
1632 min0-1,min1-1,max0-min0+1,max1-min1+1,low-(min1-min0),up-(min1-min0),
1633 ppst->param_u.fa.optwid,n0);
1636 a_res->min0 = min0-1; a_res->min1 = min1-1;
1637 a_res->max0 = max0; a_res->max1 = max1;
1639 B_ALIGN(aa0-1+min0-1,aa1-1+min1-1,max0-min0+1,max1-min1+1,
1640 low-(min1-min0),up-(min1-min0),
1642 #ifdef OLD_FASTA_GAP
1643 -(ppst->gdelval-ppst->ggapval),
1648 a_res->res,&a_res->nres,ppst->param_u.fa.optwid,n0,f_str->bss);
1654 LOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
1656 int low, int up, int **W, int G,int H,
1657 int *psi, int *psj, int *pei, int *pej, int MW,
1658 struct f_struct *f_str)
1661 register struct bdstr *bssp;
1665 int best_score, starti, startj, endi, endj;
1677 *psi = *psj = *pei = *pej;
1681 *psi = *psj = *pei = *pej;
1686 fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
1691 j = (MW + 2 + 2) * sizeof(struct bdstr);
1693 /* already done by init_work();
1694 if (f_str->bss==NULL) f_str->bss = (struct bdstr *) ckalloc(j);
1697 if (low > 0) leftd = 1;
1698 else if (up < 0) leftd = band;
1704 for (j = leftd+1; j <= rightd; j++) {
1708 bssp[rightd+1].CC = MININT;
1709 bssp[rightd+1].DD = MININT;
1713 bssp[leftd-1].CC = MININT;
1714 bssp[leftd].DD = -G;
1715 for (i = si+1; i <= ei; i++) {
1716 if (i > N-up) rightd--;
1717 if (leftd > 1) leftd--;
1719 if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
1720 if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
1722 if (ib > N) fprintf(stderr,"B[%d] out of range %d\n",ib,N);
1729 if (c > best_score) {
1734 for (curd=leftd+1; curd <= rightd; curd++) {
1735 if ((c = c-m) > (e = e-H)) e = c;
1736 if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
1738 if ((ib=curd+low-1+i) <= 0 || ib > N)
1739 fprintf(stderr,"B[%d]:%d\n",ib,B[ib]);
1741 c = bssp[curd].CC + wa[B[curd+low-1+i]];
1747 if (c > best_score) {
1750 endj = curd+low-1+i;
1755 leftd = max(1,-endi-low+1);
1756 rightd = band-(up-(endj-endi));
1757 bssp[rightd].CC = 0;
1759 for (j = rightd-1; j >= leftd; j--) {
1760 bssp[j].CC = t = t-H;
1763 for (j = rightd+1; j <= band; ++j) bssp[j].CC = MININT;
1764 bssp[leftd-1].CC = bssp[leftd-1].DD = MININT;
1765 bssp[rightd].DD = -G;
1767 for (i = endi; i >= 1; i--) {
1768 if (i+low <= 0) leftd++;
1769 if (rightd < band) rightd++;
1771 if ((c = bssp[rightd-1].CC-m) > (d = bssp[rightd-1].DD-H)) d = c;
1772 if ((ib = rightd+low-1+i) <= N) c = bssp[rightd].CC+wa[B[ib]];
1775 if (ib <= 0) fprintf(stderr,"rB[%d] <1\n",ib);
1779 bssp[rightd].DD = d;
1780 bssp[rightd].CC = c;
1781 if (c == best_score) {
1787 for (curd=rightd-1; curd >= leftd; curd--) {
1788 if ((c = c-m) > (e = e-H)) e = c;
1789 if ((c = bssp[curd-1].CC-m) > (d = bssp[curd-1].DD-H)) d = c;
1792 if ((ib=curd+low-1+i) <= 0 || ib > N)
1793 fprintf(stderr,"i: %d, B[%d]:%d\n",i,ib,B[ib]);
1795 c = bssp[curd].CC + wa[B[curd+low-1+i]];
1800 if (c == best_score) {
1802 startj = curd+low-1+i;
1807 if (flag == 1) break;
1810 if (starti < 0 || starti > M || startj < 0 || startj > N) {
1811 printf("starti=%d, startj=%d\n",starti,startj);
1812 *psi = *psj = *pei = *pej;
1822 /* A PACKAGE FOR GLOBALLY ALIGNING TWO SEQUENCES WITHIN A BAND:
1824 To invoke, call B_ALIGN(A,B,M,N,L,U,W,G,H,S,MW,MX).
1825 The parameters are explained as follows:
1826 A, B : two sequences to be aligned
1827 M : the length of sequence A
1828 N : the length of sequence B
1829 L : lower bound of the band
1830 U : upper bound of the band
1831 W : scoring table for matches and mismatches
1832 G : gap-opening penalty
1833 H : gap-extension penalty
1834 S : script for DISPLAY routine
1835 MW : maximum window size
1836 MX : maximum length sequence M to be aligned
1840 static int *MP[3]; /* save crossing points */
1841 static int *FP; /* forward dividing points */
1842 static char *MT[3]; /* 0: rep, 1: del, 2: ins */
1845 /* bg_align(A,B,M,N,up,low,tb,te) returns the cost of an optimum conversion between
1846 A[1..M] and B[1..N] and appends such a conversion to the current script.
1847 tb(te)= 1 no gap-open penalty if the conversion begins(ends) with a delete.
1848 tb(te)= 2 no gap-open penalty if the conversion begins(ends) with an insert.
1851 bg_align(const unsigned char *A, const unsigned char *B,
1853 int low, int up, int tb, int te,
1854 int **w, int g, int h,
1855 struct bdstr *bss, int **sapp, int *last)
1857 int rmid, k, l, r, v, kt;
1862 int leftd, rightd; /* for CC, DD, CP and DP */
1863 register int curd; /* current index for CC, DD CP and DP */
1865 register int c, d, e;
1866 int t, fr, *wa, ib, m;
1868 /* Boundary cases: M <= 0 , N <= 0, or up-low <= 0 */
1870 if (M > 0) { DEL(M) }
1877 if ((band = up-low+1) <= 1) {
1878 for (i = 1; i <= M; i++) { REP }
1882 /* Divide: Find all crossing points */
1884 /* Initialization */
1888 rmid = low + midd - 1;
1893 for (j = 0; j < midd; j++)
1894 bss[j].CP = bss[j].DP = -1;
1895 for (j = midd; j <= rightd; j++) {
1896 bss[j].CP = bss[j].DP = 0;
1901 MT[0][0] = MT[1][0] = MT[2][0] = 0;
1902 } else if (leftd > midd) {
1904 for (j = 0; j <= midd; j++) {
1905 bss[j].CP = bss[j].DP = fr;
1907 for (j = midd+1; j <= rightd; j++)
1908 bss[j].CP = bss[j].DP = -1;
1912 MT[0][fr] = MT[1][fr] = MT[2][fr] = 0;
1915 for (j = 0; j < midd; j++) {
1916 bss[j].CP = bss[j].DP = 0;
1918 for (j = midd; j <= rightd; j++) {
1919 bss[j].CP = bss[j].DP = 0;
1924 MT[0][0] = MT[1][0] = MT[2][0] = 0;
1930 for (j = leftd+1; j <= rightd; j++) {
1931 bss[j].CC = t = t-h;
1934 bss[rightd+1].CC = MININT;
1935 bss[rightd+1].DD = MININT;
1936 if (tb == 1) bss[leftd].DD = 0;
1937 else bss[leftd].DD = -g;
1938 bss[leftd-1].CC = MININT;
1939 for (i = 1; i <= M; i++) {
1940 if (i > N-up) rightd--;
1941 if (leftd > 1) leftd--;
1943 if ((c = bss[leftd+1].CC-m) > (d = bss[leftd+1].DD-h)) {
1945 bss[leftd].DP = bss[leftd+1].CP;
1946 } else bss[leftd].DP = bss[leftd+1].DP;
1947 if ((ib = leftd+low-1+i) > 0) c = bss[leftd].CC+wa[B[ib]];
1948 if (d > c || ib <= 0) {
1950 bss[leftd].CP = bss[leftd].DP;
1956 if (leftd == midd) bss[leftd].CP = bss[leftd].DP = IP = i;
1957 for (curd=leftd+1; curd <= rightd; curd++) {
1959 if ((c = c-m) > (e = e-h)) {
1961 IP = bss[curd-1].CP;
1962 } /* otherwise, IP is unchanged */
1963 if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
1965 bss[curd].DP = bss[curd+1].CP;
1967 bss[curd].DP = bss[curd+1].DP;
1969 c = bss[curd].CC + wa[B[curd+low-1+i]];
1970 if (c < d || c < e) {
1976 bss[curd].CP = bss[curd].DP;
1978 } /* otherwise, CP is unchanged */
1982 if ((c = c-m) > (e = e-h)) {
1984 MP[1][i] = bss[curd-1].CP;
1990 if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
1992 MP[2][i] = bss[curd+1].CP;
1995 MP[2][i] = bss[curd+1].DP;
1998 c = bss[curd].CC + wa[B[curd+low-1+i]];
1999 if (c < d || c < e) {
2002 MP[0][i] = MP[1][i];
2006 MP[0][i] = MP[2][i];
2014 MP[1][i] = MP[0][i];
2015 MT[1][i] = MT[0][i];
2018 MP[2][i] = MP[0][i];
2019 MT[2][i] = MT[0][i];
2021 bss[curd].CP = bss[curd].DP = IP = i;
2028 /* decide which path to be traced back */
2029 if (te == 1 && d+g > c) {
2032 } else if (te == 2 && e+g > c) {
2039 if (rmid > N-M) l = 2;
2040 else if (rmid < N-M) l = 1;
2043 /* Conquer: Solve subproblems recursively */
2047 for (; k > -1; r=k, k=MP[l][r], l=MT[l][r]){
2049 FT[k] = l; /* l=0,1,2 */
2051 /* forward dividing */
2052 if (r == -1) { /* optimal alignment did not cross the middle diagonal */
2054 bg_align(A,B,M,N,rmid+1,up,tb,te,w,g,h,bss, sapp, last);
2057 bg_align(A,B,M,N,low,rmid-1,tb,te,w,g,h,bss, sapp, last);
2066 bg_align(A,B,r-1,r+rmid,rmid+1,min(up,r+rmid),tb,1,w,g,h,bss,sapp,last);
2068 } else if (rmid > 0) {
2069 bg_align(A,B,r,r+rmid-1,max(-r,low),rmid-1,tb,2,w,g,h,bss,sapp,last);
2073 /* intermediate blocks */
2076 for (; l > -1; k = l, l = FP[k], kt = FT[k]) {
2077 if (kt == 0) { REP }
2078 else if (kt == 1) { /* right-hand side triangle */
2081 bg_align(A+k,B+k+rmid+1,t1,t1,0,min(t1,t2),2,1,w,g,h,bss,sapp,last);
2084 else { /* kt == 2, left-hand side triangle */
2087 bg_align(A+k+1,B+k+rmid,t1,t1,max(-t1,t3),0,1,2,w,g,h,bss,sapp,last);
2096 bg_align(A+k,B+t1,M-k,N-t1,0,min(N-t1,t2),2,te,w,g,h,bss,sapp,last);
2097 } else if (N-M < rmid) {
2100 bg_align(A+k+1,B+k+rmid,t1,N-(k+rmid),max(-t1,t3),0,1,te,w,g,h,
2107 int B_ALIGN(const unsigned char *A, const unsigned char *B,
2109 int low, int up, int **W, int G, int H, int *S, int *nS,
2110 int MW, int MX, struct bdstr *bss)
2116 int **sapp, *sapp_v, *last, last_v;
2126 low = min(max(-M, low),min(N-M,0));
2127 up = max(min(N, up),max(N-M,0));
2130 if (M > 0) { DEL(M); }
2137 if (up-low+1 <= 1) {
2139 for (i = 1; i <= M; i++) {
2148 MT[0] = (char *) ckalloc(mj);
2149 MT[1] = (char *) ckalloc(mj);
2150 MT[2] = (char *) ckalloc(mj);
2151 FT = (char *) ckalloc(mj);
2154 MP[0] = (int *) ckalloc(mj);
2155 MP[1] = (int *) ckalloc(mj);
2156 MP[2] = (int *) ckalloc(mj);
2157 FP = (int *) ckalloc(mj);
2160 c = bg_align(A,B,M,N,low,up,0,0,W,G,H,bss, sapp, last);
2162 check_score = BCHECK_SCORE(A,B,M,N,S,W,G,H,nS);
2164 free(FP); free(MP[2]); free(MP[1]); free(MP[0]);
2165 free(FT); free(MT[2]); free(MT[1]); free(MT[0]);
2168 if (check_score != c)
2169 printf("\nBCheck_score=%d != %d\n", check_score,c);
2173 int do_walign (const unsigned char *aa0, int n0,
2174 const unsigned char *aa1, int n1,
2176 struct pstruct *ppst,
2177 struct f_struct *f_str,
2178 struct a_res_str *a_res,
2181 int hoff, optflag_s, optcut_s, optwid_s, n10, score;
2182 const unsigned char *aa1p;
2186 f_str->n10 = n10=aatran(aa1,f_str->aa1x,n1,frame);
2187 do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, &rst, &hoff);
2195 a_res->res = f_str->res;
2199 return sw_walign(aa0, n0, aa1p, n10, ppst, f_str, a_res);
2201 optflag_s = ppst->param_u.fa.optflag;
2202 optcut_s = ppst->param_u.fa.optcut;
2203 optwid_s = ppst->param_u.fa.optwid;
2204 ppst->param_u.fa.optflag = 1;
2205 ppst->param_u.fa.optcut = 0;
2206 ppst->param_u.fa.optwid *= 2;
2208 do_fasta(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff);
2210 if (rst.score[0]>0) {
2211 score=bd_walign(aa0, n0, aa1p, n10, ppst, f_str, hoff, a_res);
2218 ppst->param_u.fa.optflag = optflag_s;
2219 ppst->param_u.fa.optcut = optcut_s;
2220 ppst->param_u.fa.optwid = optwid_s;
2226 pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {
2229 f_str->n10 = aatran(aa1,f_str->aa1x,n1,frame);
2233 /* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
2234 /* call from calcons, calc_id, calc_code */
2236 aln_func_vals(int frame, struct a_struct *aln) {
2246 aln->frame = 3 - frame;
2248 else aln->llrev = 0;
2250 aln->llfact = aln->qlfact = aln->llmult = 1;
2252 if (frame > 0) aln->qlrev = 1;
2253 else aln->qlrev = 0;
2260 int calcons(const unsigned char *aa0, int n0,
2261 const unsigned char *aa1, int n1,
2263 struct a_struct *aln, struct a_res_str a_res,
2265 char *seqc0, char *seqc1, char *seqca,
2266 struct f_struct *f_str)
2269 int op, lenc, nd, ns, itmp;
2270 const unsigned char *aa1p;
2271 char *sp0, *sp1, *spa, *sq;
2275 if (pst.ext_sq_set) { sq = pst.sqx; }
2276 else { sq = pst.sq; }
2286 aln->amin0 = a_res.min0;
2287 aln->amax0 = a_res.max0;
2288 aln->amin1 = a_res.min1;
2289 aln->amax1 = a_res.max1;
2290 /* will we show all the start ?*/
2291 if (min(a_res.min0,a_res.min1) < aln->llen || aln->showall==1)
2292 if (a_res.min0 >= a_res.min1) { /* aa0 extends more to left */
2294 if (aln->showall==1) mins = a_res.min0;
2295 else mins = min(a_res.min0,aln->llcntx);
2296 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2297 aln->smin0 = a_res.min0-mins;
2298 if ((mins-a_res.min1)>0) {
2299 memset(seqc1,' ',mins-a_res.min1);
2300 aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
2304 aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2305 aln->smin1 = a_res.min1-mins;
2310 if (aln->showall == 1) mins=a_res.min1;
2311 else mins = min(a_res.min1,aln->llcntx);
2312 aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
2313 aln->smin1 = a_res.min1-mins;
2314 if ((mins-a_res.min0)>0) {
2315 memset(seqc0,' ',mins-a_res.min0);
2316 aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
2320 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2321 aln->smin0 = a_res.min0-mins;
2325 mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
2327 aln->smin0=a_res.min0 - mins;
2328 aln->smin1=a_res.min1 - mins;
2329 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2330 aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2332 /* set the alignment code to zero for context */
2333 memset(seqca,0,mins);
2337 aln->smin0=a_res.min0;
2338 aln->smin1=a_res.min1;
2341 /* now get the middle */
2347 lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2351 while (i0 < a_res.max0 || i1 < a_res.max1) {
2352 if (op == 0 && *rp == 0) {
2355 if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
2356 else if (itmp == 0) { *spa = M_ZERO;}
2357 else {*spa = M_POS;}
2358 if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
2361 *sp0 = sq[aa0[i0++]];
2362 *sp1 = sq[aa1p[i1++]];
2363 if (toupper(*sp0) == toupper(*sp1)) {
2367 else if (pst.nt_align) {
2368 if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
2369 (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
2373 else if (toupper(*sp0) == 'N') aln->ngap_q++;
2374 else if (toupper(*sp1) == 'N') aln->ngap_l++;
2376 sp0++; sp1++; spa++;
2379 if (op==0) op = *rp++;
2382 *sp1++ = sq[aa1p[i1++]];
2389 *sp0++ = sq[aa0[i0++]];
2402 /* now we have the middle, get the right end */
2403 if (!aln->llcntx_flg) {
2404 ns = mins + lenc + aln->llen; /* show an extra line? */
2405 ns -= (itmp = ns %aln->llen); /* itmp = left over on last line */
2406 if (itmp>aln->llen/2) ns += aln->llen; /* more than 1/2 , use another*/
2407 nd = ns - (mins+lenc); /* this much extra */
2409 else nd = aln->llcntx;
2411 if (nd > max(n0-a_res.max0,nn1-a_res.max1))
2412 nd = max(n0-a_res.max0,nn1-a_res.max1);
2414 if (aln->showall==1) {
2415 nd = max(n0-a_res.max0,nn1-a_res.max1); /* reset for showall=1 */
2417 aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
2418 aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2419 /* fill with blanks - this is required to use one 'nc' */
2420 memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2421 memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2424 if ((nd-(n0-a_res.max0))>0) {
2425 aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
2426 memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2428 else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
2430 if ((nd-(nn1-a_res.max1))>0) {
2431 aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2432 memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2434 else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
2437 /* fprintf(stderr,"%d\n",mins+lenc+nd); */
2439 return mins+lenc+nd;
2442 int calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
2443 const unsigned char *aa1, int n1,
2445 struct a_struct *aln,
2446 struct a_res_str a_res,
2448 char *seqc0, char *seqc0a, char *seqc1, char *seqca,
2449 char *ann_arr, struct f_struct *f_str)
2452 int op, lenc, nd, ns, itmp;
2453 const unsigned char *aa1p;
2454 char *sp0, *sp0a, *sp1, *spa, *sq;
2458 if (pst.ext_sq_set) {
2473 aln->amin0 = a_res.min0;
2474 aln->amax0 = a_res.max0;
2475 aln->amin1 = a_res.min1;
2476 aln->amax1 = a_res.max1;
2477 /* will we show all the start ?*/
2478 /* will we show all the start ?*/
2479 if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1)
2480 if (a_res.min0>=a_res.min1) { /* aa0 extends more to left */
2482 if (aln->showall==1) mins = a_res.min0;
2483 else mins = min(a_res.min0,aln->llcntx);
2484 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2485 aln->smin0 = a_res.min0-mins;
2486 if ((mins-a_res.min1)>0) {
2487 memset(seqc1,' ',mins-a_res.min1);
2488 aancpy(seqc1+mins-a_res.min1,(char *)aa1p,a_res.min1,pst);
2492 aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2493 aln->smin1 = a_res.min1-mins;
2498 if (aln->showall == 1) mins=a_res.min1;
2499 else mins = min(a_res.min1,aln->llcntx);
2500 aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
2501 aln->smin1 = a_res.min1-mins;
2502 if ((mins-a_res.min0)>0) {
2503 memset(seqc0,' ',mins-a_res.min0);
2504 aancpy(seqc0+mins-a_res.min0,(char *)aa0,a_res.min0,pst);
2508 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2509 aln->smin0 = a_res.min0-mins;
2513 mins= min(aln->llcntx,min(a_res.min0,a_res.min1));
2515 aln->smin0=a_res.min0 - smins;
2516 aln->smin1=a_res.min1 - smins;
2517 aancpy(seqc0,(char *)aa0+a_res.min0-mins,mins,pst);
2518 aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
2520 /* set the alignment code to zero for context */
2521 memset(seqca,0,mins);
2522 memset(seqc0a,' ',mins);
2526 aln->smin0=a_res.min0;
2527 aln->smin1=a_res.min1;
2530 /* now get the middle */
2537 lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2541 while (i0 < a_res.max0 || i1 < a_res.max1) {
2542 if (op == 0 && *rp == 0) {
2545 if ((itmp=pst.pam2[0][aa0[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
2546 else if (itmp == 0) { *spa = M_ZERO;}
2547 else {*spa = M_POS;}
2548 if (*spa == M_POS || *spa==M_ZERO) aln->nsim++;
2550 *sp0a++ = ann_arr[aa0a[i0]];
2552 *sp0 = sq[aa0[i0++]];
2553 *sp1 = sq[aa1p[i1++]];
2555 if (toupper(*sp0) == toupper(*sp1)) {
2559 else if (pst.nt_align) {
2560 if ((toupper(*sp0) == 'T' && toupper(*sp1) == 'U') ||
2561 (toupper(*sp0)=='U' && toupper(*sp1)=='T')) {
2565 else if (toupper(*sp0) == 'N') aln->ngap_q++;
2566 else if (toupper(*sp1) == 'N') aln->ngap_l++;
2568 sp0++; sp1++; spa++;
2571 if (op==0) op = *rp++;
2574 *sp1++ = sq[aa1p[i1++]];
2582 *sp0a++ = ann_arr[aa0a[i0]];
2583 *sp0++ = sq[aa0[i0++]];
2594 *sp0a = *spa = '\0';
2596 /* now we have the middle, get the right end */
2597 if (!aln->llcntx_flg) {
2598 ns = mins + lenc + aln->llen; /* show an extra line? */
2599 ns -= (itmp = ns %aln->llen); /* itmp = left over on last line */
2600 if (itmp>aln->llen/2) ns += aln->llen; /* more than 1/2 , use another*/
2601 nd = ns - (mins+lenc); /* this much extra */
2603 else nd = aln->llcntx;
2605 if (nd > max(n0-a_res.max0,nn1-a_res.max1))
2606 nd = max(n0-a_res.max0,nn1-a_res.max1);
2608 if (aln->showall==1) {
2609 nd = max(n0-a_res.max0,nn1-a_res.max1); /* reset for showall=1 */
2611 aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,n0-a_res.max0,pst);
2612 aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2613 /* fill with blanks - this is required to use one 'nc' */
2614 memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2615 memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2618 if ((nd-(n0-a_res.max0))>0) {
2619 aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,(n0-a_res.max0),pst);
2620 memset(seqc0+mins+lenc+n0-a_res.max0,' ',nd-(n0-a_res.max0));
2622 else aancpy(seqc0+mins+lenc,(char *)aa0+a_res.max0,nd,pst);
2624 if ((nd-(nn1-a_res.max1))>0) {
2625 aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
2626 memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
2628 else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
2631 /* fprintf(stderr,"%d\n",mins+lenc+nd); */
2633 return mins+lenc+nd;
2637 update_code(char *al_str, int al_str_max, int op, int op_cnt) {
2639 char op_char[5]={"=-+*"};
2642 sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
2643 strncat(al_str,tmp_cnt,al_str_max);
2647 /* build an array of match/ins/del - length strings */
2648 int calc_code(const unsigned char *aa0, int n0,
2649 const unsigned char *aa1, int n1,
2650 struct a_struct *aln, struct a_res_str a_res,
2652 char *al_str, int al_str_n, struct f_struct *f_str)
2657 const unsigned char *aa1p;
2661 if (pst.ext_sq_set) {
2676 aln->amin0 = a_res.min0;
2677 aln->amax0 = a_res.max0;
2678 aln->amin1 = a_res.min1;
2679 aln->amax1 = a_res.max1;
2682 lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = p_op = 0;
2688 while (i0 < a_res.max0 || i1 < a_res.max1) {
2689 if (op == 0 && *rp == 0) {
2691 if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
2693 sp0 = sq[aa0[i0++]];
2694 sp1 = sq[aa1p[i1++]];
2696 if (p_op == 0 || p_op==3) {
2697 if (sp0 != '*' && sp1 != '*') {
2699 update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2700 op_cnt = 1; p_op = 0;
2705 update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2706 op_cnt = 1; p_op = 3;
2710 update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt);
2711 op_cnt = 1; p_op = 0;
2717 if (toupper(sp0) == toupper(sp1)) aln->nident++;
2718 else if (pst.nt_align) {
2719 if ((toupper(sp0) == 'T' && toupper(sp1) == 'U') ||
2720 (toupper(sp0)=='U' && toupper(sp1)=='T')) aln->nident++;
2721 else if (toupper(sp0) == 'N') aln->ngap_q++;
2722 else if (toupper(sp1) == 'N') aln->ngap_l++;
2726 if (op==0) op = *rp++;
2728 if (p_op == 1) { op_cnt++;}
2730 update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2731 op_cnt = 1; p_op = 1;
2733 op--; lenc++; i1++; aln->ngap_q++;
2736 if (p_op == 2) { op_cnt++;}
2738 update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2739 op_cnt = 1; p_op = 2;
2741 op++; lenc++; i0++; aln->ngap_l++;
2745 update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt);
2750 int calc_id(const unsigned char *aa0, int n0,
2751 const unsigned char *aa1, int n1,
2752 struct a_struct *aln,
2753 struct a_res_str a_res,
2755 struct f_struct *f_str)
2760 const unsigned char *aa1p;
2764 if (pst.ext_sq_set) { sq = pst.sqx; }
2765 else { sq = pst.sq; }
2775 aln->amin0 = a_res.min0;
2776 aln->amax0 = a_res.max0;
2777 aln->amin1 = a_res.min1;
2778 aln->amax1 = a_res.max1;
2781 lenc = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = 0;
2785 while (i0 < a_res.max0 || i1 < a_res.max1) {
2786 if (op == 0 && *rp == 0) {
2790 if (pst.pam2[0][aa0[i0]][aa1p[i1]]>=0) { aln->nsim++;}
2792 sp0 = sq[aa0[i0++]];
2793 sp1 = sq[aa1p[i1++]];
2794 if (toupper(sp0) == toupper(sp1)) {aln->nident++;}
2795 else if (pst.nt_align) {
2796 if ((toupper(sp0)=='T' && toupper(sp1)== 'U')||
2797 (toupper(sp0)=='U' && toupper(sp1)=='T')) {aln->nident++;}
2798 else if (toupper(sp0) == 'N') aln->ngap_q++;
2799 else if (toupper(sp1) == 'N') aln->ngap_l++;
2803 if (op==0) op = *rp++;
2804 if (op>0) {op--; lenc++; i1++; aln->ngap_q++;}
2805 else {op++; lenc++; i0++; aln->ngap_l++; }
2816 update_params(struct qmng_str *qm_msg, struct pstruct *ppst)
2818 ppst->n0 = qm_msg->n0;