source: tags/svn.1.5.4/ARBDB/adali.cxx

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