source: tags/ms_r17q3/ARBDB/adGene.cxx

Last change on this file was 15176, checked in by westram, 8 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.6 KB
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        ARB_calloc(pos, 1);
166        pos->parts      = parts;
167        pos->start_pos  = (size_t*)ARB_calloc<char>(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) { // split
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, long seqLength) {
321    // if 'seqLength' is != 0, it is used to check the correctness of 'pos'
322    // (otherwise this function reads the genome-sequence to detect its length)
323
324    GB_ERROR  error          = 0;
325    GBDATA   *gb_pos_joined  = GB_entry(gb_gene, "pos_joined");
326    GBDATA   *gb_pos_certain = GB_entry(gb_gene, "pos_certain");
327    GBDATA   *gb_pos_start;
328    GBDATA   *gb_pos_stop;
329    GBDATA   *gb_pos_complement;
330    int       p;
331
332    gb_assert(pos);
333
334    gb_pos_start             = GB_search(gb_gene, "pos_start", GB_STRING);
335    if (!gb_pos_start) error = GB_await_error();
336
337    if (!error) {
338        gb_pos_stop             = GB_search(gb_gene, "pos_stop", GB_STRING);
339        if (!gb_pos_stop) error = GB_await_error();
340    }
341    if (!error) {
342        gb_pos_complement             = GB_search(gb_gene, "pos_complement", GB_STRING);
343        if (!gb_pos_complement) error = GB_await_error();
344    }
345
346    if (!error) {
347        if (pos->start_uncertain) {
348            if (!gb_pos_certain) {
349                gb_pos_certain             = GB_search(gb_gene, "pos_certain", GB_STRING);
350                if (!gb_pos_certain) error = GB_await_error();
351            }
352        }
353        else {
354            if (gb_pos_certain) {
355                error          = GB_delete(gb_pos_certain);
356                gb_pos_certain = 0;
357            }
358        }
359    }
360
361    // test data
362    if (!error) {
363        size_t length;
364
365        if (seqLength) {
366            length = seqLength;
367        }
368        else { // unknown -> autodetect
369            GBDATA *gb_organism = GB_get_grandfather(gb_gene);
370            GBDATA *gb_genome   = GBT_find_sequence(gb_organism, GENOM_ALIGNMENT);
371
372            length = GB_read_count(gb_genome);
373        }
374
375        for (p = 0; p<pos->parts && !error; ++p) {
376            char c;
377
378            c = pos->complement[p]; gb_assert(c == 0 || c == 1);
379            if (c<0 || c>1) {
380                error = GBS_global_string("Illegal value %i in complement", int(c));
381            }
382            else {
383                if (pos->start_pos[p]>pos->stop_pos[p]) {
384                    error = GBS_global_string("Illegal positions (%zu>%zu)", pos->start_pos[p], pos->stop_pos[p]);
385                }
386                else if (pos->start_pos[p] == 0) {
387                    error = GBS_global_string("Illegal start position %zu", pos->start_pos[p]);
388                }
389                else if (pos->stop_pos[p] > length) {
390                    error = GBS_global_string("Illegal stop position %zu (>length(=%zu))", pos->stop_pos[p], length);
391                }
392                else {
393                    if (pos->start_uncertain) {
394                        c       = pos->start_uncertain[p];
395                        char c2 = pos->stop_uncertain[p];
396
397                        if      (!c  || strchr("<=>+", c)  == 0) error = GBS_global_string("Invalid uncertainty '%c'", c);
398                        else if (!c2 || strchr("<=>-", c2) == 0) error = GBS_global_string("Invalid uncertainty '%c'", c2);
399                        else {
400                            if (c == '+' || c2 == '-') {
401                                if (c == '+' && c2 == '-') {
402                                    if (pos->start_pos[p] != pos->stop_pos[p]-1) {
403                                        error = GBS_global_string("Invalid positions %zu^%zu for uncertainties +-", pos->start_pos[p], pos->stop_pos[p]);
404                                    }
405                                }
406                                else {
407                                    error = "uncertainties '+' and '-' can only be used together";
408                                }
409                            }
410                        }
411                    }
412                }
413            }
414        }
415    }
416
417    if (!error) {
418        if (pos->parts == 1) {
419            if (gb_pos_joined) error = GB_delete(gb_pos_joined);
420
421            if (!error) error = GB_write_string(gb_pos_start,      GBS_global_string("%zu", pos->start_pos[0]));
422            if (!error) error = GB_write_string(gb_pos_stop,       GBS_global_string("%zu", pos->stop_pos[0]));
423            if (!error) error = GB_write_string(gb_pos_complement, GBS_global_string("%c", pos->complement[0]+'0'));
424
425            if (!error && gb_pos_certain) {
426                error = GB_write_string(gb_pos_certain, GBS_global_string("%c%c", pos->start_uncertain[0], pos->stop_uncertain[0]));
427            }
428        }
429        else {
430            if (!gb_pos_joined) {
431                gb_pos_joined             = GB_search(gb_gene, "pos_joined", GB_INT);
432                if (!gb_pos_joined) error = GB_await_error();
433            }
434            if (!error) error = GB_write_int(gb_pos_joined, pos->parts * (pos->joinable ? 1 : -1)); // neg. parts means not joinable
435
436            if (!error) {
437                GBS_strstruct *start      = GBS_stropen(12*pos->parts);
438                GBS_strstruct *stop       = GBS_stropen(12*pos->parts);
439                GBS_strstruct *complement = GBS_stropen(2*pos->parts);
440                GBS_strstruct *uncertain  = GBS_stropen(3*pos->parts);
441
442                for (p = 0; p<pos->parts; ++p) {
443                    if (p>0) {
444                        GBS_chrcat(start, ',');
445                        GBS_chrcat(stop, ',');
446                        GBS_chrcat(complement, ',');
447                        GBS_chrcat(uncertain, ',');
448                    }
449                    GBS_strcat(start, GBS_global_string("%zu", pos->start_pos[p]));
450                    GBS_strcat(stop,  GBS_global_string("%zu", pos->stop_pos[p]));
451                    GBS_chrcat(complement, pos->complement[p]+'0');
452                    if (gb_pos_certain) {
453                        GBS_chrcat(uncertain, pos->start_uncertain[p]);
454                        GBS_chrcat(uncertain, pos->stop_uncertain[p]);
455                    }
456                }
457
458                char *sstart      = GBS_strclose(start);
459                char *sstop       = GBS_strclose(stop);
460                char *scomplement = GBS_strclose(complement);
461                char *suncertain  = GBS_strclose(uncertain);
462
463                error             = GB_write_string(gb_pos_start, sstart);
464                if (!error) error = GB_write_string(gb_pos_stop, sstop);
465                if (!error) error = GB_write_string(gb_pos_complement, scomplement);
466                if (!error && gb_pos_certain) error = GB_write_string(gb_pos_certain, suncertain);
467
468                free(suncertain);
469                free(scomplement);
470                free(sstop);
471                free(sstart);
472            }
473        }
474    }
475
476    return error;
477}
478
479static GEN_position *location2sort = 0;
480
481static int cmp_location_parts(const void *v1, const void *v2) {
482    int i1 = *(int*)v1;
483    int i2 = *(int*)v2;
484
485    int cmp = location2sort->start_pos[i1]-location2sort->start_pos[i2];
486    if (!cmp) {
487        cmp = location2sort->stop_pos[i1]-location2sort->stop_pos[i2];
488    }
489    return cmp;
490}
491
492void GEN_sortAndMergeLocationParts(GEN_position *location) {
493    // Note: makes location partly invalid (only start_pos + stop_pos are valid afterwards)
494    int  parts = location->parts;
495    int *idx   = ARB_alloc<int>(parts); // idx[newpos] = oldpos
496    int  i, p;
497
498    for (p = 0; p<parts; ++p) idx[p] = p;
499
500    location2sort = location;
501    qsort(idx, parts, sizeof(*idx), cmp_location_parts);
502    location2sort = 0;
503
504    for (p = 0; p<parts; ++p) {
505        i = idx[p];
506
507#define swap(a, b, type) do { type tmp = (a); (a) = (b); (b) = (tmp); } while (0)
508
509        if (i != p) {
510            swap(location->start_pos[i],  location->start_pos[p],  size_t);
511            swap(location->stop_pos[i],   location->stop_pos[p],   size_t);
512            swap(idx[i], idx[p], int);
513        }
514    }
515
516#if defined(DEBUG) && 0
517    printf("Locations sorted:\n");
518    for (p = 0; p<parts; ++p) {
519        printf("  [%i] %i - %i %i\n", p, location->start_pos[p], location->stop_pos[p], (int)(location->complement[p]));
520    }
521#endif // DEBUG
522
523    i = 0;
524    for (p = 1; p<parts; p++) {
525        if ((location->stop_pos[i]+1) >= location->start_pos[p]) {
526            // parts overlap or are directly consecutive
527
528            location->stop_pos[i]  = location->stop_pos[p];
529        }
530        else {
531            i++;
532            location->start_pos[i] = location->start_pos[p];
533            location->stop_pos[i]  = location->stop_pos[p];
534        }
535    }
536    location->parts = i+1;
537
538#if defined(DEBUG) && 0
539    parts = location->parts;
540    printf("Locations merged:\n");
541    for (p = 0; p<parts; ++p) {
542        printf("  [%i] %i - %i %i\n", p, location->start_pos[p], location->stop_pos[p], (int)(location->complement[p]));
543    }
544#endif // DEBUG
545
546    free(idx);
547}
548
549
550
551//  -----------------------------------------
552//      test if species is pseudo-species
553
554const char *GEN_origin_organism(GBDATA *gb_pseudo) {
555    GBDATA *gb_origin = GB_entry(gb_pseudo, "ARB_origin_species");
556    return gb_origin ? GB_read_char_pntr(gb_origin) : 0;
557}
558const char *GEN_origin_gene(GBDATA *gb_pseudo) {
559    GBDATA *gb_origin = GB_entry(gb_pseudo, "ARB_origin_gene");
560    return gb_origin ? GB_read_char_pntr(gb_origin) : 0;
561}
562
563bool GEN_is_pseudo_gene_species(GBDATA *gb_species) {
564    return GEN_origin_organism(gb_species) != 0;
565}
566
567//  ------------------------------------------------
568//      find organism or gene for pseudo-species
569
570GB_ERROR GEN_organism_not_found(GBDATA *gb_pseudo) {
571    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
572    gb_assert(GEN_find_origin_organism(gb_pseudo, 0) == 0);
573
574    return GB_export_errorf("The gene-species '%s' refers to an unknown organism (%s)\n"
575                            "This occurs if you rename or delete the organism or change the entry\n"
576                            "'ARB_origin_species' and will most likely cause serious problems.",
577                            GBT_read_name(gb_pseudo),
578                            GEN_origin_organism(gb_pseudo));
579}
580
581// @@@ FIXME: missing: GEN_gene_not_found (like GEN_organism_not_found)
582
583// ------------------------------
584//      search pseudo species
585
586static const char *pseudo_species_hash_key(const char *organism_name, const char *gene_name) {
587    return GBS_global_string("%s*%s", organism_name, gene_name);
588}
589
590static GBDATA *GEN_read_pseudo_species_from_hash(const GB_HASH *pseudo_hash, const char *organism_name, const char *gene_name) {
591    return (GBDATA*)GBS_read_hash(pseudo_hash, pseudo_species_hash_key(organism_name, gene_name));
592}
593
594void GEN_add_pseudo_species_to_hash(GBDATA *gb_pseudo, GB_HASH *pseudo_hash) {
595    const char *organism_name = GEN_origin_organism(gb_pseudo);
596    const char *gene_name     = GEN_origin_gene(gb_pseudo);
597
598    gb_assert(organism_name);
599    gb_assert(gene_name);
600
601    GBS_write_hash(pseudo_hash, pseudo_species_hash_key(organism_name, gene_name), (long)gb_pseudo);
602}
603
604GB_HASH *GEN_create_pseudo_species_hash(GBDATA *gb_main, int additionalSize) {
605    GB_HASH *pseudo_hash = GBS_create_hash(GBT_get_species_count(gb_main)+additionalSize, GB_IGNORE_CASE);
606    GBDATA  *gb_pseudo;
607
608    for (gb_pseudo = GEN_first_pseudo_species(gb_main);
609         gb_pseudo;
610         gb_pseudo = GEN_next_pseudo_species(gb_pseudo))
611    {
612        GEN_add_pseudo_species_to_hash(gb_pseudo, pseudo_hash);
613    }
614
615    return pseudo_hash;
616}
617
618GBDATA *GEN_find_pseudo_species(GBDATA *gb_main, const char *organism_name, const char *gene_name, const GB_HASH *pseudo_hash) {
619    // parameter pseudo_hash :
620    // 0 -> use slow direct search [if you only search one]
621    // otherwise it shall be a hash generated by GEN_create_pseudo_species_hash() [if you search several times]
622    // Note : use GEN_add_pseudo_species_to_hash to keep hash up-to-date
623    GBDATA *gb_pseudo;
624
625    if (pseudo_hash) {
626        gb_pseudo = GEN_read_pseudo_species_from_hash(pseudo_hash, organism_name, gene_name);
627    }
628    else {
629        for (gb_pseudo = GEN_first_pseudo_species(gb_main);
630             gb_pseudo;
631             gb_pseudo = GEN_next_pseudo_species(gb_pseudo))
632        {
633            const char *origin_gene_name = GEN_origin_gene(gb_pseudo);
634            if (strcmp(gene_name, origin_gene_name) == 0) {
635                const char *origin_species_name = GEN_origin_organism(gb_pseudo);
636                if (strcmp(organism_name, origin_species_name) == 0) {
637                    break; // found pseudo species
638                }
639            }
640        }
641    }
642    return gb_pseudo;
643}
644
645// -----------------------
646//      search origins
647
648GBDATA *GEN_find_origin_organism(GBDATA *gb_pseudo, const GB_HASH *organism_hash) {
649    // parameter organism_hash:
650    // 0 -> use slow direct search [if you only search one or two]
651    // otherwise it shall be a hash generated by GBT_create_organism_hash() [if you search several times]
652    // Note : use GBT_add_item_to_hash() to keep hash up-to-date
653
654    const char *origin_species_name;
655    GBDATA     *gb_organism = 0;
656    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
657
658    origin_species_name = GEN_origin_organism(gb_pseudo);
659    if (origin_species_name) {
660        if (organism_hash) {
661            gb_organism = (GBDATA*)GBS_read_hash(organism_hash, origin_species_name);
662        }
663        else {
664            gb_organism = GBT_find_species_rel_species_data(GB_get_father(gb_pseudo), origin_species_name);
665        }
666    }
667
668    return gb_organism;
669}
670
671GBDATA *GEN_find_origin_gene(GBDATA *gb_pseudo, const GB_HASH *organism_hash) {
672    const char *origin_gene_name;
673
674    gb_assert(GEN_is_pseudo_gene_species(gb_pseudo));
675
676    origin_gene_name = GEN_origin_gene(gb_pseudo);
677    if (origin_gene_name) {
678        GBDATA *gb_organism = GEN_find_origin_organism(gb_pseudo, organism_hash);
679        gb_assert(gb_organism);
680
681        return GEN_find_gene(gb_organism, origin_gene_name);
682    }
683    return 0;
684}
685
686//  --------------------------------
687//      find pseudo-species
688
689GBDATA* GEN_first_pseudo_species(GBDATA *gb_main) {
690    GBDATA *gb_species = GBT_first_species(gb_main);
691
692    if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) return gb_species;
693    return GEN_next_pseudo_species(gb_species);
694}
695
696GBDATA* GEN_next_pseudo_species(GBDATA *gb_species) {
697    if (gb_species) {
698        while (1) {
699            gb_species = GBT_next_species(gb_species);
700            if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) break;
701        }
702    }
703    return gb_species;
704}
705
706static GBDATA* GEN_next_marked_pseudo_species(GBDATA *gb_species) {
707    if (gb_species) {
708        while (1) {
709            gb_species = GBT_next_marked_species(gb_species);
710            if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) break;
711        }
712    }
713    return gb_species;
714}
715
716GBDATA *GEN_first_marked_pseudo_species(GBDATA *gb_main) {
717    GBDATA *gb_species = GBT_first_marked_species(gb_main);
718
719    if (!gb_species || GEN_is_pseudo_gene_species(gb_species)) return gb_species;
720    return GEN_next_marked_pseudo_species(gb_species);
721}
722
723// ------------------
724//      organisms
725
726bool GEN_is_organism(GBDATA *gb_species) {
727    gb_assert(GEN_is_genome_db(GB_get_root(gb_species), -1)); // assert this is a genome db
728    // otherwise it is an error to use GEN_is_organism (or its callers)!!!!
729
730    return GB_entry(gb_species, GENOM_ALIGNMENT) != 0;
731}
732
733GBDATA *GEN_find_organism(GBDATA *gb_main, const char *name) {
734    GBDATA *gb_orga = GBT_find_species(gb_main, name);
735    if (gb_orga) {
736        if (!GEN_is_organism(gb_orga)) {
737            fprintf(stderr, "ARBDB-warning: found unspecific species named '%s', but expected an 'organism' with that name\n", name);
738            gb_orga = 0;
739        }
740    }
741    return gb_orga;
742}
743
744GBDATA *GEN_first_organism(GBDATA *gb_main) {
745    GBDATA *gb_organism = GBT_first_species(gb_main);
746
747    if (!gb_organism || GEN_is_organism(gb_organism)) return gb_organism;
748    return GEN_next_organism(gb_organism);
749}
750GBDATA *GEN_next_organism(GBDATA *gb_organism) {
751    if (gb_organism) {
752        while (1) {
753            gb_organism = GBT_next_species(gb_organism);
754            if (!gb_organism || GEN_is_organism(gb_organism)) break;
755        }
756    }
757    return gb_organism;
758
759}
760
761long GEN_get_organism_count(GBDATA *gb_main) {
762    long    count       = 0;
763    GBDATA *gb_organism = GEN_first_organism(gb_main);
764    while (gb_organism) {
765        count++;
766        gb_organism = GEN_next_organism(gb_organism);
767    }
768    return count;
769}
770
771
772GBDATA *GEN_first_marked_organism(GBDATA *gb_main) {
773    GBDATA *gb_organism = GBT_first_marked_species(gb_main);
774
775    if (!gb_organism || GEN_is_organism(gb_organism)) return gb_organism;
776    return GEN_next_marked_organism(gb_organism);
777}
778GBDATA *GEN_next_marked_organism(GBDATA *gb_organism) {
779    if (gb_organism) {
780        while (1) {
781            gb_organism = GBT_next_marked_species(gb_organism);
782            if (!gb_organism || GEN_is_organism(gb_organism)) break;
783        }
784    }
785    return gb_organism;
786}
787
788char *GEN_global_gene_identifier(GBDATA *gb_gene, GBDATA *gb_organism) {
789    if (!gb_organism) {
790        gb_organism = GB_get_grandfather(gb_gene);
791        gb_assert(gb_organism);
792    }
793
794    return GBS_global_string_copy("%s/%s", GBT_read_name(gb_organism), GBT_read_name(gb_gene));
795}
796
797// --------------------------------------------------------------------------------
798
799#ifdef UNIT_TESTS
800#include <test_unit.h>
801#include <arb_unit_test.h>
802#include <arb_defs.h>
803
804static struct arb_unit_test::test_alignment_data TestAlignmentData_Genome[] = {
805    { 0, "spec", "AUCUCCUAAACCCAACCGUAGUUCGAAUUGAG" },
806};
807
808#define TEST_EXPECT_MEMBER_EQUAL(s1,s2,member) TEST_EXPECT_EQUAL((s1)->member, (s2)->member)
809
810#define TEST_EXPECT_GENPOS_EQUAL(p1,p2) do {                            \
811        TEST_EXPECT_MEMBER_EQUAL(p1, p2, parts);                        \
812        TEST_EXPECT_MEMBER_EQUAL(p1, p2, joinable);                     \
813        for (int p = 0; p<(p1)->parts; ++p) {                           \
814            TEST_EXPECT_MEMBER_EQUAL(p1, p2, start_pos[p]);             \
815            TEST_EXPECT_MEMBER_EQUAL(p1, p2, stop_pos[p]);              \
816            TEST_EXPECT_MEMBER_EQUAL(p1, p2, complement[p]);            \
817            if ((p1)->start_uncertain) {                                \
818                TEST_EXPECT_MEMBER_EQUAL(p1, p2, start_uncertain[p]);   \
819                TEST_EXPECT_MEMBER_EQUAL(p1, p2, stop_uncertain[p]);    \
820            }                                                           \
821        }                                                               \
822    } while(0)
823
824#define TEST_WRITE_READ_GEN_POSITION(pos)                               \
825    do {                                                                \
826        error = GEN_write_position(gb_gene, (pos), 0);                  \
827        if (!error) {                                                   \
828            GEN_position *rpos = GEN_read_position(gb_gene);            \
829            if (!rpos) {                                                \
830                error = GB_await_error();                               \
831            }                                                           \
832            else {                                                      \
833                TEST_EXPECT_GENPOS_EQUAL((pos), rpos);                  \
834                GEN_free_position(rpos);                                \
835            }                                                           \
836        }                                                               \
837        TEST_EXPECT_NULL(error.deliver());                              \
838    } while(0)
839
840#define TEST_WRITE_GEN_POSITION_ERROR(pos,exp_error) do {               \
841        error = GEN_write_position(gb_gene, &*(pos), 0);                \
842        TEST_EXPECT_EQUAL(error.deliver(), exp_error);                  \
843    } while(0)
844   
845#define TEST_GENPOS_FIELD(field,value) do {                             \
846        GBDATA *gb_field = GB_entry(gb_gene, (field));                  \
847        if ((value)) {                                                  \
848            TEST_REJECT_NULL(gb_field);                              \
849            TEST_EXPECT_EQUAL(GB_read_char_pntr(gb_field), (value));    \
850        }                                                               \
851        else {                                                          \
852            TEST_EXPECT_NULL(gb_field);                                 \
853        }                                                               \
854    } while(0)
855
856#define TEST_GENPOS_FIELDS(start,stop,complement,certain) do {          \
857        TEST_GENPOS_FIELD("pos_start", start);                          \
858        TEST_GENPOS_FIELD("pos_stop", stop);                            \
859        TEST_GENPOS_FIELD("pos_complement", complement);                \
860        TEST_GENPOS_FIELD("pos_certain", certain);                      \
861    } while(0)
862
863#define TEST_GENE_SEQ_AND_LENGTH(werr,wseq,wlen) do {                   \
864        size_t len;                                                     \
865        char *seq = GBT_read_gene_sequence_and_length(gb_gene, true, '-', &len); \
866        TEST_EXPECT_EQUAL(GB_have_error(), werr);                       \
867        if (seq) {                                                      \
868            TEST_EXPECT_EQUAL(len, (size_t)(wlen));                     \
869            TEST_EXPECT_EQUAL(seq, (wseq));                             \
870            free(seq);                                                  \
871        }                                                               \
872        else {                                                          \
873            GB_clear_error();                                           \
874        }                                                               \
875    } while(0)
876   
877void TEST_GEN_position() {
878    // see also ../GENOM_IMPORT/Location.cxx@TEST_gene_location
879
880    GB_shell   shell;
881    ARB_ERROR  error;
882    GBDATA    *gb_main = TEST_CREATE_DB(error, "ali_genom", TestAlignmentData_Genome, false);
883
884    TEST_EXPECT_NULL(error.deliver());
885
886    {
887        GB_transaction ta(gb_main);
888
889        GBDATA *gb_organism  = GBT_find_species(gb_main, "spec"); TEST_REJECT_NULL(gb_organism);
890        GBDATA *gb_gene_data = GEN_findOrCreate_gene_data(gb_organism); TEST_REJECT_NULL(gb_gene_data);
891        GBDATA *gb_gene      = GEN_create_nonexisting_gene_rel_gene_data(gb_gene_data, "gene"); TEST_REJECT_NULL(gb_gene);
892
893        typedef SmartCustomPtr(GEN_position, GEN_free_position) GEN_position_Ptr;
894        GEN_position_Ptr pos;
895
896        pos = GEN_new_position(1, false);
897
898        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal start position 0");
899
900        pos->start_pos[0]  = 5;
901        pos->stop_pos[0]   = 10;
902        pos->complement[0] = 1;
903
904        GEN_use_uncertainties(&*pos);
905
906        TEST_WRITE_READ_GEN_POSITION(&*pos);
907        TEST_GENPOS_FIELDS("5", "10", "1", "==");
908
909        TEST_GENE_SEQ_AND_LENGTH(false, "TTTAGG", 6);
910
911        // ----------
912
913        pos = GEN_new_position(3, false);
914
915        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal start position 0");
916
917        GEN_use_uncertainties(&*pos);
918
919        pos->start_pos[0]  = 5;   pos->start_pos[1]  = 10;  pos->start_pos[2]  = 25;
920        pos->stop_pos[0]   = 15;  pos->stop_pos[1]   = 20;  pos->stop_pos[2]   = 25;
921        pos->complement[0] = 0;   pos->complement[1] = 1;   pos->complement[2] = 0;
922
923        pos->start_uncertain[0] = '<';
924        pos->stop_uncertain[2] = '>';
925
926        TEST_WRITE_READ_GEN_POSITION(&*pos);
927        TEST_GENPOS_FIELDS("5,10,25", "15,20,25", "0,1,0", "<=,==,=>");
928
929        TEST_GENE_SEQ_AND_LENGTH(false, "CCUAAACCCAA-TACGGTTGGGT-G", 25);
930
931        pos->stop_uncertain[2] = 'x';
932        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid uncertainty 'x'"); 
933
934        pos->stop_uncertain[2] = '+';
935        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid uncertainty '+'"); // invalid for stop
936       
937        pos->start_uncertain[2] = '+';
938        pos->stop_uncertain[2]  = '-';
939        TEST_WRITE_GEN_POSITION_ERROR(pos, "Invalid positions 25^25 for uncertainties +-");
940
941        pos->stop_pos[2] = 26;
942        TEST_WRITE_GEN_POSITION_ERROR(pos, (void*)NULL);
943       
944        pos->stop_pos[0] = 100;
945        TEST_WRITE_GEN_POSITION_ERROR(pos, "Illegal stop position 100 (>length(=32))");
946    }
947
948    GB_close(gb_main);
949}
950
951#endif // UNIT_TESTS
952
Note: See TracBrowser for help on using the repository browser.