root/trunk/ARBDB/adGene.cxx

Revision 8745, 34.4 KB (checked in by westram, 20 hours ago)
  • format string fixes (32bit compile complained)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : adGene.cxx                                        //
4//   Purpose   : Basic gene access functions                       //
5//                                                                 //
6//   Coded by Ralf Westram (coder@reallysoft.de) in July 2002      //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "gb_local.h"
13
14#include <adGene.h>
15#include <arbdbt.h>
16#include <arb_strbuf.h>
17#include <arb_strarray.h>
18
19
20bool GEN_is_genome_db(GBDATA *gb_main, int default_value) {
21    // default_value ==  0 -> default to normal database
22    //               ==  1 -> default to GENOM database
23    //               == -1 -> assume that type is already defined
24    GBDATA *gb_genom_db = GB_entry(gb_main, GENOM_DB_TYPE);
25
26    if (!gb_genom_db) {         // no DB-type entry -> create one with default
27        GB_ERROR error = NULL;
28
29        assert_or_exit(default_value != -1); // first call to GEN_is_genome_db has to provide a 'default_value'
30
31        gb_genom_db             = GB_create(gb_main, GENOM_DB_TYPE, GB_INT);
32        if (!gb_genom_db) error = GB_await_error();
33        else error              = GB_write_int(gb_genom_db, default_value);
34
35        if (error) GBK_terminatef("Fatal in GEN_is_genome_db: %s", error);
36    }
37
38    return GB_read_int(gb_genom_db) != 0;
39}
40
41//  --------------
42//      genes:
43
44GBDATA* GEN_findOrCreate_gene_data(GBDATA *gb_species) {
45    GBDATA *gb_gene_data = GB_search(gb_species, "gene_data", GB_CREATE_CONTAINER);
46    gb_assert(gb_gene_data);
47    return gb_gene_data;
48}
49
50GBDATA* GEN_find_gene_data(GBDATA *gb_species) {
51    return GB_search(gb_species, "gene_data", GB_FIND);
52}
53
54GBDATA* GEN_expect_gene_data(GBDATA *gb_species) {
55    GBDATA *gb_gene_data = GB_search(gb_species, "gene_data", GB_FIND);
56    gb_assert(gb_gene_data);
57    return gb_gene_data;
58}
59
60GBDATA* GEN_find_gene_rel_gene_data(GBDATA *gb_gene_data, const char *name) {
61    GBDATA *gb_name = GB_find_string(gb_gene_data, "name", name, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
62
63    if (gb_name) return GB_get_father(gb_name); // found existing gene
64    return 0;
65}
66
67GBDATA* GEN_find_gene(GBDATA *gb_species, const char *name) {
68    // find existing gene. returns 0 if it does not exist.
69    GBDATA *gb_gene_data = GEN_find_gene_data(gb_species);
70    return gb_gene_data ? GEN_find_gene_rel_gene_data(gb_gene_data, name) : 0;
71}
72
73static GBDATA* GEN_create_nonexisting_gene_rel_gene_data(GBDATA *gb_gene_data, const char *name) {
74    GB_ERROR  error   = GB_push_transaction(gb_gene_data);
75    GBDATA   *gb_gene = 0;
76
77    gb_assert(!GEN_find_gene_rel_gene_data(gb_gene_data, name)); // don't call this function if you are not sure that the gene does not exists!
78
79    if (!error) {
80        gb_gene = GB_create_container(gb_gene_data, "gene");
81        error   = gb_gene ? GBT_write_string(gb_gene, "name", name) : GB_await_error();
82    }
83
84    gb_assert(gb_gene || error);
85    error = GB_end_transaction(gb_gene_data, error);
86    if (error) GB_export_error(error);
87
88    return gb_gene;
89}
90
91GBDATA* GEN_create_nonexisting_gene(GBDATA *gb_species, const char *name) { // needed by ../PERL_SCRIPTS/GENOME/GI.pm@create_nonexisting_gene
92    return GEN_create_nonexisting_gene_rel_gene_data(GEN_findOrCreate_gene_data(gb_species), name); 
93} 
94
95GBDATA* GEN_find_or_create_gene_rel_gene_data(GBDATA *gb_gene_data, const char *name) {
96    GBDATA *gb_gene = 0;
97
98    // Search for a gene, when gene does not exist create it
99    if (!name || !name[0]) {
100        GB_export_error("Missing gene name");
101    }
102    else {
103        GBDATA *gb_name = GB_find_string(gb_gene_data, "name", name, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
104
105        if (gb_name) {
106            gb_gene = GB_get_father(gb_name); // found existing gene
107        }
108        else {
109            GB_ERROR error = GB_push_transaction(gb_gene_data);
110
111            if (!error) {
112                gb_gene = GB_create_container(gb_gene_data, "gene");
113                error   = GBT_write_string(gb_gene, "name", name);
114            }
115            error = GB_end_transaction(gb_gene_data, error);
116            if (error) {
117                gb_gene = NULL;
118                GB_export_error(error);
119            }
120        }
121    }
122    return gb_gene;
123}
124
125GBDATA* GEN_first_gene(GBDATA *gb_species) {
126    return GB_entry(GEN_expect_gene_data(gb_species), "gene");
127}
128
129GBDATA* GEN_first_gene_rel_gene_data(GBDATA *gb_gene_data) {
130    return GB_entry(gb_gene_data, "gene");
131}
132
133GBDATA* GEN_next_gene(GBDATA *gb_gene) {
134    gb_assert(GB_has_key(gb_gene, "gene"));
135    return GB_nextEntry(gb_gene);
136}
137
138GBDATA *GEN_first_marked_gene(GBDATA *gb_species) {
139    return GB_first_marked(GEN_expect_gene_data(gb_species), "gene");
140}
141GBDATA *GEN_next_marked_gene(GBDATA *gb_gene) {
142    return GB_next_marked(gb_gene, "gene");
143}
144
145// ----------------------
146//      gene position
147
148static GEN_position *lastFreedPosition = 0;
149
150GEN_position *GEN_new_position(int parts, bool joinable) {
151    GEN_position *pos;
152
153    size_t pos_size  = parts*sizeof(pos->start_pos[0]);
154    size_t comp_size = parts*sizeof(pos->complement[0]);
155    size_t data_size = 2*pos_size+3*comp_size;
156
157    gb_assert(parts>0);
158
159    if (lastFreedPosition && lastFreedPosition->parts == parts) {
160        pos               = lastFreedPosition;
161        lastFreedPosition = 0;
162        memset(pos->start_pos, 0, data_size);
163    }
164    else {
165        pos             = (GEN_position*)GB_calloc(1, sizeof(*pos));
166        pos->parts      = parts;
167        pos->start_pos  = (size_t*)GB_calloc(1, data_size);
168        pos->stop_pos   = pos->start_pos+parts;
169        pos->complement = (unsigned char*)(pos->stop_pos+parts);
170    }
171
172    pos->joinable        = joinable;
173    pos->start_uncertain = 0;
174    pos->stop_uncertain  = 0;
175
176    return pos;
177}
178
179void GEN_use_uncertainties(GEN_position *pos) {
180    if (pos->start_uncertain == 0) {
181        // space was already allocated in GEN_new_position
182        pos->start_uncertain = pos->complement+pos->parts;
183        pos->stop_uncertain  = pos->start_uncertain+pos->parts;
184
185        size_t comp_size = pos->parts*sizeof(pos->complement[0]);
186        memset(pos->start_uncertain, '=', 2*comp_size);
187    }
188}
189
190void GEN_free_position(GEN_position *pos) {
191    if (pos) {
192        if (lastFreedPosition) {
193            free(lastFreedPosition->start_pos); // rest is allocated together with start_pos
194            free(lastFreedPosition);
195        }
196
197        lastFreedPosition = pos;
198    }
199}
200
201static struct GEN_position_mem_handler {
202    ~GEN_position_mem_handler() { GEN_free_position(NULL); }
203} GEN_position_dealloc;
204
205
206static GB_ERROR parseCSV(GBDATA *gb_gene, const char *field_name, size_t parts_expected, ConstStrArray& parseTable) {
207    // reads a field and splits the content at ','
208    // results are stored in parseTable
209
210    GB_ERROR  error      = 0;
211    GBDATA   *gb_field   = GB_entry(gb_gene, field_name);
212    if (!gb_field) error = GBS_global_string("Expected entry '%s' missing", field_name);
213    else {
214        char *content       = GB_read_string(gb_field);
215        if (!content) error = GB_await_error();
216        else {
217            parseTable.erase();
218            GBT_splitNdestroy_string(parseTable, content, ',');
219            if (parseTable.size() != parts_expected) {
220                error = GBS_global_string("Expected %zu CSV, found %zu", parts_expected, parseTable.size());
221            }
222        }
223    }
224    return error;
225}
226
227static GB_ERROR parsePositions(GBDATA *gb_gene, const char *field_name, int parts_expected, size_t *results, ConstStrArray& parseTable) {
228    GB_ERROR error = parseCSV(gb_gene, field_name, parts_expected, parseTable);
229    if (!error) {
230        int p;
231        for (p = 0; p<parts_expected && !error; p++) {
232            char *end;
233            results[p] = strtol(parseTable[p], &end, 10);
234            if (end == parseTable[p]) { // error
235                error = GBS_global_string("can't convert '%s' to number", parseTable[p]);
236            }
237        }
238    }
239    if (error) {
240        error = GBS_global_string("While parsing field '%s': %s", field_name, error);
241    }
242    return error;
243}
244
245GEN_position *GEN_read_position(GBDATA *gb_gene) {
246    int           parts         = 1;
247    bool          joinable      = false;
248    GBDATA       *gb_pos_joined = GB_entry(gb_gene, "pos_joined");
249    GEN_position *pos           = 0;
250    GB_ERROR      error         = 0;
251
252    if (gb_pos_joined) {
253        parts = GB_read_int(gb_pos_joined);
254        if (parts != 1) { // splitted
255            if (parts>1) joinable = true;
256            else if (parts<-1) parts = -parts; // neg value means "not joinable" (comes from feature location 'order(...)')
257            else error = GBS_global_string("Illegal value %i in 'pos_joined'", parts);
258        }
259    }
260
261    if (!error) {
262        pos = GEN_new_position(parts, joinable);
263
264        ConstStrArray parseTable;
265        parseTable.reserve(parts);
266
267        error =             parsePositions(gb_gene, "pos_start", parts, pos->start_pos, parseTable);
268        if (!error) error = parsePositions(gb_gene, "pos_stop",  parts, pos->stop_pos,  parseTable);
269
270        int p;
271        if (!error) {
272            error = parseCSV(gb_gene, "pos_complement",  parts, parseTable);
273            for (p = 0; p<parts && !error; p++) {
274                const char *val = parseTable[p];
275                if ((val[0] != '0' && val[0] != '1') || val[1] != 0) {
276                    error = GBS_global_string("Invalid content '%s' in 'pos_complement' (expected: \"01\")", val);
277                }
278                else {
279                    pos->complement[p] = (unsigned char)atoi(val);
280                }
281            }
282        }
283
284        if (!error) {
285            GBDATA *gb_pos_certain = GB_entry(gb_gene, "pos_certain");
286
287            if (gb_pos_certain) {
288                error = parseCSV(gb_gene, "pos_certain",  parts, parseTable);
289                GEN_use_uncertainties(pos);
290                for (p = 0; p<parts && !error; p++) {
291                    const unsigned char *val = (unsigned char *)(parseTable[p]);
292                    int                  vp;
293
294                    for (vp = 0; vp<2; vp++) {
295                        unsigned char c = val[vp];
296                        if (c != '<' && c != '=' && c != '>' && (c != "+-"[vp])) {
297                            error = GBS_global_string("Invalid content '%s' in 'pos_certain' (expected 2 from \"<=>\")", val);
298                        }
299                    }
300                    if (!error) {
301                        pos->start_uncertain[p] = val[0];
302                        pos->stop_uncertain[p]  = val[1];
303                    }
304                }
305            }
306        }
307    }
308
309    gb_assert(error || pos);
310    if (error) {
311        GB_export_error(error);
312        if (pos) {
313            GEN_free_position(pos);
314            pos = 0;
315        }
316    }
317    return pos;
318}
319
320GB_ERROR GEN_write_position(GBDATA *gb_gene, const GEN_position *pos) {
321    GB_ERROR  error          = 0;
322    GBDATA   *gb_pos_joined  = GB_entry(gb_gene, "pos_joined");
323    GBDATA   *gb_pos_certain = GB_entry(gb_gene, "pos_certain");
324    GBDATA   *gb_pos_start;
325    GBDATA   *gb_pos_stop;
326    GBDATA   *gb_pos_complement;
327    int       p;
328
329    gb_assert(pos);
330
331    gb_pos_start             = GB_search(gb_gene, "pos_start", GB_STRING);
332    if (!gb_pos_start) error = GB_await_error();
333
334    if (!error) {
335        gb_pos_stop             = GB_search(gb_gene, "pos_stop", GB_STRING);
336        if (!gb_pos_stop) error = GB_await_error();
337    }
338    if (!error) {
339        gb_pos_complement             = GB_search(gb_gene, "pos_complement", GB_STRING);
340        if (!gb_pos_complement) error = GB_await_error();
341    }
342
343    if (!error) {
344        if (pos->start_uncertain) {
345            if (!gb_pos_certain) {
346                gb_pos_certain             = GB_search(gb_gene, "pos_certain", GB_STRING);
347                if (!gb_pos_certain) error = GB_await_error();
348            }
349        }
350        else {
351            if (gb_pos_certain) {
352                error          = GB_delete(gb_pos_certain);
353                gb_pos_certain = 0;
354            }
355        }
356    }
357
358    // test data
359    if (!error) {
360        size_t length;
361        {
362            GBDATA *gb_organism = GB_get_grandfather(gb_gene);
363            GBDATA *gb_genome   = GBT_read_sequence(gb_organism, GENOM_ALIGNMENT);
364           
365            length = GB_read_count(gb_genome);
366        }
367
368        for (p = 0; p<pos->parts && !error; ++p) {
369            char c;
370
371            c = pos->complement[p]; gb_assert(c == 0 || c == 1);
372            if (c<0 || c>1) {
373                error = GBS_global_string("Illegal value %i in complement", int(c));
374            }
375            else {
376                if (pos->start_pos[p]>pos->stop_pos[p]) {
377                    error = GBS_global_string("Illegal positions (%zu>%zu)", pos->start_pos[p], pos->stop_pos[p]);
378                }
379                else if (pos->start_pos[p] == 0) {
380                    error = GBS_global_string("Illegal start position %zu", pos->start_pos[p]);
381                }
382                else if (pos->stop_pos[p] > length) {
383                    error = GBS_global_string("Illegal stop position %zu (>length(=%zu))", pos->stop_pos[p], length);
384                }
385                else {
386                    if (pos->start_uncertain) {
387                        c       = pos->start_uncertain[p];
388                        char c2 = pos->stop_uncertain[p];
389
390                        if      (!c  || strchr("<=>+", c)  == 0) error = GBS_global_string("Invalid uncertainty '%c'", c);
391                        else if (!c2 || strchr("<=>-", c2) == 0) error = GBS_global_string("Invalid uncertainty '%c'", c2);
392                        else {
393                            if (c == '+' || c2 == '-') {
394                                if (c == '+' && c2 == '-') {
395                                    if (pos->start_pos[p] != pos->stop_pos[p]-1) {
396                                        error = GBS_global_string("Invalid positions %zu^%zu for uncertainties +-", pos->start_pos[p], pos->stop_pos[p]);
397                                    }
398                                }
399                                else {
400                                    error = "uncertainties '+' and '-' can only be used together";
401                                }
402                            }
403                        }
404                    }
405                }
406            }
407        }
408    }
409
410    if (!error) {
411        if (pos->parts == 1) {
412            if (gb_pos_joined) error = GB_delete(gb_pos_joined);
413
414            if (!error) error = GB_write_string(gb_pos_start,      GBS_global_string("%zu", pos->start_pos[0]));
415            if (!error) error = GB_write_string(gb_pos_stop,       GBS_global_string("%zu", pos->stop_pos[0]));
416            if (!error) error = GB_write_string(gb_pos_complement, GBS_global_string("%c", pos->complement[0]+'0'));
417
418            if (!error && gb_pos_certain) {
419                error = GB_write_string(gb_pos_certain, GBS_global_string("%c%c", pos->start_uncertain[0], pos->stop_uncertain[0]));
420            }
421        }
422        else {
423            if (!gb_pos_joined) {
424                gb_pos_joined             = GB_search(gb_gene, "pos_joined", GB_INT);
425                if (!gb_pos_joined) error = GB_await_error();
426            }
427            if (!error) error = GB_write_int(gb_pos_joined, pos->parts * (pos->joinable ? 1 : -1)); // neg. parts means not joinable
428
429            if (!error) {
430                GBS_strstruct *start      = GBS_stropen(12*pos->parts);
431                GBS_strstruct *stop       = GBS_stropen(12*pos->parts);
432                GBS_strstruct *complement = GBS_stropen(2*pos->parts);
433                GBS_strstruct *uncertain  = GBS_stropen(3*pos->parts);
434
435                for (p = 0; p<pos->parts; ++p) {
436                    if (p>0) {
437                        GBS_chrcat(start, ',');
438                        GBS_chrcat(stop, ',');
439                        GBS_chrcat(complement, ',');
440                        GBS_chrcat(uncertain, ',');
441                    }
442                    GBS_strcat(start, GBS_global_string("%zu", pos->start_pos[p]));
443                    GBS_strcat(stop,  GBS_global_string("%zu", pos->stop_pos[p]));
444                    GBS_chrcat(complement, pos->complement[p]+'0');
445                    if (gb_pos_certain) {
446                        GBS_chrcat(uncertain, pos->start_uncertain[p]);
447                        GBS_chrcat(uncertain, pos->stop_uncertain[p]);
448                    }
449                }
450
451                char *sstart      = GBS_strclose(start);
452                char *sstop       = GBS_strclose(stop);
453                char *scomplement = GBS_strclose(complement);
454                char *suncertain  = GBS_strclose(uncertain);
455
456                error             = GB_write_string(gb_pos_start, sstart);
457                if (!error) error = GB_write_string(gb_pos_stop, sstop);
458                if (!error) error = GB_write_string(gb_pos_complement, scomplement);
459                if (!error && gb_pos_certain) error = GB_write_string(gb_pos_certain, suncertain);
460
461                free(suncertain);
462                free(scomplement);
463                free(sstop);
464                free(sstart);
465            }
466        }
467    }
468
469    return error;
470}
471
472static GEN_position *location2sort = 0;
473
474static int cmp_location_parts(const void *v1, const void *v2) {
475    int i1 = *(int*)v1;
476    int i2 = *(int*)v2;
477
478    int cmp = location2sort->start_pos[i1]-location2sort->start_pos[i2];
479    if (!cmp) {
480        cmp = location2sort->stop_pos[i1]-location2sort->stop_pos[i2];
481    }
482    return cmp;
483}
484
485void GEN_sortAndMergeLocationParts(GEN_position *location) {
486    // Note: makes location partly invalid (only start_pos + stop_pos are valid afterwards)
487    int  parts = location->parts;
488    int *idx   = (int*)malloc(parts*sizeof(*idx)); // idx[newpos] = oldpos
489    int  i, p;
490
491    for (p = 0; p<parts; ++p) idx[p] = p;
492
493    location2sort = location;
494    qsort(idx, parts, sizeof(*idx), cmp_location_parts);
495    location2sort = 0;
496
497    for (p = 0; p<parts; ++p) {
498        i = idx[p];
499
500#define swap(a, b, type) do { type tmp = (a); (a) = (b); (b) = (tmp); } while (0)
501
502        if (i != p) {
503            swap(location->start_pos[i],  location->start_pos[p],  size_t);
504            swap(location->stop_pos[i],   location->stop_pos[p],   size_t);
505            swap(idx[i], idx[p], int);
506        }
507    }
508
509#if defined(DEBUG) && 0
510    printf("Locations sorted:\n");
511    for (p = 0; p<parts; ++p) {
512        printf("  [%i] %i - %i %i\n", p, location->start_pos[p], location->stop_pos[p], (int)(location->complement[p]));
513    }
514#endif // DEBUG
515
516    i = 0;
517    for (p = 1; p<parts; p++) {
518        if ((location->stop_pos[i]+1) >= location->start_pos[p]) {
519            // parts overlap or are directly consecutive
520
521            location->stop_pos[i]  = location->stop_pos[p];
522        }
523        else {
524            i++;
525            location->start_pos[i] = location->start_pos[p];
526            location->stop_pos[i]  = location->stop_pos[p];
527        }
528    }
529    location->parts = i+1;
530
531#if defined(DEBUG) && 0
532    parts = location->parts;
533    printf("Locations merged:\n");
534    for (p = 0; p<parts; ++p) {
535        printf("  [%i] %i - %i %i\n", p, location->start_pos[p], location->stop_pos[p], (int)(location->complement[p]));
536    }
537#endif // DEBUG
538
539    free(idx);
540}
541
542
543
544//  -----------------------------------------
545//      test if species is pseudo-species
546
547const char *GEN_origin_organism(GBDATA *gb_pseudo) {
548    GBDATA *gb_origin = GB_entry(gb_pseudo, "ARB_origin_species");
549    return gb_origin ? GB_read_char_pntr(gb_origin) : 0;
550}
551const char *GEN_origin_gene(GBDATA *gb_pseudo) {
552    GBDATA *gb_origin = GB_entry(gb_pseudo, "ARB_origin_gene");
553    return gb_origin ? GB_read_char_pntr(gb_origin) : 0;
554}
555
556bool GEN_is_pseudo_gene_species(GBDATA *gb_species) {
557    return GEN_origin_organism(gb_species) != 0;
558}
559
560//  ------------------------------------------------
561//      find organism or gene for pseudo-species
562
563GB_ERROR GEN_organism_not_found(GBDATA *gb_pseudo) {
564    GB_ERROR error = 0;
565
566    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
567    gb_assert(GEN_find_origin_organism(gb_pseudo, 0) == 0);
568
569    error = GB_export_errorf("The gene-species '%s' refers to an unknown organism (%s)\n"
570                             "This occurs if you rename or delete the organism or change the entry\n"
571                             "'ARB_origin_species' and will most likely cause serious problems.",
572                             GBT_read_name(gb_pseudo),
573                             GEN_origin_organism(gb_pseudo));
574
575    return error;
576}
577
578// @@@ FIXME: missing: GEN_gene_not_found (like GEN_organism_not_found)
579
580// ------------------------------
581//      search pseudo species
582
583static const char *pseudo_species_hash_key(const char *organism_name, const char *gene_name) {
584    return GBS_global_string("%s*%s", organism_name, gene_name);
585}
586
587static GBDATA *GEN_read_pseudo_species_from_hash(GB_HASH *pseudo_hash, const char *organism_name, const char *gene_name) {
588    return (GBDATA*)GBS_read_hash(pseudo_hash, pseudo_species_hash_key(organism_name, gene_name));
589}
590
591void GEN_add_pseudo_species_to_hash(GBDATA *gb_pseudo, GB_HASH *pseudo_hash) {
592    const char *organism_name = GEN_origin_organism(gb_pseudo);
593    const char *gene_name     = GEN_origin_gene(gb_pseudo);
594
595    gb_assert(organism_name);
596    gb_assert(gene_name);
597
598    GBS_write_hash(pseudo_hash, pseudo_species_hash_key(organism_name, gene_name), (long)gb_pseudo);
599}
600
601GB_HASH *GEN_create_pseudo_species_hash(GBDATA *gb_main, int additionalSize) {
602    GB_HASH *pseudo_hash = GBS_create_hash(GBT_get_species_count(gb_main)+additionalSize, GB_IGNORE_CASE);
603    GBDATA  *gb_pseudo;
604
605    for (gb_pseudo = GEN_first_pseudo_species(gb_main);
606         gb_pseudo;
607         gb_pseudo = GEN_next_pseudo_species(gb_pseudo))
608    {
609        GEN_add_pseudo_species_to_hash(gb_pseudo, pseudo_hash);
610    }
611
612    return pseudo_hash;
613}
614
615GBDATA *GEN_find_pseudo_species(GBDATA *gb_main, const char *organism_name, const char *gene_name, GB_HASH *pseudo_hash) {
616    // parameter pseudo_hash :
617    // 0 -> use slow direct search [if you only search one]
618    // otherwise it shall be a hash generated by GEN_create_pseudo_species_hash() [if you search several times]
619    // Note : use GEN_add_pseudo_species_to_hash to keep hash up-to-date
620    GBDATA *gb_pseudo;
621
622    if (pseudo_hash) {
623        gb_pseudo = GEN_read_pseudo_species_from_hash(pseudo_hash, organism_name, gene_name);
624    }
625    else {
626        for (gb_pseudo = GEN_first_pseudo_species(gb_main);
627             gb_pseudo;
628             gb_pseudo = GEN_next_pseudo_species(gb_pseudo))
629        {
630            const char *origin_gene_name = GEN_origin_gene(gb_pseudo);
631            if (strcmp(gene_name, origin_gene_name) == 0) {
632                const char *origin_species_name = GEN_origin_organism(gb_pseudo);
633                if (strcmp(organism_name, origin_species_name) == 0) {
634                    break; // found pseudo species
635                }
636            }
637        }
638    }
639    return gb_pseudo;
640}
641
642// -----------------------
643//      search origins
644
645GBDATA *GEN_find_origin_organism(GBDATA *gb_pseudo, GB_HASH *organism_hash) {
646    // parameter organism_hash:
647    // 0 -> use slow direct search [if you only search one or two]
648    // otherwise it shall be a hash generated by GBT_create_organism_hash() [if you search several times]
649    // Note : use GBT_add_item_to_hash() to keep hash up-to-date
650
651    const char *origin_species_name;
652    GBDATA     *gb_organism = 0;
653    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
654
655    origin_species_name = GEN_origin_organism(gb_pseudo);
656    if (origin_species_name) {
657        if (organism_hash) {
658            gb_organism = (GBDATA*)GBS_read_hash(organism_hash, origin_species_name);
659        }
660        else {
661            gb_organism = GBT_find_species_rel_species_data(GB_get_father(gb_pseudo), origin_species_name);
662        }
663    }
664
665    return gb_organism;
666}
667
668GBDATA *GEN_find_origin_gene(GBDATA *gb_pseudo, GB_HASH *organism_hash) {
669    const char *origin_gene_name;
670
671    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
672
673    origin_gene_name = GEN_origin_gene(gb_pseudo);
674    if (origin_gene_name) {
675        GBDATA *gb_organism = GEN_find_origin_organism(gb_pseudo, organism_hash);
676        gb_assert(gb_organism);
677
678        return GEN_find_gene(gb_organism, origin_gene_name);
679    }
680    return 0;
681}
682
683//  --------------------------------
684//      find pseudo-species
685
686GBDATA* GEN_first_pseudo_species(GBDATA *gb_main) {
687    GBDATA *gb_species = GBT_first_species(gb_main);
688
689    if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) return gb_species;
690    return GEN_next_pseudo_species(gb_species);
691}
692
693GBDATA* GEN_next_pseudo_species(GBDATA *gb_species) {
694    if (gb_species) {
695        while (1) {
696            gb_species = GBT_next_species(gb_species);
697            if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) break;
698        }
699    }
700    return gb_species;
701}
702
703static GBDATA* GEN_next_marked_pseudo_species(GBDATA *gb_species) {
704    if (gb_species) {
705        while (1) {
706            gb_species = GBT_next_marked_species(gb_species);
707            if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) break;
708        }
709    }
710    return gb_species;
711}
712
713GBDATA *GEN_first_marked_pseudo_species(GBDATA *gb_main) {
714    GBDATA *gb_species = GBT_first_marked_species(gb_main);
715
716    if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) return gb_species;
717    return GEN_next_marked_pseudo_species(gb_species);
718}
719
720// ------------------
721//      organisms
722
723bool GEN_is_organism(GBDATA *gb_species) {
724    gb_assert(GEN_is_genome_db(GB_get_root(gb_species), -1)); // assert this is a genome db
725    // otherwise it is an error to use GEN_is_organism (or its callers)!!!!
726
727    return GB_entry(gb_species, GENOM_ALIGNMENT) != 0;
728}
729
730GBDATA *GEN_find_organism(GBDATA *gb_main, const char *name) {
731    GBDATA *gb_orga = GBT_find_species(gb_main, name);
732    if (gb_orga) {
733        if (!GEN_is_organism(gb_orga)) {
734            fprintf(stderr, "ARBDB-warning: found unspecific species named '%s', but expected an 'organism' with that name\n", name);
735            gb_orga = 0;
736        }
737    }
738    return gb_orga;
739}
740
741GBDATA *GEN_first_organism(GBDATA *gb_main) {
742    GBDATA *gb_organism = GBT_first_species(gb_main);
743
744    if (!gb_organism || GEN_is_organism(gb_organism)) return gb_organism;
745    return GEN_next_organism(gb_organism);
746}
747GBDATA *GEN_next_organism(GBDATA *gb_organism) {
748    if (gb_organism) {
749        while (1) {
750            gb_organism = GBT_next_species(gb_organism);
751            if (!gb_organism || GEN_is_organism(gb_organism)) break;
752        }
753    }
754    return gb_organism;
755
756}
757
758long GEN_get_organism_count(GBDATA *gb_main) {
759    long    count       = 0;
760    GBDATA *gb_organism = GEN_first_organism(gb_main);
761    while (gb_organism) {
762        count++;
763        gb_organism = GEN_next_organism(gb_organism);
764    }
765    return count;
766}
767
768
769GBDATA *GEN_first_marked_organism(GBDATA *gb_main) {
770    GBDATA *gb_organism = GBT_first_marked_species(gb_main);
771
772    if (!gb_organism || GEN_is_organism(gb_organism)) return gb_organism;
773    return GEN_next_marked_organism(gb_organism);
774}
775GBDATA *GEN_next_marked_organism(GBDATA *gb_organism) {
776    if (gb_organism) {
777        while (1) {
778            gb_organism = GBT_next_marked_species(gb_organism);
779            if (!gb_organism || GEN_is_organism(gb_organism)) break;
780        }
781    }
782    return gb_organism;
783}
784
785char *GEN_global_gene_identifier(GBDATA *gb_gene, GBDATA *gb_organism) {
786    if (!gb_organism) {
787        gb_organism = GB_get_grandfather(gb_gene);
788        gb_assert(gb_organism);
789    }
790
791    return GBS_global_string_copy("%s/%s", GBT_read_name(gb_organism), GBT_read_name(gb_gene));
792}
793
794// --------------------------------------------------------------------------------
795
796#ifdef UNIT_TESTS
797#include <test_unit.h>
798#include <arb_unit_test.h>
799#include <arb_defs.h>
800
801static struct arb_unit_test::test_alignment_data TestAlignmentData_Genome[] = {
802    { 0, "spec", "AUCUCCUAAACCCAACCGUAGUUCGAAUUGAG" },
803};
804
805#define TEST_ASSERT_MEMBER_EQUAL(s1,s2,member) TEST_ASSERT_EQUAL((s1)->member, (s2)->member)
806
807#define TEST_ASSERT_GENPOS_EQUAL(p1,p2) do {                            \
808        TEST_ASSERT_MEMBER_EQUAL(p1, p2, parts);                        \
809        TEST_ASSERT_MEMBER_EQUAL(p1, p2, joinable);                     \
810        for (int p = 0; p<(p1)->parts; ++p) {                           \
811            TEST_ASSERT_MEMBER_EQUAL(p1, p2, start_pos[p]);             \
812            TEST_ASSERT_MEMBER_EQUAL(p1, p2, stop_pos[p]);              \
813            TEST_ASSERT_MEMBER_EQUAL(p1, p2, complement[p]);            \
814            if ((p1)->start_uncertain) {                                \
815                TEST_ASSERT_MEMBER_EQUAL(p1, p2, start_uncertain[p]);   \
816                TEST_ASSERT_MEMBER_EQUAL(p1, p2, stop_uncertain[p]);    \
817            }                                                           \
818        }                                                               \
819    } while(0)
820
821#define TEST_WRITE_READ_GEN_POSITION(pos)                               \
822    do {                                                                \
823        error = GEN_write_position(gb_gene, (pos));                     \
824        if (!error) {                                                   \
825            GEN_position *rpos = GEN_read_position(gb_gene);            \
826            if (!rpos) {                                                \
827                error = GB_await_error();                               \
828            }                                                           \
829            else {                                                      \
830                TEST_ASSERT_GENPOS_EQUAL((pos), rpos);                  \
831                GEN_free_position(rpos);                                \
832            }                                                           \
833        }                                                               \
834        TEST_ASSERT_NULL(error.deliver());                              \
835    } while(0)
836
837#define TEST_WRITE_GEN_POSITION_ERROR(pos,exp_error) do {               \
838        error = GEN_write_position(gb_gene, &*(pos));                   \
839        TEST_ASSERT_EQUAL(error.deliver(), exp_error);                  \
840    } while(0)
841   
842#define TEST_GENPOS_FIELD(field,value) do {                             \
843        GBDATA *gb_field = GB_entry(gb_gene, (field));                  \
844        if ((value)) {                                                  \
845            TEST_ASSERT(gb_field);                                      \
846            TEST_ASSERT_EQUAL(GB_read_char_pntr(gb_field), (value));    \
847        }                                                               \
848        else {                                                          \
849            TEST_ASSERT(!gb_field);                                     \
850        }                                                               \
851    } while(0)
852
853#define TEST_GENPOS_FIELDS(start,stop,complement,certain) do {          \
854        TEST_GENPOS_FIELD("pos_start", start);                          \
855        TEST_GENPOS_FIELD("pos_stop", stop);                            \
856        TEST_GENPOS_FIELD("pos_complement", complement);                \
857        TEST_GENPOS_FIELD("pos_certain", certain);                      \
858    } while(0)
859
860#define TEST_GENE_SEQ_AND_LENGTH(werr,wseq,wlen) do {                   \
861        size_t len;                                                     \
862        char *seq = GBT_read_gene_sequence_and_length(gb_gene, true, '-', &len); \
863        TEST_ASSERT_EQUAL(GB_have_error(), werr);                       \
864        if (seq) {                                                      \
865            TEST_ASSERT_EQUAL(len, (size_t)(wlen));                     \
866            TEST_ASSERT_EQUAL(seq, (wseq));                             \
867            free(seq);                                                  \
868        }                                                               \
869        else {                                                          \
870            GB_clear_error();                                           \
871        }                                                               \
872    } while(0)
873   
874void TEST_GEN_position() {
875    // see also ../GENOM_IMPORT/Location.cxx@TEST_gene_location
876
877    GB_shell   shell;
878    ARB_ERROR  error;
879    GBDATA    *gb_main = TEST_CREATE_DB(error, "ali_genom", TestAlignmentData_Genome, false);
880
881    TEST_ASSERT_NULL(error.deliver());
882
883    {
884        GB_transaction ta(gb_main);
885
886        GBDATA *gb_organism  = GBT_find_species(gb_main, "spec"); TEST_ASSERT(gb_organism);
887        GBDATA *gb_gene_data = GEN_findOrCreate_gene_data(gb_organism); TEST_ASSERT(gb_gene_data);
888        GBDATA *gb_gene      = GEN_create_nonexisting_gene_rel_gene_data(gb_gene_data, "gene"); TEST_ASSERT(gb_gene);
889
890        typedef SmartCustomPtr(GEN_position, GEN_free_position) GEN_position_Ptr;
891        GEN_position_Ptr pos;
892
893        pos = GEN_new_position(1, false);
894
895        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal start position 0");
896
897        pos->start_pos[0]  = 5;
898        pos->stop_pos[0]   = 10;
899        pos->complement[0] = 1;
900
901        GEN_use_uncertainties(&*pos);
902
903        TEST_WRITE_READ_GEN_POSITION(&*pos);
904        TEST_GENPOS_FIELDS("5", "10", "1", "==");
905
906        TEST_GENE_SEQ_AND_LENGTH(false, "TTTAGG", 6);
907
908        // ----------
909
910        pos = GEN_new_position(3, false);
911
912        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal start position 0");
913
914        GEN_use_uncertainties(&*pos);
915
916        pos->start_pos[0]  = 5;   pos->start_pos[1]  = 10;  pos->start_pos[2]  = 25;
917        pos->stop_pos[0]   = 15;  pos->stop_pos[1]   = 20;  pos->stop_pos[2]   = 25;
918        pos->complement[0] = 0;   pos->complement[1] = 1;   pos->complement[2] = 0;
919
920        pos->start_uncertain[0] = '<';
921        pos->stop_uncertain[2] = '>';
922
923        TEST_WRITE_READ_GEN_POSITION(&*pos);
924        TEST_GENPOS_FIELDS("5,10,25", "15,20,25", "0,1,0", "<=,==,=>");
925
926        TEST_GENE_SEQ_AND_LENGTH(false, "CCUAAACCCAA-TACGGTTGGGT-G", 25);
927
928        pos->stop_uncertain[2] = 'x';
929        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid uncertainty 'x'"); 
930
931        pos->stop_uncertain[2] = '+';
932        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid uncertainty '+'"); // invalid for stop
933       
934        pos->start_uncertain[2] = '+';
935        pos->stop_uncertain[2]  = '-';
936        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid positions 25^25 for uncertainties +-");
937
938        pos->stop_pos[2] = 26;
939        TEST_WRITE_GEN_POSITION_ERROR(pos, (void*)NULL);
940       
941        pos->stop_pos[0] = 100;
942        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal stop position 100 (>length(=32))");
943    }
944
945    GB_close(gb_main);
946}
947
948#endif // UNIT_TESTS
Note: See TracBrowser for help on using the browser.