dropnfa.c

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

	    break;
	 }
      }

/*	now recalculate where the score fits */
      if (start == NULL)
	 start = &sarr[si];
      else
	 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;
	 }
      si++;
   }

   if (start != NULL)
      return (start->score);
   else
      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;
	 }
}

int dmatch (const unsigned char *aa0, int n0,
		   const unsigned char *aa1, int n1,
		   int hoff, int window, 
		   int **pam2, int gdelval, int ggapval,
		   struct f_struct *f_str)
{
   int low, up;

   window = min (n1, window);
   /* hoff is the offset found from aa1 to seq 2 by hmatch */

   low = -window/2-hoff;
   up = low+window;

   return FLOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
		      pam2,
#ifdef OLD_FASTA_GAP
		       -(gdelval-ggapval),
#else
		       -gdelval,
#endif
		       -ggapval,window,f_str);
 }


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

   To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,MW).
   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
	MW  : maximum window size
*/

#include <stdio.h>

#define MININT -9999999

int
FLOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
	     int M, int N, int low, int up,
	     int **W, int G,int H, int MW,
	     struct f_struct *f_str)
{
  int band;
  register struct bdstr *bssp;
  int i, j, si, ei;
  int c, d, e, m;
  int leftd, rightd;
  int best_score;
  int *wa, curd;
  int ib;
  
  bssp = f_str->bss;

  m = G+H;
  low = max(-M, low);
  up = min(N, up);
  
  if (N <= 0) return 0;

  if (M <= 0) return 0;

  band = up-low+1;
  if (band < 1) {
    fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
	    M, N, low, up);
    return 0;
  }

  if (low > 0) leftd = 1;
  else if (up < 0) leftd = band;
  else leftd = 1-low;
  rightd = band;
  si = max(0,-up);	/* start index -1 */
  ei = min(M,N-low);	/* end index */
  bssp[leftd].CC = 0;
  for (j = leftd+1; j <= rightd; j++) {
    bssp[j].CC = 0;
    bssp[j].DD = -G;
  }

  bssp[rightd+1].CC = MININT;
  bssp[rightd+1].DD = MININT;

  best_score = 0;
  bssp[leftd-1].CC = MININT;
  bssp[leftd].DD = -G;

  for (i = si+1; i <= ei; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = W[A[i]];
    if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
    if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];

    if (d > c) c = d;
    if (c < 0) c = 0;
    e = c-G;
    bssp[leftd].DD = d;
    bssp[leftd].CC = c;
    if (c > best_score) best_score = c;

    for (curd=leftd+1; curd <= rightd; curd++) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      if (c < 0) c = 0;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c > best_score) best_score = c;
    }
  }

  return best_score;
}

/* ckalloc - allocate space; check for success */
char *ckalloc(size_t amount)
{
  char *p;

  if ((p = malloc( (unsigned) amount)) == NULL)
    w_abort("Ran out of memory.","");
  return(p);
}

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

static int
BCHECK_SCORE(const unsigned char *A, const unsigned char *B,
		       int M, int N, int *S, int **w, int g, int h,
		       int *nres)
{ 
  register int   i,  j, op, nc;
  int *Ssave;
  int score;

  score = i = j = op = nc = 0;
  Ssave = S;
  while (i < M || j < N) {
	op = *S++;
	if (op == 0) {
	  score = w[A[++i]][B[++j]] + score;
	  nc++;
/*  	  fprintf(stderr,"op0 %4d %4d %4d %4d\n",i,j,w[A[i]][B[i]],score); */
	}
	else if (op > 0) {
	  score = score - (g+op*h);
/*  	  fprintf(stderr,"op> %4d %4d %4d %4d %4d\n",i,j,op,-(g+op*h),score); */
	  j += op;
	  nc += op;
	} else {
	  score = score - (g-op*h);
/*  	  fprintf(stderr,"op< %4d %4d %4d %4d %4d\n",i,j,op,-(g-op*h),score); */
	  i -= op;
	  nc -= op;
	}
  }
  *nres = nc;
  return score;
}


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

   To invoke, call LOCAL_ALIGN(A,B,M,N,L,U,W,G,H,S,dflag,&SI,&SJ,&EI,&EJ,MW).
   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
	dflag : 0 - no display or backward pass
	*SI : starting position of sequence A in the optimal local alignment
	*SJ : starting position of sequence B in the optimal local alignment
	*EI : ending position of sequence A in the optimal local alignment
	*EJ : ending position of sequence B in the optimal local alignment
	MW  : maximum window size
*/

int bd_walign (const unsigned char *aa0, int n0,
		 const unsigned char *aa1, int n1,
		 struct pstruct *ppst, 
		 struct f_struct *f_str, int hoff,
		 struct a_res_str *a_res)
{
   int low, up, score;
   int min0, min1, max0, max1;
   int window;

   window = min (n1, ppst->param_u.fa.optwid);
   /* hoff is the offset found from aa1 to seq 2 by hmatch */

   low = -window/2-hoff;
   up = low+window;

   score=LOCAL_ALIGN(aa0-1,aa1-1,n0,n1, low, up,
		    ppst->pam2[0],
#ifdef OLD_FASTA_GAP
		     -(ppst->gdelval-ppst->ggapval),
#else
		     -ppst->gdelval,
#endif
		     -ppst->ggapval,
		    &min0,&min1,&max0,&max1,ppst->param_u.fa.optwid,f_str);
  
   if (score <=0) {
     fprintf(stderr,"n0/n1: %d/%d hoff: %d window: %d\n",
	     n0, n1, hoff, window);
     return 0;
   }
  
/*
  fprintf(stderr," ALIGN: start0: %d start1: %d stop0: %d stop1: %d, bot: %d top: %d, win: %d MX %d\n",
  min0-1,min1-1,max0-min0+1,max1-min1+1,low-(min1-min0),up-(min1-min0),
  ppst->param_u.fa.optwid,n0);
*/

   a_res->min0 = min0-1; a_res->min1 = min1-1;
   a_res->max0 = max0; a_res->max1 = max1; 

   B_ALIGN(aa0-1+min0-1,aa1-1+min1-1,max0-min0+1,max1-min1+1,
	   low-(min1-min0),up-(min1-min0),
	   ppst->pam2[0],
#ifdef OLD_FASTA_GAP
	   -(ppst->gdelval-ppst->ggapval),
#else
	   -ppst->gdelval,
#endif
	   -ppst->ggapval,
	   a_res->res,&a_res->nres,ppst->param_u.fa.optwid,n0,f_str->bss);

   return score;
}

static int
LOCAL_ALIGN(const unsigned char *A, const unsigned char *B,
	    int M, int N,
	    int low, int up, int **W, int G,int H,
	    int *psi, int *psj, int *pei, int *pej, int MW,
	    struct f_struct *f_str)
{ 
  int band;
  register struct bdstr *bssp;
  int i, j, si, ei;
  int c, d, e, t, m;
  int leftd, rightd;
  int best_score, starti, startj, endi, endj;
  int *wa, curd;
  int ib;
  char flag;
  
  bssp = f_str->bss;

  m = G+H;
  low = max(-M, low);
  up = min(N, up);
  
  if (N <= 0) { 
    *psi = *psj = *pei = *pej;
    return 0;
  }
  if (M <= 0) {
    *psi = *psj = *pei = *pej;
    return 0;
  }
  band = up-low+1;
  if (band < 1) {
    fprintf(stderr,"low > up is unacceptable!: M: %d N: %d l/u: %d/%d\n",
	    M, N, low, up);
    return -1;
  }

  j = (MW + 2 + 2) * sizeof(struct bdstr);

  /* already done by init_work(); 
  if (f_str->bss==NULL) f_str->bss = (struct bdstr *) ckalloc(j);
  */

  if (low > 0) leftd = 1;
  else if (up < 0) leftd = band;
  else leftd = 1-low;
  rightd = band;
  si = max(0,-up);
  ei = min(M,N-low);
  bssp[leftd].CC = 0;
  for (j = leftd+1; j <= rightd; j++) {
    bssp[j].CC = 0;
    bssp[j].DD = -G;
  }
  bssp[rightd+1].CC = MININT;
  bssp[rightd+1].DD = MININT;
  best_score = 0;
  endi = si;
  endj = si+low;
  bssp[leftd-1].CC = MININT;
  bssp[leftd].DD = -G;
  for (i = si+1; i <= ei; i++) {
    if (i > N-up) rightd--;
    if (leftd > 1) leftd--;
    wa = W[A[i]];
    if ((c = bssp[leftd+1].CC-m) > (d = bssp[leftd+1].DD-H)) d = c;
    if ((ib = leftd+low-1+i ) > 0) c = bssp[leftd].CC+wa[B[ib]];
/*
    if (ib > N) fprintf(stderr,"B[%d] out of range %d\n",ib,N);
*/
    if (d > c) c = d;
    if (c < 0) c = 0;
    e = c-G;
    bssp[leftd].DD = d;
    bssp[leftd].CC = c;
    if (c > best_score) {
      best_score = c;
      endi = i;
      endj = ib;
    }
    for (curd=leftd+1; curd <= rightd; curd++) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd+1].CC-m) > (d = bssp[curd+1].DD-H)) d = c;
/*
      if ((ib=curd+low-1+i) <= 0 || ib > N)
	fprintf(stderr,"B[%d]:%d\n",ib,B[ib]);
*/
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      if (c < 0) c = 0;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c > best_score) {
	best_score = c;
	endi = i;
	endj = curd+low-1+i;
      }
    }
  }
  
  leftd = max(1,-endi-low+1);
  rightd = band-(up-(endj-endi));
  bssp[rightd].CC = 0;
  t = -G;
  for (j = rightd-1; j >= leftd; j--) {
    bssp[j].CC = t = t-H;
    bssp[j].DD = t-G;
  }
  for (j = rightd+1; j <= band; ++j) bssp[j].CC = MININT;
  bssp[leftd-1].CC = bssp[leftd-1].DD = MININT;
  bssp[rightd].DD = -G;
  flag = 0;
  for (i = endi; i >= 1; i--) {
    if (i+low <= 0) leftd++;
    if (rightd < band) rightd++;
    wa = W[A[i]];
    if ((c = bssp[rightd-1].CC-m) > (d = bssp[rightd-1].DD-H)) d = c;
    if ((ib = rightd+low-1+i) <= N) c = bssp[rightd].CC+wa[B[ib]];

/*
    if (ib <= 0) fprintf(stderr,"rB[%d] <1\n",ib);
*/
    if (d > c) c = d;
    e = c-G;
    bssp[rightd].DD = d;
    bssp[rightd].CC = c;
    if (c == best_score) {
      starti = i;
      startj = ib;
      flag = 1;
      break;
    }
    for (curd=rightd-1; curd >= leftd; curd--) {
      if ((c = c-m) > (e = e-H)) e = c;
      if ((c = bssp[curd-1].CC-m) > (d = bssp[curd-1].DD-H)) d = c;

/*
      if ((ib=curd+low-1+i) <= 0 || ib > N)
	fprintf(stderr,"i: %d, B[%d]:%d\n",i,ib,B[ib]);
*/
      c = bssp[curd].CC + wa[B[curd+low-1+i]];
      if (e > c) c = e;
      if (d > c) c = d;
      bssp[curd].CC = c;
      bssp[curd].DD = d;
      if (c == best_score) {
	starti = i;
	startj = curd+low-1+i;
	flag = 1;
	break;
      }
    }
    if (flag == 1) break;
  }
  
  if (starti < 0 || starti > M || startj < 0 || startj > N) {
    printf("starti=%d, startj=%d\n",starti,startj);
    *psi = *psj = *pei = *pej;
    exit(1);
  }
  *psi = starti;
  *psj = startj;
  *pei = endi;
  *pej = endj;
  return best_score;
}