source: branches/port5/ARBDB/adali.c

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