Context Navigation

source: branches/nameserver/GDE/CLUSTALW/trees.c

Visit:

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

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

Download in other formats: