source: tags/arb_5.1/ISLAND_HOPPING/align.c

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