source: tags/ms_r17q3/ARBDB/adali.cxx

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