source: branches/stable/DIST/DI_clustertree.cxx

Last change on this file was 11455, checked in by westram, 6 years ago
  • use is_nan etc
File size: 13.8 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : DI_clustertree.cxx                                //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Coded by Ralf Westram (coder@reallysoft.de) in October 2009   //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "di_clustertree.hxx"
13#include <arb_progress.h>
14#include <set>
15#include <cmath>
16#include <limits>
17
18using namespace std;
19
20typedef std::set<LeafRelation> SortedPairValues; // iterator runs from small to big values
21
22// ------------------------
23//      ClusterTreeRoot
24
25
26ClusterTreeRoot::ClusterTreeRoot(AliView *aliview, AP_sequence *seqTemplate_, AP_FLOAT maxDistance_, size_t minClusterSize_)
27    : ARB_seqtree_root(aliview, new ClusterTreeNodeFactory, seqTemplate_, false),
28      maxDistance(maxDistance_),
29      minClusterSize(minClusterSize_)
30{}
31
32#if defined(DEBUG)
33struct ClusterStats {
34    size_t clusters;
35    size_t subclusters;
36    size_t loadedSequences;
37
38    ClusterStats()
39        : clusters(0)
40        , subclusters(0)
41        , loadedSequences(0)
42    {}
43};
44
45static void collectStats(ClusterTree *node, ClusterStats *stats) {
46    switch (node->get_state()) {
47        case CS_IS_CLUSTER: stats->clusters++; break;
48        case CS_SUB_CLUSTER: stats->subclusters++; break;
49        default: break;
50    }
51    if (node->is_leaf) {
52        AP_sequence *seq = node->get_seq();
53        if (seq->got_sequence()) stats->loadedSequences++;
54    }
55    else {
56        collectStats(node->get_leftson(), stats);
57        collectStats(node->get_rightson(), stats);
58    }
59}
60#endif // DEBUG
61
62GB_ERROR ClusterTreeRoot::find_clusters() {
63    ClusterTree *root = get_root_node();
64
65    root->init_tree();
66#if defined(DEBUG)
67    printf("----------------------------------------\n");
68    printf("maxDistance:       %f\n", maxDistance);
69    printf("minClusterSize:    %u\n", minClusterSize);
70    printf("Leafs in tree:     %u\n", root->get_leaf_count());
71    printf("Possible clusters: %u\n", root->get_cluster_count());
72#endif // DEBUG
73
74    arb_progress cluster_progress(root->get_cluster_count());
75    cluster_progress.auto_subtitles("Cluster");
76
77    GB_ERROR error = NULL;
78
79    try {
80        root->detect_clusters(cluster_progress);
81    }
82    catch (const char *msg) { error = msg; }
83    catch (...) { cl_assert(0); }
84
85    if (error) {
86        // avoid further access after error
87        change_root(root, NULL);
88        delete root;
89    }
90    else {
91#if defined(DEBUG)
92        ClusterStats stats;
93        collectStats(root, &stats);
94
95        printf("Found clusters:::::::: %zu\n", stats.clusters);
96        printf("Found subclusters::::: %zu\n", stats.subclusters);
97        printf("Loaded sequences:::::: %zu (%5.2f%%)\n",
98               stats.loadedSequences,
99               (100.0*stats.loadedSequences)/root->get_leaf_count());
100
101#if defined(TRACE_DIST_CALC)
102        size_t distances         = root->get_calculated_distances();
103        size_t leafs             = root->get_leaf_count();
104        size_t existingDistances = (leafs*leafs)/2-1;
105
106        printf("Calculated distances:: %zu (%5.2f%%)\n",
107               distances,
108               (100.0*distances)/existingDistances);
109#endif                                              // TRACE_DIST_CALC
110
111#endif                                              // DEBUG
112    }
113
114    return error;
115}
116
117// --------------------
118//      ClusterTree
119
120ClusterTree *ClusterTree::get_cluster(size_t num) {
121    cl_assert(num < get_cluster_count());
122
123    if (num == 0) return this;
124
125    ClusterTree *lson = get_leftson();
126    if (lson->get_cluster_count() <= num) {
127        return get_rightson()->get_cluster(num - lson->get_cluster_count());
128    }
129    return lson->get_cluster(num);
130}
131
132void ClusterTree::init_tree() {
133    cl_assert(state == CS_UNKNOWN);
134
135#if defined(TRACE_DIST_CALC)
136    calculatedDistances = 0;
137#endif // TRACE_DIST_CALC
138
139    if (get_father()) {
140        depth = get_father()->get_depth()+1;
141    }
142    else {
143        depth = 1;
144        cl_assert(is_root_node());
145    }
146
147    if (is_leaf) {
148        leaf_count = 1;
149        clus_count = 0;
150        state      = CS_TOO_SMALL;
151    }
152    else {
153        ClusterTree *lson = get_leftson();
154        ClusterTree *rson = get_rightson();
155
156        lson->init_tree();
157        rson->init_tree();
158
159        leaf_count = lson->get_leaf_count() + rson->get_leaf_count();
160
161        ClusterTreeRoot *troot = get_tree_root();
162        if (leaf_count<troot->get_minClusterSize()) {
163            state      = CS_TOO_SMALL;
164            clus_count = 0;
165        }
166        else {
167            state      = CS_MAYBE_CLUSTER;
168            clus_count = lson->get_cluster_count() + rson->get_cluster_count() + 1;
169            min_bases  = numeric_limits<AP_FLOAT>::infinity();
170        }
171    }
172   
173    cl_assert(state != CS_UNKNOWN);
174}
175
176void ClusterTree::detect_clusters(arb_progress& progress) {
177    if (state == CS_MAYBE_CLUSTER) {
178        cl_assert(!is_leaf);
179
180        ClusterTree *lson = get_leftson();
181        ClusterTree *rson = get_rightson();
182
183        lson->detect_clusters(progress);
184        rson->detect_clusters(progress);
185
186#if defined(TRACE_DIST_CALC)
187        calculatedDistances += get_leftson()->calculatedDistances;
188        calculatedDistances += get_rightson()->calculatedDistances;
189#endif // TRACE_DIST_CALC
190
191        if (lson->get_state() == CS_NO_CLUSTER || rson->get_state() == CS_NO_CLUSTER) {
192            // one son is no cluster -> can't be cluster myself
193            state = CS_NO_CLUSTER;
194        }
195        else {
196            state = CS_IS_CLUSTER;                  // assume this is a cluster
197
198            // Get set of branches (sorted by branch distance)
199
200            get_branch_dists();                     // calculates branchDists
201            SortedPairValues pairsByBranchDistance; // biggest branch distances at end
202            {
203                LeafRelationCIter bd_end = branchDists->end();
204                for (LeafRelationCIter bd = branchDists->begin(); bd != bd_end; ++bd) {
205                    pairsByBranchDistance.insert(LeafRelation(bd->first, bd->second));
206                }
207            }
208
209            cl_assert(pairsByBranchDistance.size() == possible_relations());
210
211            // calculate real sequence distance
212            // * stop when distance > maxDistance is detected
213            // * calculate longest branchdistance first
214            //   (Assumption is: big branch distance <=> big sequence distance)
215
216            AP_FLOAT maxDistance       = get_tree_root()->get_maxDistance();
217            AP_FLOAT worstDistanceSeen = -1.0;
218
219            cl_assert(!sequenceDists);
220            sequenceDists = new LeafRelations;
221
222            SortedPairValues::const_reverse_iterator pair_end = pairsByBranchDistance.rend();
223            for (SortedPairValues::const_reverse_iterator pair = pairsByBranchDistance.rbegin(); pair != pair_end; ++pair) {
224                AP_FLOAT realDist = get_seqDist(pair->get_pair());
225
226                if (realDist>worstDistanceSeen) {
227                    worstDistanceSeen = realDist;
228
229                    delete worstKnownDistance;
230                    worstKnownDistance = new TwoLeafs(pair->get_pair());
231
232                    if (realDist>maxDistance) {     // big distance found -> not a cluster
233                        state = CS_NO_CLUSTER;
234                        break;
235                    }
236                }
237            }
238
239#if defined(TRACE_DIST_CALC)
240            calculatedDistances += sequenceDists->size();
241#endif // TRACE_DIST_CALC
242
243            if (state == CS_IS_CLUSTER) {
244#if defined(DEBUG)
245                // test whether all distances have been calculated
246                cl_assert(!is_leaf);
247                size_t knownDistances = sequenceDists->size() +
248                    get_leftson()->known_seqDists() +
249                    get_rightson()->known_seqDists();
250
251                cl_assert(knownDistances == possible_relations());
252#endif // DEBUG
253
254                get_leftson()->state  = CS_SUB_CLUSTER;
255                get_rightson()->state = CS_SUB_CLUSTER;
256            }
257
258        }
259
260        cl_assert(state == CS_NO_CLUSTER || state == CS_IS_CLUSTER); // detection failure
261
262        lson->oblivion(false);
263        rson->oblivion(false);
264
265        progress.inc();
266        if (progress.aborted()) throw "aborted on userrequest";
267    }
268
269    cl_assert(state != CS_MAYBE_CLUSTER);
270}
271
272void ClusterTree::oblivion(bool forgetDistances) {
273    delete branchDepths;
274    branchDepths = NULL;
275
276    delete branchDists;
277    branchDists = NULL;
278
279    if (state == CS_NO_CLUSTER) {
280        if (forgetDistances) {
281            delete sequenceDists;
282            sequenceDists = NULL;
283        }
284        forgetDistances = true;
285    }
286    else {
287        cl_assert(state & (CS_IS_CLUSTER|CS_SUB_CLUSTER|CS_TOO_SMALL));
288    }
289
290    if (!is_leaf) {
291        get_leftson()->oblivion(forgetDistances);
292        get_rightson()->oblivion(forgetDistances);
293    }
294}
295
296void ClusterTree::calc_branch_depths() {
297    cl_assert(!branchDepths);
298    branchDepths = new NodeValues;
299    if (is_leaf) {
300        (*branchDepths)[this] = 0.0;
301    }
302    else {
303        for (int s = 0; s<2; s++) {
304            const NodeValues *sonDepths = s ? get_leftson()->get_branch_depths() : get_rightson()->get_branch_depths();
305            AP_FLOAT          lenToSon  = s ? leftlen : rightlen;
306
307            NodeValues::const_iterator sd_end = sonDepths->end();
308            for (NodeValues::const_iterator sd = sonDepths->begin(); sd != sd_end; ++sd) {
309                (*branchDepths)[sd->first] = sd->second+lenToSon;
310            }
311        }
312    }
313}
314
315void ClusterTree::calc_branch_dists() {
316    cl_assert(!branchDists);
317    if (!is_leaf) {
318        branchDists = new LeafRelations;
319
320        for (int s = 0; s<2; s++) {
321            ClusterTree         *son      = s ? get_leftson() : get_rightson();
322            const LeafRelations *sonDists = son->get_branch_dists();
323
324            cl_assert(sonDists || son->is_leaf);
325            if (sonDists) branchDists->insert(sonDists->begin(), sonDists->end());
326        }
327
328        const NodeValues *ldepths = get_leftson()->get_branch_depths();
329        const NodeValues *rdepths = get_rightson()->get_branch_depths();
330
331        NodeValues::const_iterator l_end = ldepths->end();
332        NodeValues::const_iterator r_end = rdepths->end();
333
334        for (NodeValues::const_iterator ld = ldepths->begin(); ld != l_end; ++ld) {
335            AP_FLOAT llen = ld->second+leftlen;
336            for (NodeValues::const_iterator rd = rdepths->begin(); rd != r_end; ++rd) {
337                AP_FLOAT rlen = rd->second+rightlen;
338
339                (*branchDists)[TwoLeafs(ld->first, rd->first)] = llen+rlen;
340            }
341        }
342
343        cl_assert(branchDists->size() == possible_relations());
344    }
345}
346
347AP_FLOAT ClusterTree::get_seqDist(const TwoLeafs& pair) {
348    // searches or calculates the distance between the sequences of leafs in 'pair'
349
350    const AP_FLOAT *found = has_seqDist(pair);
351    if (!found) {
352        const ClusterTree *commonFather = pair.first()->commonFatherWith(pair.second());
353
354        while (commonFather != this) {
355            found = commonFather->has_seqDist(pair);
356            if (found) break;
357
358            commonFather = commonFather->get_father();
359            cl_assert(commonFather); // first commonFather was not inside this!
360        }
361    }
362
363    if (!found) {
364        // calculate the distance
365
366        const AP_sequence *seq1 = pair.first()->get_seq();
367        const AP_sequence *seq2 = pair.second()->get_seq();
368
369        seq1->lazy_load_sequence();
370        seq2->lazy_load_sequence();
371
372        {
373            AP_sequence *ancestor     = seq1->dup();
374            AP_FLOAT     mutations    = ancestor->combine(seq1, seq2);
375            AP_FLOAT     wbc1         = seq1->weighted_base_count();
376            AP_FLOAT     wbc2         = seq2->weighted_base_count();
377            AP_FLOAT     minBaseCount = wbc1 < wbc2 ? wbc1 : wbc2;
378            AP_FLOAT     dist;
379
380            if (minBaseCount>0) {
381                dist = mutations / minBaseCount;
382                if (minBaseCount < min_bases) min_bases = minBaseCount;
383            }
384            else { // got at least one empty sequence -> no distance
385                dist      = 0.0;
386                min_bases = minBaseCount;
387            }
388
389#if defined(DEBUG) && 0
390            const char *name1 = pair.first()->name;
391            const char *name2 = pair.second()->name;
392            printf("%s <-?-> %s: mutations=%5.2f wbc1=%5.2f wbc2=%5.2f mbc=%5.2f dist=%5.2f\n",
393                   name1, name2, mutations, wbc1, wbc2, minBaseCount, dist);
394#endif // DEBUG
395
396            cl_assert(!is_nan_or_inf(dist));
397           
398            (*sequenceDists)[pair] = dist;
399            delete ancestor;
400        }
401
402        found = has_seqDist(pair);
403        cl_assert(found);
404    }
405
406    cl_assert(found);
407    return *found;
408}
409
410const AP_FLOAT *ClusterTree::has_seqDist(const TwoLeafs& pair) const {
411    if (sequenceDists) {
412        LeafRelationCIter found = sequenceDists->find(pair);
413        if (found != sequenceDists->end()) {
414            return &(found->second);
415        }
416    }
417    return NULL;
418}
419
420const ClusterTree *ClusterTree::commonFatherWith(const ClusterTree *other) const {
421    int depth_diff = get_depth() - other->get_depth();
422    if (depth_diff) {
423        if (depth_diff<0) { // this is less deep than other
424            return other->get_father()->commonFatherWith(this);
425        }
426        else { // other is less deep than this
427            return get_father()->commonFatherWith(other);
428        }
429    }
430    else { // both at same depth
431        if (this == other) { // common father reached ?
432            return this;
433        }
434        else { // otherwise check fathers
435            return get_father()->commonFatherWith(other->get_father());
436        }
437    }
438}
439
Note: See TracBrowser for help on using the repository browser.