binaries/src/fasta34/dropff2.c

   1
   2 /* copyright (c) 1998, 1999 William R. Pearson and the U. of Virginia */
   3
   4 /*  - dropffa.c,v 1.1.1.1 1999/10/22 20:55:59 wrp Exp */
   5
   6 /* this code implements the "fastf" algorithm, which is designed to
   7    deconvolve mixtures of protein sequences derived from mixed-peptide
   8    Edman sequencing.  The expected input is:
   9
  10    >test | 40001 90043 | mgstm1
  11    MGCEN,
  12    MIDYP,
  13    MLLAY,
  14    MLLGY
  15
  16    Where the ','s indicate the length/end of the sequencing cycle
  17    data.  Thus, in this example, the sequence is from a mixture of 4
  18    peptides, M was found in the first position, G,I, and L(2) at the second,
  19    C,D, L(2) at the third, etc.
  20
  21    Because the sequences are derived from mixtures, there need not be
  22    any partial sequence "MGCEN", the actual deconvolved sequence might be
  23    "MLDGN".
  24 */
  25
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29 #include <stdlib.h>
  30 #include <math.h>
  31 #include <ctype.h>
  32
  33 #include "defs.h"
  34 #include "param.h"
  35 #include "structs.h"
  36 #include "tatstats.h"
  37
  38 #define EOSEQ 0
  39 #define ESS 49
  40 #define MAXHASH 32
  41 #define NMAP MAXHASH+1
  42 #define NMAP_X 23       /* re-code NMAP for 'X' */
  43 #define NMAP_Z 24       /* re-code NMAP for '*' */
  44
  45 #ifndef MAXSAV
  46 #define MAXSAV 10
  47 #endif
  48
  49 #define DROP_INTERN
  50 #include "drop_func.h"
  51
  52 static char *verstr="4.21 May 2006 (ajm/wrp)";
  53
  54 int shscore(unsigned char *aa0, const int n0, int **pam2, int nsq);
  55 void update_code(char *al_str, int al_str_max, int op, int op_cnt, int fnum);
  56 extern void aancpy(char *to, char *from, int count, struct pstruct pst);
  57
  58 #ifdef TFAST
  59 extern int aatran(const unsigned char *ntseq, unsigned char *aaseq,
  60                   const int maxs, const int frame);
  61 #endif
  62
  63 struct hlstr { int next, pos;};
  64
  65 void savemax(struct dstruct *, struct f_struct *);
  66
  67 static int m0_spam(unsigned char *, const unsigned char *, int, struct savestr *,
  68             int **, struct f_struct *);
  69 static int m1_spam(unsigned char *, int,
  70                    const unsigned char *, int,
  71                    struct savestr *, int **, int, struct f_struct *);
  72
  73 int sconn(struct savestr **v, int nsave, int cgap,
  74           struct f_struct *, struct rstruct *, struct pstruct *,
  75           const unsigned char *aa0, int n0,
  76           const unsigned char *aa1, int n1,
  77           int opt_prob);
  78
  79 void kpsort(struct savestr **, int);
  80 void kssort(struct savestr **, int);
  81 void kpsort(struct savestr **, int);
  82
  83 int
  84 sconn_a(unsigned char *, int, int, struct f_struct *,
  85         struct a_res_str *);
  86
  87 /* initialize for fasta */
  88
  89 void
  90 init_work (unsigned char *aa0, int n0,
  91            struct pstruct *ppst,
  92            struct f_struct **f_arg)
  93 {
  94    int mhv, phv;
  95    int hmax;
  96    int i0, ii0, hv;
  97    struct f_struct *f_str;
  98
  99    int maxn0;
 100    int i, j, q;
 101    struct savestr *vmptr;
 102    int *res;
 103
 104    f_str = (struct f_struct *) calloc(1, sizeof(struct f_struct));
 105    if(f_str == NULL) {
 106      fprintf(stderr, "Couldn't calloc f_str\n");
 107      exit(1);
 108    }
 109
 110    ppst->sw_flag = 0;
 111
 112    /* fastf3 cannot work with lowercase symbols as low complexity;
 113       thus, NMAP must be disabled; this depends on aascii['X']  */
 114    if (ppst->hsq[NMAP_X] == NMAP ) {ppst->hsq[NMAP_X]=1;}
 115    if (ppst->hsq[NMAP_Z] == NMAP ) {ppst->hsq[NMAP_Z]=1;}
 116
 117    /*   this does not work for share ppst structs, as in threads */
 118    /*else {fprintf(stderr," cannot find 'X'==NMAP\n");} */
 119
 120    for (i0 = 1, mhv = -1; i0 <= ppst->nsq; i0++)
 121       if (ppst->hsq[i0] < NMAP && ppst->hsq[i0] > mhv) mhv = ppst->hsq[i0];
 122
 123    if (mhv <= 0) {
 124       fprintf (stderr, " maximum hsq <=0 %d\n", mhv);
 125       exit (1);
 126    }
 127
 128    for (f_str->kshft = 0; mhv > 0; mhv /= 2)
 129       f_str->kshft++;
 130
 131 /*      kshft = 2;      */
 132    hmax = hv = (1 << f_str->kshft);
 133    f_str->hmask = (hmax >> f_str->kshft) - 1;
 134
 135    if ((f_str->aa0 = (unsigned char *) calloc(n0+1, sizeof(char))) == NULL) {
 136      fprintf (stderr, " cannot allocate f_str->aa0 array; %d\n",n0+1);
 137      exit (1);
 138    }
 139    for (i=0; i<n0; i++) f_str->aa0[i] = aa0[i];
 140    aa0 = f_str->aa0;
 141
 142    if ((f_str->aa0t = (unsigned char *) calloc(n0+1, sizeof(char))) == NULL) {
 143      fprintf (stderr, " cannot allocate f_str0->aa0t array; %d\n",n0+1);
 144      exit (1);
 145    }
 146    f_str->aa0ix = 0;
 147
 148    if ((f_str->harr = (struct hlstr *) calloc (hmax, sizeof (struct hlstr))) == NULL) {
 149      fprintf (stderr, " cannot allocate hash array; hmax: %d hmask: %d\n",
 150               hmax,f_str->hmask);
 151      exit (1);
 152    }
 153    if ((f_str->pamh1 = (int *) calloc (ppst->nsq+1, sizeof (int))) == NULL) {
 154      fprintf (stderr, " cannot allocate pamh1 array\n");
 155      exit (1);
 156    }
 157    if ((f_str->pamh2 = (int *) calloc (hmax, sizeof (int))) == NULL) {
 158      fprintf (stderr, " cannot allocate pamh2 array\n");
 159      exit (1);
 160    }
 161    if ((f_str->link = (struct hlstr *) calloc (n0, sizeof (struct hlstr))) == NULL) {
 162      fprintf (stderr, " cannot allocate hash link array");
 163      exit (1);
 164    }
 165
 166    for (i0 = 0; i0 < hmax; i0++) {
 167       f_str->harr[i0].next = -1;
 168       f_str->harr[i0].pos = -1;
 169    }
 170
 171    for (i0 = 0; i0 < n0; i0++) {
 172       f_str->link[i0].next = -1;
 173       f_str->link[i0].pos = -1;
 174    }
 175
 176    /* encode the aa0 array */
 177    /*
 178      this code has been modified to allow for mixed peptide sequences
 179       aa0[] = 5 8 9 3 4 NULL 5 12 3 7 2 NULL
 180       the 'NULL' character resets the hash position counter, to indicate that
 181       any of several residues can be in the same position.
 182       We also need to keep track of the number of times this has happened, so that
 183       we can redivide the sequence later
 184
 185       i0 counts through the sequence
 186       ii0 counts through the hashed sequence
 187
 188       */
 189
 190    f_str->nm0 = 1;
 191    f_str->nmoff = -1;
 192    phv = hv = 0;
 193    for (i0= ii0 = 0; i0 < n0; i0++, ii0++) {
 194      /* reset the counter and start hashing again */
 195      if (aa0[i0] == ESS || aa0[i0] == 0) {
 196        aa0[i0] = 0;     /* set ESS to 0 */
 197        /*       fprintf(stderr," converted ',' to 0\n");*/
 198        i0++;    /* skip over the blank */
 199        f_str->nm0++;
 200        if (f_str->nmoff < 0) f_str->nmoff = i0;
 201        phv = hv = 0;
 202        ii0 = 0;
 203      }
 204      hv = ppst->hsq[aa0[i0]];
 205      f_str->link[i0].next = f_str->harr[hv].next;
 206      f_str->link[i0].pos = f_str->harr[hv].pos;
 207      f_str->harr[hv].next = i0;
 208      f_str->harr[hv].pos = ii0;
 209      f_str->pamh2[hv] = ppst->pam2[0][aa0[i0]][aa0[i0]];
 210    }
 211    if (f_str-> nmoff < 0) f_str->nmoff = n0;
 212
 213
 214 #ifdef DEBUG
 215    /*
 216    fprintf(stderr," nmoff: %d/%d nm0: %d\n", f_str->nmoff, n0,f_str->nm0);
 217    */
 218 #endif
 219
 220 /*
 221 #ifdef DEBUG
 222    fprintf(stderr," hmax: %d\n",hmax);
 223    for ( hv=0; hv<hmax; hv++)
 224        fprintf(stderr,"%2d %c %3d %3d\n",hv,
 225                (hv > 0 && hv < ppst->nsq ) ? ppst->sq[ppst->hsq[hv]] : ' ',
 226                f_str->harr[hv].pos,f_str->harr[hv].next);
 227    fprintf(stderr,"----\n");
 228    for ( hv=0; hv<n0; hv++)
 229        fprintf(stderr,"%2d: %3d %3d\n",hv,
 230                f_str->link[hv].pos,f_str->link[hv].next);
 231 #endif
 232 */
 233
 234    f_str->maxsav = MAXSAV;
 235    if ((f_str->vmax = (struct savestr *)
 236         calloc(MAXSAV,sizeof(struct savestr)))==NULL) {
 237      fprintf(stderr, "Couldn't allocate vmax[%d].\n",f_str->maxsav);
 238      exit(1);
 239    }
 240
 241    if ((f_str->vptr = (struct savestr **)
 242         calloc(MAXSAV,sizeof(struct savestr *)))==NULL) {
 243      fprintf(stderr, "Couldn't allocate vptr[%d].\n",f_str->maxsav);
 244      exit(1);
 245    }
 246
 247    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
 248      vmptr->used = (int *) calloc(n0, sizeof(int));
 249      if(vmptr->used == NULL) {
 250        fprintf(stderr, "Couldn't alloc vmptr->used\n");
 251        exit(1);
 252      }
 253    }
 254
 255 /* this has been modified from 0..<ppst->nsq to 1..<=ppst->nsq because the
 256    pam2[0][0] is now undefined for consistency with blast
 257 */
 258
 259    for (i0 = 1; i0 <= ppst->nsq; i0++)
 260      f_str->pamh1[i0] = ppst->pam2[0][i0][i0];
 261
 262    ppst->param_u.fa.cgap = shscore(aa0,f_str->nmoff-1,ppst->pam2[0],ppst->nsq)/3;
 263    if (ppst->param_u.fa.cgap > ppst->param_u.fa.bestmax/4)
 264      ppst->param_u.fa.cgap = ppst->param_u.fa.bestmax/4;
 265
 266    f_str->ndo = 0;
 267    f_str->noff = n0-1;
 268    if (f_str->diag==NULL)
 269      f_str->diag = (struct dstruct *) calloc ((size_t)MAXDIAG,
 270                                               sizeof (struct dstruct));
 271
 272    if (f_str->diag == NULL)
 273    {
 274       fprintf (stderr, " cannot allocate diagonal arrays: %ld\n",
 275               (long) MAXDIAG * (long) (sizeof (struct dstruct)));
 276       exit (1);
 277    }
 278
 279 #ifdef TFAST
 280    if ((f_str->aa1x =(unsigned char *)calloc((size_t)ppst->maxlen+2,
 281                                              sizeof(unsigned char)))
 282        == NULL) {
 283      fprintf (stderr, "cannot allocate aa1x array %d\n", ppst->maxlen+2);
 284      exit (1);
 285    }
 286    f_str->aa1x++;
 287 #endif
 288
 289    /* allocate space for the scoring arrays */
 290    maxn0 = n0 + 4;
 291
 292    maxn0 = max(3*n0/2,MIN_RES);
 293    if ((res = (int *)calloc((size_t)maxn0,sizeof(int)))==NULL) {
 294      fprintf(stderr,"cannot allocate alignment results array %d\n",maxn0);
 295      exit(1);
 296    }
 297    f_str->res = res;
 298    f_str->max_res = maxn0;
 299
 300    /* Tatusov Statistics Setup */
 301
 302    /* initialize priors array. */
 303    if((f_str->priors = (double *)calloc(ppst->nsq+1, sizeof(double))) == NULL) {
 304      fprintf(stderr, "Couldn't allocate priors array.\n");
 305      exit(1);
 306    }
 307    calc_priors(f_str->priors, ppst, f_str, NULL, 0, ppst->pseudocts);
 308
 309    f_str->dotat = 0;
 310    f_str->shuff_cnt = ppst->shuff_node;
 311
 312    /* End of Tatusov Statistics Setup */
 313
 314    *f_arg = f_str;
 315 }
 316
 317
 318 /* pstring1 is a message to the manager, currently 512 */
 319 /* pstring2 is the same information, but in a markx==10 format */
 320 void
 321 get_param (struct pstruct *pstr, char *pstring1, char *pstring2)
 322 {
 323 #ifndef TFAST
 324   char *pg_str="FASTF";
 325 #else
 326   char *pg_str="TFASTF";
 327 #endif
 328
 329   sprintf (pstring1, "%s (%s) function [%s matrix (%d:%d)] join: %d",pg_str,verstr,
 330            pstr->pamfile, pstr->pam_h,pstr->pam_l,pstr->param_u.fa.cgap);
 331
 332   if (pstr->param_u.fa.iniflag) strcat(pstring1," init1");
 333   /*
 334   if (pstr->zsflag==0) strcat(pstring1," not-scaled");
 335   else if (pstr->zsflag==1) strcat(pstring1," reg.-scaled");
 336   */
 337
 338   if (pstring2 != NULL) {
 339     sprintf (pstring2, "; pg_name: %s\n; pg_ver: %s\n; pg_matrix: %s (%d:%d)\n\
 340 ; pg_join: %d\n",
 341              pg_str,verstr, pstr->pamfile, pstr->pam_h,pstr->pam_l,
 342              pstr->param_u.fa.cgap);
 343    }
 344 }
 345
 346 void
 347 close_work (const unsigned char *aa0, const int n0,
 348             struct pstruct *ppst,
 349             struct f_struct **f_arg)
 350 {
 351   struct f_struct *f_str;
 352   struct savestr *vmptr;
 353
 354   f_str = *f_arg;
 355
 356   if (f_str != NULL) {
 357
 358     for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
 359       free(vmptr->used);
 360
 361     free(f_str->res);
 362 #ifdef TFAST
 363     free(f_str->aa1x - 1); /* allocated, then aa1x++'ed */
 364 #endif
 365     free(f_str->diag);
 366     free(f_str->link);
 367     free(f_str->pamh2);
 368     free(f_str->pamh1);
 369     free(f_str->harr);
 370     free(f_str->aa0t);
 371     free(f_str->aa0);
 372     free(f_str->priors);
 373     free(f_str);
 374     *f_arg = NULL;
 375   }
 376 }
 377
 378 int do_fastf (unsigned char *aa0, int n0,
 379               const unsigned char *aa1, int n1,
 380               struct pstruct *ppst, struct f_struct *f_str,
 381               struct rstruct *rst, int *hoff, int opt_prob)
 382 {
 383    int     nd;          /* diagonal array size */
 384    int     lhval;
 385    int     kfact;
 386    register struct dstruct *dptr;
 387    register int tscor;
 388    register struct dstruct *diagp;
 389    struct dstruct *dpmax;
 390    register int lpos;
 391    int     tpos, npos;
 392    struct savestr *vmptr;
 393    int     scor, tmp;
 394    int     im, ib, nsave;
 395    int     cmps ();             /* comparison routine for ksort */
 396    int *hsq;
 397
 398    hsq = ppst->hsq;
 399
 400    if (n1 < 1) {
 401      rst->score[0] = rst->score[1] = rst->score[2] = 0;
 402      rst->escore = 1.0;
 403      rst->segnum = 0;
 404      rst->seglen = 0;
 405      return 1;
 406    }
 407
 408    if (n0+n1+1 >= MAXDIAG) {
 409      fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
 410      rst->score[0] = rst->score[1] = rst->score[2] = -1;
 411      rst->escore = 2.0;
 412      rst->segnum = 0;
 413      rst->seglen = 0;
 414      return -1;
 415    }
 416
 417    nd = n0 + n1;
 418
 419    dpmax = &f_str->diag[nd];
 420    for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;) {
 421       dptr->stop = -1;
 422       dptr->dmax = NULL;
 423       dptr++->score = 0;
 424    }
 425
 426    /* initialize the saved segment structures */
 427    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
 428       vmptr->score = 0;
 429       memset(vmptr->used, 0, n0 * sizeof(int));
 430    }
 431
 432    f_str->lowmax = f_str->vmax;
 433    f_str->lowscor = 0;
 434
 435    /* start hashing */
 436
 437    diagp = &f_str->diag[f_str->noff];
 438    for (lhval = lpos = 0; lpos < n1; lpos++, diagp++) {
 439      if (hsq[aa1[lpos]]>=NMAP) {
 440        lpos++ ; diagp++;
 441        while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
 442        if (lpos >= n1) break;
 443        lhval = 0;
 444      }
 445      lhval = hsq[aa1[lpos]];
 446      for (tpos = f_str->harr[lhval].pos, npos = f_str->harr[lhval].next;
 447           tpos >= 0; tpos = f_str->link[npos].pos, npos = f_str->link[npos].next) {
 448        /* tscor gets position of end of current lpos diag run */
 449        if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
 450          tscor++;                /* move forward one */
 451          if ((tscor -= lpos) <= 0) { /* check for size of gap to this hit - */
 452                                      /* includes implicit -1 mismatch penalty */
 453            scor = dptr->score;       /* current score of this run */
 454            if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 &&
 455                f_str->lowscor < scor)   /* if updating tscor makes run worse, */
 456              savemax (dptr, f_str);     /* save it */
 457
 458            if ((tscor += scor) >= kfact) {  /* add to current run if continuing */
 459                                             /* is better than restart (kfact) */
 460                dptr->score = tscor;
 461                dptr->stop = lpos;
 462              }
 463              else {
 464                dptr->score = kfact;     /* starting over is better */
 465                dptr->start = (dptr->stop = lpos);
 466              }
 467          }
 468          else {                         /* continue current run */
 469            dptr->score += f_str->pamh1[aa0[tpos]];
 470            dptr->stop = lpos;
 471          }
 472        }
 473        else {                           /* no diagonal run yet */
 474          dptr->score = f_str->pamh2[lhval];
 475          dptr->start = (dptr->stop = lpos);
 476        }
 477      }                          /* end tpos */
 478    }                            /* end lpos */
 479
 480    for (dptr = f_str->diag; dptr < dpmax;) {
 481      if (dptr->score > f_str->lowscor) savemax (dptr, f_str);
 482      dptr->stop = -1;
 483      dptr->dmax = NULL;
 484      dptr++->score = 0;
 485    }
 486    f_str->ndo = nd;
 487
 488 /*
 489         at this point all of the elements of aa1[lpos]
 490         have been searched for elements of aa0[tpos]
 491         with the results in diag[dpos]
 492 */
 493
 494    /* set up pointers for sorting */
 495
 496    for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
 497      if (vmptr->score > 0) {
 498        vmptr->score = m0_spam (aa0, aa1, n1, vmptr, ppst->pam2[0], f_str);
 499        f_str->vptr[nsave++] = vmptr;
 500      }
 501    }
 502
 503    /* sort them */
 504    kssort (f_str->vptr, nsave);
 505
 506
 507 #ifdef DEBUG
 508    /*
 509    for (ib=0; ib<nsave; ib++) {
 510      fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
 511              f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
 512              f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
 513              f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
 514              f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
 515        for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
 516          fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
 517                  ppst->sq[aa1[im]]);
 518        fputc('\n',stderr);
 519    }
 520    fprintf(stderr,"---\n");
 521    */
 522    /* now use m_spam to re-evaluate */
 523    /*
 524    for (tpos = 0; tpos < n0; tpos++) {
 525      fprintf(stderr,"%c:%2d ",ppst->sq[aa0[tpos]],aa0[tpos]);
 526      if (tpos %10 == 9) fputc('\n',stderr);
 527    }
 528    fputc('\n',stderr);
 529    */
 530 #endif
 531
 532    f_str->aa0ix = 0;
 533    for (ib=0; ib < nsave; ib++) {
 534      if ((vmptr=f_str->vptr[ib])->score > 0) {
 535        vmptr->score = m1_spam (aa0, n0, aa1, n1, vmptr,
 536                                ppst->pam2[0], ppst->pam_l, f_str);
 537      }
 538    }
 539    /* reset aa0 - modified by m1_spam */
 540    for (tpos = 0; tpos < n0; tpos++) {
 541      if (aa0[tpos] >= 32) aa0[tpos] -= 32;
 542    }
 543
 544    kssort(f_str->vptr,nsave);
 545
 546    for ( ; nsave > 0; nsave--)
 547      if (f_str->vptr[nsave-1]->score >0) break;
 548
 549    if (nsave <= 0) {
 550      f_str->nsave = 0;
 551      rst->score[0] = rst->score[1] = rst->score[2] = 0;
 552      rst->escore = 1.0;
 553
 554      return 1;
 555    }
 556    else f_str->nsave = nsave;
 557
 558
 559 #ifdef DEBUG
 560    /*
 561    fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,f_str->noff);
 562    for (ib=0; ib<nsave; ib++) {
 563      fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
 564              f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
 565              f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
 566              f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
 567              f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
 568      for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
 569        fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
 570                ppst->sq[aa1[im]]);
 571      fputc('\n',stderr);
 572    }
 573
 574    fprintf(stderr,"---\n");
 575    */
 576 #endif
 577
 578    scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap, f_str,
 579                  rst, ppst, aa0, n0, aa1, n1, opt_prob);
 580
 581    for (vmptr=f_str->vptr[0],ib=1; ib<nsave; ib++)
 582      if (f_str->vptr[ib]->score > vmptr->score) vmptr=f_str->vptr[ib];
 583
 584    rst->score[1] = vmptr->score;
 585    rst->score[0] = rst->score[2] = max (scor, vmptr->score);
 586
 587    return 1;
 588 }
 589
 590 void do_work (const unsigned char *aa0, int n0,
 591               const unsigned char *aa1, int n1,
 592               int frame,
 593               struct pstruct *ppst, struct f_struct *f_str,
 594               int qr_flg, struct rstruct *rst)
 595 {
 596   int opt_prob;
 597   int hoff, n10, i;
 598
 599   if (qr_flg==1 && f_str->shuff_cnt <= 0) {
 600     rst->escore = 2.0;
 601     rst->score[0]=rst->score[1]=rst->score[2]= -1;
 602     return;
 603   }
 604
 605   if (f_str->dotat || ppst->zsflag == 4 || ppst->zsflag == 14 ) opt_prob=1;
 606   else opt_prob = 0;
 607   if (ppst->zsflag == 2 || ppst->zsflag == 12) opt_prob = 0;
 608   if (qr_flg) {
 609     opt_prob=1;
 610     /*    if (frame==1) */
 611       f_str->shuff_cnt--;
 612   }
 613
 614   if (n1 < 1) {
 615     rst->score[0] = rst->score[1] = rst->score[2] = -1;
 616     rst->escore = 2.0;
 617     return;
 618   }
 619
 620 #ifdef TFAST
 621   n10=aatran(aa1,f_str->aa1x,n1,frame);
 622   if (ppst->debug_lib)
 623     for (i=0; i<n10; i++)
 624       if (f_str->aa1x[i]>ppst->nsq) {
 625         fprintf(stderr,
 626                 "residue[%d/%d] %d range (%d)\n",i,n1,
 627                 f_str->aa1x[i],ppst->nsq);
 628         f_str->aa1x[i]=0;
 629         n10=i-1;
 630       }
 631
 632   do_fastf (f_str->aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, opt_prob);
 633 #else   /* FASTF */
 634   do_fastf (f_str->aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, opt_prob);
 635 #endif
 636
 637   rst->comp = rst->H = -1.0;
 638
 639 }
 640
 641 void do_opt (const unsigned char *aa0, int n0,
 642              const unsigned char *aa1, int n1,
 643              int frame,
 644              struct pstruct *ppst,
 645              struct f_struct *f_str,
 646              struct rstruct *rst)
 647 {
 648   int optflag, tscore, hoff, n10;
 649
 650   optflag = ppst->param_u.fa.optflag;
 651   ppst->param_u.fa.optflag = 1;
 652
 653 #ifdef TFAST
 654   n10=aatran(aa1,f_str->aa1x,n1,frame);
 655   do_fastf (f_str->aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, 1);
 656 #else   /* FASTA */
 657   do_fastf(f_str->aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, 1);
 658 #endif
 659   ppst->param_u.fa.optflag = optflag;
 660 }
 661
 662 void
 663 savemax (dptr, f_str)
 664   register struct dstruct *dptr;
 665   struct f_struct *f_str;
 666 {
 667    register int dpos;
 668    register struct savestr *vmptr;
 669    register int i;
 670
 671    dpos = (int) (dptr - f_str->diag);
 672
 673 /* check to see if this is the continuation of a run that is already saved */
 674
 675    if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
 676          vmptr->start == dptr->start)
 677    {
 678       vmptr->stop = dptr->stop;
 679       if ((i = dptr->score) <= vmptr->score)
 680          return;
 681       vmptr->score = i;
 682       if (vmptr != f_str->lowmax)
 683          return;
 684    }
 685    else
 686    {
 687       i = f_str->lowmax->score = dptr->score;
 688       f_str->lowmax->dp = dpos;
 689       f_str->lowmax->start = dptr->start;
 690       f_str->lowmax->stop = dptr->stop;
 691       dptr->dmax = f_str->lowmax;
 692    }
 693
 694    for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
 695       if (vmptr->score < i)
 696       {
 697          i = vmptr->score;
 698          f_str->lowmax = vmptr;
 699       }
 700    f_str->lowscor = i;
 701 }
 702
 703 /* this version of spam() is designed to work with a collection of
 704    subfragments, selecting the best amino acid at each position so
 705    that, from each subfragment, each position is only used once.
 706
 707    As a result, m_spam needs to know the number of fragments.
 708
 709    In addition, it now requires a global alignment to the fragment
 710    and resets the start and stop positions
 711
 712    */
 713
 714 static int
 715 m1_spam (unsigned char *aa0, int n0,
 716          const unsigned char *aa1, int n1,
 717          struct savestr *dmax, int **pam2, int pam_l,
 718          struct f_struct *f_str)
 719 {
 720   int     tpos, lpos, im, ii, nm, ci;
 721   int     tot, ctot, pv;
 722
 723   struct {
 724     int     start, stop, score;
 725   } curv, maxv;
 726   unsigned char *aa0p;
 727   const unsigned char *aa1p;
 728
 729   lpos = dmax->start;                   /* position in library sequence */
 730    tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
 731    /* force global alignment, reset start*/
 732    if (tpos < lpos) {
 733      lpos = dmax->start -= tpos;
 734      tpos = 0;
 735    }
 736    else {
 737      tpos -= lpos;
 738      lpos = dmax->start = 0;
 739    }
 740
 741    dmax->stop = dmax->start + (f_str->nmoff -2 - tpos);
 742    if (dmax->stop > n1) dmax->stop = n1;
 743
 744    /*
 745    if (dmax->start < 0) {
 746      tpos = -dmax->start;
 747      lpos = dmax->start=0;
 748    }
 749    else tpos = 0;
 750    */
 751
 752    aa1p = &aa1[lpos];
 753    aa0p = &aa0[tpos];
 754
 755    nm = f_str->nm0;
 756
 757    tot = curv.score = maxv.score = 0;
 758    for (; lpos <= dmax->stop; lpos++,aa0p++,aa1p++) {
 759      ctot = pam_l;
 760      ci = -1;
 761      for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
 762        if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
 763          ctot = pv;
 764          ci = ii;
 765 /*       fprintf(stderr, "lpos: %d im: %d ii: %d ci: %d ctot: %d pi: %d pv: %d\n", lpos, im, ii, ci, ctot, aa0p[ii], pam2[aa0p[ii]][*aa1p]); */
 766        }
 767      }
 768      tot += ctot;
 769      if (ci >= 0 && aa0p[ci] < 32) {
 770 #ifdef DEBUG
 771 /*         fprintf(stderr, "used: lpos: %d ci: %d : %c\n", lpos, ci, sq[aa0p[ci]]); */
 772 #endif
 773        aa0p[ci] +=  32;
 774        dmax->used[&aa0p[ci] - aa0] = 1;
 775      }
 776    }
 777    return tot;
 778 }
 779
 780 int ma_spam (unsigned char *aa0, int n0, const unsigned char *aa1,
 781              struct savestr *dmax, struct pstruct *ppst,
 782              struct f_struct *f_str)
 783 {
 784   int **pam2;
 785   int     tpos, lpos, im, ii, nm, ci, lp0;
 786   int     tot, ctot, pv;
 787   struct {
 788     int     start, stop, score;
 789   } curv, maxv;
 790    const unsigned char *aa1p;
 791    unsigned char *aa0p, *aa0pt;
 792    int aa0t_flg;
 793
 794    pam2 = ppst->pam2[0];
 795    aa0t_flg = 0;
 796
 797    lpos = dmax->start;                  /* position in library sequence */
 798    tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
 799    lp0 = lpos = dmax->start;
 800    aa1p = &aa1[lpos];
 801    aa0p = &aa0[tpos];                   /* real aa0 sequence */
 802
 803                         /* the destination aa0 sequence (without nulls) */
 804    aa0pt = &f_str->aa0t[f_str->aa0ix];
 805
 806    curv.start = lpos;
 807    nm = f_str->nm0;
 808
 809    /* sometimes, tpos may be > 0, with lpos = 0 - fill with 'X' */
 810    if (lpos == 0 && tpos > 0)
 811      for (ii = 0; ii < tpos; ii++) *aa0pt++ = 31;  /* filler character */
 812
 813    tot = curv.score = maxv.score = 0;
 814    for (; lpos <= dmax->stop; lpos++) {
 815      ctot = ppst->pam_l;
 816      ci = -1;
 817      for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
 818        if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
 819          ctot = pv;
 820          ci = ii;
 821        }
 822      }
 823      tot += ctot;
 824      if (ci >= 0) {
 825        if (ci >= n0) {fprintf(stderr," warning - ci off end %d/%d\n",ci,n0);}
 826        else {
 827          *aa0pt++ = aa0p[ci];
 828          aa0p[ci] +=  32;
 829          aa0t_flg=1;
 830        }
 831      }
 832      aa0p++; aa1p++;
 833    }
 834
 835    if (aa0t_flg) {
 836      dmax->dp -= f_str->aa0ix;          /* shift ->dp for aa0t */
 837      if ((ci=(int)(aa0pt-f_str->aa0t)) > n0) {
 838        fprintf(stderr," warning - aapt off %d/%d end\n",ci,n0);
 839      }
 840      else
 841        *aa0pt++ = 0;                    /* skip over NULL */
 842
 843      aa0pt = &f_str->aa0t[f_str->aa0ix];
 844      aa1p = &aa1[lp0];
 845
 846      /*
 847      for (im = 0; im < f_str->nmoff; im++)
 848        fprintf(stderr,"%c:%c,",ppst->sq[aa0pt[im]],ppst->sq[aa1p[im]]);
 849      fprintf(stderr,"- %3d (%3d:%3d)\n",dmax->score,f_str->aa0ix,lp0);
 850      */
 851
 852      f_str->aa0ix += f_str->nmoff;      /* update offset into aa0t */
 853    }
 854    /*
 855       fprintf(stderr," ma_spam returning: %d\n",tot);
 856    */
 857    return tot;
 858 }
 859
 860 static int
 861 m0_spam (unsigned char *aa0, const unsigned char *aa1, int n1,
 862          struct savestr *dmax, int **pam2,
 863          struct f_struct *f_str)
 864 {
 865    int tpos, lpos, lend, im, ii, nm;
 866    int     tot, ctot, pv;
 867    struct {
 868      int     start, stop, score;
 869    } curv, maxv;
 870    const unsigned char *aa0p, *aa1p;
 871
 872    lpos = dmax->start;                  /* position in library sequence */
 873    tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
 874    if (tpos > 0) {
 875      if (lpos-tpos >= 0) {
 876        lpos = dmax->start -= tpos;    /* force global alignment, reset start*/
 877        tpos = 0;
 878      }
 879      else {
 880        tpos -= lpos;
 881        lpos = dmax->start = 0;
 882      }
 883    }
 884
 885    nm = f_str->nm0;
 886    lend = dmax->stop;
 887    if (n1 - (lpos + f_str->nmoff-2) < 0 ) {
 888      lend = dmax->stop = (lpos - tpos) + f_str->nmoff-2;
 889      if (lend >= n1) lend = n1-1;
 890    }
 891
 892    aa1p = &aa1[lpos];
 893    aa0p = &aa0[tpos];
 894
 895    curv.start = lpos;
 896
 897    tot = curv.score = maxv.score = 0;
 898    for (; lpos <= lend; lpos++) {
 899      ctot = -10000;
 900      for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
 901        if ((pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
 902          ctot = pv;
 903        }
 904      }
 905      tot += ctot;
 906      aa0p++; aa1p++;
 907    }
 908
 909    /* reset dmax if necessary */
 910
 911    return tot;
 912 }
 913
 914 /* sconn links up non-overlapping alignments and calculates the score */
 915
 916 int sconn (struct savestr **v, int n, int cgap, struct f_struct *f_str,
 917            struct rstruct *rst, struct pstruct *ppst,
 918            const unsigned char *aa0, int n0,
 919            const unsigned char *aa1, int n1,
 920            int opt_prob)
 921 {
 922    int     i, si, cmpp ();
 923    struct slink *start, *sl, *sj, *so, sarr[MAXSAV];
 924    int     lstart, plstop;
 925    double tatprob;
 926
 927    /* sarr[] saves each alignment score/position, and provides a link
 928       back to the previous alignment that maximizes the score */
 929
 930    /*   sort the score left to right in lib pos */
 931    kpsort (v, n);
 932
 933    start = NULL;
 934
 935    /* for the remaining runs, see if they fit */
 936    for (i = 0, si = 0; i < n; i++) {
 937
 938      /* if the score is less than the gap penalty, it never helps */
 939      if (!opt_prob && (v[i]->score < cgap) ){ continue; }
 940
 941      lstart = v[i]->start;
 942
 943      /* put the run in the group */
 944      sarr[si].vp = v[i];
 945      sarr[si].score = v[i]->score;
 946      sarr[si].next = NULL;
 947      sarr[si].prev = NULL;
 948      sarr[si].tat = NULL;
 949
 950      if(opt_prob) {
 951        sarr[si].tatprob =
 952          calc_tatusov(NULL, &sarr[si], aa0, n0, aa1, n1,
 953                       ppst->pam2[0],ppst->nsq, f_str,
 954                       ppst->pseudocts, opt_prob,ppst->zsflag);
 955        sarr[si].tat = sarr[si].newtat;
 956      }
 957
 958      /* if it fits, then increase the score */
 959      for (sl = start; sl != NULL; sl = sl->next) {
 960        plstop = sl->vp->stop;
 961        /* if end < start or start > end, add score */
 962        if (plstop < lstart ) {
 963          if(!opt_prob) {
 964            sarr[si].score = sl->score + v[i]->score;
 965            sarr[si].prev = sl;
 966            /*
 967              fprintf(stderr,"sconn %d added %d/%d getting %d; si: %d, tat: %g\n",
 968              i,v[i]->start, v[i]->score,sarr[si].score,si, 2.0);
 969            */
 970            break;
 971          } else {
 972            tatprob =
 973              calc_tatusov(sl, &sarr[si], aa0, n0, aa1, n1,
 974                           ppst->pam2[0], ppst->nsq, f_str,
 975                           ppst->pseudocts, opt_prob, ppst->zsflag);
 976            /* if our tatprob gets worse when we add this, forget it */
 977            if(tatprob > sarr[si].tatprob) {
 978              free(sarr[si].newtat->probs); /* get rid of new tat struct */
 979              free(sarr[si].newtat);
 980              continue;
 981            } else {
 982              sarr[si].tatprob = tatprob;
 983              free(sarr[si].tat->probs); /* get rid of old tat struct */
 984              free(sarr[si].tat);
 985              sarr[si].tat = sarr[si].newtat;
 986              sarr[si].prev = sl;
 987              sarr[si].score = sl->score + v[i]->score;
 988              /*
 989                fprintf(stderr,"sconn TAT %d added %d/%d getting %d; si: %d, tat: %g\n",
 990                i,v[i]->start, v[i]->score,sarr[si].score,si, tatprob);
 991              */
 992              break;
 993            }
 994          }
 995        }
 996      }
 997
 998      /* now recalculate where the score fits - resort the scores */
 999      if (start == NULL) {
1000        start = &sarr[si];
1001      } else {
1002        if(!opt_prob) { /* sort by scores */
1003          for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
1004            if (sarr[si].score > sj->score) { /* if new score > best score */
1005              sarr[si].next = sj;             /* previous best linked to best */
1006              if (so != NULL)
1007                so->next = &sarr[si];         /* old best points to new best */
1008              else
1009                start = &sarr[si];
1010              break;
1011            }
1012            so = sj;                          /* old-best saved in so */
1013          }
1014        } else { /* sort by tatprobs */
1015          for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
1016            if ( sarr[si].tatprob < sj->tatprob ||
1017                 ((sarr[si].tatprob == sj->tatprob) && sarr[si].score > sj->score) ) {
1018              sarr[si].next = sj;
1019              if (so != NULL)
1020                so->next = &sarr[si];
1021              else
1022                start = &sarr[si];
1023              break;
1024            }
1025            so = sj;
1026          }
1027        }
1028      }
1029      si++;
1030    }
1031
1032    if(opt_prob) {
1033      for (i = 0 ; i < si ; i++) {
1034        free(sarr[i].tat->probs);
1035        free(sarr[i].tat);
1036      }
1037    }
1038
1039    if (start != NULL) {
1040
1041      if(opt_prob)
1042        rst->escore = start->tatprob;
1043      else
1044        rst->escore = 2.0;
1045
1046      rst->segnum = rst->seglen = 0;
1047      for(sj = start ; sj != NULL; sj = sj->prev) {
1048        rst->segnum++;
1049        rst->seglen += sj->vp->stop - sj->vp->start + 1;
1050      }
1051      return (start->score);
1052    } else {
1053
1054      if(opt_prob)
1055        rst->escore = 1.0;
1056      else
1057        rst->escore = 2.0;
1058
1059      rst->segnum = rst->seglen = 0;
1060      return (0);
1061    }
1062 }
1063
1064 void
1065 kssort (struct savestr **v, int n)
1066 {
1067    int     gap, i, j;
1068    struct savestr *tmp;
1069
1070    for (gap = n / 2; gap > 0; gap /= 2)
1071       for (i = gap; i < n; i++)
1072          for (j = i - gap; j >= 0; j -= gap)
1073          {
1074             if (v[j]->score >= v[j + gap]->score)
1075                break;
1076             tmp = v[j];
1077             v[j] = v[j + gap];
1078             v[j + gap] = tmp;
1079          }
1080 }
1081 void
1082 kpsort (v, n)
1083 struct savestr *v[];
1084 int     n;
1085 {
1086    int     gap, i, j;
1087    struct savestr *tmp;
1088
1089    for (gap = n / 2; gap > 0; gap /= 2)
1090       for (i = gap; i < n; i++)
1091          for (j = i - gap; j >= 0; j -= gap)
1092          {
1093             if (v[j]->start <= v[j + gap]->start)
1094                break;
1095             tmp = v[j];
1096             v[j] = v[j + gap];
1097             v[j + gap] = tmp;
1098          }
1099 }
1100
1101 /* sorts alignments from right to left (back to front) based on stop */
1102
1103 void
1104 krsort (v, n)
1105 struct savestr *v[];
1106 int     n;
1107 {
1108    int     gap, i, j;
1109    struct savestr *tmp;
1110
1111    for (gap = n / 2; gap > 0; gap /= 2)
1112       for (i = gap; i < n; i++)
1113          for (j = i - gap; j >= 0; j -= gap)
1114          {
1115             if (v[j]->stop > v[j + gap]->stop)
1116                break;
1117             tmp = v[j];
1118             v[j] = v[j + gap];
1119             v[j + gap] = tmp;
1120          }
1121 }
1122
1123 int  do_walign (const unsigned char *aa0, int n0,
1124                 const unsigned char *aa1, int n1,
1125                 int frame,
1126                 struct pstruct *ppst,
1127                 struct f_struct *f_str,
1128                 struct a_res_str *a_res,
1129                 int *have_ares)
1130 {
1131   int hoff, n10;
1132   struct rstruct rst;
1133   int ib;
1134   unsigned char *aa0t;
1135   const unsigned char *aa1p;
1136
1137 #ifdef TFAST
1138   f_str->n10 = n10 = aatran(aa1,f_str->aa1x,n1,frame);
1139   aa1p = f_str->aa1x;
1140 #else
1141   n10 = n1;
1142   aa1p = aa1;
1143 #endif
1144
1145   do_fastf(f_str->aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff, 1);
1146
1147   /* the alignment portion takes advantage of the information left
1148      over in f_str after do_fastf is done.  in particular, it is
1149      easy to run a modified sconn() to produce the alignments.
1150
1151      unfortunately, the alignment display routine wants to have
1152      things encoded as with bd_align and sw_align, so we need to do that.
1153      */
1154
1155   if ((aa0t = (unsigned char *)calloc(n0+1,sizeof(unsigned char)))==NULL) {
1156     fprintf(stderr," cannot allocate aa0t %d\n",n0+1);
1157     exit(1);
1158   }
1159
1160    kssort (f_str->vptr, f_str->nsave);
1161    f_str->aa0ix = 0;
1162    if (f_str->nsave > f_str->nm0) f_str->nsave = f_str->nm0;
1163    for (ib=0; ib < f_str->nm0; ib++) {
1164      if (f_str->vptr[ib]->score > 0) {
1165        f_str->vptr[ib]->score =
1166          ma_spam (f_str->aa0, n0, aa1p, f_str->vptr[ib], ppst, f_str);
1167      }
1168    }
1169
1170    /* after ma_spam is over, we need to reset aa0 */
1171    for (ib = 0; ib < n0; ib++) {
1172      if (f_str->aa0[ib] >= 32) f_str->aa0[ib] -= 32;
1173    }
1174
1175    kssort(f_str->vptr,f_str->nsave);
1176
1177    for ( ; f_str->nsave > 0; f_str->nsave--)
1178      if (f_str->vptr[f_str->nsave-1]->score >0) break;
1179
1180   a_res->nres = sconn_a (aa0t,n0, ppst->param_u.fa.cgap, f_str,a_res);
1181   free(aa0t);
1182
1183   a_res->res = f_str->res;
1184   *have_ares = 0;
1185   return rst.score[0];
1186 }
1187
1188 /* this version of sconn is modified to provide alignment information */
1189
1190 int sconn_a (unsigned char *aa0, int n0, int cgap,
1191              struct f_struct *f_str,
1192              struct a_res_str *a_res)
1193 {
1194    int     i, si, cmpp (), n;
1195    unsigned char *aa0p;
1196    int sx, dx, doff;
1197
1198    struct savestr **v;
1199    struct slink {
1200      int     score;
1201      struct savestr *vp;
1202      struct slink *snext;
1203      struct slink *aprev;
1204    } *start, *sl, *sj, *so, sarr[MAXSAV];
1205    int     lstop, plstart;
1206    int *res, nres, tres;
1207
1208 /*      sort the score left to right in lib pos */
1209
1210    v = f_str->vptr;
1211    n = f_str->nsave;
1212
1213    krsort (v, n);       /* sort from left to right in library */
1214
1215    start = NULL;
1216
1217 /*      for each alignment, see if it fits */
1218
1219    for (i = 0, si = 0; i < n; i++) {
1220
1221 /*      if the score is less than the join threshold, skip it */
1222      if (v[i]->score < cgap) continue;
1223
1224      lstop = v[i]->stop;                /* have right-most lstart */
1225
1226 /*      put the alignment in the group */
1227
1228      sarr[si].vp = v[i];
1229      sarr[si].score = v[i]->score;
1230      sarr[si].snext = NULL;
1231      sarr[si].aprev = NULL;
1232
1233 /*      if it fits, then increase the score */
1234 /* start points to a sorted (by total score) list of candidate
1235    overlaps */
1236
1237      for (sl = start; sl != NULL; sl = sl->snext) {
1238        plstart = sl->vp->start;
1239        if (plstart > lstop ) {
1240          sarr[si].score = sl->score + v[i]->score;
1241          sarr[si].aprev = sl;
1242          break;         /* quit as soon as the alignment has been added */
1243        }
1244      }
1245
1246 /* now recalculate the list of best scores */
1247      if (start == NULL)
1248        start = &sarr[si];       /* put the first one in the list */
1249      else
1250        for (sj = start, so = NULL; sj != NULL; sj = sj->snext) {
1251          if (sarr[si].score > sj->score) { /* new score better than old */
1252            sarr[si].snext = sj;         /* snext best after new score */
1253            if (so != NULL)
1254              so->snext = &sarr[si];     /* prev_best->snext points to best */
1255            else  start = &sarr[si];     /* start points to best */
1256            break;                       /* stop looking */
1257          }
1258          so = sj;               /* previous candidate best */
1259        }
1260      si++;                              /* increment to snext alignment */
1261    }
1262
1263    /* we have the best set of alignments, write them to *res */
1264    if (start != NULL) {
1265      res = f_str->res;  /* set a destination for the alignment ops */
1266      tres = nres = 0;   /* alignment op length = 0 */
1267      aa0p = aa0;        /* point into query (needed for calcons later) */
1268      a_res->min1 = start->vp->start;    /* start in library */
1269      a_res->min0 = 0;                   /* start in query */
1270      for (sj = start; sj != NULL; sj = sj->aprev ) {
1271        doff = (int)(aa0p-aa0) - (sj->vp->start-sj->vp->dp+f_str->noff);
1272        /*
1273        fprintf(stderr,"doff: %3d\n",doff);
1274        */
1275        for (dx=sj->vp->start,sx=sj->vp->start-sj->vp->dp+f_str->noff;
1276             dx <= sj->vp->stop; dx++) {
1277          *aa0p++ = f_str->aa0t[sx++];   /* copy residue into aa0 */
1278          tres++;                        /* bump alignment counter */
1279          res[nres++] = 0;               /* put 0-op in res */
1280        }
1281        sj->vp->dp -= doff;
1282        if (sj->aprev != NULL) {
1283          if (sj->aprev->vp->start - sj->vp->stop - 1 > 0 )
1284          /* put an insert op into res to get to next aligned block */
1285            tres += res[nres++] = (sj->aprev->vp->start - sj->vp->stop - 1);
1286        }
1287        /*
1288        fprintf(stderr,"t0: %3d, tx: %3d, l0: %3d, lx: %3d, dp: %3d noff: %3d, score: %3d\n",
1289                sj->vp->start - sj->vp->dp + f_str->noff,
1290                sj->vp->stop - sj->vp->dp + f_str->noff,
1291                sj->vp->start,sj->vp->stop,sj->vp->dp,
1292                f_str->noff,sj->vp->score);
1293        fprintf(stderr,"%3d - %3d: %3d\n",
1294                sj->vp->start,sj->vp->stop,sj->vp->score);
1295        */
1296        a_res->max1 = sj->vp->stop;
1297        a_res->max0 = a_res->max1 - sj->vp->dp + f_str->noff;
1298      }
1299
1300      /*
1301      fprintf(stderr,"(%3d - %3d):(%3d - %3d)\n",
1302      a_res->min0,a_res->max0,a_res->min1,a_res->max1);
1303      */
1304
1305      /* now replace f_str->aa0t with aa0 */
1306      for (i=0; i<n0; i++) f_str->aa0t[i] = aa0[i];
1307
1308      return tres;
1309    }
1310    else return (0);
1311 }
1312
1313 /* calculate the 100% identical score */
1314 int
1315 shscore(unsigned char *aa0, int n0, int **pam2, int nsq)
1316 {
1317   int i, sum;
1318   for (i=0,sum=0; i<n0; i++)
1319     if (aa0[i]!=0 && aa0[i]<=nsq) sum += pam2[aa0[i]][aa0[i]];
1320   return sum;
1321 }
1322
1323 void
1324 pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {
1325
1326 #ifdef TFAST
1327   f_str->n10=aatran(aa1,f_str->aa1x,n1,frame);
1328 #endif
1329 }
1330
1331 /* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
1332 /* call from calcons, calc_id, calc_code */
1333 void
1334 aln_func_vals(int frame, struct a_struct *aln) {
1335
1336 #ifdef TFAST
1337   aln->qlrev = 0;
1338   aln->qlfact = 1;
1339   aln->llfact = aln->llmult = 3;
1340   aln->frame = 0;
1341   if (frame > 3) aln->llrev = 1;
1342 #else   /* FASTF */
1343   aln->llfact = aln->qlfact = aln->llmult = 1;
1344   aln->llrev = aln->qlrev = 0;
1345   aln->frame = 0;
1346 #endif
1347 }
1348
1349 #include "a_mark.h"
1350
1351 int calcons(const unsigned char *aa0, int n0,
1352             const unsigned char *aa1, int n1,
1353             int *nc,
1354             struct a_struct *aln,
1355             struct a_res_str a_res,
1356             struct pstruct pst,
1357             char *seqc0, char *seqc1, char *seqca,
1358             struct f_struct *f_str)
1359 {
1360   int i0, i1, nn1, n0t;
1361   int op, lenc, len_gap, nd, ns, itmp;
1362   const unsigned char *aa1p;
1363   char *sp0, *sp1, *sq, *spa;
1364   int *rp;
1365   int mins, smins;
1366
1367   /* do not allow low complexity */
1368   sq = pst.sq;
1369
1370 #ifndef TFAST
1371   aa1p = aa1;
1372   nn1 = n1;
1373 #else
1374   aa1p = f_str->aa1x;
1375   nn1 = f_str->n10;
1376 #endif
1377
1378   /* first fill in the ends */
1379   /*   a_res.min0--; a_res.min1--; */
1380   n0 -= (f_str->nm0-1);
1381
1382   aln->amin0 = a_res.min0;
1383   aln->amin1 = a_res.min1;
1384   aln->amax0 = a_res.max0;
1385   aln->amax1 = a_res.max1;
1386
1387   if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1) {
1388     /* will we show all the start ?*/
1389     smins=0;
1390     mins = min(a_res.min1,aln->llen/2);
1391     aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
1392     aln->smin1 = a_res.min1-mins;
1393     if ((mins-a_res.min0)>0) {
1394       memset(seqc0,' ',mins-a_res.min0);
1395       aancpy(seqc0+mins-a_res.min0,(char *)f_str->aa0t,a_res.min0,pst);
1396       aln->smin0 = 0;
1397     }
1398     else {
1399       aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
1400       aln->smin0 = a_res.min0-mins;
1401     }
1402   }
1403   else {
1404     mins= min(aln->llen/2,min(a_res.min0,a_res.min1));
1405     smins=mins;
1406     aln->smin0=a_res.min0;
1407     aln->smin1=a_res.min1;
1408     aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
1409     aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
1410   }
1411
1412   memset(seqca,M_BLANK,mins);
1413
1414 /* now get the middle */
1415
1416   spa = seqca+mins;
1417   sp0 = seqc0+mins;
1418   sp1 = seqc1+mins;
1419   rp = a_res.res;
1420   n0t = lenc = len_gap = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = op = 0;
1421   i0 = a_res.min0;
1422   i1 = a_res.min1;
1423
1424   while (i0 < a_res.max0 || i1 < a_res.max1) {
1425     if (op == 0 && *rp == 0) {
1426       op = *rp++;
1427
1428       if ((itmp=pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
1429       else if (itmp == 0) { *spa = M_ZERO;}
1430       else {*spa = M_POS;}
1431       if (*spa == M_ZERO || *spa == M_POS) { aln->nsim++;}
1432
1433       *sp0 = sq[f_str->aa0t[i0++]];
1434       *sp1 = sq[aa1p[i1++]];
1435       n0t++;
1436       lenc++;
1437       if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
1438       sp0++; sp1++; spa++;
1439     }
1440     else {
1441       if (op==0) { op = *rp++;}
1442       if (op>0) {
1443         *sp0++ = '-';
1444         *sp1++ = sq[aa1p[i1++]];
1445         *spa++ = M_DEL;
1446         op--;
1447         len_gap++;
1448         lenc++;
1449       }
1450       else {
1451         *sp0++ = sq[f_str->aa0t[i0++]];
1452         *sp1++ = '-';
1453         *spa++ = M_DEL;
1454         op++;
1455         n0t++;
1456         len_gap++;
1457         lenc++;
1458       }
1459     }
1460   }
1461
1462   *spa = '\0';
1463   *nc = lenc-len_gap;
1464
1465   /* now we have the middle, get the right end */
1466   /* ns is amount to be shown */
1467   /* nd is amount remaining to be shown */
1468   ns = mins + lenc + aln->llen;
1469   ns -= (itmp = ns %aln->llen);
1470   if (itmp>aln->llen/2) ns += aln->llen;
1471   nd = ns - (mins+lenc);
1472   if (nd > max(n0t-a_res.max0,nn1-a_res.max1)) nd = max(n0t-a_res.max0,nn1-a_res.max1);
1473
1474   if (aln->showall==1) {
1475     nd = max(n0t-a_res.max0,nn1-a_res.max1);    /* reset for showall=1 */
1476     /* get right end */
1477     /* there isn't any aa0 to get */
1478     memset(seqc0+mins+lenc,' ',n0t-a_res.max0);
1479     aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
1480     /* fill with blanks - this is required to use one 'nc' */
1481     memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
1482     memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
1483   }
1484   else {
1485     memset(seqc0+mins+lenc,' ',nd);
1486     if ((nd-(nn1-a_res.max1))>0) {
1487       aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
1488       memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
1489     }
1490     else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
1491   }
1492
1493   return mins+lenc+nd;
1494 }
1495
1496 int calcons_a(const unsigned char *aa0, unsigned char *aa0a, int n0,
1497               const unsigned char *aa1, int n1,
1498               int *nc,
1499               struct a_struct *aln,
1500               struct a_res_str a_res,
1501               struct pstruct pst,
1502               char *seqc0, char *seqc0a, char *seqc1, char *seqca,
1503               char *ann_arr, struct f_struct *f_str)
1504 {
1505   int i0, i1, nn1, n0t;
1506   int op, lenc, len_gap, nd, ns, itmp;
1507   const unsigned char *aa1p;
1508   char *sp0, *sp0a, *sp1, *sq, *spa;
1509   int *rp;
1510   int mins, smins;
1511
1512   /* do not allow low complexity */
1513   sq = pst.sq;
1514
1515 #ifndef TFAST
1516   aa1p = aa1;
1517   nn1 = n1;
1518 #else
1519   aa1p = f_str->aa1x;
1520   nn1 = f_str->n10;
1521 #endif
1522
1523   aln->amin0 = a_res.min0;
1524   aln->amin1 = a_res.min1;
1525   aln->amax0 = a_res.max0;
1526   aln->amax1 = a_res.max1;
1527
1528   /* first fill in the ends */
1529   n0 -= (f_str->nm0-1);
1530
1531   if (min(a_res.min0,a_res.min1)<aln->llen || aln->showall==1) {
1532     /* will we show all the start ?*/
1533     smins=0;
1534     mins = min(a_res.min1,aln->llen/2);
1535     aancpy(seqc1,(char *)(aa1p+a_res.min1-mins),mins,pst);
1536     aln->smin1 = a_res.min1-mins;
1537     if ((mins-a_res.min0)>0) {
1538       memset(seqc0,' ',mins-a_res.min0);
1539       aancpy(seqc0+mins-a_res.min0,(char *)f_str->aa0t,a_res.min0,pst);
1540       aln->smin0 = 0;
1541     }
1542     else {
1543       aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
1544       aln->smin0 = a_res.min0-mins;
1545     }
1546   }
1547   else {
1548     mins= min(aln->llen/2,min(a_res.min0,a_res.min1));
1549     smins=mins;
1550     aln->smin0=a_res.min0;
1551     aln->smin1=a_res.min1;
1552     aancpy(seqc0,(char *)f_str->aa0t+a_res.min0-mins,mins,pst);
1553     aancpy(seqc1,(char *)aa1p+a_res.min1-mins,mins,pst);
1554   }
1555
1556   memset(seqca,M_BLANK,mins);
1557   memset(seqc0a,' ',mins);
1558
1559 /* now get the middle */
1560
1561   spa = seqca+mins;
1562   sp0 = seqc0+mins;
1563   sp0a = seqc0a+mins;
1564   sp1 = seqc1+mins;
1565   rp = a_res.res;
1566   n0t = lenc = len_gap = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = op = 0;
1567   i0 = a_res.min0;
1568   i1 = a_res.min1;
1569
1570   while (i0 < a_res.max0 || i1 < a_res.max1) {
1571     if (op == 0 && *rp == 0) {
1572       op = *rp++;
1573
1574       if ((itmp=pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]])<0) { *spa = M_NEG; }
1575       else if (itmp == 0) { *spa = M_ZERO;}
1576       else {*spa = M_POS;}
1577       if (*spa == M_ZERO || *spa == M_POS) { aln->nsim++;}
1578
1579       *sp0a++ = ' ';
1580       *sp0 = sq[f_str->aa0t[i0++]];
1581       *sp1 = sq[aa1p[i1++]];
1582       n0t++;
1583       lenc++;
1584       if (toupper(*sp0) == toupper(*sp1)) {aln->nident++; *spa = M_IDENT;}
1585       sp0++; sp1++; spa++;
1586     }
1587     else {
1588       if (op==0) { op = *rp++;}
1589       if (op>0) {
1590         *sp0++ = '-';
1591         *sp0a++ = ' ';
1592         *sp1++ = sq[aa1p[i1++]];
1593         *spa++ = M_DEL;
1594         op--;
1595         len_gap++;
1596         lenc++;
1597       }
1598       else {
1599         *sp0++ = sq[f_str->aa0t[i0++]];
1600         *sp0a++ = ' ';
1601         *sp1++ = '-';
1602         *spa++ = M_DEL;
1603         op++;
1604         n0t++;
1605         len_gap++;
1606         lenc++;
1607       }
1608     }
1609   }
1610
1611   *sp0a = *spa = '\0';
1612   *nc = lenc-len_gap;
1613
1614   /* now we have the middle, get the right end */
1615   /* ns is amount to be shown */
1616   /* nd is amount remaining to be shown */
1617   ns = mins + lenc + aln->llen;
1618   ns -= (itmp = ns %aln->llen);
1619   if (itmp>aln->llen/2) ns += aln->llen;
1620   nd = ns - (mins+lenc);
1621   if (nd > max(n0t-a_res.max0,nn1-a_res.max1)) nd = max(n0t-a_res.max0,nn1-a_res.max1);
1622
1623   if (aln->showall==1) {
1624     nd = max(n0t-a_res.max0,nn1-a_res.max1);    /* reset for showall=1 */
1625     /* get right end */
1626     /* there isn't any aa0 to get */
1627     memset(seqc0+mins+lenc,' ',n0t-a_res.max0);
1628     aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
1629     /* fill with blanks - this is required to use one 'nc' */
1630     memset(seqc0+mins+lenc+n0t-a_res.max0,' ',nd-(n0t-a_res.max0));
1631     memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
1632   }
1633   else {
1634     memset(seqc0+mins+lenc,' ',nd);
1635     if ((nd-(nn1-a_res.max1))>0) {
1636       aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nn1-a_res.max1,pst);
1637       memset(seqc1+mins+lenc+nn1-a_res.max1,' ',nd-(nn1-a_res.max1));
1638     }
1639     else aancpy(seqc1+mins+lenc,(char *)aa1p+a_res.max1,nd,pst);
1640   }
1641
1642   return mins+lenc+nd;
1643 }
1644
1645 void aa0shuffle(unsigned char *aa0, int n0, struct f_struct *f_str) {
1646
1647   int i, j, k;
1648   unsigned char tmp;
1649
1650   for (i = f_str->nmoff-1 ; --i ; ) {
1651
1652     /* j = nrand(i); if (i == j) continue;*/       /* shuffle columns */
1653     j = (f_str->nmoff - 2) - i; if (i <= j) break; /* reverse columns */
1654
1655     /* swap all i'th column residues for all j'th column residues */
1656     for(k = 0 ; k < f_str->nm0 ; k++) {
1657       tmp = aa0[(k * (f_str->nmoff)) + i];
1658       aa0[(k * (f_str->nmoff)) + i] = aa0[(k * (f_str->nmoff)) + j];
1659       aa0[(k * (f_str->nmoff)) + j] = tmp;
1660     }
1661   }
1662 }
1663
1664 /* build an array of match/ins/del - length strings */
1665 int calc_code(const unsigned char *aa0, const int n0,
1666               const unsigned char *aa1, const int n1,
1667               struct a_struct *aln,
1668               struct a_res_str a_res,
1669               struct pstruct pst,
1670               char *al_str, int al_str_n, struct f_struct *f_str)
1671 {
1672   int i0, i1, nn1;
1673   int op, lenc, len_gap;
1674   int p_op, op_cnt;
1675   const unsigned char *aa1p;
1676   char tmp_cnt[20];
1677   char sp0, sp1, *sq;
1678   int *rp;
1679   int mins, smins;
1680   int fnum = 0;
1681
1682   if (pst.ext_sq_set) {
1683     sq = pst.sqx;
1684   }
1685   else {
1686     sq = pst.sq;
1687   }
1688
1689 #ifndef TFAST
1690   aa1p = aa1;
1691   nn1 = n1;
1692 #else
1693   aa1p = f_str->aa1x;
1694   nn1 = f_str->n10;
1695 #endif
1696
1697   rp = a_res.res;
1698   lenc = len_gap = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = aln->nfs = op = p_op = 0;
1699   op_cnt = 0;
1700
1701   i0 = a_res.min0;
1702   i1 = a_res.min1;
1703   tmp_cnt[0]='\0';
1704
1705   fnum = f_str->aa0ti[i0] + 1;
1706   while (i0 < a_res.max0 || i1 < a_res.max1) {
1707     if (op == 0 && *rp == 0) {
1708       if (p_op == 0) {  op_cnt++;}
1709       else {
1710         update_code(al_str,al_str_n-strlen(al_str),p_op,op_cnt,fnum);
1711         op_cnt = 1; p_op = 0;
1712         fnum = f_str->aa0ti[i0] + 1;
1713       }
1714       op = *rp++;
1715       lenc++;
1716       if (pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]]>=0) {aln->nsim++;}
1717       sp0 = pst.sq[f_str->aa0t[i0++]];
1718       sp1 = pst.sq[aa1p[i1++]];
1719       if (toupper(sp0) == toupper(sp1)) aln->nident++;
1720     }
1721     else {
1722       if (op==0) op = *rp++;
1723       if (op>0) {
1724         if (p_op == 1) { op_cnt++;}
1725         else {
1726           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);
1727           op_cnt = 1; p_op = 1; fnum = f_str->aa0ti[i0] + 1;
1728         }
1729         op--; lenc++; i1++; len_gap++;
1730       }
1731       else {
1732         if (p_op == 2) { op_cnt++;}
1733         else {
1734           update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);
1735           op_cnt = 1; p_op = 2; fnum = f_str->aa0ti[i0] + 1;
1736         }
1737         op++; lenc++; i0++; len_gap++;
1738       }
1739     }
1740   }
1741   update_code(al_str,al_str_n - strlen(al_str),p_op,op_cnt,fnum);
1742
1743   return lenc - len_gap;
1744 }
1745
1746 void
1747 update_code(char *al_str, int al_str_max, int op, int op_cnt, int fnum) {
1748
1749   char op_char[4]={"=-+"};
1750   char tmp_cnt[20];
1751
1752   if (op == 0)
1753     sprintf(tmp_cnt,"%c%d[%d]",op_char[op],op_cnt,fnum);
1754   else
1755     sprintf(tmp_cnt,"%c%d",op_char[op],op_cnt);
1756
1757   strncat(al_str,tmp_cnt,al_str_max);
1758 }
1759
1760 int calc_id(const unsigned char *aa0, int n0,
1761             const unsigned char *aa1, int n1,
1762             struct a_struct *aln,
1763             struct a_res_str a_res,
1764             struct pstruct pst,
1765             struct f_struct *f_str)
1766 {
1767   int i0, i1, nn1, n0t;
1768   int op, lenc, len_gap;
1769   const unsigned char *aa1p;
1770   int sp0, sp1;
1771   int *rp;
1772   int mins, smins;
1773   char *sq;
1774
1775   if (pst.ext_sq_set) {
1776     sq = pst.sqx;
1777   }
1778   else {
1779     sq = pst.sq;
1780   }
1781
1782 #ifndef TFAST
1783   aa1p = aa1;
1784   nn1 = n1;
1785 #else
1786   aa1p = f_str->aa1x;
1787   nn1 = f_str->n10;
1788 #endif
1789
1790   /* first fill in the ends */
1791   /* a_res.min0--; a_res.min1--; */
1792   n0 -= (f_str->nm0-1);
1793
1794   /* now get the middle */
1795   rp = a_res.res;
1796   n0t = lenc = len_gap = aln->nident = aln->nsim = aln->ngap_q = aln->ngap_l = op = 0;
1797   i0 = a_res.min0;
1798   i1 = a_res.min1;
1799
1800   while (i0 < a_res.max0 || i1 < a_res.max1) {
1801     if (op == 0 && *rp == 0) {
1802       op = *rp++;
1803       if (pst.pam2[0][f_str->aa0t[i0]][aa1p[i1]]>=0) {aln->nsim++;}
1804       sp0 = sq[f_str->aa0t[i0++]];
1805       sp1 = sq[aa1p[i1++]];
1806       n0t++;
1807       lenc++;
1808       if (toupper(sp0) == toupper(sp1)) aln->nident++;
1809     }
1810     else {
1811       if (op==0) { op = *rp++;}
1812       if (op>0) {
1813         i1++;
1814         op--;
1815         len_gap++;
1816         lenc++;
1817       }
1818       else {
1819         i0++;
1820         op++;
1821         n0t++;
1822         len_gap++;
1823         lenc++;
1824       }
1825     }
1826   }
1827   return lenc-len_gap;
1828 }
1829
1830 #ifdef PCOMPLIB
1831
1832 #include "structs.h"
1833 #include "p_mw.h"
1834
1835 void
1836 update_params(struct qmng_str *qm_msg,
1837               struct mngmsg *m_msg, struct pstruct *ppst)
1838 {
1839   m_msg->n0 = ppst->n0 = qm_msg->n0;
1840   m_msg->nm0 = qm_msg->nm0;
1841   m_msg->escore_flg = qm_msg->escore_flg;
1842   m_msg->qshuffle = qm_msg->qshuffle;
1843 }
1844 #endif