root/trunk/DIST/DI_matr.cxx

Revision 8748, 55.8 KB (checked in by westram, 5 hours ago)
  • use simpler way to declare attributed static/inline functions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : DI_matr.cxx                                       //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "di_protdist.hxx"
12#include "di_clusters.hxx"
13#include "dist.hxx"
14#include "di_view_matrix.hxx"
15#include "di_matr.hxx"
16#include "di_awars.hxx"
17
18#include <neighbourjoin.hxx>
19#include <AP_seq_dna.hxx>
20#include <AP_filter.hxx>
21#include <BI_helix.hxx>
22#include <CT_ctree.hxx>
23#include <ColumnStat.hxx>
24
25#include <awt.hxx>
26#include <awt_sel_boxes.hxx>
27#include <awt_filter.hxx>
28#include <awt_macro.hxx>
29
30#include <aw_preset.hxx>
31#include <aw_awars.hxx>
32#include <aw_file.hxx>
33#include <aw_msg.hxx>
34#include <arb_progress.h>
35#include <aw_root.hxx>
36
37#include <gui_aliview.hxx>
38
39#include <climits>
40#include <ctime>
41#include <cmath>
42#include <arb_sort.h>
43#include <arb_global_defs.h>
44
45// --------------------------------------------------------------------------------
46
47#define AWAR_DIST_BOOTSTRAP_COUNT AWAR_DIST_PREFIX "bootstrap/count"
48#define AWAR_DIST_CANCEL_CHARS    AWAR_DIST_PREFIX "cancel/chars"
49
50#define AWAR_DIST_COLUMN_STAT_NAME AWAR_DIST_COLUMN_STAT_PREFIX "name"
51
52#define AWAR_DIST_TREE_SORT_NAME "tmp/" AWAR_DIST_TREE_PREFIX "sort_tree_name"
53#define AWAR_DIST_TREE_COMP_NAME "tmp/" AWAR_DIST_TREE_PREFIX "compr_tree_name"
54#define AWAR_DIST_TREE_STD_NAME  AWAR_DIST_TREE_PREFIX "tree_name"
55
56#define AWAR_DIST_SAVE_MATRIX_TYPE     AWAR_DIST_SAVE_MATRIX_BASE "/type"
57#define AWAR_DIST_SAVE_MATRIX_FILENAME AWAR_DIST_SAVE_MATRIX_BASE "/file_name"
58
59// --------------------------------------------------------------------------------
60
61DI_GLOBAL_MATRIX GLOBAL_MATRIX;
62
63static DI_dmatrix *di_dmatrix = 0;
64static AP_matrix   DI_dna_matrix(AP_MAX);
65
66static void refresh_matrix_display() {
67    if (di_dmatrix) {
68        di_dmatrix->resized();
69        di_dmatrix->display(false);
70    }
71}
72
73static void delete_matrix_cb(AW_root *) {
74    GLOBAL_MATRIX.forget();
75}
76
77static AW_window *create_dna_matrix_window(AW_root *aw_root) {
78    AW_window_simple *aws = 0;
79    aws = new AW_window_simple;
80    aws->init(aw_root, "SET_DNA_MATRIX", "SET MATRIX");
81    aws->auto_increment(50, 50);
82    aws->button_length(10);
83    aws->callback(AW_POPDOWN);
84    aws->create_button("CLOSE", "CLOSE");
85
86    aws->callback(AW_POPUP_HELP, (AW_CL)"user_matrix.hlp");
87    aws->create_button("HELP", "HELP");
88
89    aws->at_newline();
90
91    DI_dna_matrix.create_input_fields(aws, AWAR_DIST_MATRIX_DNA_BASE);
92    aws->window_fit();
93    return aws;
94}
95
96void DI_create_matrix_variables(AW_root *aw_root, AW_default def, AW_default db)
97{
98    GB_transaction dummy(db);
99    DI_dna_matrix.set_description("", "");
100    DI_dna_matrix.x_description[AP_A] = strdup("A");
101    DI_dna_matrix.x_description[AP_C] = strdup("C");
102    DI_dna_matrix.x_description[AP_G] = strdup("G");
103    DI_dna_matrix.x_description[AP_T] = strdup("TU");
104    DI_dna_matrix.x_description[AP_S] = strdup("GAP");
105
106    DI_dna_matrix.y_description[AP_A] = strdup("A");
107    DI_dna_matrix.y_description[AP_C] = strdup("C");
108    DI_dna_matrix.y_description[AP_G] = strdup("G");
109    DI_dna_matrix.y_description[AP_T] = strdup("TU");
110    DI_dna_matrix.y_description[AP_S] = strdup("GAP");
111
112    DI_dna_matrix.create_awars(aw_root, AWAR_DIST_MATRIX_DNA_BASE);
113
114    aw_root->awar_int(AWAR_DIST_MATRIX_DNA_ENABLED, 0)->add_callback(delete_matrix_cb); // user matrix disabled by default
115    {
116        GBDATA *gbd = GB_search(AW_ROOT_DEFAULT, AWAR_DIST_MATRIX_DNA_BASE, GB_FIND);
117        GB_add_callback(gbd, GB_CB_CHANGED, (GB_CB)delete_matrix_cb, 0);
118    }
119
120    aw_root->awar_string("tmp/dummy_string", "0", def);
121
122    aw_root->awar_string(AWAR_DIST_WHICH_SPECIES, "marked", def)->add_callback(delete_matrix_cb);
123    {
124        char *default_ali = GBT_get_default_alignment(db);
125        aw_root->awar_string(AWAR_DIST_ALIGNMENT, "", def)     ->add_callback(delete_matrix_cb)->write_string(default_ali);
126        free(default_ali);
127    }
128    aw_root->awar_string(AWAR_DIST_FILTER_ALIGNMENT, "none", def);
129    aw_root->awar_string(AWAR_DIST_FILTER_NAME,      "none", def);
130    aw_root->awar_string(AWAR_DIST_FILTER_FILTER,    "",     def)->add_callback(delete_matrix_cb);
131    aw_root->awar_int   (AWAR_DIST_FILTER_SIMPLIFY,  0,      def)->add_callback(delete_matrix_cb);
132
133    aw_root->awar_string(AWAR_DIST_CANCEL_CHARS, ".", def)->add_callback(delete_matrix_cb);
134    aw_root->awar_int(AWAR_DIST_CORR_TRANS, (int)DI_TRANSFORMATION_SIMILARITY, def)->add_callback(delete_matrix_cb);
135
136    aw_root->awar(AWAR_DIST_FILTER_ALIGNMENT)->map(AWAR_DIST_ALIGNMENT);
137
138    AW_create_fileselection_awars(aw_root, AWAR_DIST_SAVE_MATRIX_BASE, ".", "", "infile", def);
139    aw_root->awar_int(AWAR_DIST_SAVE_MATRIX_TYPE, 0, def);
140
141    aw_root->awar_string(AWAR_DIST_TREE_STD_NAME,  "tree_nj", def);
142    {
143        char *currentTree = aw_root->awar_string(AWAR_TREE, "", db)->read_string();
144        aw_root->awar_string(AWAR_DIST_TREE_CURR_NAME, currentTree, def);
145        aw_root->awar_string(AWAR_DIST_TREE_SORT_NAME, currentTree, def)->add_callback(delete_matrix_cb);
146        free(currentTree);
147    }
148    aw_root->awar_string(AWAR_DIST_TREE_COMP_NAME, NO_TREE_SELECTED, def)->add_callback(delete_matrix_cb);
149
150    aw_root->awar_int(AWAR_DIST_BOOTSTRAP_COUNT, 1000, def);
151
152    aw_root->awar_float(AWAR_DIST_MIN_DIST, 0.0);
153    aw_root->awar_float(AWAR_DIST_MAX_DIST, 0.0);
154
155    aw_root->awar_string(AWAR_SPECIES_NAME, "", db);
156
157    DI_create_cluster_awars(aw_root, def, db);
158
159#if defined(DEBUG)
160    AWT_create_db_browser_awars(aw_root, def);
161#endif // DEBUG
162
163    {
164        GB_push_transaction(db);
165
166        GBDATA *gb_species_data = GBT_get_species_data(db);
167        GB_add_callback(gb_species_data, GB_CB_CHANGED, (GB_CB)delete_matrix_cb, 0);
168
169        GB_pop_transaction(db);
170    }
171}
172
173DI_ENTRY::DI_ENTRY(GBDATA *gbd, DI_MATRIX *phmatri) {
174    memset((char *)this, 0, sizeof(DI_ENTRY));
175    phmatrix = phmatri;
176
177    GBDATA *gb_ali = GB_entry(gbd, phmatrix->get_aliname());
178    if (gb_ali) {
179        GBDATA *gb_data = GB_entry(gb_ali, "data");
180        if (gb_data) {
181            if (phmatrix->is_AA) {
182                sequence = sequence_protein = new AP_sequence_simple_protein(phmatrix->get_aliview());
183            }
184            else {
185                sequence = sequence_parsimony = new AP_sequence_parsimony(phmatrix->get_aliview());
186            }
187            sequence->bind_to_species(gbd);
188            sequence->lazy_load_sequence(); // load sequence
189
190            name      = GBT_read_string(gbd, "name");
191            full_name = GBT_read_string(gbd, "full_name");
192        }
193    }
194}
195
196DI_ENTRY::DI_ENTRY(char *namei, DI_MATRIX *phmatri)
197{
198    memset((char *)this, 0, sizeof(DI_ENTRY));
199    phmatrix = phmatri;
200    this->name = strdup(namei);
201}
202
203DI_ENTRY::~DI_ENTRY()
204{
205    delete sequence_protein;
206    delete sequence_parsimony;
207    free(name);
208    free(full_name);
209
210}
211
212DI_MATRIX::DI_MATRIX(const AliView& aliview_, AW_root *awr) {
213    memset((char *)this, 0, sizeof(*this));
214    aw_root = awr;
215    aliview = new AliView(aliview_);
216}
217
218char *DI_MATRIX::unload() {
219    for (long i=0; i<nentries; i++) {
220        delete entries[i];
221    }
222    freenull(entries);
223    nentries = 0;
224    return 0;
225}
226
227DI_MATRIX::~DI_MATRIX()
228{
229    unload();
230    delete matrix;
231    delete aliview;
232}
233
234struct TreeOrderedSpecies {
235    GBDATA  *gbd;
236    int      order_index;
237
238    TreeOrderedSpecies(const MatrixOrder& order, GBDATA *gb_spec)
239        : gbd(gb_spec),
240          order_index(order.get_index(GBT_read_name(gbd)))
241    {}
242};
243
244MatrixOrder::MatrixOrder(GBDATA *gb_main, GB_CSTR sort_tree_name)
245    : name2pos(NULL),
246      leafs(0)
247{
248    if (sort_tree_name) {
249        int size;
250        GBT_TREE *sort_tree = GBT_read_tree_and_size(gb_main, sort_tree_name, sizeof(GBT_TREE), &size);
251
252        if (sort_tree) {
253            leafs    = size+1;
254            name2pos = GBS_create_hash(leafs, GB_IGNORE_CASE);
255
256            IF_DEBUG(int leafsLoaded = leafs);
257            leafs = 0;
258            insert_in_hash(sort_tree);
259
260            arb_assert(leafsLoaded == leafs);
261        }
262        else {
263            GB_clear_error();
264        }
265    }
266}
267static int TreeOrderedSpecies_cmp(const void *p1, const void *p2, void *) {
268    TreeOrderedSpecies *s1 = (TreeOrderedSpecies*)p1;
269    TreeOrderedSpecies *s2 = (TreeOrderedSpecies*)p2;
270
271    return s2->order_index - s1->order_index;
272}
273
274void MatrixOrder::applyTo(TreeOrderedSpecies **species_array, size_t array_size) const {
275    GB_sort((void**)species_array, 0, array_size, TreeOrderedSpecies_cmp, NULL);
276}
277
278char *DI_MATRIX::load(LoadWhat what, const MatrixOrder& order, bool show_warnings, GBDATA **species_list) {
279    GBDATA     *gb_main = get_gb_main();
280    const char *use     = get_aliname();
281
282    GB_push_transaction(gb_main);
283
284    seq_len          = GBT_get_alignment_len(gb_main, use);
285    is_AA            = GBT_is_alignment_protein(gb_main, use);
286    gb_species_data  = GBT_get_species_data(gb_main);
287    entries_mem_size = 1000;
288
289    entries = (DI_ENTRY **)calloc(sizeof(DI_ENTRY*), (size_t)entries_mem_size);
290
291    nentries = 0;
292
293    int no_of_species = -1;
294    switch (what) {
295        case DI_LOAD_ALL:
296            no_of_species = GBT_get_species_count(gb_main);
297            break;
298        case DI_LOAD_MARKED:
299            no_of_species = GBT_count_marked_species(gb_main);
300            break;
301        case DI_LOAD_LIST:
302            di_assert(species_list);
303            for (no_of_species = 0; species_list[no_of_species]; ++no_of_species) ;
304            break;
305    }
306
307    TreeOrderedSpecies *species_to_load[no_of_species];
308
309    {
310        int i = 0;
311        switch (what) {
312            case DI_LOAD_ALL: {
313                for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_species(gb_species), ++i) {
314                    species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
315                }
316                break;
317            }
318            case DI_LOAD_MARKED: {
319                for (GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_marked_species(gb_species), ++i) {
320                    species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
321                }
322                break;
323            }
324            case DI_LOAD_LIST: {
325                for (i = 0; species_list[i]; ++i) {
326                    species_to_load[i] = new TreeOrderedSpecies(order, species_list[i]);
327                }
328                break;
329            }
330        }
331        arb_assert(i == no_of_species);
332    }
333
334    if (order.defined()) {
335        order.applyTo(species_to_load, no_of_species);
336        if (show_warnings) {
337            int species_not_in_sort_tree = 0;
338            for (int i = 0; i<no_of_species; ++i) {
339                if (!species_to_load[i]->order_index) {
340                    species_not_in_sort_tree++;
341                }
342            }
343            if (species_not_in_sort_tree) {
344                aw_message(GBS_global_string("%i of the affected species are not in sort-tree", species_not_in_sort_tree));
345            }
346        }
347    }
348    else {
349        if (show_warnings) {
350            aw_message("No valid tree given to sort matrix (using default database order)");
351        }
352    }
353
354    if (no_of_species>entries_mem_size) {
355        entries_mem_size = no_of_species;
356        entries = (DI_ENTRY **)realloc(entries, (size_t)(sizeof(DI_ENTRY*)*entries_mem_size));
357    }
358
359    for (int i = 0; i<no_of_species; ++i) {
360        DI_ENTRY *phentry = new DI_ENTRY(species_to_load[i]->gbd, this);
361        if (phentry->sequence) {    // a species found
362            arb_assert(nentries<entries_mem_size);
363            entries[nentries++] = phentry;
364        }
365        else {
366            delete phentry;
367        }
368        delete species_to_load[i];
369        species_to_load[i] = NULL;
370    }
371
372    GB_pop_transaction(gb_main);
373    return 0;
374}
375
376void DI_MATRIX::clear(DI_MUT_MATR &hits)
377{
378    int i, j;
379    for (i=0; i<AP_MAX; i++) {
380        for (j=0; j<AP_MAX; j++) {
381            hits[i][j] = 0;
382        }
383    }
384}
385
386void DI_MATRIX::make_sym(DI_MUT_MATR &hits)
387{
388    int i, j;
389    for (i=AP_A; i<AP_MAX; i*=2) {
390        for (j=AP_A; j<=i; j*=2) {
391            hits[i][j] = hits[j][i] = hits[i][j] + hits[j][i];
392        }
393    }
394}
395
396void DI_MATRIX::rate_write(DI_MUT_MATR &hits, FILE *out) {
397    int i, j;
398    for (i=AP_A; i<AP_MAX; i*=2) {
399        for (j=AP_A; j<AP_MAX; j*=2) {
400            fprintf(out, "%5li ", hits[i][j]);
401        }
402        fprintf(out, "\n");
403    }
404}
405
406long *DI_MATRIX::create_helix_filter(BI_helix *helix, const AP_filter *filter) {
407    long   *result = (long *)calloc(sizeof(long), filter->get_filtered_length());
408    long   *temp   = (long *)calloc(sizeof(long), filter->get_filtered_length());
409    int     count  = 0;
410    size_t  i;
411
412    for (i=0; i<filter->get_length(); i++) {
413        if (filter->use_position(i)) {
414            temp[i] = count;
415            if (helix->pairtype(i) == HELIX_PAIR) {
416                result[count] = helix->opposite_position(i);
417            }
418            else {
419                result[count] = -1;
420            }
421            count++;
422        }
423    }
424    while (--count >= 0) {
425        if (result[count] >= 0) {
426            result[count] = temp[result[count]];
427        }
428    }
429    free(temp);
430    return result;
431}
432
433GB_ERROR DI_MATRIX::calculate_rates(DI_MUT_MATR &hrates, DI_MUT_MATR &nrates, DI_MUT_MATR &pairs, long *filter) {
434    if (nentries<=1) {
435        return "There are no species selected";
436    }
437
438    arb_progress progress("rates", (nentries*(nentries+1))/2);
439    GB_ERROR     error = NULL;
440
441    long s_len = aliview->get_length();
442
443    this->clear(hrates);
444    this->clear(nrates);
445    this->clear(pairs);
446
447    for (long row = 0; row<nentries && !error; row++) {
448        for (long col=0; col<row; col++) {
449            const unsigned char *seq1 = entries[row]->sequence_parsimony->get_usequence();
450            const unsigned char *seq2 = entries[col]->sequence_parsimony->get_usequence();
451            for (long pos = 0; pos < s_len; pos++) {
452                if (filter[pos]>=0) {
453                    hrates[*seq1][*seq2]++;
454                }
455                else {
456                    nrates[*seq1][*seq2]++;
457                }
458                seq1++; seq2++;
459            }
460            progress.inc_and_check_user_abort(error);
461        }
462    }
463    for (long row = 0; row<nentries; row++) {
464        const unsigned char *seq1 = entries[row]->sequence_parsimony->get_usequence();
465        for (long pos = 0; pos < s_len; pos++) {
466            if (filter[pos]>=0) {
467                pairs[seq1[pos]][seq1[filter[pos]]]++;
468            }
469        }
470    }
471    return error;
472}
473
474// ----------------------------------------------------------------
475//      Some test functions to check correlated base correction
476
477GB_ERROR DI_MATRIX::haeschoe(const char *path) {
478    static BI_helix *helix = 0;
479    if (!helix) helix      = new BI_helix();
480
481    GB_ERROR error = helix->init(get_gb_main());
482    if (!error) {
483        FILE *out = fopen(path, "w");
484        if (!out) {
485            GB_export_IO_error("writing", path);
486            error = GB_await_error();
487        }
488        else {
489            arb_progress progress("Calculating distance matrix");
490
491            fprintf(out, "Pairs in helical regions:\n");
492            long *filter = create_helix_filter(helix, aliview->get_filter());
493           
494            DI_MUT_MATR temp, temp2, pairs;
495            error = calculate_rates(temp, temp2, pairs, filter);
496            if (!error) {
497                rate_write(pairs, out);
498                make_sym(temp); make_sym(temp2);
499                fprintf(out, "\nRatematrix helical parts:\n");
500                rate_write(temp, out);
501                fprintf(out, "\nRatematrix non helical parts:\n");
502                rate_write(temp2, out);
503
504                long row, col, pos, s_len;
505
506                s_len = aliview->get_length();
507                fprintf(out, "\nDistance matrix (Helixdist Helixlen Nonhelixdist Nonhelixlen):");
508                fprintf(out, "\n%li", nentries);
509
510                const int MAXDISTDEBUG = 1000;
511                double    distdebug[MAXDISTDEBUG];
512
513                for (pos = 0; pos<MAXDISTDEBUG; pos++) distdebug[pos] = 0.0;
514
515                arb_progress dist_progress("distance", (nentries*(nentries+1))/2);
516                for (row = 0; row<nentries && !error; row++) {
517                    fprintf (out, "\n%s  ", entries[row]->name);
518
519                    for (col=0; col<row && !error; col++) {
520                        const unsigned char *seq1, *seq2;
521
522                        seq1 = entries[row]->sequence_parsimony->get_usequence();
523                        seq2 = entries[col]->sequence_parsimony->get_usequence();
524                        this->clear(temp);
525                        this->clear(temp2);
526
527                        for (pos = 0; pos < s_len; pos++) {
528                            if (filter[pos]>=0) temp [*seq1][*seq2]++;
529                            else                temp2[*seq1][*seq2]++;
530                            seq1++; seq2++;
531                        }
532                        long hdist = 0, hall2 = 0;
533                        int i, j;
534                        for (i=AP_A; i<AP_MAX; i*=2) {
535                            for (j=AP_A; j<AP_MAX; j*=2) {
536                                hall2 += temp[i][j];
537                                if (i!=j) hdist += temp[i][j];
538                            }
539                        }
540
541                        long dist = 0, all = 0;
542                        for (i=AP_A; i<=AP_T; i*=2) {
543                            for (j=AP_A; j<=AP_T; j*=2) {
544                                all += temp2[i][j];
545                                if (i!=j) dist += temp2[i][j];
546                            }
547                        }
548                        fprintf(out, "(%4li:%4li %4li:%4li) ", hdist, hall2, dist, all);
549                        if (all>100) {
550                            distdebug[dist*MAXDISTDEBUG/all] = (double)hdist/(double)hall2;
551                        }
552                        dist_progress.inc_and_check_user_abort(error);
553                    }
554                }
555
556                if (!error) {
557                    fprintf (out, "\n");
558                    fprintf (out, "\nhdist/dist:\n");
559
560                    for (pos = 1; pos<MAXDISTDEBUG; pos++) {
561                        if (distdebug[pos]) {
562                            fprintf(out, "%4f %5f\n", (double)pos/(double)MAXDISTDEBUG, distdebug[pos]/(double)pos*(double)MAXDISTDEBUG);
563                        }
564                    }
565                }
566            }
567           
568            fclose(out);
569        }
570    }
571    return error;
572}
573
574char *DI_MATRIX::calculate_overall_freqs(double rel_frequencies[AP_MAX], char *cancel)
575{
576    long hits2[AP_MAX];
577    long sum   = 0;
578    int  i, row;
579    int  pos;
580    int  b;
581    long s_len = aliview->get_length();
582
583    memset((char *) &hits2[0], 0, sizeof(hits2));
584    for (row = 0; row < nentries; row++) {
585        const char *seq1 = entries[row]->sequence_parsimony->get_sequence();
586        for (pos = 0; pos < s_len; pos++) {
587            b = *(seq1++);
588            if (cancel[b]) continue;
589            hits2[b]++;
590        }
591    }
592    for (i = 0; i < AP_MAX; i++)    sum += hits2[i];
593    for (i = 0; i < AP_MAX; i++)    rel_frequencies[i] = hits2[i] / (double) sum;
594    return 0;
595}
596
597double DI_MATRIX::corr(double dist, double b, double & sigma) {
598    const double eps = 0.01;
599    double ar = 1.0 - dist/b;
600    sigma = 1000.0;
601    if (ar< eps) return 3.0;
602    sigma = b/ar;
603    return - b * log(1-dist/b);
604}
605
606GB_ERROR DI_MATRIX::calculate(AW_root *awr, char *cancel, double /* alpha */, DI_TRANSFORMATION transformation, bool *aborted_flag)
607{
608    if (transformation == DI_TRANSFORMATION_HAESCH) {
609        GB_ERROR error = haeschoe("outfile");
610        if (error) return error;
611        return "Your matrices have been written to 'outfile'\nSorry I can not make a tree";
612    }
613    int user_matrix_enabled = awr->awar(AWAR_DIST_MATRIX_DNA_ENABLED)->read_int();
614    if (user_matrix_enabled) {  // set transformation Matrix
615        switch (transformation) {
616            case DI_TRANSFORMATION_NONE:
617            case DI_TRANSFORMATION_SIMILARITY:
618            case DI_TRANSFORMATION_JUKES_CANTOR:
619                break;
620            default:
621                return GB_export_error("Sorry not implemented:\n"
622                                       "    This kind of distance correction does not support\n"
623                                       "    a user defined matrix");
624        }
625        DI_dna_matrix.read_awars(awr, AWAR_DIST_MATRIX_DNA_BASE);
626        DI_dna_matrix.normize();
627    }
628
629
630
631    matrix = new AP_smatrix(nentries);
632
633    long   s_len = aliview->get_length();
634    long   hits[AP_MAX][AP_MAX];
635    int    row;
636    int    col;
637    size_t i;
638
639    if (nentries<=1) {
640        return "There are no species selected";
641    }
642    memset(&cancel_columns[0], 0, 256);
643
644    for (i=0; i<strlen(cancel); i++) {
645        int j = cancel[i];
646        j = AP_sequence_parsimony::table[j];
647        cancel_columns[j] = 1;
648    }
649    long    columns;
650    double b;
651    long frequencies[AP_MAX];
652    double rel_frequencies[AP_MAX];
653    double S_square=0;
654
655    switch (transformation) {
656        case DI_TRANSFORMATION_FELSENSTEIN:
657            this->calculate_overall_freqs(rel_frequencies, cancel_columns);
658            S_square = 0.0;
659            for (i=0; i<AP_MAX; i++) S_square += rel_frequencies[i]*rel_frequencies[i];
660            break;
661        default:    break;
662    };
663
664    arb_progress progress("Calculating distance matrix", (nentries*(nentries+1))/2);
665    GB_ERROR     error = NULL;
666    for (row = 0; row<nentries && !error; row++) {
667        for (col=0; col<=row && !error; col++) {
668            columns = 0;
669           
670            const unsigned char *seq1 = entries[row]->sequence_parsimony->get_usequence();
671            const unsigned char *seq2 = entries[col]->sequence_parsimony->get_usequence();
672
673            b = 0.0;
674            switch (transformation) {
675                case  DI_TRANSFORMATION_JUKES_CANTOR:
676                    b = 0.75;
677                    // fall-through
678                case  DI_TRANSFORMATION_NONE:
679                case  DI_TRANSFORMATION_SIMILARITY: {
680                    double  dist = 0.0;
681                    if (user_matrix_enabled) {
682                        memset((char *)hits, 0, sizeof(long) * AP_MAX * AP_MAX);
683                        int pos;
684                        for (pos = s_len; pos >= 0; pos--) {
685                            hits[*(seq1++)][*(seq2++)]++;
686                        }
687                        int x, y;
688                        double diffsum = 0.0;
689                        double all_sum = 0.001;
690                        for (x = AP_A; x < AP_MAX; x*=2) {
691                            for (y = AP_A; y < AP_MAX; y*=2) {
692                                if (x==y) {
693                                    all_sum += hits[x][y];
694                                }
695                                else {
696                                    diffsum += hits[x][y] * DI_dna_matrix.get(x, y);
697                                    all_sum += hits[x][y] * DI_dna_matrix.get(x, y);
698                                }
699                            }
700                        }
701                        dist = diffsum / all_sum;
702                    }
703                    else {
704                        int pos;
705                        int b1, b2;
706                        for (pos = s_len; pos >= 0; pos--) {
707                            b1 = *(seq1++);
708                            b2 = *(seq2++);
709                            if (cancel_columns[b1]) continue;
710                            if (cancel_columns[b2]) continue;
711                            columns++;
712                            if (b1&b2) continue;
713                            dist+=1.0;
714                        }
715                        if (columns == 0) columns = 1;
716                        dist /= columns;
717                    }
718                    if (transformation==DI_TRANSFORMATION_SIMILARITY) {
719                        dist =  (1.0-dist);
720                    }
721                    else if (b) {
722                        double sigma;
723                        dist = this->corr(dist, b, sigma);
724                    }
725                    matrix->set(row, col, dist);
726                    break;
727                }
728                case DI_TRANSFORMATION_KIMURA:
729                case DI_TRANSFORMATION_OLSEN:
730                case DI_TRANSFORMATION_FELSENSTEIN: {
731                    int    pos;
732                    double dist = 0.0;
733                    long   N, P, Q, M;
734                    double p, q;
735
736                    memset((char *)hits, 0, sizeof(long) * AP_MAX * AP_MAX);
737                    for (pos = s_len; pos >= 0; pos--) {
738                        hits[*(seq1++)][*(seq2++)]++;
739                    }
740                    switch (transformation) {
741                        case DI_TRANSFORMATION_KIMURA:
742                            P = hits[AP_A][AP_G] +
743                                hits[AP_G][AP_A] +
744                                hits[AP_C][AP_T] +
745                                hits[AP_T][AP_C];
746                            Q = hits[AP_A][AP_C] +
747                                hits[AP_A][AP_T] +
748                                hits[AP_C][AP_A] +
749                                hits[AP_T][AP_A] +
750                                hits[AP_G][AP_C] +
751                                hits[AP_G][AP_T] +
752                                hits[AP_C][AP_G] +
753                                hits[AP_T][AP_G];
754                            M = hits[AP_A][AP_A] +
755                                hits[AP_C][AP_C] +
756                                hits[AP_G][AP_G] +
757                                hits[AP_T][AP_T];
758                            N = P+Q+M;
759                            if (N==0) N=1;
760                            p = (double)P/(double)N;
761                            q = (double)Q/(double)N;
762                            dist = - .5 * log(
763                                              (1.0-2.0*p-q)*sqrt(1.0-2.0*q)
764                                              );
765                            break;
766
767                        case DI_TRANSFORMATION_OLSEN:
768                        case DI_TRANSFORMATION_FELSENSTEIN:
769
770                            memset((char *)frequencies, 0,
771                                   sizeof(long) * AP_MAX);
772
773                            N = 0;
774                            M = 0;
775
776                            for (i=0; i<AP_MAX; i++) {
777                                if (cancel_columns[i]) continue;
778                                unsigned int j;
779                                for (j=0; j<i; j++) {
780                                    if (cancel_columns[j]) continue;
781                                    frequencies[i] +=
782                                        hits[i][j]+
783                                        hits[j][i];
784                                }
785                                frequencies[i] += hits[i][i];
786                                N += frequencies[i];
787                                M += hits[i][i];
788                            }
789                            if (N==0) N=1;
790                            if (transformation == DI_TRANSFORMATION_OLSEN) { // Calc sum square freq individually for each line
791                                S_square = 0.0;
792                                for (i=0; i<AP_MAX; i++) S_square += frequencies[i]*frequencies[i];
793                                b = 1.0 - S_square/((double)N*(double)N);
794                            }
795                            else {
796                                b = 1.0 - S_square;
797                            }
798
799                            dist = ((double)(N-M)) / (double) N;
800                            double sigma;
801                            dist = this->corr(dist, b, sigma);
802                            break;
803
804                        default: return "Sorry: Transformation not implemented";
805                    }
806                    matrix->set(row, col, dist);
807                    break;
808                }
809                default:;
810            }   // switch
811            progress.inc_and_check_user_abort(error);
812        }
813    }
814    if (aborted_flag && progress.aborted()) *aborted_flag = true;
815    return error;
816}
817
818GB_ERROR DI_MATRIX::calculate_pro(DI_TRANSFORMATION transformation, bool *aborted_flag) {
819    di_cattype whichcat;
820    di_codetype whichcode = universal;
821
822    switch (transformation) {
823        case DI_TRANSFORMATION_NONE:
824            whichcat = none;
825            break;
826        case DI_TRANSFORMATION_SIMILARITY:
827            whichcat = similarity;
828            break;
829        case DI_TRANSFORMATION_KIMURA:
830            whichcat = kimura;
831            break;
832        case DI_TRANSFORMATION_PAM:
833            whichcat = pam;
834            break;
835        case DI_TRANSFORMATION_CATEGORIES_HALL:
836            whichcat = hall;
837            break;
838        case DI_TRANSFORMATION_CATEGORIES_BARKER:
839            whichcat = george;
840            break;
841        case DI_TRANSFORMATION_CATEGORIES_CHEMICAL:
842            whichcat = chemical;
843            break;
844        default:
845            return "This correction is not available for protein data";
846    }
847    matrix = new AP_smatrix(nentries);
848
849    di_protdist prodist(whichcode, whichcat, nentries, entries, aliview->get_length(), matrix);
850    return prodist.makedists(aborted_flag);
851}
852
853__ATTR__USERESULT static GB_ERROR di_calculate_matrix(AW_window *aww, const WeightedFilter *weighted_filter, bool bootstrap_flag, bool show_warnings, bool *aborted_flag) {
854    // sets 'aborted_flag' to true, if it is non-NULL and the calculation has been aborted
855    GB_push_transaction(GLOBAL_gb_main);
856    if (GLOBAL_MATRIX.exists()) {
857        GB_pop_transaction(GLOBAL_gb_main);
858        return 0;
859    }
860    AW_root *aw_root = aww->get_root();
861    char    *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
862    long ali_len = GBT_get_alignment_len(GLOBAL_gb_main, use);
863    if (ali_len<=0) {
864        GB_pop_transaction(GLOBAL_gb_main);
865        delete use;
866        aw_message("Please select a valid alignment");
867        return "Error";
868    }
869    arb_progress progress("Calculating matrix");
870   
871    char *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
872
873    AliView *aliview = weighted_filter->create_aliview(use);
874    if (bootstrap_flag) aliview->get_filter()->enable_bootstrap();
875
876    char *load_what      = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
877    char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
878
879    LoadWhat all_flag = (strcmp(load_what, "all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
880    GB_ERROR error    = 0;
881    {
882        DI_MATRIX *phm   = new DI_MATRIX(*aliview, aw_root);
883        phm->matrix_type = DI_MATRIX_FULL;
884
885        static char        *last_sort_tree_name = 0;
886        static MatrixOrder *last_order          = 0;
887
888        if (!last_sort_tree_name || !sort_tree_name || strcmp(last_sort_tree_name, sort_tree_name) != 0) {
889            last_sort_tree_name = nulldup(sort_tree_name);
890            last_order = new MatrixOrder(GLOBAL_gb_main, sort_tree_name);
891        }
892        di_assert(last_order);
893        phm->load(all_flag, *last_order, show_warnings, NULL);
894
895        free(sort_tree_name);
896        GB_pop_transaction(GLOBAL_gb_main);
897        bool aborted = false;
898        if (progress.aborted()) {
899            phm->unload();
900            error   = "Aborted by user";
901            aborted = true;
902        }
903        else {
904            DI_TRANSFORMATION trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
905
906            if (phm->is_AA) error = phm->calculate_pro(trans, &aborted);
907            else error            = phm->calculate(aw_root, cancel, 0.0, trans, &aborted);
908        }
909
910        if (aborted) {
911            di_assert(error);
912            if (aborted_flag) *aborted_flag = true;
913        }
914        if (error) {
915            delete phm;
916            GLOBAL_MATRIX.forget();
917        }
918        else {
919            GLOBAL_MATRIX.replaceBy(phm);
920        }
921    }
922    free(load_what);
923
924    free(cancel);
925    delete aliview;
926
927    free(use);
928    return error;
929}
930
931static void di_mark_by_distance(AW_window *aww, AW_CL cl_weightedFilter) {
932    AW_root *aw_root    = aww->get_root();
933    double   lowerBound = aw_root->awar(AWAR_DIST_MIN_DIST)->read_float();
934    double   upperBound = aw_root->awar(AWAR_DIST_MAX_DIST)->read_float();
935
936    GB_ERROR error = 0;
937    if (lowerBound >= upperBound) {
938        error = GBS_global_string("Lower bound (%f) has to be smaller than upper bound (%f)", lowerBound, upperBound);
939    }
940    else if (lowerBound<0.0 || lowerBound > 1.0) {
941        error = GBS_global_string("Lower bound (%f) is not in allowed range [0.0 .. 1.0]", lowerBound);
942    }
943    else if (upperBound<0.0 || upperBound > 1.0) {
944        error = GBS_global_string("Upper bound (%f) is not in allowed range [0.0 .. 1.0]", upperBound);
945    }
946    else {
947        GB_transaction ta(GLOBAL_gb_main);
948
949        char   *selected    = aw_root->awar(AWAR_SPECIES_NAME)->read_string();
950        GBDATA *gb_selected = 0;
951
952        if (!selected[0]) {
953            error = "Please select a species";
954        }
955        else {
956            gb_selected = GBT_find_species(GLOBAL_gb_main, selected);
957            if (!gb_selected) {
958                error = GBS_global_string("Couldn't find species '%s'", selected);
959            }
960            else {
961                char              *use             = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
962                char              *cancel          = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
963                DI_TRANSFORMATION  trans           = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
964                WeightedFilter    *weighted_filter = (WeightedFilter*)cl_weightedFilter;
965                AliView           *aliview         = weighted_filter->create_aliview(use);
966
967                if (!error) {
968                    DI_MATRIX *prev_global = GLOBAL_MATRIX.swap(NULL);
969
970                    size_t speciesCount   = GBT_get_species_count(GLOBAL_gb_main);
971                    bool   markedSelected = false;
972
973                    arb_progress progress("Mark species by distance", speciesCount);
974                    MatrixOrder  order(GLOBAL_gb_main, NULL);
975
976                    for (GBDATA *gb_species = GBT_first_species(GLOBAL_gb_main);
977                         gb_species && !error;
978                         gb_species = GBT_next_species(gb_species))
979                    {
980                        DI_MATRIX *phm         = new DI_MATRIX(*aliview, aw_root);
981                        phm->matrix_type       = DI_MATRIX_FULL;
982                        GBDATA *species_pair[] = { gb_selected, gb_species, NULL };
983
984                        phm->load(DI_LOAD_LIST, order, false, species_pair);
985
986                        if (phm->is_AA) {
987                            error = phm->calculate_pro(trans, NULL);
988                        }
989                        else {
990                            error = phm->calculate(aw_root, cancel, 0.0, trans, NULL);
991                        }
992
993                        double dist_value = phm->matrix->get(0, 1);                         // distance or conformance
994                        bool   mark       = (lowerBound <= dist_value && dist_value <= upperBound);
995                        GB_write_flag(gb_species, mark);
996
997                        if (!markedSelected) {
998                            dist_value = phm->matrix->get(0, 0);                                     // distance or conformance to self
999                            mark       = (lowerBound <= dist_value && dist_value <= upperBound);
1000                            GB_write_flag(gb_selected, mark);
1001
1002                            markedSelected = true;
1003                        }
1004
1005                        delete phm;
1006                        progress.inc_and_check_user_abort(error);
1007                    }
1008
1009                    di_assert(!GLOBAL_MATRIX.exists());
1010                    ASSERT_RESULT(DI_MATRIX*, NULL, GLOBAL_MATRIX.swap(prev_global));
1011                }
1012
1013                delete aliview;
1014                free(cancel);
1015                free(use);
1016            }
1017        }
1018
1019        free(selected);
1020        error = ta.close(error);
1021    }
1022
1023    if (error) {
1024        aw_message(error);
1025    }
1026}
1027
1028static void di_view_matrix_cb(AW_window *aww, AW_CL cl_sparam) {
1029    save_matrix_params *sparam = (save_matrix_params*)cl_sparam;
1030    GB_ERROR            error  = di_calculate_matrix(aww, sparam->weighted_filter, 0, true, NULL);
1031    if (error) return;
1032
1033    if (!di_dmatrix) di_dmatrix = new DI_dmatrix();
1034
1035    static AW_window *viewer = 0;
1036    if (!viewer) viewer = DI_create_view_matrix_window(aww->get_root(), di_dmatrix, sparam);
1037
1038    di_dmatrix->init();
1039    di_dmatrix->display(false);
1040
1041    GLOBAL_MATRIX.set_changed_cb(refresh_matrix_display);
1042
1043    viewer->activate();
1044}
1045
1046static void di_save_matrix_cb(AW_window *aww, AW_CL cl_weightedFilter) {
1047    // save the matrix
1048    WeightedFilter *weighted_filter = (WeightedFilter*)cl_weightedFilter;
1049    GB_ERROR error = di_calculate_matrix(aww, weighted_filter, 0, true, NULL);
1050    if (!error) {
1051        char              *filename = aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_FILENAME)->read_string();
1052        enum DI_SAVE_TYPE  type     = (enum DI_SAVE_TYPE)aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_TYPE)->read_int();
1053
1054        GLOBAL_MATRIX.get()->save(filename, type);
1055        free(filename);
1056    }
1057    AW_refresh_fileselection(aww->get_root(), AWAR_DIST_SAVE_MATRIX_BASE);
1058    aww->hide_or_notify(error);
1059}
1060
1061AW_window *DI_create_save_matrix_window(AW_root *aw_root, save_matrix_params *save_params) {
1062    static AW_window_simple *aws = 0;
1063    if (!aws) {
1064        aws = new AW_window_simple;
1065        aws->init(aw_root, "SAVE_MATRIX", "Save Matrix");
1066        aws->load_xfig("sel_box_user.fig");
1067
1068        aws->at("close"); aws->callback((AW_CB0)AW_POPDOWN);
1069        aws->create_button("CLOSE", "CANCEL", "C");
1070
1071
1072        aws->at("help"); aws->callback(AW_POPUP_HELP, (AW_CL)"save_matrix.hlp");
1073        aws->create_button("HELP", "HELP", "H");
1074
1075        aws->at("user");
1076        aws->create_option_menu(AWAR_DIST_SAVE_MATRIX_TYPE);
1077        aws->insert_default_option("Phylip Format (Lower Triangular Matrix)", "P", DI_SAVE_PHYLIP_COMP);
1078        aws->insert_option("Readable (using NDS)", "R", DI_SAVE_READABLE);
1079        aws->insert_option("Tabbed (using NDS)", "R", DI_SAVE_TABBED);
1080        aws->update_option_menu();
1081
1082        AW_create_fileselection(aws, save_params->awar_base);
1083
1084        aws->at("save2");
1085        aws->callback(di_save_matrix_cb, (AW_CL)save_params->weighted_filter);
1086        aws->create_button("SAVE", "SAVE", "S");
1087
1088        aws->callback((AW_CB0)AW_POPDOWN);
1089        aws->at("cancel2");
1090        aws->create_button("CLOSE", "CANCEL", "C");
1091    }
1092    return aws;
1093}
1094
1095static AW_window *awt_create_select_cancel_window(AW_root *aw_root)
1096{
1097    AW_window_simple *aws = new AW_window_simple;
1098    aws->init(aw_root, "SELECT_CHARS_TO_CANCEL_COLUMN", "CANCEL SELECT");
1099    aws->load_xfig("di_cancel.fig");
1100
1101    aws->at("close"); aws->callback((AW_CB0)AW_POPDOWN);
1102    aws->create_button("CLOSE", "CLOSE", "C");
1103
1104    aws->at("cancel");
1105    aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1106
1107    return (AW_window *)aws;
1108}
1109
1110static const char *enum_trans_to_string[] = {
1111    "none",
1112    "similarity",
1113    "jukes_cantor",
1114    "felsenstein",
1115
1116    "pam",
1117    "hall",
1118    "barker",
1119    "chemical",
1120
1121    "haesch",
1122    "kimura",
1123    "olsen",
1124    "felsenstein voigt",
1125    "olsen voigt",
1126    "max ml"
1127};
1128
1129static void di_calculate_tree_cb(AW_window *aww, AW_CL cl_weightedFilter, AW_CL bootstrap_flag) {
1130    AW_root  *aw_root   = aww->get_root();
1131    GB_ERROR  error     = 0;
1132    GBT_TREE *tree      = 0;
1133    StrArray *all_names = 0;
1134
1135    int loop_count      = 0;
1136    int bootstrap_count = aw_root->awar(AWAR_DIST_BOOTSTRAP_COUNT)->read_int();
1137
1138    {
1139        char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1140        error           = GBT_check_tree_name(tree_name);
1141        free(tree_name);
1142    }
1143
1144    WeightedFilter *weighted_filter = (WeightedFilter*)cl_weightedFilter;
1145
1146    SmartPtr<arb_progress>  progress;
1147    SmartPtr<ConsensusTree> ctree;
1148
1149    if (!error) {
1150        if (bootstrap_flag) {
1151            if (bootstrap_count) {
1152                progress = new arb_progress("Calculating bootstrap trees", bootstrap_count);
1153            }
1154            else {
1155                progress = new arb_progress("Calculating bootstrap trees (KILL to stop)", INT_MAX);
1156            }
1157            progress->auto_subtitles("tree");
1158        }
1159        else {
1160            progress = new arb_progress("Calculating tree");
1161        }
1162
1163        if (bootstrap_flag) {
1164            GLOBAL_MATRIX.forget();
1165            GLOBAL_MATRIX.set_changed_cb(NULL); // otherwise matrix window will repeatedly pop up/down
1166
1167            error = di_calculate_matrix(aww, weighted_filter, bootstrap_flag, true, NULL);
1168            if (!error) {
1169                DI_MATRIX *matr = GLOBAL_MATRIX.get();
1170                if (!matr) {
1171                    error = "unexpected error in di_calculate_matrix_cb (data missing)";
1172                }
1173                else {
1174                    all_names = new StrArray;
1175                    all_names->reserve(matr->nentries+2);
1176
1177                    for (long i=0; i<matr->nentries; i++) {
1178                        all_names->put(strdup(matr->entries[i]->name));
1179                    }
1180                    ctree = new ConsensusTree(*all_names);
1181                }
1182            }
1183        }
1184    }
1185
1186    do {
1187        if (error) break;
1188
1189        bool aborted = false;
1190
1191        if (bootstrap_flag && loop_count) GLOBAL_MATRIX.forget(); // in first loop we already have a valid matrix
1192        error = di_calculate_matrix(aww, weighted_filter, bootstrap_flag, !bootstrap_flag, &aborted);
1193        if (error && aborted) {
1194            error = 0;          // clear error (otherwise no tree will be read below)
1195            break;              // end of bootstrap
1196        }
1197
1198        if (!GLOBAL_MATRIX.exists()) {
1199            error = "unexpected error in di_calculate_matrix_cb (data missing)";
1200            break;
1201        }
1202
1203        DI_MATRIX  *matr  = GLOBAL_MATRIX.get();
1204        char     **names = (char **)calloc(sizeof(char *), (size_t)matr->nentries+2);
1205        for (long i=0; i<matr->nentries; i++) {
1206            names[i] = matr->entries[i]->name;
1207        }
1208        tree = neighbourjoining(names, matr->matrix->m, matr->nentries, sizeof(GBT_TREE));
1209
1210        if (bootstrap_flag) {
1211            ctree->insert(tree, 1);
1212            GBT_delete_tree(tree); tree = 0;
1213            loop_count++;
1214            progress->inc();
1215            if (!bootstrap_count) { // when waiting for kill
1216                int        t     = time(0);
1217                static int tlast = 0;
1218
1219                if (t>tlast) {
1220                    progress->force_update();
1221                    tlast = t;
1222                }
1223            }
1224        }
1225        free(names);
1226    } while (bootstrap_flag && loop_count != bootstrap_count);
1227
1228    if (!error) {
1229        if (bootstrap_flag) {
1230            tree  = ctree->get_consensus_tree();
1231            error = GBT_is_invalid(tree);
1232            di_assert(!error);
1233        }
1234
1235        char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1236        GB_begin_transaction(GLOBAL_gb_main);
1237        error = GBT_write_tree(GLOBAL_gb_main, 0, tree_name, tree);
1238
1239        if (!error) {
1240            char       *filter_name = AWT_get_combined_filter_name(aw_root, "dist");
1241            int         transr      = aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1242            const char *comment     = GBS_global_string("PRG=dnadist CORR=%s FILTER=%s PKG=ARB", enum_trans_to_string[transr], filter_name);
1243
1244            error = GBT_write_tree_rem(GLOBAL_gb_main, tree_name, comment);
1245            free(filter_name);
1246        }
1247
1248        error = GB_end_transaction(GLOBAL_gb_main, error);
1249        free(tree_name);
1250    }
1251
1252    if (tree) GBT_delete_tree(tree);
1253
1254    // aw_status(); // remove 'abort' flag (@@@ got no equiv for arb_progress yet. really needed?)
1255
1256    if (bootstrap_flag) {
1257        if (all_names) delete all_names;
1258        GLOBAL_MATRIX.forget();
1259    }
1260#if defined(DEBUG)
1261    else {
1262        di_assert(all_names == 0);
1263    }
1264#endif // DEBUG
1265
1266    if (error) aw_message(error);
1267    progress->done();
1268}
1269
1270
1271static void di_autodetect_callback(AW_window *aww)
1272{
1273    GB_push_transaction(GLOBAL_gb_main);
1274
1275    GLOBAL_MATRIX.forget();
1276
1277    AW_root *aw_root = aww->get_root();
1278    char    *use     = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1279    long     ali_len = GBT_get_alignment_len(GLOBAL_gb_main, use);
1280
1281    if (ali_len<=0) {
1282        GB_pop_transaction(GLOBAL_gb_main);
1283        delete use;
1284        aw_message("Please select a valid alignment");
1285        return;
1286    }
1287
1288    arb_progress progress("Analyzing data");
1289
1290    char *filter_str = aw_root->awar(AWAR_DIST_FILTER_FILTER)->read_string();
1291    char *cancel     = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1292
1293    AliView *aliview;
1294    {
1295        AP_filter *ap_filter = NULL;
1296        long       flen  = strlen(filter_str);
1297
1298        if (flen == ali_len) {
1299            ap_filter = new AP_filter(filter_str, "0", ali_len);
1300        }
1301        else {
1302            if (flen) {
1303                aw_message("WARNING: YOUR FILTER LEN IS NOT THE ALIGNMENT LEN\nFILTER IS TRUNCATED WITH ZEROS OR CUTTED");
1304                ap_filter = new AP_filter(filter_str, "0", ali_len);
1305            }
1306            else {
1307                ap_filter = new AP_filter(ali_len); // unfiltered
1308            }
1309        }
1310
1311        AP_weights ap_weights(ap_filter);
1312
1313        aliview = new AliView(GLOBAL_gb_main, *ap_filter, ap_weights, use);
1314        delete ap_filter;
1315    }
1316
1317    {
1318        DI_MATRIX phm(*aliview, aw_root);
1319
1320        {
1321            char *load_what      = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
1322            char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
1323
1324            LoadWhat all_flag = (strcmp(load_what, "all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
1325
1326            GLOBAL_MATRIX.forget();
1327
1328            MatrixOrder order(GLOBAL_gb_main, sort_tree_name);
1329            phm.load(all_flag, order, true, NULL);
1330
1331            free(sort_tree_name);
1332            free(load_what);
1333        }
1334
1335        GB_pop_transaction(GLOBAL_gb_main);
1336
1337        progress.subtitle("Search Correction");
1338
1339        phm.analyse();
1340    }
1341
1342    free(cancel);
1343    delete aliview;
1344
1345    free(use);
1346    free(filter_str);
1347}
1348
1349__ATTR__NORETURN static void di_exit(AW_window *aww) {
1350    if (GLOBAL_gb_main) {
1351        aww->get_root()->unlink_awars_from_DB(GLOBAL_gb_main);
1352        GB_close(GLOBAL_gb_main);
1353    }
1354    exit(0);
1355}
1356
1357static void di_calculate_full_matrix_cb(AW_window *aww, AW_CL cl_weightedFilter) {
1358    GLOBAL_MATRIX.forget_if_not_has_type(DI_MATRIX_FULL);
1359
1360    WeightedFilter *weighted_filter = (WeightedFilter*)cl_weightedFilter;
1361    GB_ERROR        error           = di_calculate_matrix(aww, weighted_filter, 0, true, NULL);
1362
1363    aw_message_if(error);
1364}
1365
1366static void di_compress_tree_cb(AW_window *aww, AW_CL cl_weightedFilter) {
1367    GB_transaction  ta(GLOBAL_gb_main);
1368    char           *treename = aww->get_root()->awar(AWAR_DIST_TREE_COMP_NAME)->read_string();
1369    GB_ERROR        error    = 0;
1370    GBT_TREE       *tree     = GBT_read_tree(GLOBAL_gb_main, treename, sizeof(GBT_TREE));
1371   
1372    if (!tree) {
1373        error = GB_await_error();
1374    }
1375    else {
1376        GLOBAL_MATRIX.forget_if_not_has_type(DI_MATRIX_COMPRESSED);
1377
1378        WeightedFilter *weighted_filter = (WeightedFilter*)cl_weightedFilter;
1379        error = di_calculate_matrix(aww, weighted_filter, 0, true, NULL);
1380        if (!error && !GLOBAL_MATRIX.exists()) error = "Failed to calculate your matrix (bug?)";
1381        if (!error) error = GLOBAL_MATRIX.get()->compress(tree);
1382
1383        GBT_delete_tree(tree);
1384    }
1385    free(treename);
1386    aw_message_if(error);
1387}
1388
1389static void di_define_sort_tree_name_cb(AW_window *aww) {
1390    AW_root *aw_root = aww->get_root();
1391    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1392    aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->write_string(tree_name);
1393    free(tree_name);
1394}
1395static void di_define_compression_tree_name_cb(AW_window *aww) {
1396    AW_root *aw_root = aww->get_root();
1397    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1398    aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->write_string(tree_name);
1399    free(tree_name);
1400}
1401
1402static void di_define_save_tree_name_cb(AW_window *aww) {
1403    AW_root *aw_root = aww->get_root();
1404    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1405    aw_root->awar(AWAR_DIST_TREE_STD_NAME)->write_string(tree_name);
1406    free(tree_name);
1407}
1408
1409
1410AW_window *DI_create_matrix_window(AW_root *aw_root) {
1411    AW_window_simple_menu *aws = new AW_window_simple_menu;
1412    aws->init(aw_root, "NEIGHBOUR JOINING", "NEIGHBOUR JOINING [ARB_DIST]");
1413    aws->load_xfig("di_ge_ma.fig");
1414    aws->button_length(10);
1415
1416    aws->at("close");
1417    aws->callback(di_exit);
1418    aws->create_button("CLOSE", "CLOSE", "C");
1419
1420    aws->at("help");
1421    aws->callback(AW_POPUP_HELP, (AW_CL)"dist.hlp");
1422    aws->create_button("HELP", "HELP", "H");
1423
1424
1425    GB_push_transaction(GLOBAL_gb_main);
1426
1427#if defined(DEBUG)
1428    AWT_create_debug_menu(aws);
1429#endif // DEBUG
1430
1431    aws->create_menu("FILE", "F", AWM_ALL);
1432    aws->insert_menu_topic("macros", "Macros  ...", "M", "macro.hlp", AWM_ALL, (AW_CB)awt_popup_macro_window, (AW_CL)"NEIGHBOUR_JOINING", (AW_CL)GLOBAL_gb_main);
1433    aws->insert_menu_topic("quit",   "Quit",        "Q", "quit.hlp",  AWM_ALL, (AW_CB)di_exit,  0,  0);
1434
1435    aws->create_menu("Properties", "P", AWM_ALL);
1436    aws->insert_menu_topic("frame_props", "Frame ...",                  "F", "props_frame.hlp", AWM_ALL, AW_POPUP,(AW_CL)AW_preset_window,   0);
1437    aws->sep______________();
1438    AW_insert_common_property_menu_entries(aws);
1439    aws->sep______________();
1440    aws->insert_menu_topic("save_props",  "Save Properties (dist.arb)", "S", "savedef.hlp",     AWM_ALL,          (AW_CB)AW_save_properties, 0, 0);
1441
1442    aws->insert_help_topic("help ...", "h", "dist.hlp", AWM_ALL, (AW_CB)AW_POPUP_HELP, (AW_CL)"dist.hlp", 0);
1443
1444    // ------------------
1445    //      left side
1446
1447    aws->at("which_species");
1448    aws->create_option_menu(AWAR_DIST_WHICH_SPECIES);
1449    aws->insert_option("all", "a", "all");
1450    aws->insert_default_option("marked",   "m", "marked");
1451    aws->update_option_menu();
1452
1453    aws->at("which_alignment");
1454    awt_create_selection_list_on_alignments(GLOBAL_gb_main, (AW_window *)aws, AWAR_DIST_ALIGNMENT, "*=");
1455
1456    // filter & weights
1457
1458    AW_awar *awar_dist_alignment = aws->get_root()->awar_string(AWAR_DIST_ALIGNMENT);
1459    WeightedFilter *weighted_filter =               // do NOT free (bound to callbacks)
1460        new WeightedFilter(GLOBAL_gb_main, aws->get_root(), AWAR_DIST_FILTER_NAME, AWAR_DIST_COLUMN_STAT_NAME, awar_dist_alignment);
1461
1462    aws->at("filter_select");
1463    aws->callback(AW_POPUP, (AW_CL)awt_create_select_filter_win, (AW_CL)(weighted_filter->get_adfiltercbstruct()));
1464    aws->create_button("SELECT_FILTER", AWAR_DIST_FILTER_NAME);
1465
1466    aws->at("weights_select");
1467    aws->sens_mask(AWM_EXP);
1468    aws->callback(AW_POPUP, (AW_CL)COLSTAT_create_selection_window, (AW_CL)weighted_filter->get_column_stat());
1469    aws->create_button("SELECT_COL_STAT", AWAR_DIST_COLUMN_STAT_NAME);
1470    aws->sens_mask(AWM_ALL);
1471
1472    aws->at("which_cancel");
1473    aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1474
1475    aws->at("cancel_select");
1476    aws->callback((AW_CB1)AW_POPUP, (AW_CL)awt_create_select_cancel_window);
1477    aws->create_button("SELECT_CANCEL_CHARS", "Info", "C");
1478
1479    aws->at("change_matrix");
1480    aws->callback(AW_POPUP, (AW_CL)create_dna_matrix_window, 0);
1481    aws->create_button("EDIT_MATRIX", "Edit Matrix");
1482
1483    aws->at("enable");
1484    aws->create_toggle(AWAR_DIST_MATRIX_DNA_ENABLED);
1485
1486    aws->at("which_correction");
1487    aws->create_option_menu(AWAR_DIST_CORR_TRANS, NULL, "");
1488    aws->insert_option("none",                    "n", (int)DI_TRANSFORMATION_NONE);
1489    aws->insert_option("similarity",              "n", (int)DI_TRANSFORMATION_SIMILARITY);
1490    aws->insert_option("jukes-cantor (dna)",      "c", (int)DI_TRANSFORMATION_JUKES_CANTOR);
1491    aws->insert_option("felsenstein (dna)",       "f", (int)DI_TRANSFORMATION_FELSENSTEIN);
1492    aws->insert_option("olsen (dna)",             "o", (int)DI_TRANSFORMATION_OLSEN);
1493    aws->insert_option("felsenstein/voigt (exp)", "1", (int)DI_TRANSFORMATION_FELSENSTEIN_VOIGT);
1494    aws->insert_option("olsen/voigt (exp)",       "2", (int)DI_TRANSFORMATION_OLSEN_VOIGT);
1495    aws->insert_option("haesch (exp)",            "f", (int)DI_TRANSFORMATION_HAESCH);
1496    aws->insert_option("kimura (pro)",            "k", (int)DI_TRANSFORMATION_KIMURA);
1497    aws->insert_option("PAM (protein)",           "c", (int)DI_TRANSFORMATION_PAM);
1498    aws->insert_option("Cat. Hall(exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_HALL);
1499    aws->insert_option("Cat. Barker(exp)",        "c", (int)DI_TRANSFORMATION_CATEGORIES_BARKER);
1500    aws->insert_option("Cat.Chem (exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_CHEMICAL);
1501    aws->insert_default_option("unknown", "u", (int)DI_TRANSFORMATION_NONE);
1502
1503    aws->update_option_menu();
1504
1505    aws->at("autodetect");   // auto
1506    aws->callback(di_autodetect_callback);
1507    aws->sens_mask(AWM_EXP);
1508    aws->create_button("AUTODETECT_CORRECTION", "AUTODETECT", "A");
1509    aws->sens_mask(AWM_ALL);
1510
1511    // -------------------
1512    //      right side
1513
1514
1515    aws->at("mark_distance");
1516    aws->callback(di_mark_by_distance, (AW_CL)weighted_filter);
1517    aws->create_autosize_button("MARK_BY_DIST", "Mark all species");
1518
1519    aws->at("mark_lower");
1520    aws->create_input_field(AWAR_DIST_MIN_DIST, 5);
1521
1522    aws->at("mark_upper");
1523    aws->create_input_field(AWAR_DIST_MAX_DIST, 5);
1524
1525    // -----------------
1526
1527    // tree selection
1528
1529    aws->at("tree_list");
1530    awt_create_selection_list_on_trees(GLOBAL_gb_main, (AW_window *)aws, AWAR_DIST_TREE_CURR_NAME);
1531
1532    aws->at("detect_clusters");
1533    aws->callback(AW_POPUP, (AW_CL)DI_create_cluster_detection_window, (AW_CL)weighted_filter);
1534    aws->create_autosize_button("DETECT_CLUSTERS", "Detect homogenous clusters in tree", "D");
1535
1536    // matrix calculation
1537
1538    aws->button_length(18);
1539    aws->at("calculate");
1540    aws->callback(di_calculate_full_matrix_cb, (AW_CL)weighted_filter);
1541    aws->create_button("CALC_FULL_MATRIX", "Calculate\nFull Matrix", "F");
1542    aws->at("compress");
1543    aws->callback(di_compress_tree_cb, (AW_CL)weighted_filter);
1544    aws->create_button("CALC_COMPRESSED_MATRIX", "Calculate\nCompressed Matrix", "C");
1545
1546    aws->button_length(13);
1547
1548    {
1549        save_matrix_params *sparams = new save_matrix_params; // do not free (bound to callbacks)
1550
1551        sparams->awar_base       = AWAR_DIST_SAVE_MATRIX_BASE;
1552        sparams->weighted_filter = weighted_filter;
1553
1554        aws->at("save_matrix");
1555        aws->callback(AW_POPUP, (AW_CL)DI_create_save_matrix_window, (AW_CL)sparams);
1556        aws->create_button("SAVE_MATRIX", "Save matrix", "M");
1557
1558        aws->at("view_matrix");
1559        aws->callback(di_view_matrix_cb, (AW_CL)sparams);
1560        aws->create_button("VIEW_MATRIX", "View matrix", "V");
1561    }
1562
1563    aws->button_length(22);
1564    aws->at("use_compr_tree");
1565    aws->callback(di_define_compression_tree_name_cb);
1566    aws->create_button("USE_COMPRESSION_TREE", "Use to compress", "");
1567    aws->at("use_sort_tree");
1568    aws->callback(di_define_sort_tree_name_cb);
1569    aws->create_button("USE_SORT_TREE", "Use to sort", "");
1570
1571    aws->at("compr_tree_name");
1572    aws->create_input_field(AWAR_DIST_TREE_COMP_NAME, 12);
1573    aws->at("sort_tree_name");
1574    aws->create_input_field(AWAR_DIST_TREE_SORT_NAME, 12);
1575
1576    // tree calculation
1577
1578    aws->button_length(18);
1579    aws->at("t_calculate");
1580    aws->callback(di_calculate_tree_cb, (AW_CL)weighted_filter, 0);
1581    aws->create_button("CALC_TREE", "Calculate \ntree", "C");
1582    aws->at("bootstrap");
1583    aws->callback(di_calculate_tree_cb, (AW_CL)weighted_filter, 1);
1584    aws->create_button("CALC_BOOTSTRAP_TREE", "Calculate \nbootstrap tree");
1585
1586    aws->button_length(22);
1587    aws->at("use_existing");
1588    aws->callback(di_define_save_tree_name_cb);
1589    aws->create_button("USE_NAME", "Use as new tree name", "");
1590
1591    aws->at("calc_tree_name");
1592    aws->create_input_field(AWAR_DIST_TREE_STD_NAME, 12);
1593
1594    aws->at("bcount");
1595    aws->create_input_field(AWAR_DIST_BOOTSTRAP_COUNT, 7);
1596
1597    GB_pop_transaction(GLOBAL_gb_main);
1598    return aws;
1599}
Note: See TracBrowser for help on using the browser.