source: tags/arb-6.0/TOOLS/arb_read_tree.cxx

Last change on this file was 11488, checked in by westram, 8 years ago
  • reintegrates 'tree' into 'trunk'
    • implements #417 (multifurcate tree)
    • tree display
      • adds MULTIFURC MODE
      • reordered modes (synchronizes NTREE and PARSIMONY)
    • branch analysis
      • display number of multifurcations in 'mark long branches'
      • display "in-tree-distance" and "per-species-distance"
    • added function to toggle '100%' bootstraps
    • document bug in GBT_remove_leafs (#452)
  • adds:
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.9 KB
Line 
1// ============================================================= //
2//                                                               //
3//   File      : arb_read_tree.cxx                               //
4//   Purpose   :                                                 //
5//                                                               //
6//   Institute of Microbiology (Technical University Munich)     //
7//   http://www.arb-home.de/                                     //
8//                                                               //
9// ============================================================= //
10
11#include <TreeRead.h>
12#include <arb_strbuf.h>
13#include <arb_defs.h>
14#include <ctime>
15
16
17static void add_bootstrap(GBT_TREE *node, double hundred) {
18    // add_bootstrap interprets the length of the branches as bootstrap value
19    // (this is needed by Phylip DNAPARS/PROTPARS with bootstrapping)
20    //
21    // 'hundred' specifies which value represents 100%
22
23    if (node->is_leaf) {
24        node->remove_remark();
25        return;
26    }
27
28    node->leftlen  /= hundred;
29    node->rightlen /= hundred;
30
31    double left_bs  = node->leftlen  * 100.0;
32    double right_bs = node->rightlen * 100.0;
33
34#if defined(DEBUG) && 0
35    fprintf(stderr, "node->leftlen  = %f left_bs  = %f\n", node->leftlen, left_bs);
36    fprintf(stderr, "node->rightlen = %f right_bs = %f\n", node->rightlen, right_bs);
37#endif // DEBUG
38
39    node->leftson ->set_bootstrap(left_bs);
40    node->rightson->set_bootstrap(right_bs);
41
42    node->leftlen  = DEFAULT_BRANCH_LENGTH; // reset branchlengths
43    node->rightlen = DEFAULT_BRANCH_LENGTH;
44
45    add_bootstrap(node->leftson, hundred);
46    add_bootstrap(node->rightson, hundred);
47}
48
49static void show_message(GBDATA *gb_main, const char *msg) {
50    if (gb_main) {
51        GBT_message(gb_main, msg);
52    }
53    else {
54        fflush(stdout);
55        printf("arb_read_tree: %s\n", msg);
56    }
57}
58static void show_error(GBDATA *gb_main, GB_ERROR error) {
59    if (error) show_message(gb_main, GBS_global_string("Error running arb_read_tree (%s)", error));
60}
61
62static void error_with_usage(GBDATA *gb_main, GB_ERROR error) {
63    fputs("Usage: arb_read_tree [options] tree_name treefile [comment]\n"
64          "Available options:\n"
65          "    -db database savename    specify database and savename (default is 'running ARB')\n"
66          "    -scale factor            scale branchlengths by 'factor'\n"
67          "    -consense numberOfTrees  reinterpret branchlengths as consense values\n"
68          "    -commentFromFile file    read tree comment from 'file'\n"
69          , stdout);
70
71    show_error(gb_main, error);
72}
73
74struct parameters {
75    const char *dbname;
76    const char *dbsavename;
77    const char *tree_name;
78    const char *treefilename;
79    const char *comment;
80    const char *commentFile;
81   
82    bool   scale;
83    double scale_factor;
84
85    bool consense;
86    int  calculated_trees;
87
88    parameters()
89        : dbname(":"),
90          dbsavename(NULL),
91          tree_name(NULL),
92          treefilename(NULL),
93          comment(NULL),
94          commentFile(NULL),
95          scale(false),
96          scale_factor(0.0),
97          consense(false),
98          calculated_trees(0)
99    {
100    }
101   
102#define SHIFT_ARGS(off) do { argc -= off; argv += off; } while (0)
103#define SHIFT_NONSWITCHES(off) do { nonSwitches -= off; nonSwitch += off; } while (0)
104
105    GB_ERROR scan(int argc, char **argv) {
106        GB_ERROR error = NULL;
107
108        const char  *nonSwitch_buf[20];
109        const char **nonSwitch   = nonSwitch_buf;
110        int          nonSwitches = 0;
111
112        SHIFT_ARGS(1);              // position onto first argument
113
114        while (argc>0 && !error) {
115            if (strcmp("-scale", argv[0]) == 0) {
116                scale = true;
117                if (argc<2) error = "-scale expects a 2nd argument (scale factor)";
118                else {
119                    scale_factor = atof(argv[1]);
120                    SHIFT_ARGS(2);
121                }
122            }
123            else if (strcmp("-consense", argv[0]) == 0) {
124                consense = true;
125                if (argc<2) error = "-consense expects a 2nd argument (number of trees)";
126                else {
127                    calculated_trees = atoi(argv[1]);
128                    if (calculated_trees < 1) {
129                        error = GBS_global_string("Illegal # of trees (%i) for -consense (minimum is 1)", calculated_trees);
130                    }
131                    else SHIFT_ARGS(2);
132                }
133            }
134            else if (strcmp("-commentFromFile", argv[0]) == 0) {
135                if (argc<2) error = "-commentFromFile expects a 2nd argument (file containing comment)";
136                else {
137                    commentFile = argv[1];
138                    SHIFT_ARGS(2);
139                }
140            }
141            else if (strcmp("-db", argv[0]) == 0) {
142                if (argc<3) error = "-db expects two arguments (database and savename)";
143                else {
144                    dbname     = argv[1];
145                    dbsavename = argv[2];
146                    SHIFT_ARGS(3);
147                }
148            }
149            else {
150                nonSwitch[nonSwitches++] = argv[0];
151                SHIFT_ARGS(1);
152            }
153        }
154
155        if (!error) {
156            if (!nonSwitches) error = "Missing argument 'tree_name'";
157            else {
158                tree_name = nonSwitch[0];
159                SHIFT_NONSWITCHES(1);
160            }
161        }
162        if (!error) {
163            if (!nonSwitches) error = "Missing argument 'treefile'";
164            else {
165                treefilename = nonSwitch[0];
166                SHIFT_NONSWITCHES(1);
167            }
168        }
169        if (!error && nonSwitches>0) {
170            comment = nonSwitch[0];
171            SHIFT_NONSWITCHES(1);
172        }
173
174        if (!error && nonSwitches>0) {
175            error = GBS_global_string("unexpected argument(s): %s ..", nonSwitch[0]);
176        }
177        return error;
178    }
179};
180
181int main(int argc, char **argv) {
182    parameters param;
183    GB_ERROR error = param.scan(argc, argv);
184
185    GBDATA   *gb_main      = NULL;
186    GBDATA   *gb_msg_main  = NULL;
187    bool      connectToArb = strcmp(param.dbname, ":") == 0;
188    GB_shell  shell;
189
190    if (!error || connectToArb) {
191        gb_main                       = GB_open(param.dbname, connectToArb ? "r" : "rw");
192        if (connectToArb) gb_msg_main = gb_main;
193    }
194
195    if (error) error_with_usage(gb_main, error);
196    else {
197        if (!gb_main) {
198            if (connectToArb) error = "you have to start an arbdb server first";
199            else error = GBS_global_string("can't open db (Reason: %s)", GB_await_error());
200        }
201
202        char     *comment_from_file     = 0;
203        char     *comment_from_treefile = 0;
204        GBT_TREE *tree                  = 0;
205
206        if (!error) {
207            if (param.commentFile) {
208                comment_from_file = GB_read_file(param.commentFile);
209                if (!comment_from_file) {
210                    comment_from_file = GBS_global_string_copy("Error reading from comment-file '%s':\n%s", param.commentFile, GB_await_error());
211                }
212            }
213
214            show_message(gb_msg_main, GBS_global_string("Reading tree from '%s' ..", param.treefilename));
215            {
216                char *warnings             = 0;
217                bool  allow_length_scaling = !param.consense && !param.scale;
218
219                tree = TREE_load(param.treefilename, GBT_TREE_NodeFactory(), &comment_from_treefile, allow_length_scaling, &warnings);
220                if (!tree) {
221                    error = GB_await_error();
222                }
223                else if (warnings) {
224                    show_message(gb_msg_main, warnings);
225                    free(warnings);
226                }
227            }
228        }
229
230        if (!error) {
231            if (param.scale) {
232                show_message(gb_msg_main, GBS_global_string("Scaling branch lengths by factor %f", param.scale_factor));
233                TREE_scale(tree, param.scale_factor, 1.0);
234            }
235
236            if (param.consense) {
237                if (param.calculated_trees < 1) {
238                    error = "Minimum for -consense is 1";
239                }
240                else {
241                    show_message(gb_msg_main, GBS_global_string("Reinterpreting branch lengths as consense values (%i trees)", param.calculated_trees));
242                    add_bootstrap(tree, param.calculated_trees);
243                }
244            }
245        }
246
247        if (!error) {
248            error = GB_begin_transaction(gb_main);
249
250            if (!error && tree->is_leaf) error = "Cannot load tree (need at least 2 leafs)";
251            if (!error) error                  = GBT_write_tree(gb_main, param.tree_name, tree);
252
253            if (!error) {
254                // write tree comment
255                const char *comments[] = {
256                    param.comment,
257                    comment_from_file,
258                    comment_from_treefile,
259                };
260
261                GBS_strstruct *buf   = GBS_stropen(5000);
262                bool           empty = true;
263
264                for (size_t c = 0; c<ARRAY_ELEMS(comments); c++) {
265                    if (comments[c]) {
266                        if (!empty) GBS_chrcat(buf, '\n');
267                        GBS_strcat(buf, comments[c]);
268                        empty = false;
269                    }
270                }
271
272                char *cmt = GBS_strclose(buf);
273
274                error = GBT_write_tree_remark(gb_main, param.tree_name, cmt);
275
276                free(cmt);
277            }
278
279            error = GB_end_transaction(gb_main, error);
280        }
281
282        if (error) show_error(gb_main, error);
283        else       show_message(gb_msg_main, GBS_global_string("Tree %s read into the database", param.tree_name));
284
285        delete tree;
286        free(comment_from_file);
287        free(comment_from_treefile);
288    }
289
290    if (gb_main) {
291        if (!error && !connectToArb) {
292            error = GB_save_as(gb_main, param.dbsavename, "a");
293            if (error) show_error(gb_main, error);
294        }
295        GB_close(gb_main);
296    }
297
298    return error ? EXIT_FAILURE : EXIT_SUCCESS;
299}
Note: See TracBrowser for help on using the repository browser.