dropfs2.c

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

   for (i = 0, si = 0; i < n; i++) {

     /* the segment is worth adding; find out where? */
     lstart = v[i]->start;
     ltmp = v[i]->stop;
     tstart = lstart - v[i]->dp + f_str->noff;
     tstop = ltmp - v[i]->dp + f_str->noff;

     /*	put the run in the group */
     sarr[si].vp = v[i];
     sarr[si].score = v[i]->score;
     sarr[si].next = NULL;
     sarr[si].prev = NULL;
     sarr[si].tat = NULL;

/*
  opt_prob for FASTS only has to do with using aa1 for priors,
  i.e. we always calculate tatprobs for segments in FASTS (unlike
  FASTF)
*/
     if(opt_prob) {
       sarr[si].tatprob = 
	 calc_tatusov(NULL, &sarr[si], aa0, n0, aa1, n1, 
		      ppst->pam2[0], ppst->nsq, f_str, 
		      ppst->pseudocts, opt_prob, ppst->zsflag);
       if (sarr[si].tatprob < 0.0) {
	 fprintf(stderr," negative tatprob: %lg\n",sarr[si].tatprob);
	 sarr[si].tatprob = 1.0;
       }
       sarr[si].tat = sarr[si].newtat;
     }

/*  if it fits, then increase the score

    start points to the highest scoring run
    -> next is the second highest, etc.
    put the segment into the highest scoring run that it fits into
*/
     for (sl = start; sl != NULL; sl = sl->next) {
       ltmp = sl->vp->start;
 /* plstop -> previous lstop */
       plstop = sl->vp->stop;
 /* ptstart -> previous t(query) start */
       ptstart = ltmp - sl->vp->dp + f_str->noff;
 /* ptstop -> previous t(query) stop */
       ptstop = plstop - sl->vp->dp + f_str->noff;
#ifndef FASTM
 /* if the previous library stop is before the current library start */
       if (plstop < lstart && ( ptstop < tstart || ptstart > tstop))
#else
 /* if the previous library stop is before the current library start */
       if (plstop < lstart && ptstop < tstart)
#endif
       {
	 if(!opt_prob) {
	    sarr[si].score = sl->score + v[i]->score;
	    sarr[si].prev = sl;
	    break;
	  } else {
	    tatprob = calc_tatusov(sl, &sarr[si], aa0, n0, aa1, n1, 
				   ppst->pam2[0], ppst->nsq, f_str, 
				   ppst->pseudocts, opt_prob, ppst->zsflag);
	    /* if our tatprob gets worse when we add this, forget it */
	    if(tatprob > sarr[si].tatprob) {
	      free(sarr[si].newtat->probs); /* get rid of new tat struct */
	      free(sarr[si].newtat);
	      continue; /* reuse this sarr[si] */
	    } else {
	      sarr[si].tatprob = tatprob;
	      free(sarr[si].tat->probs); /* get rid of old tat struct */
	      free(sarr[si].tat);
	      sarr[si].tat = sarr[si].newtat;
	      sarr[si].prev = sl;
	      sarr[si].score = sl->score + v[i]->score;
	      /*
		fprintf(stderr,"sconn %d added %d/%d getting %d; si: %d, tat: %g\n",
		i,v[i]->start, v[i]->score,sarr[si].score,si, tatprob);
	      */
	      break;
	    }
	  }
	}
      }
      
      /* now recalculate where the score fits */
      if (start == NULL) start = &sarr[si];
      else {
	if(!opt_prob) {
	  for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	    if (sarr[si].score > sj->score) {
	      sarr[si].next = sj;
	      if (so != NULL)
		so->next = &sarr[si];
	      else
		start = &sarr[si];
	      break;
	    }
	    so = sj;
	  }
	} else {
	  for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	    if ( sarr[si].tatprob < sj->tatprob ||
		 ((sarr[si].tatprob == sj->tatprob) && sarr[si].score > sj->score) ) {
	      sarr[si].next = sj;
	      if (so != NULL)
		so->next = &sarr[si];
	      else
		start = &sarr[si];
	      break;
	    }
	    so = sj;
	  }
	}
      }

      si++;
   }
      
   if(opt_prob) {
     for (i = 0 ; i < si ; i++) {
       free(sarr[i].tat->probs);
       free(sarr[i].tat);
     }
   }

   if (start != NULL) {
     if(opt_prob) {
       rst->escore = start->tatprob;
     } else {
       rst->escore = 2.0;
     }

     rst->segnum = rst->seglen = 0;
     for(sj = start ; sj != NULL; sj = sj->prev) {
       rst->segnum++;
       rst->seglen += sj->vp->stop - sj->vp->start + 1;
     }
     return (start->score);
   } else {
     rst->escore = 1.0;
   }

   rst->segnum = rst->seglen = 0;
   return (0);
}

void
kssort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->score >= v[j + gap]->score)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

void
kpsort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->start <= v[j + gap]->start)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

/* calculate the 100% identical score */
int
shscore(const unsigned char *aa0, const int n0, int **pam2, int nsq)
{
  int i, sum;
  for (i=0,sum=0; i<n0; i++)
    if (aa0[i] != EOSEQ && aa0[i]<=nsq) sum += pam2[aa0[i]][aa0[i]];
  return sum;
}

/* sorts alignments from right to left (back to front) based on stop */

void
krsort (v, n)
struct savestr *v[];
int     n;
{
   int     gap, i, j;
   struct savestr *tmp;

   for (gap = n / 2; gap > 0; gap /= 2)
      for (i = gap; i < n; i++)
	 for (j = i - gap; j >= 0; j -= gap)
	 {
	    if (v[j]->stop > v[j + gap]->stop)
	       break;
	    tmp = v[j];
	    v[j] = v[j + gap];
	    v[j + gap] = tmp;
	 }
}

void
do_walign (const unsigned char *aa0, int n0,
	   const unsigned char *aa1, int n1,
	   int frame,
	   struct pstruct *ppst, 
	   struct f_struct *f_str, 
	   struct a_res_str *a_res,
	   int *have_ares)
{
  int hoff, n10;
  struct rstruct rst;
  int ib, i;
  unsigned char *aa0t;
  const unsigned char *aa1p;
  struct savestr *vmptr;

#ifdef TFAST
  f_str->n10 = n10 = aatran(aa1,f_str->aa1x,n1,frame);
  aa1p = f_str->aa1x;
#else
  n10 = n1;
  aa1p = aa1;
#endif

  do_fasts(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff, 1, f_str->maxsav_w);

  /* the alignment portion takes advantage of the information left
     over in f_str after do_fasts is done.  in particular, it is
     easy to run a modified sconn() to produce the alignments.

     unfortunately, the alignment display routine wants to have
     things encoded as with bd_align and sw_align, so we need to do that.
     */

  /* unnecessary; do_fasts just did this */
  /*  kssort(f_str->vptr,f_str->nsave);  */

   /* at some point, we want one best score for each of the segments */

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

   if ((aa0t = (unsigned char *)calloc(n0+1,sizeof(unsigned char)))==NULL) {
     fprintf(stderr," cannot allocate aa0t %d\n",n0+1);
     exit(1);
   }

   /* copy aa0[] into f_str->aa0t[] */
   for (i=0; i<n0; i++) f_str->aa0t[i] = aa0t[i] = aa0[i];
   f_str->aa0t[i] = aa0t[i] = '\0';

   a_res->nres = sconn_a (aa0t,n0,aa1p,n10,f_str, a_res, ppst);

   free(aa0t);

   a_res->res = f_str->res;
   *have_ares = 0;
   a_res->sw_score = rst.score[0];
}

/* this version of sconn is modified to provide alignment information */
/* in addition, it needs to know whether a segment has been used before */

/* sconn_a fills in the res[nres] array, but this is passed implicitly
   through f_str->res[f_str->nres] */

int sconn_a (unsigned char *aa0, int n0,
	     const unsigned char *aa1, int n1,
	     struct f_struct *f_str,
	     struct a_res_str *a_res,
	     struct pstruct *ppst)
{
   int     i, si, cmpp (), n;
   unsigned char *aa0p;
   int sx, dx, doff, *aa0tip;

   struct savestr **v;
   struct slink *start, *sl, *sj, *so, *sarr;
   int     lstart, lstop, ltmp, plstart, tstart, plstop, ptstop, ptstart, tstop;

   int *res, nres, tres;

   double tatprob;

/*	sort the score left to right in lib pos */

   v = f_str->vptr;
   n = f_str->nsave;
   sarr = f_str->sarr;

   /* set things up in case nothing fits */
   if (n <=0 || v[0]->score <= 0) return 0;

   if (v[0]->score < 0) {
     sarr[0].vp = v[0];
     sarr[0].score = v[0]->score;
     sarr[0].next = NULL;
     sarr[0].prev = NULL;
     start = &sarr[0];
   }
   else {

     krsort (v, n);	/* sort from left to right in library */

     start = NULL;

     /*	for each alignment, see if it fits */


     for (i = 0, si = 0; i < n; i++) {
       /*	if the score is less than the join threshold, skip it */

       if (v[i]->score < 0) continue;

       lstart = v[i]->start;
       lstop = v[i]->stop;
       tstart = lstart - v[i]->dp + f_str->noff;
       tstop = lstop - v[i]->dp + f_str->noff;

       /*	put the alignment in the group */

       sarr[si].vp = v[i];
       sarr[si].score = v[i]->score;
       sarr[si].next = NULL;
       sarr[si].prev = NULL;
       sarr[si].tat = NULL;

       sarr[si].tatprob = 
	 calc_tatusov(NULL, &sarr[si], aa0, n0, aa1, n1, 
		      ppst->pam2[0], ppst->nsq, f_str, 
		      ppst->pseudocts, 1, ppst->zsflag);
       sarr[si].tat = sarr[si].newtat;


       /* 	if it fits, then increase the score */
       /* start points to a sorted (by total score) list of candidate
	  overlaps */

       for (sl = start; sl != NULL; sl = sl->next) { 
	 plstart = sl->vp->start;
	 plstop = sl->vp->stop;
	 ptstart = plstart - sl->vp->dp + f_str->noff;
	 ptstop = plstop - sl->vp->dp + f_str->noff;
#ifndef FASTM
	 if (plstart > lstop && (ptstop < tstart || ptstart > tstop)) {
#else
         if (plstop > lstart && ptstart > tstop) {
#endif
	   /* alignment always uses probabilistic scoring ... */
	   /*   sarr[si].score = sl->score + v[i]->score;
		sarr[si].prev = sl;
		break; */		/* quit as soon as the alignment has been added */

	   tatprob = calc_tatusov(sl, &sarr[si], aa0, n0, aa1, n1, 
				  ppst->pam2[0], ppst->nsq, f_str, 
				  ppst->pseudocts, 1, ppst->zsflag);
	   /* if our tatprob gets worse when we add this, forget it */
	   if(tatprob > sarr[si].tatprob) {
	     free(sarr[si].newtat->probs); /* get rid of new tat struct */
	     free(sarr[si].newtat);
	     continue; /* reuse this sarr[si] */
	   } else {
	     sarr[si].tatprob = tatprob;
	     free(sarr[si].tat->probs); /* get rid of old tat struct */
	     free(sarr[si].tat);
	     sarr[si].tat = sarr[si].newtat;
	     sarr[si].prev = sl;
	     sarr[si].score = sl->score + v[i]->score;
	     /*
	       fprintf(stderr,"sconn %d added %d/%d getting %d; si: %d, tat: %g\n",
	       i,v[i]->start, v[i]->score,sarr[si].score,si, tatprob);
	     */
	     break;
	   }
	 }
       }

       /* now recalculate the list of best scores */
       if (start == NULL)
	 start = &sarr[si];	/* put the first one in the list */
       else
	 for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	   /* if (sarr[si].score > sj->score) { */ /* new score better than old */
	   if ( sarr[si].tatprob < sj->tatprob ||
		((sarr[si].tatprob == sj->tatprob) && sarr[si].score > sj->score) ) {
	     sarr[si].next = sj;		/* next best after new score */
	     if (so != NULL)
	       so->next = &sarr[si];	/* prev_best->next points to best */
	     else  start = &sarr[si];	/* start points to best */
	     break;			/* stop looking */
	   }
	   so = sj;		/* previous candidate best */
	 }
       si++;				/* increment to next alignment */
     }
   }

   for (i = 0 ; i < si ; i++) {
     free(sarr[i].tat->probs);
     free(sarr[i].tat);
   }

   res = f_str->res;
   tres = nres = 0;
   aa0p = aa0;
   aa0tip = f_str->aa0ti;	/* point to temporary index */
   a_res->min1 = start->vp->start;
   a_res->min0 = 0;

   for (sj = start; sj != NULL; sj = sj->prev ) {
     doff = (int)(aa0p-aa0) - (sj->vp->start-sj->vp->dp+f_str->noff);
     
     /* fprintf(stderr,"doff: %3d\n",doff); */
     
     for (dx=sj->vp->start,sx=sj->vp->start-sj->vp->dp+f_str->noff;
	  dx <= sj->vp->stop; dx++) {
       *aa0tip++ = f_str->aa0i[sx];	/* save index */
       *aa0p++ = f_str->aa0t[sx++];	/* save sequence at index */
       tres++;
       res[nres++] = 0;
     }
     sj->vp->dp -= doff;
     if (sj->prev != NULL) {
       if (sj->prev->vp->start - sj->vp->stop - 1 > 0 )
	 tres += res[nres++] = (sj->prev->vp->start - sj->vp->stop - 1);
     }

     /*
     fprintf(stderr,"t0: %3d, tx: %3d, l0: %3d, lx: %3d, dp: %3d noff: %3d, score: %3d\n",
       sj->vp->start - sj->vp->dp + f_str->noff,
       sj->vp->stop - sj->vp->dp + f_str->noff,
       sj->vp->start,sj->vp->stop,sj->vp->dp,
       f_str->noff,sj->vp->score);

       fprintf(stderr,"%3d - %3d: %3d\n",
       sj->vp->start,sj->vp->stop,sj->vp->score);
     */
     a_res->max1 = sj->vp->stop+1;
     a_res->max0 = a_res->max1 - sj->vp->dp + f_str->noff;
   }

   /*
   fprintf(stderr,"(%3d - %3d):(%3d - %3d)\n",
	   a_res->min0,a_res->max0,a_res->min1,a_res->max1);
   */
   
   /* now replace f_str->aa0t with aa0
      (f_str->aa0t is permanent, aa0 is not)*/
   for (i=0; i<n0; i++) f_str->aa0t[i] = aa0[i];

   return tres;
}

/* for fasts (and fastf), pre_cons needs to set up f_str as well as do
   necessary translations - for right now, simply do do_walign */

void
pre_cons(const unsigned char *aa1, int n1, int frame, struct f_struct *f_str) {

#ifdef TFAST
  f_str->n10=aatran(aa1,f_str->aa1x,n1,frame);
#endif

}

/* aln_func_vals - set up aln.qlfact, qlrev, llfact, llmult, frame, llrev */
/* call from calcons, calc_id, calc_code */
void 
aln_func_vals(int frame, struct a_struct *aln) {

#ifdef TFAST
  aln->qlrev = 0;
  aln->qlfact= 1;
  aln->llfact = aln->llmult = 3;
  if (frame > 3) aln->llrev = 1;
  else aln->llrev = 0;
  aln->frame = 0;
#else	/* FASTS */
  aln->llfact = aln->llmult = aln->qlfact = 1;
  aln->llrev = aln->qlrev = 0;
  aln->frame = 0;
#endif
}

void aaptrshuffle(unsigned char *res, int n) {

  int i, j;
  unsigned char tmp;

  for( i = n; --i; ) {

    /* j = nrand(i); if (i == j) continue; */ /* shuffle */
    j = (n - 1) - i; if (i <= j ) break; /* reverse */

    tmp = res[i];
    res[i] = res[j];
    res[j] = tmp;
  }
}

void aa0shuffle(unsigned char *aa0, int n0, struct f_struct *f_str) {

  int i;
  int j;

  for(i = 0 ; i < f_str->nm0 ; i++) { /* for each fragment */

    aaptrshuffle(&(aa0[f_str->nmoff[i]]), 
		 f_str->nmoff[i+1] - f_str->nmoff[i] - 1 );

  }

}

#ifdef PCOMPLIB

#include "structs.h"
#include "p_mw.h"

void
update_params(struct qmng_str *qm_msg,
	      struct mngmsg *m_msg, struct pstruct *ppst)
{
  m_msg->n0 = ppst->n0 = qm_msg->n0;
  m_msg->nm0 = qm_msg->nm0;
  m_msg->escore_flg = qm_msg->escore_flg;
  m_msg->qshuffle = qm_msg->qshuffle;
}
#endif