source: branches/port5/SEQ_QUALITY/SQ_helix.h

Last change on this file was 5675, checked in by westram, 16 years ago
  • removed automatic timestamps (the best they were good for, were vc-conflicts)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.3 KB
Line 
1//  ==================================================================== //
2//                                                                       //
3//    File      : SQ_helix.h                                             //
4//    Purpose   : Class used for calculation of helix layout             //
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#ifndef _CPP_STRING
16#include <string>
17#endif
18#ifndef _CPP_MAP
19#include <map>
20#endif
21
22#ifndef BI_HELIX_HXX
23#include <BI_helix.hxx>
24#endif
25
26class SQ_helix {
27public:
28    SQ_helix(int size);
29    ~SQ_helix();
30    void SQ_calc_helix_layout(const char *sequence, GBDATA * gb_main,
31            char *alignment_name, GBDATA * gb_quality, AP_filter * filter);
32    int SQ_get_no_helix() const {
33        return count_no_helix;
34    }
35    int SQ_get_weak_helix() const {
36        return count_weak_helix;
37    }
38    int SQ_get_strong_helix() const {
39        return count_strong_helix;
40    }
41private:
42    const char *sequence;
43    int count_strong_helix;
44    int count_weak_helix;
45    int count_no_helix;
46    int size;
47
48    static BI_helix *helix;
49    static GBDATA *helix_gb_main;
50    static std::string helix_ali_name;
51    static std::map<int, int> filterMap;
52    static bool has_filterMap;
53
54    static BI_helix & getHelix(GBDATA * gb_main, const char *ali_name);
55};
56
57// global data
58BI_helix *SQ_helix::helix = 0;
59GBDATA *SQ_helix::helix_gb_main = 0;
60std::string SQ_helix::helix_ali_name;
61std::map<int, int> SQ_helix::filterMap;
62bool SQ_helix::has_filterMap = false;
63
64// SQ_helix implementation
65BI_helix & SQ_helix::getHelix(GBDATA * gb_main, const char *ali_name) {
66    if (!helix || gb_main != helix_gb_main || strcmp(helix_ali_name.c_str(),
67            ali_name) != 0) {
68        delete helix;
69        helix = new BI_helix;
70
71        helix->init(gb_main, ali_name);
72
73        helix_gb_main = gb_main;
74        helix_ali_name = ali_name;
75    }
76    return *helix;
77}
78
79SQ_helix::SQ_helix(int size_) :
80    sequence(0), count_strong_helix(0), count_weak_helix(0), count_no_helix(0),
81            size(size_) {
82}
83
84SQ_helix::~SQ_helix() {
85}
86
87void SQ_helix::SQ_calc_helix_layout(const char *seq, GBDATA * gb_main,
88        char *alignment_name, GBDATA * gb_quality, AP_filter * filter) {
89    getHelix(gb_main, alignment_name);
90
91    // one call should be enough here (alignment does not change during the whole evaluation)
92    if (!has_filterMap) {
93        filterMap.clear();
94
95        for (int filter_pos = 0; filter_pos < filter->real_len; filter_pos++) {
96            filterMap[filter->filterpos_2_seqpos[filter_pos]] = filter_pos;
97        }
98
99        has_filterMap = true;
100    }
101
102    if (!helix->has_entries()) {
103        count_strong_helix = 1;
104        count_weak_helix = 1;
105        count_no_helix = 1;
106    } else {
107        // calculate the number of strong, weak and no helixes
108        std::map<int,int>::iterator it;
109
110        for (int filter_pos = 0; filter_pos < filter->real_len; filter_pos++) {
111            int seq_pos = filter->filterpos_2_seqpos[filter_pos];
112
113            BI_PAIR_TYPE pair_type = helix->pairtype(seq_pos);
114            if (pair_type == HELIX_PAIR) {
115                int v_seq_pos = helix->opposite_position(seq_pos);
116
117                if (v_seq_pos > seq_pos) { // ignore right helix positions
118                    it = filterMap.find(v_seq_pos);
119
120                    if (it != filterMap.end()) {
121                        char left = seq[filter_pos];
122                        char right = seq[it->second];
123                        int check = helix->check_pair(left, right, pair_type);
124
125                        switch (check) {
126                        case 2:
127                            count_strong_helix++;
128                            break;
129                        case 1:
130                            count_weak_helix++;
131                            break;
132                        case 0:
133                            if (!((left == '-') && (right == '-')))
134                                count_no_helix++;
135                            break;
136                        }
137                    }
138                }
139            }
140        }
141    }
142
143    //     /*if (count_strong_helix != 0)*/ count_strong_helix = count_strong_helix / 2;
144    //     /*if (count_weak_helix   != 0)*/ count_weak_helix = count_weak_helix / 2;
145    //     /*if (count_no_helix     != 0)*/ count_no_helix = count_no_helix / 2;
146
147    GBDATA *gb_result1 = GB_search(gb_quality, "number_of_no_helix", GB_INT);
148    seq_assert(gb_result1);
149    GB_write_int(gb_result1, count_no_helix);
150
151    GBDATA *gb_result2 = GB_search(gb_quality, "number_of_weak_helix", GB_INT);
152    seq_assert(gb_result2);
153    GB_write_int(gb_result2, count_weak_helix);
154
155    GBDATA *gb_result3 =
156            GB_search(gb_quality, "number_of_strong_helix", GB_INT);
157    seq_assert(gb_result3);
158    GB_write_int(gb_result3, count_strong_helix);
159}
Note: See TracBrowser for help on using the repository browser.