source: tags/ms_r16q2/TOOLS/arb_read_tree.cxx

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