ARB
SQ_helix.h
Go to the documentation of this file.
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 _GLIBCXX_STRING
16 #include <string>
17 #endif
18 #ifndef _GLIBCXX_MAP
19 #include <map>
20 #endif
21 
22 #ifndef BI_HELIX_HXX
23 #include <BI_helix.hxx>
24 #endif
25 #ifndef AP_FILTER_HXX
26 #include <AP_filter.hxx>
27 #endif
28 #ifndef ARBDB_H
29 #include <arbdb.h>
30 #endif
31 
32 
33 class SQ_helix : virtual Noncopyable {
34  const char *sequence;
35  int count_strong_helix;
36  int count_weak_helix;
37  int count_no_helix;
38  int size;
39 
40  static SmartPtr<BI_helix> helix;
41  static GBDATA *helix_gb_main;
42  static std::string helix_ali_name;
43  static std::map<int, int> filterMap;
44  static bool has_filterMap;
45 
46  static BI_helix & getHelix(GBDATA * gb_main, const char *ali_name);
47 
48 public:
49  SQ_helix(int size);
50  ~SQ_helix();
51  void SQ_calc_helix_layout(const char *sequence, GBDATA * gb_main,
52  char *alignment_name, GBDATA * gb_quality, AP_filter * filter);
53  int SQ_get_no_helix() const {
54  return count_no_helix;
55  }
56  int SQ_get_weak_helix() const {
57  return count_weak_helix;
58  }
59  int SQ_get_strong_helix() const {
60  return count_strong_helix;
61  }
62 };
63 
64 // global data
65 SmartPtr<BI_helix> SQ_helix::helix;
66 GBDATA *SQ_helix::helix_gb_main = NULp;
67 std::string SQ_helix::helix_ali_name;
68 std::map<int, int> SQ_helix::filterMap;
69 bool SQ_helix::has_filterMap = false;
70 
71 // SQ_helix implementation
72 BI_helix & SQ_helix::getHelix(GBDATA * gb_main, const char *ali_name) {
73  if (helix.isNull() ||
74  gb_main != helix_gb_main ||
75  strcmp(helix_ali_name.c_str(), ali_name) != 0)
76  {
77  helix = new BI_helix;
78  helix->init(gb_main, ali_name);
79 
80  helix_gb_main = gb_main;
81  helix_ali_name = ali_name;
82  }
83  return *helix;
84 }
85 
86 SQ_helix::SQ_helix(int size_) :
87  sequence(NULp),
88  count_strong_helix(0),
89  count_weak_helix(0),
90  count_no_helix(0),
91  size(size_)
92 {}
93 
95 }
96 
97 void SQ_helix::SQ_calc_helix_layout(const char *seq, GBDATA * gb_main,
98  char *alignment_name, GBDATA * gb_quality, AP_filter * filter) {
99  getHelix(gb_main, alignment_name);
100 
101  size_t filterLen = filter->get_filtered_length();
102  const size_t *filterpos_2_seqpos = filter->get_filterpos_2_seqpos();
103 
104  // one call should be enough here (alignment does not change during the whole evaluation)
105  if (!has_filterMap) {
106  filterMap.clear();
107 
108  for (size_t filter_pos = 0; filter_pos < filterLen; filter_pos++) {
109  filterMap[filterpos_2_seqpos[filter_pos]] = filter_pos;
110  }
111 
112  has_filterMap = true;
113  }
114 
115  if (!helix->has_entries()) {
116  count_strong_helix = 1;
117  count_weak_helix = 1;
118  count_no_helix = 1;
119  }
120  else {
121  // calculate the number of strong, weak and no helixes
122  std::map<int, int>::iterator it;
123 
124  for (size_t filter_pos = 0; filter_pos < filterLen; filter_pos++) {
125  int seq_pos = filterpos_2_seqpos[filter_pos];
126 
127  BI_PAIR_TYPE pair_type = helix->pairtype(seq_pos);
128  if (pair_type == HELIX_PAIR) {
129  int v_seq_pos = helix->opposite_position(seq_pos);
130 
131  if (v_seq_pos > seq_pos) { // ignore right helix positions
132  it = filterMap.find(v_seq_pos);
133 
134  if (it != filterMap.end()) {
135  char left = seq[filter_pos];
136  char right = seq[it->second];
137  int check = helix->check_pair(left, right, pair_type);
138 
139  switch (check) {
140  case 2:
141  count_strong_helix++;
142  break;
143  case 1:
144  count_weak_helix++;
145  break;
146  case 0:
147  if (!((left == '-') && (right == '-')))
148  count_no_helix++;
149  break;
150  }
151  }
152  }
153  }
154  }
155  }
156 
157  GBDATA *gb_result1 = GB_search(gb_quality, "number_of_no_helix", GB_INT);
158  seq_assert(gb_result1);
159  GB_write_int(gb_result1, count_no_helix);
160 
161  GBDATA *gb_result2 = GB_search(gb_quality, "number_of_weak_helix", GB_INT);
162  seq_assert(gb_result2);
163  GB_write_int(gb_result2, count_weak_helix);
164 
165  GBDATA *gb_result3 =
166  GB_search(gb_quality, "number_of_strong_helix", GB_INT);
167  seq_assert(gb_result3);
168  GB_write_int(gb_result3, count_strong_helix);
169 }
const size_t * get_filterpos_2_seqpos() const
Definition: AP_filter.hxx:91
return string(buffer, length)
int SQ_get_no_helix() const
Definition: SQ_helix.h:53
int SQ_get_strong_helix() const
Definition: SQ_helix.h:59
BI_PAIR_TYPE
Definition: BI_helix.hxx:18
size_t opposite_position(size_t pos) const
Definition: BI_helix.hxx:96
static char * alignment_name
int SQ_get_weak_helix() const
Definition: SQ_helix.h:56
bool isNull() const
test if SmartPtr is NULp
Definition: smartptr.h:248
FILE * seq
Definition: rns.c:46
~SQ_helix()
Definition: SQ_helix.h:94
void SQ_calc_helix_layout(const char *sequence, GBDATA *gb_main, char *alignment_name, GBDATA *gb_quality, AP_filter *filter)
Definition: SQ_helix.h:97
#define seq_assert(bed)
GB_ERROR GB_write_int(GBDATA *gbd, long i)
Definition: arbdb.cxx:1244
size_t get_filtered_length() const
Definition: AP_filter.hxx:82
SQ_helix(int size)
Definition: SQ_helix.h:86
bool has_entries() const
Definition: BI_helix.hxx:89
int check_pair(char left, char right, BI_PAIR_TYPE pair_type)
Definition: BI_helix.cxx:322
BI_PAIR_TYPE pairtype(size_t pos) const
Definition: BI_helix.hxx:101
#define NULp
Definition: cxxforward.h:97
GBDATA * gb_main
Definition: adname.cxx:33
GBDATA * GB_search(GBDATA *gbd, const char *fieldpath, GB_TYPES create)
Definition: adquery.cxx:531
Definition: arbdb.h:66