source: tags/initial/GDE/CLUSTALW/pairalign.c

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 KB
Line 
1/* Change int h to int gh everywhere  DES June 1994 */
2#include <stdio.h>
3#include <stdlib.h>
4#include <string.h>
5#include <math.h>
6#include "clustalw.h"
7
8#define MIN(a,b) ((a)<(b)?(a):(b))
9#define MAX(a,b) ((a)>(b)?(a):(b))
10
11#define gap(k)  ((k) <= 0 ? 0 : g + gh * (k))
12#define tbgap(k)  ((k) <= 0 ? 0 : tb + gh * (k))
13#define tegap(k)  ((k) <= 0 ? 0 : te + gh * (k))
14
15/*
16*       Prototypes
17*/
18static void add(sint v);
19static sint calc_score(sint iat, sint jat, sint v1, sint v2);
20static float tracepath(sint tsb1,sint tsb2);
21static void forward_pass(char *ia, char *ib, sint n, sint m);
22static void reverse_pass(char *ia, char *ib);
23static sint diff(sint A, sint B, sint M, sint N, sint tb, sint te);
24static void del(sint k);
25
26/*
27 *   Global variables
28 */
29#ifdef MAC
30#define pwint   short
31#else
32#define pwint   int
33#endif
34static sint             int_scale;
35
36extern double   **tmat;
37extern float    pw_go_penalty;
38extern float    pw_ge_penalty;
39extern float    transition_weight;
40extern sint     nseqs;
41extern sint     max_aa;
42extern sint     gap_pos1,gap_pos2;
43extern sint     max_aln_length;
44extern sint     *seqlen_array;
45extern sint     debug;
46extern sint     mat_avscore;
47extern short    blosum30mt[],pam350mt[],idmat[],pw_usermat[],pw_userdnamat[];
48extern short    clustalvdnamt[],swgapdnamt[];
49extern short    gon250mt[];
50extern short    def_dna_xref[],def_aa_xref[],pw_dna_xref[],pw_aa_xref[];
51extern Boolean  dnaflag;
52extern char     **seq_array;
53extern char     *amino_acid_codes;
54extern char     pw_mtrxname[];
55extern char     pw_dnamtrxname[];
56
57static float    mm_score;
58static sint     print_ptr,last_print;
59static sint     *displ;
60static pwint    *HH, *DD, *RR, *SS;
61static sint     g, gh;
62static sint     seq1, seq2;
63static sint     matrix[NUMRES][NUMRES];
64static pwint    maxscore;
65static sint     sb1, sb2, se1, se2;
66
67
68sint pairalign(sint istart, sint iend, sint jstart, sint jend)
69{
70  short  *mat_xref;
71  static sint    si, sj, i;
72  static sint    n,m,len1,len2;
73  static sint    maxres;
74  static short    *matptr;
75  static char   c;
76  static float gscale,ghscale;
77
78        displ = (sint *)ckalloc((2*max_aln_length+1) * sizeof(sint));
79        HH = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
80        DD = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
81        RR = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
82        SS = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
83               
84#ifdef MAC
85       int_scale = 10;
86#else
87       int_scale = 100;
88#endif
89  gscale=ghscale=1.0;
90  if (dnaflag)
91     {
92if (debug>1) fprintf(stdout,"matrix %s\n",pw_dnamtrxname);
93       if (strcmp(pw_dnamtrxname, "iub") == 0)
94          { 
95             matptr = swgapdnamt;
96             mat_xref = def_dna_xref;
97          }
98       else if (strcmp(pw_dnamtrxname, "clustalw") == 0)
99          { 
100             matptr = clustalvdnamt;
101             mat_xref = def_dna_xref;
102             gscale=0.6667;
103             ghscale=0.751;
104          }
105       else
106          {
107             matptr = pw_userdnamat;
108             mat_xref = pw_dna_xref;
109          }
110            maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
111            if (maxres == 0) return((sint)-1);
112
113            matrix[0][4]=transition_weight*matrix[0][0];
114            matrix[4][0]=transition_weight*matrix[0][0];
115            matrix[2][11]=transition_weight*matrix[0][0];
116            matrix[11][2]=transition_weight*matrix[0][0];
117            matrix[2][12]=transition_weight*matrix[0][0];
118            matrix[12][2]=transition_weight*matrix[0][0];
119    }
120  else
121    {
122if (debug>1) fprintf(stdout,"matrix %s\n",pw_mtrxname);
123       if (strcmp(pw_mtrxname, "blosum") == 0)
124          {
125             matptr = blosum30mt;
126             mat_xref = def_aa_xref;
127          }
128       else if (strcmp(pw_mtrxname, "pam") == 0)
129          {
130             matptr = pam350mt;
131             mat_xref = def_aa_xref;
132          }
133       else if (strcmp(pw_mtrxname, "gonnet") == 0)
134          {
135             matptr = gon250mt;
136             int_scale /= 10;
137             mat_xref = def_aa_xref;
138          }
139       else if (strcmp(pw_mtrxname, "id") == 0)
140          {
141             matptr = idmat;
142             mat_xref = def_aa_xref;
143          }
144       else
145          {
146             matptr = pw_usermat;
147             mat_xref = pw_aa_xref;
148          }
149
150       maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
151       if (maxres == 0) return((sint)-1);
152    }
153
154
155  for (si=MAX(0,istart);si<nseqs && si<iend;si++)
156   {
157     n = seqlen_array[si+1];
158     len1 = 0;
159     for (i=1;i<=n;i++) {
160                c = seq_array[si+1][i];
161                if ((c!=gap_pos1) && (c != gap_pos2)) len1++;
162     }
163
164     for (sj=MAX(si+1,jstart+1);sj<nseqs && sj<jend;sj++)
165      {
166        m = seqlen_array[sj+1];
167                len2 = 0;
168                for (i=1;i<=m;i++) {
169                        c = seq_array[sj+1][i];
170                        if ((c!=gap_pos1) && (c != gap_pos2)) len2++;
171                }
172
173        if (dnaflag) {
174           g = 200 * (float)pw_go_penalty*gscale;
175           gh = pw_ge_penalty * int_scale*ghscale;
176        }
177        else {
178           if (mat_avscore <= 0)
179              g = 200 * (float)(pw_go_penalty + log((double)(MIN(n,m))));
180           else
181              g = 2 * mat_avscore * (float)(pw_go_penalty +
182                    log((double)(MIN(n,m))));
183           gh = pw_ge_penalty * int_scale;
184        }
185
186if (debug>1) fprintf(stdout,"go %d ge %d\n",(pint)g,(pint)gh);
187
188/*
189   align the sequences
190*/
191        seq1 = si+1;
192        seq2 = sj+1;
193
194        forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
195           n, m);
196
197        reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0]);
198
199        last_print = 0;
200        print_ptr = 1;
201/*
202        sb1 = sb2 = 1;
203        se1 = n-1;
204        se2 = m-1;
205*/
206
207/* use Myers and Miller to align two sequences */
208
209        maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
210        (sint)0, (sint)0);
211 
212/* calculate percentage residue identity */
213
214        mm_score = tracepath(sb1,sb2);
215
216                if(len1==0 || len2==0) mm_score=0;
217                else
218                        mm_score /= (float)MIN(len1,len2);
219
220        tmat[si+1][sj+1] = ((float)100.0 - mm_score)/(float)100.0;
221        tmat[sj+1][si+1] = ((float)100.0 - mm_score)/(float)100.0;
222
223if (debug>1)
224{
225        fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
226                           (pint)si+1,(pint)sj+1, 
227                           (pint)mm_score, 
228                           (pint)maxscore/(MIN(len1,len2)*100));
229}
230else
231{
232        info("Sequences (%d:%d) Aligned. Score:  %d",
233                                      (pint)si+1,(pint)sj+1, 
234                                      (pint)mm_score);
235}
236
237   }
238  }
239   displ=ckfree((void *)displ);
240   HH=ckfree((void *)HH);
241   DD=ckfree((void *)DD);
242   RR=ckfree((void *)RR);
243   SS=ckfree((void *)SS);
244
245
246  return((sint)1);
247}
248
249static void add(sint v)
250{
251
252        if(last_print<0) {
253                displ[print_ptr-1] = v;
254                displ[print_ptr++] = last_print;
255        }
256        else
257                last_print = displ[print_ptr++] = v;
258}
259
260static sint calc_score(sint iat,sint jat,sint v1,sint v2)
261{
262        sint ipos,jpos;
263                sint ret;
264
265        ipos = v1 + iat;
266        jpos = v2 + jat;
267
268        ret=matrix[(int)seq_array[seq1][ipos]][(int)seq_array[seq2][jpos]];
269
270        return(ret);
271}
272
273
274static float tracepath(sint tsb1,sint tsb2)
275{
276        char c1,c2;
277    sint  i1,i2;
278    sint i,k,pos,to_do;
279        sint count;
280        float score;
281/*      char *s1, *s2;
282*/
283        to_do=print_ptr-1;
284        i1 = tsb1;
285        i2 = tsb2;
286
287        pos = 0;
288        count = 0;
289        for(i=1;i<=to_do;++i) {
290
291if (debug>1) fprintf(stdout,"%d ",(pint)displ[i]);
292                if(displ[i]==0) {
293                        c1 = seq_array[seq1][i1];
294                        c2 = seq_array[seq2][i2];
295/*
296if (debug>1)
297{
298if (c1>max_aa) s1[pos] = '-';
299else s1[pos]=amino_acid_codes[c1];
300if (c2>max_aa) s2[pos] = '-';
301else s2[pos]=amino_acid_codes[c2];
302}
303*/
304                        if ((c1!=gap_pos1) && (c1 != gap_pos2) &&
305                                    (c1 == c2)) count++;
306                        ++i1;
307                        ++i2;
308                        ++pos;
309                }
310                else {
311                        if((k=displ[i])>0) {
312/*
313if (debug>1)
314for (r=0;r<k;r++)
315{
316s1[pos+r]='-';
317if (seq_array[seq2][i2+r]>max_aa) s2[pos+r] = '-';
318else s2[pos+r]=amino_acid_codes[seq_array[seq2][i2+r]];
319}
320*/
321                                i2 += k;
322                                pos += k;
323                                count--;
324                        }
325                        else {
326/*
327if (debug>1)
328for (r=0;r<(-k);r++)
329{
330s2[pos+r]='-';
331if (seq_array[seq1][i1+r]>max_aa) s1[pos+r] = '-';
332else s1[pos+r]=amino_acid_codes[seq_array[seq1][i1+r]];
333}
334*/
335                                i1 -= k;
336                                pos -= k;
337                                count--;
338                        }
339                }
340        }
341/*
342if (debug>1) fprintf(stdout,"\n");
343if (debug>1)
344{
345for (i=0;i<pos;i++) fprintf(stdout,"%c",s1[i]);
346fprintf(stdout,"\n");
347for (i=0;i<pos;i++) fprintf(stdout,"%c",s2[i]);
348fprintf(stdout,"\n");
349}
350*/
351        if (count <= 0) count = 1;
352        score = 100.0 * (float)count;
353        return(score);
354}
355
356
357static void forward_pass(char *ia, char *ib, sint n, sint m)
358{
359
360  sint i,j;
361  pwint f,hh,p,t;
362
363  maxscore = 0;
364  se1 = se2 = 0;
365  for (i=0;i<=m;i++)
366    {
367       HH[i] = 0;
368       DD[i] = -g;
369    }
370
371  for (i=1;i<=n;i++)
372     {
373        hh = p = 0;
374                f = -g;
375
376        for (j=1;j<=m;j++)
377           {
378
379              f -= gh; 
380              t = hh - g - gh;
381              if (f<t) f = t;
382
383              DD[j] -= gh;
384              t = HH[j] - g - gh;
385              if (DD[j]<t) DD[j] = t;
386
387              hh = p + matrix[(int)ia[i]][(int)ib[j]];
388              if (hh<f) hh = f;
389              if (hh<DD[j]) hh = DD[j];
390              if (hh<0) hh = 0;
391
392              p = HH[j];
393              HH[j] = hh;
394
395              if (hh > maxscore)
396                {
397                   maxscore = hh;
398                   se1 = i;
399                   se2 = j;
400                }
401           }
402     }
403
404}
405
406
407static void reverse_pass(char *ia, char *ib)
408{
409
410  sint i,j;
411  pwint f,hh,p,t;
412  pwint cost;
413
414  cost = 0;
415  sb1 = sb2 = 0;
416  for (i=se2;i>0;i--)
417    {
418       HH[i] = -1;
419       DD[i] = -1;
420    }
421
422  for (i=se1;i>0;i--)
423     {
424        hh = f = -1;
425        if (i == se1) p = 0;
426        else p = -1;
427
428        for (j=se2;j>0;j--)
429           {
430
431              f -= gh; 
432              t = hh - g - gh;
433              if (f<t) f = t;
434
435              DD[j] -= gh;
436              t = HH[j] - g - gh;
437              if (DD[j]<t) DD[j] = t;
438
439              hh = p + matrix[(int)ia[i]][(int)ib[j]];
440              if (hh<f) hh = f;
441              if (hh<DD[j]) hh = DD[j];
442
443              p = HH[j];
444              HH[j] = hh;
445
446              if (hh > cost)
447                {
448                   cost = hh;
449                   sb1 = i;
450                   sb2 = j;
451                   if (cost >= maxscore) break;
452                }
453           }
454        if (cost >= maxscore) break;
455     }
456
457}
458
459static int diff(sint A,sint B,sint M,sint N,sint tb,sint te)
460{
461                sint type;
462        sint midi,midj,i,j;
463        int midh;
464        static pwint f, hh, e, s, t;
465
466        if(N<=0)  {
467                if(M>0) {
468                        del(M);
469                }
470
471                return(-(int)tbgap(M));
472        }
473
474        if(M<=1) {
475                if(M<=0) {
476                        add(N);
477                        return(-(int)tbgap(N));
478                }
479
480                midh = -(tb+gh) - tegap(N);
481                hh = -(te+gh) - tbgap(N);
482                if (hh>midh) midh = hh;
483                midj = 0;
484                for(j=1;j<=N;j++) {
485                        hh = calc_score(1,j,A,B)
486                                    - tegap(N-j) - tbgap(j-1);
487                        if(hh>midh) {
488                                midh = hh;
489                                midj = j;
490                        }
491                }
492
493                if(midj==0) {
494                        del(1);
495                        add(N);
496                }
497                else {
498                        if(midj>1)
499                                add(midj-1);
500                        displ[print_ptr++] = last_print = 0;
501                        if(midj<N)
502                                add(N-midj);
503                }
504                return midh;
505        }
506
507/* Divide: Find optimum midpoint (midi,midj) of cost midh */
508
509        midi = M / 2;
510        HH[0] = 0.0;
511        t = -tb;
512        for(j=1;j<=N;j++) {
513                HH[j] = t = t-gh;
514                DD[j] = t-g;
515        }
516
517        t = -tb;
518        for(i=1;i<=midi;i++) {
519                s=HH[0];
520                HH[0] = hh = t = t-gh;
521                f = t-g;
522                for(j=1;j<=N;j++) {
523                        if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
524                        if ((hh=HH[j]-g-gh) > (e=DD[j]-gh)) e=hh;
525                        hh = s + calc_score(i,j,A,B);
526                        if (f>hh) hh = f;
527                        if (e>hh) hh = e;
528
529                        s = HH[j];
530                        HH[j] = hh;
531                        DD[j] = e;
532                }
533        }
534
535        DD[0]=HH[0];
536
537        RR[N]=0;
538        t = -te;
539        for(j=N-1;j>=0;j--) {
540                RR[j] = t = t-gh;
541                SS[j] = t-g;
542        }
543
544        t = -te;
545        for(i=M-1;i>=midi;i--) {
546                s = RR[N];
547                RR[N] = hh = t = t-gh;
548                f = t-g;
549
550                for(j=N-1;j>=0;j--) {
551
552                        if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
553                        if ((hh=RR[j]-g-gh) > (e=SS[j]-gh)) e=hh;
554                        hh = s + calc_score(i+1,j+1,A,B);
555                        if (f>hh) hh = f;
556                        if (e>hh) hh = e;
557
558                        s = RR[j];
559                        RR[j] = hh;
560                        SS[j] = e;
561
562                }
563        }
564
565        SS[N]=RR[N];
566
567        midh=HH[0]+RR[0];
568        midj=0;
569        type=1;
570        for(j=0;j<=N;j++) {
571                hh = HH[j] + RR[j];
572                if(hh>=midh)
573                        if(hh>midh || (HH[j]!=DD[j] && RR[j]==SS[j])) {
574                                midh=hh;
575                                midj=j;
576                        }
577        }
578
579        for(j=N;j>=0;j--) {
580                hh = DD[j] + SS[j] + g;
581                if(hh>midh) {
582                        midh=hh;
583                        midj=j;
584                        type=2;
585                }
586        }
587
588        /* Conquer recursively around midpoint  */
589
590
591        if(type==1) {             /* Type 1 gaps  */
592                diff(A,B,midi,midj,tb,g);
593                diff(A+midi,B+midj,M-midi,N-midj,g,te);
594        }
595        else {
596                diff(A,B,midi-1,midj,tb,0.0);
597                del(2);
598                diff(A+midi+1,B+midj,M-midi-1,N-midj,0.0,te);
599        }
600
601        return midh;       /* Return the score of the best alignment */
602}
603
604static void del(sint k)
605{
606        if(last_print<0)
607                last_print = displ[print_ptr-1] -= k;
608        else
609                last_print = displ[print_ptr++] = -(k);
610}
611
612
Note: See TracBrowser for help on using the repository browser.