dropnfa.c

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

/* A PACKAGE FOR GLOBALLY ALIGNING TWO SEQUENCES WITHIN A BAND:

   To invoke, call B_ALIGN(A,B,M,N,L,U,W,G,H,S,MW,MX).
   The parameters are explained as follows:
	A, B : two sequences to be aligned
	M : the length of sequence A
	N : the length of sequence B
	L : lower bound of the band
	U : upper bound of the band
	W : scoring table for matches and mismatches
	G : gap-opening penalty
	H : gap-extension penalty
	S : script for DISPLAY routine
	MW : maximum window size
	MX : maximum length sequence M to be aligned
*/

static int IP;
static int *MP[3];		/* save crossing points */
static int *FP;			/* forward dividing points */
static char *MT[3];		/* 0: rep, 1: del, 2: ins */
static char *FT;

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

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

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

#define REP { *last = (*sapp)[0] = 0; (*sapp)++;} /* Append "Replace" op */

/* bg_align(A,B,M,N,up,low,tb,te) returns the cost of an optimum conversion between
  A[1..M] and B[1..N] and appends such a conversion to the current script.
  tb(te)= 1  no gap-open penalty if the conversion begins(ends) with a delete.
  tb(te)= 2  no gap-open penalty if the conversion begins(ends) with an insert.
*/
static int
bg_align(const unsigned char *A, const unsigned char *B, 
	 int M, int N,
	 int low, int up, int tb, int te,
	 int **w, int g, int h,
	 struct bdstr *bss, int **sapp, int *last)
{
  int rmid, k, l, r, v, kt;
  int t1, t2, t3;

  {
  int band, midd;
  int leftd, rightd;	/* for CC, DD, CP and DP */
  register int curd;	/* current index for CC, DD CP and DP */
  register int i, j;
  register int c, d, e;
  int t, fr, *wa, ib, m;

  /* Boundary cases: M <= 0 , N <= 0, or up-low <= 0 */
  if (N <= 0) { 
    if (M > 0) { DEL(M) }
    return 0;
  }
  if (M <= 0) {
    INS(N)
    return 0;
  }
  if ((band = up-low+1) <= 1) {
    for (i = 1; i <= M; i++) { REP }
    return 0;
  }

  /* Divide: Find all crossing points */

  /* Initialization */
  m = g + h;

  midd = band/2 + 1;
  rmid = low + midd - 1;
  leftd = 1-low;
  rightd = up-low+1;
  if (leftd < midd) {
    fr = -1;
    for (j = 0; j < midd; j++) 
      bss[j].CP = bss[j].DP = -1;
    for (j = midd; j <= rightd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    MP[0][0] = -1;
    MP[1][0] = -1;
    MP[2][0] = -1;
    MT[0][0] = MT[1][0] = MT[2][0] = 0;
  } else if (leftd > midd) {
    fr = leftd-midd;
    for (j = 0; j <= midd; j++) {
      bss[j].CP = bss[j].DP = fr;
    }
    for (j = midd+1; j <= rightd; j++) 
      bss[j].CP = bss[j].DP = -1;
    MP[0][fr] = -1;
    MP[1][fr] = -1;
    MP[2][fr] = -1;
    MT[0][fr] = MT[1][fr] = MT[2][fr] = 0;
  } else {
    fr = 0;
    for (j = 0; j < midd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    for (j = midd; j <= rightd; j++) {
      bss[j].CP = bss[j].DP = 0;
    }
    MP[0][0] = -1;
    MP[1][0] = -1;
    MP[2][0] = -1;
    MT[0][0] = MT[1][0] = MT[2][0] = 0;
  }

  bss[leftd].CC = 0;
  if (tb == 2) t = 0;
  else t = -g;
  for (j = leftd+1; j <= rightd; j++) {
    bss[j].CC = t = t-h;
    bss[j].DD = t-g;
  }
  bss[rightd+1].CC = MININT;
  bss[rightd+1].DD = MININT;
  if (tb == 1) bss[leftd].DD = 0;
  else bss[leftd].DD = -g;
  bss[leftd-1].CC = MININT;
  for (i = 1; i <= M; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = w[A[i]];
    if ((c = bss[leftd+1].CC-m) > (d = bss[leftd+1].DD-h)) {
      d = c;
      bss[leftd].DP = bss[leftd+1].CP;
    } else bss[leftd].DP = bss[leftd+1].DP;
    if ((ib = leftd+low-1+i) > 0) c = bss[leftd].CC+wa[B[ib]];
    if (d > c || ib <= 0) {
      c = d;
      bss[leftd].CP = bss[leftd].DP;
    }
    e = c-g;
    bss[leftd].DD = d;
    bss[leftd].CC = c;
    IP = bss[leftd].CP;
    if (leftd == midd) bss[leftd].CP = bss[leftd].DP = IP = i;
    for (curd=leftd+1; curd <= rightd; curd++) {
      if (curd != midd) {
	if ((c = c-m) > (e = e-h)) {
	  e = c;
	  IP = bss[curd-1].CP;
	}  /* otherwise, IP is unchanged */
	if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
	  d = c;
	  bss[curd].DP = bss[curd+1].CP;
	} else {
	  bss[curd].DP = bss[curd+1].DP;
	}
	c = bss[curd].CC + wa[B[curd+low-1+i]];
	if (c < d || c < e) {
	  if (e > d) {
	    c = e;
	    bss[curd].CP = IP;
	  } else {
	    c = d;
	    bss[curd].CP = bss[curd].DP;
	  }
	} /* otherwise, CP is unchanged */
	bss[curd].CC = c;
	bss[curd].DD = d;
      } else {
	if ((c = c-m) > (e = e-h)) {
	  e = c;
	  MP[1][i] = bss[curd-1].CP;
	  MT[1][i] = 2;
	} else {
	  MP[1][i] = IP;
	  MT[1][i] = 2;
	}
	if ((c = bss[curd+1].CC-m) > (d = bss[curd+1].DD-h)) {
	  d = c;
	  MP[2][i] = bss[curd+1].CP;
	  MT[2][i] = 1;
	} else {
	  MP[2][i] = bss[curd+1].DP;
	  MT[2][i] = 1;
	}
	c = bss[curd].CC + wa[B[curd+low-1+i]];
	if (c < d || c < e) {
	  if (e > d) {
	    c = e;
	    MP[0][i] = MP[1][i];
	    MT[0][i] = 2;
	  } else {
	    c = d;
	    MP[0][i] = MP[2][i];
	    MT[0][i] = 1;
	  }
	} else {
	  MP[0][i] = i-1;
	  MT[0][i] = 0;
	}
	if (c-g > e) {
	  MP[1][i] = MP[0][i];
	  MT[1][i] = MT[0][i];
	}
	if (c-g > d) {
	  MP[2][i] = MP[0][i];
	  MT[2][i] = MT[0][i];
	}
	bss[curd].CP = bss[curd].DP = IP = i;
	bss[curd].CC = c;
	bss[curd].DD = d;
      }
    }
  }

  /* decide which path to be traced back */
  if (te == 1 && d+g > c) {
    k = bss[rightd].DP;
    l = 2;
  } else if (te == 2 && e+g > c) {
    k = IP;
    l = 1;
  } else {
    k = bss[rightd].CP;
    l = 0;
  }
  if (rmid > N-M) l = 2;
  else if (rmid < N-M) l = 1;
  v = c;
  }
  /* Conquer: Solve subproblems recursively */

  /* trace back */
  r = -1;	
  for (; k > -1; r=k, k=MP[l][r], l=MT[l][r]){
    FP[k] = r;
    FT[k] = l;	/* l=0,1,2 */
  }
  /* forward dividing */
  if (r == -1) { /* optimal alignment did not cross the middle diagonal */
    if (rmid < 0) {
      bg_align(A,B,M,N,rmid+1,up,tb,te,w,g,h,bss, sapp, last);
    }
    else {
      bg_align(A,B,M,N,low,rmid-1,tb,te,w,g,h,bss, sapp, last);
    }
  } else {
    k = r;
    l = FP[k];
    kt = FT[k];

    /* first block */
    if (rmid < 0) {
      bg_align(A,B,r-1,r+rmid,rmid+1,min(up,r+rmid),tb,1,w,g,h,bss,sapp,last);
      DEL(1)
    } else if (rmid > 0) {
      bg_align(A,B,r,r+rmid-1,max(-r,low),rmid-1,tb,2,w,g,h,bss,sapp,last);
      INS(1)
    }

    /* intermediate blocks */
    t2 = up-rmid-1;
    t3 = low-rmid+1;
    for (; l > -1; k = l, l = FP[k], kt = FT[k]) {
      if (kt == 0) { REP }
      else if (kt == 1) { /* right-hand side triangle */
	INS(1)
	t1 = l-k-1;
	bg_align(A+k,B+k+rmid+1,t1,t1,0,min(t1,t2),2,1,w,g,h,bss,sapp,last);
	DEL(1)
      }
      else { /* kt == 2, left-hand side triangle */
	DEL(1)
	t1 = l-k-1;
	bg_align(A+k+1,B+k+rmid,t1,t1,max(-t1,t3),0,1,2,w,g,h,bss,sapp,last);
	INS(1)
      }
    }

    /* last block */
    if (N-M > rmid) {
      INS(1)
      t1 = k+rmid+1;
      bg_align(A+k,B+t1,M-k,N-t1,0,min(N-t1,t2),2,te,w,g,h,bss,sapp,last);
    } else if (N-M < rmid) {
      DEL(1)
      t1 = M-(k+1);
      bg_align(A+k+1,B+k+rmid,t1,N-(k+rmid),max(-t1,t3),0,1,te,w,g,h,
	       bss,sapp,last);
    }
  }
  return(v);
}

int B_ALIGN(const unsigned char *A, const unsigned char *B,
	    int M, int N,
	    int low, int up, int **W, int G, int H, int *S, int *nS,
	    int MW, int MX, struct bdstr *bss)
{ 
  int c, i, j;
  int g, h;
  size_t mj;
  int check_score;
  int **sapp, *sapp_v, *last, last_v;

  g = G;
  h = H;
  sapp_v = S;
  sapp = &sapp_v;

  last_v = 0;
  last = &last_v;

  low = min(max(-M, low),min(N-M,0));
  up = max(min(N, up),max(N-M,0));

  if (N <= 0) { 
    if (M > 0) { DEL(M); }
    return -gap(M);
  }
  if (M <= 0) {
    INS(N);
    return -gap(N);
  }
  if (up-low+1 <= 1) {
    c = 0;
    for (i = 1; i <= M; i++) {
      REP;
      c += W[A[i]][B[i]];
    }
    return c;
  }

  if (MT[0]==NULL) {
    mj = MX+1;
    MT[0] = (char *) ckalloc(mj);
    MT[1] = (char *) ckalloc(mj);
    MT[2] = (char *) ckalloc(mj);
    FT = (char *) ckalloc(mj);

    mj *= sizeof(int);
    MP[0] = (int *) ckalloc(mj);
    MP[1] = (int *) ckalloc(mj);
    MP[2] = (int *) ckalloc(mj);
    FP = (int *) ckalloc(mj);
  }

  c = bg_align(A,B,M,N,low,up,0,0,W,G,H,bss, sapp, last);

  check_score = BCHECK_SCORE(A,B,M,N,S,W,G,H,nS);

  free(FP); free(MP[2]); free(MP[1]); free(MP[0]);
  free(FT); free(MT[2]); free(MT[1]); free(MT[0]);
  MT[0]=NULL;

  if (check_score != c)
    printf("\nBCheck_score=%d != %d\n", check_score,c);
  return c;
}

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, optflag_s, optcut_s, optwid_s, n10, score;
  const unsigned char *aa1p;
  struct rstruct rst;

#ifdef TFASTA  
  f_str->n10 = n10=aatran(aa1,f_str->aa1x,n1,frame);
  do_fasta (aa0, n0, f_str->aa1x, n10, ppst, f_str, &rst, &hoff);
  aa1p = f_str->aa1x;

#else
  n10 = n1;
  aa1p = aa1;
#endif

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

  if (ppst->sw_flag) {
  a_res->sw_score = sw_walign(f_str->pam2p[0], n0, aa1, n1, 
#ifdef OLD_FASTA_GAP
			      -(ppst->gdelval - ppst->ggapval),
#else
			      -ppst->gdelval,
#endif
			      -ppst->ggapval,
			      f_str->ss, a_res);
  }
  else {
    optflag_s = ppst->param_u.fa.optflag;
    optcut_s = ppst->param_u.fa.optcut;
    optwid_s = ppst->param_u.fa.optwid;
    ppst->param_u.fa.optflag = 1;
    ppst->param_u.fa.optcut = 0;
    ppst->param_u.fa.optwid *= 2;

    do_fasta(aa0, n0, aa1p, n10, ppst, f_str, &rst, &hoff);

    if (rst.score[0]>0) {
      score=bd_walign(aa0, n0, aa1p, n10, ppst, f_str, hoff, a_res);
    }
    else {
      a_res->nres = 0;
      score=0;
    }

    ppst->param_u.fa.optflag = optflag_s;
    ppst->param_u.fa.optcut = optcut_s;
    ppst->param_u.fa.optwid = optwid_s;
    a_res->sw_score = score;
  }
}

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

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

#ifdef PCOMPLIB

#include "w_mw.h"

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