source: branches/profile/SEQ_QUALITY/SQ_functions.cxx

Last change on this file was 10113, checked in by westram, 13 years ago
  • changed dynamic value-initializations into default-initializations
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 38.9 KB
Line 
1//  ==================================================================== //
2//                                                                       //
3//    File      : SQ_functions.cxx                                       //
4//    Purpose   : Implementation of SQ_functions.h                       //
5//                                                                       //
6//                                                                       //
7//  Coded by Juergen Huber in July 2003 - February 2004                  //
8//  Coded by Kai Bader (baderk@in.tum.de) in 2007 - 2008                 //
9//  Copyright Department of Microbiology (Technical University Munich)   //
10//                                                                       //
11//  Visit our web site at: http://www.arb-home.de/                       //
12//                                                                       //
13//  ==================================================================== //
14
15#include "SQ_ambiguities.h"
16#include "SQ_helix.h"
17#include "SQ_physical_layout.h"
18#include "SQ_functions.h"
19
20#include <aw_preset.hxx>
21#include <arb_progress.h>
22#include <arbdbt.h>
23
24using namespace std;
25
26typedef GBDATA *(*species_iterator)(GBDATA *);
27
28static SQ_GroupDataDictionary group_dict;
29
30enum {
31    CS_CLEAR, CS_PASS1
32};
33
34void SQ_clear_group_dictionary() {
35    SQ_GroupDataDictionary tmp;
36    swap(tmp, group_dict);
37}
38
39static GB_ERROR no_data_error(GBDATA * gb_species, const char *ali_name) {
40    GBDATA *gb_name = GB_entry(gb_species, "name");
41    const char *name = "<unknown>";
42    if (gb_name)
43        name = GB_read_char_pntr(gb_name);
44    return GBS_global_string("Species '%s' has no data in alignment '%s'",
45            name, ali_name);
46}
47
48static int sq_round(double value) {
49    int x;
50
51    value += 0.5;
52    x = (int) floor(value);
53    return x;
54}
55
56GB_ERROR SQ_remove_quality_entries(GBDATA *gb_main) {
57    GB_push_transaction(gb_main);
58    GB_ERROR error = NULL;
59
60    for (GBDATA *gb_species = GBT_first_species(gb_main); gb_species && !error; gb_species
61    = GBT_next_species(gb_species)) {
62        GBDATA *gb_quality = GB_search(gb_species, "quality", GB_FIND);
63        if (gb_quality)
64            error = GB_delete(gb_quality);
65    }
66    if (error)
67        GB_abort_transaction(gb_main);
68    else
69        GB_pop_transaction(gb_main);
70    return error;
71}
72
73GB_ERROR SQ_evaluate(GBDATA * gb_main, const SQ_weights & weights, bool marked_only) {
74    GB_push_transaction(gb_main);
75    char *alignment_name = GBT_get_default_alignment(gb_main);
76    seq_assert(alignment_name);
77
78    species_iterator getFirst = marked_only ? GBT_first_marked_species : GBT_first_species;
79    species_iterator getNext  = marked_only ? GBT_next_marked_species  : GBT_next_species;
80
81    GB_ERROR error = 0;
82    for (GBDATA *gb_species = getFirst(gb_main); gb_species && !error; gb_species = getNext(gb_species)) {
83        GBDATA *gb_name = GB_entry(gb_species, "name");
84        if (!gb_name) {
85            error = GB_get_error();
86        }
87        else {
88            GBDATA *gb_quality = GB_entry(gb_species, "quality");
89            if (gb_quality) {
90                GBDATA *gb_quality_ali = GB_entry(gb_quality, alignment_name);
91
92                if (!gb_quality_ali) {
93                    error = GBS_global_string("No alignment entry '%s' in quality data", alignment_name);
94                }
95                else {
96
97                    // evaluate the percentage of bases the actual sequence consists of
98                    GBDATA *gb_result1 = GB_search(gb_quality_ali, "percent_of_bases", GB_INT);
99                    int     bases      = GB_read_int(gb_result1);
100                   
101                    double result = bases<4 ? 0 : (bases<6 ? 1 : 2);
102
103                    if (result != 0) result = (result * weights.bases) / 2;
104
105                    double value = 0;
106                    value += result;
107
108                    // evaluate the difference in number of bases from sequence to group
109                    GBDATA *gb_result2 = GB_search(gb_quality_ali, "percent_base_deviation", GB_INT);
110                    int     dfa        = GB_read_int(gb_result2);
111                    if (abs(dfa) < 2) {
112                        result = 5;
113                    }
114                    else {
115                        if (abs(dfa) < 4) {
116                            result = 4;
117                        }
118                        else {
119                            if (abs(dfa) < 6) {
120                                result = 3;
121                            }
122                            else {
123                                if (abs(dfa) < 8) {
124                                    result = 2;
125                                }
126                                else {
127                                    if (abs(dfa) < 10) {
128                                        result = 1;
129                                    }
130                                    else {
131                                        result = 0;
132                                    }
133                                }
134                            }
135                        }
136                    }
137                    if (result != 0) result = (result * weights.diff_from_average) / 5;
138                    value += result;
139
140                    // evaluate the number of positions where no helix can be built
141                    GBDATA *gb_result3 = GB_search(gb_quality_ali, "number_of_no_helix", GB_INT);
142                    int     noh        = GB_read_int(gb_result3);
143                    if (noh < 20) {
144                        result = 5;
145                    }
146                    else {
147                        if (noh < 50) {
148                            result = 4;
149                        }
150                        else {
151                            if (noh < 125) {
152                                result = 3;
153                            }
154                            else {
155                                if (noh < 250) {
156                                    result = 2;
157                                }
158                                else {
159                                    if (noh < 500) {
160                                        result = 1;
161                                    }
162                                    else {
163                                        result = 0;
164                                    }
165                                }
166                            }
167                        }
168                    }
169                    if (result != 0) result = (result * weights.helix) / 5;
170                    value += result;
171
172                    // evaluate the consensus
173                    GBDATA *gb_result4 = GB_search(gb_quality_ali, "consensus_evaluated", GB_INT);
174                    int     cos        = GB_read_int(gb_result4);
175
176                    result                   = cos;
177                    if (result != 0) result  = (result * weights.consensus) / 12;
178                    value                   += result;
179
180                    // evaluate the number of iupacs in a sequence
181                    GBDATA *gb_result5 = GB_search(gb_quality_ali, "iupac_value", GB_INT);
182                    int     iupv       = GB_read_int(gb_result5);
183                   
184                    if (iupv < 1) {
185                        result = 3;
186                    }
187                    else {
188                        if (iupv < 5) {
189                            result = 2;
190                        }
191                        else {
192                            if (iupv < 10) {
193                                result = 1;
194                            }
195                            else {
196                                result = 0;
197                            }
198                        }
199                    }
200                    if (result != 0) result = (result * weights.iupac) / 3;
201                    value += result;
202
203                    // evaluate the difference in the GC proportion from sequence to group
204                    GBDATA *gb_result6 = GB_search(gb_quality_ali, "percent_GC_difference", GB_INT);
205                    int     gcprop     = GB_read_int(gb_result6);
206                   
207                    if (abs(gcprop) < 1) {
208                        result = 5;
209                    }
210                    else {
211                        if (abs(gcprop) < 2) {
212                            result = 4;
213                        }
214                        else {
215                            if (abs(gcprop) < 4) {
216                                result = 3;
217                            }
218                            else {
219                                if (abs(gcprop) < 8) {
220                                    result = 2;
221                                }
222                                else {
223                                    if (abs(gcprop) < 16)
224                                        result = 1;
225                                    else {
226                                        result = 0;
227                                    }
228                                }
229                            }
230                        }
231                    }
232                    if (result != 0) result = (result * weights.gc) / 5;
233                    value += result;
234
235                    // write the final value of the evaluation
236                    int     value2     = sq_round(value);
237                    GBDATA *gb_result7 = GB_search(gb_quality_ali, "evaluation", GB_INT);
238                    seq_assert(gb_result7);
239                    GB_write_int(gb_result7, value2);
240                }
241            }
242        }
243    }
244    free(alignment_name);
245
246    error = GB_end_transaction(gb_main, error);
247    return error;
248}
249
250static char *SQ_fetch_filtered_sequence(GBDATA * read_sequence, AP_filter * filter) {
251    char *filteredSequence = 0;
252    if (read_sequence) {
253        const char   *rawSequence        = GB_read_char_pntr(read_sequence);
254        int           filteredLength     = filter->get_filtered_length();
255        const size_t *filterpos_2_seqpos = filter->get_filterpos_2_seqpos();
256
257        filteredSequence = (char*)malloc(filteredLength * sizeof(char));
258        if (filteredSequence) {
259            for (int i = 0; i < filteredLength; ++i) {
260                filteredSequence[i] = rawSequence[filterpos_2_seqpos[i]];
261            }
262        }
263    }
264    return filteredSequence;
265}
266
267static GB_ERROR SQ_pass1(SQ_GroupData * globalData, GBDATA * gb_main, GBT_TREE * node, AP_filter * filter) {
268    GB_push_transaction(gb_main);
269
270    GB_ERROR  error          = 0;
271    char     *alignment_name = GBT_get_default_alignment(gb_main); seq_assert(alignment_name);
272    GBDATA   *gb_species     = node->gb_node;
273    GBDATA   *gb_name        = GB_entry(gb_species, "name");
274
275    if (!gb_name) error = GB_get_error();
276    else {
277        GBDATA *gb_ali = GB_entry(gb_species, alignment_name);
278
279        if (!gb_ali) {
280            error = no_data_error(gb_species, alignment_name);
281        }
282        else {
283            GBDATA *gb_quality = GB_search(gb_species, "quality", GB_CREATE_CONTAINER);
284
285            if (!gb_quality) {
286                error = GB_get_error();
287            }
288
289            GBDATA *read_sequence  = GB_entry(gb_ali, "data");
290            GBDATA *gb_quality_ali = GB_search(gb_quality, alignment_name, GB_CREATE_CONTAINER);
291           
292            if (!gb_quality_ali) error = GB_get_error();
293
294            // real calculations start here
295            if (read_sequence) {
296                char *rawSequence = SQ_fetch_filtered_sequence(read_sequence, filter);
297                int sequenceLength = filter->get_filtered_length();
298
299                // calculate physical layout of sequence
300                {
301                    SQ_physical_layout ps_chan;
302                    ps_chan.SQ_calc_physical_layout(rawSequence, sequenceLength, gb_quality_ali);
303
304                    // calculate the average number of bases in group
305                    globalData->SQ_count_sequences();
306                    globalData->SQ_set_avg_bases(ps_chan.
307                    SQ_get_number_of_bases());
308                    globalData->SQ_set_avg_gc(ps_chan.
309                    SQ_get_gc_proportion());
310                }
311
312                // get values for ambiguities
313                {
314                    SQ_ambiguities ambi_chan;
315                    ambi_chan.SQ_count_ambiguities(rawSequence, sequenceLength, gb_quality_ali);
316                }
317
318                // calculate the number of strong, weak and no helixes
319                {
320                    SQ_helix heli_chan(sequenceLength);
321                    heli_chan.SQ_calc_helix_layout(rawSequence, gb_main, alignment_name, gb_quality_ali, filter);
322                }
323
324                // calculate consensus sequence
325                {
326                    if (!globalData->SQ_is_initialized()) {
327                        globalData->SQ_init_consensus(sequenceLength);
328                    }
329                    globalData->SQ_add_sequence(rawSequence);
330                }
331
332                free(rawSequence);
333            }
334        }
335    }
336
337    free(alignment_name);
338
339    if (error)
340        GB_abort_transaction(gb_main);
341    else
342        GB_pop_transaction(gb_main);
343
344    return error;
345}
346
347GB_ERROR SQ_pass1_no_tree(SQ_GroupData * globalData, GBDATA * gb_main, AP_filter * filter, arb_progress& progress) {
348    GBDATA *read_sequence = 0;
349
350    GB_push_transaction(gb_main);
351
352    char     *alignment_name = GBT_get_default_alignment(gb_main);
353    GB_ERROR  error          = 0;
354
355    seq_assert(alignment_name);
356
357    // first pass operations
358    for (GBDATA *gb_species = GBT_first_species(gb_main); gb_species && !error; gb_species = GBT_next_species(gb_species)) {
359        GBDATA *gb_name = GB_entry(gb_species, "name");
360
361        if (!gb_name) {
362            error = GB_get_error();
363        }
364        else {
365            GBDATA *gb_ali = GB_entry(gb_species, alignment_name);
366
367            if (!gb_ali) {
368                error = no_data_error(gb_species, alignment_name);
369            }
370            else {
371                GBDATA *gb_quality = GB_search(gb_species, "quality", GB_CREATE_CONTAINER);
372                if (!gb_quality) {
373                    error = GB_get_error();
374                }
375
376                read_sequence = GB_entry(gb_ali, "data");
377
378                GBDATA *gb_quality_ali = GB_search(gb_quality, alignment_name, GB_CREATE_CONTAINER);
379                if (!gb_quality_ali)
380                    error = GB_get_error();
381
382                // real calculations start here
383                if (read_sequence) {
384                    char *rawSequence = SQ_fetch_filtered_sequence(read_sequence, filter);
385                    int sequenceLength = filter->get_filtered_length();
386
387                    // calculate physical layout of sequence
388                    SQ_physical_layout *ps_chan = new SQ_physical_layout;
389                    ps_chan->SQ_calc_physical_layout(rawSequence, sequenceLength, gb_quality_ali);
390
391                    // calculate the average number of bases in group
392                    globalData->SQ_count_sequences();
393                    globalData->SQ_set_avg_bases(ps_chan->SQ_get_number_of_bases());
394                    globalData->SQ_set_avg_gc(ps_chan->SQ_get_gc_proportion());
395                    delete ps_chan;
396
397                    // get values for ambiguities
398                    SQ_ambiguities *ambi_chan = new SQ_ambiguities;
399                    ambi_chan->SQ_count_ambiguities(rawSequence, sequenceLength, gb_quality_ali);
400                    delete ambi_chan;
401
402                    // calculate the number of strong, weak and no helixes
403                    SQ_helix *heli_chan = new SQ_helix(sequenceLength);
404                    heli_chan->SQ_calc_helix_layout(rawSequence, gb_main, alignment_name, gb_quality_ali, filter);
405                    delete heli_chan;
406
407                    // calculate consensus sequence
408                    {
409                        if (!globalData->SQ_is_initialized()) {
410                            globalData->SQ_init_consensus(sequenceLength);
411                        }
412                        globalData->SQ_add_sequence(rawSequence);
413                    }
414                    delete(rawSequence);
415                }
416            }
417        }
418        progress.inc_and_check_user_abort(error);
419    }
420
421    free(alignment_name);
422
423    if (error)
424        GB_abort_transaction(gb_main);
425    else
426        GB_pop_transaction(gb_main);
427
428    return error;
429}
430
431static GB_ERROR SQ_pass2(const SQ_GroupData * globalData, GBDATA * gb_main, GBT_TREE * node, AP_filter * filter) {
432    GB_push_transaction(gb_main);
433
434    char *alignment_name = GBT_get_default_alignment(gb_main);
435    seq_assert(alignment_name);
436
437    GBDATA   *gb_species = node->gb_node;
438    GBDATA   *gb_name    = GB_entry(gb_species, "name");
439    GB_ERROR  error      = 0;
440
441    if (!gb_name) error = GB_get_error();
442    else {
443        GBDATA *gb_ali = GB_entry(gb_species, alignment_name);
444
445        if (!gb_ali) {
446            error = no_data_error(gb_species, alignment_name);
447        }
448        else {
449            GBDATA *gb_quality = GB_search(gb_species, "quality", GB_CREATE_CONTAINER);
450            if (!gb_quality) error = GB_get_error();
451
452            GBDATA *gb_quality_ali = GB_search(gb_quality, alignment_name, GB_CREATE_CONTAINER);
453            if (!gb_quality_ali) error = GB_get_error();
454
455            GBDATA *read_sequence = GB_entry(gb_ali, "data");
456
457            // real calculations start here
458            if (read_sequence) {
459                char *rawSequence = SQ_fetch_filtered_sequence(read_sequence, filter);
460
461                /*
462                  calculate the average number of bases in group, and the difference of
463                  a single sequence in group from it
464                */
465                {
466                    GBDATA *gb_result1   = GB_search(gb_quality_ali, "number_of_bases", GB_INT);
467                    int     bases        = GB_read_int(gb_result1);
468                    int     avg_bases    = globalData->SQ_get_avg_bases();
469                    int     diff_percent = 0;
470
471                    if (avg_bases != 0) {
472                        double diff  = bases - avg_bases;
473                        diff         = (100 * diff) / avg_bases;
474                        diff_percent = sq_round(diff);
475                    }
476
477                    GBDATA *gb_result2 = GB_search(gb_quality_ali, "percent_base_deviation", GB_INT);
478                    seq_assert(gb_result2);
479                    GB_write_int(gb_result2, diff_percent);
480                }
481
482                /*
483                  calculate the average gc proportion in group, and the difference of
484                  a single sequence in group from it
485                */
486                {
487                    GBDATA *gb_result6   = GB_search(gb_quality_ali, "GC_proportion", GB_FLOAT);
488                    double  gcp          = GB_read_float(gb_result6);
489                    double  avg_gc       = globalData->SQ_get_avg_gc();
490                    int     diff_percent = 0;
491
492                    if (avg_gc != 0) {
493                        double diff  = gcp - avg_gc;
494                        diff         = (100 * diff) / avg_gc;
495                        diff_percent = sq_round(diff);
496                    }
497
498                    GBDATA *gb_result7 = GB_search(gb_quality_ali, "percent_GC_difference", GB_INT);
499                    seq_assert(gb_result7);
500                    GB_write_int(gb_result7, diff_percent);
501                }
502
503                /*
504                  get groupnames of visited groups
505                  search for name in group dictionary
506                  evaluate sequence with group consensus
507                */
508                GBDATA *gb_con     = GB_search(gb_quality_ali, "consensus_conformity", GB_CREATE_CONTAINER);
509                if (!gb_con) error = GB_get_error();
510
511                GBDATA *gb_dev     = GB_search(gb_quality_ali, "consensus_deviation", GB_CREATE_CONTAINER);
512                if (!gb_dev) error = GB_get_error();
513
514                GBT_TREE *backup       = node; // needed?
515                int       whilecounter = 0;
516                double    eval         = 0;
517
518                while (backup->father) {
519                    if (backup->name) {
520                        SQ_GroupDataDictionary::iterator GDI = group_dict.find(backup->name);
521                        if (GDI != group_dict.end()) {
522                            SQ_GroupDataPtr GD_ptr = GDI->second;
523
524                            consensus_result cr = GD_ptr->SQ_calc_consensus(rawSequence);
525
526                            double value1 = cr.conformity;
527                            double value2 = cr.deviation;
528                            int    value3 = GD_ptr->SQ_get_nr_sequences();
529
530                            GBDATA *gb_node_entry = GB_search(gb_con, "name", GB_STRING);
531                            seq_assert(gb_node_entry);
532                            GB_write_string(gb_node_entry, backup->name);
533
534                            gb_node_entry = GB_search(gb_con, "value", GB_FLOAT); seq_assert(gb_node_entry);
535                            GB_write_float(gb_node_entry, value1);
536
537                            gb_node_entry = GB_search(gb_con, "num_species", GB_INT); seq_assert(gb_node_entry);
538                            GB_write_int(gb_node_entry, value3);
539
540                            gb_node_entry = GB_search(gb_dev, "name", GB_STRING); seq_assert(gb_node_entry);
541                            GB_write_string(gb_node_entry, backup->name);
542
543                            gb_node_entry = GB_search(gb_dev, "value", GB_FLOAT); seq_assert(gb_node_entry);
544                            GB_write_float(gb_node_entry, value2);
545
546                            gb_node_entry = GB_search(gb_dev, "num_species", GB_INT); seq_assert(gb_node_entry);
547                            GB_write_int(gb_node_entry, value3);
548
549                            // if you parse the upper two values in the evaluate() function cut the following out
550                            // for time reasons i do the evaluation here, as i still have the upper two values
551                            // -------------cut this-----------------
552                            if (value1 > 0.95) {
553                                eval += 5;
554                            }
555                            else {
556                                if (value1 > 0.8) {
557                                    eval += 4;
558                                }
559                                else {
560                                    if (value1 > 0.6) {
561                                        eval += 3;
562                                    }
563                                    else {
564                                        if (value1 > 0.4) {
565                                            eval += 2;
566                                        }
567                                        else {
568                                            if (value1 > 0.25) {
569                                                eval += 1;
570                                            }
571                                            else {
572                                                eval += 0;
573                                            }
574                                        }
575                                    }
576                                }
577                            }
578                            if (value2 > 0.6) {
579                                eval += 0;
580                            }
581                            else {
582                                if (value2 > 0.4) {
583                                    eval += 1;
584                                }
585                                else {
586                                    if (value2 > 0.2) {
587                                        eval += 2;
588                                    }
589                                    else {
590                                        if (value2 > 0.1) {
591                                            eval += 3;
592                                        }
593                                        else {
594                                            if (value2 > 0.05) {
595                                                eval += 4;
596                                            }
597                                            else {
598                                                if (value2 > 0.025) {
599                                                    eval += 5;
600                                                }
601                                                else {
602                                                    if (value2 > 0.01) {
603                                                        eval += 6;
604                                                    }
605                                                    else {
606                                                        eval += 7;
607                                                    }
608                                                }
609                                            }
610                                        }
611                                    }
612                                }
613                            }
614                            whilecounter++;
615                            // ---------to this and scroll down--------
616                        }
617                    }
618                    backup = backup->father;
619                }
620
621                // --------also cut this------
622                int evaluation = 0;
623                if (eval != 0) {
624                    eval = eval / whilecounter;
625                    evaluation = sq_round(eval);
626                }
627                GBDATA *gb_result5 = GB_search(gb_quality_ali, "consensus_evaluated", GB_INT);
628                seq_assert(gb_result5);
629                GB_write_int(gb_result5, evaluation);
630                // --------end cut this-------
631
632                free(rawSequence);
633            }
634        }
635    }
636
637    free(alignment_name);
638
639    if (error)
640        GB_abort_transaction(gb_main);
641    else
642        GB_pop_transaction(gb_main);
643
644    return error;
645}
646
647GB_ERROR SQ_pass2_no_tree(const SQ_GroupData * globalData, GBDATA * gb_main, AP_filter * filter, arb_progress& progress) {
648    GBDATA *read_sequence = 0;
649
650    GB_push_transaction(gb_main);
651
652    char *alignment_name = GBT_get_default_alignment(gb_main);
653    seq_assert(alignment_name);
654
655    // second pass operations
656    GB_ERROR error = 0;
657    for (GBDATA *gb_species = GBT_first_species(gb_main); gb_species && !error; gb_species = GBT_next_species(gb_species)) {
658        GBDATA *gb_name = GB_entry(gb_species, "name");
659
660        if (!gb_name) {
661            error = GB_get_error();
662        }
663        else {
664            GBDATA *gb_ali = GB_entry(gb_species, alignment_name);
665            if (!gb_ali) {
666                error = no_data_error(gb_species, alignment_name);
667            }
668            else {
669                GBDATA *gb_quality = GB_search(gb_species, "quality", GB_CREATE_CONTAINER);
670                if (!gb_quality) error = GB_get_error();
671
672                GBDATA *gb_quality_ali = GB_search(gb_quality, alignment_name, GB_CREATE_CONTAINER);
673                if (!gb_quality_ali) error = GB_get_error();
674
675                read_sequence = GB_entry(gb_ali, "data");
676
677                // real calculations start here
678                if (read_sequence) {
679                    const char *rawSequence = SQ_fetch_filtered_sequence(read_sequence, filter);
680
681                    /*
682                      calculate the average number of bases in group, and the difference of
683                      a single sequence in group from it
684                    */
685                    {
686                        GBDATA *gb_result1   = GB_search(gb_quality_ali, "number_of_bases", GB_INT);
687                        int     bases        = GB_read_int(gb_result1);
688                        int     avg_bases    = globalData->SQ_get_avg_bases();
689                        int     diff_percent = 0;
690
691                        if (avg_bases != 0) {
692                            double diff  = bases - avg_bases;
693                            diff         = (100 * diff) / avg_bases;
694                            diff_percent = sq_round(diff);
695                        }
696
697                        GBDATA *gb_result2 = GB_search(gb_quality_ali, "percent_base_deviation", GB_INT);
698                        seq_assert(gb_result2);
699                        GB_write_int(gb_result2, diff_percent);
700                    }
701
702                    /*
703                     calculate the average gc proportion in group, and the difference of
704                     a single sequence in group from it
705                    */
706                    {
707                        GBDATA *gb_result6   = GB_search(gb_quality_ali, "GC_proportion", GB_FLOAT);
708                        double  gcp          = GB_read_float(gb_result6);
709                        double  avg_gc       = globalData->SQ_get_avg_gc();
710                        int     diff_percent = 0;
711
712                        if (avg_gc != 0) {
713                            double diff  = gcp - avg_gc;
714                            diff         = (100 * diff) / avg_gc;
715                            diff_percent = sq_round(diff);
716                        }
717
718                        GBDATA *gb_result7 = GB_search(gb_quality_ali, "percent_GC_difference", GB_INT);
719                        seq_assert(gb_result7);
720                        GB_write_int(gb_result7, diff_percent);
721                    }
722                    /*
723                      get groupnames of visited groups
724                      search for name in group dictionary
725                      evaluate sequence with group consensus
726                    */
727                    GBDATA *gb_con     = GB_search(gb_quality_ali, "consensus_conformity", GB_CREATE_CONTAINER);
728                    if (!gb_con) error = GB_get_error();
729
730                    GBDATA *gb_dev     = GB_search(gb_quality_ali, "consensus_deviation", GB_CREATE_CONTAINER);
731                    if (!gb_dev) error = GB_get_error();
732
733                    consensus_result cr = globalData->SQ_calc_consensus(rawSequence);
734
735                    double value1 = cr.conformity;
736                    double value2 = cr.deviation;
737                    int    value3 = globalData->SQ_get_nr_sequences();
738
739                    GBDATA *gb_node_entry = GB_search(gb_con, "name", GB_STRING);
740                    seq_assert(gb_node_entry);
741                    GB_write_string(gb_node_entry, "one global consensus");
742
743                    gb_node_entry = GB_search(gb_con, "value", GB_FLOAT); seq_assert(gb_node_entry);
744                    GB_write_float(gb_node_entry, value1);
745
746                    gb_node_entry = GB_search(gb_con, "num_species", GB_INT); seq_assert(gb_node_entry);
747                    GB_write_int(gb_node_entry, value3);
748
749                    gb_node_entry = GB_search(gb_dev, "name", GB_STRING); seq_assert(gb_node_entry);
750                    GB_write_string(gb_node_entry, "one global consensus");
751
752                    gb_node_entry = GB_search(gb_dev, "value", GB_FLOAT); seq_assert(gb_node_entry);
753                    GB_write_float(gb_node_entry, value2);
754
755                    gb_node_entry = GB_search(gb_dev, "num_species", GB_INT); seq_assert(gb_node_entry);
756                    GB_write_int(gb_node_entry, value3);
757
758                    double eval = 0;
759                   
760                    // if you parse the upper two values in the evaluate() function cut the following out
761                    // for time reasons i do the evaluation here, as i still have the upper two values
762                    // -------------cut this-----------------
763                    if (value1 > 0.95) {
764                        eval += 5;
765                    }
766                    else {
767                        if (value1 > 0.8) {
768                            eval += 4;
769                        }
770                        else {
771                            if (value1 > 0.6) {
772                                eval += 3;
773                            }
774                            else {
775                                if (value1 > 0.4) {
776                                    eval += 2;
777                                }
778                                else {
779                                    if (value1 > 0.25) {
780                                        eval += 1;
781                                    }
782                                    else {
783                                        eval += 0;
784                                    }
785                                }
786                            }
787                        }
788                    }
789                    if (value2 > 0.6) {
790                        eval += 0;
791                    }
792                    else {
793                        if (value2 > 0.4) {
794                            eval += 1;
795                        }
796                        else {
797                            if (value2 > 0.2) {
798                                eval += 2;
799                            }
800                            else {
801                                if (value2 > 0.1) {
802                                    eval += 3;
803                                }
804                                else {
805                                    if (value2 > 0.05) {
806                                        eval += 4;
807                                    }
808                                    else {
809                                        if (value2 > 0.025) {
810                                            eval += 5;
811                                        }
812                                        else {
813                                            if (value2 > 0.01) {
814                                                eval += 6;
815                                            }
816                                            else {
817                                                eval += 7;
818                                            }
819                                        }
820                                    }
821                                }
822                            }
823                        }
824                    }
825
826                    {
827                        int evaluation            = 0;
828                        if (eval != 0) evaluation = sq_round(eval);
829
830                        GBDATA *gb_result5 = GB_search(gb_quality_ali, "consensus_evaluated", GB_INT);
831                        seq_assert(gb_result5);
832                        GB_write_int(gb_result5, evaluation);
833                    }
834                    // --------end cut this-------
835                    delete(rawSequence);
836                }
837            }
838        }
839        progress.inc_and_check_user_abort(error);
840    }
841    free(alignment_name);
842
843    if (error)
844        GB_abort_transaction(gb_main);
845    else
846        GB_pop_transaction(gb_main);
847
848    return error;
849}
850
851int SQ_count_nodes(GBT_TREE *node) {
852    // calculate number of nodes in tree
853    return GBT_count_leafs(node)*2-1;
854}
855
856static void create_multi_level_consensus(GBT_TREE * node, SQ_GroupData * data) {
857    SQ_GroupData *newData = data->clone(); // save actual consensus
858    *newData = *data;
859    group_dict[node->name] = newData; // and link it with an name
860}
861
862void SQ_calc_and_apply_group_data(GBT_TREE * node, GBDATA * gb_main, SQ_GroupData * data, AP_filter * filter, arb_progress& progress) {
863    if (node->is_leaf) {
864        if (node->gb_node) {
865            SQ_pass1(data, gb_main, node, filter);
866            seq_assert(data->getSize()> 0);
867        }
868    }
869    else {
870        GBT_TREE *node1 = node->leftson;
871        GBT_TREE *node2 = node->rightson;
872
873        if (node->name) {
874            SQ_GroupData *leftData      = NULL;
875            bool          parentIsEmpty = false;
876
877            if (data->getSize() == 0) {
878                parentIsEmpty = true;
879                SQ_calc_and_apply_group_data(node1, gb_main, data, filter, progress); // process left branch with empty data
880                seq_assert(data->getSize()> 0);
881            }
882            else {
883                leftData = data->clone(); // create new empty SQ_GroupData
884                SQ_calc_and_apply_group_data(node1, gb_main, leftData, filter, progress); // process left branch
885                seq_assert(leftData->getSize()> 0);
886            }
887
888            SQ_GroupData *rightData = data->clone(); // create new empty SQ_GroupData
889            SQ_calc_and_apply_group_data(node2, gb_main, rightData, filter, progress); // process right branch
890            seq_assert(rightData->getSize()> 0);
891
892            if (!parentIsEmpty) {
893                data->SQ_add(*leftData);
894                delete leftData;
895            }
896
897            data->SQ_add(*rightData);
898            delete rightData;
899
900            create_multi_level_consensus(node, data);
901        }
902        else {
903            SQ_calc_and_apply_group_data(node1, gb_main, data, filter, progress); // enter left branch
904            seq_assert(data->getSize()> 0);
905
906            SQ_calc_and_apply_group_data(node2, gb_main, data, filter, progress); // enter right branch
907            seq_assert(data->getSize()> 0);
908        }
909    }
910    progress.inc();
911}
912
913void SQ_calc_and_apply_group_data2(GBT_TREE * node, GBDATA * gb_main, const SQ_GroupData * data, AP_filter * filter, arb_progress& progress) {
914    if (node->is_leaf) {
915        if (node->gb_node) {
916            SQ_pass2(data, gb_main, node, filter);
917        }
918    }
919    else {
920        GBT_TREE *node1 = node->leftson;
921        GBT_TREE *node2 = node->rightson;
922
923        if (node1) SQ_calc_and_apply_group_data2(node1, gb_main, data, filter, progress);
924        if (node2) SQ_calc_and_apply_group_data2(node2, gb_main, data, filter, progress);
925    }
926    progress.inc();
927}
928
929// marks species that are below threshold "evaluation"
930GB_ERROR SQ_mark_species(GBDATA * gb_main, int condition, bool marked_only) {
931    char *alignment_name;
932    int result = 0;
933
934    GBDATA *read_sequence = 0;
935    GBDATA *gb_species;
936    GB_ERROR error = 0;
937
938    GB_push_transaction(gb_main);
939    alignment_name = GBT_get_default_alignment(gb_main); seq_assert(alignment_name);
940
941    species_iterator  getFirst      = 0;
942    species_iterator  getNext       = 0;
943
944    if (marked_only) {
945        getFirst = GBT_first_marked_species;
946        getNext = GBT_next_marked_species;
947    }
948    else {
949        getFirst = GBT_first_species;
950        getNext = GBT_next_species;
951    }
952
953    for (gb_species = getFirst(gb_main); gb_species; gb_species = getNext(gb_species)) {
954        GBDATA *gb_ali = GB_entry(gb_species, alignment_name);
955        bool marked = false;
956        if (gb_ali) {
957            GBDATA *gb_quality = GB_search(gb_species, "quality", GB_CREATE_CONTAINER);
958            if (gb_quality) {
959                read_sequence = GB_entry(gb_ali, "data");
960                if (read_sequence) {
961                    GBDATA *gb_quality_ali = GB_search(gb_quality, alignment_name, GB_CREATE_CONTAINER);
962                    if (gb_quality_ali) {
963                        GBDATA *gb_result1 = GB_search(gb_quality_ali, "evaluation", GB_INT);
964                        result = GB_read_int(gb_result1);
965
966                        if (result < condition) marked = true;
967                    }
968                }
969            }
970        }
971
972        if (GB_read_flag(gb_species) != marked) {
973            GB_write_flag(gb_species, marked);
974        }
975    }
976    free(alignment_name);
977
978    if (error)
979        GB_abort_transaction(gb_main);
980    else
981        GB_pop_transaction(gb_main);
982
983    return error;
984}
985
986SQ_TREE_ERROR SQ_check_tree_structure(GBT_TREE * node) {
987    SQ_TREE_ERROR retval = NONE;
988
989    if (!node)
990        return MISSING_NODE;
991
992    if (node->is_leaf) {
993        if (!node->gb_node)
994            retval = ZOMBIE;
995    }
996    else {
997        retval = SQ_check_tree_structure(node->leftson);
998        if (retval == NONE)
999            retval = SQ_check_tree_structure(node->rightson);
1000    }
1001
1002    return retval;
1003}
Note: See TracBrowser for help on using the repository browser.