source: tags/arb-6.0/MERGE/MG_preserves.cxx

Last change on this file was 12267, checked in by westram, 10 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.5 KB
Line 
1//  ==================================================================== //
2//                                                                       //
3//    File      : MG_preserves.cxx                                       //
4//    Purpose   : find candidates for alignment preservation             //
5//                                                                       //
6//                                                                       //
7//  Coded by Ralf Westram (coder@reallysoft.de) in July 2003             //
8//  Copyright Department of Microbiology (Technical University Munich)   //
9//                                                                       //
10//  Visit our web site at: http://www.arb-home.de/                       //
11//                                                                       //
12//                                                                       //
13//  ==================================================================== //
14
15#include "merge.hxx"
16#include "MG_adapt_ali.hxx"
17
18#include <aw_awar.hxx>
19#include <aw_root.hxx>
20#include <aw_msg.hxx>
21#include <aw_window.hxx>
22#include <aw_select.hxx>
23#include <arb_progress.h>
24#include <arb_global_defs.h>
25#include <arbdbt.h>
26
27#include <set>
28#include <string>
29
30using namespace std;
31
32// find species/SAIs to preserve alignment
33
34#define AWAR_REMAP_CANDIDATE     AWAR_MERGE_TMP "remap_candidates"
35#define AWAR_REMAP_ALIGNMENT     AWAR_MERGE_TMP "remap_alignment"
36#define AWAR_REMAP_SEL_REFERENCE AWAR_MERGE_TMP "remap_reference"
37
38struct preserve_para {
39    AW_selection_list *alignmentList;
40    AW_selection_list *refCandidatesList; // reference candidates
41    AW_selection_list *usedRefsList;
42};
43
44static void get_global_alignments(ConstStrArray& ali_names) {
45    // get all alignment names available in both databases
46    GBT_get_alignment_names(ali_names, GLOBAL_gb_src);
47    GBDATA *gb_presets = GBT_get_presets(GLOBAL_gb_dst);
48
49    int i;
50    for (i = 0; ali_names[i]; ++i) {
51        GBDATA *gb_ali_name = GB_find_string(gb_presets, "alignment_name", ali_names[i], GB_IGNORE_CASE, SEARCH_GRANDCHILD);
52        if (!gb_ali_name) ali_names.remove(i--);
53    }
54}
55
56static void init_alignments(preserve_para *para) {
57    // initialize the alignment selection list
58    ConstStrArray ali_names;
59    get_global_alignments(ali_names);
60    para->alignmentList->init_from_array(ali_names, "All");
61}
62
63static void clear_candidates(preserve_para *para) {
64    // clear the candidate list
65    AW_selection_list *candList = para->refCandidatesList;
66
67    candList->clear();
68    candList->insert_default(DISPLAY_NONE, NO_ALI_SELECTED);
69    candList->update();
70}
71
72static long count_bases(const char *data, long len) {
73    long count = 0;
74    for (long i = 0; i<len; ++i) {
75        if (data[i] != '-' && data[i] != '.') {
76            ++count;
77        }
78    }
79    return count;
80}
81
82// count real bases
83// (gb_data should point to ali_xxx/data)
84static long count_bases(GBDATA *gb_data) {
85    long count = 0;
86    switch (GB_read_type(gb_data)) {
87        case GB_STRING:
88            count = count_bases(GB_read_char_pntr(gb_data), GB_read_count(gb_data));
89            break;
90        case GB_BITS: {
91            char *bitstring = GB_read_as_string(gb_data);
92            long  len       = strlen(bitstring);
93            count           = count_bases(bitstring, len);
94            free(bitstring);
95            break;
96        }
97        default:
98            GB_export_errorf("Data type %s is not supported", GB_get_type_name(gb_data));
99            break;
100    }
101    return count;
102}
103
104// -------------------------
105//      class Candidate
106
107class Candidate {
108    string name;                // species/SAI name
109    double score;
110    int    found_alignments;    // counts alignments with data in both databases
111    long   base_count_diff;
112
113public:
114    Candidate(bool is_species, const char *name_, GBDATA *gb_src, GBDATA *gb_dst, const CharPtrArray& ali_names)
115        : name(is_species ? name_ : string("SAI:")+name_)
116    {
117        found_alignments = 0;
118        score            = 0.0;
119        base_count_diff  = 0;
120        bool valid       = true;
121
122        mg_assert(!GB_have_error());
123
124        for (int i = 0; valid && ali_names[i]; ++i) {
125            if (GBDATA *gb_src_data = GBT_read_sequence(gb_src, ali_names[i])) {
126                if (GBDATA *gb_dst_data = GBT_read_sequence(gb_dst, ali_names[i])) {
127                    ++found_alignments;
128
129                    long src_bases  = count_bases(gb_src_data);
130                    long dst_bases  = count_bases(gb_dst_data);
131                    base_count_diff = labs(src_bases-dst_bases);
132
133                    if (GB_have_error()) valid = false;
134
135                    double factor;
136                    if (base_count_diff>5)  factor = 0.2;
137                    else                    factor = (6-base_count_diff)*(1.0/6);
138
139                    long min_bases  = src_bases<dst_bases ? src_bases : dst_bases;
140                    score          += min_bases*factor;
141                }
142            }
143        }
144
145        if (!valid) found_alignments = 0;
146    }
147    ~Candidate() {}
148
149    int has_alignments() const { return found_alignments; }
150    double get_score() const { return score; }
151
152    const char *get_name() const { return name.c_str(); }
153    const char *get_entry() const {
154        return GBS_global_string("%24s  %i %li %f", get_name(), found_alignments, base_count_diff, score);
155    }
156
157    bool operator < (const Candidate& other) const { // sort highest score first
158        int ali_diff = found_alignments-other.found_alignments;
159        if (ali_diff) return ali_diff>0;
160        return score > other.score;
161    }
162};
163
164// use SmartPtr to memory-manage Candidate's
165static bool operator < (const SmartPtr<Candidate>& c1, const SmartPtr<Candidate>& c2) {
166    return c1->operator<(*c2);
167}
168typedef set< SmartPtr<Candidate> > Candidates;
169
170static void find_species_candidates(Candidates& candidates, const CharPtrArray& ali_names) {
171    // collect names of all species in source database
172    GB_HASH      *src_species = GBT_create_species_hash(GLOBAL_gb_src);
173    long          src_count   = GBS_hash_count_elems(src_species);
174    arb_progress  progress("Examining species", src_count);
175    bool          aborted     = false;
176
177    // find existing species in destination database
178    for (GBDATA *gb_dst_species = GBT_first_species(GLOBAL_gb_dst);
179         gb_dst_species && !aborted;
180         gb_dst_species = GBT_next_species(gb_dst_species))
181    {
182        const char *dst_name = GBT_read_name(gb_dst_species);
183
184        if (GBDATA *gb_src_species = (GBDATA*)GBS_read_hash(src_species, dst_name)) {
185            Candidate *cand = new Candidate(true, dst_name, gb_src_species, gb_dst_species, ali_names);
186
187            if (cand->has_alignments() && cand->get_score()>0.0) {
188                candidates.insert(cand);
189            }
190            else {
191                if (GB_have_error()) {
192                    aw_message(GBS_global_string("Invalid adaption candidate '%s' (%s)", dst_name, GB_await_error()));
193                }
194                delete cand;
195            }
196
197            progress.inc();
198            aborted = progress.aborted();
199        }
200    }
201
202    GBS_free_hash(src_species);
203    progress.done();
204}
205
206static void find_SAI_candidates(Candidates& candidates, const CharPtrArray& ali_names) {
207    // add all candidate SAIs to 'candidates'
208    GB_HASH      *src_SAIs  = GBT_create_SAI_hash(GLOBAL_gb_src);
209    long          src_count = GBS_hash_count_elems(src_SAIs);
210    arb_progress  progress("Examining SAIs", src_count);
211
212    // find existing SAIs in destination database
213    for (GBDATA *gb_dst_SAI = GBT_first_SAI(GLOBAL_gb_dst);
214         gb_dst_SAI;
215         gb_dst_SAI = GBT_next_SAI(gb_dst_SAI))
216    {
217        const char *dst_name = GBT_read_name(gb_dst_SAI);
218
219        if (GBDATA *gb_src_SAI = (GBDATA*)GBS_read_hash(src_SAIs, dst_name)) {
220            Candidate *cand = new Candidate(false, dst_name, gb_src_SAI, gb_dst_SAI, ali_names);
221
222            if (cand->has_alignments() && cand->get_score()>0.0) {
223                candidates.insert(cand);
224            }
225            else {
226                if (GB_have_error()) {
227                    aw_message(GBS_global_string("Invalid adaption candidate 'SAI:%s' (%s)", dst_name, GB_await_error()));
228                }
229                delete cand;
230            }
231
232            progress.inc();
233        }
234    }
235
236    GBS_free_hash(src_SAIs);
237}
238
239static void calculate_preserves_cb(AW_window *aww, preserve_para *para) {
240    // FIND button (rebuild candidates list)
241
242    GB_transaction ta_src(GLOBAL_gb_src);
243    GB_transaction ta_dst(GLOBAL_gb_dst);
244
245    clear_candidates(para);
246
247    const char *ali = aww->get_root()->awar(AWAR_REMAP_ALIGNMENT)->read_char_pntr();
248    Candidates  candidates;
249
250    arb_progress("Searching candidates");
251
252    // add candidates
253    {
254        ConstStrArray ali_names;
255        if (0 == strcmp(ali, "All")) {
256            get_global_alignments(ali_names);
257        }
258        else {
259            ali_names.put(ali);
260        }
261        find_SAI_candidates(candidates, ali_names);
262        find_species_candidates(candidates, ali_names);
263    }
264
265    int                   count       = 0;
266    Candidates::iterator  e           = candidates.end();
267    AW_selection_list    *refCandList = para->refCandidatesList;
268
269    for (Candidates::iterator i = candidates.begin();
270         i != e && count<5000;
271         ++i, ++count)
272    {
273        string name  = (*i)->get_name();
274        string shown = (*i)->get_entry();
275
276        refCandList->insert(shown.c_str(), name.c_str());
277    }
278
279    refCandList->update();
280}
281
282
283
284static void read_references(ConstStrArray& refs, AW_root *aw_root)  {
285    char *ref_string = aw_root->awar(AWAR_REMAP_SPECIES_LIST)->read_string();
286    GBT_splitNdestroy_string(refs, ref_string, " \n,;", true);
287}
288static void write_references(AW_root *aw_root, const CharPtrArray& ref_array) {
289    char *ref_string = GBT_join_names(ref_array, '\n');
290    aw_root->awar(AWAR_REMAP_SPECIES_LIST)->write_string(ref_string);
291    aw_root->awar(AWAR_REMAP_ENABLE)->write_int(ref_string[0] != 0);
292    free(ref_string);
293}
294static void select_reference(AW_root *aw_root, const char *ref_to_select) {
295    aw_root->awar(AWAR_REMAP_SEL_REFERENCE)->write_string(ref_to_select);
296}
297static char *get_selected_reference(AW_root *aw_root) {
298    return aw_root->awar(AWAR_REMAP_SEL_REFERENCE)->read_string();
299}
300
301static void refresh_reference_list_cb(AW_root *aw_root, preserve_para *para) {
302    ConstStrArray  refs;
303    read_references(refs, aw_root);
304    para->usedRefsList->init_from_array(refs, "");
305}
306
307static void add_selected_cb(AW_window *aww, preserve_para *para) {
308    // ADD button (add currently selected candidate to references)
309    AW_root       *aw_root = aww->get_root();
310    ConstStrArray  refs;
311    read_references(refs, aw_root);
312
313    char *candidate  = aw_root->awar(AWAR_REMAP_CANDIDATE)->read_string();
314    char *selected   = get_selected_reference(aw_root);
315    int   cand_index = GBT_names_index_of(refs, candidate);
316    int   sel_index  = GBT_names_index_of(refs, selected);
317
318    if (cand_index == -1) GBT_names_add(refs, sel_index+1, candidate);
319    else                  GBT_names_move(refs, cand_index, sel_index);
320
321    write_references(aw_root, refs);
322    select_reference(aw_root, candidate);
323
324    free(selected);
325    free(candidate);
326
327    para->refCandidatesList->move_selection(1); 
328}
329
330static void clear_references_cb(AW_window *aww) {
331    // CLEAR button (clear references)
332    aww->get_root()->awar(AWAR_REMAP_SPECIES_LIST)->write_string("");
333}
334
335static void del_reference_cb(AW_window *aww) {
336    AW_root       *aw_root = aww->get_root();
337    ConstStrArray  refs;
338    read_references(refs, aw_root);
339
340    char *selected  = get_selected_reference(aw_root);
341    int   sel_index = GBT_names_index_of(refs, selected);
342
343    if (sel_index >= 0) {
344        select_reference(aw_root, refs[sel_index+1]);
345        GBT_names_erase(refs, sel_index);
346        write_references(aw_root, refs);
347    }
348
349    free(selected);
350}
351
352static void lower_reference_cb(AW_window *aww) {
353    AW_root       *aw_root = aww->get_root();
354    ConstStrArray  refs;
355    read_references(refs, aw_root);
356
357    char *selected  = get_selected_reference(aw_root);
358    int   sel_index = GBT_names_index_of(refs, selected);
359
360    if (sel_index >= 0) {
361        GBT_names_move(refs, sel_index, sel_index+1);
362        write_references(aw_root, refs);
363    }
364
365    free(selected);
366}
367static void raise_reference_cb(AW_window *aww) {
368    AW_root       *aw_root = aww->get_root();
369    ConstStrArray  refs;
370    read_references(refs, aw_root);
371
372    char *selected  = get_selected_reference(aw_root);
373    int   sel_index = GBT_names_index_of(refs, selected);
374
375    if (sel_index > 0) {
376        GBT_names_move(refs, sel_index, sel_index-1);
377        write_references(aw_root, refs);
378    }
379
380    free(selected);
381}
382
383static void test_references_cb(AW_window *aww) {
384    char           *reference_species_names = aww->get_root()->awar(AWAR_REMAP_SPECIES_LIST)->read_string();
385    GB_transaction  tm(GLOBAL_gb_src);
386    GB_transaction  td(GLOBAL_gb_dst);
387   
388    MG_remaps test_mapping(GLOBAL_gb_src, GLOBAL_gb_dst, true, reference_species_names); // will raise aw_message's in case of problems
389
390    free(reference_species_names);
391}
392
393static void init_preserve_awars(AW_root *aw_root) {
394    aw_root->awar_string(AWAR_REMAP_ALIGNMENT,     "", GLOBAL_gb_dst);
395    aw_root->awar_string(AWAR_REMAP_CANDIDATE,     "", GLOBAL_gb_dst);
396    aw_root->awar_string(AWAR_REMAP_SEL_REFERENCE, "", GLOBAL_gb_dst);
397}
398
399AW_window *MG_create_preserves_selection_window(AW_root *aw_root) {
400    // SELECT PRESERVES window
401    init_preserve_awars(aw_root);
402
403    AW_window_simple *aws = new AW_window_simple;
404
405    aws->init(aw_root, "SELECT_PRESERVES", "Select adaption candidates");
406    aws->load_xfig("merge/preserves.fig");
407
408    aws->at("close"); aws->callback((AW_CB0)AW_POPDOWN);
409    aws->create_button("CLOSE", "CLOSE", "C");
410
411    aws->at("help");
412    aws->callback(makeHelpCallback("mg_preserve.hlp"));
413    aws->create_button("HELP", "HELP", "H");
414
415    // ----------
416
417    preserve_para *para = new preserve_para; // do not free (is passed to callback)
418
419    aws->at("ali");
420    para->alignmentList = aws->create_selection_list(AWAR_REMAP_ALIGNMENT, 10, 30, true);
421
422    // ----------
423
424    aws->at("adapt");
425    aws->label("Adapt alignments");
426    aws->create_toggle(AWAR_REMAP_ENABLE);
427
428    aws->at("reference");
429    para->usedRefsList = aws->create_selection_list(AWAR_REMAP_SEL_REFERENCE, 10, 30, true);
430
431    aws->button_length(8);
432
433    aws->at("clear");
434    aws->callback(clear_references_cb);
435    aws->create_button("CLEAR", "Clear", "C");
436
437    aws->at("del");
438    aws->callback(del_reference_cb);
439    aws->create_button("DEL", "Del", "L");
440
441    aws->at("up");
442    aws->callback(raise_reference_cb);
443    aws->create_button("UP", "Up", "U");
444
445    aws->at("down");
446    aws->callback(lower_reference_cb);
447    aws->create_button("DOWN", "Down", "D");
448
449    aws->at("test");
450    aws->callback(test_references_cb);
451    aws->create_button("TEST", "Test", "T");
452   
453    // ----------
454
455    aws->at("find");
456    aws->callback(makeWindowCallback(calculate_preserves_cb, para));
457    aws->create_autosize_button("FIND", "Find candidates", "F", 1);
458
459    aws->at("add");
460    aws->callback(makeWindowCallback(add_selected_cb, para));
461    aws->create_button("ADD", "Add", "A");
462
463    aws->at("candidate");
464    para->refCandidatesList = aws->create_selection_list(AWAR_REMAP_CANDIDATE, 10, 30, true);
465
466    {
467        GB_transaction ta_src(GLOBAL_gb_src);
468        GB_transaction ta_dst(GLOBAL_gb_dst);
469
470        init_alignments(para);
471        clear_candidates(para);
472    }
473
474    {
475        AW_awar *awar_list = aw_root->awar(AWAR_REMAP_SPECIES_LIST);
476        awar_list->add_callback(makeRootCallback(refresh_reference_list_cb, para));
477        awar_list->touch(); // trigger callback
478    }
479
480    return aws;
481}
482
483
484
485
Note: See TracBrowser for help on using the repository browser.