source: tags/ms_r18q1/TOOLS/arb_read_tree.cxx

Last change on this file was 17110, checked in by westram, 6 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.8 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 GBDATA *gb_msg_main = NULp;
51
52static void show_message(const char *msg) {
53    if (gb_msg_main) {
54        GBT_message(gb_msg_main, msg);
55    }
56    else {
57        fflush(stdout);
58        printf("arb_read_tree: %s\n", msg);
59    }
60}
61static void show_error(GB_ERROR error) {
62    if (error) show_message(GBS_global_string("Error running arb_read_tree (%s)", error));
63}
64
65static void error_with_usage(GB_ERROR error) {
66    fputs("Usage: arb_read_tree [options] tree_name treefile [comment]\n"
67          "Available options:\n"
68          "    -db database savename    specify database and savename (default is 'running ARB')\n"
69          "    -scale factor            scale branchlengths by 'factor'\n"
70          "    -consense numberOfTrees  reinterpret branchlengths as consense values\n"
71          "    -commentFromFile file    read tree comment from 'file'\n"
72          , stdout);
73
74    show_error(error);
75}
76
77struct parameters {
78    const char *dbname;
79    const char *dbsavename;
80    const char *tree_name;
81    const char *treefilename;
82    const char *comment;
83    const char *commentFile;
84   
85    bool   scale;
86    double scale_factor;
87
88    bool consense;
89    int  calculated_trees;
90
91    parameters()
92        : dbname(":"),
93          dbsavename(NULp),
94          tree_name(NULp),
95          treefilename(NULp),
96          comment(NULp),
97          commentFile(NULp),
98          scale(false),
99          scale_factor(0.0),
100          consense(false),
101          calculated_trees(0)
102    {}
103   
104#define SHIFT_ARGS(off) do { argc -= off; argv += off; } while (0)
105#define SHIFT_NONSWITCHES(off) do { nonSwitches -= off; nonSwitch += off; } while (0)
106
107    GB_ERROR scan(int argc, char **argv) {
108        GB_ERROR error = NULp;
109
110        const char  *nonSwitch_buf[20];
111        const char **nonSwitch   = nonSwitch_buf;
112        int          nonSwitches = 0;
113
114        SHIFT_ARGS(1);              // position onto first argument
115
116        while (argc>0 && !error) {
117            if (strcmp("-scale", argv[0]) == 0) {
118                scale = true;
119                if (argc<2) error = "-scale expects a 2nd argument (scale factor)";
120                else {
121                    scale_factor = atof(argv[1]);
122                    SHIFT_ARGS(2);
123                }
124            }
125            else if (strcmp("-consense", argv[0]) == 0) {
126                consense = true;
127                if (argc<2) error = "-consense expects a 2nd argument (number of trees)";
128                else {
129                    calculated_trees = atoi(argv[1]);
130                    if (calculated_trees < 1) {
131                        error = GBS_global_string("Illegal # of trees (%i) for -consense (minimum is 1)", calculated_trees);
132                    }
133                    else SHIFT_ARGS(2);
134                }
135            }
136            else if (strcmp("-commentFromFile", argv[0]) == 0) {
137                if (argc<2) error = "-commentFromFile expects a 2nd argument (file containing comment)";
138                else {
139                    commentFile = argv[1];
140                    SHIFT_ARGS(2);
141                }
142            }
143            else if (strcmp("-db", argv[0]) == 0) {
144                if (argc<3) error = "-db expects two arguments (database and savename)";
145                else {
146                    dbname     = argv[1];
147                    dbsavename = argv[2];
148                    SHIFT_ARGS(3);
149                }
150            }
151            else {
152                nonSwitch[nonSwitches++] = argv[0];
153                SHIFT_ARGS(1);
154            }
155        }
156
157        if (!error) {
158            if (!nonSwitches) error = "Missing argument 'tree_name'";
159            else {
160                tree_name = nonSwitch[0];
161                SHIFT_NONSWITCHES(1);
162            }
163        }
164        if (!error) {
165            if (!nonSwitches) error = "Missing argument 'treefile'";
166            else {
167                treefilename = nonSwitch[0];
168                SHIFT_NONSWITCHES(1);
169            }
170        }
171        if (!error && nonSwitches>0) {
172            comment = nonSwitch[0];
173            SHIFT_NONSWITCHES(1);
174        }
175
176        if (!error && nonSwitches>0) {
177            error = GBS_global_string("unexpected argument(s): %s ..", nonSwitch[0]);
178        }
179        return error;
180    }
181};
182
183int main(int argc, char **argv) {
184    parameters param;
185    GB_ERROR error = param.scan(argc, argv);
186
187    GBDATA   *gb_main      = NULp;
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(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     = NULp;
204        char     *comment_from_treefile = NULp;
205        TreeNode *tree                  = NULp;
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(GBS_global_string("Reading tree from '%s' ..", param.treefilename));
216            {
217                char *warnings             = NULp;
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(warnings);
226                    free(warnings);
227                }
228            }
229        }
230
231        if (!error) {
232            if (param.scale) {
233                show_message(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(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 tree_remark(5000);
263                bool          empty = true;
264
265                for (size_t c = 0; c<ARRAY_ELEMS(comments); c++) {
266                    if (comments[c]) {
267                        if (!empty) tree_remark.put('\n');
268                        tree_remark.cat(comments[c]);
269                        empty = false;
270                    }
271                }
272
273                error = GBT_write_tree_remark(gb_main, param.tree_name, tree_remark.get_data());
274            }
275
276            error = GB_end_transaction(gb_main, error);
277        }
278
279        if (error) show_error(error);
280        else       show_message(GBS_global_string("Tree %s read into the database", param.tree_name));
281
282        UNCOVERED();
283        destroy(tree);
284        free(comment_from_file);
285        free(comment_from_treefile);
286    }
287
288    if (gb_main) {
289        if (!error && !connectToArb) {
290            error = GB_save_as(gb_main, param.dbsavename, "a");
291            if (error) show_error(error);
292        }
293        GB_close(gb_main);
294    }
295
296    return error ? EXIT_FAILURE : EXIT_SUCCESS;
297}
Note: See TracBrowser for help on using the repository browser.