1 | #include <stdio.h> |
---|
2 | #include <string.h> |
---|
3 | #include <stdlib.h> |
---|
4 | #include "clustalv.h" |
---|
5 | |
---|
6 | /* |
---|
7 | * Prototypes |
---|
8 | */ |
---|
9 | |
---|
10 | extern Boolean read_tree(int *,double *,int *,int *,int *); |
---|
11 | extern void *ckalloc(size_t); |
---|
12 | |
---|
13 | void init_myers(void); |
---|
14 | void myers(int); |
---|
15 | |
---|
16 | static void group_gap(int,int,char *); |
---|
17 | static void add_ggaps(void); |
---|
18 | static void add(int); |
---|
19 | static int calc_weight(int,int,int,int); |
---|
20 | static void fill_pam(void); |
---|
21 | static int diff(int,int,int,int,int,int); |
---|
22 | static void do_align(int,int *); |
---|
23 | static int *alist; |
---|
24 | static int ** naa1,**naa2,**naas; |
---|
25 | |
---|
26 | /* |
---|
27 | * Global variables |
---|
28 | */ |
---|
29 | |
---|
30 | extern int nseqs, next, profile1_nseqs; /* DES */ |
---|
31 | extern int *seqlen_array; |
---|
32 | extern char **seq_array; |
---|
33 | extern int pamo[]; |
---|
34 | extern Boolean dnaflag,percent; |
---|
35 | extern int xover,big_pam; |
---|
36 | extern int nblocks; |
---|
37 | extern FILE *clustal_outfile,*tree; |
---|
38 | extern char treename[],seqname[],mtrxnam[],**names; |
---|
39 | extern char *matptr,idmat[],pam100mt[],pam250mt[]; |
---|
40 | extern int ktup,wind_gap,window; |
---|
41 | extern int signif; |
---|
42 | |
---|
43 | |
---|
44 | int gap_open, gap_extend; |
---|
45 | int dna_gap_open, dna_gap_extend; |
---|
46 | int prot_gap_open, prot_gap_extend; |
---|
47 | |
---|
48 | Boolean is_weight; |
---|
49 | |
---|
50 | extern int *zza,*zzb,*zzc,*zzd; /* allocated in show_pair.c */ |
---|
51 | static int *group; |
---|
52 | static int print_ptr,last_print,pos1,pos2; |
---|
53 | extern int * displ; /* allocated in show_pair.c */ |
---|
54 | static int pam[21][21]; |
---|
55 | int weights[21][21]; |
---|
56 | static int *fst_list,*snd_list; |
---|
57 | |
---|
58 | #define gap(x) ((x) <= 0 ? 0 : gap_open + gap_extend * (x)) |
---|
59 | |
---|
60 | void init_myers() |
---|
61 | { |
---|
62 | register int i; |
---|
63 | |
---|
64 | group = (int *)ckalloc( (MAXN+1) * sizeof (int)); |
---|
65 | |
---|
66 | alist = (int *)ckalloc( (MAXN+1) * sizeof (int)); |
---|
67 | fst_list = (int *)ckalloc( (MAXN+1) * sizeof (int)); |
---|
68 | snd_list = (int *)ckalloc( (MAXN+1) * sizeof (int)); |
---|
69 | |
---|
70 | naa1 = (int **)ckalloc(21 * sizeof (int *) ); |
---|
71 | for(i=0;i<21;i++) |
---|
72 | naa1[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int)); |
---|
73 | naa2 = (int **)ckalloc(21 * sizeof (int *) ); |
---|
74 | for(i=0;i<21;i++) |
---|
75 | naa2[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int)); |
---|
76 | naas = (int **)ckalloc(21 * sizeof (int *) ); |
---|
77 | for(i=0;i<21;i++) |
---|
78 | naas[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int)); |
---|
79 | |
---|
80 | dna_gap_open = 10; /* Default gap penalties for DNA */ |
---|
81 | dna_gap_extend = 10; |
---|
82 | |
---|
83 | prot_gap_open = 10; /* Default gap penalties for protein */ |
---|
84 | prot_gap_extend = 10; |
---|
85 | |
---|
86 | is_weight = TRUE; |
---|
87 | } |
---|
88 | |
---|
89 | |
---|
90 | |
---|
91 | static void group_gap(int len,int sclass,char *seq) |
---|
92 | { |
---|
93 | int i,j,k,xtra; |
---|
94 | |
---|
95 | for(i=1;i<=nseqs;++i) |
---|
96 | if(group[i] == sclass) { |
---|
97 | xtra = len - seqlen_array[i]; |
---|
98 | if(xtra>0) |
---|
99 | for(j=1;j<=xtra;++j) |
---|
100 | seq_array[i][seqlen_array[i]+j] = -1; |
---|
101 | for(j=1;j<=len;++j) |
---|
102 | if(seq[j] == '-') { |
---|
103 | for(k=len;k>=j+1;--k) |
---|
104 | seq_array[i][k] = seq_array[i][k-1]; |
---|
105 | seq_array[i][j] = -1; |
---|
106 | } |
---|
107 | seqlen_array[i] = len; |
---|
108 | } |
---|
109 | } |
---|
110 | |
---|
111 | |
---|
112 | static void add_ggaps() |
---|
113 | { |
---|
114 | int i,j,k,pos,to_do,len; |
---|
115 | char str1[MAXLEN+1],str2[MAXLEN+1]; |
---|
116 | |
---|
117 | pos=1; |
---|
118 | to_do=print_ptr-1; |
---|
119 | |
---|
120 | for(i=1;i<=to_do;++i) { |
---|
121 | if(displ[i]==0) { |
---|
122 | str1[pos]=str2[pos]='*'; |
---|
123 | ++pos; |
---|
124 | } |
---|
125 | else { |
---|
126 | if((k=displ[i])>0) { |
---|
127 | for(j=0;j<=k-1;++j) { |
---|
128 | str2[pos+j]='*'; |
---|
129 | str1[pos+j]='-'; |
---|
130 | } |
---|
131 | pos += k; |
---|
132 | } |
---|
133 | else { |
---|
134 | k = (displ[i]<0) ? displ[i] * -1 : displ[i]; |
---|
135 | for(j=0;j<=k-1;++j) { |
---|
136 | str1[pos+j]='*'; |
---|
137 | str2[pos+j]='-'; |
---|
138 | } |
---|
139 | pos += k; |
---|
140 | } |
---|
141 | } |
---|
142 | } |
---|
143 | |
---|
144 | len = --pos; |
---|
145 | group_gap(len,1,str1); |
---|
146 | group_gap(len,2,str2); |
---|
147 | } |
---|
148 | |
---|
149 | |
---|
150 | static void add(int v) |
---|
151 | { |
---|
152 | |
---|
153 | if(last_print<0) { |
---|
154 | displ[print_ptr-1] = v; |
---|
155 | displ[print_ptr++] = last_print; |
---|
156 | } |
---|
157 | else |
---|
158 | last_print = displ[print_ptr++] = v; |
---|
159 | } |
---|
160 | |
---|
161 | |
---|
162 | static int calc_weight(int iat,int jat,int v1,int v2) |
---|
163 | { |
---|
164 | int sum,i,j,lookn,ret; |
---|
165 | int ipos,jpos; |
---|
166 | |
---|
167 | sum=lookn=0; |
---|
168 | ipos = v1 + iat -1; |
---|
169 | jpos = v2 + jat -1; |
---|
170 | |
---|
171 | ret=0; |
---|
172 | if(pos1>=pos2) { |
---|
173 | for(i=1;i<=pos2;++i) { |
---|
174 | j=seq_array[alist[i]][jpos]; |
---|
175 | if(j>0) { |
---|
176 | sum += naas[j][ipos]; |
---|
177 | ++lookn; |
---|
178 | } |
---|
179 | } |
---|
180 | } |
---|
181 | else { |
---|
182 | for(i=1;i<=pos1;++i) { |
---|
183 | j=seq_array[alist[i]][ipos]; |
---|
184 | if(j>0) { |
---|
185 | sum += naas[j][jpos]; |
---|
186 | ++lookn; |
---|
187 | } |
---|
188 | } |
---|
189 | } |
---|
190 | |
---|
191 | if(sum > 0 ) ret = sum / lookn; |
---|
192 | return ret; |
---|
193 | } |
---|
194 | |
---|
195 | |
---|
196 | static void fill_pam() |
---|
197 | { |
---|
198 | register int i,j,pos; |
---|
199 | |
---|
200 | pos=0; |
---|
201 | |
---|
202 | for(i=0;i<20;++i) |
---|
203 | for(j=0;j<=i;++j) |
---|
204 | pam[i][j]=pamo[pos++]; |
---|
205 | |
---|
206 | for(i=0;i<20;++i) |
---|
207 | for(j=0;j<=i;++j) |
---|
208 | pam[j][i]=pam[i][j]; |
---|
209 | |
---|
210 | if(dnaflag) { |
---|
211 | xover=4; |
---|
212 | big_pam=8; |
---|
213 | for(i=0;i<5;++i) |
---|
214 | for(j=0;j<5;++j) { |
---|
215 | if(i==j) |
---|
216 | weights[i][j]=0; |
---|
217 | else |
---|
218 | weights[j][i]=10; |
---|
219 | } |
---|
220 | if(is_weight) { |
---|
221 | weights[1][3]=4; |
---|
222 | weights[3][1]=4; |
---|
223 | weights[2][4]=4; |
---|
224 | weights[4][2]=4; |
---|
225 | } |
---|
226 | } |
---|
227 | else { |
---|
228 | /* |
---|
229 | fprintf(stdout,"\nxover = %d; big_pam = %d\n",xover,big_pam); |
---|
230 | */ |
---|
231 | for(i=1;i<21;++i) |
---|
232 | for(j=1;j<21;++j) { |
---|
233 | weights[i][j] = big_pam - pam[i-1][j-1]; |
---|
234 | /* |
---|
235 | fprintf(stdout,"\n%2d vs %2d: %5d",i,j,weights[i][j]); |
---|
236 | */ |
---|
237 | } |
---|
238 | for(i=0;i<21;++i) { |
---|
239 | weights[0][i] = xover; |
---|
240 | weights[i][0] = xover; |
---|
241 | } |
---|
242 | } |
---|
243 | } |
---|
244 | |
---|
245 | static int diff(int v1,int v2,int v3,int v4,int st,int en) |
---|
246 | { |
---|
247 | int ctrc,ctri,ctrj,i,j,k,l,m,n,p,q,flag; |
---|
248 | |
---|
249 | q = gap_open + gap_extend; |
---|
250 | |
---|
251 | if(v4<=0) { |
---|
252 | if(v3>0) { |
---|
253 | if(last_print<0) |
---|
254 | last_print = displ[print_ptr-1] -= v3; |
---|
255 | else |
---|
256 | last_print = displ[print_ptr++] = -(v3); |
---|
257 | } |
---|
258 | |
---|
259 | return gap(v3); |
---|
260 | } |
---|
261 | |
---|
262 | if(v3<=1) { |
---|
263 | if(v3<=0) { |
---|
264 | add(v4); |
---|
265 | return gap(v4); |
---|
266 | } |
---|
267 | if(st>en) |
---|
268 | st=en; |
---|
269 | |
---|
270 | /***************if(!v4)*********BUG********************************/ |
---|
271 | |
---|
272 | ctrc = (st+gap_extend) + gap(v4); |
---|
273 | ctrj = 0; |
---|
274 | for(j=1;j<=v4;++j) { |
---|
275 | k = calc_weight(1,j,v1,v2) + gap(v4-j) + gap(j-1); |
---|
276 | if(k<ctrc) { |
---|
277 | ctrc = k; |
---|
278 | ctrj = j; |
---|
279 | } |
---|
280 | } |
---|
281 | |
---|
282 | if(!ctrj) { |
---|
283 | add(v4); |
---|
284 | if(last_print<0) |
---|
285 | last_print = displ[print_ptr-1] -= 1; |
---|
286 | else |
---|
287 | last_print = displ[print_ptr++] = -1; |
---|
288 | } |
---|
289 | else { |
---|
290 | if(ctrj>1) |
---|
291 | add(ctrj-1); |
---|
292 | displ[print_ptr++] = last_print = 0; |
---|
293 | if(ctrj<v4) |
---|
294 | add(v4-ctrj); |
---|
295 | } |
---|
296 | return ctrc; |
---|
297 | } |
---|
298 | |
---|
299 | |
---|
300 | ctri = v3 / 2; |
---|
301 | zza[0] = 0; |
---|
302 | p = gap_open; |
---|
303 | for(j=1;j<=v4;++j) { |
---|
304 | p += gap_extend; |
---|
305 | zza[j] = p; |
---|
306 | zzb[j] = p + gap_open; |
---|
307 | } |
---|
308 | |
---|
309 | p=st; |
---|
310 | for(i=1;i<=ctri;++i) { |
---|
311 | n=zza[0]; |
---|
312 | p += gap_extend; |
---|
313 | k = p; |
---|
314 | zza[0]=k; |
---|
315 | l = p+gap_open; |
---|
316 | for(j=1;j<=v4;++j) { |
---|
317 | k += q; |
---|
318 | l += gap_extend; |
---|
319 | if(k<l) |
---|
320 | l=k; |
---|
321 | k = zza[j] + q; |
---|
322 | m = zzb[j] + gap_extend; |
---|
323 | if(k<m) |
---|
324 | m=k; |
---|
325 | k = n + calc_weight(i,j,v1,v2); |
---|
326 | if(l<k) |
---|
327 | k=l; |
---|
328 | if(m<k) |
---|
329 | k=m; |
---|
330 | n=zza[j]; |
---|
331 | zza[j]=k; |
---|
332 | zzb[j]=m; |
---|
333 | } |
---|
334 | } |
---|
335 | |
---|
336 | zzb[0]=zza[0]; |
---|
337 | zzc[v4]=0; |
---|
338 | p=gap_open; |
---|
339 | for(j=v4-1;j>-1;--j) { |
---|
340 | p += gap_extend; |
---|
341 | zzc[j]=p; |
---|
342 | zzd[j]=p+gap_open; |
---|
343 | } |
---|
344 | p=en; |
---|
345 | for(i=v3-1;i>=ctri;--i) { |
---|
346 | n=zzc[v4]; |
---|
347 | p += gap_extend; |
---|
348 | k = p; |
---|
349 | zzc[v4] = k; |
---|
350 | l = p+gap_open; |
---|
351 | for(j=v4-1;j>=0;--j) { |
---|
352 | k += q; |
---|
353 | l += gap_extend; |
---|
354 | if(k<l) |
---|
355 | l=k; |
---|
356 | k = zzc[j] + q; |
---|
357 | m = zzd[j] + gap_extend; |
---|
358 | if(k<m) |
---|
359 | m=k; |
---|
360 | k = n + calc_weight(i+1,j+1,v1,v2); |
---|
361 | if(l<k) |
---|
362 | k=l; |
---|
363 | if(m<k) |
---|
364 | k=m; |
---|
365 | n=zzc[j]; |
---|
366 | zzc[j]=k; |
---|
367 | zzd[j]=m; |
---|
368 | } |
---|
369 | } |
---|
370 | |
---|
371 | zzd[v4]=zzc[v4]; |
---|
372 | ctrc=zza[0]+zzc[0]; |
---|
373 | ctrj=0; |
---|
374 | flag=1; |
---|
375 | for(j=0;j<=v4;++j) { |
---|
376 | k = zza[j] + zzc[j]; |
---|
377 | if(k<=ctrc) |
---|
378 | if(k<ctrc || ((zza[j]!=zzb[j]) && (zzc[j]==zzd[j]))) { |
---|
379 | ctrc=k; |
---|
380 | ctrj=j; |
---|
381 | } |
---|
382 | } |
---|
383 | |
---|
384 | for(j=v4;j>=0;--j) { |
---|
385 | k = zzb[j] + zzd[j] - gap_open; |
---|
386 | if(k<ctrc) { |
---|
387 | ctrc=k; |
---|
388 | ctrj=j; |
---|
389 | flag=2; |
---|
390 | } |
---|
391 | } |
---|
392 | |
---|
393 | /* Conquer recursively around midpoint */ |
---|
394 | |
---|
395 | if(flag==1) { /* Type 1 gaps */ |
---|
396 | diff(v1,v2,ctri,ctrj,st,gap_open); |
---|
397 | diff(v1+ctri,v2+ctrj,v3-ctri,v4-ctrj,gap_open,en); |
---|
398 | } |
---|
399 | else { |
---|
400 | diff(v1,v2,ctri-1,ctrj,st,0); |
---|
401 | if(last_print<0) /* Delete 2 */ |
---|
402 | last_print = displ[print_ptr-1] -= 2; |
---|
403 | else |
---|
404 | last_print = displ[print_ptr++] = -2; |
---|
405 | diff(v1+ctri+1,v2+ctrj,v3-ctri-1,v4-ctrj,0,en); |
---|
406 | } |
---|
407 | return ctrc; /* Return the score of the best alignment */ |
---|
408 | } |
---|
409 | |
---|
410 | |
---|
411 | |
---|
412 | |
---|
413 | |
---|
414 | |
---|
415 | static void do_align(int v1,int *score) |
---|
416 | { |
---|
417 | int i,j,k,l1,l2,n; |
---|
418 | int t_arr[21]; |
---|
419 | |
---|
420 | l1=l2=pos1=pos2=0; |
---|
421 | |
---|
422 | for(i=1;i<=MAXLEN;++i) { |
---|
423 | for(j=0;j<21;++j) |
---|
424 | naa1[j][i]=naa2[j][i]=0; |
---|
425 | for(j=1;j<21;++j) |
---|
426 | naas[j][i]=0; |
---|
427 | } |
---|
428 | |
---|
429 | for(i=1;i<=nseqs;++i) { |
---|
430 | if(group[i]==1) { |
---|
431 | fst_list[++pos1]=i; |
---|
432 | for(j=1;j<=seqlen_array[i];++j) |
---|
433 | if(seq_array[i][j]>0) { |
---|
434 | ++naa1[seq_array[i][j]][j]; |
---|
435 | ++naa1[0][j]; |
---|
436 | } |
---|
437 | if(seqlen_array[i]>l1) |
---|
438 | l1=seqlen_array[i]; |
---|
439 | } |
---|
440 | else if(group[i]==2) { |
---|
441 | snd_list[++pos2]=i; |
---|
442 | for(j=1;j<=seqlen_array[i];++j) |
---|
443 | if(seq_array[i][j]>0) { |
---|
444 | ++naa2[seq_array[i][j]][j]; |
---|
445 | ++naa2[0][j]; |
---|
446 | } |
---|
447 | if(seqlen_array[i]>l2) l2=seqlen_array[i]; |
---|
448 | } |
---|
449 | } |
---|
450 | |
---|
451 | if(pos1>=pos2) { |
---|
452 | for(i=1;i<=pos2;++i) |
---|
453 | alist[i]=snd_list[i]; |
---|
454 | for(n=1;n<=l1;++n) { |
---|
455 | for(i=1;i<21;++i) |
---|
456 | t_arr[i]=0; |
---|
457 | for(i=1;i<21;++i) |
---|
458 | if(naa1[i][n]>0) |
---|
459 | for(j=1;j<21;++j) |
---|
460 | t_arr[j] += (weights[i][j]*naa1[i][n]); |
---|
461 | k = naa1[0][n]; |
---|
462 | if(k>0) |
---|
463 | for(i=1;i<21;++i) |
---|
464 | naas[i][n]=t_arr[i]/k; |
---|
465 | } |
---|
466 | } |
---|
467 | else { |
---|
468 | for(i=1;i<=pos1;++i) |
---|
469 | alist[i]=fst_list[i]; |
---|
470 | for(n=1;n<=l2;++n) { |
---|
471 | for(i=1;i<21;++i) |
---|
472 | t_arr[i]=0; |
---|
473 | for(i=1;i<21;++i) |
---|
474 | if(naa2[i][n]>0) |
---|
475 | for(j=1;j<21;++j) |
---|
476 | t_arr[j] += (weights[i][j]*naa2[i][n]); |
---|
477 | k = naa2[0][n]; |
---|
478 | if(k>0) |
---|
479 | for(i=1;i<21;++i) |
---|
480 | naas[i][n]=t_arr[i]/k; |
---|
481 | } |
---|
482 | } |
---|
483 | |
---|
484 | *score=diff(1,1,l1,l2,v1,v1); /* Myers and Miller alignment now */ |
---|
485 | } |
---|
486 | |
---|
487 | |
---|
488 | void myers(int align_type) /* align_type = 0 for full progressive alignment*/ |
---|
489 | /* align_type = 1 for a profile alignment */ |
---|
490 | { |
---|
491 | /* static char *nbases = "XACGT"; |
---|
492 | static char seq1[MAXLEN+1]; */ |
---|
493 | char temp[MAXLINE]; |
---|
494 | int i,set; |
---|
495 | int sets,entries,idummy,score; |
---|
496 | double dummy; |
---|
497 | |
---|
498 | nblocks=0; |
---|
499 | |
---|
500 | fprintf(stdout,"\nStart of Multiple Alignment\n"); |
---|
501 | |
---|
502 | fill_pam(); |
---|
503 | |
---|
504 | if(align_type == 0) { /* a full progressive alignment */ |
---|
505 | sets=0; |
---|
506 | tree=fopen(treename,"r"); |
---|
507 | while(fgets(temp,MAXLINE,tree)!=NULL) ++sets; |
---|
508 | fseek(tree,0,0); |
---|
509 | fprintf(stdout,"There are %d sets\n",sets); |
---|
510 | } |
---|
511 | else /* just one set (a profile alignment) */ |
---|
512 | sets = 1; |
---|
513 | |
---|
514 | fprintf(stdout,"Aligning...\n"); |
---|
515 | |
---|
516 | for(set=1;set<=sets;++set) { |
---|
517 | if(align_type == 0) |
---|
518 | read_tree(&entries,&dummy,&idummy,&idummy,group); |
---|
519 | else { |
---|
520 | for(i=1; i<=profile1_nseqs; ++i) |
---|
521 | group[i] = 1; |
---|
522 | for(i=profile1_nseqs+1; i<=nseqs; ++i) |
---|
523 | group[i] = 2; |
---|
524 | entries = nseqs; |
---|
525 | } |
---|
526 | last_print=0; |
---|
527 | print_ptr=1; |
---|
528 | do_align(gap_open,&score); |
---|
529 | fprintf(stdout,"Set %d: Entries:%d Score:%d\n",set,entries,score); |
---|
530 | add_ggaps(); |
---|
531 | } |
---|
532 | if(align_type == 0) fclose(tree); |
---|
533 | |
---|
534 | /* make the rest (output stuff) into routine clustal_out in file amenu.c */ |
---|
535 | |
---|
536 | } |
---|
537 | |
---|
538 | |
---|