source: tags/arb-6.0/NALIGNER/ali_profile.hxx

Last change on this file was 11002, checked in by westram, 11 years ago
  • 'class { public' → struct
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.6 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : ali_profile.hxx                                   //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#ifndef ALI_PROFILE_HXX
12#define ALI_PROFILE_HXX
13
14#ifndef ALI_ARBDB_HXX
15#include "ali_arbdb.hxx"
16#endif
17#ifndef ALI_PT_HXX
18#include "ali_pt.hxx"
19#endif
20
21#define ALI_PROFILE_BORDER_LEFT '['
22#define ALI_PROFILE_BORDER_RIGHT ']'
23
24typedef void ALI_MAP_; // make module independent
25
26
27struct ALI_PROFILE_CONTEXT {
28    ALI_ARBDB *arbdb;
29    ALI_PT *pt;
30
31    int find_family_mode;
32
33    int exclusive_flag;
34    int mark_family_flag;
35    int mark_family_extension_flag;
36
37    int min_family_size;
38    int max_family_size;
39
40    float min_weight;
41    float ext_max_weight;
42
43    float multi_gap_factor;
44
45    float insert_factor;
46    float multi_insert_factor;
47
48    double substitute_matrix[5][5]; // a c g u -
49    double binding_matrix[5][5]; // a c g u -
50};
51
52
53// Class for a family member
54struct ali_family_member {
55    ALI_SEQUENCE *sequence;
56    float matches;
57    float weight;
58
59    ali_family_member(ALI_SEQUENCE *seq, float m, float w = 0.0) {
60        sequence = seq;
61        matches = m;
62        weight = w;
63    }
64};
65
66
67// Class for the profiling
68class ALI_PROFILE : virtual Noncopyable {
69    ALI_NORM_SEQUENCE *norm_sequence;
70    unsigned long prof_len;
71
72    long **helix;                    // base to base connection
73    char **helix_borders;            // borders of the helix '[' and ']'
74    unsigned long helix_len;
75
76    float (**base_weights)[4];         // relative weight of base i
77    float (**sub_costs)[6];            // costs to substitute with base i
78    float sub_costs_maximum;
79
80    float insert_cost;
81    float multi_insert_cost;
82
83    float multi_gap_factor;
84
85    float (*binding_costs)[5][5];       // Matrix for binding costs a c g u -
86    float w_bind_maximum;
87
88    unsigned long *lmin, *lmax;
89    float ***gap_costs;
90    float ***gap_percents;
91
92    int is_binding_marker(char c);
93
94    ALI_TLIST<ali_family_member *> *find_family(ALI_SEQUENCE *sequence,
95                                                ALI_PROFILE_CONTEXT *context);
96    void calculate_costs(ALI_TLIST<ali_family_member *> *family_list,
97                         ALI_PROFILE_CONTEXT *context);
98
99    int find_next_helix(char h[], unsigned long h_len, unsigned long pos,
100                        unsigned long *helix_nr,
101                        unsigned long *start, unsigned long *end);
102    int find_comp_helix(char h[], unsigned long h_len, unsigned long pos,
103                        unsigned long helix_nr,
104                        unsigned long *start, unsigned long *end);
105    void delete_comp_helix(char h1[], char h2[], unsigned long h_len,
106                           unsigned long start, unsigned long end);
107    // int map_helix(char h[], unsigned long h_len,
108                  // unsigned long start1, unsigned long end1,
109                  // unsigned long start2, unsigned long end2);
110    void initialize_helix(ALI_PROFILE_CONTEXT *context);
111
112public:
113
114    ALI_PROFILE(ALI_SEQUENCE *sequence, ALI_PROFILE_CONTEXT *context);
115    ~ALI_PROFILE();
116
117    void print(int start = -1, int end = -1)
118    {
119        int i;
120        size_t j;
121
122        if (start < 0 || (unsigned long)(start) > prof_len - 1)
123            start = 0;
124
125        if ((unsigned long)(end) > prof_len - 1 || end < 0)
126            end = (int)prof_len - 1;
127
128        if (start > end) {
129            i = start;
130            start = end;
131            end = i;
132        }
133
134
135        printf("\nProfile from %d to %d\n", start, end);
136        printf("Substitutions kosten:\n");
137        for (i = start; i <= end; i++) {
138            printf("%2d: ", i);
139            for (j = 0; j < 6; j++)
140                printf("%4.2f ", (*sub_costs)[i][j]);
141            printf("\n");
142        }
143
144        printf("\nGap Bereiche:\n");
145        for (i = start; i <= end; i++) {
146            printf("%2d: %3ld %3ld\n", i, lmin[i], lmax[i]);
147            printf("  : ");
148            for (j = 0; j <= lmax[i] - lmin[i] + 1; j++)
149                printf("%4.2f ", (*gap_percents)[i][j]);
150            printf("\n  : ");
151            for (j = 0; j <= lmax[i] - lmin[i] + 1; j++)
152                printf("%4.2f ", (*gap_costs)[i][j]);
153            printf("\n");
154        }
155
156    }
157
158    int is_in_helix(unsigned long pos,
159                    unsigned long *first, unsigned long *last);
160    int is_outside_helix(unsigned long pos,
161                         unsigned long *first, unsigned long *last);
162    int complement_position(unsigned long pos) {
163        if (pos >= helix_len)
164            return -1;
165        else
166            return (*helix)[pos];
167    }
168    int is_in_helix(unsigned long pos) {
169        long l;
170        if ((*helix_borders)[pos] == ALI_PROFILE_BORDER_LEFT ||
171            (*helix_borders)[pos] == ALI_PROFILE_BORDER_RIGHT)
172            return 1;
173        for (l = (long) pos - 1; l >= 0; l--)
174            switch ((*helix_borders)[l]) {
175                case ALI_PROFILE_BORDER_LEFT:
176                    return 1;
177                case ALI_PROFILE_BORDER_RIGHT:
178                    return 0;
179            }
180        for (l = (long) pos + 1; l < (long) prof_len; l++)
181            switch ((*helix_borders)[l]) {
182                case ALI_PROFILE_BORDER_LEFT:
183                    return 0;
184                case ALI_PROFILE_BORDER_RIGHT:
185                    return 1;
186            }
187        return 0;
188    }
189
190    int is_internal_last(unsigned long pos) {
191        return norm_sequence->is_begin(pos + 1);
192    }
193
194    char *cheapest_sequence();
195    char *borders_sequence() {
196        unsigned long i;
197        char *str, *str_buffer;
198
199        str = (char *) calloc((unsigned int) prof_len, sizeof(char));
200        if (str == 0) ali_fatal_error("Out of memory");
201
202        str_buffer = str;
203        for (i = 0; i < prof_len; i++)
204            switch ((*helix_borders)[i]) {
205                case ALI_PROFILE_BORDER_LEFT: *str++ = '['; break;
206                case ALI_PROFILE_BORDER_RIGHT: *str++ = ']'; break;
207                default: *str++ = ' ';
208            }
209        return str_buffer;
210    }
211
212    unsigned long length() {
213        return prof_len;
214    }
215    ALI_NORM_SEQUENCE *sequence() {
216        return norm_sequence;
217    }
218    unsigned long sequence_length() {
219        return norm_sequence->length();
220    }
221
222    float base_weight(unsigned long pos, unsigned char c) {
223        if (c > 3)
224            ali_fatal_error("Out of range", "ALI_PROFILE::base_weight()");
225        return (*base_weights)[pos][c];
226    }
227
228    float w_sub(unsigned long positiony, unsigned long positionx) {
229        return (*sub_costs)[positiony][norm_sequence->base(positionx)];
230    }
231    float w_sub_gap(unsigned long positiony) {
232        return (*sub_costs)[positiony][4];
233    }
234    float w_sub_gap_cheap(unsigned long positiony) {
235        return (*sub_costs)[positiony][4] * multi_gap_factor;
236    }
237    float w_sub_multi_gap_cheap(unsigned long start, unsigned long end) {
238        float costs = 0.0;
239        unsigned long i;
240        for (i = start; i <= end; i++)
241            costs += w_sub_gap(i);
242        return costs * multi_gap_factor;
243    }
244
245    float w_sub_maximum() {
246        return sub_costs_maximum;
247    }
248
249    float w_del(unsigned long begin, unsigned long end) {
250        if (begin <= lmin[end])
251            return (*gap_costs)[end][0];
252        else {
253            if (begin <= lmax[end])
254                return (*gap_costs)[end][begin - lmin[end]];
255            else
256                return (*gap_costs)[end][lmax[end] - lmin[end] + 1];
257        }
258    }
259    float w_del_cheap(unsigned long position) {
260        return (*gap_costs)[position][0];
261    }
262    float w_del_multi(unsigned long start, unsigned long end) {
263        float costs = 0.0;
264        unsigned long i;
265        for (i = start; i <= end; i++)
266            costs += w_del(start, i);
267        return costs;
268    }
269    float w_del_multi_cheap(unsigned long start, unsigned long end) {
270        float costs = 0.0;
271        unsigned long i;
272        for (i = start; i <= end; i++)
273            costs += w_del_cheap(i);
274        return costs;
275    }
276    float w_del_multi_unweighted(unsigned long start, unsigned long end) {
277        float costs = 0.0;
278        unsigned long i;
279        for (i = start; i <= end; i++)
280            costs += w_del(i, i);
281        return costs;
282    }
283
284    float w_ins(unsigned long /* positionx */, unsigned long /* positiony */) {
285        return insert_cost;
286    }
287    float w_ins_cheap(unsigned long /* positionx */, unsigned long /* positiony */) {
288        return multi_insert_cost;
289    }
290    float w_ins_multi_unweighted(unsigned long startx, unsigned long endx) {
291        return (endx - startx + 1) * insert_cost;
292    }
293    float w_ins_multi(unsigned long startx, unsigned long endx) {
294        return insert_cost + (endx - startx)*multi_insert_cost;
295    }
296    float w_ins_multi_cheap(unsigned long startx, unsigned long endx) {
297        return (endx - startx + 1)*multi_insert_cost;
298    }
299
300
301    float gap_percent(unsigned long begin, unsigned long end) {
302        if (begin <= lmin[end])
303            return (*gap_percents)[end][0];
304        else {
305            if (begin <= lmax[end])
306                return (*gap_percents)[end][begin - lmin[end]];
307            else
308                return (*gap_percents)[end][lmax[end] - lmin[end] + 1];
309        }
310    }
311
312    float w_bind(unsigned long pos_1, unsigned char base_1,
313                 unsigned long pos_2, unsigned char base_2) {
314        int i, j;
315        float cost_in, cost = 0;
316        if (ali_is_real_base_or_gap(base_1)) {
317            if (ali_is_real_base_or_gap(base_2)) {
318                return (*binding_costs)[base_1][base_2];
319            }
320            for (i = 0; i < 4; i++)
321                cost += (*base_weights)[pos_2][i] * (*binding_costs)[base_1][i];
322            return cost;
323        }
324        if (ali_is_real_base_or_gap(base_2)) {
325            for (i = 0; i < 4; i++)
326                cost += (*base_weights)[pos_1][i] * (*binding_costs)[i][base_2];
327            return cost;
328        }
329        for (i = 0; i < 4; i++) {
330            for (j = 0, cost_in = 0; j < 4; j++)
331                cost_in += (*base_weights)[pos_2][j] * (*binding_costs)[i][j];
332            cost += (*base_weights)[pos_1][i] * cost_in;
333        }
334        return cost;
335    }
336
337    float w_binding(unsigned long first_position_of_sequence,
338                    ALI_SEQUENCE *sequence);
339};
340
341#else
342#error ali_profile.hxx included twice
343#endif // ALI_PROFILE_HXX
Note: See TracBrowser for help on using the repository browser.