1 | // ============================================================= |
---|
2 | /* */ |
---|
3 | // File : align.c |
---|
4 | // Purpose : |
---|
5 | /* */ |
---|
6 | // Institute of Microbiology (Technical University Munich) |
---|
7 | // http://www.arb-home.de/ |
---|
8 | /* */ |
---|
9 | // ============================================================= |
---|
10 | |
---|
11 | #include "mem.h" |
---|
12 | #include "trnsprob.h" |
---|
13 | |
---|
14 | #define EXTERN |
---|
15 | #include "i-hopper.h" |
---|
16 | |
---|
17 | #include <stdio.h> |
---|
18 | #include <math.h> |
---|
19 | |
---|
20 | #define MSIZE 512 |
---|
21 | #define ESIZE 36 |
---|
22 | |
---|
23 | //============================================================================ |
---|
24 | |
---|
25 | #define MINDIST 0.001 |
---|
26 | #define MAXDIST 1.000 |
---|
27 | |
---|
28 | #define MAXGAP 16 |
---|
29 | #define NOWHERE (MAXGAP+1) |
---|
30 | |
---|
31 | #define MINLENGTH 5 |
---|
32 | |
---|
33 | typedef struct score { |
---|
34 | double score; |
---|
35 | int up,down; |
---|
36 | } Score; |
---|
37 | |
---|
38 | |
---|
39 | typedef struct fragment { |
---|
40 | int beginX,beginY,length; |
---|
41 | struct fragment *next; |
---|
42 | } Fragment; |
---|
43 | |
---|
44 | typedef struct island { |
---|
45 | double score,upperScore,lowerScore; |
---|
46 | int beginX,beginY,endX,endY; |
---|
47 | int hasUpper,hasLower; |
---|
48 | struct island *next,*nextBeginIndex,*nextEndIndex,*nextSelected,*upper,*lower; |
---|
49 | struct fragment *fragments; |
---|
50 | } Island; |
---|
51 | |
---|
52 | static double **GP,**GS; |
---|
53 | |
---|
54 | // #define TEST |
---|
55 | |
---|
56 | #ifdef TEST |
---|
57 | |
---|
58 | #define TELEN 193 |
---|
59 | static double TE[TELEN][TELEN]; |
---|
60 | |
---|
61 | static int te=TRUE; |
---|
62 | |
---|
63 | #endif |
---|
64 | |
---|
65 | static Score **U; |
---|
66 | |
---|
67 | static double Thres,LogThres,Supp,GapA,GapB,GapC,expectedScore,**GE; |
---|
68 | |
---|
69 | static Island *Z,**ZB,**ZE; |
---|
70 | |
---|
71 | //============================================================================ |
---|
72 | |
---|
73 | static void initScore(double ***pP,double ***pS) { |
---|
74 | |
---|
75 | *pP=newDoubleMatrix(N,N); |
---|
76 | *pS=newDoubleMatrix(N1,N1); |
---|
77 | |
---|
78 | } |
---|
79 | |
---|
80 | //............................................................................ |
---|
81 | |
---|
82 | static void uninitScore(double ***pP,double ***pS) { |
---|
83 | |
---|
84 | freeMatrix(pS); |
---|
85 | freeMatrix(pP); |
---|
86 | |
---|
87 | } |
---|
88 | |
---|
89 | //............................................................................ |
---|
90 | |
---|
91 | static void updateScore( |
---|
92 | double **P,double **S, |
---|
93 | double F[N],double X[6],double dist |
---|
94 | ) { |
---|
95 | int i,j; |
---|
96 | double ***SS,s,smin,smax; |
---|
97 | |
---|
98 | P[T][T]=F[T]*F[T]; P[T][C]=F[T]*F[C]; P[T][A]=F[T]*F[A]; P[T][G]=F[T]*F[G]; |
---|
99 | P[C][T]=F[C]*F[T]; P[C][C]=F[C]*F[C]; P[C][A]=F[C]*F[A]; P[C][G]=F[C]*F[G]; |
---|
100 | P[A][T]=F[A]*F[T]; P[A][C]=F[A]*F[C]; P[A][A]=F[A]*F[A]; P[A][G]=F[A]*F[G]; |
---|
101 | P[G][T]=F[G]*F[T]; P[G][C]=F[G]*F[C]; P[G][A]=F[G]*F[A]; P[G][G]=F[G]*F[G]; |
---|
102 | |
---|
103 | initTrnsprob(&SS); |
---|
104 | updateTrnsprob(SS,F,X,REV); |
---|
105 | getTrnsprob(S,SS,dist); |
---|
106 | uninitTrnsprob(&SS); |
---|
107 | S[T][T]/=F[T]; S[T][C]/=F[C]; S[T][A]/=F[A]; S[T][G]/=F[G]; |
---|
108 | S[C][T]/=F[T]; S[C][C]/=F[C]; S[C][A]/=F[A]; S[C][G]/=F[G]; |
---|
109 | S[A][T]/=F[T]; S[A][C]/=F[C]; S[A][A]/=F[A]; S[A][G]/=F[G]; |
---|
110 | S[G][T]/=F[T]; S[G][C]/=F[C]; S[G][A]/=F[A]; S[G][G]/=F[G]; |
---|
111 | |
---|
112 | for(i=0;i<N;i++) for(j=0;j<N;j++) S[i][j]=log(S[i][j]); |
---|
113 | |
---|
114 | S[T][N]=F[T]*S[T][T]+F[C]*S[T][C]+F[A]*S[T][A]+F[G]*S[T][G]; |
---|
115 | S[C][N]=F[T]*S[C][T]+F[C]*S[C][C]+F[A]*S[C][A]+F[G]*S[C][G]; |
---|
116 | S[A][N]=F[T]*S[A][T]+F[C]*S[A][C]+F[A]*S[A][A]+F[G]*S[A][G]; |
---|
117 | S[G][N]=F[T]*S[G][T]+F[C]*S[G][C]+F[A]*S[G][A]+F[G]*S[G][G]; |
---|
118 | S[N][T]=F[T]*S[T][T]+F[C]*S[C][T]+F[A]*S[A][T]+F[G]*S[G][T]; |
---|
119 | S[N][C]=F[T]*S[T][C]+F[C]*S[C][C]+F[A]*S[A][C]+F[G]*S[G][C]; |
---|
120 | S[N][A]=F[T]*S[T][A]+F[C]*S[C][A]+F[A]*S[A][A]+F[G]*S[G][A]; |
---|
121 | S[N][G]=F[T]*S[T][G]+F[C]*S[C][G]+F[A]*S[A][G]+F[G]*S[G][G]; |
---|
122 | S[N][N]=F[T]*S[T][N]+F[C]*S[C][N]+F[A]*S[A][N]+F[G]*S[G][N]; |
---|
123 | |
---|
124 | smin=0.; smax=0.; |
---|
125 | for(i=0;i<N;i++) |
---|
126 | for(j=0;j<N;j++) { |
---|
127 | s=S[i][j]; |
---|
128 | if(s<smin) smin=s; else |
---|
129 | if(s>smax) smax=s; |
---|
130 | } |
---|
131 | |
---|
132 | } |
---|
133 | |
---|
134 | //============================================================================ |
---|
135 | |
---|
136 | static void initEntropy(double ***EE) { |
---|
137 | *EE=newDoubleMatrix(ESIZE+1,MSIZE); |
---|
138 | } |
---|
139 | |
---|
140 | //............................................................................ |
---|
141 | |
---|
142 | static void uninitEntropy(double ***EE) { |
---|
143 | freeMatrix(EE); |
---|
144 | } |
---|
145 | |
---|
146 | //............................................................................ |
---|
147 | |
---|
148 | static double getEntropy(double **E,double m,int l) { |
---|
149 | int k,ia,ib,ja,jb; |
---|
150 | double *M,mmin,idm,la,lb,ma,mb; |
---|
151 | |
---|
152 | M=E[0]; ma=0.5*(M[1]-M[0]); mmin=M[0]-ma; idm=0.5/ma; |
---|
153 | |
---|
154 | k=floor((m-mmin)*idm); |
---|
155 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
156 | |
---|
157 | if(k==0) {ja=k; jb=k+1;} else |
---|
158 | if(k==MSIZE-1) {ja=k-1; jb=k; } else |
---|
159 | if(m>M[k]) {ja=k; jb=k+1;} else |
---|
160 | {ja=k-1; jb=k; } |
---|
161 | |
---|
162 | ma=M[ja]; mb=M[jb]; |
---|
163 | |
---|
164 | if(l<=16) { |
---|
165 | return( |
---|
166 | ((mb-m)*E[l][ja]+(m-ma)*E[l][jb])/(mb-ma) |
---|
167 | ); |
---|
168 | } |
---|
169 | else { |
---|
170 | if(l<=32) { |
---|
171 | if(l<=18) {la=16; lb=18; ia=16; ib=17;} else |
---|
172 | if(l<=20) {la=18; lb=20; ia=17; ib=18;} else |
---|
173 | if(l<=22) {la=20; lb=22; ia=18; ib=19;} else |
---|
174 | if(l<=24) {la=22; lb=24; ia=19; ib=20;} else |
---|
175 | if(l<=26) {la=24; lb=26; ia=20; ib=21;} else |
---|
176 | if(l<=28) {la=26; lb=28; ia=21; ib=22;} else |
---|
177 | if(l<=30) {la=28; lb=30; ia=22; ib=23;} else |
---|
178 | if(l<=32) {la=30; lb=32; ia=23; ib=24;} |
---|
179 | } else |
---|
180 | if(l<=64) { |
---|
181 | if(l<=36) {la=32; lb=36; ia=24; ib=25;} else |
---|
182 | if(l<=40) {la=36; lb=40; ia=25; ib=26;} else |
---|
183 | if(l<=44) {la=40; lb=44; ia=26; ib=27;} else |
---|
184 | if(l<=48) {la=44; lb=48; ia=27; ib=28;} else |
---|
185 | if(l<=52) {la=48; lb=52; ia=28; ib=29;} else |
---|
186 | if(l<=56) {la=52; lb=56; ia=29; ib=30;} else |
---|
187 | if(l<=60) {la=56; lb=60; ia=30; ib=31;} else |
---|
188 | if(l<=64) {la=60; lb=64; ia=31; ib=32;} |
---|
189 | } |
---|
190 | else { |
---|
191 | if(l<= 128) {la= 64; lb= 128; ia=32; ib=33;} else |
---|
192 | if(l<= 256) {la= 128; lb= 256; ia=33; ib=34;} else |
---|
193 | if(l<= 512) {la= 256; lb= 512; ia=34; ib=35;} else |
---|
194 | {la= 512; lb=1024; ia=35; ib=36;} |
---|
195 | } |
---|
196 | return( |
---|
197 | ( |
---|
198 | +(lb-l)*((mb-m)*E[ia][ja]+(m-ma)*E[ia][jb]) |
---|
199 | +(l-la)*((mb-m)*E[ib][ja]+(m-ma)*E[ib][jb]) |
---|
200 | ) |
---|
201 | /((lb-la)*(mb-ma)) |
---|
202 | ); |
---|
203 | } |
---|
204 | |
---|
205 | } |
---|
206 | |
---|
207 | //............................................................................ |
---|
208 | |
---|
209 | static void updateEntropy(double **P,double **S,double **E) { |
---|
210 | int i,j,k,l,ll,lll; |
---|
211 | double *M,m,dm,idm,mmin,mmax; |
---|
212 | |
---|
213 | mmin=0.; mmax=0.; |
---|
214 | for(i=0;i<N;i++) |
---|
215 | for(j=0;j<N;j++) {m=S[i][j]; if(m<mmin) mmin=m; else if(m>mmax) mmax=m;} |
---|
216 | dm=(mmax-mmin)/MSIZE; |
---|
217 | idm=1./dm; |
---|
218 | |
---|
219 | M=E[0]; |
---|
220 | for(i=0,m=mmin+0.5*dm;i<MSIZE;i++,m+=dm) M[i]=m; |
---|
221 | |
---|
222 | for(i=0;i<MSIZE;i++) E[1][i]=0.; |
---|
223 | for(i=0;i<N;i++) |
---|
224 | for(j=0;j<N;j++) { |
---|
225 | k=floor((S[i][j]-mmin)*idm); |
---|
226 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
227 | E[1][k]+=P[i][j]; |
---|
228 | } |
---|
229 | |
---|
230 | for(k=0;k<MSIZE;k++) E[2][k]=0.; |
---|
231 | for(i=0;i<MSIZE;i++) |
---|
232 | for(j=0;j<MSIZE;j++) { |
---|
233 | m=(M[i]+M[j])/2; |
---|
234 | k=floor((m-mmin)*idm); |
---|
235 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
236 | E[2][k]+=E[1][i]*E[1][j]; |
---|
237 | } |
---|
238 | |
---|
239 | for(ll=2,lll=1;ll<=8;ll++,lll++) { |
---|
240 | l=ll+lll; |
---|
241 | for(k=0;k<MSIZE;k++) E[l][k]=0.; |
---|
242 | for(i=0;i<MSIZE;i++) |
---|
243 | for(j=0;j<MSIZE;j++) { |
---|
244 | m=(ll*M[i]+lll*M[j])/l; |
---|
245 | k=floor((m-mmin)*idm); |
---|
246 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
247 | E[l][k]+=E[ll][i]*E[lll][j]; |
---|
248 | } |
---|
249 | l=ll+ll; |
---|
250 | for(k=0;k<MSIZE;k++) E[l][k]=0.; |
---|
251 | for(i=0;i<MSIZE;i++) |
---|
252 | for(j=0;j<MSIZE;j++) { |
---|
253 | m=(M[i]+M[j])/2; |
---|
254 | k=floor((m-mmin)*idm); |
---|
255 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
256 | E[l][k]+=E[ll][i]*E[ll][j]; |
---|
257 | } |
---|
258 | } |
---|
259 | |
---|
260 | for(l=17,ll=l-8;l<=24;l++,ll++) { |
---|
261 | for(k=0;k<MSIZE;k++) E[l][k]=0.; |
---|
262 | for(i=0;i<MSIZE;i++) |
---|
263 | for(j=0;j<MSIZE;j++) { |
---|
264 | m=(M[i]+M[j])/2; |
---|
265 | k=floor((m-mmin)*idm); |
---|
266 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
267 | E[l][k]+=E[ll][i]*E[ll][j]; |
---|
268 | } |
---|
269 | } |
---|
270 | |
---|
271 | for(l=25,ll=l-8;l<=32;l++,ll++) { |
---|
272 | for(k=0;k<MSIZE;k++) E[l][k]=0.; |
---|
273 | for(i=0;i<MSIZE;i++) |
---|
274 | for(j=0;j<MSIZE;j++) { |
---|
275 | m=(M[i]+M[j])/2; |
---|
276 | k=floor((m-mmin)*idm); |
---|
277 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
278 | E[l][k]+=E[ll][i]*E[ll][j]; |
---|
279 | } |
---|
280 | } |
---|
281 | |
---|
282 | for(l=33,ll=l-1;l<=36;l++,ll++) { |
---|
283 | for(k=0;k<MSIZE;k++) E[l][k]=0.; |
---|
284 | for(i=0;i<MSIZE;i++) |
---|
285 | for(j=0;j<MSIZE;j++) { |
---|
286 | m=(M[i]+M[j])/2; |
---|
287 | k=floor((m-mmin)*idm); |
---|
288 | if(k>=MSIZE) k=MSIZE-1; else if(k<0) k=0; |
---|
289 | E[l][k]+=E[ll][i]*E[ll][j]; |
---|
290 | } |
---|
291 | } |
---|
292 | |
---|
293 | for(l=1;l<=ESIZE;l++) {m=0.; for(k=MSIZE-1;k>=0;k--) m=E[l][k]+=m;} |
---|
294 | |
---|
295 | for(i=1;i<=ESIZE;i++) |
---|
296 | for(j=0;j<MSIZE;j++) E[i][j]=(E[i][j]>REAL_MIN)?log(E[i][j]):LOG_REAL_MIN; |
---|
297 | |
---|
298 | { |
---|
299 | FILE *fp; |
---|
300 | int nbp; |
---|
301 | double m0,dm2; |
---|
302 | |
---|
303 | nbp=22; |
---|
304 | fp=fopen("distrib.txt","wt"); |
---|
305 | fprintf(fp,"\n(* score per basepair distribution for gap-less random path 1=Smin 100=Smax *)\nListPlot3D[({"); |
---|
306 | m0=M[0]; dm2=(M[MSIZE-1]-M[0])/99.; |
---|
307 | for(i=1;i<=nbp;i++) { |
---|
308 | fprintf(fp,"\n{%f",exp(getEntropy(E,m0,i))); |
---|
309 | for(j=1;j<=99;j++) fprintf(fp,",%f",exp(getEntropy(E,m0+dm2*j,i))); |
---|
310 | fprintf(fp,"}"); if(i<nbp) fprintf(fp,","); |
---|
311 | } |
---|
312 | fprintf(fp,"}),PlotRange->{{1,100},{1,%d},{0,1}},ViewPoint->{1.5,1.1,0.8}]\n",nbp); |
---|
313 | fclose(fp); |
---|
314 | |
---|
315 | } |
---|
316 | |
---|
317 | } |
---|
318 | |
---|
319 | //============================================================================ |
---|
320 | |
---|
321 | static Island *newIsland(char *X,char *Y,int i,int j,int d) { |
---|
322 | Island *p; Fragment *f,**ff,*L; int k,ii,jj,iii,jjj,l; double s; |
---|
323 | |
---|
324 | p=(Island*)newBlock(sizeof(Island)); |
---|
325 | |
---|
326 | ii=i; jj=j; l=0; s=0.; |
---|
327 | |
---|
328 | if(d>0) { |
---|
329 | ff=&L; |
---|
330 | for(;;) { |
---|
331 | s+=GS[(int)X[ii]][(int)Y[jj]]; |
---|
332 | if(++l==1) {iii=ii; jjj=jj;} |
---|
333 | k=U[ii][jj].up; |
---|
334 | if(k) { |
---|
335 | f=(Fragment*)newBlock(sizeof(Fragment)); |
---|
336 | f->beginX=iii; f->beginY=jjj; f->length=l; |
---|
337 | l=0; *ff=f; ff=&f->next; |
---|
338 | } |
---|
339 | if(k==NOWHERE) break; |
---|
340 | ii++; jj++; if(k>=0) ii+=k; else jj-=k; |
---|
341 | } |
---|
342 | *ff=NULL; |
---|
343 | p->beginX= i; p->beginY= j; |
---|
344 | p->endX =ii; p->endY =jj; |
---|
345 | } |
---|
346 | else { |
---|
347 | L=NULL; |
---|
348 | for(;;) { |
---|
349 | s+=GS[(int)X[ii]][(int)Y[jj]]; |
---|
350 | if(++l==1) {iii=ii; jjj=jj;} |
---|
351 | k=U[ii][jj].down; |
---|
352 | if(k) { |
---|
353 | f=(Fragment*)newBlock(sizeof(Fragment)); |
---|
354 | f->beginX=iii-l+1; f->beginY=jjj-l+1; f->length=l; |
---|
355 | l=0; f->next=L; L=f; |
---|
356 | } |
---|
357 | if(k==NOWHERE) break; |
---|
358 | ii--; jj--; if(k>=0) ii-=k; else jj+=k; |
---|
359 | } |
---|
360 | p->beginX=ii; p->beginY=jj; |
---|
361 | p->endX = i; p->endY = j; |
---|
362 | } |
---|
363 | |
---|
364 | p->fragments=L; |
---|
365 | |
---|
366 | p->score=s+GapC*expectedScore; |
---|
367 | |
---|
368 | return(p); |
---|
369 | } |
---|
370 | |
---|
371 | //............................................................................ |
---|
372 | |
---|
373 | static void freeIsland(Island **pp) { |
---|
374 | Island *p; Fragment *q; |
---|
375 | |
---|
376 | p=*pp; |
---|
377 | |
---|
378 | while(p->fragments) {q=p->fragments; p->fragments=q->next; freeBlock(&q);} |
---|
379 | |
---|
380 | freeBlock(pp); |
---|
381 | |
---|
382 | } |
---|
383 | |
---|
384 | //............................................................................ |
---|
385 | |
---|
386 | static void registerIsland(Island *f) { |
---|
387 | Island **pli; |
---|
388 | |
---|
389 | f->next=Z; Z=f; |
---|
390 | |
---|
391 | pli=&ZB[f->beginX+f->beginY]; |
---|
392 | f->nextBeginIndex=*pli; *pli=f; |
---|
393 | |
---|
394 | pli=&ZE[f->endX+f->endY]; |
---|
395 | f->nextEndIndex=*pli; *pli=f; |
---|
396 | |
---|
397 | } |
---|
398 | |
---|
399 | //............................................................................ |
---|
400 | |
---|
401 | static Island *selectUpperIslands(Island *f,int nX,int nY,int *incomplete) { |
---|
402 | int i,j; Island *l,*g; |
---|
403 | |
---|
404 | l=NULL; |
---|
405 | |
---|
406 | for(i=f->endX+f->endY+2,j=nX+nY-2;i<=j;i++) |
---|
407 | for(g=ZB[i];g;g=g->nextBeginIndex) |
---|
408 | if(g->beginX>f->endX&&g->beginY>f->endY) { |
---|
409 | if(!g->hasUpper) {*incomplete=TRUE; return(NULL);} |
---|
410 | g->nextSelected=l; l=g; |
---|
411 | } |
---|
412 | |
---|
413 | *incomplete=FALSE; |
---|
414 | |
---|
415 | return(l); |
---|
416 | } |
---|
417 | |
---|
418 | //............................................................................ |
---|
419 | |
---|
420 | static Island *selectLowerIslands(Island *f,int *incomplete) { |
---|
421 | int i; Island *l,*g; |
---|
422 | |
---|
423 | l=NULL; |
---|
424 | |
---|
425 | for(i=f->beginX+f->beginY-2;i>=0;i--) |
---|
426 | for(g=ZE[i];g;g=g->nextEndIndex) |
---|
427 | if(g->endX<f->beginX&&g->endY<f->beginY) { |
---|
428 | if(!g->hasLower) {*incomplete=TRUE; return(NULL);} |
---|
429 | g->nextSelected=l; l=g; |
---|
430 | } |
---|
431 | |
---|
432 | *incomplete=FALSE; |
---|
433 | |
---|
434 | return(l); |
---|
435 | } |
---|
436 | |
---|
437 | //............................................................................ |
---|
438 | |
---|
439 | static int areEqual(Island *a,Island *b) { |
---|
440 | Fragment *fa,*fb; |
---|
441 | |
---|
442 | if( a->beginX != b->beginX ) return(FALSE); |
---|
443 | if( a->beginY != b->beginY ) return(FALSE); |
---|
444 | if( a->endX != b->endX ) return(FALSE); |
---|
445 | if( a->endY != b->endY ) return(FALSE); |
---|
446 | |
---|
447 | fa=a->fragments; fb=b->fragments; |
---|
448 | while(fa&&fb) { |
---|
449 | if( fa->beginX != fb->beginX ) return(FALSE); |
---|
450 | if( fa->beginY != fb->beginY ) return(FALSE); |
---|
451 | if( fa->length != fb->length ) return(FALSE); |
---|
452 | fa=fa->next; fb=fb->next; |
---|
453 | } |
---|
454 | if(fa||fb) return(FALSE); |
---|
455 | |
---|
456 | return(TRUE); |
---|
457 | } |
---|
458 | |
---|
459 | //............................................................................ |
---|
460 | |
---|
461 | static int isUnique(Island *f) { |
---|
462 | Island *v; |
---|
463 | |
---|
464 | for(v=ZB[f->beginX+f->beginY];v;v=v->nextBeginIndex) |
---|
465 | if(areEqual(f,v)) return(FALSE); |
---|
466 | |
---|
467 | return(TRUE); |
---|
468 | } |
---|
469 | |
---|
470 | //............................................................................ |
---|
471 | |
---|
472 | static int isSignificant(Island *f) { |
---|
473 | int l; |
---|
474 | Fragment *q; |
---|
475 | |
---|
476 | l=0; for(q=f->fragments;q;q=q->next) l+=q->length; |
---|
477 | |
---|
478 | if(l<MINLENGTH) return(FALSE); |
---|
479 | |
---|
480 | if(getEntropy(GE,f->score/l,l)<=LogThres) return(TRUE); |
---|
481 | |
---|
482 | return(FALSE); |
---|
483 | } |
---|
484 | |
---|
485 | //............................................................................ |
---|
486 | |
---|
487 | static int I,J,K; |
---|
488 | |
---|
489 | //.... |
---|
490 | |
---|
491 | static void drawLowerPath(Island *f,int nX,char *X,char *XX,int nY,char *Y,char *YY) { |
---|
492 | int k; |
---|
493 | Fragment *q; |
---|
494 | |
---|
495 | if(f->lower) drawLowerPath(f->lower,nX,X,XX,nY,Y,YY); |
---|
496 | |
---|
497 | for(q=f->fragments;q;q=q->next) { |
---|
498 | while(I<q->beginX) {XX[K]=decodeBase(X[I++]); YY[K]='-'; K++;} |
---|
499 | while(J<q->beginY) {XX[K]='-'; YY[K]=decodeBase(Y[J++]); K++;} |
---|
500 | for(k=0;k<q->length;k++) |
---|
501 | {XX[K]=decodeBase(X[I++]); YY[K]=decodeBase(Y[J++]); K++;} |
---|
502 | } |
---|
503 | |
---|
504 | } |
---|
505 | |
---|
506 | //.... |
---|
507 | |
---|
508 | static void drawPath(Island *f,int nX,char *X,char *XX,int nY,char *Y,char *YY) { |
---|
509 | int k; |
---|
510 | Island *p; |
---|
511 | Fragment *q; |
---|
512 | |
---|
513 | I=0; J=0; K=0; |
---|
514 | |
---|
515 | if(f->lower) drawLowerPath(f->lower,nX,X,XX,nY,Y,YY); |
---|
516 | |
---|
517 | for(p=f;p;p=p->upper) |
---|
518 | for(q=p->fragments;q;q=q->next) { |
---|
519 | while(I<q->beginX) {XX[K]=decodeBase(X[I++]); YY[K]='-'; K++;} |
---|
520 | while(J<q->beginY) {XX[K]='-'; YY[K]=decodeBase(Y[J++]); K++;} |
---|
521 | for(k=0;k<q->length;k++) |
---|
522 | {XX[K]=decodeBase(X[I++]); YY[K]=decodeBase(Y[J++]); K++;} |
---|
523 | } |
---|
524 | |
---|
525 | while(I<nX) {XX[K]=decodeBase(X[I++]); YY[K]='-'; K++;} |
---|
526 | while(J<nY) {XX[K]='-'; YY[K]=decodeBase(Y[J++]); K++;} |
---|
527 | |
---|
528 | XX[K]='\0'; |
---|
529 | YY[K]='\0'; |
---|
530 | |
---|
531 | } |
---|
532 | |
---|
533 | //............................................................................ |
---|
534 | |
---|
535 | static void drawIsland(Island *f) { |
---|
536 | int i,j,k; |
---|
537 | Fragment *q; |
---|
538 | |
---|
539 | #ifdef TEST |
---|
540 | double score; |
---|
541 | score=f->score; |
---|
542 | #endif |
---|
543 | |
---|
544 | for(q=f->fragments;q;q=q->next) { |
---|
545 | for(i=q->beginX,j=q->beginY,k=0;k<q->length;k++,i++,j++) { |
---|
546 | |
---|
547 | #ifdef TEST |
---|
548 | if( |
---|
549 | i>=0&&i<TELEN&&j>=0&&j<TELEN // &&score>TE[i][j] |
---|
550 | ) { |
---|
551 | TE[i][j]=score; |
---|
552 | } |
---|
553 | #endif |
---|
554 | |
---|
555 | } |
---|
556 | } |
---|
557 | |
---|
558 | } |
---|
559 | |
---|
560 | //============================================================================ |
---|
561 | |
---|
562 | static double secS(int i,int j,char X[],int secX[],char Y[],int secY[]) { |
---|
563 | |
---|
564 | if(secX[i]||secY[j]) { |
---|
565 | if(secX[i]&&secY[j]) return(GS[(int)X[i]][(int)Y[j]]-Supp*expectedScore); |
---|
566 | return(-1.e34); |
---|
567 | } |
---|
568 | |
---|
569 | return(GS[(int)X[i]][(int)Y[j]]); |
---|
570 | } |
---|
571 | |
---|
572 | //============================================================================ |
---|
573 | |
---|
574 | static void AlignTwo( |
---|
575 | int nX,char X[],int secX[],char XX[],int nY,char Y[],int secY[],char YY[] |
---|
576 | ) { |
---|
577 | int r,changed,i,j,k,ii,jj,startx,stopx,starty,stopy; |
---|
578 | double s,ss; |
---|
579 | Island *z,*zz,*zzz,*best; |
---|
580 | |
---|
581 | //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO |
---|
582 | |
---|
583 | Z=NULL; for(i=nX+nY-2;i>=0;i--) {ZB[i]=NULL; ZE[i]=NULL;} |
---|
584 | |
---|
585 | startx=0; stopx=nX-1; |
---|
586 | for(i=startx,ii=i-1;i<=stopx;i++,ii++) { |
---|
587 | starty=0; stopy=nY-1; |
---|
588 | for(j=starty,jj=j-1;j<=stopy;j++,jj++) { |
---|
589 | s=0.; r=NOWHERE; |
---|
590 | if(i>startx&&j>starty) { |
---|
591 | ss=U[ii][jj].score; if(ss>s) {s=ss; r=0;} |
---|
592 | for(k=1;k<=MAXGAP&&ii-k>=0;k++) { |
---|
593 | ss=U[ii-k][jj].score+(GapA+k*GapB)*expectedScore; if(ss>s) {s=ss; r= k;} |
---|
594 | } |
---|
595 | for(k=1;k<=MAXGAP&&jj-k>=0;k++) { |
---|
596 | ss=U[ii][jj-k].score+(GapA+k*GapB)*expectedScore; if(ss>s) {s=ss; r=-k;} |
---|
597 | } |
---|
598 | } |
---|
599 | s+=secS(i,j,X,secX,Y,secY); if(s<0.) {s=0.; r=NOWHERE;} |
---|
600 | U[i][j].score=s; |
---|
601 | U[i][j].down=r; |
---|
602 | } |
---|
603 | } |
---|
604 | |
---|
605 | startx=0; stopx=nX-1; |
---|
606 | for(i=startx;i<=stopx;i++) { |
---|
607 | starty=0; stopy=nY-1; |
---|
608 | for(j=starty;j<=stopy;j++) { |
---|
609 | ii=i; jj=j; s=U[ii][jj].score; U[ii][jj].up=NOWHERE; |
---|
610 | for(;;) { |
---|
611 | k=U[ii][jj].down; if(k==NOWHERE) break; |
---|
612 | ii--; jj--; if(k>=0) ii-=k; else jj+=k; |
---|
613 | if(U[ii][jj].score>s) break; |
---|
614 | U[ii][jj].score=s; U[ii][jj].up=k; |
---|
615 | } |
---|
616 | } |
---|
617 | } |
---|
618 | |
---|
619 | startx=0; stopx=nX-1; |
---|
620 | for(i=startx;i<=stopx;i++) { |
---|
621 | starty=0; stopy=nY-1; |
---|
622 | for(j=starty;j<=stopy;j++) { |
---|
623 | if(U[i][j].down!=NOWHERE) continue; |
---|
624 | if(U[i][j].up==NOWHERE) continue; |
---|
625 | z=newIsland(X,Y,i,j,1); |
---|
626 | if(isUnique(z)&&isSignificant(z)) registerIsland(z); |
---|
627 | else freeIsland(&z); |
---|
628 | } |
---|
629 | } |
---|
630 | |
---|
631 | //---- |
---|
632 | |
---|
633 | startx=nX-1; stopx=0; |
---|
634 | for(i=startx,ii=i+1;i>=stopx;i--,ii--) { |
---|
635 | starty=nY-1; stopy=0; |
---|
636 | for(j=starty,jj=j+1;j>=stopy;j--,jj--) { |
---|
637 | s=0.; r=NOWHERE; |
---|
638 | if(i<startx&&j<starty) { |
---|
639 | ss=U[ii][jj].score; if(ss>s) {s=ss; r=0;} |
---|
640 | for(k=1;k<=MAXGAP&&ii+k<nX;k++) { |
---|
641 | ss=U[ii+k][jj].score+(GapA+k*GapB)*expectedScore; if(ss>s) {s=ss; r= k;} |
---|
642 | } |
---|
643 | for(k=1;k<=MAXGAP&&jj+k<nY;k++) { |
---|
644 | ss=U[ii][jj+k].score+(GapA+k*GapB)*expectedScore; if(ss>s) {s=ss; r=-k;} |
---|
645 | } |
---|
646 | } |
---|
647 | s+=secS(i,j,X,secX,Y,secY); if(s<0.) {s=0.; r=NOWHERE;} |
---|
648 | U[i][j].score=s; |
---|
649 | U[i][j].up=r; |
---|
650 | } |
---|
651 | } |
---|
652 | |
---|
653 | startx=nX-1; stopx=0; |
---|
654 | for(i=startx;i>=stopx;i--) { |
---|
655 | starty=nY-1; stopy=0; |
---|
656 | for(j=starty;j>=stopy;j--) { |
---|
657 | ii=i; jj=j; s=U[ii][jj].score; U[ii][jj].down=NOWHERE; |
---|
658 | for(;;) { |
---|
659 | k=U[ii][jj].up; if(k==NOWHERE) break; |
---|
660 | ii++; jj++; if(k>=0) ii+=k; else jj-=k; |
---|
661 | if(U[ii][jj].score>s) break; |
---|
662 | U[ii][jj].score=s; U[ii][jj].down=k; |
---|
663 | } |
---|
664 | } |
---|
665 | } |
---|
666 | |
---|
667 | startx=nX-1; stopx=0; |
---|
668 | for(i=startx;i>=stopx;i--) { |
---|
669 | starty=nY-1; stopy=0; |
---|
670 | for(j=starty;j>=stopy;j--) { |
---|
671 | if(U[i][j].up!=NOWHERE) continue; |
---|
672 | if(U[i][j].down==NOWHERE) continue; |
---|
673 | z=newIsland(X,Y,i,j,-1); |
---|
674 | if(isUnique(z)&&isSignificant(z)) registerIsland(z); |
---|
675 | else freeIsland(&z); |
---|
676 | } |
---|
677 | } |
---|
678 | |
---|
679 | //***************** |
---|
680 | |
---|
681 | for(z=Z;z;z=z->next) {z->hasUpper=FALSE; z->upper=NULL; z->upperScore=0.;} |
---|
682 | do { changed=FALSE; |
---|
683 | for(z=Z;z;z=z->next) { |
---|
684 | if(z->hasUpper) continue; |
---|
685 | zz=selectUpperIslands(z,nX,nY,&i); if(i) continue; |
---|
686 | if(zz) { |
---|
687 | s=0.; best=NULL; |
---|
688 | for(zzz=zz;zzz;zzz=zzz->nextSelected) { |
---|
689 | if(zzz->upperScore+zzz->score>s) {s=zzz->upperScore+zzz->score; best=zzz;} |
---|
690 | } |
---|
691 | if(best) {z->upper=best; z->upperScore=s;} |
---|
692 | } |
---|
693 | z->hasUpper=TRUE; |
---|
694 | changed=TRUE; |
---|
695 | } |
---|
696 | } while(changed); |
---|
697 | |
---|
698 | for(z=Z;z;z=z->next) {z->hasLower=FALSE; z->lower=NULL; z->lowerScore=0.;} |
---|
699 | do { changed=FALSE; |
---|
700 | for(z=Z;z;z=z->next) { |
---|
701 | if(z->hasLower) continue; |
---|
702 | zz=selectLowerIslands(z,&i); if(i) continue; |
---|
703 | if(zz) { |
---|
704 | s=0.; best=NULL; |
---|
705 | for(zzz=zz;zzz;zzz=zzz->nextSelected) { |
---|
706 | if(zzz->lowerScore+zzz->score>s) {s=zzz->lowerScore+zzz->score; best=zzz;} |
---|
707 | } |
---|
708 | if(best) {z->lower=best; z->lowerScore=s;} |
---|
709 | } |
---|
710 | z->hasLower=TRUE; |
---|
711 | changed=TRUE; |
---|
712 | } |
---|
713 | } while(changed); |
---|
714 | |
---|
715 | //***************** |
---|
716 | |
---|
717 | s=0.; best=NULL; |
---|
718 | for(z=Z;z;z=z->next) { |
---|
719 | ss=z->score+z->upperScore+z->lowerScore; |
---|
720 | if(ss>s) {best=z; s=ss;} |
---|
721 | } |
---|
722 | |
---|
723 | #ifdef TEST |
---|
724 | for(i=0;i<TELEN;i++) for(j=0;j<TELEN;j++) TE[i][j]=0.; |
---|
725 | #endif |
---|
726 | |
---|
727 | if(best) { drawPath(best,nX,X,XX,nY,Y,YY); |
---|
728 | drawIsland(best); |
---|
729 | for(z=best->upper;z;z=z->upper) drawIsland(z); |
---|
730 | for(z=best->lower;z;z=z->lower) drawIsland(z); |
---|
731 | } |
---|
732 | |
---|
733 | //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO |
---|
734 | |
---|
735 | #ifdef TEST |
---|
736 | |
---|
737 | if(te) {FILE *fp; te=FALSE; |
---|
738 | |
---|
739 | fp=fopen("subst.txt","wt"); |
---|
740 | fprintf(fp,"\n(* substitution matrix *)\nListDensityPlot[{"); |
---|
741 | for(i=0;i<N;i++) { |
---|
742 | fprintf(fp,"\n{%f",GS[i][0]); |
---|
743 | for(j=1;j<N;j++) fprintf(fp,",%f",GS[i][j]); |
---|
744 | fprintf(fp,"}"); if(i<N-1) fprintf(fp,","); |
---|
745 | } |
---|
746 | fprintf(fp,"}]\n"); |
---|
747 | fclose(fp); |
---|
748 | |
---|
749 | fp=fopen("correl.txt","w"); |
---|
750 | fprintf(fp,"\n(* correlation matrix *)\n{"); |
---|
751 | for(i=0;i<TELEN&&i<nX;i++) { |
---|
752 | fprintf(fp,"\n{"); |
---|
753 | fprintf(fp,"%f",secS(i,0,X,secX,Y,secY)); |
---|
754 | for(j=1;j<TELEN&&j<nY;j++) fprintf(fp,",%f",secS(i,j,X,secX,Y,secY)); |
---|
755 | if(i==TELEN-1||i==nX-1)fprintf(fp,"}"); else fprintf(fp,"},"); |
---|
756 | } |
---|
757 | fprintf(fp,"}//ListDensityPlot\n"); |
---|
758 | fclose(fp); |
---|
759 | |
---|
760 | fp=fopen("path.txt","w"); |
---|
761 | fprintf(fp,"\n(* pairwise alignment *)\n{"); |
---|
762 | for(i=0;i<TELEN&&i<nX;i++) { |
---|
763 | fprintf(fp,"\n{"); |
---|
764 | fprintf(fp,"%f\n",TE[i][0]); |
---|
765 | for(j=1;j<TELEN&&j<nY;j++) fprintf(fp,",%f",TE[i][j]); |
---|
766 | if(i==TELEN-1||i==nX-1)fprintf(fp,"}"); else fprintf(fp,"},"); |
---|
767 | } |
---|
768 | fprintf(fp,"}//ListDensityPlot\n"); |
---|
769 | fclose(fp); |
---|
770 | |
---|
771 | for(i=0;i<TELEN;i++) for(j=0;j<TELEN;j++) TE[i][j]=0.; |
---|
772 | for(z=Z;z;z=z->next) drawIsland(z); |
---|
773 | |
---|
774 | fp=fopen("islands.txt","w"); |
---|
775 | fprintf(fp,"\n(* smith-waterman-islands *)\n{"); |
---|
776 | for(i=0;i<TELEN&&i<nX;i++) { |
---|
777 | fprintf(fp,"\n{"); |
---|
778 | fprintf(fp,"%f\n",TE[i][0]); |
---|
779 | for(j=1;j<TELEN&&j<nY;j++) fprintf(fp,",%f",TE[i][j]); |
---|
780 | if(i==TELEN-1||i==nX-1)fprintf(fp,"}"); else fprintf(fp,"},"); |
---|
781 | } |
---|
782 | fprintf(fp,"}//ListDensityPlot\n"); |
---|
783 | fclose(fp); |
---|
784 | |
---|
785 | } |
---|
786 | |
---|
787 | #endif |
---|
788 | |
---|
789 | //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO |
---|
790 | |
---|
791 | while(Z) {z=Z; Z=Z->next; freeIsland(&z);} |
---|
792 | |
---|
793 | } |
---|
794 | |
---|
795 | //............................................................................ |
---|
796 | |
---|
797 | void Align( |
---|
798 | int nX,char X[],int secX[],char **XX,int nY,char Y[],int secY[],char **YY, |
---|
799 | int freqs,double fT,double fC,double fA,double fG, |
---|
800 | double rTC,double rTA,double rTG,double rCA,double rCG,double rAG, |
---|
801 | double dist,double supp,double gapA,double gapB,double gapC,double thres |
---|
802 | ) { |
---|
803 | double F[N],R[6]; |
---|
804 | char *s; |
---|
805 | int i,j,maxlen; |
---|
806 | |
---|
807 | *XX = (char*)newBlock((nX+nY+1)*sizeof(char)); |
---|
808 | *YY = (char*)newBlock((nX+nY+1)*sizeof(char)); |
---|
809 | |
---|
810 | Supp=supp; |
---|
811 | GapA=gapA; |
---|
812 | GapB=gapB; |
---|
813 | GapC=gapC; |
---|
814 | Thres=thres; |
---|
815 | |
---|
816 | if(dist>MAXDIST||dist<MINDIST) {Error="Bad argument"; return;} |
---|
817 | |
---|
818 | if(GapA<0.||GapB<0.) {Error="Bad argument"; return;} |
---|
819 | |
---|
820 | if(Thres>1.||Thres<=0.) {Error="Bad argument"; return;} |
---|
821 | LogThres=log(Thres); |
---|
822 | |
---|
823 | if(freqs) { |
---|
824 | if(fT<=0.||fC<=0.||fA<=0.||fG<=0.) {Error="Bad argument"; return;} |
---|
825 | } |
---|
826 | else { |
---|
827 | fT=0.; fC=0.; fA=0.; fG=0.; |
---|
828 | for(s=X;*s;s++) { |
---|
829 | switch(*s) { |
---|
830 | case 'T': fT++; break; case 'C': fC++; break; |
---|
831 | case 'A': fA++; break; case 'G': fG++; break; |
---|
832 | default: fT+=0.25; fC+=0.25; fA+=0.25; fG+=0.25; |
---|
833 | } |
---|
834 | } |
---|
835 | for(s=Y;*s;s++) { |
---|
836 | switch(*s) { |
---|
837 | case 'T': fT++; break; case 'C': fC++; break; |
---|
838 | case 'A': fA++; break; case 'G': fG++; break; |
---|
839 | default: fT+=0.25; fC+=0.25; fA+=0.25; fG+=0.25; |
---|
840 | } |
---|
841 | } |
---|
842 | } |
---|
843 | |
---|
844 | if(rTC<=0.||rTA<=0.||rTG<=0.||rCA<=0.||rCG<=0.||rAG<=0.) {Error="Bad argument"; return;} |
---|
845 | |
---|
846 | normalizeBaseFreqs(F,fT,fC,fA,fG); |
---|
847 | normalizeRateParams(R,rTC,rTA,rTG,rCA,rCG,rAG); |
---|
848 | |
---|
849 | maxlen=nX>nY?nX:nY; |
---|
850 | |
---|
851 | for(i=0;i<nX;i++) X[i]=encodeBase(X[i]); |
---|
852 | for(i=0;i<nY;i++) Y[i]=encodeBase(Y[i]); |
---|
853 | |
---|
854 | U=(Score **)newMatrix(maxlen,maxlen,sizeof(Score)); |
---|
855 | ZB=(Island **)newVector(2*maxlen,sizeof(Island *)); |
---|
856 | ZE=(Island **)newVector(2*maxlen,sizeof(Island *)); |
---|
857 | |
---|
858 | initScore(&GP,&GS); |
---|
859 | initEntropy(&GE); |
---|
860 | |
---|
861 | updateScore(GP,GS,F,R,dist); |
---|
862 | updateEntropy(GP,GS,GE); |
---|
863 | |
---|
864 | expectedScore=0.; |
---|
865 | for(i=0;i<N;i++) for(j=0;j<N;j++) expectedScore+=GP[i][j]*GS[i][j]; |
---|
866 | |
---|
867 | AlignTwo(nX,X,secX,*XX,nY,Y,secY,*YY); |
---|
868 | |
---|
869 | uninitEntropy(&GE); |
---|
870 | uninitScore(&GP,&GS); |
---|
871 | |
---|
872 | freeBlock(&ZE); |
---|
873 | freeBlock(&ZB); |
---|
874 | freeMatrix(&U); |
---|
875 | |
---|
876 | for(i=0;i<nX;i++) X[i]=decodeBase(X[i]); |
---|
877 | for(i=0;i<nY;i++) Y[i]=decodeBase(Y[i]); |
---|
878 | |
---|
879 | } |
---|