dropnnw.c

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

     (do_karlin(aa1, n1, ppst->pam2[0], ppst,f_str->aa0_f, 
		f_str->kar_p, &lambda, &H)>0)) {
    rst->comp = 1.0/lambda;
    rst->H = H;
  }
  else {rst->comp = rst->H = -1.0;}

}

#define gap(k)  ((k) <= 0 ? 0 : g+h*(k))	/* k-symbol indel cost */

static int
FGLOBAL_ALIGN(int *waa, const unsigned char *aa1,
	      int n0, int n1,
	      int q, int r,
	      struct swstr *ss)
{
  const unsigned char *aa1p;
  register int *pwaa;
  int i,j;
  struct swstr *ssj;
  int t;
  int e, f, h, p;
  int qr;
  int score;
  int ij, max_col, max_ij;

  /* q - gap open is positve */
  /* r - gap extend is positive */

  qr = q+r;

  score = -BIGNUM;

  /* initialize 0th row */
  ss[0].H = 0;
  ss[0].E = t = -q;
  for (ssj = ss+1; ssj <= ss+n0 ; ssj++) {
    ssj->H = t = t - r;
    ssj->E = t - q;
  }

  aa1p = aa1;
  t = -q;
  while (*aa1p) {
    p = ss[0].H;
#if defined(GLOBAL_GLOBAL)
    ss[0].H = h = t = t - r;
#else
    ss[0].H = h = t = 0;
#endif
    f = t - q;
    pwaa = waa + (*aa1p++ * n0);
    for (ssj = ss+1; ssj <= ss+n0; ssj++) {	/* go across query */
      if ((h =   h    - qr) > (f =   f    - r)) f = h;
      if ((h = ssj->H - qr) > (e = ssj->E - r)) e = h;
      h = p + *pwaa++;
      if (h < f) h = f;
      if (h < e) h = e;
      p = ssj->H;
      ssj->H = h;
      ssj->E = e;
    }
#if !defined(GLOBAL_GLOBAL)
    if (h > score) {
      score = h;	/* at end of query, update score */
    }
#endif
  }				/* done with forward pass */
#ifdef GLOBAL_GLOBAL
  score = h;
#endif
  return score;
}

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)
{
}

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)
{
  const unsigned char *aa0p, *aa1p;
  register int *pwaa;
  register int i, j;
  register struct swstr *ssj;
  struct swstr *ss;
  int *res, *waa, **pam2p;
  int e, f, h, p;
  int q, r, qr, t;
  int     score;
  int cost, I, J, K, L;

  /* now we need alignment storage - get it */
  if ((a_res->res = (int *)calloc((size_t)f_str->max_res,sizeof(int)))==NULL) {
    fprintf(stderr," *** cannot allocate alignment results array %d\n",f_str->max_res);
    exit(1);
  }
  
  *have_ares = 0x3;	/* set 0x2 bit to indicate local copy */
  a_res->next = NULL;

  ss = f_str->ss;	/* row vector */
  waa = f_str->waa_a;	/* this time use universal pam2[0] */
  pam2p = f_str->pam2p[0];	/* pam2p profile */

#ifdef OLD_FASTA_GAP
  q = -(ppst->gdelval - ppst->ggapval);
#else
  q = -ppst->gdelval;
#endif

  r = -ppst->ggapval;
  qr = q + r;

  score = -BIGNUM;
  J = n0-1; L = 0;	/* alignments are global in aa0[n0] */

  /* initialize 0th row */
  ss[0].H = 0;
  ss[0].E = t = -q;	/* must be re-initialized because it was
			   filled in the reverse direction in previous
			   invocations */
  for (ssj=ss+1; ssj <= ss+n0 ; ssj++) {
    ssj->H = t = t - r;
    ssj->E = t - q;
  }

  aa1p = aa1;
  i = 0;
  t = -q;
  while (*aa1p) {
    p = ss[0].H;
#ifndef GLOBAL_GLOBAL
    ss[0].H = h = t = 0;
#else
    ss[0].H = h = t = t - r;
#endif
    f = t - q;
    pwaa = waa + (*aa1p++ * n0);
    for (ssj = ss+1; ssj <= ss+n0; ssj++) {
      if ((h =   h     - qr) > /* gap open from left best */
	  /* gap extend from left gapped */
	  (f =   f     - r)) f = h;	/* if better, use new gap opened */
      if ((h = ssj->H - qr) >	/* gap open from up best */
	  /* gap extend from up gap */
	  (e = ssj->E - r)) e = h;	/* if better, use new gap opened */
      h = p + *pwaa++;		/* diagonal match */
      if (h < f ) h = f;	/* left gap better, reset */
      if (h < e ) h = e;	/* up gap better, reset */
      p = ssj->H;		/* save previous best score */
      ssj->H = h;		/* save (new) up diag-matched */
      ssj->E = e;		/* save upper gap opened */
    }
#ifndef GLOBAL_GLOBAL
    if (h > score) {		/* ? new best score at the end of each row */
      score = h;		/* save best */
      I = i;			/* row */
    }
#endif
    /*
    fprintf(stderr," r %d - score: %d ssj[]: %d\n", i,score,
	    ss[(i <= n0) ? i-1 : n0-1].H);
    */
    i++;
  }				/* done with forward pass */

#ifdef GLOBAL_GLOBAL
  cost = score = h;
  K = 0;
  I = n1 - 1;
#else
  /*   fprintf(stderr, " r: %d - score: %d\n", I, score); */

  /* to get the start point, go backwards */
  
  cost = -BIGNUM;
  K = 0;
  ss[n0].H = 0;
  t = -q;
  for (ssj=ss+n0-1; ssj>=ss; ssj--) {
    ssj->H = t = t - r;
    ssj->E= t - q;
  }

  t = 0;
  for (i=I; i>=0; i--) {
    p = ss[n0].H;
    ss[n0].H = h = t = t-r;
    f = t-q;
    for (ssj=ss+J, j= J; j>=0; ssj--, j--) {
      if ((h =   h     - qr) > /* gap open from left best */
	  /* gap extend from left gapped */
	  (f =   f     - r)) f = h;	/* if better, use new gap opened */
      if ((h = ssj->H - qr) >	/* gap open from up best */
	  /* gap extend from up gap */
	  (e = ssj->E - r)) e = h;	/* if better, use new gap opened */
      h = p + pam2p[j][aa1[i]];		/* diagonal match */
      if (h < f ) h = f;	/* left gap better, reset */
      if (h < e ) h = e;	/* up gap better, reset */
      p = ssj->H;		/* save previous best score */
      ssj->H = h;		/* save (new) up diag-matched */
      ssj->E = e;		/* save upper gap opened */

      /*
      f = max (f,h-q)-r;
      ssj->E = max(ssj->E,ssj->H-q)-r;
      h = max(max(ssj->E,f),p+f_str->pam2p[0][j][aa1[i]]);
      p = ssj->H;
      ssj->H=h;
      */
    }
    if (h > cost) {
      cost = h;
      K = i;
      if (cost >= score) goto found;
    }
  }
  /* at this point, ss[0].E has a very high value for good alignments */
 found:
#endif /* not GLOBAL_GLOBAL */

/*   fprintf(stderr," *** %d=%d: L: %3d-%3d/%3d; K: %3d-%3d/%3d\n",score,cost,L,J+1,n0,K,I+1,n1); */

/* in the f_str version, the *res array is already allocated at 4*n0/3 */

  a_res->max0 = J+1; a_res->min0 = L; a_res->max1 = I+1; a_res->min1 = K;
  
/* this code no longer refers to aa0[], it uses pam2p[0][L] instead */
  NW_ALIGN(L,&aa1[K-1],J-L+1,I-K+1,f_str->pam2p[0],q,r,
	   a_res->res,&a_res->nres,f_str);

/*  DISPLAY(&aa0[L-1],&aa1[K-1],J-L+1,I-K+1,res,L,K,ppst->sq); */

/* return *res and nres */

  a_res->sw_score = score;
}


/*
#define XTERNAL
#include "upam.h"

void
print_seq_prof(unsigned char *A, int M,
	       unsigned char *B, int N,
	       int **w, int iw, int dir) {
  char c_max;
  int i_max, j_max, i,j;

  char *c_dir="LRlr";

  for (i=1; i<=min(60,M); i++) {
    fprintf(stderr,"%c",aa[A[i]]);
  }
  fprintf(stderr, - %d\n,M);

  for (i=0; i<min(60,M); i++) {
    i_max = -1;
    for (j=1; j<21; j++) {
      if (w[iw+i][j]> i_max) {
	i_max = w[iw+i][j]; 
	j_max = j;
      }
    }
    fprintf(stderr,"%c",aa[j_max]);
  }
  fputc(':',stderr);

  for (i=1; i<=min(60,N); i++) {
    fprintf(stderr,"%c",aa[B[i]]);
  }

  fprintf(stderr," -%c: %d,%d\n",c_dir[dir],M,N);
}
*/

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

  aln->llfact = aln->llmult = aln->qlfact = 1;
  aln->qlrev = aln->llrev = 0;
  aln->frame = 0;
}

#ifdef PCOMPLIB
#include "p_mw.h"
void
update_params(struct qmng_str *qm_msg, struct pstruct *ppst)
{
  ppst->n0 = qm_msg->n0;
}
#endif