1
2 #include <stdio.h>
3 #include <stdlib.h>
4
5
6 struct sx_s {int C1, C2, C3, I1, I2, I3, flag; };
7
8 static struct sx_s *cur=NULL;
9
10 #define max(x,y) ((x) > (y) ? (x) : (y))
11 #define min(x,y) ((x) < (y) ? (x) : (y))
12
13 static init_row(row, sp)
14 struct sx_s *row;
15 int sp;
16 {
17 int i;
18 for (i = 0; i < sp; i++) {
19 row[i].C1 = row[i].I1 = 0;
20 row[i].C2 = row[i].I2 = 0;
21 row[i].C3 = row[i].I3 = 0;
22 row[i].flag = 0;
23 }
24 }
25
lx_align(char * prot_seq,int len_prot,char * dna_prot_seq,int len_dna_prot,int ** pam_matrix,int gopen,int gext,int gshift,int start_diag,int width)26 extern lx_align(char *prot_seq, /* array with protein sequence numbers*/
27 int len_prot, /* length of prot. seq */
28 char *dna_prot_seq, /* translated DNA sequence numbers*/
29 int len_dna_prot, /* length trans. seq. */
30 int **pam_matrix, /* scoring matrix */
31 int gopen, int gext, /* gap open, gap extend penalties */
32 int gshift, /* frame-shift penalty */
33 int start_diag, /* start diagonal of band */
34 int width) /* width for band alignment */
35 {
36 char *ckalloc();
37 int i, j, bd, bd1, x1, x2, sp, p1=0, p2=0;
38 int sc, del, best = 0, cd,ci, e1, e2, e3, cd1, cd2, cd3, f, gg;
39 register int *wt;
40 register char *dp;
41 register struct sx_s *ap, *aq;
42
43 sp = width+7;
44 gg = gopen+gext;
45 /* sp = sp/3; */
46 if (cur == NULL)
47 cur = (struct sx_s *) ckalloc(sizeof(struct sx_s)*sp);
48
49 init_row(cur, sp);
50
51 /*
52 if (start_diag %3 !=0) start_diag = start_diag/3-1;
53 else start_diag = start_diag/3;
54 */
55
56 /*
57 if (width % 3 != 0) width = width/3+1;
58 else width = width /3;
59 */
60
61 x1 = start_diag; /* x1 = lower bound of DNA */
62 x2 = 1; /* the amount of position shift from last row*/
63
64 /* i counts through protein sequence, x1 through DNAp */
65
66 for (i = max(0,-width-start_diag), x1+=i; i < len_prot; i++, x1++) {
67 bd = min(x1+width, len_dna_prot/3); /* upper bound of band */
68 bd1 = max(0,x1); /* lower bound of band */
69 wt = pam_matrix[prot_seq[i]];
70 del = 1-x1; /*adjustment*/
71 bd += del;
72 bd1 +=del;
73
74 ap = &cur[bd1];
75 aq = ap+1;
76 e1 = cur[bd1-1].C3;
77 e2 = ap->C1;
78 cd1 = cd2= cd3= 0;
79
80 for (dp = &dna_prot_seq[(bd1-del)*3]; ap < &cur[bd]; ap++) {
81 sc = max(max(e1, (e3=ap->C2))-gshift, e2)+wt[*dp++];
82 if (cd1 > sc) sc = cd1;
83 cd1 -= gext;
84 if ((ci = aq->I1) > 0) {
85 if (sc < ci) { ap->C1 = ci; ap->I1 = ci-gext;}
86 else {
87 ap->C1 = sc;
88 sc -= gg;
89 if (sc > 0) {
90 if (sc > best) best =sc;
91 if (cd1 < sc) cd1 = sc;
92 ap->I1 = max(ci-gext, sc);
93 } else ap->I1 = ci-gext;
94 }
95 } else {
96 if (sc <= 0) {
97 ap->I1 = ap->C1 = 0;
98 } else {
99 ap->C1 = sc; sc-=gg;
100 if (sc >0) {
101 if (sc > best) best =sc;
102 if (cd1 < sc) cd1 = sc;
103 ap->I1 = sc;
104 } else ap->I1 = 0;
105 }
106 }
107 sc = max(max(e2, (e1=ap->C3))-gshift, e3)+wt[*dp++];
108 if (cd2 > sc) sc = cd2;
109 cd2 -= gext;
110 if ((ci = aq->I2) > 0) {
111 if (sc < ci) { ap->C2 = ci; ap->I2 = ci-gext;}
112 else {
113 ap->C2 = sc;
114 sc -= gg;
115 if (sc > 0) {
116 if (sc > best) best =sc;
117 if (cd2 < sc) cd2 = sc;
118 ap->I2 = max(ci-gext, sc);
119 }
120 }
121 } else {
122 if (sc <= 0) {
123 ap->I2 = ap->C2 = 0;
124 } else {
125 ap->C2 = sc; sc-=gg;
126 if (sc >0) {
127 if (sc > best) best =sc;
128 if (cd2 < sc) cd2 = sc;
129 ap->I2 = sc;
130 } else ap->I2 = 0;
131 }
132 }
133 sc = max(max(e3, (e2=aq->C1))-gshift, e1)+wt[*dp++];
134 if (cd3 > sc) sc = cd3;
135 cd3 -= gext;
136 if ((ci = aq++->I3) > 0) {
137 if (sc < ci) { ap->C3 = ci; ap->I3 = ci-gext;}
138 else {
139 ap->C3 = sc;
140 sc -= gg;
141 if (sc > 0) {
142 if (sc > best) best =sc;
143 if (cd3 < sc) cd3 = sc;
144 ap->I3 = max(ci-gext, sc);
145 }
146 }
147 } else {
148 if (sc <= 0) {
149 ap->I3 = ap->C3 = 0;
150 } else {
151 ap->C3 = sc; sc-=gg;
152 if (sc >0) {
153 if (sc > best) best =sc;
154 if (cd3 < sc) cd3 = sc;
155 ap->I3 = sc;
156 } else ap->I3 = 0;
157 }
158 }
159 }
160 }
161 /* printf("The best score is %d\n", best); */
162 return best+gopen+gext;
163 }
164