root/trunk/ARBDB/adali.cxx

Revision 8607, 34.5 KB (checked in by westram, 5 weeks ago)

merge from e4fix [8135] [8136] [8137] [8138] [8139] [8140] [8141] [8142] [8143] [8144] [8222]
this revives the reverted patches [8129] [8130] [8131] [8132]

  • fixes
    • some free/delete mismatches
    • wrong definition of ORF objects (Level was no bit value)
    • amino consensus (failed for columns only containing 'C')
  • rename
    • AA_sequence_term -> orf_term
    • ED4_sequence_terminal_basic -> ED4_abstract_sequence_terminal
  • cleaned up hierarchy dumps
  • tweaked is_terminal()/to_terminal()
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : adali.cxx                                         //
4//   Purpose   : alignments                                        //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include <arbdbt.h>
12#include <adGene.h>
13#include <arb_strarray.h>
14
15#include "gb_local.h"
16
17static long check_for_species_without_data(const char *species_name, long value, void *counterPtr) {
18    if (value == 1) {
19        long cnt = *((long*)counterPtr);
20        if (cnt<40) {
21            GB_warningf("Species '%s' has no data in any alignment", species_name);
22        }
23        *((long*)counterPtr) = cnt+1;
24    }
25    return value; // new hash value
26}
27
28#if defined(WARN_TODO)
29#warning GBT_check_data ignores given 'alignment_name' if we have a default alignment. seems wrong!
30#endif
31
32GBDATA *GBT_get_presets(GBDATA *gb_main) {
33    return GBT_find_or_create(gb_main, "presets", 7);
34}
35
36int GBT_count_alignments(GBDATA *gb_main) {
37    int     count      = 0;
38    GBDATA *gb_presets = GBT_get_presets(gb_main);
39    for (GBDATA *gb_ali = GB_entry(gb_presets, "alignment");
40         gb_ali;
41         gb_ali = GB_nextEntry(gb_ali))
42    {
43        ++count;
44    }
45    return count;
46}
47
48static GB_ERROR GBT_check_alignment(GBDATA *gb_main, GBDATA *preset_alignment, GB_HASH *species_name_hash) {
49    /* check
50     * - whether alignment has correct size and
51     * - whether all data is present.
52     *
53     * Sets the security deletes and writes.
54     *
55     * If 'species_name_hash' is not NULL,
56     * - it initially has to contain value == 1 for each existing species.
57     * - afterwards it will contain value  == 2 for each species where an alignment has been found.
58     */
59
60    GBDATA *gb_species_data  = GBT_get_species_data(gb_main);
61    GBDATA *gb_extended_data = GBT_get_SAI_data(gb_main);
62
63    GB_ERROR  error      = 0;
64    char     *ali_name   = GBT_read_string(preset_alignment, "alignment_name");
65    if (!ali_name) error = "Alignment w/o 'alignment_name'";
66
67    if (!error) {
68        long    security_write = -1;
69        long    stored_ali_len = -1;
70        long    found_ali_len  = -1;
71        long    aligned        = 1;
72        GBDATA *gb_ali_len     = 0;
73
74        {
75            GBDATA *gb_ali_wsec = GB_entry(preset_alignment, "alignment_write_security");
76            if (!gb_ali_wsec) {
77                error = "has no 'alignment_write_security' entry";
78            }
79            else {
80                security_write = GB_read_int(gb_ali_wsec);
81            }
82        }
83
84
85        if (!error) {
86            gb_ali_len = GB_entry(preset_alignment, "alignment_len");
87            if (!gb_ali_len) {
88                error = "has no 'alignment_len' entry";
89            }
90            else {
91                stored_ali_len = GB_read_int(gb_ali_len);
92            }
93        }
94
95        if (!error) {
96            GBDATA *gb_species;
97            for (gb_species = GBT_first_species_rel_species_data(gb_species_data);
98                 gb_species && !error;
99                 gb_species = GBT_next_species(gb_species))
100            {
101                GBDATA     *gb_name        = GB_entry(gb_species, "name");
102                const char *name           = 0;
103                int         alignment_seen = 0;
104                GBDATA     *gb_ali         = 0;
105
106                if (!gb_name) {
107                    // fatal: name is missing -> create a unique name
108                    char *unique = GBT_create_unique_species_name(gb_main, "autoname.");
109                    error        = GBT_write_string(gb_species, "name", unique);
110
111                    if (!error) {
112                        gb_name = GB_entry(gb_species, "name");
113                        GBS_write_hash(species_name_hash, unique, 1); // not seen before
114                        GB_warningf("Seen unnamed species (gave name '%s')", unique);
115                    }
116                    free(unique);
117                }
118
119                if (!error) {
120                    name = GB_read_char_pntr(gb_name);
121                    if (species_name_hash) {
122                        int seen = GBS_read_hash(species_name_hash, name);
123
124                        gb_assert(seen != 0); // species_name_hash not initialized correctly
125                        if (seen == 2) alignment_seen = 1; // already seen an alignment
126                    }
127                }
128
129                if (!error) {
130                    GB_push_my_security(gb_name);
131
132                    error             = GB_write_security_delete(gb_name, 7);
133                    if (!error) error = GB_write_security_write(gb_name, 6);
134
135                    if (!error) {
136                        gb_ali = GB_entry(gb_species, ali_name);
137                        if (gb_ali) {
138                            GBDATA *gb_data = GB_entry(gb_ali, "data");
139                            if (!gb_data) {
140                                error = GBT_write_string(gb_ali, "data", "Error: entry 'data' was missing and therefore was filled with this text.");
141                                GB_warningf("No '%s/data' entry for species '%s' (has been filled with dummy data)", ali_name, name);
142                            }
143                            else {
144                                if (GB_read_type(gb_data) != GB_STRING) {
145                                    GB_delete(gb_data);
146                                    error = GBS_global_string("'%s/data' of species '%s' had wrong DB-type (%s) and has been deleted!",
147                                                              ali_name, name, GB_read_key_pntr(gb_data));
148                                }
149                                else {
150                                    long data_len = GB_read_string_count(gb_data);
151                                    if (found_ali_len != data_len) {
152                                        if (found_ali_len>0)        aligned       = 0;
153                                        if (found_ali_len<data_len) found_ali_len = data_len;
154                                    }
155
156                                    error = GB_write_security_delete(gb_data, 7);
157
158                                    if (!alignment_seen && species_name_hash) { // mark as seen
159                                        GBS_write_hash(species_name_hash, name, 2); // 2 means "species has data in at least 1 alignment"
160                                    }
161                                }
162                            }
163                        }
164                    }
165
166                    if (!error) error = GB_write_security_delete(gb_species, security_write);
167
168                    GB_pop_my_security(gb_name);
169                }
170            }
171        }
172
173        if (!error) {
174            GBDATA *gb_sai;
175            for (gb_sai = GBT_first_SAI_rel_SAI_data(gb_extended_data);
176                 gb_sai && !error;
177                 gb_sai = GBT_next_SAI(gb_sai))
178            {
179                GBDATA *gb_sai_name = GB_entry(gb_sai, "name");
180                GBDATA *gb_ali;
181
182                if (!gb_sai_name) continue;
183
184                GB_write_security_delete(gb_sai_name, 7);
185
186                gb_ali = GB_entry(gb_sai, ali_name);
187                if (gb_ali) {
188                    GBDATA *gb_sai_data;
189                    for (gb_sai_data = GB_child(gb_ali);
190                         gb_sai_data;
191                         gb_sai_data = GB_nextChild(gb_sai_data))
192                    {
193                        long type = GB_read_type(gb_sai_data);
194                        long data_len;
195
196                        if (type == GB_DB || type < GB_BITS) continue;
197                        if (GB_read_key_pntr(gb_sai_data)[0] == '_') continue; // e.g. _STRUCT (of secondary structure)
198
199                        data_len = GB_read_count(gb_sai_data);
200
201                        if (found_ali_len != data_len) {
202                            if (found_ali_len>0)        aligned       = 0;
203                            if (found_ali_len<data_len) found_ali_len = data_len;
204                        }
205                    }
206                }
207            }
208        }
209
210        if (!error && stored_ali_len != found_ali_len) error = GB_write_int(gb_ali_len, found_ali_len);
211        if (!error) error = GBT_write_int(preset_alignment, "aligned", aligned);
212
213        if (error) {
214            error = GBS_global_string("Error checking alignment '%s':\n%s\n",  ali_name, error);
215        }
216    }
217
218    free(ali_name);
219
220    return error;
221}
222
223GB_ERROR GBT_check_data(GBDATA *Main, const char *alignment_name) {
224    /* alignment_name
225     *  == 0     -> check all existing alignments
226     * otherwise -> check only one alignment
227     */
228    GB_ERROR  error             = 0;
229    GBDATA   *gb_sd             = GBT_get_species_data(Main);
230    GBDATA   *gb_presets        = GBT_get_presets(Main);
231    GB_HASH  *species_name_hash = 0;
232
233    // create rest of main containers
234    GBT_get_SAI_data(Main);
235    GBT_get_tree_data(Main);
236
237    if (alignment_name) {
238        GBDATA *gb_ali_name = GB_find_string(gb_presets, "alignment_name", alignment_name, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
239        if (!gb_ali_name) {
240            error = GBS_global_string("Alignment '%s' does not exist - it can't be checked.", alignment_name);
241        }
242    }
243
244    if (!error) {
245        // check whether we have an default alignment
246        GBDATA *gb_use = GB_entry(gb_presets, "use");
247        if (!gb_use) {
248            // if we have no default alignment -> look for any alignment
249            GBDATA *gb_ali_name = GB_find_string(gb_presets, "alignment_name", alignment_name, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
250
251            error = gb_ali_name ?
252                GBT_write_string(gb_presets, "use", GB_read_char_pntr(gb_ali_name)) :
253                "No alignment defined";
254        }
255    }
256
257    if (!alignment_name && !error) {
258        // if all alignments are checked -> use species_name_hash to detect duplicated species and species w/o data
259        GBDATA *gb_species;
260        long    duplicates = 0;
261        species_name_hash  = GBS_create_hash(GBT_get_species_count(Main), GB_IGNORE_CASE);
262
263        if (!error) {
264            for (gb_species = GBT_first_species_rel_species_data(gb_sd);
265                 gb_species;
266                 gb_species = GBT_next_species(gb_species))
267            {
268                const char *name = GBT_read_name(gb_species);
269
270                if (GBS_read_hash(species_name_hash, name)) duplicates++;
271                GBS_incr_hash(species_name_hash, name);
272            }
273        }
274
275        if (duplicates) {
276            error = GBS_global_string("Database is corrupted:\n"
277                                      "Found %li duplicated species with identical names!\n"
278                                      "Fix the problem using\n"
279                                      "   'Search For Equal Fields and Mark Duplicates'\n"
280                                      "in ARB_NTREE search tool, save DB and restart ARB."
281                                      , duplicates);
282        }
283    }
284
285    if (!error) {
286        GBDATA *gb_ali;
287
288        for (gb_ali = GB_entry(gb_presets, "alignment");
289             gb_ali && !error;
290             gb_ali = GB_nextEntry(gb_ali))
291        {
292            error = GBT_check_alignment(Main, gb_ali, species_name_hash);
293        }
294    }
295
296    if (species_name_hash) {
297        if (!error) {
298            long counter = 0;
299            GBS_hash_do_loop(species_name_hash, check_for_species_without_data, &counter);
300            if (counter>0) {
301                GB_warningf("Found %li species without alignment data (only some were listed)", counter);
302            }
303        }
304
305        GBS_free_hash(species_name_hash);
306    }
307
308    return error;
309}
310
311void GBT_get_alignment_names(ConstStrArray& names, GBDATA *gbd) {
312    /* Get names of existing alignments from database.
313     *
314     * Returns: array of strings, the last element is NULL
315     * (Note: use GBT_free_names() to free result)
316     */
317
318    GBDATA *presets = GBT_get_presets(gbd);
319    for (GBDATA *ali = GB_entry(presets, "alignment"); ali; ali = GB_nextEntry(ali)) {
320        GBDATA *name = GB_entry(ali, "alignment_name");
321        names.put(name ? GB_read_char_pntr(name) : "<unnamed alignment>");
322    }
323}
324
325static char *gbt_nonexisting_alignment(GBDATA *gbMain) {
326    char  *ali_other = 0;
327    int    counter;
328
329    for (counter = 1; !ali_other; ++counter) {
330        ali_other = GBS_global_string_copy("ali_x%i", counter);
331        if (GBT_get_alignment(gbMain, ali_other) != 0) freenull(ali_other); // exists -> continue
332    }
333
334    return ali_other;
335}
336
337GB_ERROR GBT_check_alignment_name(const char *alignment_name)
338{
339    GB_ERROR error;
340    if ((error = GB_check_key(alignment_name))) return error;
341    if (strncmp(alignment_name, "ali_", 4)) {
342        return GB_export_errorf("your alignment_name '%s' must start with 'ali_'",
343                                alignment_name);
344    }
345    return 0;
346}
347
348static GB_ERROR create_ali_strEntry(GBDATA *gb_ali, const char *field, const char *strval, long write_protection) {
349    GB_ERROR  error  = 0;
350    GBDATA   *gb_sub = GB_create(gb_ali, field, GB_STRING);
351
352    if (!gb_sub) error = GB_await_error();
353    else {
354        error             = GB_write_string(gb_sub, strval);
355        if (!error) error = GB_write_security_delete(gb_sub, 7);
356        if (!error) error = GB_write_security_write(gb_sub, write_protection);
357    }
358
359    if (error) {
360        error = GBS_global_string("failed to create alignment subentry '%s'\n"
361                                  "(Reason: %s)", field, error);
362    }
363
364    return error;
365}
366static GB_ERROR create_ali_intEntry(GBDATA *gb_ali, const char *field, int intval, long write_protection) {
367    GB_ERROR  error  = 0;
368    GBDATA   *gb_sub = GB_create(gb_ali, field, GB_INT);
369
370    if (!gb_sub) error = GB_await_error();
371    else {
372        error             = GB_write_int(gb_sub, intval);
373        if (!error) error = GB_write_security_delete(gb_sub, 7);
374        if (!error) error = GB_write_security_write(gb_sub, write_protection);
375    }
376
377    if (error) {
378        error = GBS_global_string("failed to create alignment subentry '%s'\n"
379                                  "(Reason: %s)", field, error);
380    }
381
382    return error;
383}
384
385GBDATA *GBT_create_alignment(GBDATA *gbd, const char *name, long len, long aligned, long security, const char *type) {
386    /* create alignment
387     *
388     * returns pointer to alignment or
389     * NULL (in this case an error has been exported)
390     */
391    GB_ERROR  error      = NULL;
392    GBDATA   *gb_presets = GBT_get_presets(gbd);
393    GBDATA   *result     = NULL;
394
395    if (!gb_presets) {
396        error = GBS_global_string("can't find/create 'presets' (Reason: %s)", GB_await_error());
397    }
398    else {
399        error = GBT_check_alignment_name(name);
400        if (!error && (security<0 || security>6)) {
401            error = GBS_global_string("Illegal security value %li (allowed 0..6)", security);
402        }
403        if (!error) {
404            const char *allowed_types = ":dna:rna:ami:usr:";
405            int         tlen          = strlen(type);
406            const char *found         = strstr(allowed_types, type);
407            if (!found || found == allowed_types || found[-1] != ':' || found[tlen] != ':') {
408                error = GBS_global_string("Invalid alignment type '%s'", type);
409            }
410        }
411
412        if (!error) {
413            GBDATA *gb_name = GB_find_string(gb_presets, "alignment_name", name, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
414
415            if (gb_name) error = GBS_global_string("Alignment '%s' already exists", name);
416            else {
417                GBDATA *gb_ali     = GB_create_container(gb_presets, "alignment");
418                if (!gb_ali) error = GB_await_error();
419                else {
420                    error = GB_write_security_delete(gb_ali, 6);
421                    if (!error) error = create_ali_strEntry(gb_ali, "alignment_name",           name, 6);
422                    if (!error) error = create_ali_intEntry(gb_ali, "alignment_len",            len, 0);
423                    if (!error) error = create_ali_intEntry(gb_ali, "aligned",                  aligned <= 0 ? 0 : 1, 0);
424                    if (!error) error = create_ali_intEntry(gb_ali, "alignment_write_security", security, 6);
425                    if (!error) error = create_ali_strEntry(gb_ali, "alignment_type",           type, 0);
426                }
427
428                if (!error) result = gb_ali;
429            }
430        }
431    }
432
433    if (!result) {
434        gb_assert(error);
435        GB_export_errorf("in GBT_create_alignment: %s", error);
436    }
437#if defined(DEBUG)
438    else gb_assert(!error);
439#endif // DEBUG
440
441    return result;
442}
443
444
445static GB_ERROR gbt_rename_alignment_of_item(GBDATA *gb_item_container, const char *item_name, const char *item_entry_name,
446                                             const char *source, const char *dest, int copy, int dele)
447{
448    GB_ERROR  error = 0;
449    GBDATA   *gb_item;
450
451    for (gb_item = GB_entry(gb_item_container, item_entry_name);
452         gb_item && !error;
453         gb_item = GB_nextEntry(gb_item))
454    {
455        GBDATA *gb_ali = GB_entry(gb_item, source);
456        if (!gb_ali) continue;
457
458        if (copy) {
459            GBDATA *gb_new = GB_entry(gb_item, dest);
460            if (gb_new) {
461                error = GBS_global_string("Entry '%s' already exists", dest);
462            }
463            else {
464                gb_new             = GB_create_container(gb_item, dest);
465                if (!gb_new) error = GB_await_error();
466                else error         = GB_copy(gb_new, gb_ali);
467            }
468        }
469        if (dele) error = GB_delete(gb_ali);
470    }
471
472    if (error && gb_item) {
473        error = GBS_global_string("%s\n(while renaming alignment for %s '%s')", error, item_name, GBT_read_name(gb_item));
474    }
475
476    return error;
477}
478
479GB_ERROR GBT_rename_alignment(GBDATA *gbMain, const char *source, const char *dest, int copy, int dele)
480{
481    /* if copy == 1 then create a copy
482     * if dele == 1 then delete old
483     */
484
485    GB_ERROR  error            = 0;
486    int       is_case_error    = 0;
487    GBDATA   *gb_presets       = GBT_get_presets(gbMain);
488    GBDATA   *gb_species_data  = GBT_get_species_data(gbMain);
489    GBDATA   *gb_extended_data = GBT_get_SAI_data(gbMain);
490
491    if (!gb_presets || !gb_species_data || !gb_extended_data) error = GB_await_error();
492
493    // create copy and/or delete old alignment description
494    if (!error) {
495        GBDATA *gb_old_alignment = GBT_get_alignment(gbMain, source);
496
497        if (!gb_old_alignment) {
498            error = GB_await_error();
499        }
500        else {
501            if (copy) {
502                GBDATA *gbh = GBT_get_alignment(gbMain, dest);
503                if (gbh) {
504                    error         = GBS_global_string("destination alignment '%s' already exists", dest);
505                    is_case_error = (strcasecmp(source, dest) == 0); // test for case-only difference
506                }
507                else {
508                    GB_clear_error();
509                    error = GBT_check_alignment_name(dest);
510                    if (!error) {
511                        GBDATA *gb_new_alignment = GB_create_container(gb_presets, "alignment");
512                        error                    = GB_copy(gb_new_alignment, gb_old_alignment);
513                        if (!error) error        = GBT_write_string(gb_new_alignment, "alignment_name", dest);
514                    }
515                }
516            }
517
518            if (dele && !error) {
519                error = GB_delete(gb_old_alignment);
520            }
521        }
522    }
523
524    // change default alignment
525    if (!error && dele && copy) {
526        error = GBT_write_string(gb_presets, "use", dest);
527    }
528
529    // copy and/or delete alignment entries in species
530    if (!error) {
531        error = gbt_rename_alignment_of_item(gb_species_data, "Species", "species", source, dest, copy, dele);
532    }
533
534    // copy and/or delete alignment entries in SAIs
535    if (!error) {
536        error = gbt_rename_alignment_of_item(gb_extended_data, "SAI", "extended", source, dest, copy, dele);
537    }
538
539    if (is_case_error) {
540        // alignments source and dest only differ in case
541        char *ali_other = gbt_nonexisting_alignment(gbMain);
542        gb_assert(copy);
543
544        printf("Renaming alignment '%s' -> '%s' -> '%s' (to avoid case-problem)\n", source, ali_other, dest);
545
546        error             = GBT_rename_alignment(gbMain, source, ali_other, 1, dele);
547        if (!error) error = GBT_rename_alignment(gbMain, ali_other, dest, 1, 1);
548
549        free(ali_other);
550    }
551
552    return error;
553}
554
555// -----------------------------------------
556//      alignment related item functions
557
558NOT4PERL GBDATA *GBT_add_data(GBDATA *species, const char *ali_name, const char *key, GB_TYPES type) {
559    // goes to header: __ATTR__DEPRECATED_TODO("better use GBT_create_sequence_data()")
560
561    /* replace this function by GBT_create_sequence_data
562     * the same as GB_search(species, 'ali_name/key', GB_CREATE)
563     *
564     * Note: The behavior is weird, cause it does sth special for GB_STRING (write default content "...")
565     */
566
567    GBDATA *gb_gb;
568    GBDATA *gb_data;
569    if (GB_check_key(ali_name)) {
570        return NULL;
571    }
572    if (GB_check_hkey(key)) {
573        return NULL;
574    }
575    gb_gb = GB_entry(species, ali_name);
576    if (!gb_gb) gb_gb = GB_create_container(species, ali_name);
577
578    if (type == GB_STRING) {
579        gb_data = GB_search(gb_gb, key, GB_FIND);
580        if (!gb_data) {
581            gb_data = GB_search(gb_gb, key, GB_STRING);
582            GB_write_string(gb_data, "...");
583        }
584    }
585    else {
586        gb_data = GB_search(gb_gb, key, type);
587    }
588    return gb_data;
589}
590
591NOT4PERL GBDATA *GBT_create_sequence_data(GBDATA *species, const char *ali_name, const char *key, GB_TYPES type, int security_write) {
592    GBDATA *gb_data = GBT_add_data(species, ali_name, key, type);
593    if (gb_data) {
594        GB_ERROR error = GB_write_security_write(gb_data, security_write);
595        if (error) {
596            GB_export_error(error);
597            gb_data = 0;
598        }
599    }
600    return gb_data;
601}
602
603
604GB_ERROR GBT_write_sequence(GBDATA *gb_data, const char *ali_name, long ali_len, const char *sequence) {
605    /* writes a sequence which is generated by GBT_add_data,
606     * cuts sequence after alignment len only if bases e ".-nN"
607     */
608    int slen = strlen(sequence);
609    int old_char = 0;
610    GB_ERROR error = 0;
611    if (slen > ali_len) {
612        int i;
613        for (i = slen -1; i>=ali_len; i--) {
614            if (!strchr("-.nN", sequence[i])) break;    // real base after end of alignment
615        }
616        i++;                            // points to first 0 after alignment
617        if (i > ali_len) {
618            GBDATA *gb_main = GB_get_root(gb_data);
619            ali_len = GBT_get_alignment_len(gb_main, ali_name);
620            if (slen > ali_len) {               // check for modified alignment len
621                GBT_set_alignment_len(gb_main, ali_name, i);
622                ali_len = i;
623            }
624        }
625        if (slen > ali_len) {
626            old_char = sequence[ali_len];
627            ((char*)sequence)[ali_len] = 0;
628        }
629    }
630    error = GB_write_string(gb_data, sequence);
631    if (slen> ali_len) ((char*)sequence)[ali_len] = old_char;
632    return error;
633}
634
635
636GBDATA *GBT_gen_accession_number(GBDATA *gb_species, const char *ali_name) {
637    GBDATA *gb_acc = GB_entry(gb_species, "acc");
638    if (!gb_acc) {
639        GBDATA *gb_data = GBT_read_sequence(gb_species, ali_name);
640        if (gb_data) {                                     // found a valid alignment
641            GB_CSTR     sequence = GB_read_char_pntr(gb_data);
642            long        id       = GBS_checksum(sequence, 1, ".-");
643            const char *acc      = GBS_global_string("ARB_%lX", id);
644            GB_ERROR    error    = GBT_write_string(gb_species, "acc", acc);
645
646            if (error) GB_export_error(error);
647        }
648    }
649    return gb_acc;
650}
651
652
653int GBT_is_partial(GBDATA *gb_species, int default_value, int define_if_undef) {
654    // checks whether a species has a partial or full sequence
655    //
656    // Note: partial sequences should not be used for tree calculations
657    //
658    // returns: 0 if sequence is full
659    //          1 if sequence is partial
660    //          -1 in case of error (which is exported in this case)
661    //
662    // if the sequence has no 'ARB_partial' entry it returns 'default_value'
663    // if 'define_if_undef' is true then create an 'ARB_partial'-entry with the default value
664
665    int       result     = -1;
666    GB_ERROR  error      = 0;
667    GBDATA   *gb_partial = GB_entry(gb_species, "ARB_partial");
668
669    if (gb_partial) {
670        result = GB_read_int(gb_partial);
671        if (result != 0 && result != 1) {
672            error = "Illegal value for 'ARB_partial' (only 1 or 0 allowed)";
673        }
674    }
675    else {
676        if (define_if_undef) {
677            error = GBT_write_int(gb_species, "ARB_partial", default_value);
678        }
679        result = default_value;
680    }
681
682    if (error) {
683        GB_export_error(error);
684        return -1;
685    }
686    return result;
687}
688
689#if defined(WARN_TODO)
690#warning rename GBT_read_sequence - it does not read the sequence itself
691#endif
692GBDATA *GBT_read_sequence(GBDATA *gb_species, const char *aliname) {
693    GBDATA *gb_ali = GB_entry(gb_species, aliname);
694    return gb_ali ? GB_entry(gb_ali, "data") : 0;
695}
696
697char *GBT_get_default_alignment(GBDATA *gb_main) {
698    return GBT_read_string(gb_main, "presets/use");
699}
700
701GB_ERROR GBT_set_default_alignment(GBDATA *gb_main, const char *alignment_name) {
702    return GBT_write_string(gb_main, "presets/use", alignment_name);
703}
704
705GBDATA *GBT_get_alignment(GBDATA *gb_main, const char *aliname) {
706    GBDATA *gb_presets        = GBT_get_presets(gb_main);
707    GBDATA *gb_alignment_name = GB_find_string(gb_presets, "alignment_name", aliname, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
708
709    if (!gb_alignment_name) {
710        GB_export_errorf("alignment '%s' not found", aliname);
711        return NULL;
712    }
713    return GB_get_father(gb_alignment_name);
714}
715
716#if defined(WARN_TODO)
717#warning recode and change result type to long* ?
718#endif
719long GBT_get_alignment_len(GBDATA *gb_main, const char *aliname) {
720    GBDATA *gb_alignment = GBT_get_alignment(gb_main, aliname);
721    return gb_alignment ? *GBT_read_int(gb_alignment, "alignment_len") : -1;
722}
723
724GB_ERROR GBT_set_alignment_len(GBDATA *gb_main, const char *aliname, long new_len) {
725    GB_ERROR  error        = 0;
726    GBDATA   *gb_alignment = GBT_get_alignment(gb_main, aliname);
727
728    if (gb_alignment) {
729        GB_push_my_security(gb_main);
730        error             = GBT_write_int(gb_alignment, "alignment_len", new_len); // write new len
731        if (!error) error = GBT_write_int(gb_alignment, "aligned", 0);             // mark as unaligned
732        GB_pop_my_security(gb_main);
733    }
734    else error = GB_export_errorf("Alignment '%s' not found", aliname);
735    return error;
736}
737
738char *GBT_get_alignment_type_string(GBDATA *gb_main, const char *aliname) {
739    char   *result       = NULL;
740    GBDATA *gb_alignment = GBT_get_alignment(gb_main, aliname);
741    if (gb_alignment) {
742        result = GBT_read_string(gb_alignment, "alignment_type");
743        gb_assert(result);
744    }
745    return result;
746}
747
748GB_alignment_type GBT_get_alignment_type(GBDATA *gb_main, const char *aliname) {
749    char              *ali_type = GBT_get_alignment_type_string(gb_main, aliname);
750    GB_alignment_type  at       = GB_AT_UNKNOWN;
751
752    if (ali_type) {
753        switch (ali_type[0]) {
754            case 'r': if (strcmp(ali_type, "rna")==0) at = GB_AT_RNA; break;
755            case 'd': if (strcmp(ali_type, "dna")==0) at = GB_AT_DNA; break;
756            case 'a': if (strcmp(ali_type, "ami")==0) at = GB_AT_AA; break;
757            case 'p': if (strcmp(ali_type, "pro")==0) at = GB_AT_AA; break;
758            default: gb_assert(0); break;
759        }
760        free(ali_type);
761    }
762    return at;
763}
764
765bool GBT_is_alignment_protein(GBDATA *gb_main, const char *alignment_name) {
766    return GBT_get_alignment_type(gb_main, alignment_name) == GB_AT_AA;
767}
768
769// -----------------------
770//      gene sequence
771
772static const char *gb_cache_genome(GBDATA *gb_genome) {
773    static GBDATA *gb_last_genome = 0;
774    static char   *last_genome    = 0;
775
776    if (gb_genome != gb_last_genome) {
777        free(last_genome);
778
779        last_genome    = GB_read_string(gb_genome);
780        gb_last_genome = gb_genome;
781    }
782
783    return last_genome;
784}
785
786struct gene_part_pos {
787    int            parts;       // initialized for parts
788    unsigned char *certain;     // contains parts "=" chars
789    char           offset[256];
790};
791
792static gene_part_pos *gpp = 0;
793
794static void init_gpp(int parts) {
795    if (!gpp) {
796        int i;
797        gpp          = (gene_part_pos*)malloc(sizeof(*gpp));
798        gpp->certain = 0;
799
800        for (i = 0; i<256; ++i) gpp->offset[i] = 0;
801
802        gpp->offset[(int)'+'] = 1;
803        gpp->offset[(int)'-'] = -1;
804    }
805    else {
806        if (parts>gpp->parts) freenull(gpp->certain);
807    }
808
809    if (!gpp->certain) {
810        int forParts           = parts+10;
811        gpp->certain           = (unsigned char *)malloc(forParts+1);
812        memset(gpp->certain, '=', forParts);
813        gpp->certain[forParts] = 0;
814        gpp->parts             = forParts;
815    }
816}
817
818static void getPartPositions(const GEN_position *pos, int part, size_t *startPos, size_t *stopPos) {
819    // returns 'startPos' and 'stopPos' of one part of a gene
820    gb_assert(part<pos->parts);
821    *startPos = pos->start_pos[part]+gpp->offset[(pos->start_uncertain ? pos->start_uncertain : gpp->certain)[part]];
822    *stopPos  = pos->stop_pos [part]+gpp->offset[(pos->stop_uncertain  ? pos->stop_uncertain  : gpp->certain)[part]];
823}
824
825NOT4PERL char *GBT_read_gene_sequence_and_length(GBDATA *gb_gene, bool use_revComplement, char partSeparator, size_t *gene_length) {
826    // return the sequence data of a gene
827    //
828    // if use_revComplement is true -> use data from complementary strand (if complement is set for gene)
829    //                    otherwise -> use data from primary strand (sort+merge parts by position)
830    //
831    // if partSeparator not is 0 -> insert partSeparator between single (non-merged) parts
832    //
833    // returns sequence as result (and length of sequence if 'gene_length' points to something)
834    //
835    // if 'pos_certain' contains '+', start behind position (or end at position)
836    //                           '-', start at position (or end before position)
837    //
838    // For zero-length genes (e.g. "711^712") this function returns an empty string.
839
840    GB_ERROR      error      = 0;
841    char         *result     = 0;
842    GBDATA       *gb_species = GB_get_grandfather(gb_gene);
843    GEN_position *pos        = GEN_read_position(gb_gene);
844
845    if (!pos) error = GB_await_error();
846    else {
847        GBDATA        *gb_seq        = GBT_read_sequence(gb_species, "ali_genom");
848        unsigned long  seq_length    = GB_read_count(gb_seq);
849        int            p;
850        int            parts         = pos->parts;
851        int            resultlen     = 0;
852        int            separatorSize = partSeparator ? 1 : 0;
853
854        init_gpp(parts);
855
856        // test positions and calculate overall result length
857        for (p = 0; p<parts && !error; p++) {
858            size_t start;
859            size_t stop;
860            getPartPositions(pos, p, &start, &stop);
861
862            if (start<1 || start>(stop+1) || stop > seq_length) { // do not reject zero-length genes (start == stop+1)
863                error = GBS_global_string("Illegal gene position(s): start=%zu, end=%zu, seq.length=%li",
864                                          start, stop, seq_length);
865            }
866            else {
867                resultlen += stop-start+1;
868            }
869        }
870
871        if (separatorSize) resultlen += (parts-1)*separatorSize;
872
873        if (!error) {
874            char T_or_U = 0;
875            if (use_revComplement) {
876                error = GBT_determine_T_or_U(GB_AT_DNA, &T_or_U, "reverse-complement");
877            }
878            else if (parts>1) {
879                GEN_sortAndMergeLocationParts(pos);
880                parts = pos->parts; // may have changed
881            }
882
883            if (!error) {
884                const char *seq_data = gb_cache_genome(gb_seq);
885                char       *resultpos;
886
887                result    = (char*)malloc(resultlen+1);
888                resultpos = result;
889
890                if (gene_length) *gene_length = resultlen;
891
892                for (p = 0; p<parts; ++p) {
893                    size_t start;
894                    size_t stop;
895                    getPartPositions(pos, p, &start, &stop);
896
897                    int len = stop-start+1;
898
899                    if (separatorSize && p>0) *resultpos++ = partSeparator;
900
901                    memcpy(resultpos, seq_data+start-1, len);
902                    if (T_or_U && pos->complement[p]) {
903                        GBT_reverseComplementNucSequence(resultpos, len, T_or_U);
904                    }
905                    resultpos += len;
906                }
907
908                resultpos[0] = 0;
909            }
910        }
911        GEN_free_position(pos);
912    }
913
914    gb_assert(result || error);
915
916    if (error) {
917        char *id = GEN_global_gene_identifier(gb_gene, gb_species);
918        error    = GB_export_errorf("Can't read sequence of '%s' (Reason: %s)", id, error);
919        free(id);
920        free(result);
921        result   = 0;
922    }
923
924    return result;
925}
926
927char *GBT_read_gene_sequence(GBDATA *gb_gene, bool use_revComplement, char partSeparator) {
928    return GBT_read_gene_sequence_and_length(gb_gene, use_revComplement, partSeparator, 0);
929}
930
931// --------------------------------------------------------------------------------
932
933#ifdef UNIT_TESTS
934#include <test_unit.h>
935
936void TEST_alignment() {
937    GB_shell  shell;
938    GBDATA   *gb_main = GB_open("TEST_prot.arb", "r");
939
940    {
941        GB_transaction ta(gb_main);
942
943        TEST_ASSERT_EQUAL(GBT_count_alignments(gb_main), 2);
944
945        char *def_ali_name = GBT_get_default_alignment(gb_main);
946        TEST_ASSERT_EQUAL(def_ali_name, "ali_tuf_dna");
947
948        {
949            ConstStrArray names;
950            GBT_get_alignment_names(names, gb_main);
951            {
952                char *joined = GBT_join_names(names, '*');
953                TEST_ASSERT_EQUAL(joined, "ali_tuf_pro*ali_tuf_dna");
954                free(joined);
955            }
956
957            for (int i = 0; names[i]; ++i) {
958                long len = GBT_get_alignment_len(gb_main, names[i]);
959                TEST_ASSERT_EQUAL(len, !i ? 473 : 1426);
960
961                char *type_name = GBT_get_alignment_type_string(gb_main, names[i]);
962                TEST_ASSERT_EQUAL(type_name, !i ? "ami" : "dna");
963                free(type_name);
964
965                GB_alignment_type type = GBT_get_alignment_type(gb_main, names[i]);
966                TEST_ASSERT_EQUAL(type, !i ? GB_AT_AA : GB_AT_DNA);
967                TEST_ASSERT(!i == GBT_is_alignment_protein(gb_main, names[i]));
968            }
969        }
970
971        free(def_ali_name);
972    }
973
974    GB_close(gb_main);
975}
976
977#endif // UNIT_TESTS
Note: See TracBrowser for help on using the browser.