source: trunk/GDE/CLUSTAL/showpair.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.7 KB
Line 
1#include <stdio.h>
2#include <string.h>
3#include <stdlib.h>
4#include <math.h>
5#include "clustalv.h"   
6
7
8/*
9*       Prototypes
10*/
11
12extern void *ckalloc(size_t);
13extern void error(const char *,...);
14extern int *seqlen_array;
15extern char **seq_array;
16
17void init_show_pair(void);
18void show_pair(void);
19static void make_p_ptrs(int *,int *,int,int);
20static void make_n_ptrs(int *,int *,int,int);
21static void put_frag(int,int,int,int);
22static int frag_rel_pos(int,int,int,int);
23static void pair_align(int,int,int);
24static void des_quick_sort(int *, int *, int);
25
26
27/*
28*        Global variables
29*/
30
31extern int next,nseqs;
32extern Boolean dnaflag;
33extern double **tmat;
34
35int ktup,window,wind_gap,signif;                      /* Pairwise aln. params */
36int  dna_ktup, dna_window, dna_wind_gap, dna_signif;  /* params for DNA */
37int prot_ktup,prot_window,prot_wind_gap,prot_signif;  /* params for prots */
38
39
40Boolean percent;
41
42static int curr_frag,maxsf,vatend;
43static int **accum;
44int *displ;                                             /* also used in myers.c     */
45int *zza, *zzb, *zzc, *zzd;             /* also used in myers.c     */
46static int *diag_index;
47static char *slopes;
48
49void init_show_pair(void)
50{
51        register int i;
52
53        accum = (int **)ckalloc( 5*sizeof (int *) );
54        for (i=0;i<5;i++)
55                accum[i] = (int *) ckalloc(FSIZE * sizeof (int) );
56
57        displ      = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
58        slopes     = (char *)ckalloc( (2*MAXLEN +1) * sizeof (char));
59        diag_index = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
60
61        zza = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
62        zzb = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
63
64        zzc = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
65        zzd = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
66
67        dna_ktup      = 2;   /* default parameters for DNA */
68        dna_wind_gap  = 5;
69        dna_signif    = 4;
70        dna_window    = 4;
71
72    prot_ktup     = 1;   /* default parameters for proteins */
73        prot_wind_gap = 3;
74        prot_signif   = 5;
75        prot_window   = 5;
76
77        percent=TRUE;
78}
79
80static void make_p_ptrs(int *tptr,int *pl,int naseq,int l)
81{
82        static int a[]={ 0, 1, 20, 400 };
83        int i,j,limit,code,flag;
84        char residue;
85       
86        limit = (int) pow((double)20,(double)ktup);
87        for(i=1;i<=limit;++i)
88                pl[i]=0;
89        for(i=1;i<=l;++i)
90                tptr[i]=0;
91       
92        for(i=1;i<=(l-ktup+1);++i) {
93                code=0;
94                flag=FALSE;
95                for(j=1;j<=ktup;++j) {
96                        residue = seq_array[naseq][i+j-1];
97                        if(residue<=0) {
98                                flag=TRUE;
99                                break;
100                        }
101                        code += ((residue-1) * a[j]);
102                }
103                if(flag)
104                        continue;
105                ++code;
106                if(pl[code]!=0)
107                        tptr[i]=pl[code];
108                pl[code]=i;
109        }
110}
111
112
113static void make_n_ptrs(int *tptr,int *pl,int naseq,int len)
114{
115        static int pot[]={ 0, 1, 4, 16, 64, 256, 1024, 4096 };
116        int i,j,limit,code,flag;
117        char residue;
118       
119        limit = (int) pow((double)4,(double)ktup);
120       
121        for(i=1;i<=limit;++i)
122                pl[i]=0;
123        for(i=1;i<=len;++i)
124                tptr[i]=0;
125       
126        for(i=1;i<=len-ktup+1;++i) {
127                code=0;
128                flag=FALSE;
129                for(j=1;j<=ktup;++j) {
130                        residue = seq_array[naseq][i+j-1];
131                        if(residue<=0) {
132                                flag=TRUE;
133                                break;
134                        }
135                        code += ((residue-1) * pot[j]);
136                }
137                if(flag)
138                        continue;
139                ++code;
140                if(pl[code]!=0)
141                        tptr[i]=pl[code];
142                pl[code]=i;
143        }
144}
145
146
147static void put_frag(int fs,int v1,int v2,int flen)
148{
149        int end;
150       
151        accum[0][curr_frag]=fs;
152        accum[1][curr_frag]=v1;
153        accum[2][curr_frag]=v2;
154        accum[3][curr_frag]=flen;
155       
156        if(!maxsf) {
157                maxsf=1;
158                accum[4][curr_frag]=0;
159                return;
160        }
161       
162        if(fs >= accum[0][maxsf]) {
163                accum[4][curr_frag]=maxsf;
164                maxsf=curr_frag;
165                return;
166        }
167        else {
168                next=maxsf;
169                while(TRUE) {
170                        end=next;
171                        next=accum[4][next];
172                        if(fs>=accum[0][next])
173                                break;
174                }
175                accum[4][curr_frag]=next;
176                accum[4][end]=curr_frag;
177        }
178}
179
180
181static int frag_rel_pos(int a1,int b1,int a2,int b2)
182{
183        int ret;
184       
185        ret=FALSE;
186        if(a1-b1==a2-b2) {
187                if(a2<a1)
188                        ret=TRUE;
189        }
190        else {
191                if(a2+ktup-1<a1 && b2+ktup-1<b1)
192                        ret=TRUE;
193        }
194        return ret;
195}
196
197
198static void des_quick_sort(int *array1, int *array2, int array_size)
199/*  */
200/* Quicksort routine, adapted from chapter 4, page 115 of software tools */
201/* by Kernighan and Plauger, (1986) */
202/* Sort the elements of array1 and sort the */
203/* elements of array2 accordingly */
204/*  */
205{
206        int temp1, temp2;
207        int p, pivlin;
208        int i, j;
209        int lst[50], ust[50];       /* the maximum no. of elements must be*/
210                                                                /* < log(base2) of 50 */
211
212        lst[1] = 1;
213        ust[1] = array_size;
214        p = 1;
215
216        while(p > 0) {
217                if(lst[p] >= ust[p])
218                        p--;
219                else {
220                        i = lst[p] - 1;
221                        j = ust[p];
222                        pivlin = array1[j];
223                        while(i < j) {
224                                for(i=i+1; array1[i] < pivlin; i++)
225                                        ;
226                                for(j=j-1; j > i; j--)
227                                        if(array1[j] <= pivlin) break;
228                                if(i < j) {
229                                        temp1     = array1[i];
230                                        array1[i] = array1[j];
231                                        array1[j] = temp1;
232                                       
233                                        temp2     = array2[i];
234                                        array2[i] = array2[j];
235                                        array2[j] = temp2;
236                                }
237                        }
238                       
239                        j = ust[p];
240
241                        temp1     = array1[i];
242                        array1[i] = array1[j];
243                        array1[j] = temp1;
244
245                        temp2     = array2[i];
246                        array2[i] = array2[j];
247                        array2[j] = temp2;
248
249                        if(i-lst[p] < ust[p] - i) {
250                                lst[p+1] = lst[p];
251                                ust[p+1] = i - 1;
252                                lst[p]   = i + 1;
253                        }
254                        else {
255                                lst[p+1] = i + 1;
256                                ust[p+1] = ust[p];
257                                ust[p]   = i - 1;
258                        }
259                        p = p + 1;
260                }
261        }
262        return;
263
264}
265
266
267
268
269
270static void pair_align(int seq_no,int l1,int l2)
271{
272        int pot[8],i,j,l,m,flag,limit,pos,tl1,vn1,vn2,flen,osptr,fs;
273        int tv1,tv2,encrypt,subt1,subt2,rmndr;
274        char residue;
275       
276        if(dnaflag) {
277                for(i=1;i<=ktup;++i)
278                        pot[i] = (int) pow((double)4,(double)(i-1));
279                limit = (int) pow((double)4,(double)ktup);
280        }
281        else {
282                pot[1]=1;
283                pot[2]=20;
284                pot[3]=400;
285                limit = (int) pow(20.0,(double)ktup);
286        }
287       
288        tl1 = (l1+l2)-1;
289       
290        for(i=1;i<=tl1;++i) {
291                slopes[i]=displ[i]=0;
292                diag_index[i] = i;
293        }
294       
295
296/* increment diagonal score for each k_tuple match */
297
298        for(i=1;i<=limit;++i) {
299                vn1=zzc[i];
300                while(TRUE) {
301                        if(!vn1) break;
302                        vn2=zzd[i];
303                        while(vn2 != 0) {
304                                osptr=vn1-vn2+l2;
305                                ++displ[osptr];
306                                vn2=zzb[vn2];
307                        }
308                        vn1=zza[vn1];
309                }
310        }
311
312/* choose the top SIGNIF diagonals */
313
314        des_quick_sort(displ, diag_index, tl1);
315
316        j = tl1 - signif + 1;
317        if(j < 1) j = 1;
318 
319/* flag all diagonals within WINDOW of a top diagonal */
320
321        for(i=tl1; i>=j; i--) 
322                if(displ[i] > 0) {
323                        pos = diag_index[i];
324                        l = (1  >pos-window) ? 1   : pos-window;
325                        m = (tl1<pos+window) ? tl1 : pos+window;
326                        for(; l <= m; l++) 
327                                slopes[l] = 1;
328                }
329
330        for(i=1; i<=tl1; i++)  displ[i] = 0;
331
332       
333        curr_frag=maxsf=0;
334       
335        for(i=1;i<=(l1-ktup+1);++i) {
336                encrypt=flag=0;
337                for(j=1;j<=ktup;++j) {
338                        residue = seq_array[seq_no][i+j-1];
339                        if(residue<=0) {
340                                flag=TRUE;
341                                break;
342                        }
343                        encrypt += ((residue-1)*pot[j]);
344                }
345                if(flag) continue;
346                ++encrypt;
347       
348                vn2=zzd[encrypt];
349       
350                flag=FALSE;
351                while(TRUE) {
352                        if(!vn2) {
353                                flag=TRUE;
354                                break;
355                        }
356                        osptr=i-vn2+l2;
357                        if(slopes[osptr]!=1) {
358                                vn2=zzb[vn2];
359                                continue;
360                        }
361                        flen=0;
362                        fs=ktup;
363                        next=maxsf;
364               
365               
366                /*
367                * A-loop
368                */
369               
370                        while(TRUE) {
371                                if(!next) {
372                                        ++curr_frag;
373                                        if(curr_frag>=FSIZE) {
374                                                fprintf(stdout,"(Partial alignment)");
375                                                vatend=1;
376                                                return;
377                                        }
378                                        displ[osptr]=curr_frag;
379                                        put_frag(fs,i,vn2,flen);
380                                }
381                                else {
382                                        tv1=accum[1][next];
383                                        tv2=accum[2][next];
384                                        if(frag_rel_pos(i,vn2,tv1,tv2)) {
385                                                if(i-vn2==accum[1][next]-accum[2][next]) {
386                                                        if(i>accum[1][next]+(ktup-1))
387                                                                fs=accum[0][next]+ktup;
388                                                        else {
389                                                                rmndr=i-accum[1][next];
390                                                                fs=accum[0][next]+rmndr;
391                                                        }
392                                                        flen=next;
393                                                        next=0;
394                                                        continue;
395                                                }
396                                                else {
397                                                        if(displ[osptr]==0)
398                                                                subt1=ktup;
399                                                        else {
400                                                                if(i>accum[1][displ[osptr]]+(ktup-1))
401                                                                        subt1=accum[0][displ[osptr]]+ktup;
402                                                                else {
403                                                                        rmndr=i-accum[1][displ[osptr]];
404                                                                        subt1=accum[0][displ[osptr]]+rmndr;
405                                                                }
406                                                        }
407                                                        subt2=accum[0][next]-wind_gap+ktup;
408                                                        if(subt2>subt1) {
409                                                                flen=next;
410                                                                fs=subt2;
411                                                        }
412                                                        else {
413                                                                flen=displ[osptr];
414                                                                fs=subt1;
415                                                        }
416                                                        next=0;
417                                                        continue;
418                                                }
419                                        }
420                                        else {
421                                                next=accum[4][next];
422                                                continue;
423                                        }
424                                }
425                                break;
426                        }
427                /*
428                * End of Aloop
429                */
430               
431                        vn2=zzb[vn2];
432                }
433        }
434        vatend=0;
435}               
436
437void show_pair()
438{
439        int i,j,dsr;
440        double calc_score;
441       
442        fprintf(stdout,"\n\n");
443       
444        for(i=1;i<=nseqs;++i) {
445                if(dnaflag)
446                        make_n_ptrs(zza,zzc,i,seqlen_array[i]);
447                else
448                        make_p_ptrs(zza,zzc,i,seqlen_array[i]);
449                for(j=i+1;j<=nseqs;++j) {
450                        if(dnaflag)
451                                make_n_ptrs(zzb,zzd,j,seqlen_array[j]);
452                        else
453                                make_p_ptrs(zzb,zzd,j,seqlen_array[j]);
454                        pair_align(i,seqlen_array[i],seqlen_array[j]);
455                        if(!maxsf)
456                                calc_score=0.0;
457                        else {
458                                calc_score=(double)accum[0][maxsf];
459                                if(percent) {
460                                        dsr=(seqlen_array[i]<seqlen_array[j]) ?
461                                                        seqlen_array[i] : seqlen_array[j];
462                                        calc_score = (calc_score/(double)dsr) * 100.0;
463                                }
464                        }
465                        tmat[i][j]=calc_score;
466                        tmat[j][i]=calc_score;
467                        if(calc_score>0.1)
468                                fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %lg\n",
469               i,j,calc_score);
470                        else
471                                fprintf(stdout,"Sequences (%d:%d) Not Aligned\n",i,j);
472                }
473        }
474}
475
Note: See TracBrowser for help on using the repository browser.