source: tags/ms_r16q3/NALIGNER/ali_main.cxx

Last change on this file was 9899, checked in by epruesse, 11 years ago

pass down argc/argv to AW_root::AW_root (for GTK)
GTK will pick its parameters (e.g. —g-fatal-warnings) from argv and
remove them, hence argv must not be const and AW_root() should
be called before command line parsing (if possible).

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.7 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : ali_main.cxx                                      //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "ali_global.hxx"
12
13
14#define HELIX_PAIRS   "helix_pairs"
15#define HELIX_LINE    "helix_line"
16#define ALI_CONSENSUS "ALI_CON"
17#define ALI_ERROR     "ALI_ERR"
18#define ALI_INTERVALS "ALI_INT"
19
20static ALI_GLOBAL aligs;
21
22
23void message(char *errortext);
24
25
26static const char *ali_version =
27    "\nALIGNER   V2.0  (Boris Reichel 5/95)\n";
28
29static const char *ali_man_line[] = {
30    "Parameter:",
31    "-s<species>       aligne the sequence of <species>",
32    "-f<species_1>,...,<species_n>[;<extension_1>,...,<extension_k>]  use specified family and family extension",
33    "-P<pt_server>     use the PT_server <pt_server>",
34    "[-D<db_server>]   use the DB_server <db_server>",
35    "[-nx]             not exclusive mode (profile)",
36    "[-ms]             mark species (profile)",
37    "[-mf]             mark used family (profile)",
38    "[-mfe]            mark used family extension (profile)",
39    "[-mgf]            multi gap factor (0.1) (profile)",
40    "[-if]             insert factor (2.0) (profile)",
41    "[-mif]            multi insert factor (0.5) (profile)",
42    "[-m]              mark all (-m  == -ms -mf -mfe) (profile)",
43    "[-d<filename>]    use <filename> for default values",
44    "[-msubX1,X2,...,X25] use X1,...,X25 for the substitute matrix: (profile)",
45    "                     a   c   g   u   -",
46    "                a   X1  X2  X3  X4  X5",
47    "                c   X6  X7  X8  X9  X10",
48    "                g   X11 X12 X13 X14 X15",
49    "                u   X16 X17 X18 X19 X20",
50    "                -   X21 X22 X23 X24 X25",
51    "[-mbindX1,X2,...,X25] use X1,...,X25 for the binding matrix: (profile)",
52    "                     a   c   g   u   -",
53    "                a   X1  X2  X3  X4  X5",
54    "                c   X6  X7  X8  X9  X10",
55    "                g   X11 X12 X13 X14 X15",
56    "                u   X16 X17 X18 X19 X20",
57    "                -   X21 X22 X23 X24 X25",
58    "[-maxf]           maximal number of family members (10) (profile)",
59    "[-minf]           minimal number of family members (5) (profile)",
60    "[-minw]           minimal weight for family members (0.7) (profile)",
61    "[-maxew]          maximal weight for family extension members (0.2) (profile)",
62    "unused [-cl]             cost threshold (low) (0.25) ALI_ERR = ','",
63    "unused [-cm]             cost threshold (middle) (0.5) ALI_ERR = '-'",
64    "unused [-ch]             cost threshold (high) (0.8) ALI_ERR = '='",
65    "[-mm]             maximal number of maps for solution (prealigner) (1000)",
66    "[-mma]            maximal number of maps for aligner (prealigner) (2)",
67    "[-csub]           cost threshold for substitution (prealigner) (0.5)",
68    "[-chel]           cost threshold for helix binding (prealigner) (2.0)",
69    "[-ec]             error count (prealigner) (2)",
70    "[-ib]             interval border (prealigner) (5)",
71    "[-ic]             interval center (prealigner) (5)",
72    0
73};
74
75
76static void print_man() {
77    // Print a short parameter description
78    int i;
79
80    for (i = 0; ali_man_line[i] != 0; i++)
81        fprintf(stderr, "%s\n", ali_man_line[i]);
82}
83
84
85void ali_fatal_error(const char *message, const char *func) {
86    fprintf(stderr, "FATAL ERROR %s: %s\n", func, message);
87    exit(-1);
88}
89
90void ali_error(const char *message, const char *func) {
91    fprintf(stderr, "ERROR %s: %s\n", func, message);
92    exit(-1);
93}
94
95
96
97static int get_species(char *species_string, unsigned int species_number, char *buffer) {
98    // Get one species of a list
99    while (species_number > 0 && *species_string != '\0') {
100        while (*species_string != '\0' && *species_string != ',')
101            species_string++;
102        if (*species_string != '\0')
103            species_string++;
104        species_number--;
105    }
106
107    if (*species_string != '\0') {
108        while (*species_string != '\0' && *species_string != ',')
109            *buffer++ = *species_string++;
110        *buffer = '\0';
111    }
112    else {
113        return 0;
114    }
115
116    return 1;
117}
118
119
120static int check_base_invariance(char *seq1, char *seq2)
121{
122    while (*seq1 != '\0' && !ali_is_base(*seq1))
123        seq1++;
124    while (*seq2 != '\0' && !ali_is_base(*seq2))
125        seq2++;
126    while (*seq1 != '\0' && *seq2 != '\0') {
127        if (*seq1 != *seq2)
128            return 0;
129        seq1++;
130        seq2++;
131        while (*seq1 != '\0' && !ali_is_base(*seq1))
132            seq1++;
133        while (*seq2 != '\0' && !ali_is_base(*seq2))
134            seq2++;
135    }
136
137    if (*seq1 == *seq2)
138        return 1;
139
140    return 0;
141}
142
143
144static int convert_for_back_write(char *seq_new, char *seq_orig) {
145    // Convert the working sequenz into the original bases
146
147    while (*seq_new != '\0' && (ali_is_dot(*seq_new) || ali_is_gap(*seq_new)))
148        seq_new++;
149    while (*seq_orig != '\0' && (ali_is_dot(*seq_orig) || ali_is_gap(*seq_orig)))
150        seq_orig++;
151
152    while (*seq_new != '\0' && *seq_orig != '\0') {
153        if (*seq_new != *seq_orig) {
154            switch (*seq_new) {
155                case 'a':
156                    switch (*seq_orig) {
157                        case 'A':
158                            *seq_new = 'A';
159                            break;
160                        default:
161                            ali_error("Unexpected character in original sequence");
162                    }
163                    break;
164                case 'c':
165                    switch (*seq_orig) {
166                        case 'C':
167                            *seq_new = 'C';
168                            break;
169                        default:
170                            ali_error("Unexpected character in original sequence");
171                    }
172                    break;
173                case 'g':
174                    switch (*seq_orig) {
175                        case 'G':
176                            *seq_new = 'G';
177                            break;
178                        default:
179                            ali_error("Unexpected character in original sequence");
180                    }
181                    break;
182                case 'u':
183                    switch (*seq_orig) {
184                        case 'U':
185                            *seq_new = 'U';
186                            break;
187                        case 't':
188                            *seq_new = 't';
189                            break;
190                        case 'T':
191                            *seq_new = 'T';
192                            break;
193                    }
194                    break;
195                case 'n':
196                    *seq_new = *seq_orig;
197                    break;
198                default:
199                    ali_fatal_error("Unexpected character in generated sequence");
200            }
201        }
202        seq_new++;
203        seq_orig++;
204        while (*seq_new != '\0' && (ali_is_dot(*seq_new) || ali_is_gap(*seq_new)))
205            seq_new++;
206        while (*seq_orig != '\0' && (ali_is_dot(*seq_orig) || ali_is_gap(*seq_orig)))
207            seq_orig++;
208    }
209
210    if (*seq_new == *seq_orig)
211        return 1;
212
213    return 0;
214}
215
216
217
218int ARB_main(int argc, char *argv[]) {
219    int                            i;
220    char                           message_buffer[100];
221    char                           species_name[100], species_number;
222    ALI_PREALIGNER                *align_prealigner;
223    ali_prealigner_approx_element *approx_elem;
224    ALI_SEQUENCE                  *sequence;
225
226
227    ali_message(ali_version);
228
229    aligs.init(&argc, (const char**)argv);
230
231    if (!aligs.species_name || argc > 1) {
232        printf("Unknowen : ");
233        for (i = 1; i < argc; i++)
234            printf("%s ", argv[i]);
235        printf("\n\n");
236        print_man();
237        exit (-1);
238    }
239
240    // Main loop
241    species_number = 0;
242    while (get_species(aligs.species_name, species_number, species_name)) {
243        species_number++;
244
245        sprintf(message_buffer,
246                "\nStarting alignment of sequence: %s", species_name);
247        ali_message(message_buffer);
248
249        // Get all information of the sequence
250        aligs.arbdb.begin_transaction();
251        ALI_SEQUENCE *align_sequence;
252        align_sequence = aligs.arbdb.get_sequence(species_name,
253                                                  aligs.mark_species_flag);
254        if (align_sequence == 0) {
255            ali_error("Can't read sequence from database");
256            ali_message("Aborting alignment of sequence");
257        }
258        else {
259            char *align_string;
260            char *align_string_original;
261
262            align_string = align_sequence->string();
263            align_string_original = aligs.arbdb.get_sequence_string(species_name,
264                                                                    aligs.mark_species_flag);
265            aligs.arbdb.commit_transaction();
266
267            if (align_sequence == 0)
268                ali_warning("Can't read sequence from database");
269            else {
270                // make profile for sequence
271                ALI_PROFILE *align_profile;
272                align_profile = new ALI_PROFILE(align_sequence, &aligs.prof_context);
273
274                // write information about the profile to the database
275                aligs.arbdb.begin_transaction();
276                char *String = align_profile->cheapest_sequence();
277                aligs.arbdb.put_SAI("ALI_CON", String);
278                freeset(String, align_profile->borders_sequence());
279                free(String);
280                aligs.arbdb.commit_transaction();
281
282                // make prealignment
283                align_prealigner = new ALI_PREALIGNER(&aligs.preali_context,
284                                                      align_profile,
285                                                      0, align_profile->sequence_length() - 1,
286                                                      0, align_profile->length() - 1);
287                ALI_SEQUENCE *align_pre_sequence_i, *align_pre_sequence;
288                ALI_SUB_SOLUTION *align_pre_solution;
289                ALI_TLIST<ali_prealigner_approx_element *> *align_pre_approx;
290                align_pre_sequence_i = align_prealigner->sequence();
291                align_pre_sequence = align_prealigner->sequence_without_inserts();
292                align_pre_solution = align_prealigner->solution();
293                align_pre_approx = align_prealigner->approximation();
294                delete align_prealigner;
295
296                align_pre_solution->print();
297
298                // write result of alignment into database
299                aligs.arbdb.begin_transaction();
300                String = align_pre_sequence_i->string();
301                aligs.arbdb.put_SAI("ALI_PRE_I", String);
302                freeset(String, align_pre_sequence->string());
303                aligs.arbdb.put_SAI("ALI_PRE", String);
304                free(String);
305                aligs.arbdb.commit_transaction();
306
307
308                sprintf(message_buffer, "%d solutions generated (taking the first)",
309                        align_pre_approx->cardinality());
310                ali_message(message_buffer);
311
312                if (align_pre_approx->is_empty())
313                    ali_fatal_error("List of approximations is empty");
314
315                // Write result back to the database
316                approx_elem = align_pre_approx->first();
317
318                sequence = approx_elem->map->sequence(align_profile->sequence());
319                String = sequence->string();
320
321                if (!check_base_invariance(String, align_string))
322                    ali_error("Bases changed in output sequence");
323
324                if (!convert_for_back_write(String, align_string_original))
325                    ali_fatal_error("Can't convert correctly");
326
327                aligs.arbdb.begin_transaction();
328                aligs.arbdb.put_sequence_string(species_name, String);
329                aligs.arbdb.put_SAI("ALI_INSERTS", approx_elem->ins_marker);
330                aligs.arbdb.commit_transaction();
331                delete sequence;
332
333                // Delete all Objects
334                free(align_string);
335                free(align_string_original);
336                delete align_pre_solution;
337                delete align_pre_approx;
338                delete align_profile;
339                delete align_pre_sequence;
340                delete align_pre_sequence_i;
341            }
342            delete align_sequence;
343        }
344    }
345
346    ali_message("Aligner terminated\n");
347    return 0;
348}
Note: See TracBrowser for help on using the repository browser.