dropfx.c

序列对齐 Compare a protein sequence to a protein sequence database or a DNA sequence to a DNA sequenc

    free(f_str->harr);
    free(f_str);
    *f_arg = NULL;
  }
}

void do_fastx (const unsigned char *aa0, int n0,
	       const unsigned char *aa1, int n1,
	       struct pstruct *ppst, struct f_struct *f_str,
	       struct rstruct *rst, int *hoff)
{
   int     nd;		/* diagonal array size */
   int     lhval;
   int     kfact;
   int i;
   int my_hoff;
   register struct dstruct *dptr;
   register int tscor;

#ifndef ALLOCN0
   register struct dstruct *diagp;
#else
   register int dpos;
   int     lposn0;
#endif
   struct dstruct *dpmax;
   register int lpos;
   int     tpos;
   struct savestr *vmptr;
   int     scor, tmp;
   int     im, ib, nsave;
   int ktup, kt1, *hsq, ip, lkt;
#ifndef TFAST
   int n0x31, n0x32;
   n0x31 = (n0-2)/3;
   n0x32 = n0x31+1+(n0-n0x31-1)/2;
#else
   const unsigned char *fs;
   unsigned char *fd;
   int n1x31, n1x32, last_n1, itemp;
   n1x31 = (n1-2)/3;
   n1x32 = n1x31+1+(n1-n1x31-1)/2;
#endif

  if (ppst->ext_sq_set) {
    ip = 1;
    hsq = ppst->hsqx;
  }
  else {
    ip = 0;
    hsq = ppst->hsq;
  }

   ktup = ppst->param_u.fa.ktup;
   kt1 = ktup-1;

   if (n1 < ktup) {
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     return;
   }

   if (n0+n1+1 >= MAXDIAG) {
     fprintf(stderr,"n0,n1 too large: %d, %d\n",n0,n1);
     rst->score[0] = rst->score[1] = rst->score[2] = -1;
     return;
   }

   f_str->noff = n0 - 1;

#ifdef ALLOCN0
   nd = n0;
#endif

#ifndef ALLOCN0
   nd = n0 + n1;
#endif

   dpmax = &f_str->diag[nd];
   for (dptr = &f_str->diag[f_str->ndo]; dptr < dpmax;)
   {
      dptr->stop = -1;
      dptr->dmax = NULL;
      dptr++->score = 0;
   }

   for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
      vmptr->score = 0;
   f_str->lowmax = f_str->vmax;
   f_str->lowscor = 0;

   /* start hashing */
   lhval = 0;
   lkt = kt1;
   for (lpos = 0; (lpos < lkt || hsq[aa1[lpos]]>=NMAP) && lpos<n1; lpos++) {
     if (hsq[aa1[lpos]]>=NMAP) {
       lhval = 0; lkt=lpos+ktup; continue;
#ifdef ALLOCN0		/* reinitialize dptr */
       dptr = &f_str->diag[lpos % nd];
       dptr->stop = -1;
       dptr->dmax = NULL;
       dptr->score = 0;
#endif
     }
     lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
   }

#ifndef ALLOCN0
   diagp = &f_str->diag[f_str->noff + lkt];
   for (; lpos < n1; lpos++, diagp++) {
     /*     if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; continue;} */
     if (hsq[aa1[lpos]]>=NMAP) {
       lpos++ ; diagp++;
       while (lpos < n1 && hsq[aa1[lpos]]>=NMAP) {lpos++; diagp++;}
       if (lpos >= n1) break;
       lhval = 0;
     }
     lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
     for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
       if ((tscor = (dptr = &diagp[-tpos])->stop) >= 0) {
#else
   lposn0 = f_str->noff + lpos;
   for (; lpos < n1; lpos++, lposn0++) {
     if (hsq[aa1[lpos]]>=NMAP) {lhval = 0; goto loopl;}
     lhval = ((lhval & f_str->hmask) << f_str->kshft) + hsq[aa1[lpos]];
     for (tpos = f_str->harr[lhval]; tpos >= 0; tpos = f_str->link[tpos]) {
       dpos = lposn0 - tpos;
       if ((tscor = (dptr = &f_str->diag[dpos % nd])->stop) >= 0) {
#endif
	 tscor += ktup;
	 if ((tscor -= lpos) <= 0) {	/* better to start over */
	   scor = dptr->score;
	   if ((tscor += (kfact = f_str->pamh2[lhval])) < 0 && f_str->lowscor < scor)
#ifdef ALLOCN0
	     savemax (dptr, dpos, f_str);
#else
	     savemax (dptr, f_str);
#endif
	     if ((tscor += scor) >= kfact) {
	       dptr->score = tscor;
	       dptr->stop = lpos;
	     }
	     else {
	       dptr->score = kfact;
	       dptr->start = (dptr->stop = lpos) - kt1;
	     }
	 }				/* continue current run in diagonal */
	 else {
	   dptr->score += f_str->pamh1[aa0[tpos]];
	   dptr->stop = lpos;
	 }
       }
       else {
	 dptr->score = f_str->pamh2[lhval];
	 dptr->start = (dptr->stop = lpos) - kt1;
       }
     }				/* end tpos */

#ifdef ALLOCN0
      /* reinitialize diag structure */
   loopl:
     if ((dptr = &f_str->diag[lpos % nd])->score > f_str->lowscor) {
	   savemax (dptr, lpos, f_str);
     }
     dptr->stop = -1;
     dptr->dmax = NULL;
     dptr->score = 0;
#endif
   }				/* end lpos */

#ifdef ALLOCN0
   for (tpos = 0, dpos = f_str->noff + n1 - 1; tpos < n0; tpos++, dpos--) {
     if ((dptr = &f_str->diag[dpos % nd])->score > f_str->lowscor)
       savemax (dptr, dpos, f_str);
   }
#else
   for (dptr = f_str->diag; dptr < dpmax;) {
     if (dptr->score > f_str->lowscor) savemax (dptr, f_str);
     dptr->stop = -1;
     dptr->dmax = NULL;
     dptr++->score = 0;
   }
   f_str->ndo = nd;
#endif

/*
        at this point all of the elements of aa1[lpos]
        have been searched for elements of aa0[tpos]
        with the results in diag[dpos]
*/

   for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
   {
     /*
       fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	       f_str->noff+vmptr->start-vmptr->dp,
	       f_str->noff+vmptr->stop-vmptr->dp,
	       vmptr->start,vmptr->stop,
	       vmptr->dp,vmptr->score);
     */
      if (vmptr->score > 0) {
	 vmptr->score = spam (aa0, aa1, vmptr, ppst->pam2[ip], f_str);
	 f_str->vptr[nsave++] = vmptr;
      }
   }

   if (nsave <= 0) {
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     return;
   }
       
#ifndef TFAST
   /* FASTX code here to modify the start, stop points for 
      the three phases of the translated protein sequence
      */
   /*
     fprintf(stderr,"n0x: %d; n0x31:%d; n0x32: %d\n",n0,n0x31,n0x32);
     for (ib=0; ib<nsave; ib++) {
       fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	       f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	       f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
	       f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
	       f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
     }

     fprintf(stderr,"---\n");
   */
   for (ib=0; ib<nsave; ib++) {
     if (f_str->noff-f_str->vptr[ib]->dp+f_str->vptr[ib]->start >= n0x32)
       f_str->vptr[ib]->dp += n0x32;
     if (f_str->noff-f_str->vptr[ib]->dp +f_str->vptr[ib]->start >= n0x31)
       f_str->vptr[ib]->dp += n0x31;
   }
	    
   /*
     for (ib=0; ib<nsave; ib++) {
       fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	       f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	       f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
	       f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
	       f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
     }
     */
#else

   /* TFASTX code here to modify the start, stop points for 
	     the three phases of the translated protein sequence
	     TFASTX modifies library start points, rather than 
	     query start points
	     */

     /*
   fprintf(stderr,"n0: %d; noff: %d; n1: %d; n1x31: %d n1x32 %d\n",n0, f_str->noff,n1,n1x31,n1x32);
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	     f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	     f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
	     f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
	     f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
   }

   fprintf(stderr,"---\n");
   */

   for (ib=0; ib<nsave; ib++) {
     if (f_str->vptr[ib]->start >= n1x32) {
       f_str->vptr[ib]->start -= n1x32;
       f_str->vptr[ib]->stop -= n1x32;
       f_str->vptr[ib]->dp -= n1x32;
     }
     if (f_str->vptr[ib]->start >= n1x31) {
       f_str->vptr[ib]->start -= n1x31;
       f_str->vptr[ib]->stop -= n1x31;
       f_str->vptr[ib]->dp -= n1x31;
     }
   }
	    
   /*
   for (ib=0; ib<nsave; ib++) {
     fprintf(stderr,"0: %4d-%4d  1: %4d-%4d  dp: %d score: %d\n",
	     f_str->noff+f_str->vptr[ib]->start-f_str->vptr[ib]->dp,
	     f_str->noff+f_str->vptr[ib]->stop-f_str->vptr[ib]->dp,
	     f_str->vptr[ib]->start,f_str->vptr[ib]->stop,
	     f_str->vptr[ib]->dp,f_str->vptr[ib]->score);
   }
   */

#endif /* TFASTX */

   scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap, 
		 ppst->param_u.fa.pgap, f_str);

   for (vmptr=f_str->vptr[0],ib=1; ib<nsave; ib++)
     if (f_str->vptr[ib]->score > vmptr->score) vmptr=f_str->vptr[ib];

/*  kssort (f_str->vptr, nsave); */

   rst->score[1] = vmptr->score;	/* best single score - init1*/
   rst->score[0] = max (scor, vmptr->score);  /* initn */
   rst->score[2] = rst->score[0];		/* initn */

   my_hoff=f_str->noff - vmptr->dp;

   /*
   if (n1 > 5000) {
     fprintf(stderr," Long n1: %d\n",n1);
   }
   */

   if (ppst->param_u.fa.optflag) {
     if (rst->score[0] > ppst->param_u.fa.optcut) {
#ifndef TFAST
       rst->score[2] = dmatchx(aa0, n0,aa1,n1,my_hoff,
			     ppst->param_u.fa.optwid, ppst->pam2[ip],
			     ppst->gdelval,ppst->ggapval,ppst->gshift,f_str);
#else /* TFASTX */
     /* generate f_str->aa1y */
/*
     for (i=0; i<n1; i++) {
       fputc(ppst->sq[aa1[i]],stderr);
       if (i%60==59) fputc('\n',stderr);
     }
     fprintf(stderr,"\n-----\n");
*/
     for (fs=aa1,itemp=0; itemp <3; itemp++,fs++) {
       for (fd= &f_str->aa1y[itemp]; *fs!=EOSEQ; fd += 3, fs++) *fd = *fs;
       *fd=EOSEQ;
     }

/*
     for (i=0; i<n1; i++) {
       fputc(ppst->sq[f_str->aa1y[i]],stderr);
       if (i%60==59) fputc('\n',stderr);
     }
*/
     rst->score[2] = dmatchx(aa0, n0, aa1, n1, my_hoff=vmptr->dp-f_str->noff,
			     ppst->param_u.fa.optwid, ppst->pam2[ip],
			     ppst->gdelval,ppst->ggapval,ppst->gshift,f_str);
#endif /* TFASTX */
     }
   }
   *hoff = my_hoff;
}

/* returns rst.score[0] - initn
   	   rst.score[1] - init1
	   rst.score[2] - opt
*/

void do_work (const unsigned char *aa0, int n0,
	      const unsigned char *aa1, int n1,
	      int frame,
	      struct pstruct *ppst, struct f_struct *f_str,
	      int qr_flg, struct rstruct *rst)
{
  int hoff;
  int last_n1, itx, itt, n10, i;

#ifdef TFAST
  unsigned char *aa1x;
  /* aa0 has a protein sequence */
  /* aa1 has a raw DNA sequence */

  itt = frame;
  last_n1 = 0;
  aa1x = f_str->aa1x;
  for (itx= itt*3; itx< itt*3+3; itx++) {
    n10  = saatran(aa1,&aa1x[last_n1],n1,itx);
    /*
    fprintf(stderr," itt %d itx: %d\n",itt,itx);
    for (i=0; i<n10; i++) {
      fprintf(stderr,"%c",aa[f_str->aa1x[last_n1+i]]);
      if ((i%60)==59) fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
    */
    last_n1 += n10+1;
  }
  n10 = last_n1-1;
#endif

  rst->score[0] = rst->score[1] = rst->score[2] = 0;
  rst->escore = 1.0;
  rst->segnum = rst->seglen = 1;

#ifndef TFAST
  do_fastx (f_str->aa0x, n0, aa1, n1, ppst, f_str, rst, &hoff);
#else /* tfastx */
  do_fastx (aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff);
#endif
}

void do_opt (const unsigned char *aa0, int n0,
	     const unsigned char *aa1, int n1,
	     int frame,
	     struct pstruct *ppst,
	     struct f_struct *f_str,
	     struct rstruct *rst)
{
  int optflag, tscore, hoff;

  optflag = ppst->param_u.fa.optflag;
  ppst->param_u.fa.optflag = 1;

#ifndef TFAST
  do_fastx (f_str->aa0x, n0, aa1, n1, ppst, f_str, rst, &hoff);
#else
  do_fastx (aa0, n0, aa1, n1, ppst, f_str, rst, &hoff);
#endif

  ppst->param_u.fa.optflag = optflag;
}

#ifdef ALLOCN0
void
savemax (dptr, dpos, f_str)
  register struct dstruct *dptr;
  int  dpos;
  struct f_struct *f_str;
{
   register struct savestr *vmptr;
   register int i;

#else
void
savemax (dptr, f_str)
  register struct dstruct *dptr;
  struct f_struct *f_str;
{
   register int dpos;
   register struct savestr *vmptr;
   register int i;

   dpos = (int) (dptr - f_str->diag);

#endif

/* check to see if this is the continuation of a run that is already saved */

   if ((vmptr = dptr->dmax) != NULL && vmptr->dp == dpos &&
	 vmptr->start == dptr->start)
   {
      vmptr->stop = dptr->stop;
      if ((i = dptr->score) <= vmptr->score)
	 return;
      vmptr->score = i;
      if (vmptr != f_str->lowmax)
	 return;
   }
   else
   {
      i = f_str->lowmax->score = dptr->score;
      f_str->lowmax->dp = dpos;
      f_str->lowmax->start = dptr->start;
      f_str->lowmax->stop = dptr->stop;
      dptr->dmax = f_str->lowmax;
   }

   for (vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++)
      if (vmptr->score < i)
      {
	 i = vmptr->score;
	 f_str->lowmax = vmptr;
      }
   f_str->lowscor = i;
}

int spam (const unsigned char *aa0, const unsigned char *aa1,
	  struct savestr *dmax, int **pam2,
	  struct f_struct *f_str)
{
   int     lpos;
   int     tot, mtot;
   struct {
     int     start, stop, score;
   } curv, maxv;
   const unsigned char *aa0p, *aa1p;

   aa1p = &aa1[lpos = dmax->start];
   aa0p = &aa0[lpos - dmax->dp + f_str->noff];
   curv.start = lpos;

   tot = curv.score = maxv.score = 0;
   for (; lpos <= dmax->stop; lpos++) {
     tot += pam2[*aa0p++][*aa1p++];
     if (tot > curv.score) {
       curv.stop = lpos;
       curv.score = tot;
      }
      else if (tot < 0) {
	if (curv.score > maxv.score) {
	  maxv.start = curv.start;
	  maxv.stop = curv.stop;
	  maxv.score = curv.score;
	}
	tot = curv.score = 0;
	curv.start = lpos+1;
      }
   }

   if (curv.score > maxv.score) {
     maxv.start = curv.start;
     maxv.stop = curv.stop;
     maxv.score = curv.score;
   }