1 /*
2 * consensus.c
3 *
4 *
5 * Part of TREE-PUZZLE 5.2 (July 2004)
6 *
7 * (c) 2003-2004 by Heiko A. Schmidt, Korbinian Strimmer, and Arndt von Haeseler
8 * (c) 1999-2003 by Heiko A. Schmidt, Korbinian Strimmer,
9 * M. Vingron, and Arndt von Haeseler
10 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
11 *
12 * All parts of the source except where indicated are distributed under
13 * the GNU public licence. See http://www.opensource.org for details.
14 *
15 * ($Id$)
16 *
17 */
18
19
20 #define EXTERN extern
21 #include "consensus.h"
22 #ifdef PARALLEL
23 # include "ppuzzle.h"
24 #endif
25
26 int ZZZ=0;
27 /* fprintf(stderr, "ZZZ: %d (%s:%d)\n", ZZZ++, __FILE__, __LINE__); */
28
29
30 /******************************************************************************/
31 /* functions for computing the consensus tree */
32 /******************************************************************************/
33
34 /* prepare for consensus tree analysis */
initconsensus()35 void initconsensus()
36 {
37 # if ! PARALLEL
38 biparts = new_cmatrix(Maxspc-3, Maxspc);
39 # endif /* PARALLEL */
40
41 if (Maxspc % 32 == 0)
42 splitlength = Maxspc/32;
43 else splitlength = (Maxspc + 32 - (Maxspc % 32))/32;
44 numbiparts = 0; /* no pattern stored so far */
45 maxbiparts = 0; /* no memory reserved so far */
46 splitfreqs = NULL;
47 splitpatterns = NULL;
48 splitsizes = NULL;
49 splitcomptemp = (uli *) calloc((size_t) splitlength, sizeof(uli) );
50 if (splitcomptemp == NULL) maerror("splitcomptemp in initconsensus");
51 } /* initconsensus */
52
53 /******************/
54
55
56 /* de-trueID-ed (HAS) */
57 /* recursive function to get bipartitions */
58 /* traversal should be optimazable (HAS) */
makepart(int i,int curribrnch,ONEEDGE * edge,int * edgeofleaf,cmatrix biparts,int Maxspc)59 void makepart(int i, /* */
60 int curribrnch, /* in: current branch in traversal */
61 ONEEDGE *edge, /* in: tree topology */
62 int *edgeofleaf, /* in: ext. edge list */
63 cmatrix biparts, /* out: split strings, edge by edge */
64 int Maxspc) /* in: number of species in tree */
65 {
66 if ( edge[i].downright == NULL || edge[i].downleft == NULL) {
67 biparts[curribrnch][edge[i].taxon] = '*';
68 return;
69 } else { /* still on inner branch */
70 makepart(edge[i].downleft->numedge, curribrnch, edge, edgeofleaf, biparts, Maxspc);
71 makepart(edge[i].downright->numedge, curribrnch, edge, edgeofleaf, biparts, Maxspc);
72 }
73 } /* makepart */
74
75 /******************/
76
77
78 /* de-trueID-ed (HAS) */
79 /* compute bipartitions of current puzzling step tree */
80 /* traversal should be optimazable (HAS) */
computebiparts(ONEEDGE * edge,int * edgeofleaf,int psteprootleaf,cmatrix biparts,int outgroup,int Maxspc)81 void computebiparts(ONEEDGE *edge, /* in: tree topology */
82 int *edgeofleaf, /* in: ext. edge list */
83 int psteprootleaf, /* virtual root from pstep */
84 cmatrix biparts, /* out: splits */
85 int outgroup, /* in: outgroup of tree */
86 int Maxspc) /* in: number of taxa in tree */
87 {
88 int i, j, curribrnch;
89
90 curribrnch = -1;
91
92 for (i = 0; i < Maxspc - 3; i++)
93 for (j = 0; j < Maxspc; j++)
94 biparts[i][j] = '.';
95
96 for (i = 0; i < Maxbrnch; i++) { /* (HAS) */
97 if (!(edgeofleaf[psteprootleaf] == i
98 || edge[i].downright == NULL
99 || edge[i].downleft == NULL) ) { /* check only inner branches */
100 curribrnch++;
101 makepart(i, curribrnch, edge, edgeofleaf, biparts, Maxspc);
102
103 /* make sure that the root is always a '*' */
104 if (biparts[curribrnch][outgroup] == '.') {
105 for (j = 0; j < Maxspc; j++) {
106 if (biparts[curribrnch][j] == '.')
107 biparts[curribrnch][j] = '*';
108 else
109 biparts[curribrnch][j] = '.';
110 }
111 }
112 }
113 }
114 } /* computebiparts */
115
116 /******************/
117
118
119 /* print out the bipartition n of all different splitpatterns */
fprintfsplit(FILE * fp,uli n,uli * splitpatterns,int splitlength,int Maxspc)120 void fprintfsplit(FILE *fp, /* outputfile stream */
121 uli n, /* split number */
122 uli *splitpatterns, /* splits */
123 int splitlength, /* number of ulis per split */
124 int Maxspc) /* number of taxa */
125 {
126 int i, j; /* counter variables */
127 int col; /* current output column */
128 uli z; /* temporary variable for bit shifting */
129
130 col = 0;
131 for (i = 0; i < splitlength; i++) {
132 z = splitpatterns[n*splitlength + i];
133 for (j = 0; j < 32 && col < Maxspc; j++) {
134 if (col % 10 == 0 && col != 0) fprintf(fp, " ");
135 if (z & 1) fprintf(fp, ".");
136 else fprintf(fp, "*");
137 z = (z >> 1);
138 col++;
139 }
140 }
141 } /* printsplit */
142
143 /******************/
144
145 /* general remarks:
146
147 - every entry in consbiparts is one node of the consensus tree
148 - for each node one has to know which taxa and which other nodes
149 are *directly* descending from it
150 - for every taxon/node number there is a flag that shows
151 whether it descends from the node or not
152 - '0' means that neither a taxon nor another node with the
153 corresponding number decends from the node
154 '1' means that the corresponding taxon descends from the node
155 '2' means that the corresponding node descends from the node
156 '3' means that the corresponding taxon and node descends from the node
157 */
158
159 /* make new entries for new different bipartitions and count frequencies */
160 /* internal use: splitcomptemp */
makenewsplitentries(cmatrix bip,int numspl,uli ** in_splitpatterns,int ** in_splitsizes,uli ** in_splitfreqs,uli * in_numbiparts,uli * in_maxbiparts,int Maxspc)161 void makenewsplitentries(cmatrix bip, /* in: split string vector */
162 int numspl, /* in: no. of new splits */
163 uli **in_splitpatterns,/* io: known compr splits */
164 int **in_splitsizes, /* io: kn. split sizes: '.'*/
165 uli **in_splitfreqs, /* io: kn. split frequences*/
166 uli *in_numbiparts, /* io: no. of splits so far*/
167 uli *in_maxbiparts, /* io: reserved memory */
168 int Maxspc) /* in: no. of species */
169 {
170 int i, j; /* count variables */
171 int bpc; /* split count */
172 int identical; /* split identical flag */
173 int idflag; /* split size identical flag */
174 int bpsize; /* split size: count the '.' */
175 uli nextentry; /* temp vaiable for split amount numbiparts */
176 uli obpc; /* split count - old splits */
177 uli *splitpatternstmp; /* temp pointer, saves derefs. */
178 int *splitsizestmp; /* temp pointer, saves derefs. */
179 uli *splitfreqstmp; /* temp pointer, saves derefs. */
180 uli maxbipartstmp; /* split memory reserved so far */
181
182 splitpatternstmp = *in_splitpatterns;
183 splitsizestmp = *in_splitsizes;
184 splitfreqstmp = *in_splitfreqs;
185 maxbipartstmp = *in_maxbiparts;
186
187 /* where the next entry would be in splitpatterns */
188 nextentry = *in_numbiparts;
189
190 for (bpc = 0; bpc < numspl; bpc++) { /* for every new bipartition */
191 /* convert bipartition into a more compact format */
192 bpsize = 0;
193 for (i = 0; i < splitlength; i++) {
194 splitcomptemp[i] = 0;
195 for (j = 0; j < 32; j++) {
196 splitcomptemp[i] = splitcomptemp[i] >> 1;
197 if (i*32 + j < Maxspc)
198 if (bip[bpc][i*32 + j] == '.') {
199 /* set highest bit */
200 splitcomptemp[i] = (splitcomptemp[i] | 2147483648UL);
201 bpsize++; /* count the '.' */
202 }
203 }
204 }
205 /* compare to the *old* patterns */
206 identical = FALSE;
207 for (obpc = 0; (obpc < *in_numbiparts) && (!identical); obpc++) {
208 /* compare first partition size */
209 if (splitsizestmp[obpc] == bpsize) idflag = TRUE;
210 else idflag = FALSE;
211 /* if size is identical compare whole partition */
212 for (i = 0; (i < splitlength) && idflag; i++)
213 if (splitcomptemp[i] != splitpatternstmp[obpc*splitlength + i])
214 idflag = FALSE;
215 if (idflag) identical = TRUE;
216 }
217 if (identical) { /* if identical increase frequency */
218 splitfreqstmp[2*(obpc-1)]++;
219 } else { /* create new entry */
220 if (nextentry == maxbipartstmp) { /* reserve more memory */
221 maxbipartstmp = maxbipartstmp + 2*Maxspc;
222 splitfreqstmp = (uli *) myrealloc(splitfreqstmp,
223 2*maxbipartstmp * sizeof(uli) );
224 /* 2x: splitfreqstmp contains also an index (sorting!) */
225 if (splitfreqstmp == NULL) maerror("splitfreqstmp in makenewsplitentries");
226 splitpatternstmp = (uli *) myrealloc(splitpatternstmp,
227 splitlength*maxbipartstmp * sizeof(uli) );
228 if (splitpatternstmp == NULL) maerror("splitpatternstmp in makenewsplitentries");
229 splitsizestmp = (int *) myrealloc(splitsizestmp,
230 maxbipartstmp * sizeof(int) );
231 if (splitsizestmp == NULL) maerror("splitsizestmp in makenewsplitentries");
232 }
233 splitfreqstmp[2*nextentry] = 1; /* frequency */
234 splitfreqstmp[2*nextentry+1] = nextentry; /* index for sorting */
235 for (i = 0; i < splitlength; i++)
236 splitpatternstmp[nextentry*splitlength + i] = splitcomptemp[i];
237 splitsizestmp[nextentry] = bpsize;
238 nextentry++;
239 }
240 }
241
242 /* export new values */
243 *in_numbiparts = nextentry;
244 *in_splitpatterns = splitpatternstmp;
245 *in_splitsizes = splitsizestmp;
246 *in_splitfreqs = splitfreqstmp;
247 *in_maxbiparts = maxbipartstmp;
248
249 } /* makenewsplitentries - new */
250
251
252 /***************************************************************************/
253
254 /* copy bipartition n of all different splitpatterns to consbiparts[k] */
copysplit(uli n,uli * splitpatterns,int k,cmatrix consbipartsvec)255 void copysplit(uli n, uli *splitpatterns, int k, cmatrix consbipartsvec)
256 {
257 int i, j, col;
258 uli z;
259
260 col = 0;
261 for (i = 0; i < splitlength; i++) {
262 z = splitpatterns[n*splitlength + i];
263 for (j = 0; j < 32 && col < Maxspc; j++) {
264 if (z & 1) consbipartsvec[k][col] = '1';
265 else consbipartsvec[k][col] = '0';
266 z = (z >> 1);
267 col++;
268 consbipartsvec[k][Maxspc] = '\0';
269 }
270 }
271 } /* copysplit */
272
273 /******************/
274
275
276
277 /* <consbipartsvec> and the quartet topologies from the ML step. Values are */
278 /* stored in <qsupparr[splitnum]> */
279 /* missing parameter: quartetinfo */
quartsupport(int splitnum,cmatrix consbipartsvec,qsupportarr_t * qsupparr)280 void quartsupport(int splitnum, cmatrix consbipartsvec, qsupportarr_t *qsupparr)
281
282 {
283 int n;
284 int *clusterA;
285 int *clusterB;
286 int clustA = 0;
287 int clustB = 0;
288
289 int a1, a2, b1, b2;
290 int aa1, aa2, bb1, bb2;
291 uli fullres_pro = 0;
292 uli fullres_con = 0;
293 uli partres_pro = 0;
294 uli partres_con = 0;
295 uli unres = 0;
296 uli missing = 0;
297 uli qsum = 0;
298
299 cvector splitstr;
300 unsigned char tmpweight;
301
302 clusterA = calloc((size_t) Maxspc, (sizeof(int)));
303 if (clusterA == NULL) maerror("clusterA in quartsupport");
304 clusterB = calloc((size_t) Maxspc, (sizeof(int)));
305 if (clusterA == NULL) maerror("clusterB in quartsupport");
306
307 splitstr = consbipartsvec[splitnum];
308
309 for(n=0; n<Maxspc; n++) {
310 if (splitstr[n] == '0') clusterA[clustA++] = n;
311 else clusterB[clustB++] = n;
312 }
313
314 if ((clustA <= 1) || (clustB <= 1)) {
315 /* split in outer edge */
316 fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR HS1 TO DEVELOPERS\n\n\n");
317 # if PARALLEL
318 PP_Finalize();
319 # endif
320 exit(1);
321 }
322
323 for (aa1 = 0; aa1 < clustA-1; aa1++)
324 for (aa2 = aa1+1; aa2 < clustA; aa2++)
325 for (bb1 = 0; bb1 < clustB-1; bb1++)
326 for (bb2 = bb1+1; bb2 < clustB; bb2++) {
327 a1 = clusterA[aa1];
328 a2 = clusterA[aa2];
329 b1 = clusterB[bb1];
330 b2 = clusterB[bb2];
331
332 if (a1 < b1) {
333 if (b1 < a2) {
334 if (a2 < b2) { /* a1<b1<a2<b2 */
335 tmpweight = readquartet(a1, b1, a2, b2);
336 qsum++;
337 switch (tmpweight) {
338 case 1: fullres_con++; break;
339 case 2: fullres_pro++; break;
340 case 3: partres_pro++; break;
341 case 4: fullres_con++; break;
342 case 5: partres_con++; break;
343 case 6: partres_pro++; break;
344 case 7: unres++; break;
345 case 0: missing++; break;
346 }
347 } else { /* a1<b1<b2<a2 */
348 tmpweight = readquartet(a1, b1, b2, a2);
349 qsum++;
350 switch (tmpweight) {
351 case 1: fullres_con++; break;
352 case 2: fullres_con++; break;
353 case 3: partres_con++; break;
354 case 4: fullres_pro++; break;
355 case 5: partres_pro++; break;
356 case 6: partres_pro++; break;
357 case 7: unres++; break;
358 case 0: missing++; break;
359 }
360 }
361 } else { /* a1<a2<b1<b2 */
362 tmpweight = readquartet(a1, a2, b1, b2);
363 qsum++;
364 switch (tmpweight) {
365 case 1: fullres_pro++; break;
366 case 2: fullres_con++; break;
367 case 3: partres_pro++; break;
368 case 4: fullres_con++; break;
369 case 5: partres_pro++; break;
370 case 6: partres_con++; break;
371 case 7: unres++; break;
372 case 0: missing++; break;
373 }
374 }
375 } else {
376 if (a1 < b2) {
377 if (b2 < a2) { /* b1<a1<b2<a2 */
378 tmpweight = readquartet(b1, a1, b2, a2);
379 qsum++;
380 switch (tmpweight) {
381 case 1: fullres_con++; break;
382 case 2: fullres_pro++; break;
383 case 3: partres_pro++; break;
384 case 4: fullres_con++; break;
385 case 5: partres_con++; break;
386 case 6: partres_pro++; break;
387 case 7: unres++; break;
388 case 0: missing++; break;
389 }
390 } else { /* b1<a1<a2<b2 */
391 tmpweight = readquartet(b1, a1, a2, b2);
392 qsum++;
393 switch (tmpweight) {
394 case 1: fullres_con++; break;
395 case 2: fullres_con++; break;
396 case 3: partres_con++; break;
397 case 4: fullres_pro++; break;
398 case 5: partres_pro++; break;
399 case 6: partres_pro++; break;
400 case 7: unres++; break;
401 case 0: missing++; break;
402 }
403 }
404 } else { /* b1<b2<a1<a2 */
405 tmpweight = readquartet(b1, b2, a1, a2);
406 qsum++;
407 switch (tmpweight) {
408 case 1: fullres_pro++; break;
409 case 2: fullres_con++; break;
410 case 3: partres_pro++; break;
411 case 4: fullres_con++; break;
412 case 5: partres_pro++; break;
413 case 6: partres_con++; break;
414 case 7: unres++; break;
415 case 0: missing++; break;
416 }
417 }
418 }
419 }
420
421 qsupparr[splitnum].fullres_pro = fullres_pro;
422 qsupparr[splitnum].fullres_con = fullres_con;
423 qsupparr[splitnum].partres_pro = partres_pro;
424 qsupparr[splitnum].partres_con = partres_con;
425 qsupparr[splitnum].unres= unres;
426 qsupparr[splitnum].missing = missing;
427 qsupparr[splitnum].qsum = qsum;
428
429 free(clusterA);
430 free(clusterB);
431 } /* quartsupport */
432
433 /******************/
434
435
436
437 /* compute majority rule consensus tree */
makeconsensus(uli maxnumtrees,int qsupp_optn)438 void makeconsensus(uli maxnumtrees, int qsupp_optn)
439 {
440 int i, j, k, n, size, subnode;
441 char chari, charj;
442 int done;
443 int count01, count10, count00, count11, scount;
444
445 /* sort bipartition frequencies */
446 qsort(splitfreqs, numbiparts, 2*sizeof(uli), ulicmp);
447 /* how many bipartitions are included in the consensus tree */
448 consincluded = 0;
449 for (i = 0; (uli)i < numbiparts && i == consincluded; i++) {
450 if (2*splitfreqs[2*i] > maxnumtrees) consincluded = i + 1;
451 }
452
453 /* collect all info about majority rule consensus tree */
454 /* the +1 is due to the edge with the root */
455 /* consincluded changed to (Maxspc-3) to collect sub-50%-splits */
456 /* a maximum of (Maxspc-3) is possible in a Maxspc-tree (HAS) */
457 consconfid = new_ivector( (Maxspc-3) + 1);
458 conssizes = new_ivector(2*(Maxspc-3) + 2);
459 consbiparts = new_cmatrix( (Maxspc-3) + 1, Maxspc + 1);
460
461 /* old routine (HAS)
462 consconfid = new_ivector(consincluded + 1);
463 conssizes = new_ivector(2*consincluded + 2);
464 consbiparts = new_cmatrix(consincluded + 1, Maxspc);
465 */
466
467 for (i = 0; i < consincluded; i++) {
468 /* copy partition to consbiparts */
469 copysplit(splitfreqs[2*i+1], splitpatterns, i, consbiparts);
470 /* frequency in percent (rounded to integer) */
471 consconfid[i] = (int) floor(100.0*splitfreqs[2*i]/maxnumtrees + 0.5);
472 /* size of partition */
473 conssizes[2*i] = splitsizes[splitfreqs[2*i+1]];
474 conssizes[2*i+1] = i;
475 }
476
477 consfifty = consincluded;
478 done = FALSE;
479 i = consincluded;
480
481 /* until maximum splits (Maxspc-3) incorporated or */
482 /* first incongruence found */
483 while ((i < (Maxspc-3)) && !done) {
484
485 /* copy partition to consbiparts */
486 copysplit(splitfreqs[2*i+1], splitpatterns, i, consbiparts);
487
488 for (scount = 0; (scount<i) && (!done); scount++) {
489 count01 = 0;
490 count10 = 0;
491 count00 = 0;
492 count11 = 0;
493
494 for (n = 0; (n<Maxspc) && (!done); n++) {
495 if (consbiparts[i][n] != consbiparts[scount][n]) {
496 if (consbiparts[i][n] == '0') count01++;
497 else count10++;
498 } else {
499 if (consbiparts[i][n] == '0') count00++;
500 else count11++;
501 }
502 if ((count01 > 0) && (count10 > 0) && (count00 > 0) && (count11 > 0)) done = TRUE;
503 }
504 }
505 if (!done) {
506 /* frequency in percent (rounded to integer) */
507 consconfid[i] = (int) floor(100.0*splitfreqs[2*i]/maxnumtrees+ 0.5);
508 /* size of partition */
509 conssizes[2*i] = splitsizes[splitfreqs[2*i+1]];
510 conssizes[2*i+1] = i;
511 i++;
512 }
513 }
514
515 /* last split has same percentage as the first unincorporated, remove */
516 if (((uli)(i + 1) < numbiparts) && (splitfreqs[2*i] == splitfreqs[2*(i+1)]))
517 done = TRUE;
518
519 /* remove all splits with same percentage as first incongruence */
520 if (done) {
521 while ((i > 0) && (splitfreqs[2*i] == splitfreqs[2*(i-1)]))
522 i--;
523 }
524
525 if (conssub50_optn) {
526 consincluded = i;
527 conscongruent = i;
528 } else {
529 conscongruent = i;
530 }
531
532 for (i = 0; i < Maxspc; i++) consbiparts[consincluded][i] = '1';
533 consbiparts[consincluded][outgroup] = '0';
534 consconfid[consincluded] = 100;
535 conssizes[2*consincluded] = Maxspc - 1;
536 conssizes[2*consincluded + 1] = consincluded;
537
538 /* sort bipartitions according to cluster size */
539 qsort(conssizes, consincluded + 1, 2*sizeof(int), intcmp);
540
541 qsupportarr = malloc((size_t) (sizeof(qsupportarr_t) * consincluded));
542 if (qsupp_optn) {
543 for (i = 0; i < consincluded; i++) {
544 quartsupport(i, consbiparts, qsupportarr);
545 }
546 }
547
548 /******************************/
549 /* reconstruct consensus tree */
550 /******************************/
551 for (i = 0; i < consincluded; i++) { /* try every node */
552 size = conssizes[2*i]; /* size of current node */
553 for (j = i + 1; j < consincluded + 1; j++) {
554
555 /* compare only with nodes with more descendants */
556 if (size == conssizes[2*j]) continue;
557
558 /* check whether node i is a subnode of j */
559 subnode = FALSE;
560 for (k = 0; k < Maxspc && !subnode; k++) {
561 chari = consbiparts[ conssizes[2*i+1] ][k];
562 if (chari != '0') {
563 charj = consbiparts[ conssizes[2*j+1] ][k];
564 if (chari == charj || charj == '3') subnode = TRUE;
565 }
566 }
567
568 /* if i is a subnode of j change j accordingly */
569 if (subnode) {
570 /* remove subnode i from j */
571 for (k = 0; k < Maxspc; k++) {
572 chari = consbiparts[ conssizes[2*i+1] ][k];
573 if (chari != '0') {
574 charj = consbiparts[ conssizes[2*j+1] ][k];
575 if (chari == charj)
576 consbiparts[ conssizes[2*j+1] ][k] = '0';
577 else if (charj == '3') {
578 if (chari == '1')
579 consbiparts[ conssizes[2*j+1] ][k] = '2';
580 else if (chari == '2')
581 consbiparts[ conssizes[2*j+1] ][k] = '1';
582 else {
583 /* Consensus tree [1] */
584 fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR H TO DEVELOPERS\n\n\n");
585 # if PARALLEL
586 PP_Finalize();
587 # endif
588 exit(1);
589 }
590 } else {
591 /* Consensus tree [2] */
592 fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR I TO DEVELOPERS\n\n\n");
593 # if PARALLEL
594 PP_Finalize();
595 # endif
596 exit(1);
597 }
598 }
599 }
600 /* add link to subnode i in node j */
601 charj = consbiparts[ conssizes[2*j+1] ][ conssizes[2*i+1] ];
602 if (charj == '0')
603 consbiparts[ conssizes[2*j+1] ][ conssizes[2*i+1] ] = '2';
604 else if (charj == '1')
605 consbiparts[ conssizes[2*j+1] ][ conssizes[2*i+1] ] = '3';
606 else {
607 /* Consensus tree [3] */
608 fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR J TO DEVELOPERS\n\n\n");
609 # if PARALLEL
610 PP_Finalize();
611 # endif
612 exit(1);
613 }
614 }
615 }
616 }
617 } /* makeconsensus */
618
619 /******************/
620
621
622 /* prototype for recursion */
623 void writenode(FILE *treefile, /* in: output stream */
624 int node, /* current node */
625 int qsupp_optn, /* 'print quartet support' flag */
626 qsupportarr_t *qsupparr, /* quartet support values */
627 int *column); /* current line position */
628
629 /* write node (writeconsensustree) */
630 /* missing: column */
writenode(FILE * treefile,int node,int qsupp_optn,qsupportarr_t * qsupparr,int * column)631 void writenode(FILE *treefile, /* in: output stream */
632 int node, /* current node */
633 int qsupp_optn, /* 'print quartet support' flag */
634 qsupportarr_t *qsupparr, /* quartet support values */
635 int *column) /* current line position */
636 {
637 int i, first;
638
639 fprintf(treefile, "(");
640 (*column)++;
641 /* write descending nodes */
642 first = TRUE;
643 for (i = 0; i < Maxspc; i++) {
644 if (consbiparts[node][i] == '2' ||
645 consbiparts[node][i] == '3') {
646 if (first) first = FALSE;
647 else {
648 fprintf(treefile, ",");
649 (*column)++;
650 }
651 if ((*column) > 60) {
652 (*column) = 2;
653 fprintf(treefile, "\n");
654 }
655 /* write node i */
656 writenode(treefile, i, qsupp_optn, qsupparr, column);
657
658 if (qsupp_optn) {
659 /* reliability value as internal label */
660 fprintf(treefile, "%d", consconfid[i]);
661 (*column) += 3;
662 #if 0
663 fprintf(treefile, "\"bs=%d|", consconfid[i]);
664 (*column) += 7;
665
666 fprintf(treefile, "f+=%.0f|", (100.0*qsupparr[i].fullres_pro/qsupparr[i].qsum));
667 (*column) += 7;
668 fprintf(treefile, "p+=%.0f|", (100.0*qsupparr[i].partres_pro/qsupparr[i].qsum));
669 (*column) += 7;
670 fprintf(treefile, "f-=%.0f|", (100.0*qsupparr[i].fullres_con/qsupparr[i].qsum));
671 (*column) += 7;
672 fprintf(treefile, "p-=%.0f|", (100.0*qsupparr[i].partres_con/qsupparr[i].qsum));
673 (*column) += 7;
674 fprintf(treefile, "ur=%.0f", (100.0*qsupparr[i].unres/qsupparr[i].qsum));
675 (*column) += 6;
676 if (qsupparr[i].missing > 0) {
677 fprintf(treefile, "|md=%.0f", (100.0*qsupparr[i].missing/qsupparr[i].qsum));
678 (*column) += 7;
679 }
680 fprintf(treefile, "/%ld\"", qsupparr[i].qsum);
681 (*column) += 5;
682 #endif
683 } else {
684 /* reliability value as internal label */
685 fprintf(treefile, "%d", consconfid[i]);
686 (*column) += 3;
687 }
688 }
689 }
690 /* write descending taxa */
691 for (i = 0; i < Maxspc; i++) {
692 if (consbiparts[node][i] == '1' ||
693 consbiparts[node][i] == '3') {
694 if (first) first = FALSE;
695 else {
696 fprintf(treefile, ",");
697 (*column)++;
698 }
699 if ((*column) > 60) {
700 (*column) = 2;
701 fprintf(treefile, "\n");
702 }
703 (*column) += fputid(treefile, i);
704 }
705 }
706 fprintf(treefile, ")");
707 (*column)++;
708 } /* writenode */
709
710 /******************/
711
712
713
714 /* write consensus tree */
715 /* missing: consbiparts, consincluded, consconfid, */
writeconsensustree(FILE * treefile,int qsupp_optn,qsupportarr_t * qsupparr)716 void writeconsensustree(FILE *treefile, /* in: output stream */
717 int qsupp_optn, /* 'print quartsupp' flag */
718 qsupportarr_t *qsupparr) /* quartet support values */
719 {
720 int i, first;
721 int column;
722
723
724 column = 1;
725 fprintf(treefile, "(");
726 column += fputid(treefile, outgroup) + 2;
727 fprintf(treefile, ",");
728 /* write descending nodes */
729 first = TRUE;
730 for (i = 0; i < Maxspc; i++) {
731 if (consbiparts[consincluded][i] == '2' || consbiparts[consincluded][i] == '3') {
732 if (first) first = FALSE;
733 else {
734 fprintf(treefile, ","); /* first (outgroup) has ',' already */
735 column++;
736 }
737 if (column > 60) {
738 column = 2;
739 fprintf(treefile, "\n");
740 }
741 /* write node i */
742 writenode(treefile, i, qsupp_optn, qsupparr, &column);
743
744 if (qsupp_optn) {
745 /* reliability value as internal label */
746 fprintf(treefile, "%d", consconfid[i]);
747 column += 3;
748 #if 0
749 fprintf(treefile, "\"bs=%d|", consconfid[i]);
750 column += 7;
751
752 fprintf(treefile, "f+=%.0f|", (100.0*qsupparr[i].fullres_pro/qsupparr[i].qsum));
753 column += 7;
754 fprintf(treefile, "p+=%.0f|", (100.0*qsupparr[i].partres_pro/qsupparr[i].qsum));
755 column += 7;
756 fprintf(treefile, "f-=%.0f|", (100.0*qsupparr[i].fullres_con/qsupparr[i].qsum));
757 column += 7;
758 fprintf(treefile, "p-=%.0f|", (100.0*qsupparr[i].partres_con/qsupparr[i].qsum));
759 column += 7;
760 fprintf(treefile, "ur=%.0f", (100.0*qsupparr[i].unres/qsupparr[i].qsum));
761 column += 6;
762 if (qsupparr[i].missing > 0) {
763 fprintf(treefile, "|md=%.0f", (100.0*qsupparr[i].missing/qsupparr[i].qsum));
764 column += 7;
765 }
766 fprintf(treefile, "/%ld\"", qsupparr[i].qsum);
767 column += 5;
768 #endif
769 } else {
770 /* reliability value as internal label */
771 fprintf(treefile, "%d", consconfid[i]);
772 column += 3;
773 }
774 }
775 }
776 /* write descending taxa */
777 for (i = 0; i < Maxspc; i++) {
778 if (consbiparts[consincluded][i] == '1' ||
779 consbiparts[consincluded][i] == '3') {
780 if (first) first = FALSE;
781 else {
782 fprintf(treefile, ",");
783 column++;
784 }
785 if (column > 60) {
786 column = 2;
787 fprintf(treefile, "\n");
788 }
789 column += fputid(treefile, i);
790 }
791 }
792 fprintf(treefile, ");\n");
793 } /* writeconsensustree */
794
795 /******************/
796
797
798 /* prototype for recursion */
799 void nodecoordinates(int node);
800
801 /* establish node coordinates (plotconsensustree) */
nodecoordinates(int node)802 void nodecoordinates(int node)
803 {
804 int i, ymin, ymax, xcoordinate;
805
806 /* first establish coordinates of descending nodes */
807 for (i = 0; i < Maxspc; i++) {
808 if (consbiparts[node][i] == '2' ||
809 consbiparts[node][i] == '3')
810 nodecoordinates(i);
811 }
812
813 /* then establish coordinates of descending taxa */
814 for (i = 0; i < Maxspc; i++) {
815 if (consbiparts[node][i] == '1' ||
816 consbiparts[node][i] == '3') {
817 /* y-coordinate of taxon i */
818 ycortax[i] = ytaxcounter;
819 ytaxcounter = ytaxcounter - 2;
820 }
821 }
822
823 /* then establish coordinates of this node */
824 ymin = 2*Maxspc - 2;
825 ymax = 0;
826 xcoordinate = 0;
827 for (i = 0; i < Maxspc; i++) {
828 if (consbiparts[node][i] == '2' ||
829 consbiparts[node][i] == '3') {
830 if (ycor[i] > ymax) ymax = ycor[i];
831 if (ycor[i] < ymin) ymin = ycor[i];
832 if (xcor[i] > xcoordinate) xcoordinate = xcor[i];
833 }
834 }
835 for (i = 0; i < Maxspc; i++) {
836 if (consbiparts[node][i] == '1' ||
837 consbiparts[node][i] == '3') {
838 if (ycortax[i] > ymax) ymax = ycortax[i];
839 if (ycortax[i] < ymin) ymin = ycortax[i];
840 }
841 }
842 ycormax[node] = ymax;
843 ycormin[node] = ymin;
844 ycor[node] = (int) floor(0.5*(ymax + ymin) + 0.5);
845 if (xcoordinate == 0) xcoordinate = 9;
846 xcor[node] = xcoordinate + 4;
847 } /* nodecoordinates */
848
849 /******************/
850
851
852 /* prototype for recursion */
853 void drawnode(int node, int xold);
854
855 /* drawnode (plotconsensustree) */
drawnode(int node,int xold)856 void drawnode(int node, int xold)
857 {
858 int i, j;
859 char buf[4];
860
861 /* first draw vertical line */
862 for (i = ycormin[node] + 1; i < ycormax[node]; i++)
863 treepict[xcor[node]][i] = ':';
864
865 /* then draw descending nodes */
866 for (i = 0; i < Maxspc; i++) {
867 if (consbiparts[node][i] == '2' ||
868 consbiparts[node][i] == '3')
869 drawnode(i, xcor[node]);
870 }
871
872 /* then draw descending taxa */
873 for (i = 0; i < Maxspc; i++) {
874 if (consbiparts[node][i] == '1' ||
875 consbiparts[node][i] == '3') {
876 treepict[xcor[node]][ycortax[i]] = ':';
877 for (j = xcor[node] + 1; j < xsize-10; j++)
878 treepict[j][ycortax[i]] = '-';
879 for (j = 0; j < 10; j++)
880 treepict[xsize-10+j][ycortax[i]] = Identif[i][j];
881 }
882 }
883
884 /* then draw internal edge with consensus value */
885 treepict[xold][ycor[node]] = ':';
886 treepict[xcor[node]][ycor[node]] = ':';
887 for (i = xold + 1; i < xcor[node]-3; i++)
888 treepict[i][ycor[node]] = '-';
889 sprintf(buf, "%d", consconfid[node]);
890 if (consconfid[node] == 100) {
891 treepict[xcor[node]-3][ycor[node]] = buf[0];
892 treepict[xcor[node]-2][ycor[node]] = buf[1];
893 treepict[xcor[node]-1][ycor[node]] = buf[2];
894 } else {
895 treepict[xcor[node]-3][ycor[node]] = '-';
896 treepict[xcor[node]-2][ycor[node]] = buf[0];
897 treepict[xcor[node]-1][ycor[node]] = buf[1];
898 }
899 } /* drawnode */
900
901 /******************/
902
903
904 /* plot consensus tree */
plotconsensustree(FILE * plotfp)905 void plotconsensustree(FILE *plotfp)
906 {
907 int i, j, yroot, startree;
908
909 /* star tree or no star tree */
910 if (consincluded == 0) {
911 startree = TRUE;
912 consincluded = 1; /* avoids problems with malloc */
913 } else
914 startree = FALSE;
915
916 /* memory for x-y-coordinates of each bipartition */
917 xcor = new_ivector(consincluded);
918 ycor = new_ivector(consincluded);
919 ycormax = new_ivector(consincluded);
920 ycormin = new_ivector(consincluded);
921 if (startree) consincluded = 0; /* avoids problems with malloc */
922
923 /* y-coordinates of each taxon */
924 ycortax = new_ivector(Maxspc);
925 ycortax[outgroup] = 0;
926
927 /* establish coordinates */
928 ytaxcounter = 2*Maxspc - 2;
929
930 /* first establish coordinates of descending nodes */
931 for (i = 0; i < Maxspc; i++) {
932 if (consbiparts[consincluded][i] == '2' ||
933 consbiparts[consincluded][i] == '3')
934 nodecoordinates(i);
935 }
936
937 /* then establish coordinates of descending taxa */
938 for (i = 0; i < Maxspc; i++) {
939 if (consbiparts[consincluded][i] == '1' ||
940 consbiparts[consincluded][i] == '3') {
941 /* y-coordinate of taxon i */
942 ycortax[i] = ytaxcounter;
943 ytaxcounter = ytaxcounter - 2;
944 }
945 }
946
947 /* then establish length of root edge and size of whole tree */
948 yroot = 0;
949 xsize = 0;
950 for (i = 0; i < Maxspc; i++) {
951 if (consbiparts[consincluded][i] == '2' ||
952 consbiparts[consincluded][i] == '3') {
953 if (ycor[i] > yroot) yroot = ycor[i];
954 if (xcor[i] > xsize) xsize = xcor[i];
955 }
956 }
957 for (i = 0; i < Maxspc; i++) {
958 if (consbiparts[consincluded][i] == '1' ||
959 consbiparts[consincluded][i] == '3') {
960 if (ycortax[i] > yroot) yroot = ycortax[i];
961 }
962 }
963 if (xsize == 0) xsize = 9;
964 /* size in x direction inclusive one blank on the left */
965 xsize = xsize + 6;
966
967 /* change all x-labels so that (0,0) is down-left */
968 for (i = 0; i < consincluded; i++)
969 xcor[i] = xsize-1-xcor[i];
970
971 /* draw tree */
972 treepict = new_cmatrix(xsize, 2*Maxspc-1);
973 for (i = 0; i < xsize; i++)
974 for (j = 0; j < 2*Maxspc-1; j++)
975 treepict[i][j] = ' ';
976
977 /* draw root */
978 for (i = 1; i < yroot; i++)
979 treepict[1][i] = ':';
980 treepict[1][0] = ':';
981 for (i = 2; i < xsize - 10; i++)
982 treepict[i][0] = '-';
983 for (i = 0; i < 10; i++)
984 treepict[xsize-10+i][0] = Identif[outgroup][i];
985
986 /* then draw descending nodes */
987 for (i = 0; i < Maxspc; i++) {
988 if (consbiparts[consincluded][i] == '2' ||
989 consbiparts[consincluded][i] == '3')
990 drawnode(i, 1);
991 }
992
993 /* then draw descending taxa */
994 for (i = 0; i < Maxspc; i++) {
995 if (consbiparts[consincluded][i] == '1' ||
996 consbiparts[consincluded][i] == '3') {
997 treepict[1][ycortax[i]] = ':';
998 for (j = 2; j < xsize-10; j++)
999 treepict[j][ycortax[i]] = '-';
1000 for (j = 0; j < 10; j++)
1001 treepict[xsize-10+j][ycortax[i]] = Identif[i][j];
1002 }
1003 }
1004
1005 /* plot tree */
1006 for (i = 2*Maxspc-2; i > -1; i--) {
1007 for (j = 0; j < xsize; j++)
1008 fputc(treepict[j][i], plotfp);
1009 fputc('\n', plotfp);
1010 }
1011
1012 free_ivector(xcor);
1013 free_ivector(ycor);
1014 free_ivector(ycormax);
1015 free_ivector(ycormin);
1016 free_ivector(ycortax);
1017 free_cmatrix(treepict);
1018 } /* plotconsensustree */
1019
1020
1021
1022
1023