dropfx.c

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

      if (j >= 3) {
	sc = -score;
	e3 = e2-shift; e2 = last[j-3].C;
	e1 = last[j-2].C-shift; 
	if (e1 > sc) {sc = e1; del = 2;}
	if (e2 > sc) {sc = e2; del = 3;}
	if (e3 > sc) {sc = e3; del = 4;} 
      } else {
	sc = e2  = 0;
	if (sc < -score) sc=-score;
	else del = 3;
      }
      sc += score;
      if (sc < (ci=last[j].I)) {
	sc = ci; del = 0;
      }
      if (sc < (cd=cur[j].D)) {
	sc = cd; del = 5;
      }
      cur[j].C = sc;
      e = sc  - gop;
      if (e > cd) {
	cur[j+3].D = e-gext;
	cur_d_st[j+3] = 3;
      } else {
	cur[j+3].D = cd-gext;
	cur_d_st[j+3] = cur_d_st[j]+3;
      }
      switch(del) {
      case 5:
	e1 = cur_d_st[j];
	cur_st[j].i = cur_st[j-e1].i;
	cur_st[j].j = cur_st[j-e1].j;
	break;
      case 0:
	cur_st[j].i = cur_i_st[j].i;
	cur_st[j].j = cur_i_st[j].j;
	break;
      case 2:
      case 3:
      case 4:
	if (i) {
	  if (j-del >= 0) {
	    cur_st[j].i = last_st[j-del].i;
	    cur_st[j].j = last_st[j-del].j;
	  } else {
	    cur_st[j].i = i;
	    cur_st[j].j = 0;
	  }
	} else {
	  cur_st[j].i = 0;
	  cur_st[j].j = max(0, j-del+1);
	}
	break;
      case -1:
	cur_st[j].i = i+1;
	cur_st[j].j = j+1;
	break;
      }
      if (e > ci) {
	cur[j].I  = e -gext;
	cur_i_st[j].i = cur_st[j].i;
	cur_i_st[j].j = cur_st[j].j;
      } else {
	cur[j].I  = ci- gext;
      }
      if (sc > best) {
	x1 = cur_st[j].i;
	x2 = cur_st[j].j;
	best =sc;
	x3 = i;
	x4 = j;
      }
    }
    swap((void **)&last, (void **)&cur);
    swap((void **)&cur_st, (void **)&last_st);
  }
  /*	printf("The best score is %d\n", best); */
  *x = x1; *y = x2; *ex = x3; *ey = x4;
  free(cur_st); free(last_st); free(cur_i_st); 
  free(st_up);
  return best;
}

/* 
   Both global_up and global_down do linear space score only global 
   alignments on subsequence pro[x]...pro[ex], and dna[y]...dna[ey].
   global_up do the algorithm upwards, from row x towards row y.
   global_down do the algorithm downwards, from row y towards x.
*/

static void
global_up(st_ptr *row1, st_ptr *row2,
	  int x, int y, int ex, int ey, 
	  int **wgts, int gop, int gext, int shift,
	  unsigned char *dnap,
	  unsigned char *pro,
	  int N) {
  int i, j, k, sc, e, e1, e2, e3, t, ci, cd, score, *wt;
  st_ptr cur, last;

  cur = *row1; last = *row2;
  sc = -gop-gext;
  for (j = 1; j <= ey-y+1; j++) {
    if (j % 3 == 0) {last[j].C = sc; sc -= gext; last[j].I = sc-gop;}
    else { last[j].I = last[j].C = -10000;}
    cur[j].I = -10000;
  }  
  last[0].C = 0; cur[0].D = cur[1].D = cur[2].D = -10000;
  last[0].D = last[1].D = last[2].D = -10000;
  if (N) last[0].I = -gext; else last[0].I = -gop-gext;
  for (i = 1; i <= ex-x+1; i++) {
    wt = &wgts[pro[i+x-1]][0]; e2 = last[0].C; e1 = -10000;
    for (j = 0; j <= ey-y+1; j++) {
      t = j+y;
      sc = -10000; 
      if (t < 3) score = -10000;
      else score = wt[dnap[t-3]]; 
      if (j < 4) {
	if (j == 3) sc = e2;
	else if (j == 2) sc = e2-shift;
      } else {
	e3 = e2; e2 = e1;
	e1 = last[j-2].C;
	sc = max(max(e1, e3)-shift, e2);
      }
      sc += score;
      sc = max(sc, max(ci=last[j].I, cd = cur[j].D));
      cur[j].C = sc;
      cur[j+3].D = max(cd, sc-gop)-gext;
      cur[j].I = max(ci, sc-gop)-gext;
    }
    swap((void **)&last, (void **)&cur);
  }
  for (i = 0; i <= ey-y+1; i++) last[i].I = cur[i].I;
  if (*row1 != last) swap((void **)row1, (void **)row2);
}

static void
global_down(st_ptr *row1, st_ptr *row2,
	    int x, int y, int ex, int ey,
	    int **wgts, int gop, int gext, int shift,
	    unsigned char *dnap, unsigned char *pro,
	    int N) {
  int i, j, k, sc, del, *tmp, e,  t, e1,e2,e3, ci,cd, s1, s2, s3, *wt;
  st_ptr cur, last;

  cur = (*row1); last = *row2;
  sc = -gop-gext;
  for (j = ey-y; j >= 0; j--) {
    if ((ey-y+1-j) % 3) {last[j].C = sc; sc-=gext; last[j].I = sc-gop;}
    else  last[j].I =  last[j].C = -10000;
  } 
  last[ey-y+1].C = 0;
  cur[ey-y+1].D = cur[ey-y].D = cur[ey-y-1].D = -10000;
  last[ey-y+1].D = last[ey-y].D = last[ey-y-1].D = -10000;
  if (N) last[ey-y+1].I = -gext; else last[ey-y+1].I = -gop-gext;
  for (i = ex-x; i >= 0; i--) {
    wt = &wgts[pro[i+x]][0]; e2 = last[ey-y+1].C; 
    e1 = s2 = s3 = -10000; 
    for (j = ey-y+1; j >= 0; j--) {
      t = j+y;
      s1 = wt[dnap[t-1]];
      sc = -10000;
      if (t+3 > ey) {
	if (t+2==ey) sc = e2+s2;
	else if (t+1==ey) sc = e2-shift+s1;
      } else {
	e3 = e2; e2 = e1;
	e1 = last[j+2].C;
	sc = max(max(e1+s1, e3+s3)-shift, e2+s2);
      }
      if (sc < (cd= cur[j].D)) {
	sc = cd; 
	cur[j-3].D = cd-gext;
      } else cur[j-3].D =max(cd, sc-gop)-gext;
      if (sc < (ci= last[j].I)) {
	sc = ci; del = 0;
	cur[j].I = ci - gext;
      } else cur[j].I = max(sc-gop,ci)-gext;
      cur[j].C = sc;
      s3 = s2; s2 = s1;
    }
    swap((void **)&last, (void **)&cur);
  }
  for (i = 0; i <= ey-y+1; i++) last[i].I = cur[i].I;
  if (*row1 != last) swap((void **)row1, (void **)row2);
}

static void
init_row2(int *row, int ld) {
  int i;
  for (i = 0; i < ld; i++) row[i] = 0;
}

static void
init_ROW(st_ptr row, int ld) {
  int i;
  for (i = 0; i < ld; i++) row[i].I = row[i].D = row[i].C = 0;
}

static match_ptr
combine(match_ptr x1, match_ptr x2, int st) {
  match_ptr x;

  if (x1 == NULL) return x2;
  for (x = x1; x->next; x = x->next);
  x->next = x2;
  if (st) {
    for (x = x2; x; x = x->next) {
      x->j++;
      if (x->l == 3 || x->l == 4) break;
    }
    x->l--;
  }
  return x1;
}

/*
   global use the two upwards and downwards score only linear
   space global alignment subroutine to recursively build the
   alignment.
*/

match_ptr
global(int x, int y, int ex, int ey, 
       int **wgts, int gop, int gext, int shift,
       unsigned char *dnap, 
       unsigned char *pro,
       int N1, int N2,
       st_ptr *up_stp, st_ptr *dn_stp, st_ptr *tp_stp
       )
{
  int m;
  int m1, m2;
  match_ptr x1, x2, mm1, mm2;
  /*printf("%d %d %d %d\n", x,y, ex, ey);*/
  /*
    if the space required is limited, we can do a quadratic space
    algorithm to find the alignment.
  */
  if (ex <= x) {
    mm1  = NULL; mm2= NULL;
    for (m = y+3; m <= ey; m+=3) {
      x1 = (match_ptr) ckalloc(sizeof(match_node));
      x1->l = 5; x1->next = mm1; 
      if (mm1== NULL) mm2 = x1;
      mm1 = x1;
    }
    if (ex == x) {
      if ((ey-y) % 3 != 0) {
	x1 = (match_ptr) ckalloc(sizeof(match_node));
	x1->l = ((ey-y) % 3) +1; x1->next = NULL;
	if (mm2) mm2->next = x1;
	else mm1 = x1;
      } else {
	if (mm2) mm2->l = 4;
      }
    }
    return mm1;
  }
  if (ey <= y) {
    mm1  = NULL;
    for (m = x; m <= ex; m++) {
      x1 = (match_ptr) ckalloc(sizeof(match_node));
      x1->l = 0; x1->next = mm1; mm1 = x1;
    }
    return mm1;
  }
  if (ex -x < SGW1-1 && ey-y < SGW2-1) 
    return small_global(x,y,ex,ey,
			wgts, gop, gext, shift,
			dnap, pro, N1, N2);
  m = (x+ex)/2;
  /*     
	 Do the score only global alignment from row x to row m, m is
	 the middle row of x and ex. Store the information of row m in
	 upC, upD, and upI.
  */
  global_up(up_stp, tp_stp, x, y, m, ey, 
	    wgts, gop, gext, shift,
	    dnap, pro, N1);

  /* 
     Do the score only global alignment downwards from row ex
     to row m+1, store information of row m+1 in downC downI and downD
  */
  global_down(dn_stp, tp_stp,  m+1, y, ex, ey, 
	      wgts, gop, gext, shift,
	      dnap, pro, N2);

  /*
    Use these information of row m and m+1, to find the crossing
    point of the best alignment with the middle row. The crossing
    point is given by m1 and m2. Then we recursively call global
    itself to compute alignments in two smaller regions found by
    the crossing point and combine the two alignments to form a
    whole alignment. Return that alignment.
  */
  if (find_best(*up_stp, *dn_stp, &m1, &m2, ey-y+1, y, gop)) {
    x1 = global(x, y, m, m1, wgts, gop, gext, shift, dnap, pro, N1, 0,
		up_stp, dn_stp, tp_stp);
    x2 = global(m+1, m2, ex, ey, wgts, gop, gext, shift, dnap, pro, 0, N2,
		up_stp, dn_stp, tp_stp);
    if (m1 == m2) x1 = combine(x1,x2,1);
    else x1 = combine(x1, x2,0);
  } else {
    x1 = global(x, y, m-1, m1, wgts, gop, gext, shift,  dnap, pro, N1, 1,
		up_stp, dn_stp, tp_stp);
    x2 = global(m+2, m2, ex, ey, wgts, gop, gext, shift, dnap, pro, 1, N2,
		up_stp, dn_stp, tp_stp);
    mm1 = (match_ptr) ckalloc(sizeof(match_node));
    mm1->i = m; mm1->l = 0; mm1->j = m1;
    mm2 = (match_ptr) ckalloc(sizeof(match_node));
    mm2->i = m+1; mm2->l = 0; mm2->j = m1;
    mm1->next = mm2; mm2->next = x2;
    x1 = combine(x1, mm1, 0);
  }
  return x1;
}

static int
find_best(st_ptr up, st_ptr down,
	  int *m1, int *m2,
	  int ld, int y, int gop) {
  int i, best = -100000, j = 0, s1, s2, s3, s4, st;
  up++;
  for (i = 1; i < ld; i++) {
    s2 = up[i-1].C + down[i].C;
    s4 = up[i-1].I + down[i].I + gop;
    if (best < s2) {
      best = s2; j = i; st = 1;
    }
    if (best < s4) {
      best = s4; j = i; st = 0;
    }
  }
  *m1 = j-1+y;
  *m2 = j+y;
  /*printf("find best score =%d\n", best);*/
  return st;
} 

/*
   An alignment is represented as a linked list whose element
   is of type match_node. Each element represent an edge in the
   path of the alignment graph. The fields of match_node are
   l ---  gives the type of the edge.
   i, j --- give the end position.
*/

static match_ptr
small_global(int x, int y, int ex, int ey,
	     int **wgts, int gop, int gext, int shift,
	     unsigned char *dnap, unsigned char *pro,
	     int N1, int N2) {
  static int C[SGW1+1][SGW2+1], st[SGW1+1][SGW2+1], D[SGW2+7], I[SGW2+1];
  int i, j, e, sc, score, del, k, t, *wt, ci, cd;
  int *cI, *cD, *cC, *lC, *cst, e2, e3, e4;
  match_ptr mp, first;

  /*printf("small_global %d %d %d %d\n", x, y, ex, ey);*/
  sc = -gop-gext; C[0][0] = 0;
  if (N1) I[0] = -gext; else I[0] = sc;
  for (j = 1; j <= ey-y+1; j++) {
    if (j % 3== 0) {
      C[0][j] = sc; sc -= gext; I[j] = sc-gop;
    } else I[j] = C[0][j] = -10000;
    st[0][j] = 5;
  }
  lC = &C[0][0]; cD = D; D[0] = D[1] = D[2] = -10000;
  cI = I;
  for (i = 1; i <= ex-x+1; i++) {
    cC = &C[i][0];	
    wt = &wgts[pro[i+x-1]][0]; cst = &st[i][0];
    for (j = 0; j <=ey-y+1; j++) {
      sc = -10000; del = 0;
      ci = cI[j];
      cd= cD[j];
      t = j+y;
      if (t < 3) score = -10000;
      else score = wt[dnap[t-3]];
      if (j >= 4) {
	e2 = lC[j-2]-shift; sc = lC[j-3]; e4 = lC[j-4]-shift;
	del = 3;
	if (e2 > sc) { sc = e2; del = 2;}
	if (e4 >= sc) { sc = e4; del = 4;}
      } else {
	if (j ==3) {sc= lC[0]; del = 3;}
	else if (j == 2) {sc = lC[0]-shift; del = 2;}
      }
      sc = sc+score;
      if (sc < ci) {
	sc = ci; del = 0; 
      }
      if (sc <= cd) {
	sc = cd;
	del = 5;
      }
      cC[j] = sc;
      sc -= gop;
      if (sc < cd) {
	del += 10;
	cD[j+3] = cd - gext;
      } else cD[j+3] = sc -gext;
      if (sc < ci) {
	del += 20;
	cI[j] = ci-gext;
      } else cI[j] = sc-gext;
      *(cst++) = del;
    }
    lC = cC;
  }
  if (N2 && ci +gop > cC[ey-y+1]) {
    st[ex-x+1][ey-y+1] = 0;
    /*printf("small score = %d\n", ci+gop);*/
  } /*else printf("small score =%d\n", cC[ey-y+1]);*/
  first = NULL; e = 1;
  for (i = ex+1, j = ey+1; i > x || j > y; i--) {
    mp = (match_ptr) ckalloc(sizeof(match_node));
    mp->i = i-1;
    k  = (t=st[i-x][j-y])%10;
    mp->j = j-1;
    if (e == 5 && (t/10)%2 == 1) k = 5;
    if (e == 0 && (t/20)== 1) k = 0;
    if (k == 5) { j -= 3; i++; e=5;}
    else {j -= k;if (k==0) e= 0; else e = 1;}
    mp->l = k;
    mp->next = first;
    first = mp;
  }

  /*	for (i = 0; i <= ex-x; i++) {
	for (j = 0; j <= ey-y; j++) 
	printf("%d ", C[i][j]);
	printf("\n");
	}
  */
  return first;	
}


#define XTERNAL
#include "upam.h"

extern void display_alig(a, dna, pro,length, ld)
int *a;
unsigned char *dna, *pro;
int length, ld;
{
	int len = 0, i, j, x, y, lines, k;
	static char line1[100], line2[100], line3[100],
		 tmp[10] = "         ";
	unsigned char *dna1, c1, c2, c3, *st;

	dna1 = ckalloc((size_t)ld);
	for (st = dna, i = 0; i < ld; i++, st++) dna1[i] = aa[*st];
	line1[0] = line2[0] = line3[0] = '\0'; x= a[0]; y = a[1]-1;
 
	for (len = 0, j = 2, lines = 0; j < length; j++) {
		i = a[j];
		/*printf("%d %d %d\n", i, len, b->j);*/
		if (i > 0 && i < 5) tmp[i-2] = aa[pro[x++]];
		if (i == 5) {
		    i = 3; tmp[0] = tmp[1] = tmp[2] = '-';
		    if (a[j+1] == 2) tmp[2] = ' ';
		}
		if (i > 0) {
		    strncpy(&line1[len], (const char *)&dna1[y], i); y+=i;
		} else {line1[len] = '-'; i = 1; tmp[0] = aa[pro[x++]];}
		strncpy(&line2[len], tmp, i);
		for (k = 0; k < i; k++) {
			if (tmp[k] != ' ' && tmp[k] != '-') {
				if (k == 2) tmp[k] = '\\';
				else if (k == 1) tmp[k] = '|';
				else tmp[k] = '/';
			} else tmp[k] = ' ';
		}
		if (i == 1) tmp[0] = ' ';
		strncpy(&line3[len], tmp, i); 
		tmp[0] = tmp[1] =  tmp[2] = ' ';
		len += i;
		line1[len] = line2[len] =line3[len]  = '\0'; 
		if (len >= WIDTH) {
		    printf("\n%5d", WIDTH*lines++);
		    for (k = 10; k <= WIDTH; k+=10) 
			printf("    .    :");
		    if (k-5 < WIDTH) printf("    .");
		    c1 = line1[WIDTH]; c2 = line2[WIDTH]; c3 = line3[WIDTH];
		    line1[WIDTH] = line2[WIDTH] = line3[WIDTH] = '\0';
		    printf("\n     %s\n     %s\n     %s\n", line1, line3, line2);
		    line1[WIDTH] = c1; line2[WIDTH] = c2; line3[WIDTH] = c3;
		    strncpy(line1, &line1[WIDTH], sizeof(line1)-1);
		    strncpy(line2, &line2[WIDTH], sizeof(line2)-1);
		    strncpy(line3, &line3[WIDTH], sizeof(line3)-1);
		    len = len - WIDTH;
		}
        }
	printf("\n%5d", WIDTH*lines);
	for (k = 10; k < len; k+=10) 
	    printf("    .    :");
	if (k-5 < len) printf("    .");
	printf("\n     %s\n     %s\n     %s\n", line1, line3, line2);
}