source: branches/profile/ISLAND_HOPPING/align.cxx

Last change on this file was 8607, checked in by westram, 7 years ago

merge from e4fix [8135] [8136] [8137] [8138] [8139] [8140] [8141] [8142] [8143] [8144] [8222]
(this revives the reverted patches [8129] [8130] [8131] [8132]; see [8133])

  • fixes
    • some free/delete mismatches
    • wrong definition of ORF objects (Level was no bit value)
    • amino consensus (failed for columns only containing 'C')
  • rename
    • AA_sequence_term → orf_term
    • ED4_sequence_terminal_basic → ED4_abstract_sequence_terminal
  • cleaned up hierarchy dumps
  • tweaked is_terminal()/to_terminal()
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.5 KB
Line 
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
33typedef struct score {
34    double score;
35    int up,down;
36} Score;
37
38
39typedef struct fragment {
40    int beginX,beginY,length;
41    struct fragment *next;
42} Fragment;
43
44typedef 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
52static double **GP,**GS;
53
54// #define TEST
55
56#ifdef TEST
57
58#define TELEN 193
59static double TE[TELEN][TELEN];
60
61static int te=TRUE;
62
63#endif
64
65static Score **U;
66
67static double Thres,LogThres,Supp,GapA,GapB,GapC,expectedScore,**GE;
68
69static Island *Z,**ZB,**ZE;
70
71//============================================================================
72
73static void initScore(double ***pP,double ***pS) {
74
75    *pP=newDoubleMatrix(N,N);
76    *pS=newDoubleMatrix(N1,N1);
77
78}
79
80//............................................................................
81
82static void uninitScore(double ***pP,double ***pS) {
83
84    freeMatrix(pS);
85    freeMatrix(pP);
86
87}
88
89//............................................................................
90
91static 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
136static void initEntropy(double ***EE) {
137    *EE=newDoubleMatrix(ESIZE+1,MSIZE);
138}
139
140//............................................................................
141
142static void uninitEntropy(double ***EE) {
143    freeMatrix(EE);
144}
145
146//............................................................................
147
148static 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
209static 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
321static 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
373static 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
386static 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
401static 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
420static 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
439static 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
461static 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
472static 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
487static int I,J,K;
488
489//....
490
491static 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
508static 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
535static 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
562static 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
574static 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
797void 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}
Note: See TracBrowser for help on using the repository browser.