dropff2.c

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

   /* sarr[] saves each alignment score/position, and provides a link
      back to the previous alignment that maximizes the score */

   /*	sort the score left to right in lib pos */
   kpsort (v, n);

   start = NULL;

   /* for the remaining runs, see if they fit */
   for (i = 0, si = 0; i < n; i++) {

     /* if the score is less than the gap penalty, it never helps */
     if (!opt_prob && (v[i]->score < cgap) ){ continue; }

     lstart = v[i]->start;

     /* 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;
     
     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);
       sarr[si].tat = sarr[si].newtat;
     }

     /* if it fits, then increase the score */
     for (sl = start; sl != NULL; sl = sl->next) {
       plstop = sl->vp->stop;
       /* if end < start or start > end, add score */
       if (plstop < lstart ) {
	 if(!opt_prob) {
	   sarr[si].score = sl->score + v[i]->score;
	   sarr[si].prev = sl;
	   /*
	     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, 2.0);
	   */
	   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;
	   } 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 TAT %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 - resort the scores */
     if (start == NULL) {
       start = &sarr[si];
     } else {
       if(!opt_prob) { /* sort by scores */
	 for (sj = start, so = NULL; sj != NULL; sj = sj->next) {
	   if (sarr[si].score > sj->score) { /* if new score > best score */
	     sarr[si].next = sj;	     /* previous best linked to best */
	     if (so != NULL)		
	       so->next = &sarr[si];	     /* old best points to new best */
	     else
	       start = &sarr[si];
	     break;
	   }
	   so = sj;			     /* old-best saved in so */
	 }
       } else { /* sort by tatprobs */
	 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 {

     if(opt_prob)
       rst->escore = 1.0;
     else
       rst->escore = 2.0;

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

void
kssort (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;
	 }
}

/* 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;
  unsigned char *aa0t;
  const unsigned char *aa1p;

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

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

  /* the alignment portion takes advantage of the information left
     over in f_str after do_fastf 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.
     */

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

   kssort (f_str->vptr, f_str->nsave);
   f_str->aa0ix = 0;
   if (f_str->nsave > f_str->nm0) f_str->nsave = f_str->nm0;
   for (ib=0; ib < f_str->nm0; ib++) {
     if (f_str->vptr[ib]->score > 0) {
       f_str->vptr[ib]->score = 
	 ma_spam (f_str->aa0, n0, aa1p, f_str->vptr[ib], ppst, f_str);
     }
   }

   /* after ma_spam is over, we need to reset aa0 */
   for (ib = 0; ib < n0; ib++) {
     if (f_str->aa0[ib] >= 32) f_str->aa0[ib] -= 32;
   }

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

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

  a_res->nres = sconn_a (aa0t,n0, ppst->param_u.fa.cgap, f_str,a_res);
  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 */

int sconn_a (unsigned char *aa0, int n0, int cgap, 
	     struct f_struct *f_str,
	     struct a_res_str *a_res)
{
   int     i, si, cmpp (), n;
   unsigned char *aa0p;
   int sx, dx, doff;

   struct savestr **v;
   struct slink {
     int     score;
     struct savestr *vp;
     struct slink *snext;
     struct slink *aprev;
   } *start, *sl, *sj, *so, sarr[MAXSAV];
   int     lstop, plstart;
   int *res, nres, tres;

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

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

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

     lstop = v[i]->stop;		/* have right-most lstart */

/*	put the alignment in the group */

     sarr[si].vp = v[i];
     sarr[si].score = v[i]->score;
     sarr[si].snext = NULL;
     sarr[si].aprev = NULL;

/* 	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->snext) { 
       plstart = sl->vp->start;
       if (plstart > lstop ) {
	 sarr[si].score = sl->score + v[i]->score;
	 sarr[si].aprev = sl;
	 break;		/* quit as soon as the alignment has been added */
       }
     }

/* 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->snext) {
	 if (sarr[si].score > sj->score) { /* new score better than old */
	   sarr[si].snext = sj;		/* snext best after new score */
	   if (so != NULL)
	     so->snext = &sarr[si];	/* prev_best->snext points to best */
	   else  start = &sarr[si];	/* start points to best */
	   break;			/* stop looking */
	 }
	 so = sj;		/* previous candidate best */
       }
     si++;				/* increment to snext alignment */
   }

   /* we have the best set of alignments, write them to *res */
   if (start != NULL) {
     res = f_str->res;	/* set a destination for the alignment ops */
     tres = nres = 0;	/* alignment op length = 0 */
     aa0p = aa0;	/* point into query (needed for calcons later) */
     a_res->min1 = start->vp->start;	/* start in library */
     a_res->min0 = 0;			/* start in query */
     for (sj = start; sj != NULL; sj = sj->aprev ) {
       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++) {
	 *aa0p++ = f_str->aa0t[sx++];	/* copy residue into aa0 */
	 tres++;			/* bump alignment counter */
	 res[nres++] = 0;		/* put 0-op in res */
       }
       sj->vp->dp -= doff;
       if (sj->aprev != NULL) {
	 if (sj->aprev->vp->start - sj->vp->stop - 1 > 0 )
	 /* put an insert op into res to get to next aligned block */
	   tres += res[nres++] = (sj->aprev->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;
       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 */
     for (i=0; i<n0; i++) f_str->aa0t[i] = aa0[i];

     return tres;
   }
   else return (0);
}

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

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;
  aln->frame = 0;
  if (frame > 3) aln->llrev = 1;
#else	/* FASTF */
  aln->llfact = aln->qlfact = aln->llmult = 1;
  aln->llrev = aln->qlrev = 0;
  aln->frame = 0;
#endif
}

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

  int i, j, k;
  unsigned char tmp;

  for (i = f_str->nmoff-1 ; --i ; ) {

    /* j = nrand(i); if (i == j) continue;*/       /* shuffle columns */ 
    j = (f_str->nmoff - 2) - i; if (i <= j) break; /* reverse columns */

    /* swap all i'th column residues for all j'th column residues */
    for(k = 0 ; k < f_str->nm0 ; k++) {
      tmp = aa0[(k * (f_str->nmoff)) + i];
      aa0[(k * (f_str->nmoff)) + i] = aa0[(k * (f_str->nmoff)) + j];
      aa0[(k * (f_str->nmoff)) + j] = tmp;
    }
  }
}

#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