source: branches/profile/DIST/DI_foundclusters.cxx

Last change on this file was 12877, checked in by westram, 5 years ago
File size: 15.5 KB
Line 
1// ================================================================ //
2//                                                                  //
3//   File      : DI_foundclusters.cxx                               //
4//   Purpose   : Store results of cluster detection                 //
5//                                                                  //
6//   Coded by Ralf Westram (coder@reallysoft.de) in November 2009   //
7//   Institute of Microbiology (Technical University Munich)        //
8//   http://www.arb-home.de/                                        //
9//                                                                  //
10// ================================================================ //
11
12#include "di_foundclusters.hxx"
13#include "di_clustertree.hxx"
14
15#include <aw_window.hxx>
16#include <aw_select.hxx>
17
18#include <cmath>
19
20using namespace std;
21
22ID Cluster::unused_id = 1;
23
24static AP_FLOAT calc_mean_dist_to(const ClusterTree *distance_to, const LeafRelations *distances) {
25    size_t   count = 0;
26    AP_FLOAT sum   = 0.0;
27
28    LeafRelationCIter dist_end = distances->end();
29    for (LeafRelationCIter dist = distances->begin(); dist != dist_end; ++dist) {
30        const TwoLeafs& leafs = dist->first;
31
32        if (leafs.first() == distance_to || leafs.second() == distance_to) {
33            count++;
34            sum += dist->second;
35        }
36    }
37
38    cl_assert(count);                               // 'distance_to' does not occur in 'distances'
39
40    return sum/count;
41}
42
43Cluster::Cluster(ClusterTree *ct)
44    : next_desc(NULL)
45{
46    cl_assert(ct->get_state() == CS_IS_CLUSTER);
47
48    id = unused_id++;
49
50    const LeafRelations *distances = ct->get_sequence_dists();
51    cl_assert(distances);
52
53    LeafRelationCIter dist     = distances->begin();
54    LeafRelationCIter dist_end = distances->end();
55
56    cl_assert(dist != dist_end); // oops - empty cluster
57
58    min_dist  = max_dist = dist->second;
59    mean_dist = 0.0;
60
61    AP_FLOAT min_mean_dist = 9999999999.0;
62    representative         = NULL;
63
64    for (; dist != dist_end; ++dist) {
65        const TwoLeafs& leafs = dist->first;
66        AP_FLOAT        di    = dist->second;
67
68        if (di<min_dist) min_dist = di;
69        if (di>max_dist) max_dist = di;
70
71        mean_dist += di;
72
73        for (int i = 0; i<2; ++i) {
74            const ClusterTree *member     = i ? leafs.first() : leafs.second();
75            GBDATA            *gb_species = member->gb_node;
76
77            cl_assert(gb_species);
78            if (members.find(gb_species) == members.end()) { // not in set
79                members.insert(gb_species);
80                AP_FLOAT mean_dist_to_member = calc_mean_dist_to(member, distances);
81
82                if (mean_dist_to_member<min_mean_dist) {
83                    representative = gb_species;
84                    min_mean_dist  = mean_dist_to_member;
85                }
86            }
87        }
88
89    }
90
91    cl_assert(members.size() == ct->get_leaf_count());
92
93    min_bases = ct->get_min_bases();
94    {
95        ClusterTree *root    = ct->get_root_node();
96        size_t       allSize = root->get_leaf_count();
97        size_t       size    = ct->get_leaf_count();
98        size_t       first   = ct->relative_position_in(root);
99        size_t       last    = first+size-1;
100
101        rel_tree_pos = ((double(first)+last)/2) / allSize;
102    }
103    mean_dist /= distances->size();
104    cl_assert(representative);
105    desc = create_description(ct);
106}
107
108static string *get_downgroups(const ARB_countedTree *ct, const ARB_tree_predicate& keep_group_name, size_t& group_members) {
109    string *result = NULL;
110    if (ct->is_leaf) {
111        group_members = 0;
112    }
113    else {
114        const char *group_name = ct->get_group_name();
115
116        if (group_name && keep_group_name.selects(*ct)) {
117            group_members = ct->get_leaf_count();
118            result        = new string(ct->name);
119        }
120        else {
121            string *leftgroups  = get_downgroups(ct->get_leftson(), keep_group_name, group_members);
122            size_t  right_members;
123            string *rightgroups = get_downgroups(ct->get_rightson(), keep_group_name, right_members);
124
125            group_members += right_members;
126
127            if (leftgroups || rightgroups) {
128                if (!leftgroups)       { result = rightgroups; rightgroups = NULL; }
129                else if (!rightgroups) { result = leftgroups; leftgroups = NULL; }
130                else                     result = new string(*leftgroups+"+"+*rightgroups);
131
132                delete leftgroups;
133                delete rightgroups;
134            }
135        }
136    }
137    return result;
138}
139static const char *get_upgroup(const ARB_countedTree *ct, const ARB_tree_predicate& keep_group_name, const ARB_countedTree*& upgroupPtr, bool reuse_original_name) {
140    const char *group_name = ct->get_group_name();
141    char       *original   = (reuse_original_name && group_name) ? originalGroupName(group_name) : NULL;
142
143    if (original) {
144        static SmartCharPtr result;
145        result = original;
146        group_name = &*result;
147        upgroupPtr = ct;
148    }
149    else if (group_name && keep_group_name.selects(*ct)) {
150        upgroupPtr = ct;
151    }
152    else {
153        const ARB_countedTree *father = ct->get_father();
154        if (father) {
155            group_name = get_upgroup(father, keep_group_name, upgroupPtr, reuse_original_name);
156        }
157        else {
158            upgroupPtr = ct;
159            group_name = "WHOLE_TREE";
160        }
161    }
162
163    cl_assert(!upgroupPtr || ct->is_inside(upgroupPtr));
164    cl_assert(group_name);
165    return group_name;
166}
167
168string Cluster::create_description(const ARB_countedTree *ct) {
169    // creates cluster description (using up- and down-groups, percentages, existing groups, ...)
170    cl_assert(!ct->is_leaf);
171
172    string name;
173
174    UseAnyTree use_any_upgroup;
175
176    const ARB_countedTree *upgroup         = NULL;
177    const char            *upgroup_name    = get_upgroup(ct, use_any_upgroup, upgroup, false);
178    size_t                 upgroup_members = upgroup->get_leaf_count();
179
180    size_t cluster_members = ct->get_leaf_count();
181
182    cl_assert(upgroup_members>0); // if no group, whole tree should have been returned
183
184    if (upgroup_members == cluster_members) { // use original or existing name
185        name = string("[G] ")+upgroup_name;
186    }
187    else {
188        size_t  downgroup_members;
189        string *downgroup_names = get_downgroups(ct, use_any_upgroup, downgroup_members);
190
191        AP_FLOAT up_rel   = (AP_FLOAT)cluster_members/upgroup_members;   // used for: xx% of upgroup (big is good)
192        AP_FLOAT down_rel = (AP_FLOAT)downgroup_members/cluster_members; // used for: downgroup(s) + (1-xx)% outgroup (big is good)
193
194        if (up_rel>down_rel) { // prefer up_percent
195            name = string(GBS_global_string("%4.1f%% of %s", up_rel*100.0, upgroup_name));
196        }
197        else {
198            if (downgroup_members == cluster_members) {
199                name = downgroup_names->c_str();
200            }
201            else {
202                name = string(GBS_global_string("%s + %4.1f%% outgroup", downgroup_names->c_str(), (1-down_rel)*100.0));
203            }
204        }
205
206        delete downgroup_names;
207    }
208    return name;
209}
210
211string Cluster::get_upgroup_info(const ARB_countedTree *ct, const ARB_tree_predicate& keep_group_name) {
212    // generates a string indicating relative position in up-group
213    cl_assert(!ct->is_leaf);
214
215    const ARB_countedTree *upgroup      = NULL;
216    const char            *upgroup_name = get_upgroup(ct, keep_group_name, upgroup, true);
217
218    size_t upgroup_members = upgroup->get_leaf_count();
219    size_t cluster_members = ct->get_leaf_count();
220
221    const char *upgroup_info;
222    if (upgroup_members == cluster_members) {
223        upgroup_info = upgroup_name;
224    }
225    else {
226        size_t cluster_pos1 = ct->relative_position_in(upgroup)+1; // range [1..upgroup_members]
227        size_t cluster_posN = cluster_pos1+cluster_members-1;
228
229        upgroup_info = GBS_global_string("%zu-%zu/%zu_%s", cluster_pos1, cluster_posN, upgroup_members, upgroup_name);
230    }
231    return upgroup_info;
232}
233
234// --------------------------
235//      DisplayFormat
236
237class DisplayFormat : virtual Noncopyable {
238    size_t   max_count;
239    AP_FLOAT max_dist;
240    AP_FLOAT max_minBases;
241
242    mutable char *format_string;
243
244    static int calc_digits(long val) {
245        int digits;
246        if (val>0) digits      = log10(val)+1;
247        else if (val<0) digits = calc_digits(-val);
248        else digits            = 1;
249
250        cl_assert(digits>0);
251       
252        return digits;
253    }
254
255    static char *make_format(size_t val) {
256        int digits = calc_digits(val);
257        return GBS_global_string_copy("%%%izu", digits);
258    }
259
260    static void calc_float_format(AP_FLOAT val, int& all, int& afterdot) {
261        int digits   = calc_digits(long(val));
262        afterdot = digits <=3 ? 5-digits-1 : 0;
263        cl_assert(afterdot >= 0);
264        all = digits+(afterdot ? (afterdot+1) : 0);
265    }
266
267    static char *make_format(AP_FLOAT val) {
268        int all, afterdot;
269        calc_float_format(val, all, afterdot);
270        return GBS_global_string_copy("%%%i.%if", all, afterdot);
271    }
272
273public:
274
275    DisplayFormat()
276        : max_count(0),
277          max_dist(0.0),
278          max_minBases(0.0),
279          format_string(NULL)
280    {}
281
282    ~DisplayFormat() {
283        free(format_string);
284    }
285
286    void announce(size_t count, AP_FLOAT maxDist, AP_FLOAT minBases) {
287        max_count    = std::max(max_count, count);
288        max_dist     = std::max(max_dist, maxDist);
289        max_minBases = std::max(max_minBases, minBases);
290
291        freenull(format_string);
292    }
293
294    const char *get_header() const {
295        int countSize = calc_digits(max_count) + 2; // allow to use two of the spaces behind count
296        int distSize, baseSize;
297        {
298            int afterdot;
299            calc_float_format(max_dist,     distSize, afterdot);
300            calc_float_format(max_minBases, baseSize, afterdot);
301        }
302
303        return GBS_global_string("%-*s%*s  %*s  Groupname",
304                                 countSize, countSize<5 ? "#" : "Count",
305                                 distSize, "Mean [min-max] dist",
306                                 baseSize, "Bases");
307    }
308
309    const char *get_format() const {
310        if (!format_string) {
311            // create format string for description
312            char *count_part = make_format(max_count);
313            char *dist_part  = make_format(max_dist);
314            char *bases_part = make_format(max_minBases);
315
316            format_string = GBS_global_string_copy("%s  %s [%s-%s]  %s  %%s",
317                                                   count_part,
318                                                   dist_part, dist_part, dist_part,
319                                                   bases_part);
320
321            free(bases_part);
322            free(dist_part);
323            free(count_part);
324        }
325        return format_string;
326    }
327};
328
329
330
331void Cluster::scan_display_widths(DisplayFormat& format) const {
332    AP_FLOAT max_of_all_dists = std::max(max_dist, std::max(min_dist, mean_dist));
333    format.announce(get_member_count(), max_of_all_dists*100.0, min_bases);
334}
335
336const char *Cluster::get_list_display(const DisplayFormat *format) const {
337    const char *format_string;
338    if (format) format_string = format->get_format();
339    else        format_string = "%zu  %.3f [%.3f-%.3f]  %.1f  %s";
340
341    return GBS_global_string(format_string,
342                             get_member_count(),
343                             mean_dist*100.0,
344                             min_dist*100.0,
345                             max_dist*100.0,
346                             min_bases,
347                             desc.c_str());
348}
349
350
351// ---------------------
352//      ClustersData
353
354void ClustersData::add(ClusterPtr clus, ClusterSubset subset) {
355    known_clusters[clus->get_ID()] = clus;
356    get_subset(subset).push_back(clus->get_ID());
357}
358
359void ClustersData::remove(ClusterPtr clus, ClusterSubset subset) {
360    int pos = get_pos(clus, subset);
361    if (pos != -1) {
362        ClusterIDs&           ids      = get_subset(subset);
363        ClusterIDs::iterator  toRemove = ids.begin();
364        advance(toRemove, pos);
365        ids.erase(toRemove);
366        known_clusters.erase(clus->get_ID());
367    }
368}
369
370void ClustersData::clear(ClusterSubset subset) {
371    ClusterIDs&     ids    = get_subset(subset);
372    ClusterIDsIter  id_end = ids.end();
373    for (ClusterIDsIter  id = ids.begin(); id != id_end; ++id) {
374        known_clusters.erase(*id);
375    }
376    ids.clear();
377}
378
379void ClustersData::store(ID id) {
380    ClusterPtr toStore = clusterWithID(id);
381    remove(toStore, SHOWN_CLUSTERS);
382    add(toStore, STORED_CLUSTERS);
383}
384
385void ClustersData::store_all() {
386    copy(shown.begin(), shown.end(), back_insert_iterator<ClusterIDs>(stored));
387    shown.clear();
388    sort_needed = false;
389}
390
391void ClustersData::restore_all() {
392    copy(stored.begin(), stored.end(), back_insert_iterator<ClusterIDs>(shown));
393    stored.clear();
394    sort_needed = true;
395}
396
397void ClustersData::swap_all() {
398    swap(stored, shown);
399    sort_needed = true;
400}
401
402
403class SortClusterIDsBy {
404    const KnownClusters& known;
405    ClusterOrder         primary_order;
406    ClusterOrder         secondary_order;
407public:
408    SortClusterIDsBy(const KnownClusters& known_, ClusterOrder primary, ClusterOrder secondary)
409        : known(known_)
410        , primary_order(primary)
411        , secondary_order(secondary)
412    {}
413
414     bool operator()(ID id1, ID id2) const {
415         ClusterPtr cl1 = known.find(id1)->second;
416         ClusterPtr cl2 = known.find(id2)->second;
417
418         cl_assert(!cl1.isNull() && !cl2.isNull());
419
420         bool less = cl1->lessByOrder(*cl2, primary_order);
421         if (!less) {
422             bool equal = !cl2->lessByOrder(*cl1, primary_order);
423             if (equal) {
424                 less = cl1->lessByOrder(*cl2, secondary_order);
425             }
426         }
427         return less;
428     }
429};
430
431void ClustersData::update_cluster_selection_list() {
432    cl_assert(clusterList);
433    clusterList->clear();
434
435    if (shown.empty()) {
436        clusterList->insert_default("<No clusters detected>", ID(0));
437    }
438    else {
439        {
440            SortClusterIDsBy sortBy(known_clusters, criteria[0], criteria[1]);
441            sort(shown.begin(), shown.end(), sortBy);
442        }
443        {
444            ClusterIDsIter    cl_end = shown.end();
445            DisplayFormat format;
446
447            for (int pass = 1; pass <= 2; ++pass) {
448                if (pass == 2) {
449                    clusterList->insert(format.get_header(), ID(-1));
450                }
451                for (ClusterIDsIter cl = shown.begin(); cl != cl_end; ++cl) {
452                    ID         id      = *cl;
453                    ClusterPtr cluster = clusterWithID(id);
454
455                    if (pass == 1) {
456                        cluster->scan_display_widths(format);
457                    }
458                    else {
459                        clusterList->insert(cluster->get_list_display(&format), cluster->get_ID());
460                    }
461                }
462            }
463            clusterList->insert_default("<select no cluster>", ID(0));
464        }
465    }
466    clusterList->update();
467}
468
469void ClustersData::free() {
470    // delete calculated/stored data
471    shown.clear();
472    stored.clear();
473    update_cluster_selection_list();
474    known_clusters.clear();
475}
476
477char *originalGroupName(const char *groupname) {
478    char *original = NULL;
479    const char *was = strstr(groupname, " {was:");
480    const char *closing = strchr(groupname, '}');
481    if (was && closing) {
482        original = GB_strpartdup(was+6, closing-1);
483        if (original[0]) {
484            char *sub = originalGroupName(original);
485            if (sub) freeset(original, sub);
486        }
487        else {
488            freenull(original); // empty -> invalid
489        }
490    }
491    return original;
492}
493
494
495
Note: See TracBrowser for help on using the repository browser.