source: tags/ms_r16q3/ARBDB/adali.cxx

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