source: branches/profile/SEQ_QUALITY/SQ_helix.h

Last change on this file was 11002, checked in by westram, 6 years ago
  • 'class { public' → struct
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.5 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 _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
33class 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
48public:
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
65SmartPtr<BI_helix>  SQ_helix::helix;
66GBDATA             *SQ_helix::helix_gb_main = 0;
67std::string         SQ_helix::helix_ali_name;
68std::map<int, int>  SQ_helix::filterMap;
69bool                SQ_helix::has_filterMap = false;
70
71// SQ_helix implementation
72BI_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
86SQ_helix::SQ_helix(int size_)
87    : sequence(0),
88      count_strong_helix(0),
89      count_weak_helix(0),
90      count_no_helix(0),
91      size(size_)
92{
93}
94
95SQ_helix::~SQ_helix() {
96}
97
98void SQ_helix::SQ_calc_helix_layout(const char *seq, GBDATA * gb_main,
99                                    char *alignment_name, GBDATA * gb_quality, AP_filter * filter) {
100    getHelix(gb_main, alignment_name);
101
102    size_t        filterLen          = filter->get_filtered_length();
103    const size_t *filterpos_2_seqpos = filter->get_filterpos_2_seqpos();
104
105    // one call should be enough here (alignment does not change during the whole evaluation)
106    if (!has_filterMap) {
107        filterMap.clear();
108
109        for (size_t filter_pos = 0; filter_pos < filterLen; filter_pos++) {
110            filterMap[filterpos_2_seqpos[filter_pos]] = filter_pos;
111        }
112
113        has_filterMap = true;
114    }
115
116    if (!helix->has_entries()) {
117        count_strong_helix = 1;
118        count_weak_helix = 1;
119        count_no_helix = 1;
120    }
121    else {
122        // calculate the number of strong, weak and no helixes
123        std::map<int, int>::iterator it;
124
125        for (size_t filter_pos = 0; filter_pos < filterLen; filter_pos++) {
126            int seq_pos = filterpos_2_seqpos[filter_pos];
127
128            BI_PAIR_TYPE pair_type = helix->pairtype(seq_pos);
129            if (pair_type == HELIX_PAIR) {
130                int v_seq_pos = helix->opposite_position(seq_pos);
131
132                if (v_seq_pos > seq_pos) { // ignore right helix positions
133                    it = filterMap.find(v_seq_pos);
134
135                    if (it != filterMap.end()) {
136                        char left = seq[filter_pos];
137                        char right = seq[it->second];
138                        int check = helix->check_pair(left, right, pair_type);
139
140                        switch (check) {
141                        case 2:
142                            count_strong_helix++;
143                            break;
144                        case 1:
145                            count_weak_helix++;
146                            break;
147                        case 0:
148                            if (!((left == '-') && (right == '-')))
149                                count_no_helix++;
150                            break;
151                        }
152                    }
153                }
154            }
155        }
156    }
157
158    GBDATA *gb_result1 = GB_search(gb_quality, "number_of_no_helix", GB_INT);
159    seq_assert(gb_result1);
160    GB_write_int(gb_result1, count_no_helix);
161
162    GBDATA *gb_result2 = GB_search(gb_quality, "number_of_weak_helix", GB_INT);
163    seq_assert(gb_result2);
164    GB_write_int(gb_result2, count_weak_helix);
165
166    GBDATA *gb_result3 =
167            GB_search(gb_quality, "number_of_strong_helix", GB_INT);
168    seq_assert(gb_result3);
169    GB_write_int(gb_result3, count_strong_helix);
170}
Note: See TracBrowser for help on using the repository browser.