source: tags/initial/GDE/CLUSTAL/trees.c

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

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1        /* Phyle of filogenetic tree calculating functions for CLUSTAL V */
2#include <stdio.h>
3#include <string.h>
4#include <stdlib.h>
5#include <math.h>
6#include "clustalv.h"
7
8/*
9 *   Prototypes
10 */
11
12extern void *ckalloc(size_t);
13extern int getint(char *, int, int, int);
14extern void get_path(char *, char *);
15extern FILE * open_output_file(char *, char *, char *, char *);
16
17
18void init_trees(void);
19void phylogenetic_tree(void);
20void bootstrap_tree(void);
21void compare_tree(char **, char **, int *, int);
22void tree_gap_delete(void); /*flag all positions in alignment that have a gap */
23void dna_distance_matrix(FILE *);
24void prot_distance_matrix(FILE *);
25void nj_tree(char **, FILE *);
26void print_tree(char **, FILE *, int *);
27void root_tree(char **, int);
28
29Boolean transition(int,int);
30
31/*
32 *   Global variables
33 */
34
35
36extern double **tmat;     /* general nxn array of reals; allocated from main */
37                          /* this is used as a distance matrix */
38extern double **smat;
39extern Boolean dnaflag;   /* TRUE for DNA seqs; FALSE for proteins */
40extern Boolean tossgaps;  /* Ignore places in align. where ANY seq. has a gap*/
41extern Boolean kimura;    /* Use correction for multiple substitutions */
42extern Boolean empty;     /* any sequences in memory? */
43extern Boolean usemenu;   /* interactive (TRUE) or command line (FALSE) */
44extern void error(char *,...);  /* error reporting */
45extern int nseqs;         /* total no. of seqs. */
46extern int *seqlen_array; /* the lengths of the sequences */
47extern char **seq_array;   /* the sequences */
48extern char **names;       /* the seq. names */
49extern char seqname[];          /* name of input file */
50
51static double *av;
52static int *kill;
53static int ran_factor;
54int boot_ntrials;                       /* number of bootstrap trials */
55unsigned int boot_ran_seed;             /* random number generator seed */
56static int root_sequence;
57static int *boot_positions;
58static int *boot_totals;
59static char **standard_tree;
60static char **sample_tree;
61static FILE *phy_tree_file;
62static char phy_tree_name[FILENAMELEN+1];
63static Boolean verbose;
64static char *tree_gaps;    /* array of weights; 1 for use this posn.; 0 don't */
65
66
67void init_trees()
68{
69        register int i;
70       
71        tree_gaps = (char *)ckalloc( (MAXLEN+1) * sizeof (char) );
72       
73        boot_positions = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
74        boot_totals    = (int *)ckalloc( (MAXN+1) * sizeof (int) );
75
76        kill = (int *) ckalloc( (MAXN+1) * sizeof (int) );
77        av   = (double *) ckalloc( (MAXN+1) * sizeof (double)   );
78
79        standard_tree = (char **) ckalloc( (MAXN+1) * sizeof (char *) );
80        for(i=0; i<MAXN+1; i++) 
81                standard_tree[i] = (char *) ckalloc( (MAXN+1) * sizeof(char) );
82
83        sample_tree   = (char **) ckalloc( (MAXN+1) * sizeof (char *) );
84        for(i=0; i<MAXN+1; i++) 
85                sample_tree[i]   = (char *) ckalloc( (MAXN+1) * sizeof(char) );
86
87        boot_ntrials  = 1000;
88        boot_ran_seed = 111;
89        kimura   = FALSE;
90        tossgaps = FALSE;       
91}
92
93
94
95void phylogenetic_tree()
96{       char path[FILENAMELEN+1];
97        int j;
98
99        if(empty) {
100                error("You must load an alignment first");
101                return;
102        }
103
104        get_path(seqname,path);
105       
106        if((phy_tree_file = open_output_file(
107                "\nEnter name for tree output file  ",path,
108                phy_tree_name,"nj")) == NULL) return;
109
110        for(j=1; j<=seqlen_array[1]; ++j) 
111                boot_positions[j] = j;         
112
113        verbose = TRUE;                  /* Turn on screen output */
114        if(dnaflag)
115                dna_distance_matrix(phy_tree_file);
116        else 
117                prot_distance_matrix(phy_tree_file);
118
119        verbose = TRUE;                   /* Turn on output */
120        nj_tree(standard_tree,phy_tree_file);
121        fclose(phy_tree_file); 
122/*
123        print_tree(standard_tree,phy_tree_file);
124*/
125        fprintf(stdout,"\nPhylogenetic tree file created:   [%s]",phy_tree_name);
126}
127
128
129
130
131
132Boolean transition(int base1, int base2) /* TRUE if transition; else FALSE */
133/*
134
135   assumes that the bases of DNA sequences have been translated as
136   a,A = 1;   c,C = 2;   g,G = 3;   t,T,u,U = 4;  X or N = 0;  "-" < 0;
137
138   A <--> G  and  T <--> C  are transitions;  all others are transversions.
139
140*/
141{
142        if( ((base1 == 1) && (base2 == 3)) || ((base1 == 3) && (base2 == 1)) )
143                return TRUE;                                     /* A <--> G */
144        if( ((base1 == 4) && (base2 == 2)) || ((base1 == 2) && (base2 == 4)) )
145                return TRUE;                                     /* T <--> C */
146    return FALSE;
147}
148
149
150void tree_gap_delete()   /* flag all positions in alignment that have a gap */
151{                         /* in ANY sequence */
152        int seqn, posn;
153
154        for(posn=1; posn<=seqlen_array[1]; ++posn) {
155                tree_gaps[posn] = 0;
156        for(seqn=1; seqn<=nseqs; ++seqn)  {
157                        if(seq_array[seqn][posn] <= 0) {
158                           tree_gaps[posn] = 1;
159                                break;
160                        }
161                }
162        }
163}
164
165
166
167void nj_tree(char **tree_description, FILE *tree)
168{
169        register int i;
170        int l[4],nude,k;
171        int nc,mini,minj,j,j1,ii,jj;
172        double fnseqs,fnseqs2,sumd;
173        double diq,djq,dij,d2r,dr,dio,djo,da;
174        double tmin,total,dmin;
175        double bi,bj,b1,b2,b3,branch[4];
176        int typei,typej,type[4];             /* 0 = node; 1 = OTU */
177       
178        fnseqs = (double)nseqs;
179
180/*********************** First initialisation ***************************/
181       
182        if(verbose)  {
183                fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
184                fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
185                fprintf(tree," The Neighbor-joining Method:");
186                fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
187                fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
188                fprintf(tree,"\n\n This is an UNROOTED tree\n");
189                fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
190        }       
191
192        mini = minj = 0;
193
194        for(i=1;i<=nseqs;++i) 
195                {
196                tmat[i][i] = av[i] = 0.0;
197                kill[i] = 0;
198                }
199
200/*********************** Enter The Main Cycle ***************************/
201
202 /*     for(nc=1; nc<=(nseqs-3); ++nc) {  */                    /**start main cycle**/
203        for(nc=1; nc<=(nseqs-3); ++nc) {
204                sumd = 0.0;
205                for(j=2; j<=nseqs; ++j)
206                        for(i=1; i<j; ++i) {
207                                tmat[j][i] = tmat[i][j];
208                                sumd = sumd + tmat[i][j];
209                                smat[i][j] = smat[j][i] = 0.0;
210                        }
211
212                tmin = 99999.0;
213
214/*.................compute SMATij values and find the smallest one ........*/
215
216                for(jj=2; jj<=nseqs; ++jj) 
217                        if(kill[jj] != 1) 
218                                for(ii=1; ii<jj; ++ii)
219                                        if(kill[ii] != 1) {
220                                                diq = djq = 0.0;
221
222                                                for(i=1; i<=nseqs; ++i) {
223                                                        diq = diq + tmat[i][ii];
224                                                        djq = djq + tmat[i][jj];
225                                                }
226
227                                                dij = tmat[ii][jj];
228                                                d2r = diq + djq - (2.0*dij);
229                                                dr  = sumd - dij -d2r;
230                                                fnseqs2 = fnseqs - 2.0;
231                                                total= d2r+ fnseqs2*dij +dr*2.0;
232                                                total= total / (2.0*fnseqs2);
233                                                smat[ii][jj] = total;
234
235                                                if(total < tmin) {
236                                                        tmin = total;
237                                                        mini = ii;
238                                                        minj = jj;
239                                                }
240                                        }
241               
242
243/*.................compute branch lengths and print the results ........*/
244
245
246                dio = djo = 0.0;
247                for(i=1; i<=nseqs; ++i) {
248                        dio = dio + tmat[i][mini];
249                        djo = djo + tmat[i][minj];
250                }
251
252                dmin = tmat[mini][minj];
253                dio = (dio - dmin) / fnseqs2;
254                djo = (djo - dmin) / fnseqs2;
255                bi = (dmin + dio - djo) * 0.5;
256                bj = dmin - bi;
257                bi = bi - av[mini];
258                bj = bj - av[minj];
259
260                if( av[mini] > 0.0 )
261                        typei = 0;
262                else
263                        typei = 1;
264                if( av[minj] > 0.0 )
265                        typej = 0;
266                else
267                        typej = 1;
268
269                if(verbose) {
270                    fprintf(tree,"\n Cycle%4d     = ",nc);
271                    if(typei == 0)
272                        fprintf(tree,"Node:%4d (%9.5f) joins ",mini,bi);
273                    else
274                        fprintf(tree," SEQ:%4d (%9.5f) joins ",mini,bi);
275                    if(typej == 0) 
276                        fprintf(tree,"Node:%4d (%9.5f)",minj,bj);
277                    else
278                        fprintf(tree," SEQ:%4d (%9.5f)",minj,bj);
279                        fprintf(tree,"\n");
280                }
281
282                for(i=1; i<=nseqs; i++)
283                        tree_description[nc][i] = 0;
284
285                if(typei == 0) { 
286                        for(i=nc-1; i>=1; i--)
287                                if(tree_description[i][mini] == 1) {
288                                        for(j=1; j<=nseqs; j++) 
289                                             if(tree_description[i][j] == 1)
290                                                    tree_description[nc][j] = 1;
291                                        break;
292                                }
293                }
294                else
295                        tree_description[nc][mini] = 1;
296
297                if(typej == 0) {
298                        for(i=nc-1; i>=1; i--) 
299                                if(tree_description[i][minj] == 1) {
300                                        for(j=1; j<=nseqs; j++) 
301                                             if(tree_description[i][j] == 1)
302                                                    tree_description[nc][j] = 1;
303                                        break;
304                                }
305                }
306                else
307                        tree_description[nc][minj] = 1;
308                       
309                if(dmin <= 0.0) dmin = 0.0001;
310                av[mini] = dmin * 0.5;
311
312/*........................Re-initialisation................................*/
313
314                fnseqs = fnseqs - 1.0;
315                kill[minj] = 1;
316
317                for(j=1; j<=nseqs; ++j) 
318                        if( kill[j] != 1 ) {
319                                da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
320                                if( (mini - j) < 0 ) 
321                                        tmat[mini][j] = da;
322                                if( (mini - j) > 0)
323                                        tmat[j][mini] = da;
324                        }
325
326                for(j=1; j<=nseqs; ++j)
327                        tmat[minj][j] = tmat[j][minj] = 0.0;
328
329
330/****/  }                                               /**end main cycle**/
331
332/******************************Last Cycle (3 Seqs. left)********************/
333
334        nude = 1;
335
336        for(i=1; i<=nseqs; ++i)
337                if( kill[i] != 1 ) {
338                        l[nude] = i;
339                        nude = nude + 1;
340                }
341
342        b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
343        b2 =  tmat[l[1]][l[2]] - b1;
344        b3 =  tmat[l[1]][l[3]] - b1;
345        branch[1] = b1 - av[l[1]];
346        branch[2] = b2 - av[l[2]];
347        branch[3] = b3 - av[l[3]];
348
349
350        for(i=1; i<=nseqs; i++)
351                tree_description[nseqs-2][i] = 0;
352
353        if(verbose)
354                fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",nc);
355
356        for(i=1; i<=3; ++i) {
357           if( av[l[i]] > 0.0) {
358                if(verbose)
359                    fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",l[i],branch[i]);
360                for(k=nseqs-3; k>=1; k--)
361                        if(tree_description[k][l[i]] == 1) {
362                                for(j=1; j<=nseqs; j++)
363                                        if(tree_description[k][j] == 1)
364                                            tree_description[nseqs-2][j] = i;
365                                break;
366                        }
367           }
368           else  {
369                if(verbose)
370                    fprintf(tree,"\n\t\t  SEQ:%4d (%9.5f) ",l[i],branch[i]);
371                tree_description[nseqs-2][l[i]] = i;
372           }
373           if(i < 3) {
374                if(verbose)
375                    fprintf(tree,"joins");
376           }
377        }
378
379        if(verbose)
380                fprintf(tree,"\n");
381
382}
383
384
385
386
387void bootstrap_tree()
388{
389        int i,j,ranno;
390        int k,l,m,p;
391        unsigned int num;
392        char lin2[MAXLINE],path[MAXLINE+1];
393
394        if(empty) {
395                error("You must load an alignment first");
396                return;
397        }
398
399        get_path(seqname, path);
400       
401        if((phy_tree_file = open_output_file(
402                "\nEnter name for bootstrap output file  ",path,
403                phy_tree_name,"njb")) == NULL) return;
404
405        for(i=0;i<MAXN+1;i++)
406                boot_totals[i]=0;
407               
408        for(j=1; j<=seqlen_array[1]; ++j)  /* First select all positions for */
409                boot_positions[j] = j;     /* the "standard" tree */
410
411        verbose = TRUE;                    /* Turn on screen output */
412        if(dnaflag)
413                dna_distance_matrix(phy_tree_file);
414        else 
415                prot_distance_matrix(phy_tree_file);
416
417        verbose = TRUE;                         /* Turn on screen output */
418        nj_tree(standard_tree, phy_tree_file);  /* compute the standard tree */
419
420        fprintf(phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n");
421
422        ran_factor = RAND_MAX / seqlen_array[1];
423
424        if(usemenu) 
425                boot_ran_seed = 
426getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed);
427
428        srand(boot_ran_seed);
429        fprintf(phy_tree_file,"\n Random number generator seed = %7u\n",
430        boot_ran_seed);
431
432        if(usemenu) 
433                boot_ntrials = 
434getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials);
435
436        fprintf(phy_tree_file,"\n Number of bootstrap trials   = %7d\n",
437        boot_ntrials);
438
439        fprintf(phy_tree_file,
440        "\n\n Diagrammatic representation of the above tree: \n");
441        fprintf(phy_tree_file,"\n Each row represents 1 tree cycle;");
442        fprintf(phy_tree_file," defining 2 groups.\n");
443        fprintf(phy_tree_file,"\n Each column is 1 sequence; ");
444        fprintf(phy_tree_file,"the stars in each line show 1 group; ");
445        fprintf(phy_tree_file,"\n the dots show the other\n");
446        fprintf(phy_tree_file,"\n Numbers show occurences in bootstrap samples.");
447/*
448        print_tree(standard_tree, phy_tree_file, boot_totals);
449*/
450        verbose = FALSE;                   /* Turn OFF screen output */
451
452        fprintf(stdout,"\n\nEach dot represents 10 trials\n\n");
453        for(i=1; i<=boot_ntrials; ++i) {
454                for(j=1; j<=seqlen_array[1]; ++j) {       /* select alignment */
455                        ranno = ( rand() / ran_factor ) + 1; /* positions for */
456                        boot_positions[j] = ranno;        /* bootstrap sample */
457                }
458                if(dnaflag)
459                        dna_distance_matrix(phy_tree_file);
460                else
461                        prot_distance_matrix(phy_tree_file);
462                nj_tree(sample_tree, phy_tree_file); /* compute 1 sample tree */
463                compare_tree(standard_tree, sample_tree, boot_totals, nseqs);
464                if(i % 10  == 0) fprintf(stdout,".");
465                if(i % 100 == 0) fprintf(stdout,"\n");
466        }
467
468/*
469        fprintf(phy_tree_file,"\n\n Bootstrap totals for each group\n");
470*/
471        print_tree(standard_tree, phy_tree_file, boot_totals);
472
473        fclose(phy_tree_file);
474
475        fprintf(stdout,"\n\nBootstrap output file completed       [%s]"
476                ,phy_tree_name);
477}
478
479
480void compare_tree(char **tree1, char **tree2, int *hits, int n)
481{       
482        int i,j,k,l;
483        int nhits1, nhits2;
484
485        for(i=1; i<=n-3; i++)  {
486                for(j=1; j<=n-3; j++)  {
487                        nhits1 = 0;
488                        nhits2 = 0;
489                        for(k=1; k<=n; k++) {
490                                if(tree1[i][k] == tree2[j][k]) nhits1++;
491                                if(tree1[i][k] != tree2[j][k]) nhits2++;
492                        }
493                        if((nhits1 == nseqs) || (nhits2 == nseqs)) hits[i]++;
494                }
495        }
496}
497
498
499
500
501void print_tree(char **tree_description, FILE *tree, int *totals)
502{
503        int row,col;
504
505        fprintf(tree,"\n");
506
507        for(row=1; row<=nseqs-3; row++)  {
508                fprintf(tree," \n");
509                for(col=1; col<=nseqs; col++) { 
510                        if(tree_description[row][col] == 0)
511                                fprintf(tree,"*");
512                        else
513                                fprintf(tree,".");
514                }
515                if(totals[row] > 0)
516                        fprintf(tree,"%7d",totals[row]);
517        }
518        fprintf(tree," \n");
519        for(col=1; col<=nseqs; col++) 
520                fprintf(tree,"%1d",tree_description[nseqs-2][col]);
521        fprintf(tree,"\n");
522}
523
524
525
526void dna_distance_matrix(FILE *tree)
527{   
528        int m,n,j,i;
529        int res1, res2;
530        double p,q,e,a,b,k;     
531
532        tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
533       
534        if(verbose) {
535                fprintf(tree,"\n");
536                fprintf(tree,"\n DIST   = percentage divergence (/100)");
537                fprintf(tree,"\n p      = rate of transition (A <-> G; C <-> T)");
538                fprintf(tree,"\n q      = rate of transversion");
539                fprintf(tree,"\n Length = number of sites used in comparison");
540                fprintf(tree,"\n");
541            if(tossgaps) {
542                fprintf(tree,"\n All sites with gaps (in any sequence) deleted!");
543                fprintf(tree,"\n");
544            }
545            if(kimura) {
546                fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:");
547                fprintf(tree,"\n\n Kimura, M. (1980)");
548                fprintf(tree," A simple method for estimating evolutionary ");
549                fprintf(tree,"rates of base");
550                fprintf(tree,"\n substitutions through comparative studies of ");
551                fprintf(tree,"nucleotide sequences.");
552                fprintf(tree,"\n J. Mol. Evol., 16, 111-120.");
553                fprintf(tree,"\n\n");
554            }
555        }
556
557        for(m=1;   m<nseqs;  ++m)     /* for every pair of sequence */
558        for(n=m+1; n<=nseqs; ++n) {
559                p = q = e = 0.0;
560                tmat[m][n] = tmat[n][m] = 0.0;
561                for(i=1; i<=seqlen_array[1]; ++i) {
562                        j = boot_positions[i];
563                        if(tossgaps && (tree_gaps[j] > 0) ) 
564                                goto skip;          /* gap position */
565                        res1 = seq_array[m][j];
566                        res2 = seq_array[n][j];
567                        if( (res1 < 1) || (res2 < 1) ) 
568                                goto skip;          /* gap in a seq*/
569                        e = e + 1.0;
570                        if(res1 != res2) {
571                                if(transition(res1,res2))
572                                        p = p + 1.0;
573                                else
574                                        q = q + 1.0;
575                        }
576                        skip:;
577                }
578
579
580        /* Kimura's 2 parameter correction for multiple substitutions */
581
582                if(!kimura) {
583                        k = (p+q)/e;
584                        if(p > 0.0)
585                                p = p/e;
586                        else
587                                p = 0.0;
588                        if(q > 0.0)
589                                q = q/e;
590                        else
591                                q = 0.0;
592                        tmat[m][n] = tmat[n][m] = k;
593                        if(verbose)                    /* if screen output */
594                                fprintf(tree,       
595             "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
596                                 ,m,n,k,p,q,e);
597                }
598                else {
599                        if(p > 0.0)
600                                p = p/e;
601                        else
602                                p = 0.0;
603                        if(q > 0.0)
604                                q = q/e;
605                        else
606                                q = 0.0;
607                        a = 1.0/(1.0-2.0*p-q);
608                        b = 1.0/(1.0-2.0*q);
609                        k = 0.5*log(a) + 0.25*log(b);
610                        tmat[m][n] = tmat[n][m] = k;
611                        if(verbose)                      /* if screen output */
612                                fprintf(tree,
613             "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
614                                ,m,n,k,p,q,e);
615
616                }
617        }
618}
619
620
621
622void prot_distance_matrix(FILE *tree)
623{
624        int m,n,j,i;
625        int res1, res2;
626        double p,e,a,b,k;       
627
628        tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
629       
630        if(verbose) {
631                fprintf(tree,"\n");
632                fprintf(tree,"\n DIST   = percentage divergence (/100)");
633                fprintf(tree,"\n Length = number of sites used in comparison");
634                fprintf(tree,"\n\n");
635                if(tossgaps) {
636                        fprintf(tree,"\n All sites with gaps (in any sequence) deleted");
637                        fprintf(tree,"\n");
638                }
639                if(kimura) {
640                        fprintf(tree,"\n Distances corrected by Kimura's empirical method:");
641                        fprintf(tree,"\n\n Kimura, M. (1983)");
642                        fprintf(tree," The Neutral Theory of Molecular Evolution.");
643                        fprintf(tree,"\n Cambridge University Press, Cambridge, England.");
644                        fprintf(tree,"\n\n");
645                }
646        }
647
648        for(m=1;   m<nseqs;  ++m)     /* for every pair of sequence */
649        for(n=m+1; n<=nseqs; ++n) {
650                p = e = 0.0;
651                tmat[m][n] = tmat[n][m] = 0.0;
652                for(i=1; i<=seqlen_array[1]; ++i) {
653                        j = boot_positions[i];
654                        if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */
655                        res1 = seq_array[m][j];
656                        res2 = seq_array[n][j];
657                        if( (res1 < 1) || (res2 < 1) )  goto skip;   /* gap in a seq*/
658                        e = e + 1.0;
659                        if(res1 != res2) p = p + 1.0;
660                        skip:;
661                }
662
663                if(p <= 0.0) 
664                        k = 0.0;
665                else
666                        k = p/e;
667
668                if(kimura) 
669                        if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) );
670
671                tmat[m][n] = tmat[n][m] = k;
672                    if(verbose)                    /* if screen output */
673                        fprintf(tree,       
674                         "%4d vs.%4d  DIST = %6.4f;  length = %6.0f\n",m,n,k,e);
675        }
676}
677
678
Note: See TracBrowser for help on using the repository browser.