Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

trees.c

Visit:

Last change on this file was 19480, checked in by westram, 2 years ago
reintegrates 'gde' into 'trunk' adds: log:branches/gde@19460:19479
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 62.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	void print_nexus_tree(char *tree_description, FILE tree, sint bootstrap);
22	sint two_way_split(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap);
23	sint two_way_split_nexus(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap);
24	void print_tree(char *tree_description, FILE tree, sint *totals);
25	static Boolean is_ambiguity(char c);
26	static void overspill_message(sint overspill,sint total_dists);
27
28
29	/*
30	* Global variables
31	*/
32
33	extern sint max_names;
34
35	extern double *tmat; / general nxn array of reals; allocated from main */
36	/* this is used as a distance matrix */
37	extern Boolean dnaflag; /* TRUE for DNA seqs; FALSE for proteins */
38	extern Boolean tossgaps; /* Ignore places in align. where ANY seq. has a gap*/
39	extern Boolean kimura; /* Use correction for multiple substitutions */
40	extern Boolean output_tree_clustal; /* clustal text output for trees */
41	extern Boolean output_tree_phylip; /* phylip nested parentheses format */
42	extern Boolean output_tree_distances; /* phylip distance matrix */
43	extern Boolean output_tree_nexus; /* nexus format tree */
44	extern Boolean output_pim; /* perc identity matrix output Ramu */
45
46	extern sint bootstrap_format; /* bootstrap file format */
47	extern Boolean empty; /* any sequences in memory? */
48	extern Boolean usemenu; /* interactive (TRUE) or command line (FALSE) */
49	extern sint nseqs;
50	extern sint max_aln_length;
51	extern sint seqlen_array; / the lengths of the sequences */
52	extern char *seq_array; / the sequences */
53	extern char *names; / the seq. names */
54	extern char seqname[]; /* name of input file */
55	extern sint gap_pos1,gap_pos2;
56	extern Boolean use_ambiguities;
57	extern char *amino_acid_codes;
58
59	static double *av;
60	static double left_branch, right_branch;
61	static double save_left_branch, save_right_branch;
62	static sint *boot_totals;
63	static sint *tkill;
64	/*
65	The next line is a fossil from the days of using the cc ran()
66	static int ran_factor;
67	*/
68	static sint *boot_positions;
69	static FILE *phylip_phy_tree_file;
70	static FILE *clustal_phy_tree_file;
71	static FILE *distances_phy_tree_file;
72	static FILE *nexus_phy_tree_file;
73	static FILE pim_file; / Ramu */
74	static Boolean verbose;
75	static char *tree_gaps;
76	static sint first_seq, last_seq;
77	/* array of weights; 1 for use this posn.; 0 don't */
78
79	extern sint boot_ntrials; /* number of bootstrap trials */
80	extern unsigned sint boot_ran_seed; /* random number generator seed */
81
82	void phylogenetic_tree(char phylip_name, char clustal_name, char dist_name, char nexus_name, char *pim_name)
83	/*
84	Calculate a tree using the distances in the nseqs*nseqs array tmat.
85	This is the routine for getting the REAL trees after alignment.
86	*/
87	{ char path[FILENAMELEN+1];
88	sint i, j;
89	sint overspill = 0;
90	sint total_dists;
91	static char **standard_tree;
92	static char **save_tree;
93
94	if(empty) {
95	error("You must load an alignment first");
96	return;
97	}
98
99	if(nseqs<2) {
100	error("Alignment has only %d sequences",nseqs);
101	return;
102	}
103	first_seq=1;
104	last_seq=nseqs;
105
106	get_path(seqname,path);
107
108	if(output_tree_clustal) {
109	if (clustal_name[0]!=EOS) {
110	if((clustal_phy_tree_file = open_explicit_file(
111	clustal_name))==NULL) return;
112	}
113	else {
114	if((clustal_phy_tree_file = open_output_file(
115	"\nEnter name for CLUSTAL tree output file ",path,
116	clustal_name,"nj")) == NULL) return;
117	}
118	}
119
120	if(output_tree_phylip) {
121	if (phylip_name[0]!=EOS) {
122	if((phylip_phy_tree_file = open_explicit_file(
123	phylip_name))==NULL) return;
124	}
125	else {
126	if((phylip_phy_tree_file = open_output_file(
127	"\nEnter name for PHYLIP tree output file ",path,
128	phylip_name,"ph")) == NULL) return;
129	}
130	}
131
132	if(output_tree_distances)
133	{
134	if (dist_name[0]!=EOS) {
135	if((distances_phy_tree_file = open_explicit_file(
136	dist_name))==NULL) return;
137	}
138	else {
139	if((distances_phy_tree_file = open_output_file(
140	"\nEnter name for distance matrix output file ",path,
141	dist_name,"dst")) == NULL) return;
142	}
143	}
144
145	if(output_tree_nexus)
146	{
147	if (nexus_name[0]!=EOS) {
148	if((nexus_phy_tree_file = open_explicit_file(
149	nexus_name))==NULL) return;
150	}
151	else {
152	if((nexus_phy_tree_file = open_output_file(
153	"\nEnter name for NEXUS tree output file ",path,
154	nexus_name,"tre")) == NULL) return;
155	}
156	}
157
158	if(output_pim)
159	{
160	if (pim_name[0]!=EOS) {
161	if((pim_file = open_explicit_file(
162	pim_name))==NULL) return;
163	}
164	else {
165	if((pim_file = open_output_file(
166	"\nEnter name for % Identity matrix output file ",path,
167	pim_name,"pim")) == NULL) return;
168	}
169	}
170
171	boot_positions = (sint )ckalloc( (seqlen_array[first_seq]+2) sizeof (sint) );
172
173	for(j=1; j<=seqlen_array[first_seq]; ++j)
174	boot_positions[j] = j;
175
176	if(output_tree_clustal) {
177	verbose = TRUE; /* Turn on file output */
178	if(dnaflag)
179	overspill = dna_distance_matrix(clustal_phy_tree_file);
180	else
181	overspill = prot_distance_matrix(clustal_phy_tree_file);
182	}
183
184	if(output_tree_phylip) {
185	verbose = FALSE; /* Turn off file output */
186	if(dnaflag)
187	overspill = dna_distance_matrix(phylip_phy_tree_file);
188	else
189	overspill = prot_distance_matrix(phylip_phy_tree_file);
190	}
191
192	if(output_tree_nexus) {
193	verbose = FALSE; /* Turn off file output */
194	if(dnaflag)
195	overspill = dna_distance_matrix(nexus_phy_tree_file);
196	else
197	overspill = prot_distance_matrix(nexus_phy_tree_file);
198	}
199
200	if(output_pim) { /* Ramu */
201	verbose = FALSE; /* Turn off file output */
202	if(dnaflag)
203	calc_percidentity(pim_file);
204	else
205	calc_percidentity(pim_file);
206	}
207
208
209	if(output_tree_distances) {
210	verbose = FALSE; /* Turn off file output */
211	if(dnaflag)
212	overspill = dna_distance_matrix(distances_phy_tree_file);
213	else
214	overspill = prot_distance_matrix(distances_phy_tree_file);
215	distance_matrix_output(distances_phy_tree_file);
216	}
217
218	/* check if any distances overflowed the distance corrections */
219	if ( overspill > 0 ) {
220	total_dists = (nseqs*(nseqs-1))/2;
221	overspill_message(overspill,total_dists);
222	}
223
224	if(output_tree_clustal) verbose = TRUE; /* Turn on file output */
225
226	standard_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
227	for(i=0; i<nseqs+1; i++)
228	standard_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
229	save_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
230	for(i=0; i<nseqs+1; i++)
231	save_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
232
233	if(output_tree_clustal \|\| output_tree_phylip \|\| output_tree_nexus)
234	nj_tree(standard_tree,clustal_phy_tree_file);
235
236	for(i=1; i<nseqs+1; i++)
237	for(j=1; j<nseqs+1; j++)
238	save_tree[i][j] = standard_tree[i][j];
239
240	if(output_tree_phylip)
241	print_phylip_tree(standard_tree,phylip_phy_tree_file,0);
242
243	for(i=1; i<nseqs+1; i++)
244	for(j=1; j<nseqs+1; j++)
245	standard_tree[i][j] = save_tree[i][j];
246
247	if(output_tree_nexus)
248	print_nexus_tree(standard_tree,nexus_phy_tree_file,0);
249
250	/*
251	print_tree(standard_tree,phy_tree_file);
252	*/
253	tree_gaps=ckfree((void *)tree_gaps);
254	boot_positions=ckfree((void *)boot_positions);
255	if (left_branch != NULL) left_branch=ckfree((void *)left_branch);
256	if (right_branch != NULL) right_branch=ckfree((void *)right_branch);
257	if (tkill != NULL) tkill=ckfree((void *)tkill);
258	if (av != NULL) av=ckfree((void *)av);
259	for (i=0;i<nseqs+1;i++)
260	standard_tree[i]=ckfree((void *)standard_tree[i]);
261	standard_tree=ckfree((void *)standard_tree);
262
263	for (i=0;i<nseqs+1;i++)
264	save_tree[i]=ckfree((void *)save_tree[i]);
265	save_tree=ckfree((void *)save_tree);
266
267	if(output_tree_clustal) {
268	fclose(clustal_phy_tree_file);
269	info("Phylogenetic tree file created: [%s]",clustal_name);
270	}
271
272	if(output_tree_phylip) {
273	fclose(phylip_phy_tree_file);
274	info("Phylogenetic tree file created: [%s]",phylip_name);
275	}
276
277	if(output_tree_distances) {
278	fclose(distances_phy_tree_file);
279	info("Distance matrix file created: [%s]",dist_name);
280	}
281
282	if(output_tree_nexus) {
283	fclose(nexus_phy_tree_file);
284	info("Nexus tree file created: [%s]",nexus_name);
285	}
286
287	if(output_pim) {
288	fclose(pim_file);
289	info(" perc identity matrix file created: [%s]",pim_name);
290	}
291
292	}
293
294	static void overspill_message(sint overspill,sint total_dists)
295	{
296	char err_mess[1024]="";
297
298	sprintf(err_mess,"%d of the distances out of a total of %d",
299	(pint)overspill,(pint)total_dists);
300	strcat(err_mess,"\n were out of range for the distance correction.");
301	strcat(err_mess,"\n");
302	strcat(err_mess,"\n SUGGESTIONS: 1) remove the most distant sequences");
303	strcat(err_mess,"\n or 2) use the PHYLIP package");
304	strcat(err_mess,"\n or 3) turn off the correction.");
305	strcat(err_mess,"\n Note: Use option 3 with caution! With this degree");
306	strcat(err_mess,"\n of divergence you will have great difficulty");
307	strcat(err_mess,"\n getting robust and reliable trees.");
308	strcat(err_mess,"\n\n");
309	warning(err_mess);
310	}
311
312
313
314	Boolean transition(sint base1, sint base2) /* TRUE if transition; else FALSE */
315	/*
316
317	assumes that the bases of DNA sequences have been translated as
318	a,A = 0; c,C = 1; g,G = 2; t,T,u,U = 3; N = 4;
319	a,A = 0; c,C = 2; g,G = 6; t,T,u,U =17;
320
321	A <--> G and T <--> C are transitions; all others are transversions.
322
323	*/
324	{
325	if( ((base1 == 0) && (base2 == 6)) \|\| ((base1 == 6) && (base2 == 0)) )
326	return TRUE; /* A <--> G */
327	if( ((base1 ==17) && (base2 == 2)) \|\| ((base1 == 2) && (base2 ==17)) )
328	return TRUE; /* T <--> C */
329	return FALSE;
330	}
331
332
333	void tree_gap_delete(void) /* flag all positions in alignment that have a gap */
334	{ /* in ANY sequence */
335	sint seqn;
336	sint posn;
337
338	tree_gaps = (char )ckalloc( (max_aln_length+1) sizeof (char) );
339
340	for(posn=1; posn<=seqlen_array[first_seq]; ++posn) {
341	tree_gaps[posn] = 0;
342	for(seqn=1; seqn<=last_seq-first_seq+1; ++seqn) {
343	if((seq_array[seqn+first_seq-1][posn] == gap_pos1) \|\|
344	(seq_array[seqn+first_seq-1][posn] == gap_pos2)) {
345	tree_gaps[posn] = 1;
346	break;
347	}
348	}
349	}
350
351	}
352
353	void distance_matrix_output(FILE *ofile)
354	{
355	sint i,j;
356
357	fprintf(ofile,"%6d",(pint)last_seq-first_seq+1);
358	for(i=1;i<=last_seq-first_seq+1;i++) {
359	fprintf(ofile,"\n%-*s ",max_names,names[i]);
360	for(j=1;j<=last_seq-first_seq+1;j++) {
361	fprintf(ofile,"%6.3f ",tmat[i][j]);
362	if(j % 8 == 0) {
363	if(j!=last_seq-first_seq+1) fprintf(ofile,"\n");
364	if(j != last_seq-first_seq+1 ) fprintf(ofile," ");
365	}
366	}
367	}
368	}
369
370
371
372	#ifdef ORIGINAL_NJ_TREE
373	void nj_tree(char *tree_description, FILE tree)
374	{
375	register int i;
376	sint l[4],nude,k;
377	sint nc,mini,minj,j,ii,jj;
378	double fnseqs,fnseqs2=0,sumd;
379	double diq,djq,dij,d2r,dr,dio,djo,da;
380	double tmin,total,dmin;
381	double bi,bj,b1,b2,b3,branch[4];
382	sint typei,typej; /* 0 = node; 1 = OTU */
383
384	fnseqs = (double)last_seq-first_seq+1;
385
386	/********************* First initialisation *************************/
387
388	if(verbose) {
389	fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
390	fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
391	fprintf(tree," The Neighbor-joining Method:");
392	fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
393	fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
394	fprintf(tree,"\n\n This is an UNROOTED tree\n");
395	fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
396	}
397
398	if (fnseqs == 2) {
399	if (verbose) fprintf(tree,"Cycle 1 = SEQ: 1 (%9.5f) joins SEQ: 2 (%9.5f)",tmat[first_seq][first_seq+1],tmat[first_seq][first_seq+1]);
400	return;
401	}
402
403	mini = minj = 0;
404
405	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
406	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
407	tkill = (sint ) ckalloc( (nseqs+1) sizeof (sint) );
408	av = (double ) ckalloc( (nseqs+1) sizeof (double) );
409
410	for(i=1;i<=last_seq-first_seq+1;++i)
411	{
412	tmat[i][i] = av[i] = 0.0;
413	tkill[i] = 0;
414	}
415
416	/********************* Enter The Main Cycle *************************/
417
418	/* for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) { / /start main cycle*/
419	for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) {
420	sumd = 0.0;
421	for(j=2; j<=last_seq-first_seq+1; ++j)
422	for(i=1; i<j; ++i) {
423	tmat[j][i] = tmat[i][j];
424	sumd = sumd + tmat[i][j];
425	}
426
427	tmin = 99999.0;
428
429	/.................compute SMATij values and find the smallest one ......../
430
431	for(jj=2; jj<=last_seq-first_seq+1; ++jj)
432	if(tkill[jj] != 1)
433	for(ii=1; ii<jj; ++ii)
434	if(tkill[ii] != 1) {
435	diq = djq = 0.0;
436
437	for(i=1; i<=last_seq-first_seq+1; ++i) {
438	diq = diq + tmat[i][ii];
439	djq = djq + tmat[i][jj];
440	}
441
442	dij = tmat[ii][jj];
443	d2r = diq + djq - (2.0*dij);
444	dr = sumd - dij -d2r;
445	fnseqs2 = fnseqs - 2.0;
446	total= d2r+ fnseqs2dij +dr2.0;
447	total= total / (2.0*fnseqs2);
448
449	if(total < tmin) {
450	tmin = total;
451	mini = ii;
452	minj = jj;
453	}
454	}
455
456
457	/.................compute branch lengths and print the results ......../
458
459
460	dio = djo = 0.0;
461	for(i=1; i<=last_seq-first_seq+1; ++i) {
462	dio = dio + tmat[i][mini];
463	djo = djo + tmat[i][minj];
464	}
465
466	dmin = tmat[mini][minj];
467	dio = (dio - dmin) / fnseqs2;
468	djo = (djo - dmin) / fnseqs2;
469	bi = (dmin + dio - djo) * 0.5;
470	bj = dmin - bi;
471	bi = bi - av[mini];
472	bj = bj - av[minj];
473
474	if( av[mini] > 0.0 )
475	typei = 0;
476	else
477	typei = 1;
478	if( av[minj] > 0.0 )
479	typej = 0;
480	else
481	typej = 1;
482
483	if(verbose)
484	fprintf(tree,"\n Cycle%4d = ",(pint)nc);
485
486	/*
487	set negative branch lengths to zero. Also set any tiny positive
488	branch lengths to zero.
489	*/ if( fabs(bi) < 0.0001) bi = 0.0;
490	if( fabs(bj) < 0.0001) bj = 0.0;
491
492	if(verbose) {
493	if(typei == 0)
494	fprintf(tree,"Node:%4d (%9.5f) joins ",(pint)mini,bi);
495	else
496	fprintf(tree," SEQ:%4d (%9.5f) joins ",(pint)mini,bi);
497
498	if(typej == 0)
499	fprintf(tree,"Node:%4d (%9.5f)",(pint)minj,bj);
500	else
501	fprintf(tree," SEQ:%4d (%9.5f)",(pint)minj,bj);
502
503	fprintf(tree,"\n");
504	}
505
506
507	left_branch[nc] = bi;
508	right_branch[nc] = bj;
509
510	for(i=1; i<=last_seq-first_seq+1; i++)
511	tree_description[nc][i] = 0;
512
513	if(typei == 0) {
514	for(i=nc-1; i>=1; i--)
515	if(tree_description[i][mini] == 1) {
516	for(j=1; j<=last_seq-first_seq+1; j++)
517	if(tree_description[i][j] == 1)
518	tree_description[nc][j] = 1;
519	break;
520	}
521	}
522	else
523	tree_description[nc][mini] = 1;
524
525	if(typej == 0) {
526	for(i=nc-1; i>=1; i--)
527	if(tree_description[i][minj] == 1) {
528	for(j=1; j<=last_seq-first_seq+1; j++)
529	if(tree_description[i][j] == 1)
530	tree_description[nc][j] = 1;
531	break;
532	}
533	}
534	else
535	tree_description[nc][minj] = 1;
536
537
538	/*
539	Here is where the -0.00005 branch lengths come from for 3 or more
540	identical seqs.
541	*/
542	/* if(dmin <= 0.0) dmin = 0.0001; */
543	if(dmin <= 0.0) dmin = 0.000001;
544	av[mini] = dmin * 0.5;
545
546	/........................Re-initialisation................................/
547
548	fnseqs = fnseqs - 1.0;
549	tkill[minj] = 1;
550
551	for(j=1; j<=last_seq-first_seq+1; ++j)
552	if( tkill[j] != 1 ) {
553	da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
554	if( (mini - j) < 0 )
555	tmat[mini][j] = da;
556	if( (mini - j) > 0)
557	tmat[j][mini] = da;
558	}
559
560	for(j=1; j<=last_seq-first_seq+1; ++j)
561	tmat[minj][j] = tmat[j][minj] = 0.0;
562
563
564	/**/ } /end main cycle**/
565
566	/****************************Last Cycle (3 Seqs. left)******************/
567
568	nude = 1;
569
570	for(i=1; i<=last_seq-first_seq+1; ++i)
571	if( tkill[i] != 1 ) {
572	l[nude] = i;
573	nude = nude + 1;
574	}
575
576	b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
577	b2 = tmat[l[1]][l[2]] - b1;
578	b3 = tmat[l[1]][l[3]] - b1;
579
580	branch[1] = b1 - av[l[1]];
581	branch[2] = b2 - av[l[2]];
582	branch[3] = b3 - av[l[3]];
583
584	/* Reset tiny negative and positive branch lengths to zero */
585	if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
586	if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
587	if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
588
589	left_branch[last_seq-first_seq+1-2] = branch[1];
590	left_branch[last_seq-first_seq+1-1] = branch[2];
591	left_branch[last_seq-first_seq+1] = branch[3];
592
593	for(i=1; i<=last_seq-first_seq+1; i++)
594	tree_description[last_seq-first_seq+1-2][i] = 0;
595
596	if(verbose)
597	fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",(pint)nc);
598
599	for(i=1; i<=3; ++i) {
600	if( av[l[i]] > 0.0) {
601	if(verbose)
602	fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",(pint)l[i],branch[i]);
603	for(k=last_seq-first_seq+1-3; k>=1; k--)
604	if(tree_description[k][l[i]] == 1) {
605	for(j=1; j<=last_seq-first_seq+1; j++)
606	if(tree_description[k][j] == 1)
607	tree_description[last_seq-first_seq+1-2][j] = i;
608	break;
609	}
610	}
611	else {
612	if(verbose)
613	fprintf(tree,"\n\t\t SEQ:%4d (%9.5f) ",(pint)l[i],branch[i]);
614	tree_description[last_seq-first_seq+1-2][l[i]] = i;
615	}
616	if(i < 3) {
617	if(verbose)
618	fprintf(tree,"joins");
619	}
620	}
621
622	if(verbose)
623	fprintf(tree,"\n");
624
625	}
626
627	#else /* ORIGINAL_NJ_TREE */
628
629	void nj_tree(char *tree_description, FILE tree) {
630	void fast_nj_tree();
631
632	/fprintf(stderr, "*** call fast_nj_tree() !!!! ***\n");/
633	fast_nj_tree(tree_description, tree);
634	}
635
636
637	/****************************************************************************
638	* [ Improvement ideas in fast_nj_tree() ] by DDBJ & FUJITSU Limited.
639	* written by Tadashi Koike
640	* (takoike@genes.nig.ac.jp)
641	*******************
642	* <IMPROVEMENT 1> : Store the value of sum of the score to temporary array,
643	* and use again and again.
644	*
645	* In the main cycle, these are calculated again and again :
646	* diq = sum of tmat[n][ii] (n:1 to last_seq-first_seq+1),
647	* djq = sum of tmat[n][jj] (n:1 to last_seq-first_seq+1),
648	* dio = sum of tmat[n][mini] (n:1 to last_seq-first_seq+1),
649	* djq = sum of tmat[n][minj] (n:1 to last_seq-first_seq+1)
650	* // 'last_seq' and 'first_seq' are both constant values //
651	* and the result of above calculations is always same until
652	* a best pair of neighbour nodes is joined.
653	*
654	* So, we change the logic to calculate the sum[i] (=sum of tmat[n][i]
655	* (n:1 to last_seq-first_seq+1)) and store it to array, before
656	* beginning to find a best pair of neighbour nodes, and after that
657	* we use them again and again.
658	*
659	* tmat[i][j]
660	* 1 2 3 4 5
661	* +---+---+---+---+---+
662	* 1 \| \| \| \| \| \|
663	* +---+---+---+---+---+
664	* 2 \| \| \| \| \| \| 1) calculate sum of tmat[n][i]
665	* +---+---+---+---+---+ (n: 1 to last_seq-first_seq+1)
666	* 3 \| \| \| \| \| \| 2) store that sum value to sum[i]
667	* +---+---+---+---+---+
668	* 4 \| \| \| \| \| \| 3) use sum[i] during finding a best
669	* +---+---+---+---+---+ pair of neibour nodes.
670	* 5 \| \| \| \| \| \|
671	* +---+---+---+---+---+
672	* \| \| \| \| \|
673	* V V V V V Calculate sum , and store it to sum[i]
674	* +---+---+---+---+---+
675	* sum[i] \| \| \| \| \| \|
676	* +---+---+---+---+---+
677	*
678	* At this time, we thought that we use upper triangle of the matrix
679	* because tmat[i][j] is equal to tmat[j][i] and tmat[i][i] is equal
680	* to zero. Therefore, we prepared sum_rows[i] and sum_cols[i] instead
681	* of sum[i] for storing the sum value.
682	*
683	* tmat[i][j]
684	* 1 2 3 4 5 sum_cols[i]
685	* +---+---+---+---+---+ +---+
686	* 1 \| # \| # \| # \| # \| --> \| \| ... sum of tmat[1][2..5]
687	* + - +---+---+---+---+ +---+
688	* 2 \| # \| # \| # \| --> \| \| ... sum of tmat[2][3..5]
689	* + - + - +---+---+---+ +---+
690	* 3 \| # \| # \| --> \| \| ... sum of tmat[3][4..5]
691	* + - + - + - +---+---+ +---+
692	* 4 \| # \| --> \| \| ... sum of tmat[4][5]
693	* + - + - + - + - +---+ +---+
694	* 5 \| --> \| \| ... zero
695	* + - + - + - + - + - + +---+
696	* \| \| \| \| \|
697	* V V V V V Calculate sum , sotre to sum[i]
698	* +---+---+---+---+---+
699	* sum_rows[i] \| \| \| \| \| \|
700	* +---+---+---+---+---+
701	* \| \| \| \| \|
702	* \| \| \| \| +----- sum of tmat[1..4][5]
703	* \| \| \| +--------- sum of tmat[1..3][4]
704	* \| \| +------------- sum of tmat[1..2][3]
705	* \| +----------------- sum of tmat[1][2]
706	* +--------------------- zero
707	*
708	* And we use (sum_rows[i] + sum_cols[i]) instead of sum[i].
709	*
710	*******************
711	* <IMPROVEMENT 2> : We manage valid nodes with chain list, instead of
712	* tkill[i] flag array.
713	*
714	* In original logic, invalid(killed?) nodes after nodes-joining
715	* are managed with tkill[i] flag array (set to 1 when killed).
716	* By this method, it is conspicuous to try next node but skip it
717	* at the latter of finding a best pair of neighbor nodes.
718	*
719	* So, we thought that we managed valid nodes by using a chain list
720	* as below:
721	*
722	* 1) declare the list structure.
723	* struct {
724	* sint n; // entry number of node.
725	* void *prev; // pointer to previous entry.
726	* void *next; // pointer to next entry.
727	* }
728	* 2) construct a valid node list.
729	*
730	* +-----+ +-----+ +-----+ +-----+ +-----+
731	* NULL<-\|prev \|<---\|prev \|<---\|prev \|<---\|prev \|<- - - -\|prev \|
732	* \| 0 \| \| 1 \| \| 2 \| \| 3 \| \| n \|
733	* \| next\|--->\| next\|--->\| next\|--->\| next\|- - - ->\| next\|->NULL
734	* +-----+ +-----+ +-----+ +-----+ +-----+
735	*
736	* 3) when finding a best pair of neighbor nodes, we use
737	* this chain list as loop counter.
738	*
739	* 4) If an entry was killed by node-joining, this chain list is
740	* modified to remove that entry.
741	*
742	* EX) remove the entry No 2.
743	* +-----+ +-----+ +-----+ +-----+
744	* NULL<-\|prev \|<---\|prev \|<--------------\|prev \|<- - - -\|prev \|
745	* \| 0 \| \| 1 \| \| 3 \| \| n \|
746	* \| next\|--->\| next\|-------------->\| next\|- - - ->\| next\|->NULL
747	* +-----+ +-----+ +-----+ +-----+
748	* +-----+
749	* NULL<-\|prev \|
750	* \| 2 \|
751	* \| next\|->NULL
752	* +-----+
753	*
754	* By this method, speed is up at the latter of finding a best pair of
755	* neighbor nodes.
756	*
757	*******************
758	* <IMPROVEMENT 3> : Cut the frequency of division.
759	*
760	* At comparison between 'total' and 'tmin' in the main cycle, total is
761	* divided by (2.0*fnseqs2) before comparison. If N nodes are available,
762	* that division happen (N*(N-1))/2 order.
763	*
764	* We thought that the comparison relation between tmin and total/(2.0*fnseqs2)
765	* is equal to the comparison relation between (tmin2.0fnseqs2) and total.
766	* Calculation of (tmin2.0fnseqs2) is only one time. so we stop dividing
767	* a total value and multiply tmin and (tmin2.0fnseqs2) instead.
768	*
769	*******************
770	* <IMPROVEMENT 4> : some transformation of the equation (to cut operations).
771	*
772	* We transform an equation of calculating 'total' in the main cycle.
773	*
774	*/
775
776
777	void fast_nj_tree(char *tree_description, FILE tree)
778	{
779	register int i;
780	sint l[4],nude,k;
781	sint nc,mini,minj,j,ii,jj;
782	double fnseqs,fnseqs2=0,sumd;
783	double diq,djq,dij,dio,djo,da;
784	double tmin,total,dmin;
785	double bi,bj,b1,b2,b3,branch[4];
786	sint typei,typej; /* 0 = node; 1 = OTU */
787
788	/* IMPROVEMENT 1, STEP 0 : declare variables */
789	double sum_cols, sum_rows, *join;
790
791	/* IMPROVEMENT 2, STEP 0 : declare variables */
792	sint loop_limit;
793	typedef struct _ValidNodeID {
794	sint n;
795	struct _ValidNodeID *prev;
796	struct _ValidNodeID *next;
797	} ValidNodeID;
798	ValidNodeID tvalid, lpi, lpj, lpii, lpjj, lp_prev, *lp_next;
799
800	/*
801	* correspondence of the loop counter variables.
802	* i .. lpi->n, ii .. lpii->n
803	* j .. lpj->n, jj .. lpjj->n
804	*/
805
806	fnseqs = (double)last_seq-first_seq+1;
807
808	/********************* First initialisation *************************/
809
810	if(verbose) {
811	fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
812	fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
813	fprintf(tree," The Neighbor-joining Method:");
814	fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
815	fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
816	fprintf(tree,"\n\n This is an UNROOTED tree\n");
817	fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
818	}
819
820	if (fnseqs == 2) {
821	if (verbose) fprintf(tree,"Cycle 1 = SEQ: 1 (%9.5f) joins SEQ: 2 (%9.5f)",tmat[first_seq][first_seq+1],tmat[first_seq][first_seq+1]);
822	return;
823	}
824
825	mini = minj = 0;
826
827	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
828	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
829	tkill = (sint ) ckalloc( (nseqs+1) sizeof (sint) );
830	av = (double ) ckalloc( (nseqs+1) sizeof (double) );
831
832	/* IMPROVEMENT 1, STEP 1 : Allocate memory */
833	sum_cols = (double ) ckalloc( (nseqs+1) sizeof (double) );
834	sum_rows = (double ) ckalloc( (nseqs+1) sizeof (double) );
835	join = (double ) ckalloc( (nseqs+1) sizeof (double) );
836
837	/* IMPROVEMENT 2, STEP 1 : Allocate memory */
838	tvalid = (ValidNodeID ) ckalloc( (nseqs+1) sizeof (ValidNodeID) );
839	/* tvalid[0] is special entry in array. it points a header of valid entry list */
840	tvalid[0].n = 0;
841	tvalid[0].prev = NULL;
842	tvalid[0].next = &tvalid[1];
843
844	/* IMPROVEMENT 2, STEP 2 : Construct and initialize the entry chain list */
845	for(i=1, loop_limit = last_seq-first_seq+1,
846	lpi=&tvalid[1], lp_prev=&tvalid[0], lp_next=&tvalid[2] ;
847	i<=loop_limit ;
848	++i, ++lpi, ++lp_prev, ++lp_next)
849	{
850	tmat[i][i] = av[i] = 0.0;
851	tkill[i] = 0;
852	lpi->n = i;
853	lpi->prev = lp_prev;
854	lpi->next = lp_next;
855
856	/* IMPROVEMENT 1, STEP 2 : Initialize arrays */
857	sum_cols[i] = sum_rows[i] = join[i] = 0.0;
858	}
859	tvalid[loop_limit].next = NULL;
860
861	/*
862	* IMPROVEMENT 1, STEP 3 : Calculate the sum of score value that
863	* is sequence[i] to others.
864	*/
865	sumd = 0.0;
866	for (lpj=tvalid[0].next ; lpj!=NULL ; lpj = lpj->next) {
867	double tmp_sum = 0.0;
868	j = lpj->n;
869	/* calculate sum_rows[j] */
870	for (lpi=tvalid[0].next ; lpi->n < j ; lpi = lpi->next) {
871	i = lpi->n;
872	tmp_sum += tmat[i][j];
873	/* tmat[j][i] = tmat[i][j]; */
874	}
875	sum_rows[j] = tmp_sum;
876
877	tmp_sum = 0.0;
878	/* Set lpi to that lpi->n is greater than j */
879	if ((lpi != NULL) && (lpi->n == j)) {
880	lpi = lpi->next;
881	}
882	/* calculate sum_cols[j] */
883	for( ; lpi!=NULL ; lpi = lpi->next) {
884	i = lpi->n;
885	tmp_sum += tmat[j][i];
886	/* tmat[i][j] = tmat[j][i]; */
887	}
888	sum_cols[j] = tmp_sum;
889	}
890
891	/********************* Enter The Main Cycle *************************/
892
893	for(nc=1, loop_limit = (last_seq-first_seq+1-3); nc<=loop_limit; ++nc) {
894
895	sumd = 0.0;
896	/* IMPROVEMENT 1, STEP 4 : use sum value */
897	for(lpj=tvalid[0].next ; lpj!=NULL ; lpj = lpj->next) {
898	sumd += sum_cols[lpj->n];
899	}
900
901	/* IMPROVEMENT 3, STEP 0 : multiply tmin and 2fnseqs2 /
902	fnseqs2 = fnseqs - 2.0; /* Set fnseqs2 at this point. */
903	tmin = 99999.0 * 2.0 * fnseqs2;
904
905
906	/.................compute SMATij values and find the smallest one ......../
907
908	mini = minj = 0;
909
910	/* jj must starts at least 2 */
911	if ((tvalid[0].next != NULL) && (tvalid[0].next->n == 1)) {
912	lpjj = tvalid[0].next->next;
913	} else {
914	lpjj = tvalid[0].next;
915	}
916
917	for( ; lpjj != NULL; lpjj = lpjj->next) {
918	jj = lpjj->n;
919	for(lpii=tvalid[0].next ; lpii->n < jj ; lpii = lpii->next) {
920	ii = lpii->n;
921	diq = djq = 0.0;
922
923	/* IMPROVEMENT 1, STEP 4 : use sum value */
924	diq = sum_cols[ii] + sum_rows[ii];
925	djq = sum_cols[jj] + sum_rows[jj];
926	/*
927	* always ii < jj in this point. Use upper
928	* triangle of score matrix.
929	*/
930	dij = tmat[ii][jj];
931
932	/*
933	* IMPROVEMENT 3, STEP 1 : fnseqs2 is
934	* already calculated.
935	*/
936	/* fnseqs2 = fnseqs - 2.0 */
937
938	/* IMPROVEMENT 4 : transform the equation */
939	/-------------------------------------------------------------------
940	* OPTIMIZE of expression 'total = d2r + fnseqs2dij + dr2.0' *
941	* total = d2r + fnseq2dij + 2.0dr *
942	* = d2r + fnseq2dij + 2(sumd - dij - d2r)
943	* = d2r + fnseq2dij + 2sumd - 2dij - 2d2r *
944	* = fnseq2dij + 2sumd - 2dij - 2d2r + d2r *
945	* = fnseq2dij + 2sumd - 2dij - d2r
946	* = fnseq2dij + 2sumd - 2dij - (diq + djq - 2dij) *
947	* = fnseq2dij + 2sumd - 2dij - diq - djq + 2dij *
948	* = fnseq2dij + 2sumd - 2dij + 2dij - diq - djq *
949	* = fnseq2dij + 2sumd - diq - djq *
950	-------------------------------------------------------------------/
951	total = fnseqs2dij + 2.0sumd - diq - djq;
952
953	/*
954	* IMPROVEMENT 3, STEP 2 : abbrevlate
955	* the division on comparison between
956	* total and tmin.
957	*/
958	/* total = total / (2.0fnseqs2); /
959
960	if(total < tmin) {
961	tmin = total;
962	mini = ii;
963	minj = jj;
964	}
965	}
966	}
967
968	/* MEMO: always ii < jj in avobe loop, so mini < minj */
969
970	/.................compute branch lengths and print the results ......../
971
972
973	dio = djo = 0.0;
974
975	/* IMPROVEMENT 1, STEP 4 : use sum value */
976	dio = sum_cols[mini] + sum_rows[mini];
977	djo = sum_cols[minj] + sum_rows[minj];
978
979	dmin = tmat[mini][minj];
980	dio = (dio - dmin) / fnseqs2;
981	djo = (djo - dmin) / fnseqs2;
982	bi = (dmin + dio - djo) * 0.5;
983	bj = dmin - bi;
984	bi = bi - av[mini];
985	bj = bj - av[minj];
986
987	if( av[mini] > 0.0 )
988	typei = 0;
989	else
990	typei = 1;
991	if( av[minj] > 0.0 )
992	typej = 0;
993	else
994	typej = 1;
995
996	if(verbose)
997	fprintf(tree,"\n Cycle%4d = ",(pint)nc);
998
999	/*
1000	set negative branch lengths to zero. Also set any tiny positive
1001	branch lengths to zero.
1002	*/ if( fabs(bi) < 0.0001) bi = 0.0;
1003	if( fabs(bj) < 0.0001) bj = 0.0;
1004
1005	if(verbose) {
1006	if(typei == 0)
1007	fprintf(tree,"Node:%4d (%9.5f) joins ",(pint)mini,bi);
1008	else
1009	fprintf(tree," SEQ:%4d (%9.5f) joins ",(pint)mini,bi);
1010
1011	if(typej == 0)
1012	fprintf(tree,"Node:%4d (%9.5f)",(pint)minj,bj);
1013	else
1014	fprintf(tree," SEQ:%4d (%9.5f)",(pint)minj,bj);
1015
1016	fprintf(tree,"\n");
1017	}
1018
1019
1020	left_branch[nc] = bi;
1021	right_branch[nc] = bj;
1022
1023	for(i=1; i<=last_seq-first_seq+1; i++)
1024	tree_description[nc][i] = 0;
1025
1026	if(typei == 0) {
1027	for(i=nc-1; i>=1; i--)
1028	if(tree_description[i][mini] == 1) {
1029	for(j=1; j<=last_seq-first_seq+1; j++)
1030	if(tree_description[i][j] == 1)
1031	tree_description[nc][j] = 1;
1032	break;
1033	}
1034	}
1035	else
1036	tree_description[nc][mini] = 1;
1037
1038	if(typej == 0) {
1039	for(i=nc-1; i>=1; i--)
1040	if(tree_description[i][minj] == 1) {
1041	for(j=1; j<=last_seq-first_seq+1; j++)
1042	if(tree_description[i][j] == 1)
1043	tree_description[nc][j] = 1;
1044	break;
1045	}
1046	}
1047	else
1048	tree_description[nc][minj] = 1;
1049
1050
1051	/*
1052	Here is where the -0.00005 branch lengths come from for 3 or more
1053	identical seqs.
1054	*/
1055	/* if(dmin <= 0.0) dmin = 0.0001; */
1056	if(dmin <= 0.0) dmin = 0.000001;
1057	av[mini] = dmin * 0.5;
1058
1059	/........................Re-initialisation................................/
1060
1061	fnseqs = fnseqs - 1.0;
1062	tkill[minj] = 1;
1063
1064	/* IMPROVEMENT 2, STEP 3 : Remove tvalid[minj] from chain list. */
1065	/* [ Before ]
1066	* +---------+ +---------+ +---------+
1067	* \|prev \|<-------\|prev \|<-------\|prev \|<---
1068	* \| n \| \| n(=minj)\| \| n \|
1069	* \| next\|------->\| next\|------->\| next\|----
1070	* +---------+ +---------+ +---------+
1071	*
1072	* [ After ]
1073	* +---------+ +---------+
1074	* \|prev \|<--------------------------\|prev \|<---
1075	* \| n \| \| n \|
1076	* \| next\|-------------------------->\| next\|----
1077	* +---------+ +---------+
1078	* +---------+
1079	* NULL---\|prev \|
1080	* \| n(=minj)\|
1081	* \| next\|---NULL
1082	* +---------+
1083	*/
1084	(tvalid[minj].prev)->next = tvalid[minj].next;
1085	if (tvalid[minj].next != NULL) {
1086	(tvalid[minj].next)->prev = tvalid[minj].prev;
1087	}
1088	tvalid[minj].prev = tvalid[minj].next = NULL;
1089
1090	/* IMPROVEMENT 1, STEP 5 : re-calculate sum values. */
1091	for(lpj=tvalid[0].next ; lpj != NULL ; lpj = lpj->next) {
1092	double tmp_di = 0.0;
1093	double tmp_dj = 0.0;
1094	j = lpj->n;
1095
1096	/*
1097	* subtrace a score value related with 'minj' from
1098	* sum arrays .
1099	*/
1100	if (j < minj) {
1101	tmp_dj = tmat[j][minj];
1102	sum_cols[j] -= tmp_dj;
1103	} else if (j > minj) {
1104	tmp_dj = tmat[minj][j];
1105	sum_rows[j] -= tmp_dj;
1106	} /* nothing to do when j is equal to minj. */
1107
1108
1109	/*
1110	* subtrace a score value related with 'mini' from
1111	* sum arrays .
1112	*/
1113	if (j < mini) {
1114	tmp_di = tmat[j][mini];
1115	sum_cols[j] -= tmp_di;
1116	} else if (j > mini) {
1117	tmp_di = tmat[mini][j];
1118	sum_rows[j] -= tmp_di;
1119	} /* nothing to do when j is equal to mini. */
1120
1121	/*
1122	* calculate a score value of the new inner node.
1123	* then, store it temporary to join[] array.
1124	*/
1125	join[j] = (tmp_dj + tmp_di) * 0.5;
1126	}
1127
1128	/*
1129	* 1)
1130	* Set the score values (stored in join[]) into the matrix,
1131	* row/column position is 'mini'.
1132	* 2)
1133	* Add a score value of the new inner node to sum arrays.
1134	*/
1135	for(lpj=tvalid[0].next ; lpj != NULL; lpj = lpj->next) {
1136	j = lpj->n;
1137	if (j < mini) {
1138	tmat[j][mini] = join[j];
1139	sum_cols[j] += join[j];
1140	} else if (j > mini) {
1141	tmat[mini][j] = join[j];
1142	sum_rows[j] += join[j];
1143	} /* nothing to do when j is equal to mini. */
1144	}
1145
1146	/* Re-calculate sum_rows[mini],sum_cols[mini]. */
1147	sum_cols[mini] = sum_rows[mini] = 0.0;
1148
1149	/* calculate sum_rows[mini] */
1150	da = 0.0;
1151	for(lpj=tvalid[0].next ; lpj->n < mini ; lpj = lpj->next) {
1152	da += join[lpj->n];
1153	}
1154	sum_rows[mini] = da;
1155
1156	/* skip if 'lpj->n' is equal to 'mini' */
1157	if ((lpj != NULL) && (lpj->n == mini)) {
1158	lpj = lpj->next;
1159	}
1160
1161	/* calculate sum_cols[mini] */
1162	da = 0.0;
1163	for( ; lpj != NULL; lpj = lpj->next) {
1164	da += join[lpj->n];
1165	}
1166	sum_cols[mini] = da;
1167
1168	/*
1169	* Clean up sum_rows[minj], sum_cols[minj] and score matrix
1170	* related with 'minj'.
1171	*/
1172	sum_cols[minj] = sum_rows[minj] = 0.0;
1173	for(j=1; j<=last_seq-first_seq+1; ++j)
1174	tmat[minj][j] = tmat[j][minj] = join[j] = 0.0;
1175
1176
1177	/**/ } /end main cycle**/
1178
1179	/****************************Last Cycle (3 Seqs. left)******************/
1180
1181	nude = 1;
1182
1183	for(lpi=tvalid[0].next; lpi != NULL; lpi = lpi->next) {
1184	l[nude] = lpi->n;
1185	++nude;
1186	}
1187
1188	b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
1189	b2 = tmat[l[1]][l[2]] - b1;
1190	b3 = tmat[l[1]][l[3]] - b1;
1191
1192	branch[1] = b1 - av[l[1]];
1193	branch[2] = b2 - av[l[2]];
1194	branch[3] = b3 - av[l[3]];
1195
1196	/* Reset tiny negative and positive branch lengths to zero */
1197	if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
1198	if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
1199	if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
1200
1201	left_branch[last_seq-first_seq+1-2] = branch[1];
1202	left_branch[last_seq-first_seq+1-1] = branch[2];
1203	left_branch[last_seq-first_seq+1] = branch[3];
1204
1205	for(i=1; i<=last_seq-first_seq+1; i++)
1206	tree_description[last_seq-first_seq+1-2][i] = 0;
1207
1208	if(verbose)
1209	fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",(pint)nc);
1210
1211	for(i=1; i<=3; ++i) {
1212	if( av[l[i]] > 0.0) {
1213	if(verbose)
1214	fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",(pint)l[i],branch[i]);
1215	for(k=last_seq-first_seq+1-3; k>=1; k--)
1216	if(tree_description[k][l[i]] == 1) {
1217	for(j=1; j<=last_seq-first_seq+1; j++)
1218	if(tree_description[k][j] == 1)
1219	tree_description[last_seq-first_seq+1-2][j] = i;
1220	break;
1221	}
1222	}
1223	else {
1224	if(verbose)
1225	fprintf(tree,"\n\t\t SEQ:%4d (%9.5f) ",(pint)l[i],branch[i]);
1226	tree_description[last_seq-first_seq+1-2][l[i]] = i;
1227	}
1228	if(i < 3) {
1229	if(verbose)
1230	fprintf(tree,"joins");
1231	}
1232	}
1233
1234	if(verbose)
1235	fprintf(tree,"\n");
1236
1237
1238	/* IMPROVEMENT 1, STEP 6 : release memory area */
1239	ckfree(sum_cols);
1240	ckfree(sum_rows);
1241	ckfree(join);
1242
1243	/* IMPROVEMENT 2, STEP 4 : release memory area */
1244	ckfree(tvalid);
1245
1246	}
1247	#endif /* ORIGINAL_NJ_TREE */
1248
1249
1250
1251	void bootstrap_tree(char phylip_name,char clustal_name, char *nexus_name)
1252	{
1253	sint i,j;
1254	int ranno;
1255	char path[MAXLINE+1];
1256	char dummy[10];
1257	static char **sample_tree;
1258	static char **standard_tree;
1259	static char **save_tree;
1260	sint total_dists, overspill = 0, total_overspill = 0;
1261	sint nfails = 0;
1262
1263	if(empty) {
1264	error("You must load an alignment first");
1265	return;
1266	}
1267
1268	if(nseqs<4) {
1269	error("Alignment has only %d sequences",nseqs);
1270	return;
1271	}
1272
1273	if(!output_tree_clustal && !output_tree_phylip && !output_tree_nexus) {
1274	error("You must select either clustal or phylip or nexus tree output format");
1275	return;
1276	}
1277	get_path(seqname, path);
1278
1279	if (output_tree_clustal) {
1280	if (clustal_name[0]!=EOS) {
1281	if((clustal_phy_tree_file = open_explicit_file(
1282	clustal_name))==NULL) return;
1283	}
1284	else {
1285	if((clustal_phy_tree_file = open_output_file(
1286	"\nEnter name for bootstrap output file ",path,
1287	clustal_name,"njb")) == NULL) return;
1288	}
1289	}
1290
1291	first_seq=1;
1292	last_seq=nseqs;
1293
1294	if (output_tree_phylip) {
1295	if (phylip_name[0]!=EOS) {
1296	if((phylip_phy_tree_file = open_explicit_file(
1297	phylip_name))==NULL) return;
1298	}
1299	else {
1300	if((phylip_phy_tree_file = open_output_file(
1301	"\nEnter name for bootstrap output file ",path,
1302	phylip_name,"phb")) == NULL) return;
1303	}
1304	}
1305
1306	if (output_tree_nexus) {
1307	if (nexus_name[0]!=EOS) {
1308	if((nexus_phy_tree_file = open_explicit_file(
1309	nexus_name))==NULL) return;
1310	}
1311	else {
1312	if((nexus_phy_tree_file = open_output_file(
1313	"\nEnter name for bootstrap output file ",path,
1314	nexus_name,"treb")) == NULL) return;
1315	}
1316	}
1317
1318	boot_totals = (sint )ckalloc( (nseqs+1) sizeof (sint) );
1319	for(i=0;i<nseqs+1;i++)
1320	boot_totals[i]=0;
1321
1322	boot_positions = (sint )ckalloc( (seqlen_array[first_seq]+2) sizeof (sint) );
1323
1324	for(j=1; j<=seqlen_array[first_seq]; ++j) /* First select all positions for */
1325	boot_positions[j] = j; /* the "standard" tree */
1326
1327	if(output_tree_clustal) {
1328	verbose = TRUE; /* Turn on file output */
1329	if(dnaflag)
1330	overspill = dna_distance_matrix(clustal_phy_tree_file);
1331	else
1332	overspill = prot_distance_matrix(clustal_phy_tree_file);
1333	}
1334
1335	if(output_tree_phylip) {
1336	verbose = FALSE; /* Turn off file output */
1337	if(dnaflag)
1338	overspill = dna_distance_matrix(phylip_phy_tree_file);
1339	else
1340	overspill = prot_distance_matrix(phylip_phy_tree_file);
1341	}
1342
1343	if(output_tree_nexus) {
1344	verbose = FALSE; /* Turn off file output */
1345	if(dnaflag)
1346	overspill = dna_distance_matrix(nexus_phy_tree_file);
1347	else
1348	overspill = prot_distance_matrix(nexus_phy_tree_file);
1349	}
1350
1351	/* check if any distances overflowed the distance corrections */
1352	if ( overspill > 0 ) {
1353	total_dists = (nseqs*(nseqs-1))/2;
1354	overspill_message(overspill,total_dists);
1355	}
1356
1357	tree_gaps=ckfree((void *)tree_gaps);
1358
1359	if (output_tree_clustal) verbose = TRUE; /* Turn on screen output */
1360
1361	standard_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
1362	for(i=0; i<nseqs+1; i++)
1363	standard_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
1364
1365	/* compute the standard tree */
1366
1367	if(output_tree_clustal \|\| output_tree_phylip \|\| output_tree_nexus)
1368	nj_tree(standard_tree,clustal_phy_tree_file);
1369
1370	if (output_tree_clustal)
1371	fprintf(clustal_phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n");
1372
1373	/* save the left_branch and right_branch for phylip output */
1374	save_left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1375	save_right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1376	for (i=1;i<=nseqs;i++) {
1377	save_left_branch[i] = left_branch[i];
1378	save_right_branch[i] = right_branch[i];
1379	}
1380	/*
1381	The next line is a fossil from the days of using the cc ran()
1382	ran_factor = RAND_MAX / seqlen_array[first_seq];
1383	*/
1384
1385	if(usemenu)
1386	boot_ran_seed =
1387	getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed);
1388
1389	/* do not use the native cc ran()
1390	srand(boot_ran_seed);
1391	*/
1392	addrandinit((unsigned long) boot_ran_seed);
1393
1394	if (output_tree_clustal)
1395	fprintf(clustal_phy_tree_file,"\n Random number generator seed = %7u\n",
1396	boot_ran_seed);
1397
1398	if(usemenu)
1399	boot_ntrials =
1400	getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials);
1401
1402	if (output_tree_clustal) {
1403	fprintf(clustal_phy_tree_file,"\n Number of bootstrap trials = %7d\n",
1404	(pint)boot_ntrials);
1405
1406	fprintf(clustal_phy_tree_file,
1407	"\n\n Diagrammatic representation of the above tree: \n");
1408	fprintf(clustal_phy_tree_file,"\n Each row represents 1 tree cycle;");
1409	fprintf(clustal_phy_tree_file," defining 2 groups.\n");
1410	fprintf(clustal_phy_tree_file,"\n Each column is 1 sequence; ");
1411	fprintf(clustal_phy_tree_file,"the stars in each line show 1 group; ");
1412	fprintf(clustal_phy_tree_file,"\n the dots show the other\n");
1413	fprintf(clustal_phy_tree_file,"\n Numbers show occurrences in bootstrap samples.");
1414	}
1415	/*
1416	print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
1417	*/
1418	verbose = FALSE; /* Turn OFF screen output */
1419
1420	left_branch=ckfree((void *)left_branch);
1421	right_branch=ckfree((void *)right_branch);
1422	tkill=ckfree((void *)tkill);
1423	av=ckfree((void *)av);
1424
1425	sample_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
1426	for(i=0; i<nseqs+1; i++)
1427	sample_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
1428
1429	if (usemenu)
1430	fprintf(stdout,"\n\nEach dot represents 10 trials\n\n");
1431	total_overspill = 0;
1432	nfails = 0;
1433	for(i=1; i<=boot_ntrials; ++i) {
1434	for(j=1; j<=seqlen_array[first_seq]; ++j) { /* select alignment */
1435	/* positions for */
1436	ranno = addrand( (unsigned long) seqlen_array[1]) + 1;
1437	boot_positions[j] = ranno; /* bootstrap sample */
1438	}
1439	if(output_tree_clustal) {
1440	if(dnaflag)
1441	overspill = dna_distance_matrix(clustal_phy_tree_file);
1442	else
1443	overspill = prot_distance_matrix(clustal_phy_tree_file);
1444	}
1445
1446	if(output_tree_phylip) {
1447	if(dnaflag)
1448	overspill = dna_distance_matrix(phylip_phy_tree_file);
1449	else
1450	overspill = prot_distance_matrix(phylip_phy_tree_file);
1451	}
1452
1453	if(output_tree_nexus) {
1454	if(dnaflag)
1455	overspill = dna_distance_matrix(nexus_phy_tree_file);
1456	else
1457	overspill = prot_distance_matrix(nexus_phy_tree_file);
1458	}
1459
1460	if( overspill > 0) {
1461	total_overspill = total_overspill + overspill;
1462	nfails++;
1463	}
1464
1465	tree_gaps=ckfree((void *)tree_gaps);
1466
1467	if(output_tree_clustal \|\| output_tree_phylip \|\| output_tree_nexus)
1468	nj_tree(sample_tree,clustal_phy_tree_file);
1469
1470	left_branch=ckfree((void *)left_branch);
1471	right_branch=ckfree((void *)right_branch);
1472	tkill=ckfree((void *)tkill);
1473	av=ckfree((void *)av);
1474
1475	compare_tree(standard_tree, sample_tree, boot_totals, last_seq-first_seq+1);
1476	if (usemenu) {
1477	if(i % 10 == 0) fprintf(stdout,".");
1478	if(i % 100 == 0) fprintf(stdout,"\n");
1479	}
1480	}
1481
1482	/* check if any distances overflowed the distance corrections */
1483	if ( nfails > 0 ) {
1484	total_dists = (nseqs*(nseqs-1))/2;
1485	fprintf(stdout,"\n");
1486	fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld times %ld",
1487	(long)total_overspill,(long)total_dists,(long)boot_ntrials);
1488	fprintf(stdout,"\n were out of range for the distance correction.");
1489	fprintf(stdout,"\n This affected %d out of %d bootstrap trials.",
1490	(pint)nfails,(pint)boot_ntrials);
1491	fprintf(stdout,"\n This may not be fatal but you have been warned!");
1492	fprintf(stdout,"\n");
1493	fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
1494	fprintf(stdout,"\n or 2) remove the most distant sequences");
1495	fprintf(stdout,"\n or 3) use the PHYLIP package.");
1496	fprintf(stdout,"\n\n");
1497	if (usemenu)
1498	getstr("Press [RETURN] to continue",dummy, 10);
1499	}
1500
1501
1502	boot_positions=ckfree((void *)boot_positions);
1503
1504	for (i=1;i<nseqs+1;i++)
1505	sample_tree[i]=ckfree((void *)sample_tree[i]);
1506	sample_tree=ckfree((void *)sample_tree);
1507	/*
1508	fprintf(clustal_phy_tree_file,"\n\n Bootstrap totals for each group\n");
1509	*/
1510	if (output_tree_clustal)
1511	print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
1512
1513	save_tree = (char *) ckalloc( (nseqs+1) sizeof (char *) );
1514	for(i=0; i<nseqs+1; i++)
1515	save_tree[i] = (char ) ckalloc( (nseqs+1) sizeof(char) );
1516
1517	for(i=1; i<nseqs+1; i++)
1518	for(j=1; j<nseqs+1; j++)
1519	save_tree[i][j] = standard_tree[i][j];
1520
1521	if(output_tree_phylip) {
1522	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1523	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1524	for (i=1;i<=nseqs;i++) {
1525	left_branch[i] = save_left_branch[i];
1526	right_branch[i] = save_right_branch[i];
1527	}
1528	print_phylip_tree(standard_tree,phylip_phy_tree_file,
1529	bootstrap_format);
1530	left_branch=ckfree((void *)left_branch);
1531	right_branch=ckfree((void *)right_branch);
1532	}
1533
1534	for(i=1; i<nseqs+1; i++)
1535	for(j=1; j<nseqs+1; j++)
1536	standard_tree[i][j] = save_tree[i][j];
1537
1538	if(output_tree_nexus) {
1539	left_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1540	right_branch = (double ) ckalloc( (nseqs+2) sizeof (double) );
1541	for (i=1;i<=nseqs;i++) {
1542	left_branch[i] = save_left_branch[i];
1543	right_branch[i] = save_right_branch[i];
1544	}
1545	print_nexus_tree(standard_tree,nexus_phy_tree_file,
1546	bootstrap_format);
1547	left_branch=ckfree((void *)left_branch);
1548	right_branch=ckfree((void *)right_branch);
1549	}
1550
1551	boot_totals=ckfree((void *)boot_totals);
1552	save_left_branch=ckfree((void *)save_left_branch);
1553	save_right_branch=ckfree((void *)save_right_branch);
1554
1555	for (i=1;i<nseqs+1;i++)
1556	standard_tree[i]=ckfree((void *)standard_tree[i]);
1557	standard_tree=ckfree((void *)standard_tree);
1558
1559	for (i=0;i<nseqs+1;i++)
1560	save_tree[i]=ckfree((void *)save_tree[i]);
1561	save_tree=ckfree((void *)save_tree);
1562
1563	if (output_tree_clustal)
1564	fclose(clustal_phy_tree_file);
1565
1566	if (output_tree_phylip)
1567	fclose(phylip_phy_tree_file);
1568
1569	if (output_tree_nexus)
1570	fclose(nexus_phy_tree_file);
1571
1572	if (output_tree_clustal)
1573	info("Bootstrap output file completed [%s]"
1574	,clustal_name);
1575	if (output_tree_phylip)
1576	info("Bootstrap output file completed [%s]"
1577	,phylip_name);
1578	if (output_tree_nexus)
1579	info("Bootstrap output file completed [%s]"
1580	,nexus_name);
1581	}
1582
1583
1584	void compare_tree(char tree1, char tree2, sint *hits, sint n)
1585	{
1586	sint i,j,k;
1587	sint nhits1, nhits2;
1588
1589	for(i=1; i<=n-3; i++) {
1590	for(j=1; j<=n-3; j++) {
1591	nhits1 = 0;
1592	nhits2 = 0;
1593	for(k=1; k<=n; k++) {
1594	if(tree1[i][k] == tree2[j][k]) nhits1++;
1595	if(tree1[i][k] != tree2[j][k]) nhits2++;
1596	}
1597	if((nhits1 == last_seq-first_seq+1) \|\| (nhits2 == last_seq-first_seq+1)) hits[i]++;
1598	}
1599	}
1600	}
1601
1602
1603	void print_nexus_tree(char *tree_description, FILE tree, sint bootstrap)
1604	{
1605	sint i;
1606	sint old_row;
1607
1608	fprintf(tree,"#NEXUS\n\n");
1609
1610	fprintf(tree,"BEGIN TREES;\n\n");
1611	fprintf(tree,"\tTRANSLATE\n");
1612	for(i=1;i<nseqs;i++) {
1613	fprintf(tree,"\t\t%d %s,\n",(pint)i,names[i]);
1614	}
1615	fprintf(tree,"\t\t%d %s\n",(pint)nseqs,names[nseqs]);
1616	fprintf(tree,"\t\t;\n");
1617
1618	fprintf(tree,"\tUTREE PAUP_1= ");
1619
1620	if(last_seq-first_seq+1==2) {
1621	fprintf(tree,"(%d:%7.5f,%d:%7.5f);",first_seq,tmat[first_seq][first_seq+1],first_seq+1,tmat[first_seq][first_seq+1]);
1622	}
1623	else {
1624
1625	fprintf(tree,"(");
1626
1627	old_row=two_way_split_nexus(tree_description, tree, last_seq-first_seq+1-2,1,bootstrap);
1628	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-2]);
1629	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1630	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1631	fprintf(tree,",");
1632
1633	old_row=two_way_split_nexus(tree_description, tree, last_seq-first_seq+1-2,2,bootstrap);
1634	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-1]);
1635	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1636	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1637	fprintf(tree,",");
1638
1639	old_row=two_way_split_nexus(tree_description, tree, last_seq-first_seq+1-2,3,bootstrap);
1640	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1]);
1641	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1642	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1643	fprintf(tree,")");
1644	if (bootstrap==BS_NODE_LABELS) fprintf(tree,"TRICHOTOMY");
1645	fprintf(tree,";");
1646	}
1647	fprintf(tree,"\nENDBLOCK;\n");
1648	}
1649
1650
1651	sint two_way_split_nexus
1652	(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap)
1653	{
1654	sint row, new_row = 0, old_row, col, test_col = 0;
1655	Boolean single_seq;
1656
1657	if(start_row != last_seq-first_seq+1-2) fprintf(tree,"(");
1658
1659	for(col=1; col<=last_seq-first_seq+1; col++) {
1660	if(tree_description[start_row][col] == flag) {
1661	test_col = col;
1662	break;
1663	}
1664	}
1665
1666	single_seq = TRUE;
1667	for(row=start_row-1; row>=1; row--)
1668	if(tree_description[row][test_col] == 1) {
1669	single_seq = FALSE;
1670	new_row = row;
1671	break;
1672	}
1673
1674	if(single_seq) {
1675	tree_description[start_row][test_col] = 0;
1676	fprintf(tree,"%d",test_col+first_seq-1);
1677	if(start_row == last_seq-first_seq+1-2) {
1678	return(0);
1679	}
1680
1681	fprintf(tree,":%7.5f,",left_branch[start_row]);
1682	}
1683	else {
1684	for(col=1; col<=last_seq-first_seq+1; col++) {
1685	if((tree_description[start_row][col]==1)&&
1686	(tree_description[new_row][col]==1))
1687	tree_description[start_row][col] = 0;
1688	}
1689	old_row=two_way_split_nexus(tree_description, tree, new_row, (sint)1, bootstrap);
1690	if(start_row == last_seq-first_seq+1-2) {
1691	return(new_row);
1692	}
1693
1694	fprintf(tree,":%7.5f",left_branch[start_row]);
1695	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
1696	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1697
1698	fprintf(tree,",");
1699	}
1700
1701
1702	for(col=1; col<=last_seq-first_seq+1; col++)
1703	if(tree_description[start_row][col] == flag) {
1704	test_col = col;
1705	break;
1706	}
1707
1708	single_seq = TRUE;
1709	new_row = 0;
1710	for(row=start_row-1; row>=1; row--)
1711	if(tree_description[row][test_col] == 1) {
1712	single_seq = FALSE;
1713	new_row = row;
1714	break;
1715	}
1716
1717	if(single_seq) {
1718	tree_description[start_row][test_col] = 0;
1719	fprintf(tree,"%d",test_col+first_seq-1);
1720	fprintf(tree,":%7.5f)",right_branch[start_row]);
1721	}
1722	else {
1723	for(col=1; col<=last_seq-first_seq+1; col++) {
1724	if((tree_description[start_row][col]==1)&&
1725	(tree_description[new_row][col]==1))
1726	tree_description[start_row][col] = 0;
1727	}
1728	old_row=two_way_split_nexus(tree_description, tree, new_row, (sint)1, bootstrap);
1729	fprintf(tree,":%7.5f",right_branch[start_row]);
1730	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
1731	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1732
1733	fprintf(tree,")");
1734	}
1735	if ((bootstrap==BS_NODE_LABELS) && (boot_totals[start_row]>0))
1736	fprintf(tree,"%d",(pint)boot_totals[start_row]);
1737
1738	return(start_row);
1739	}
1740
1741
1742	void print_phylip_tree(char *tree_description, FILE tree, sint bootstrap)
1743	{
1744	sint old_row;
1745
1746	if(last_seq-first_seq+1==2) {
1747	fprintf(tree,"(%s:%7.5f,%s:%7.5f);",names[first_seq],tmat[first_seq][first_seq+1],names[first_seq+1],tmat[first_seq][first_seq+1]);
1748	return;
1749	}
1750
1751	fprintf(tree,"(\n");
1752
1753	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,1,bootstrap);
1754	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-2]);
1755	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1756	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1757	fprintf(tree,",\n");
1758
1759	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,2,bootstrap);
1760	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1-1]);
1761	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1762	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1763	fprintf(tree,",\n");
1764
1765	old_row=two_way_split(tree_description, tree, last_seq-first_seq+1-2,3,bootstrap);
1766	fprintf(tree,":%7.5f",left_branch[last_seq-first_seq+1]);
1767	if ((bootstrap==BS_BRANCH_LABELS) && (old_row>0) && (boot_totals[old_row]>0))
1768	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1769	fprintf(tree,")");
1770	if (bootstrap==BS_NODE_LABELS) fprintf(tree,"TRICHOTOMY");
1771	fprintf(tree,";\n");
1772	}
1773
1774
1775	sint two_way_split
1776	(char *tree_description, FILE tree, sint start_row, sint flag, sint bootstrap)
1777	{
1778	sint row, new_row = 0, old_row, col, test_col = 0;
1779	Boolean single_seq;
1780
1781	if(start_row != last_seq-first_seq+1-2) fprintf(tree,"(\n");
1782
1783	for(col=1; col<=last_seq-first_seq+1; col++) {
1784	if(tree_description[start_row][col] == flag) {
1785	test_col = col;
1786	break;
1787	}
1788	}
1789
1790	single_seq = TRUE;
1791	for(row=start_row-1; row>=1; row--)
1792	if(tree_description[row][test_col] == 1) {
1793	single_seq = FALSE;
1794	new_row = row;
1795	break;
1796	}
1797
1798	if(single_seq) {
1799	tree_description[start_row][test_col] = 0;
1800	fprintf(tree,"%.*s",max_names,names[test_col+first_seq-1]);
1801	if(start_row == last_seq-first_seq+1-2) {
1802	return(0);
1803	}
1804
1805	fprintf(tree,":%7.5f,\n",left_branch[start_row]);
1806	}
1807	else {
1808	for(col=1; col<=last_seq-first_seq+1; col++) {
1809	if((tree_description[start_row][col]==1)&&
1810	(tree_description[new_row][col]==1))
1811	tree_description[start_row][col] = 0;
1812	}
1813	old_row=two_way_split(tree_description, tree, new_row, (sint)1, bootstrap);
1814	if(start_row == last_seq-first_seq+1-2) {
1815	return(new_row);
1816	}
1817
1818	fprintf(tree,":%7.5f",left_branch[start_row]);
1819	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
1820	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1821
1822	fprintf(tree,",\n");
1823	}
1824
1825
1826	for(col=1; col<=last_seq-first_seq+1; col++)
1827	if(tree_description[start_row][col] == flag) {
1828	test_col = col;
1829	break;
1830	}
1831
1832	single_seq = TRUE;
1833	new_row = 0;
1834	for(row=start_row-1; row>=1; row--)
1835	if(tree_description[row][test_col] == 1) {
1836	single_seq = FALSE;
1837	new_row = row;
1838	break;
1839	}
1840
1841	if(single_seq) {
1842	tree_description[start_row][test_col] = 0;
1843	fprintf(tree,"%.*s",max_names,names[test_col+first_seq-1]);
1844	fprintf(tree,":%7.5f)\n",right_branch[start_row]);
1845	}
1846	else {
1847	for(col=1; col<=last_seq-first_seq+1; col++) {
1848	if((tree_description[start_row][col]==1)&&
1849	(tree_description[new_row][col]==1))
1850	tree_description[start_row][col] = 0;
1851	}
1852	old_row=two_way_split(tree_description, tree, new_row, (sint)1, bootstrap);
1853	fprintf(tree,":%7.5f",right_branch[start_row]);
1854	if ((bootstrap==BS_BRANCH_LABELS) && (boot_totals[old_row]>0))
1855	fprintf(tree,"[%d]",(pint)boot_totals[old_row]);
1856
1857	fprintf(tree,")\n");
1858	}
1859	if ((bootstrap==BS_NODE_LABELS) && (boot_totals[start_row]>0))
1860	fprintf(tree,"%d",(pint)boot_totals[start_row]);
1861
1862	return(start_row);
1863	}
1864
1865
1866
1867	void print_tree(char *tree_description, FILE tree, sint *totals)
1868	{
1869	sint row,col;
1870
1871	fprintf(tree,"\n");
1872
1873	for(row=1; row<=last_seq-first_seq+1-3; row++) {
1874	fprintf(tree," \n");
1875	for(col=1; col<=last_seq-first_seq+1; col++) {
1876	if(tree_description[row][col] == 0)
1877	fprintf(tree,"*");
1878	else
1879	fprintf(tree,".");
1880	}
1881	if(totals[row] > 0)
1882	fprintf(tree,"%7d",(pint)totals[row]);
1883	}
1884	fprintf(tree," \n");
1885	for(col=1; col<=last_seq-first_seq+1; col++)
1886	fprintf(tree,"%1d",(pint)tree_description[last_seq-first_seq+1-2][col]);
1887	fprintf(tree,"\n");
1888	}
1889
1890
1891
1892	sint dna_distance_matrix(FILE *tree)
1893	{
1894	sint m,n;
1895	sint j,i;
1896	sint res1, res2;
1897	sint overspill = 0;
1898	double p,q,e,a,b,k;
1899
1900	tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */
1901
1902	if(verbose) {
1903	fprintf(tree,"\n");
1904	fprintf(tree,"\n DIST = percentage divergence (/100)");
1905	fprintf(tree,"\n p = rate of transition (A <-> G; C <-> T)");
1906	fprintf(tree,"\n q = rate of transversion");
1907	fprintf(tree,"\n Length = number of sites used in comparison");
1908	fprintf(tree,"\n");
1909	if(tossgaps) {
1910	fprintf(tree,"\n All sites with gaps (in any sequence) deleted!");
1911	fprintf(tree,"\n");
1912	}
1913	if(kimura) {
1914	fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:");
1915	fprintf(tree,"\n\n Kimura, M. (1980)");
1916	fprintf(tree," A simple method for estimating evolutionary ");
1917	fprintf(tree,"rates of base");
1918	fprintf(tree,"\n substitutions through comparative studies of ");
1919	fprintf(tree,"nucleotide sequences.");
1920	fprintf(tree,"\n J. Mol. Evol., 16, 111-120.");
1921	fprintf(tree,"\n\n");
1922	}
1923	}
1924
1925	for(m=1; m<last_seq-first_seq+1; ++m) /* for every pair of sequence */
1926	for(n=m+1; n<=last_seq-first_seq+1; ++n) {
1927	p = q = e = 0.0;
1928	tmat[m][n] = tmat[n][m] = 0.0;
1929	for(i=1; i<=seqlen_array[first_seq]; ++i) {
1930	j = boot_positions[i];
1931	if(tossgaps && (tree_gaps[j] > 0) )
1932	goto skip; /* gap position */
1933	res1 = seq_array[m+first_seq-1][j];
1934	res2 = seq_array[n+first_seq-1][j];
1935	if( (res1 == gap_pos1) \|\| (res1 == gap_pos2) \|\|
1936	(res2 == gap_pos1) \|\| (res2 == gap_pos2))
1937	goto skip; /* gap in a seq*/
1938	if(!use_ambiguities)
1939	if( is_ambiguity(res1) \|\| is_ambiguity(res2))
1940	goto skip; /* ambiguity code in a seq*/
1941	e = e + 1.0;
1942	if(res1 != res2) {
1943	if(transition(res1,res2))
1944	p = p + 1.0;
1945	else
1946	q = q + 1.0;
1947	}
1948	skip:;
1949	}
1950
1951
1952	/* Kimura's 2 parameter correction for multiple substitutions */
1953
1954	if(!kimura) {
1955	if (e == 0) {
1956	fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
1957	k = 0.0;
1958	p = 0.0;
1959	q = 0.0;
1960	}
1961	else {
1962	k = (p+q)/e;
1963	if(p > 0.0)
1964	p = p/e;
1965	else
1966	p = 0.0;
1967	if(q > 0.0)
1968	q = q/e;
1969	else
1970	q = 0.0;
1971	}
1972	tmat[m][n] = tmat[n][m] = k;
1973	if(verbose) /* if screen output */
1974	fprintf(tree,
1975	"%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
1976	,(pint)m,(pint)n,k,p,q,e);
1977	}
1978	else {
1979	if (e == 0) {
1980	fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
1981	p = 0.0;
1982	q = 0.0;
1983	}
1984	else {
1985	if(p > 0.0)
1986	p = p/e;
1987	else
1988	p = 0.0;
1989	if(q > 0.0)
1990	q = q/e;
1991	else
1992	q = 0.0;
1993	}
1994
1995	if( ((2.0*p)+q) == 1.0 )
1996	a = 0.0;
1997	else
1998	a = 1.0/(1.0-(2.0*p)-q);
1999
2000	if( q == 0.5 )
2001	b = 0.0;
2002	else
2003	b = 1.0/(1.0-(2.0*q));
2004
2005	/* watch for values going off the scale for the correction. */
2006	if( (a<=0.0) \|\| (b<=0.0) ) {
2007	overspill++;
2008	k = 3.5; /* arbitrary high score */
2009	}
2010	else
2011	k = 0.5log(a) + 0.25log(b);
2012	tmat[m][n] = tmat[n][m] = k;
2013	if(verbose) /* if screen output */
2014	fprintf(tree,
2015	"%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
2016	,(pint)m,(pint)n,k,p,q,e);
2017
2018	}
2019	}
2020	return overspill; /* return the number of off-scale values */
2021	}
2022
2023
2024	sint prot_distance_matrix(FILE *tree)
2025	{
2026	sint m,n;
2027	sint j,i;
2028	sint res1, res2;
2029	sint overspill = 0;
2030	double p,e,k, table_entry;
2031
2032
2033	tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */
2034
2035	if(verbose) {
2036	fprintf(tree,"\n");
2037	fprintf(tree,"\n DIST = percentage divergence (/100)");
2038	fprintf(tree,"\n Length = number of sites used in comparison");
2039	fprintf(tree,"\n\n");
2040	if(tossgaps) {
2041	fprintf(tree,"\n All sites with gaps (in any sequence) deleted");
2042	fprintf(tree,"\n");
2043	}
2044	if(kimura) {
2045	fprintf(tree,"\n Distances up tp 0.75 corrected by Kimura's empirical method:");
2046	fprintf(tree,"\n\n Kimura, M. (1983)");
2047	fprintf(tree," The Neutral Theory of Molecular Evolution.");
2048	fprintf(tree,"\n Page 75. Cambridge University Press, Cambridge, England.");
2049	fprintf(tree,"\n\n");
2050	}
2051	}
2052
2053	for(m=1; m<nseqs; ++m) /* for every pair of sequence */
2054	for(n=m+1; n<=nseqs; ++n) {
2055	p = e = 0.0;
2056	tmat[m][n] = tmat[n][m] = 0.0;
2057	for(i=1; i<=seqlen_array[1]; ++i) {
2058	j = boot_positions[i];
2059	if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */
2060	res1 = seq_array[m][j];
2061	res2 = seq_array[n][j];
2062	if( (res1 == gap_pos1) \|\| (res1 == gap_pos2) \|\|
2063	(res2 == gap_pos1) \|\| (res2 == gap_pos2))
2064	goto skip; /* gap in a seq*/
2065	e = e + 1.0;
2066	if(res1 != res2) p = p + 1.0;
2067	skip:;
2068	}
2069
2070	if(p <= 0.0)
2071	k = 0.0;
2072	else
2073	k = p/e;
2074
2075	/* DES debug */
2076	/* fprintf(stdout,"Seq1=%4d Seq2=%4d k =%7.4f \n",(pint)m,(pint)n,k); */
2077	/* DES debug */
2078
2079	if(kimura) {
2080	if(k < 0.75) { /* use Kimura's formula */
2081	if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) );
2082	}
2083	else {
2084	if(k > 0.930) {
2085	overspill++;
2086	k = 10.0; /* arbitrarily set to 1000% */
2087	}
2088	else {
2089	table_entry = (k*1000.0) - 750.0;
2090	k = (double)dayhoff_pams[(int)table_entry];
2091	k = k/100.0;
2092	}
2093	}
2094	}
2095
2096	tmat[m][n] = tmat[n][m] = k;
2097	if(verbose) /* if screen output */
2098	fprintf(tree,
2099	"%4d vs.%4d DIST = %6.4f; length = %6.0f\n",
2100	(pint)m,(pint)n,k,e);
2101	}
2102	return overspill;
2103	}
2104
2105
2106	void guide_tree(FILE *tree,sint firstseq,sint numseqs)
2107	/*
2108	Routine for producing unrooted NJ trees from seperately aligned
2109	pairwise distances. This produces the GUIDE DENDROGRAMS in
2110	PHYLIP format.
2111	*/
2112	{
2113	static char **standard_tree;
2114	sint i;
2115	float dist;
2116
2117	phylip_phy_tree_file=tree;
2118	verbose = FALSE;
2119	first_seq=firstseq;
2120	last_seq=first_seq+numseqs-1;
2121
2122	if(numseqs==2) {
2123	dist=tmat[firstseq][firstseq+1]/2.0;
2124	fprintf(tree,"(%s:%0.5f,%s:%0.5f);\n",
2125	names[firstseq],dist,names[firstseq+1],dist);
2126	}
2127	else {
2128	standard_tree = (char *) ckalloc( (last_seq-first_seq+2) sizeof (char *) );
2129	for(i=0; i<last_seq-first_seq+2; i++)
2130	standard_tree[i] = (char ) ckalloc( (last_seq-first_seq+2) sizeof(char));
2131
2132	nj_tree(standard_tree,clustal_phy_tree_file);
2133
2134	print_phylip_tree(standard_tree,phylip_phy_tree_file,0);
2135
2136	if(left_branch != NULL) left_branch=ckfree((void *)left_branch);
2137	if(right_branch != NULL) right_branch=ckfree((void *)right_branch);
2138	if(tkill != NULL) tkill=ckfree((void *)tkill);
2139	if(av != NULL) av=ckfree((void *)av);
2140	for (i=1;i<last_seq-first_seq+2;i++)
2141	standard_tree[i]=ckfree((void *)standard_tree[i]);
2142	standard_tree=ckfree((void *)standard_tree);
2143	}
2144	fclose(phylip_phy_tree_file);
2145
2146	}
2147
2148	static Boolean is_ambiguity(char c)
2149	{
2150	int i;
2151	char codes[]="ACGTU";
2152
2153	if(use_ambiguities==TRUE)
2154	{
2155	return FALSE;
2156	}
2157
2158	for(i=0;i<5;i++)
2159	if(amino_acid_codes[c]==codes[i])
2160	return FALSE;
2161
2162	return TRUE;
2163	}
2164

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

Context Navigation

source: trunk/GDE/CLUSTALW/trees.c

Download in other formats: