1 /* ======================================================= *
2 * Copyright 1998-2008 Stephen C. Grubb *
3 * http://ploticus.sourceforge.net *
4 * Covered by GPL; see the file ./Copyright for details. *
5 * ======================================================= */
6
7 /* PROC TREE - render a tree diagram from a newick file */
8
9 #include "pl.h"
10 #define MAXNEWICKTOKENS 1000 /* max # tokens in input newick file.. indiv. tokens include: labels, commas, colons, distances */
11 #define HALFPI 1.5707963
12
13 /* static int angpos(); */
14
15 int
PLP_tree()16 PLP_tree()
17 {
18 char attr[NAMEMAXLEN], *line, *lineval;
19 int lvp, first;
20
21 char *newickfile, *symbol, *linedetails, *textdetails, *rootsym;
22 double adjx, adjy;
23 int i, j, align, c, nrows, ntok, newstring;
24
25 char *label[MAXNEWICKTOKENS];
26 int nest[MAXNEWICKTOKENS];
27 char func[MAXNEWICKTOKENS];
28 double xval[MAXNEWICKTOKENS], yval[MAXNEWICKTOKENS];
29 double dist[MAXNEWICKTOKENS];
30 char *buf;
31
32 double accum;
33 char *toklist[MAXNEWICKTOKENS];
34 int nestlevel, tagcount, naccum;
35 FILE *fp;
36
37 char symcode[80];
38 double radius, x, y, y1, y2, yofs;
39 int squaredoff, mustparen, debug;
40
41 newickfile = "";
42 symbol = "";
43 linedetails = "color=blue width=0.5";
44 textdetails = "";
45 rootsym = "shape=pixcircle color=black radius=0.02 style=filled";
46 squaredoff = 1;
47 mustparen = 1; /* if 1, the 1st char in newick file must be an open paren */
48 debug = 0;
49 TDH_setvar( "TREEDONE", "0" );
50
51 first = 1;
52 while( 1 ) {
53 line = getnextattr( first, attr, &lvp );
54 if( line == NULL ) break;
55 first = 0;
56 lineval = &line[lvp];
57
58 if( strcmp( attr, "newickfile" )==0 ) newickfile = lineval;
59 else if( strcmp( attr, "symbol" )==0 ) symbol = lineval;
60 else if( strcmp( attr, "linedetails" )==0 ) linedetails = lineval;
61 else if( strcmp( attr, "textdetails" )==0 ) textdetails = lineval;
62 else if( strcmp( attr, "rootsym" )==0 ) rootsym = lineval;
63 else if( strcmp( attr, "squaredoff" )==0 ) squaredoff = getyn( lineval );
64 else if( strcmp( attr, "mustparen" )==0 ) mustparen = getyn( lineval );
65 else if( strcmp( attr, "debug" )==0 ) debug = getyn( lineval );
66 else Eerr( 1, "attribute not recognized", attr );
67 }
68
69 if( newickfile[0] == '\0' ) return( Eerr( 673, "newickfile not specified", "" ));
70 if( !scalebeenset() ) return( Eerr( 51, "No scaled plotting area has been defined yet w/ proc areadef", "" ) );
71
72
73 /* read in the newick file..
74 *
75 * Newick grammar taken from: http://evolution.genetics.washington.edu/phylip/newick_doc.html ....
76 * ... exception: to print a ' in a label use eg. 'steve\'s kludge'
77 *
78 * Example:
79 * (NOD/LtJ:0.00000,129S1SvImJ:0.00000,((KK/HlJ:0.40085,((CAST/EiJ:0.70121,(MOLF/EiJ:0.18484,PWD/PhJ:0.03214):0.55524):0.49231,
80 * (WSB/EiJ:0.73653):0.21783):0.21668,AKR/J:0.01357,(NZW/LacJ:0.00000,BTBRT+tf/J:0.02142,FVB/NJ:0.00000,C3H/HeJ:0.00000,
81 * BALB/cByJ:0.00000,A/J:0.00000):0.01392):0.01188,DBA/2J:0.01071)[0.3333];
82 */
83
84 fp = fopen( newickfile, "r" );
85 if( fp == NULL ) return( Eerr( 674, "cannot open newickfile", newickfile ));
86
87 buf = PL_bigbuf;
88
89 /* read in the newick format data and parse into tokens.. */
90 ntok = 0;
91 newstring = 0;
92 for( i = 0; i < (MAXBIGBUF-2); i++ ) {
93 if( ntok > MAXNEWICKTOKENS-2 ) { fclose(fp); return( Eerr( 675, "sorry, too many tokens in newickfile, no tree produced", "MAXNEWICKTOKENS" )); }
94 c = getc( fp );
95 if( i == 0 && c != '(' && mustparen ) { fclose(fp); return( Eerr( 678, "first char in newickfile not a paren, quitting", "" )); }
96 if( c == EOF ) break;
97 else if( isspace( c )) { i--; continue; }
98 else if( c == ')' || c == '(' || c == ',' || c == ';' ) {
99 buf[i++] = '\0'; buf[i++] = c; buf[i] = '\0';
100 toklist[ntok++] = &buf[i-1];
101 newstring = 1;
102 }
103 else if( c == ':' ) {
104 buf[i++] = '\0';
105 buf[i] = c;
106 toklist[ntok++] = &buf[i];
107 newstring = 0;
108 }
109 else if( c == '[' ) {
110 i--;
111 while( 1 ) { /* skip [comments] ... */
112 c = getc( fp );
113 if( c == ']' || c == EOF ) break;
114 }
115 }
116 else if( c == '\'' ) { /* handle 'quoted labels' ... */
117 toklist[ntok++] = &buf[i];
118 for( j = 0; ; j++, i++ ) {
119 c = getc( fp );
120 if( ( c == '\'' && buf[i] != '\\' ) || c == EOF ) break;
121 else if( c == ' ' ) buf[i] = '_';
122 else if( j == 0 && GL_member( c, ":()" ) ) buf[i] = '_';
123 else buf[i] = c;
124 }
125 i--;
126 }
127 else {
128 if( newstring ) toklist[ntok++] = &buf[i];
129 newstring = 0;
130 buf[i] = c;
131 }
132 }
133
134 if( debug ) {
135 fprintf( PLS.diagfp, "\nTOKENS: " );
136 for( i = 0; i < ntok; i++ ) fprintf( PLS.diagfp, "[%s]", toklist[i] );
137 fprintf( PLS.diagfp, "\n(N=%d)\n\n", ntok );
138 }
139
140 fclose( fp );
141
142
143 /* now do some processing on the token list.. */
144 nestlevel = 0;
145 nrows = 0;
146 tagcount = 0;
147 for( i = 0; i < ntok; i++ ) {
148 if( debug ) fprintf( PLS.diagfp, "[%s]", toklist[i] );
149
150 if( strcmp( toklist[i], ";" )==0 ) break;
151
152 else if( strcmp( toklist[i], "(" )==0 ) {
153 label[nrows] = ""; /* was: = toklist[i]; */
154 func[nrows] = '+';
155 nestlevel++;
156 xval[nrows] = 0.0;
157 yval[nrows] = 0.0;
158 nest[nrows++] = nestlevel;
159 continue;
160 }
161 else if( strcmp( toklist[i], ")" )==0 ) {
162 if( nestlevel <= 0 ) return( Eerr( 677, "too many closing parens in newickfile, no tree produced", "" ));
163 nest[nrows] = nestlevel;
164 func[nrows] = '-';
165 nestlevel--;
166 label[nrows] = ""; /* was: = toklist[i]; */
167 xval[nrows] = 0.0;
168 yval[nrows] = 0.0;
169 /* now look ahead to next toks.. there could be an internal node label and/or a distance.. */
170 if( GL_smember( toklist[i+1], "( ) , ;" ));
171 else if( toklist[i+1][0] != ':' ) { label[nrows] = toklist[i+1]; i++; } /* capture internal node label.. */
172 if( GL_smember( toklist[i+1], "( ) , ;" ));
173 else if( toklist[i+1][0] == ':' ) {
174 dist[nrows] = atof( &toklist[i+1][1] );
175 for( j = nrows; j >= 0; j-- ) {
176 /* if( nest[j] < nest[nrows] ) break; */
177 if( nest[j] == nest[nrows] && func[j] == '+' ) break;
178 xval[j] += dist[nrows];
179 if( func[j] == 'T' && nest[j] == nest[nrows] ) { accum += yval[j]; naccum++; }
180 }
181 i++;
182 }
183
184 /* also compute midpoint for fork */
185 accum = 0.0; naccum = 0;
186 for( j = nrows; j >= 0; j-- ) {
187 /* if( nest[j] < nest[nrows] ) break; */
188 if( nest[j] == nest[nrows] && func[j] == '+' ) break;
189 else if( nest[j] == nest[nrows] && func[j] == 'T' ) { accum += yval[j]; naccum++; }
190 else if( nest[j] == nest[nrows]+1 && func[j] == '-' ) { accum += yval[j]; naccum++; }
191 }
192 if( naccum > 0 ) yval[nrows] = accum / (double) naccum;
193
194 nrows++;
195 continue;
196 }
197
198 else if( strcmp( toklist[i], "," )==0 ) continue;
199
200 else if( toklist[i][0] != ':' ) {
201 nest[nrows] = nestlevel;
202 label[nrows] = toklist[i];
203 func[nrows] = 'T';
204 xval[nrows] = 0.0;
205 yval[nrows] = 0.0;
206 tagcount++;
207 yval[nrows] = tagcount; /* rooted */
208 /* now look ahead to next tok.. there could be a distance.. */
209 if( toklist[i+1][0] == ':' ) {
210 xval[nrows] = atof( &toklist[i+1][1] );
211 i++;
212 }
213 nrows++;
214 }
215
216 else if( toklist[i][0] == ':' ) { /* no label, only a distance.. */
217 nest[nrows] = nestlevel;
218 label[nrows] = "<nolabel>";
219 func[nrows] = 'D';
220 xval[nrows] = atof( &toklist[i+1][1] );
221 tagcount++;
222 yval[nrows] = tagcount;
223 nrows++;
224 }
225 }
226
227 if( nestlevel > 0 ) return( Eerr( 676, "too many open parens in newickfile, no tree produced", "" ));
228
229
230
231 if( debug ) {
232 fprintf( PLS.diagfp, "\n\n" );
233 for( i = 0; i < nrows; i++ ) fprintf( PLS.diagfp, "%2d. %c %d %s %g %g\n", i, func[i], nest[i], label[i], xval[i], yval[i] );
234 }
235
236
237
238
239
240
241
242 /* now draw the tree.. */
243 yofs = Ecurtextheight * 0.3;
244
245 linedet( "linedetails", linedetails, 1.0 );
246
247
248 /* lines connecting leafs to forks.. */
249 for( i = 0; i < nrows; i++ ) {
250 if( func[i] == 'T' ) {
251 for( j = i; j < nrows; j++ ) {
252 if( func[j] == '-' && nest[j] == nest[i] ) {
253 y1 = Eay( EDYhi - yval[i] ) + yofs; /* flip & adjust upward a bit */
254 y2 = Eay( EDYhi - yval[j] );
255 Emov( Eax( xval[i] ), y1 );
256 if( squaredoff ) Elin( Eax( xval[j] ), y1 );
257 Elin( Eax( xval[j] ), y2 );
258 break;
259 }
260 }
261 }
262 }
263
264 /* lines connecting forks to parent forks.. */
265 for( i = 0; i < nrows; i++ ) {
266 if( func[i] == '-' && nest[i] > 1 ) {
267 for( j = i; j < nrows; j++ ) {
268 if( func[j] == '-' && nest[j] == nest[i]-1 ) {
269 y1 = Eay( EDYhi - yval[i] ); /* flip */
270 y2 = Eay( EDYhi - yval[j] );
271 Emov( Eax( xval[i] ), y1 );
272 if( squaredoff ) Elin( Eax( xval[j] ), y1 );
273 Elin( Eax( xval[j] ), y2 );
274 break;
275 }
276 }
277 }
278 }
279
280 /* dots at forks.. */
281 radius = 0.05;
282 if( strcmp( symbol, "" ) != 0 ) {
283 symdet( "symbol", symbol, symcode, &radius );
284 for( i = 0; i < nrows; i++ ) {
285 if( func[i] == '-' ) {
286 x = Eax( xval[i] );
287 y = Eay( EDYhi - yval[i] ); /* flip */
288 Emark( x, y, symcode, radius );
289 }
290 }
291 }
292
293 /* leaf labels.. */
294 textdet( "textdetails", textdetails, &align, &adjx, &adjy, -1, "R", 1.0 );
295 for( i = 0; i < nrows; i++ ) {
296 if( label[i][0] != '\0' ) {
297 x = Eax( xval[i] );
298 y = Eay( EDYhi - yval[i] ); /* flip */
299 Emov( x+0.01, y );
300 Etext( label[i] );
301 }
302 }
303
304 /* do root dot */
305 if( GL_smember( rootsym, "no none" ) == 0 && squaredoff ) {
306 i = nrows-1;
307 symdet( "root", rootsym, symcode, &radius );
308 /* Emark( Eax( xval[i] ), Eay( EDYhi-yval[i] ), symcode, radius ); */
309 Emark( Eax( xval[i] ), Eay( EDYhi-((double)tagcount/2.0) ), symcode, radius );
310 }
311
312 TDH_setvar( "TREEDONE", "1" );
313 return( 0 );
314 }
315
316
317 #ifdef HOLD
318 /* ========================== */
319 /* ANGPOS - get new position based on current position & distance and angle in rads */
320 static int
angpos(oldx,oldy,dist,ang,newx,newy)321 angpos( oldx, oldy, dist, ang, newx, newy )
322 double oldx, oldy, dist, ang, *newx, *newy;
323 {
324 ang = (ang-HALFPI) * -1.0;
325 *newx = oldx + (dist * cos( ang ));
326 *newy = oldy + (dist * sin( ang ));
327 return( 0 );
328 }
329
330 #endif
331