source: tags/svn.1.5.4/SEQ_QUALITY/SQ_functions.cxx

Last change on this file was 8313, checked in by westram, 12 years ago
  • removed dead code

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