Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

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

Visit:

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:keywords set to `Author Date Id Revision`
File size: 34.3 KB

Line
1	/* Phyle of filogenetic tree calculating functions for CLUSTAL W */
2	/* DES was here FEB. 1994 */
3
4	#include <stdio.h>
5	#include <string.h>
6	#include <stdlib.h>
7	#include <math.h>
8	#include "clustalw.h"
9	#include "dayhoff.h" /* set correction for amino acid distances >= 75% */
10
11
12	/*
13	* Prototypes
14	*/
15	Boolean transition(sint base1, sint base2);
16	void tree_gap_delete(void);
17	void distance_matrix_output(FILE *ofile);
18	void nj_tree(char *tree_description, FILE tree);
19	void compare_tree(char tree1, char tree2, sint *hits, sint n);
20	void print_phylip_tree(char *tree_description, FILE tree, sint bootstrap);
21	sint two_way_split(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap);
22	void print_tree(char *tree_description, FILE tree, sint *totals);
23
24	/*
25	* Global variables
26	*/
27
28	extern sint max_names;
29
30	extern double *tmat; / general nxn array of reals; allocated from main */
31	/* this is used as a distance matrix */
32	extern Boolean dnaflag; /* TRUE for DNA seqs; FALSE for proteins */
33	extern Boolean tossgaps; /* Ignore places in align. where ANY seq. has a gap*/
34	extern Boolean kimura; /* Use correction for multiple substitutions */
35	extern Boolean output_tree_clustal; /* clustal text output for trees */
36	extern Boolean output_tree_phylip; /* phylip nested parentheses format */
37	extern Boolean output_tree_distances; /* phylip distance matrix */
38	extern sint bootstrap_format; /* bootstrap file format */
39	extern Boolean empty; /* any sequences in memory? */
40	extern Boolean usemenu; /* interactive (TRUE) or command line (FALSE) */
41	extern sint nseqs;
42	extern sint max_aln_length;
43	extern sint seqlen_array; / the lengths of the sequences */
44	extern char *seq_array; / the sequences */
45	extern char *names; / the seq. names */
46	extern char seqname[]; /* name of input file */
47	extern sint gap_pos1,gap_pos2;
48
49	static double *av;
50	static double left_branch, right_branch;
51	static double save_left_branch, save_right_branch;
52	static sint *boot_totals;
53	static sint *tkill;
54	/*
55	The next line is a fossil from the days of using the cc ran()
56	static int ran_factor;
57	*/
58	static sint *boot_positions;
59	static FILE *phylip_phy_tree_file;
60	static FILE *clustal_phy_tree_file;
61	static FILE *distances_phy_tree_file;
62	static Boolean verbose;
63	static char *tree_gaps;
64	static sint first_seq, last_seq;
65	/* array of weights; 1 for use this posn.; 0 don't */
66
67	extern sint boot_ntrials; /* number of bootstrap trials */
68	extern unsigned sint boot_ran_seed; /* random number generator seed */
69
70	void phylogenetic_tree(char phylip_name,char clustal_name,char *dist_name)
71	/*
72	Calculate a tree using the distances in the nseqs*nseqs array tmat.
73	This is the routine for getting the REAL trees after alignment.
74	*/
75	{ char path[FILENAMELEN+1];
76	sint i, j;
77	sint overspill = 0;
78	sint total_dists;
79	static char **standard_tree;
80	char lin2[10];
81
82	if(empty) {
83	error("You must load an alignment first");
84	return;
85	}
86
87	if(nseqs<=3) {
88	error("Alignment has only %d sequences",nseqs);
89	return;
90	}
91	first_seq=1;
92	last_seq=nseqs;
93
94	get_path(seqname,path);
95
96	if(output_tree_clustal) {
97	if (clustal_name[0]!=EOS) {
98	if((clustal_phy_tree_file = open_explicit_file(
99	clustal_name))==NULL) return;
100	}
101	else {
102	if((clustal_phy_tree_file = open_output_file(
103	"\nEnter name for CLUSTAL tree output file ",path,
104	clustal_name,"nj")) == NULL) return;
105	}
106	}
107
108	if(output_tree_phylip) {
109	if (phylip_name[0]!=EOS) {
110	if((phylip_phy_tree_file = open_explicit_file(
111	phylip_name))==NULL) return;
112	}
113	else {
114	if((phylip_phy_tree_file = open_output_file(
115	"\nEnter name for PHYLIP tree output file ",path,
116	phylip_name,"ph")) == NULL) return;
117	}
118	}
119
120	if(output_tree_distances)
121	{
122	if (dist_name[0]!=EOS) {
123	if((distances_phy_tree_file = open_explicit_file(
124	dist_name))==NULL) return;
125	}
126	else {
127	if((distances_phy_tree_file = open_output_file(
128	"\nEnter name for distance matrix output file ",path,
129	dist_name,"dst")) == NULL) return;
130	}
131	}
132
133	boot_positions = (sint )ckalloc( (seqlen_array[first_seq]+2) sizeof (sint) );
134
135	for(j=1; j<=seqlen_array[first_seq]; ++j)
136	boot_positions[j] = j;
137
138	if(output_tree_clustal) {
139	verbose = TRUE; /* Turn on file output */
140	if(dnaflag)
141	overspill = dna_distance_matrix(clustal_phy_tree_file);
142	else
143	overspill = prot_distance_matrix(clustal_phy_tree_file);
144	}
145
146	if(output_tree_phylip) {
147	verbose = FALSE; /* Turn off file output */
148	if(dnaflag)
149	overspill = dna_distance_matrix(phylip_phy_tree_file);
150	else
151	overspill = prot_distance_matrix(phylip_phy_tree_file);
152	}
153
154	if(output_tree_distances) {
155	verbose = FALSE; /* Turn off file output */
156	if(dnaflag)
157	overspill = dna_distance_matrix(distances_phy_tree_file);
158	else
159	overspill = prot_distance_matrix(distances_phy_tree_file);
160	distance_matrix_output(distances_phy_tree_file);
161	}
162
163	/* check if any distances overflowed the distance corrections */
164	if ( overspill > 0 ) {
165	total_dists = (nseqs*(nseqs-1))/2;
166	fprintf(stdout,"\n");
167	fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld",
168	(long)overspill,(long)total_dists);
169	fprintf(stdout,"\n were out of range for the distance correction.");
170	fprintf(stdout,"\n This may not be fatal but you have been warned!");
171	fprintf(stdout,"\n");
172	fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
173	fprintf(stdout,"\n or 2) remove the most distant sequences");
174	fprintf(stdout,"\n or 3) use the PHYLIP package.");
175	fprintf(stdout,"\n\n");
176	if (usemenu)
177	getstr("Press [RETURN] to continue",lin2);
178	}
179
180	if(output_tree_clustal) verbose = TRUE; /* Turn on file output */
181
182	standard_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
183	for(i=0; i<nseqs+1; i++)
184	standard_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
185
186	if(output_tree_clustal \|\| output_tree_phylip)
187	nj_tree(standard_tree,clustal_phy_tree_file);
188
189	if(output_tree_phylip)
190	print_phylip_tree(standard_tree,phylip_phy_tree_file,0);
191
192	/*
193	print_tree(standard_tree,phy_tree_file);
194	*/
195	tree_gaps=ckfree((void *)tree_gaps);
196	boot_positions=ckfree((void *)boot_positions);
197	if (left_branch != NULL) left_branch=ckfree((void *)left_branch);
198	if (right_branch != NULL) right_branch=ckfree((void *)right_branch);
199	if (tkill != NULL) tkill=ckfree((void *)tkill);
200	if (av != NULL) av=ckfree((void *)av);
201	for (i=1;i<nseqs+1;i++)
202	standard_tree[i]=ckfree((void *)standard_tree[i]);
203	standard_tree=ckfree((void *)standard_tree);
204
205	if(output_tree_clustal) {
206	fclose(clustal_phy_tree_file);
207	info("Phylogenetic tree file created: [%s]",clustal_name);
208	}
209
210	if(output_tree_phylip) {
211	fclose(phylip_phy_tree_file);
212	info("Phylogenetic tree file created: [%s]",phylip_name);
213	}
214
215	if(output_tree_distances) {
216	fclose(distances_phy_tree_file);
217	info("Distance matrix file created: [%s]",dist_name);
218	}
219
220
221	}
222
223
224
225
226
227	Boolean transition(sint base1, sint base2) /* TRUE if transition; else FALSE */
228	/*
229
230	assumes that the bases of DNA sequences have been translated as
231	a,A = 0; c,C = 1; g,G = 2; t,T,u,U = 3; N = 4;
232
233	A <--> G and T <--> C are transitions; all others are transversions.
234
235	*/
236	{
237	if( ((base1 == 0) && (base2 == 2)) \|\| ((base1 == 2) && (base2 == 0)) )
238	return TRUE; /* A <--> G */
239	if( ((base1 == 3) && (base2 == 1)) \|\| ((base1 == 1) && (base2 == 3)) )
240	return TRUE; /* T <--> C */
241	return FALSE;
242	}
243
244
245	void tree_gap_delete(void) /* flag all positions in alignment that have a gap */
246	{ /* in ANY sequence */
247	sint seqn;
248	sint posn;
249
250	tree_gaps = (char )ckalloc( (max_aln_length+1) sizeof (char) );
251
252	for(posn=1; posn<=seqlen_array[first_seq]; ++posn) {
253	tree_gaps[posn] = 0;
254	for(seqn=1; seqn<=last_seq-first_seq+1; ++seqn) {
255	if((seq_array[seqn+first_seq-1][posn] == gap_pos1) \|\|
256	(seq_array[seqn+first_seq-1][posn] == gap_pos2)) {
257	tree_gaps[posn] = 1;
258	break;
259	}
260	}
261	}
262	}
263
264	void distance_matrix_output(FILE *ofile)
265	{
266	sint i,j;
267
268	fprintf(ofile,"%6d",(pint)last_seq-first_seq+1);
269	for(i=1;i<=last_seq-first_seq+1;i++) {
270	fprintf(ofile,"\n%-*s ",max_names,names[i]);
271	for(j=1;j<=last_seq-first_seq+1;j++) {
272	fprintf(ofile,"%6.3f ",tmat[i][j]);
273	if(j % 8 == 0) {
274	if(j!=last_seq-first_seq+1) fprintf(ofile,"\n");
275	if(j != last_seq-first_seq+1 ) fprintf(ofile," ");
276	}
277	}
278	}
279	}
280
281
282
283	void nj_tree(char *tree_description, FILE tree)
284	{
285	register int i;
286	sint l[4],nude,k;
287	sint nc,mini,minj,j,ii,jj;
288	double fnseqs,fnseqs2=0,sumd;
289	double diq,djq,dij,d2r,dr,dio,djo,da;
290	double tmin,total,dmin;
291	double bi,bj,b1,b2,b3,branch[4];
292	sint typei,typej; /* 0 = node; 1 = OTU */
293
294	fnseqs = (double)last_seq-first_seq+1;
295
296	/********************* First initialisation *************************/
297
298	if(verbose) {
299	fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
300	fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
301	fprintf(tree," The Neighbor-joining Method:");
302	fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
303	fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
304	fprintf(tree,"\n\n This is an UNROOTED tree\n");
305	fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
306	}
307
308	mini = minj = 0;
309
310	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
311	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
312	tkill = (sint ) ckalloc( (nseqs+1) sizeof (sint) );
313	av = (double ) ckalloc( (nseqs+1) sizeof (double) );
314
315	for(i=1;i<=last_seq-first_seq+1;++i)
316	{
317	tmat[i][i] = av[i] = 0.0;
318	tkill[i] = 0;
319	}
320
321	/********************* Enter The Main Cycle *************************/
322
323	/* for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) { / /start main cycle*/
324	for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) {
325	sumd = 0.0;
326	for(j=2; j<=last_seq-first_seq+1; ++j)
327	for(i=1; i<j; ++i) {
328	tmat[j][i] = tmat[i][j];
329	sumd = sumd + tmat[i][j];
330	}
331
332	tmin = 99999.0;
333
334	/.................compute SMATij values and find the smallest one ......../
335
336	for(jj=2; jj<=last_seq-first_seq+1; ++jj)
337	if(tkill[jj] != 1)
338	for(ii=1; ii<jj; ++ii)
339	if(tkill[ii] != 1) {
340	diq = djq = 0.0;
341
342	for(i=1; i<=last_seq-first_seq+1; ++i) {
343	diq = diq + tmat[i][ii];
344	djq = djq + tmat[i][jj];
345	}
346
347	dij = tmat[ii][jj];
348	d2r = diq + djq - (2.0*dij);
349	dr = sumd - dij -d2r;
350	fnseqs2 = fnseqs - 2.0;
351	total= d2r+ fnseqs2dij +dr2.0;
352	total= total / (2.0*fnseqs2);
353
354	if(total < tmin) {
355	tmin = total;
356	mini = ii;
357	minj = jj;
358	}
359	}
360
361
362	/.................compute branch lengths and print the results ......../
363
364
365	dio = djo = 0.0;
366	for(i=1; i<=last_seq-first_seq+1; ++i) {
367	dio = dio + tmat[i][mini];
368	djo = djo + tmat[i][minj];
369	}
370
371	dmin = tmat[mini][minj];
372	dio = (dio - dmin) / fnseqs2;
373	djo = (djo - dmin) / fnseqs2;
374	bi = (dmin + dio - djo) * 0.5;
375	bj = dmin - bi;
376	bi = bi - av[mini];
377	bj = bj - av[minj];
378
379	if( av[mini] > 0.0 )
380	typei = 0;
381	else
382	typei = 1;
383	if( av[minj] > 0.0 )
384	typej = 0;
385	else
386	typej = 1;
387
388	if(verbose)
389	fprintf(tree,"\n Cycle%4d = ",(pint)nc);
390
391	/*
392	set negative branch lengths to zero. Also set any tiny positive
393	branch lengths to zero.
394	*/ if( fabs(bi) < 0.0001) bi = 0.0;
395	if( fabs(bj) < 0.0001) bj = 0.0;
396
397	if(verbose) {
398	if(typei == 0)
399	fprintf(tree,"Node:%4d (%9.5f) joins ",(pint)mini,bi);
400	else
401	fprintf(tree," SEQ:%4d (%9.5f) joins ",(pint)mini,bi);
402
403	if(typej == 0)
404	fprintf(tree,"Node:%4d (%9.5f)",(pint)minj,bj);
405	else
406	fprintf(tree," SEQ:%4d (%9.5f)",(pint)minj,bj);
407
408	fprintf(tree,"\n");
409	}
410
411
412	left_branch[nc] = bi;
413	right_branch[nc] = bj;
414
415	for(i=1; i<=last_seq-first_seq+1; i++)
416	tree_description[nc][i] = 0;
417
418	if(typei == 0) {
419	for(i=nc-1; i>=1; i--)
420	if(tree_description[i][mini] == 1) {
421	for(j=1; j<=last_seq-first_seq+1; j++)
422	if(tree_description[i][j] == 1)
423	tree_description[nc][j] = 1;
424	break;
425	}
426	}
427	else
428	tree_description[nc][mini] = 1;
429
430	if(typej == 0) {
431	for(i=nc-1; i>=1; i--)
432	if(tree_description[i][minj] == 1) {
433	for(j=1; j<=last_seq-first_seq+1; j++)
434	if(tree_description[i][j] == 1)
435	tree_description[nc][j] = 1;
436	break;
437	}
438	}
439	else
440	tree_description[nc][minj] = 1;
441
442
443	/*
444	Here is where the -0.00005 branch lengths come from for 3 or more
445	identical seqs.
446	*/
447	/* if(dmin <= 0.0) dmin = 0.0001; */
448	if(dmin <= 0.0) dmin = 0.000001;
449	av[mini] = dmin * 0.5;
450
451	/........................Re-initialisation................................/
452
453	fnseqs = fnseqs - 1.0;
454	tkill[minj] = 1;
455
456	for(j=1; j<=last_seq-first_seq+1; ++j)
457	if( tkill[j] != 1 ) {
458	da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
459	if( (mini - j) < 0 )
460	tmat[mini][j] = da;
461	if( (mini - j) > 0)
462	tmat[j][mini] = da;
463	}
464
465	for(j=1; j<=last_seq-first_seq+1; ++j)
466	tmat[minj][j] = tmat[j][minj] = 0.0;
467
468
469	/**/ } /end main cycle**/
470
471	/****************************Last Cycle (3 Seqs. left)******************/
472
473	nude = 1;
474
475	for(i=1; i<=last_seq-first_seq+1; ++i)
476	if( tkill[i] != 1 ) {
477	l[nude] = i;
478	nude = nude + 1;
479	}
480
481	b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
482	b2 = tmat[l[1]][l[2]] - b1;
483	b3 = tmat[l[1]][l[3]] - b1;
484
485	branch[1] = b1 - av[l[1]];
486	branch[2] = b2 - av[l[2]];
487	branch[3] = b3 - av[l[3]];
488
489	/* Reset tiny negative and positive branch lengths to zero */
490	if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
491	if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
492	if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
493
494	left_branch[last_seq-first_seq+1-2] = branch[1];
495	left_branch[last_seq-first_seq+1-1] = branch[2];
496	left_branch[last_seq-first_seq+1] = branch[3];
497
498	for(i=1; i<=last_seq-first_seq+1; i++)
499	tree_description[last_seq-first_seq+1-2][i] = 0;
500
501	if(verbose)
502	fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",(pint)nc);
503
504	for(i=1; i<=3; ++i) {
505	if( av[l[i]] > 0.0) {
506	if(verbose)
507	fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",(pint)l[i],branch[i]);
508	for(k=last_seq-first_seq+1-3; k>=1; k--)
509	if(tree_description[k][l[i]] == 1) {
510	for(j=1; j<=last_seq-first_seq+1; j++)
511	if(tree_description[k][j] == 1)
512	tree_description[last_seq-first_seq+1-2][j] = i;
513	break;
514	}
515	}
516	else {
517	if(verbose)
518	fprintf(tree,"\n\t\t SEQ:%4d (%9.5f) ",(pint)l[i],branch[i]);
519	tree_description[last_seq-first_seq+1-2][l[i]] = i;
520	}
521	if(i < 3) {
522	if(verbose)
523	fprintf(tree,"joins");
524	}
525	}
526
527	if(verbose)
528	fprintf(tree,"\n");
529
530	}
531
532
533
534
535	void bootstrap_tree(char phylip_name,char clustal_name)
536	{
537	sint i,j;
538	int ranno;
539	char path[MAXLINE+1];
540	char dummy[10];
541	static char **sample_tree;
542	static char **standard_tree;
543	sint total_dists, overspill = 0, total_overspill = 0;
544	sint nfails = 0;
545
546	if(empty) {
547	error("You must load an alignment first");
548	return;
549	}
550
551	if(nseqs<=3) {
552	error("Alignment has only %d sequences",nseqs);
553	return;
554	}
555
556	if(!output_tree_clustal && !output_tree_phylip) {
557	error("You must select either clustal or phylip tree output format");
558	return;
559	}
560	get_path(seqname, path);
561
562	if (output_tree_clustal) {
563	if (clustal_name[0]!=EOS) {
564	if((clustal_phy_tree_file = open_explicit_file(
565	clustal_name))==NULL) return;
566	}
567	else {
568	if((clustal_phy_tree_file = open_output_file(
569	"\nEnter name for bootstrap output file ",path,
570	clustal_name,"njb")) == NULL) return;
571	}
572	}
573
574	first_seq=1;
575	last_seq=nseqs;
576
577	if (output_tree_phylip) {
578	if (phylip_name[0]!=EOS) {
579	if((phylip_phy_tree_file = open_explicit_file(
580	phylip_name))==NULL) return;
581	}
582	else {
583	if((phylip_phy_tree_file = open_output_file(
584	"\nEnter name for bootstrap output file ",path,
585	phylip_name,"phb")) == NULL) return;
586	}
587	}
588
589	boot_totals = (sint )ckalloc( (nseqs+1) sizeof (sint) );
590	for(i=0;i<nseqs+1;i++)
591	boot_totals[i]=0;
592
593	boot_positions = (sint )ckalloc( (seqlen_array[first_seq]+2) sizeof (sint) );
594
595	for(j=1; j<=seqlen_array[first_seq]; ++j) /* First select all positions for */
596	boot_positions[j] = j; /* the "standard" tree */
597
598	if(output_tree_clustal) {
599	verbose = TRUE; /* Turn on file output */
600	if(dnaflag)
601	overspill = dna_distance_matrix(clustal_phy_tree_file);
602	else
603	overspill = prot_distance_matrix(clustal_phy_tree_file);
604	}
605
606	if(output_tree_phylip) {
607	verbose = FALSE; /* Turn off file output */
608	if(dnaflag)
609	overspill = dna_distance_matrix(phylip_phy_tree_file);
610	else
611	overspill = prot_distance_matrix(phylip_phy_tree_file);
612	}
613
614	/* check if any distances overflowed the distance corrections */
615	if ( overspill > 0 ) {
616	total_dists = (nseqs*(nseqs-1))/2;
617	fprintf(stdout,"\n");
618	fprintf(stdout,"\n WARNING: %d of the distances out of a total of %d",
619	(pint)overspill,(pint)total_dists);
620	fprintf(stdout,"\n were out of range for the distance correction.");
621	fprintf(stdout,"\n This may not be fatal but you have been warned!");
622	fprintf(stdout,"\n");
623	fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
624	fprintf(stdout,"\n or 2) remove the most distant sequences");
625	fprintf(stdout,"\n or 3) use the PHYLIP package.");
626	fprintf(stdout,"\n\n");
627	if (usemenu)
628	getstr("Press [RETURN] to continue",dummy);
629	}
630
631	tree_gaps=ckfree((void *)tree_gaps);
632
633	if (output_tree_clustal) verbose = TRUE; /* Turn on screen output */
634
635	standard_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
636	for(i=0; i<nseqs+1; i++)
637	standard_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
638
639	/* compute the standard tree */
640
641	if(output_tree_clustal \|\| output_tree_phylip)
642	nj_tree(standard_tree,clustal_phy_tree_file);
643
644	if (output_tree_clustal)
645	fprintf(clustal_phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n");
646
647	/* save the left_branch and right_branch for phylip output */
648	save_left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
649	save_right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
650	for (i=1;i<=nseqs;i++) {
651	save_left_branch[i] = left_branch[i];
652	save_right_branch[i] = right_branch[i];
653	}
654	/*
655	The next line is a fossil from the days of using the cc ran()
656	ran_factor = RAND_MAX / seqlen_array[first_seq];
657	*/
658
659	if(usemenu)
660	boot_ran_seed =
661	getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed);
662
663	/* do not use the native cc ran()
664	srand(boot_ran_seed);
665	*/
666	addrandinit((unsigned long) boot_ran_seed);
667
668	if (output_tree_clustal)
669	fprintf(clustal_phy_tree_file,"\n Random number generator seed = %7u\n",
670	boot_ran_seed);
671
672	if(usemenu)
673	boot_ntrials =
674	getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials);
675
676	if (output_tree_clustal) {
677	fprintf(clustal_phy_tree_file,"\n Number of bootstrap trials = %7d\n",
678	(pint)boot_ntrials);
679
680	fprintf(clustal_phy_tree_file,
681	"\n\n Diagrammatic representation of the above tree: \n");
682	fprintf(clustal_phy_tree_file,"\n Each row represents 1 tree cycle;");
683	fprintf(clustal_phy_tree_file," defining 2 groups.\n");
684	fprintf(clustal_phy_tree_file,"\n Each column is 1 sequence; ");
685	fprintf(clustal_phy_tree_file,"the stars in each line show 1 group; ");
686	fprintf(clustal_phy_tree_file,"\n the dots show the other\n");
687	fprintf(clustal_phy_tree_file,"\n Numbers show occurences in bootstrap samples.");
688	}
689	/*
690	print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
691	*/
692	verbose = FALSE; /* Turn OFF screen output */
693
694	left_branch=ckfree((void *)left_branch);
695	right_branch=ckfree((void *)right_branch);
696	tkill=ckfree((void *)tkill);
697	av=ckfree((void *)av);
698
699	sample_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
700	for(i=0; i<nseqs+1; i++)
701	sample_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
702
703	if (usemenu)
704	fprintf(stdout,"\n\nEach dot represents 10 trials\n\n");
705	total_overspill = 0;
706	nfails = 0;
707	for(i=1; i<=boot_ntrials; ++i) {
708	for(j=1; j<=seqlen_array[first_seq]; ++j) { /* select alignment */
709	/* positions for */
710	ranno = addrand( (unsigned long) seqlen_array[1]) + 1;
711	boot_positions[j] = ranno; /* bootstrap sample */
712	}
713	if(output_tree_clustal) {
714	if(dnaflag)
715	overspill = dna_distance_matrix(clustal_phy_tree_file);
716	else
717	overspill = prot_distance_matrix(clustal_phy_tree_file);
718	}
719
720	if(output_tree_phylip) {
721	if(dnaflag)
722	overspill = dna_distance_matrix(phylip_phy_tree_file);
723	else
724	overspill = prot_distance_matrix(phylip_phy_tree_file);
725	}
726
727	if( overspill > 0) {
728	total_overspill = total_overspill + overspill;
729	nfails++;
730	}
731
732	tree_gaps=ckfree((void *)tree_gaps);
733
734	if(output_tree_clustal \|\| output_tree_phylip)
735	nj_tree(sample_tree,clustal_phy_tree_file);
736
737	left_branch=ckfree((void *)left_branch);
738	right_branch=ckfree((void *)right_branch);
739	tkill=ckfree((void *)tkill);
740	av=ckfree((void *)av);
741
742	compare_tree(standard_tree, sample_tree, boot_totals, last_seq-first_seq+1);
743	if (usemenu) {
744	if(i % 10 == 0) fprintf(stdout,".");
745	if(i % 100 == 0) fprintf(stdout,"\n");
746	}
747	}
748
749	/* check if any distances overflowed the distance corrections */
750	if ( nfails > 0 ) {
751	total_dists = (nseqs*(nseqs-1))/2;
752	fprintf(stdout,"\n");
753	fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld times %ld",
754	(long)total_overspill,(long)total_dists,(long)boot_ntrials);
755	fprintf(stdout,"\n were out of range for the distance correction.");
756	fprintf(stdout,"\n This affected %d out of %d bootstrap trials.",
757	(pint)nfails,(pint)boot_ntrials);
758	fprintf(stdout,"\n This may not be fatal but you have been warned!");
759	fprintf(stdout,"\n");
760	fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
761	fprintf(stdout,"\n or 2) remove the most distant sequences");
762	fprintf(stdout,"\n or 3) use the PHYLIP package.");
763	fprintf(stdout,"\n\n");
764	if (usemenu)
765	getstr("Press [RETURN] to continue",dummy);
766	}
767
768
769	boot_positions=ckfree((void *)boot_positions);
770
771	for (i=1;i<nseqs+1;i++)
772	sample_tree[i]=ckfree((void *)sample_tree[i]);
773	sample_tree=ckfree((void *)sample_tree);
774	/*
775	fprintf(clustal_phy_tree_file,"\n\n Bootstrap totals for each group\n");
776	*/
777	if (output_tree_clustal)
778	print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
779
780	if(output_tree_phylip) {
781	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
782	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
783	for (i=1;i<=nseqs;i++) {
784	left_branch[i] = save_left_branch[i];
785	right_branch[i] = save_right_branch[i];
786	}
787	print_phylip_tree(standard_tree,phylip_phy_tree_file,
788	bootstrap_format);
789	left_branch=ckfree((void *)left_branch);
790	right_branch=ckfree((void *)right_branch);
791	}
792
793	boot_totals=ckfree((void *)boot_totals);
794	save_left_branch=ckfree((void *)save_left_branch);
795	save_right_branch=ckfree((void *)save_right_branch);
796
797	for (i=1;i<nseqs+1;i++)
798	standard_tree[i]=ckfree((void *)standard_tree[i]);
799	standard_tree=ckfree((void *)standard_tree);
800
801	if (output_tree_clustal)
802	fclose(clustal_phy_tree_file);
803
804	if (output_tree_phylip)
805	fclose(phylip_phy_tree_file);
806
807	if (output_tree_clustal)
808	info("Bootstrap output file completed [%s]"
809	,clustal_name);
810	if (output_tree_phylip)
811	info("Bootstrap output file completed [%s]"
812	,phylip_name);
813	}
814
815
816	void compare_tree(char tree1, char tree2, sint *hits, sint n)
817	{
818	sint i,j,k;
819	sint nhits1, nhits2;
820
821	for(i=1; i<=n-3; i++) {
822	for(j=1; j<=n-3; j++) {
823	nhits1 = 0;
824	nhits2 = 0;
825	for(k=1; k<=n; k++) {
826	if(tree1[i][k] == tree2[j][k]) nhits1++;
827	if(tree1[i][k] != tree2[j][k]) nhits2++;
828	}
829	if((nhits1 == last_seq-first_seq+1) \|\| (nhits2 == last_seq-first_seq+1)) hits[i]++;
830	}
831	}
832	}
833
834
835	void print_phylip_tree(char *tree_description, FILE tree, sint bootstrap)
836	{
837	sint old_row;
838
839	fprintf(tree,"(\n");
840
841	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,1,bootstrap);
842	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-2]);
843	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
844	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
845	fprintf(tree,",\n");
846
847	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,2,bootstrap);
848	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-1]);
849	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
850	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
851	fprintf(tree,",\n");
852
853	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,3,bootstrap);
854	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1]);
855	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
856	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
857	fprintf(tree,")");
858	if (bootstrap==BS_NODE_LABELS) fprintf(tree,"TRICHOTOMY");
859	fprintf(tree,";\n");
860	}
861
862
863	sint two_way_split
864	(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap)
865	{
866	sint row, new_row = 0, old_row, col, test_col = 0;
867	Boolean single_seq;
868
869	if(start_row != last_seq-first_seq+1-2) fprintf(tree,"(\n");
870
871	for(col=1; col<=last_seq-first_seq+1; col++) {
872	if(tree_description[start_row][col] == flag) {
873	test_col = col;
874	break;
875	}
876	}
877
878	single_seq = TRUE;
879	for(row=start_row-1; row>=1; row--)
880	if(tree_description[row][test_col] == 1) {
881	single_seq = FALSE;
882	new_row = row;
883	break;
884	}
885
886	if(single_seq) {
887	tree_description[start_row][test_col] = 0;
888	fprintf(tree,"%.*s",max_names,names[test_col+first_seq-1]);
889	if(start_row == last_seq-first_seq+1-2) {
890	return(0);
891	}
892
893	fprintf(tree,":%7.5f,\n",left_branch[start_row]);
894	}
895	else {
896	for(col=1; col<=last_seq-first_seq+1; col++) {
897	if((tree_description[start_row][col]==1)&&
898	(tree_description[new_row][col]==1))
899	tree_description[start_row][col] = 0;
900	}
901	old_row=two_way_split(tree_description, tree, new_row, (sint)1, bootstrap);
902	if(start_row == last_seq-first_seq+1-2) {
903	return(new_row);
904	}
905
906	fprintf(tree,":%7.5f",left_branch[start_row]);
907	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
908	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
909
910	fprintf(tree,",\n");
911	}
912
913
914	for(col=1; col<=last_seq-first_seq+1; col++)
915	if(tree_description[start_row][col] == flag) {
916	test_col = col;
917	break;
918	}
919
920	single_seq = TRUE;
921	new_row = 0;
922	for(row=start_row-1; row>=1; row--)
923	if(tree_description[row][test_col] == 1) {
924	single_seq = FALSE;
925	new_row = row;
926	break;
927	}
928
929	if(single_seq) {
930	tree_description[start_row][test_col] = 0;
931	fprintf(tree,"%.*s",max_names,names[test_col+first_seq-1]);
932	fprintf(tree,":%7.5f)\n",right_branch[start_row]);
933	}
934	else {
935	for(col=1; col<=last_seq-first_seq+1; col++) {
936	if((tree_description[start_row][col]==1)&&
937	(tree_description[new_row][col]==1))
938	tree_description[start_row][col] = 0;
939	}
940	old_row=two_way_split(tree_description, tree, new_row, (sint)1, bootstrap);
941	fprintf(tree,":%7.5f",right_branch[start_row]);
942	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
943	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
944
945	fprintf(tree,")\n");
946	}
947	if ((bootstrap==BS_NODE_LABELS) && (boot_totals[start_row]>0))
948	fprintf(tree,"%d",(pint)boot_totals[start_row]);
949
950	return(start_row);
951	}
952
953
954
955	void print_tree(char *tree_description, FILE tree, sint *totals)
956	{
957	sint row,col;
958
959	fprintf(tree,"\n");
960
961	for(row=1; row<=last_seq-first_seq+1-3; row++) {
962	fprintf(tree," \n");
963	for(col=1; col<=last_seq-first_seq+1; col++) {
964	if(tree_description[row][col] == 0)
965	fprintf(tree,"*");
966	else
967	fprintf(tree,".");
968	}
969	if(totals[row] > 0)
970	fprintf(tree,"%7d",(pint)totals[row]);
971	}
972	fprintf(tree," \n");
973	for(col=1; col<=last_seq-first_seq+1; col++)
974	fprintf(tree,"%1d",(pint)tree_description[last_seq-first_seq+1-2][col]);
975	fprintf(tree,"\n");
976	}
977
978
979
980	sint dna_distance_matrix(FILE *tree)
981	{
982	sint m,n;
983	sint j,i;
984	sint res1, res2;
985	sint overspill = 0;
986	double p,q,e,a,b,k;
987
988	tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */
989
990	if(verbose) {
991	fprintf(tree,"\n");
992	fprintf(tree,"\n DIST = percentage divergence (/100)");
993	fprintf(tree,"\n p = rate of transition (A <-> G; C <-> T)");
994	fprintf(tree,"\n q = rate of transversion");
995	fprintf(tree,"\n Length = number of sites used in comparison");
996	fprintf(tree,"\n");
997	if(tossgaps) {
998	fprintf(tree,"\n All sites with gaps (in any sequence) deleted!");
999	fprintf(tree,"\n");
1000	}
1001	if(kimura) {
1002	fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:");
1003	fprintf(tree,"\n\n Kimura, M. (1980)");
1004	fprintf(tree," A simple method for estimating evolutionary ");
1005	fprintf(tree,"rates of base");
1006	fprintf(tree,"\n substitutions through comparative studies of ");
1007	fprintf(tree,"nucleotide sequences.");
1008	fprintf(tree,"\n J. Mol. Evol., 16, 111-120.");
1009	fprintf(tree,"\n\n");
1010	}
1011	}
1012
1013	for(m=1; m<last_seq-first_seq+1; ++m) /* for every pair of sequence */
1014	for(n=m+1; n<=last_seq-first_seq+1; ++n) {
1015	p = q = e = 0.0;
1016	tmat[m][n] = tmat[n][m] = 0.0;
1017	for(i=1; i<=seqlen_array[first_seq]; ++i) {
1018	j = boot_positions[i];
1019	if(tossgaps && (tree_gaps[j] > 0) )
1020	goto skip; /* gap position */
1021	res1 = seq_array[m+first_seq-1][j];
1022	res2 = seq_array[n+first_seq-1][j];
1023	if( (res1 == gap_pos1) \|\| (res1 == gap_pos2) \|\|
1024	(res2 == gap_pos1) \|\| (res2 == gap_pos2))
1025	goto skip; /* gap in a seq*/
1026	e = e + 1.0;
1027	if(res1 != res2) {
1028	if(transition(res1,res2))
1029	p = p + 1.0;
1030	else
1031	q = q + 1.0;
1032	}
1033	skip:;
1034	}
1035
1036
1037	/* Kimura's 2 parameter correction for multiple substitutions */
1038
1039	if(!kimura) {
1040	if (e == 0) {
1041	fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
1042	k = 0.0;
1043	p = 0.0;
1044	q = 0.0;
1045	}
1046	else {
1047	k = (p+q)/e;
1048	if(p > 0.0)
1049	p = p/e;
1050	else
1051	p = 0.0;
1052	if(q > 0.0)
1053	q = q/e;
1054	else
1055	q = 0.0;
1056	}
1057	tmat[m][n] = tmat[n][m] = k;
1058	if(verbose) /* if screen output */
1059	fprintf(tree,
1060	"%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
1061	,(pint)m,(pint)n,k,p,q,e);
1062	}
1063	else {
1064	if (e == 0) {
1065	fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
1066	p = 0.0;
1067	q = 0.0;
1068	}
1069	else {
1070	if(p > 0.0)
1071	p = p/e;
1072	else
1073	p = 0.0;
1074	if(q > 0.0)
1075	q = q/e;
1076	else
1077	q = 0.0;
1078	}
1079
1080	if( ((2.0*p)+q) == 1.0 )
1081	a = 0.0;
1082	else
1083	a = 1.0/(1.0-(2.0*p)-q);
1084
1085	if( q == 0.5 )
1086	b = 0.0;
1087	else
1088	b = 1.0/(1.0-(2.0*q));
1089
1090	/* watch for values going off the scale for the correction. */
1091	if( (a<=0.0) \|\| (b<=0.0) ) {
1092	overspill++;
1093	k = 3.5; /* arbitrary high score */
1094	}
1095	else
1096	k = 0.5log(a) + 0.25log(b);
1097	tmat[m][n] = tmat[n][m] = k;
1098	if(verbose) /* if screen output */
1099	fprintf(tree,
1100	"%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
1101	,(pint)m,(pint)n,k,p,q,e);
1102
1103	}
1104	}
1105	return overspill; /* return the number of off-scale values */
1106	}
1107
1108
1109	sint prot_distance_matrix(FILE *tree)
1110	{
1111	sint m,n;
1112	sint j,i;
1113	sint res1, res2;
1114	sint overspill = 0;
1115	double p,e,k, table_entry;
1116
1117
1118	tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */
1119
1120	if(verbose) {
1121	fprintf(tree,"\n");
1122	fprintf(tree,"\n DIST = percentage divergence (/100)");
1123	fprintf(tree,"\n Length = number of sites used in comparison");
1124	fprintf(tree,"\n\n");
1125	if(tossgaps) {
1126	fprintf(tree,"\n All sites with gaps (in any sequence) deleted");
1127	fprintf(tree,"\n");
1128	}
1129	if(kimura) {
1130	fprintf(tree,"\n Distances up tp 0.75 corrected by Kimura's empirical method:");
1131	fprintf(tree,"\n\n Kimura, M. (1983)");
1132	fprintf(tree," The Neutral Theory of Molecular Evolution.");
1133	fprintf(tree,"\n Page 75. Cambridge University Press, Cambridge, England.");
1134	fprintf(tree,"\n\n");
1135	}
1136	}
1137
1138	for(m=1; m<nseqs; ++m) /* for every pair of sequence */
1139	for(n=m+1; n<=nseqs; ++n) {
1140	p = e = 0.0;
1141	tmat[m][n] = tmat[n][m] = 0.0;
1142	for(i=1; i<=seqlen_array[1]; ++i) {
1143	j = boot_positions[i];
1144	if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */
1145	res1 = seq_array[m][j];
1146	res2 = seq_array[n][j];
1147	if( (res1 == gap_pos1) \|\| (res1 == gap_pos2) \|\|
1148	(res2 == gap_pos1) \|\| (res2 == gap_pos2))
1149	goto skip; /* gap in a seq*/
1150	e = e + 1.0;
1151	if(res1 != res2) p = p + 1.0;
1152	skip:;
1153	}
1154
1155	if(p <= 0.0)
1156	k = 0.0;
1157	else
1158	k = p/e;
1159
1160	/* DES debug */
1161	/* fprintf(stdout,"Seq1=%4d Seq2=%4d k =%7.4f \n",(pint)m,(pint)n,k); */
1162	/* DES debug */
1163
1164	if(kimura) {
1165	if(k < 0.75) { /* use Kimura's formula */
1166	if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) );
1167	}
1168	else {
1169	if(k > 0.930) {
1170	overspill++;
1171	k = 10.0; /* arbitrarily set to 1000% */
1172	}
1173	else {
1174	table_entry = (k*1000.0) - 750.0;
1175	k = (double)dayhoff_pams[(int)table_entry];
1176	k = k/100.0;
1177	}
1178	}
1179	}
1180
1181	tmat[m][n] = tmat[n][m] = k;
1182	if(verbose) /* if screen output */
1183	fprintf(tree,
1184	"%4d vs.%4d DIST = %6.4f; length = %6.0f\n",
1185	(pint)m,(pint)n,k,e);
1186	}
1187	return overspill;
1188	}
1189
1190
1191	void guide_tree(FILE *tree,sint firstseq,sint numseqs)
1192	/*
1193	Routine for producing unrooted NJ trees from seperately aligned
1194	pairwise distances. This produces the GUIDE DENDROGRAMS in
1195	PHYLIP format.
1196	*/
1197	{
1198	static char **standard_tree;
1199	sint i;
1200
1201	phylip_phy_tree_file=tree;
1202	verbose = FALSE;
1203	first_seq=firstseq;
1204	last_seq=first_seq+numseqs-1;
1205
1206	standard_tree = (char *) ckalloc( (last_seq-first_seq+2) sizeof (char *) );
1207	for(i=0; i<last_seq-first_seq+2; i++)
1208	standard_tree[i] = (char ) ckalloc( (last_seq-first_seq+2) sizeof(char));
1209
1210	nj_tree(standard_tree,clustal_phy_tree_file);
1211
1212	print_phylip_tree(standard_tree,phylip_phy_tree_file,0);
1213
1214	left_branch=ckfree((void *)left_branch);
1215	right_branch=ckfree((void *)right_branch);
1216	tkill=ckfree((void *)tkill);
1217	av=ckfree((void *)av);
1218	for (i=1;i<last_seq-first_seq+2;i++)
1219	standard_tree[i]=ckfree((void *)standard_tree[i]);
1220	standard_tree=ckfree((void *)standard_tree);
1221
1222	fclose(phylip_phy_tree_file);
1223
1224	}
1225

Note: See TracBrowser for help on using the repository browser.

Download in other formats: