dropff2.c

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

	       dptr->start = (dptr->stop = lpos);
	     }
	 }
	 else {				/* continue current run */
	   dptr->score += f_str->pamh1[aa0[tpos]]; 
	   dptr->stop = lpos;
	 }
       }
       else {				/* no diagonal run yet */
	 dptr->score = f_str->pamh2[lhval];
	 dptr->start = (dptr->stop = lpos);
       }
     }				/* end tpos */
   }				/* end lpos */

   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;

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

   /* set up pointers for sorting */

   for (nsave = 0, vmptr = f_str->vmax; vmptr < &f_str->vmax[MAXSAV]; vmptr++) {
     if (vmptr->score > 0) {
       vmptr->score = m0_spam (aa0, aa1, n1, vmptr, ppst->pam2[0], f_str);
       f_str->vptr[nsave++] = vmptr;
     }
   }

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

   
#ifdef DEBUG
   /*
   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);
       for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
         fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
		 ppst->sq[aa1[im]]);
       fputc('\n',stderr);
   }
   fprintf(stderr,"---\n");
   */
   /* now use m_spam to re-evaluate */
   /*
   for (tpos = 0; tpos < n0; tpos++) {
     fprintf(stderr,"%c:%2d ",ppst->sq[aa0[tpos]],aa0[tpos]);
     if (tpos %10 == 9) fputc('\n',stderr);
   }
   fputc('\n',stderr);
   */
#endif   
   
   f_str->aa0ix = 0;
   for (ib=0; ib < nsave; ib++) {
     if ((vmptr=f_str->vptr[ib])->score > 0) {
       vmptr->score = m1_spam (aa0, n0, aa1, n1, vmptr,
			       ppst->pam2[0], ppst->pam_l, f_str);
     }
   }
   /* reset aa0 - modified by m1_spam */
   for (tpos = 0; tpos < n0; tpos++) {
     if (aa0[tpos] >= 32) aa0[tpos] -= 32;
   }

   kssort(f_str->vptr,nsave);

   for ( ; nsave > 0; nsave--) 
     if (f_str->vptr[nsave-1]->score >0) break;

   if (nsave <= 0) {
     f_str->nsave = 0;
     rst->score[0] = rst->score[1] = rst->score[2] = 0;
     rst->escore = 1.0;

     return 1;
   }
   else f_str->nsave = nsave;

   
#ifdef DEBUG
   /*
   fprintf(stderr,"n0: %d; n1: %d; noff: %d\n",n0,n1,f_str->noff);
   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);
     for (im=f_str->vptr[ib]->start; im<=f_str->vptr[ib]->stop; im++)
       fprintf(stderr," %c:%c",ppst->sq[aa0[f_str->noff+im-f_str->vptr[ib]->dp]],
  	       ppst->sq[aa1[im]]);
     fputc('\n',stderr);
   }

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

   scor = sconn (f_str->vptr, nsave, ppst->param_u.fa.cgap, f_str,
		 rst, ppst, aa0, n0, aa1, n1, opt_prob);

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

   rst->score[1] = vmptr->score;
   rst->score[0] = rst->score[2] = max (scor, vmptr->score);

   return 1;
}

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 opt_prob;
  int hoff, n10, i;

  if (qr_flg==1 && f_str->shuff_cnt <= 0) {
    rst->escore = 2.0;
    rst->score[0]=rst->score[1]=rst->score[2]= -1;
    return;
  }

  if (f_str->dotat || ppst->zsflag == 4 || ppst->zsflag == 14 ) opt_prob=1;
  else opt_prob = 0;
  if (ppst->zsflag == 2 || ppst->zsflag == 12) opt_prob = 0;
  if (qr_flg) {
    opt_prob=1;
    /*    if (frame==1) */
      f_str->shuff_cnt--;
  }

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

#ifdef TFAST 
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  if (ppst->debug_lib)
    for (i=0; i<n10; i++)
      if (f_str->aa1x[i]>ppst->nsq) {
	fprintf(stderr,
		"residue[%d/%d] %d range (%d)\n",i,n1,
		f_str->aa1x[i],ppst->nsq);
	f_str->aa1x[i]=0;
	n10=i-1;
      }

  do_fastf (f_str->aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, opt_prob);
#else	/* FASTF */
  do_fastf (f_str->aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, opt_prob);
#endif

  rst->comp = rst->H = -1.0;

}

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, n10;

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

#ifdef TFAST  
  n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fastf (f_str->aa0, n0, f_str->aa1x, n10, ppst, f_str, rst, &hoff, 1);
#else	/* FASTA */
  do_fastf(f_str->aa0, n0, aa1, n1, ppst, f_str, rst, &hoff, 1);
#endif
  ppst->param_u.fa.optflag = optflag;
}

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);

/* 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;
}

/* this version of spam() is designed to work with a collection of
   subfragments, selecting the best amino acid at each position so
   that, from each subfragment, each position is only used once.

   As a result, m_spam needs to know the number of fragments.

   In addition, it now requires a global alignment to the fragment
   and resets the start and stop positions

   */

static int
m1_spam (unsigned char *aa0, int n0,
	 const unsigned char *aa1, int n1,
	 struct savestr *dmax, int **pam2, int pam_l,
	 struct f_struct *f_str)
{
  int     tpos, lpos, im, ii, nm, ci;
  int     tot, ctot, pv;

  struct {
    int     start, stop, score;
  } curv, maxv;
  unsigned char *aa0p;
  const unsigned char *aa1p;

  lpos = dmax->start;                   /* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   /* force global alignment, reset start*/
   if (tpos < lpos) {
     lpos = dmax->start -= tpos;
     tpos = 0;
   }
   else {
     tpos -= lpos;
     lpos = dmax->start = 0;
   }

   dmax->stop = dmax->start + (f_str->nmoff -2 - tpos);
   if (dmax->stop > n1) dmax->stop = n1;

   /*
   if (dmax->start < 0) {
     tpos = -dmax->start;
     lpos = dmax->start=0;
   }
   else tpos = 0;
   */

   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];

   nm = f_str->nm0;

   tot = curv.score = maxv.score = 0;
   for (; lpos <= dmax->stop; lpos++,aa0p++,aa1p++) {
     ctot = pam_l;
     ci = -1;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
	 ci = ii;
/*  	 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]); */
       }
     }
     tot += ctot;
     if (ci >= 0 && aa0p[ci] < 32) {
#ifdef DEBUG
/*         fprintf(stderr, "used: lpos: %d ci: %d : %c\n", lpos, ci, sq[aa0p[ci]]); */
#endif
       aa0p[ci] +=  32;
       dmax->used[&aa0p[ci] - aa0] = 1;
     }
   }
   return tot;
}

int ma_spam (unsigned char *aa0, int n0, const unsigned char *aa1,
	     struct savestr *dmax, struct pstruct *ppst,
	     struct f_struct *f_str)
{
  int **pam2;
  int     tpos, lpos, im, ii, nm, ci, lp0;
  int     tot, ctot, pv;
  struct {
    int     start, stop, score;
  } curv, maxv;
   const unsigned char *aa1p;
   unsigned char *aa0p, *aa0pt;
   int aa0t_flg;

   pam2 = ppst->pam2[0];
   aa0t_flg = 0;

   lpos = dmax->start;			/* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   lp0 = lpos = dmax->start;
   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];			/* real aa0 sequence */

			/* the destination aa0 sequence (without nulls) */
   aa0pt = &f_str->aa0t[f_str->aa0ix];

   curv.start = lpos;
   nm = f_str->nm0;

   /* sometimes, tpos may be > 0, with lpos = 0 - fill with 'X' */
   if (lpos == 0 && tpos > 0)
     for (ii = 0; ii < tpos; ii++) *aa0pt++ = 31;  /* filler character */

   tot = curv.score = maxv.score = 0;
   for (; lpos <= dmax->stop; lpos++) {
     ctot = ppst->pam_l;
     ci = -1;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if (aa0p[ii] < 32 && (pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
	 ci = ii;
       }
     }
     tot += ctot;
     if (ci >= 0) {
       if (ci >= n0) {fprintf(stderr," warning - ci off end %d/%d\n",ci,n0);}
       else {
	 *aa0pt++ = aa0p[ci];
	 aa0p[ci] +=  32;
	 aa0t_flg=1;
       }
     }
     aa0p++; aa1p++;
   }

   if (aa0t_flg) {
     dmax->dp -= f_str->aa0ix;		/* shift ->dp for aa0t */
     if ((ci=(int)(aa0pt-f_str->aa0t)) > n0) {
       fprintf(stderr," warning - aapt off %d/%d end\n",ci,n0);
     }
     else 
       *aa0pt++ = 0;			/* skip over NULL */

     aa0pt = &f_str->aa0t[f_str->aa0ix];
     aa1p = &aa1[lp0];

     /*
     for (im = 0; im < f_str->nmoff; im++)
       fprintf(stderr,"%c:%c,",ppst->sq[aa0pt[im]],ppst->sq[aa1p[im]]);
     fprintf(stderr,"- %3d (%3d:%3d)\n",dmax->score,f_str->aa0ix,lp0);
     */

     f_str->aa0ix += f_str->nmoff;	/* update offset into aa0t */
   }
   /*
      fprintf(stderr," ma_spam returning: %d\n",tot);
   */
   return tot;
}

static int
m0_spam (unsigned char *aa0, const unsigned char *aa1, int n1,
	 struct savestr *dmax, int **pam2,
	 struct f_struct *f_str)
{
   int tpos, lpos, lend, im, ii, nm;
   int     tot, ctot, pv;
   struct {
     int     start, stop, score;
   } curv, maxv;
   const unsigned char *aa0p, *aa1p;

   lpos = dmax->start;			/* position in library sequence */
   tpos = lpos - dmax->dp + f_str->noff; /* position in query sequence */
   if (tpos > 0) {
     if (lpos-tpos >= 0) {
       lpos = dmax->start -= tpos;    /* force global alignment, reset start*/
       tpos = 0;
     }
     else {
       tpos -= lpos;
       lpos = dmax->start = 0;
     }
   }

   nm = f_str->nm0;
   lend = dmax->stop;
   if (n1 - (lpos + f_str->nmoff-2) < 0 ) {
     lend = dmax->stop = (lpos - tpos) + f_str->nmoff-2;
     if (lend >= n1) lend = n1-1;
   }

   aa1p = &aa1[lpos];
   aa0p = &aa0[tpos];

   curv.start = lpos;

   tot = curv.score = maxv.score = 0;
   for (; lpos <= lend; lpos++) {
     ctot = -10000;
     for (im = 0, ii=0; im < nm; im++,ii+=f_str->nmoff) {
       if ((pv = pam2[aa0p[ii]][*aa1p]) > ctot) {
	 ctot = pv;
       }
     }
     tot += ctot;
     aa0p++; aa1p++;
   }

   /* reset dmax if necessary */

   return tot;
}

/* sconn links up non-overlapping alignments and calculates the score */

int sconn (struct savestr **v, int n, int cgap, struct f_struct *f_str,
	   struct rstruct *rst, struct pstruct *ppst,
	   const unsigned char *aa0, int n0,
	   const unsigned char *aa1, int n1,
	   int opt_prob)
{
   int     i, si, cmpp ();
   struct slink *start, *sl, *sj, *so, sarr[MAXSAV];
   int     lstart, plstop;
   double tatprob;