source: tags/svn.1.5.4/NALIGNER/ali_profile.hxx

Last change on this file was 7812, checked in by westram, 14 years ago

merge from dev [7751] [7752]

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