lsim4.c

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

	  pj = l_ptr->c_ss[l_ptr->n1].JJ;
	  ci = fi = l_ptr->n1;
	  limit = l_ptr->n1 - 1;
	}

	for ( i = limit; i >= l_ptr->m1 ; i-- ) {
	  f = f - R;
	  c = c - qr;
	  ORDER(f, fi, fj, c, ci, cj)
	  c = l_ptr->c_ss[i].HH - qr; 
	  ci = l_ptr->c_ss[i].II;
	  cj = l_ptr->c_ss[i].JJ;
	  d = l_ptr->c_ss[i].WW - R;
	  di = l_ptr->c_ss[i].XX;
	  dj = l_ptr->c_ss[i].YY;
	  ORDER(d, di, dj, c, ci, cj)
	  c = 0;
	  DIAG(i, l_ptr->n1, c, p+va[A[i]])
	  if ( c <= 0 ) { c = 0; ci = i; cj = l_ptr->n1; }
	  else { ci = pi; cj = pj; }
	  ORDER1(c, ci, cj, d, di, dj)
	  ORDER1(c, ci, cj, f, fi, fj)
	  p = l_ptr->c_ss[i].HH;
	  l_ptr->c_ss[i].HH = c;
	  pi = l_ptr->c_ss[i].II;
	  pj = l_ptr->c_ss[i].JJ;
	  l_ptr->c_ss[i].II = ci;
	  l_ptr->c_ss[i].JJ = cj;
	  l_ptr->c_ss[i].WW = d;
	  l_ptr->c_ss[i].XX = di;
	  l_ptr->c_ss[i].YY = dj;
	  if ( c > mini_score ) *flag_p = 1;
	  if ( ! cflag && ( (ci > l_ptr->rl && cj > l_ptr->cl) || (di > l_ptr->rl && dj > l_ptr->cl)
			    || (fi > l_ptr->rl && fj > l_ptr->cl) ) )  cflag = 1;
	}

	l_ptr->CCC[l_ptr->n1].CC = l_ptr->c_ss[l_ptr->m1].HH;
	l_ptr->CCC[l_ptr->n1].RR = l_ptr->c_ss[l_ptr->m1].II;
	l_ptr->CCC[l_ptr->n1].EE = l_ptr->c_ss[l_ptr->m1].JJ;
	l_ptr->CCC[l_ptr->n1].DD = f;
	l_ptr->CCC[l_ptr->n1].SS = fi;
	l_ptr->CCC[l_ptr->n1].FF = fj;
	if ( ! rflag && ( (ci > l_ptr->rl && cj > l_ptr->cl) || (di > l_ptr->rl && dj > l_ptr->cl)
			  || (fi > l_ptr->rl && fj > l_ptr->cl )) ) rflag = 1;
      }
    }
    if (( l_ptr->m1 == 1 && l_ptr->n1 == 1) || no_cross(flag_p, v_ptr->LIST, l_ptr) ) break;
  }
  l_ptr->m1--;
  l_ptr->n1--;
}

/* recompute the area on forward pass */
static void
small_pass(const unsigned char *A,
	   const unsigned char *B,
	   int mini_score,
	   int **pam2, int Q, int R,
	   int nseq,
	   struct vert_str *v_ptr,
	   struct l_struct *l_ptr) {

  int  i, j;		/* row and column indices */
  int  c;		/* best score at current point */
  int  f;		/* best score ending with insertion */
  int  d;		/* best score ending with deletion */
  int  p;		/* best score at (i-1, j-1) */
  int  ci, cj;		/* end-point associated with c */ 
  int  fi, fj;		/* end-point associated with f */
  int  pi, pj;		/* end-point associated with p */
  space_ptr sp;
  pair_ptr z;
  int q, r, qr;
  int  *va;		/* pointer to pam2(A[i], B[j]) */
  
  int  limit;		/* lower bound on j */

  q = Q; r = R; qr = q + r;

  for ( sp = &l_ptr->CCC[l_ptr->n1 + 1], j = l_ptr->n1+1; sp <= &l_ptr->CCC[l_ptr->nn] ; sp++, j++ ) {
    sp->CC = 0;
    sp->RR = l_ptr->m1;
    sp->EE = j;
    sp->DD = - (qr);
    sp->SS = l_ptr->m1+1;
    sp->FF = j;
  }

  for ( i = l_ptr->m1 + 1; i <= l_ptr->mm; i++) {
    c = 0;				/* Initialize column 0 */
    f = - (qr);
    ci = fi = i;
    va = pam2[A[i]];
    if ( nseq == 2 || i <= l_ptr->n1 ) {
      p = 0;
      pi = i - 1;
      cj = fj = pj = l_ptr->n1;
      limit = l_ptr->n1 + 1;
    }
    else {
      p = l_ptr->CCC[i].CC;
      pi = l_ptr->CCC[i].RR;
      pj = l_ptr->CCC[i].EE;
      cj = fj = i;
      limit = i + 1;
    }

    for ( j = limit, sp = &l_ptr->CCC[j] ; j <= l_ptr->nn ; j++, sp++ )  {  
      d = sp->DD;
      c = -1;
      DIAG(i, j, c, p+va[B[j]])		/* diagonal */
      if (c < 0) {
	p = sp->CC; pi = sp->RR; pj = sp->EE;
	if (f >= 0) {
	  c = f; ci = fi; cj = fj;
	  ORDER1(c, ci, cj,  d, sp->SS, sp->FF)
	  sp->CC = c; sp->RR = ci; sp->EE = cj;
	  sp->DD -= r; f-=r;
	} 
	else if (d >= 0) {
	  sp->CC = d; sp->RR = sp->SS; sp->EE = sp->FF;  
	  sp->DD -= r;
	}
	else {
	  sp->CC = 0;
	  sp->RR=i;
	  sp->EE = j;
	}
      }
      else { 
	ci = pi; cj = pj;
	ORDER1(c, ci, cj,  f, fi, fj)
	  ORDER1(c, ci, cj,  d, sp->SS, sp->FF)
	  p = sp->CC;
	sp->CC = c;
	pi = sp->RR;
	sp->RR = ci;
	pj = sp->EE;
	sp->EE = cj;
	f-=r;
	if (c >= qr) {
	  if ( c > mini_score ) /* add the score into list */
	    addnode(c, ci, cj, i, j, v_ptr);
	  d -= r; c-=qr;
	  ORDER1(f, fi, fj, c, ci, cj)
          ORDER1(d, sp->SS, sp->FF, c, ci, cj)
	  sp->DD = d;
	}
	else {
	  sp->DD -= r;
	}
      }
    }
  }
}

/* Add a new node into list.  */

static void
addnode(int c, int ci, int cj, int i, int j, struct vert_str *v_ptr) {

  bool found;	/* 1 if the node is in LIST */

  found = 0;
  if ( v_ptr->most != NULL && v_ptr->most->STARI == ci && v_ptr->most->STARJ == cj)
    found = 1;
  else {
    for ( v_ptr->most = v_ptr->LIST; v_ptr->most; v_ptr->most = v_ptr->most->next ) {
      if ( v_ptr->most->STARI == ci && v_ptr->most->STARJ == cj) {
	found = 1;
	break;
      }
    }
  }
  if ( found ) {
    if ( v_ptr->most->SCORE < c ) {
      v_ptr->most->SCORE = c;
      v_ptr->most->ENDI = i;
      v_ptr->most->ENDJ = j;
    }
    if ( v_ptr->most->TOP > i ) v_ptr->most->TOP = i;
    if ( v_ptr->most->BOT < i ) v_ptr->most->BOT = i;
    if ( v_ptr->most->LEFT > j ) v_ptr->most->LEFT = j;
    if ( v_ptr->most->RIGHT < j ) v_ptr->most->RIGHT = j;
  }
  else { 
    v_ptr->numnode++;
    v_ptr->most = (vertex_p) ckalloc(1,sizeof(vertex));
    v_ptr->most->SCORE = c;
    v_ptr->most->STARI = ci;
    v_ptr->most->STARJ = cj;
    v_ptr->most->ENDI = i;
    v_ptr->most->ENDJ = j;
    v_ptr->most->TOP = v_ptr->most->BOT = i;
    v_ptr->most->LEFT = v_ptr->most->RIGHT = j;
    v_ptr->most->next = v_ptr->LIST;
    v_ptr->LIST = v_ptr->most;
  }
}

/* Find and remove the largest score in list */

static vertex_p
findmax(struct vert_str *v_ptr) {
  vertex_p  ap, cur;
  register int score;

  for ( score = (v_ptr->LIST)->SCORE, cur = NULL, ap = (v_ptr->LIST); ap->next; ap = ap->next) {
    if ( ap->next->SCORE > score ) {
      cur = ap; score = ap->next->SCORE;
    }
  }
  if (cur) {ap = cur->next; cur->next = ap->next; }
  else { ap = v_ptr->LIST; v_ptr->LIST = (v_ptr->LIST)->next;}
  v_ptr->numnode--;
  v_ptr->most = v_ptr->LIST;
  return ( ap );
}

/* return 1 if no node in LIST share vertices with the area */

static bool
no_cross(int *flag_p, vertex_p LIST, struct l_struct *l_ptr) {

  vertex_p  cur;

  for ( cur = LIST; cur; cur = cur->next ) { 
    if ( cur->STARI <= l_ptr->mm && cur->STARJ <= l_ptr->nn && cur->BOT >= l_ptr->m1-1 && 
	 cur->RIGHT >= l_ptr->n1-1 && (cur->STARI < l_ptr->rl || cur->STARJ < l_ptr->cl)) { 
      if ( cur->STARI < l_ptr->rl ) l_ptr->rl = cur->STARI;
      if ( cur->STARJ < l_ptr->cl ) l_ptr->cl = cur->STARJ;
      *flag_p = 1;
      break;
    }
  }
  return !cur;
}

/* The following definitions are for function diff() */

#define gap(k)  ((k) <= 0 ? 0 : Q+R*(k)) /* k-symbol indel score */

/* Append "Delete k" op */
#define DEL(k)				\
{ l_ptr->I += k;			\
  if (*last < 0)			\
    *last = (*sapp)[-1] -= (k);		\
  else {				\
    *last = (*sapp)[0] = -(k);		\
    (*sapp)++;				\
  }					\
}

/* Append "Insert k" op */
#define INS(k)				\
{ l_ptr->J += k;			\
  if (*last < 0) {			\
   (*sapp)[-1] = (k);			\
   (*sapp)[0] = *last;			\
   (*sapp)++;				\
  }					\
  else {				\
    *last = (*sapp)[0] = (k);		\
    (*sapp)++;				\
  }					\
}

/* diff(A,B,M,N,tb,te) returns the score of an optimum conversion between
   A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
   and appends such a conversion to the current script.                   */

static int
diff(const unsigned char *A,
     const unsigned char *B,
     int M, int N,
     int tb, int te,
     int **pam2, int Q, int R,
     int **sapp, int *last,
     struct l_struct *l_ptr) {

  int   midi, midj, type;	/* Midpoint, type, and cost */
  int midc;

  register int i, j;
  register int c, e, d, s;

  pair_ptr z;

  int t;
  int *va;
  int qr;

  bool tt;

  qr = Q + R;

  /* Boundary cases: M <= 1 or N == 0 */

  if (N <= 0){
    if (M > 0) {DEL(M)}
    return - gap(M);
  }

  if (M <= 1) {
    if (M <= 0) {
      INS(N)
      return - gap(N);
    }

    if (tb > te) tb = te;
    midc = - (tb + R + gap(N) );
    midj = 0;

    va = pam2[A[1]];
    for (j = 1; j <= N; j++) {
      for ( tt = 1, z = l_ptr->row[l_ptr->I+1]; z != 0; z = z->NEXT ) {
	if ( z->COL == j+l_ptr->J ) { tt = 0; break; }		
      }
      if (tt) {
	c = va[B[j]] - ( gap(j-1) + gap(N-j) );
	if (c > midc) { midc = c;  midj = j; }
      }
    }

    if (midj == 0) { INS(N) DEL(1) }
    else  {
      if (midj > 1) INS(midj-1)
      *last = (*sapp)[0] = 0;
      (*sapp)++;

      /* mark (A[I],B[J]) as used: put J into list row[I] */	
      l_ptr->I++; l_ptr->J++;
      z = ( pair_ptr ) ckalloc(1,sizeof(pair));
      z->COL = l_ptr->J;			
      z->NEXT = l_ptr->row[l_ptr->I];				
      l_ptr->row[l_ptr->I] = z;
      if (midj < N) INS(N-midj)
    }
    return midc;
  }

  /* Divide: Find optimum midpoint (midi,midj) of cost midc */

  midi = M/2;		/* Forward phase:                          */
  l_ptr->r_ss[0].CC = 0;		/*   Compute C(M/2,k) & D(M/2,k) for all k */
  t = -Q;
  for (j = 1; j <= N; j++) {
    l_ptr->r_ss[j].CC = t = t-R;
    l_ptr->r_ss[j].DD = t-Q;
  }
  t = -tb;
  for (i = 1; i <= midi; i++) {
    s = l_ptr->r_ss[0].CC;
    l_ptr->r_ss[0].CC = c = t = t-R;
    e = t-Q;
    va = pam2[A[i]];
    for (j = 1; j <= N; j++) {
      if ((c = c - qr) > (e = e - R)) e = c;
      if ((c = l_ptr->r_ss[j].CC - qr) > (d = l_ptr->r_ss[j].DD - R)) d = c;
      DIAG(i+l_ptr->I, j+l_ptr->J, c, s+va[B[j]])
      if (c < d) c = d;
      if (c < e) c = e;
      s = l_ptr->r_ss[j].CC;
      l_ptr->r_ss[j].CC = c;
      l_ptr->r_ss[j].DD = d;
    }
  }
  l_ptr->r_ss[0].DD = l_ptr->r_ss[0].CC;

  l_ptr->r_ss[N].RR = 0;			/* Reverse phase:                          */
  t = -Q;			/*   Compute R(M/2,k) & S(M/2,k) for all k */
  for (j = N-1; j >= 0; j--) {
    l_ptr->r_ss[j].RR = t = t-R;
    l_ptr->r_ss[j].SS = t-Q;
  }
  t = -te;

  for (i = M-1; i >= midi; i--) {
    s = l_ptr->r_ss[N].RR;
    l_ptr->r_ss[N].RR = c = t = t-R;
    e = t-Q;
    va = pam2[A[i+1]];
    for (j = N-1; j >= 0; j--) {
      if ((c = c - qr) > (e = e - R)) e = c;
      if ((c = l_ptr->r_ss[j].RR - qr) > (d = l_ptr->r_ss[j].SS - R)) d = c;
      DIAG(i+1+l_ptr->I, j+1+l_ptr->J, c, s+va[B[j+1]])
      if (c < d) c = d;
      if (c < e) c = e;
      s = l_ptr->r_ss[j].RR;
      l_ptr->r_ss[j].RR = c;
      l_ptr->r_ss[j].SS = d;
    }
  }
  l_ptr->r_ss[N].SS = l_ptr->r_ss[N].RR;

  midc = l_ptr->r_ss[0].CC+l_ptr->r_ss[0].RR;		/* Find optimal midpoint */
  midj = 0;
  type = 1;
  for (j = 0; j <= N; j++) {
    if ((c = l_ptr->r_ss[j].CC + l_ptr->r_ss[j].RR) >= midc) {
      if (c > midc ||
	  (l_ptr->r_ss[j].CC != l_ptr->r_ss[j].DD && l_ptr->r_ss[j].RR == l_ptr->r_ss[j].SS)) {
	midc = c;
	midj = j;
      }
    }
  }
  for (j = N; j >= 0; j--)
    if ((c = l_ptr->r_ss[j].DD + l_ptr->r_ss[j].SS + Q) > midc)
      { midc = c;
      midj = j;
      type = 2;
      }

/* Conquer: recursively around midpoint */

  if (type == 1) {
    (void)diff(A,B,midi,midj,tb,Q,pam2,Q,R, sapp, last, l_ptr);
    (void)diff(A+midi,B+midj,M-midi,N-midj,Q,te,pam2,Q,R, sapp, last, l_ptr);
  }
  else {
    (void)diff(A,B,midi-1,midj,tb,0,pam2,Q,R, sapp, last, l_ptr);
    DEL(2);
    (void)diff(A+midi+1,B+midj,M-midi-1,N-midj,0,te,pam2,Q,R, sapp, last, l_ptr);
  }
  return midc;
}

/* CHECK_SCORE - return the score of the alignment stored in S */

static int CHECK_SCORE(const unsigned char *A, const unsigned char *B,
		       int M, int N,
		       int *S, int **w,
		       int qq, int rr, int *NC)
{ 
  register int   i,  j, op, nc;
  int score;

  /*  print_seq_prof(A,M,w,iw); */

  score = i = j = op = nc = 0;
  while (i < M || j < N) {
    op = *S++;
    if (op == 0) {
      score = w[A[++i]][B[++j]] + score;
      nc++;
    }
    else if (op > 0) {
      score = score - (qq + op*rr);
      j += op;
      nc += op;
    } else {	/* op < 0 */
      score = score - (qq - op*rr);
      i -= op;	/* i increased */
      nc -= op;	/* nc increased */
    }
  }
  *NC = nc;
  return score;
}

/* ckalloc - allocate space; check for success */
void *ckalloc(size_t amount, size_t size)
{
  void *p;
  static size_t mtotal;

  mtotal += amount * size;

  if ((p = malloc( amount * size )) == NULL) {
    fprintf(stderr,"Ran out of near memory: %ld*%ld/%ld\n",amount,size,mtotal);
    exit(1);
  }
  return(p);
}