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 | |
---|
26 | class SQ_helix { |
---|
27 | public: |
---|
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 | } |
---|
41 | private: |
---|
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 |
---|
58 | BI_helix *SQ_helix::helix = 0; |
---|
59 | GBDATA *SQ_helix::helix_gb_main = 0; |
---|
60 | std::string SQ_helix::helix_ali_name; |
---|
61 | std::map<int, int> SQ_helix::filterMap; |
---|
62 | bool SQ_helix::has_filterMap = false; |
---|
63 | |
---|
64 | // SQ_helix implementation |
---|
65 | BI_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 | |
---|
79 | SQ_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 | |
---|
84 | SQ_helix::~SQ_helix() { |
---|
85 | } |
---|
86 | |
---|
87 | void 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 | } |
---|