source: tags/arb-6.0/TOOLS/arb_gene_probe.cxx

Last change on this file was 11888, checked in by westram, 8 years ago
  • highlight differences in production code triggered by build-flag 'UNIT_TESTS'
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.2 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : arb_gene_probe.cxx                                //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include <arbdbt.h>
12#include <adGene.h>
13
14#include <map>
15#include <list>
16#include <set>
17#include <string>
18
19#include <unistd.h>
20#include <sys/types.h>
21
22#define gp_assert(cond) arb_assert(cond)
23
24using namespace std;
25
26#if defined(DEBUG)
27// #define CREATE_DEBUG_FILES
28// #define DUMP_OVERLAP_CALC
29#endif // DEBUG
30
31// --------------------------------------------------------------------------------
32
33static int gene_counter          = 0; // pre-incremented counters
34static int splitted_gene_counter = 0;
35static int intergene_counter     = 0;
36
37struct nameOrder {
38    bool operator()(const char *name1, const char *name2) const {
39        // Normally it is sufficient to have any order, as long as it is strict.
40        // But for UNIT_TESTS we need a reproducable order, which does not
41        // depend on memory layout of DB elements.
42#if defined(UNIT_TESTS) // UT_DIFF
43        return strcmp(name1, name2)<0; // slow, determined by species names
44#else
45        return (name1-name2)<0;        // fast, but depends on memory layout (e.g. on MEMORY_TEST in gb_memory.h)
46#endif
47    }
48};
49
50typedef map<const char *, string, nameOrder> FullNameMap;
51static FullNameMap names;
52
53// --------------------------------------------------------------------------------
54
55struct PositionPair {
56    int begin;                  // these positions are in range [0 .. genome_length-1]
57    int end;
58
59    static int genome_length;
60
61#if defined(DEBUG)
62    void check_legal() const {
63        gp_assert(begin >= 0);
64        gp_assert(begin <= end);
65        gp_assert(end < genome_length);
66    }
67#endif // DEBUG
68
69    PositionPair() : begin(-1), end(-1) {}
70    PositionPair(int begin_, int end_) : begin(begin_), end(end_) {
71#if defined(DEBUG)
72        check_legal();
73#endif // DEBUG
74    }
75
76    int length() const { return end-begin+1; }
77
78    bool overlapsWith(const PositionPair& other) const {
79#if defined(DEBUG)
80        check_legal();
81        other.check_legal();
82#endif // DEBUG
83        return ! ((end < other.begin) || (other.end < begin));
84    }
85
86#if defined(DUMP_OVERLAP_CALC)
87    void dump(const char *note) const {
88        printf("%s begin=%i end=%i\n", note, begin, end);
89    }
90#endif // DUMP_OVERLAP_CALC
91};
92
93int PositionPair::genome_length = 0;
94
95typedef list<PositionPair> PositionPairList;
96
97struct ltNonOverlap {
98    // sorting with this operator identifies all overlapping PositionPair's as "equal"
99    bool operator ()(const PositionPair& p1, const PositionPair& p2) const {
100        return p1.end < p2.begin;
101    }
102};
103
104class GenePositionMap {
105    typedef set<PositionPair, ltNonOverlap> OverlappingGeneSet;
106
107    OverlappingGeneSet usedRanges;
108    unsigned long      overlapSize;
109    unsigned long      geneSize;
110public:
111    GenePositionMap() : overlapSize(0), geneSize(0) {}
112
113    void announceGene(PositionPair gene);
114    GB_ERROR buildIntergeneList(const PositionPair& wholeGenome, PositionPairList& intergeneList) const;
115    unsigned long getOverlap() const { return overlapSize; }
116    unsigned long getAllGeneSize() const { return geneSize; }
117
118#if defined(DUMP_OVERLAP_CALC)
119    void dump() const;
120#endif // DUMP_OVERLAP_CALC
121};
122
123// ____________________________________________________________
124// start of implementation of class GenePositionMap:
125
126void GenePositionMap::announceGene(PositionPair gene)
127{
128    OverlappingGeneSet::iterator found = usedRanges.find(gene);
129    if (found == usedRanges.end()) { // gene does not overlap with currently known ranges
130        usedRanges.insert(gene); // add to known ranges
131    }
132    else {
133        // 'found' overlaps with 'gene'
134        int gene_length = gene.length();
135
136        do {
137            gp_assert(gene.overlapsWith(*found));
138
139            gene                = PositionPair(min(found->begin, gene.begin), max(found->end, gene.end)); // calc combined range
140            int combined_length = gene.length();
141
142            size_t overlap  = (found->length()+gene_length)-combined_length;
143            overlapSize    += overlap;
144            geneSize       += gene_length;
145
146            usedRanges.erase(found);
147
148            gene_length = combined_length;
149            found       = usedRanges.find(gene); // search for further overlaps
150        } while (found != usedRanges.end());
151
152        usedRanges.insert(gene); // insert the combined range
153    }
154}
155
156GB_ERROR GenePositionMap::buildIntergeneList(const PositionPair& wholeGenome, PositionPairList& intergeneList) const
157{
158    OverlappingGeneSet::iterator end  = usedRanges.end();
159    OverlappingGeneSet::iterator curr = usedRanges.begin();
160    OverlappingGeneSet::iterator prev = end;
161
162    if (curr == end) { // nothing defined -> use whole genome as one big intergene
163        intergeneList.push_back(wholeGenome);
164    }
165    else {
166        if (curr->begin > wholeGenome.begin) { // intergene before first gene range ?
167            intergeneList.push_back(PositionPair(wholeGenome.begin, curr->begin-1));
168        }
169
170        prev = curr; ++curr;
171
172        while (curr != end) {
173            if (prev->end < curr->begin) {
174                if (prev->end != (curr->begin-1)) { // not directly adjacent
175                    intergeneList.push_back(PositionPair(prev->end+1, curr->begin-1));
176                }
177            }
178            else {
179                return "Internal error: Overlapping gene ranges";
180            }
181
182            prev = curr; ++curr;
183        }
184
185        if (prev != end && prev->end < wholeGenome.end) {
186            intergeneList.push_back(PositionPair(prev->end+1, wholeGenome.end));
187        }
188    }
189    return 0;
190}
191
192#if defined(DUMP_OVERLAP_CALC)
193void GenePositionMap::dump() const
194{
195    printf("List of ranges used by genes:\n");
196    for (OverlappingGeneSet::iterator g = usedRanges.begin(); g != usedRanges.end(); ++g) {
197        g->dump("- ");
198    }
199    printf("Overlap: %lu bases\n", getOverlap());
200}
201#endif // DUMP_OVERLAP_CALC
202
203// -end- of implementation of class GenePositionMap.
204
205static GB_ERROR create_data_entry(GBDATA *gb_species2, const char *sequence, int seqlen) {
206    GB_ERROR  error         = 0;
207    char     *gene_sequence = new char[seqlen+1];
208
209    memcpy(gene_sequence, sequence, seqlen);        // @@@ FIXME: avoid this copy!
210    gene_sequence[seqlen] = 0;
211
212    GBDATA *gb_ali     = GB_create_container(gb_species2, "ali_ptgene");
213    if (!gb_ali) error = GB_await_error();
214    else    error      = GBT_write_string(gb_ali, "data", gene_sequence);
215
216    delete [] gene_sequence;
217    return error;
218}
219
220#if defined(DEBUG)
221static void CHECK_SEMI_ESCAPED(const char *name) {
222    // checks whether all ";\\" are escaped
223    while (*name) {
224        gp_assert(*name != ';'); // oops, unescaped ';'
225        if (*name == '\\') ++name;
226        ++name;
227    }
228}
229#else
230#define CHECK_SEMI_ESCAPED(s)
231#endif // DEBUG
232
233
234static GBDATA *create_gene_species(GBDATA *gb_species_data2, const char *internal_name, const char *long_name, int abspos, const char *sequence, int length) {
235    // Note: 'sequence' is not necessarily 0-terminated!
236
237#if defined(DEBUG)
238    const char *firstSem = strchr(long_name, ';');
239    gp_assert(firstSem);
240    CHECK_SEMI_ESCAPED(firstSem+1);
241#endif // DEBUG
242
243    GB_ERROR  error       = GB_push_transaction(gb_species_data2);
244    GBDATA   *gb_species2 = 0;
245
246    if (!error) {
247        gb_species2 = GB_create_container(gb_species_data2, "species");
248        if (!gb_species2) error = GB_await_error();
249    }
250
251    if (!error) {
252        GBDATA *gb_name = GB_create(gb_species2, "name", GB_STRING);
253
254        if (!gb_name) error = GB_await_error();
255        else {
256            error = GB_write_string(gb_name, internal_name);
257            if (!error) {
258                const char *static_internal_name = GB_read_char_pntr(gb_name); // use static copy from db as map-index (internal_name is temporary)
259                error                            = create_data_entry(gb_species2, sequence, length);
260                if (!error) {
261                    names[static_internal_name] = long_name;
262                    error = GBT_write_int(gb_species2, "abspos", abspos);
263                }
264            }
265        }
266    }
267
268    error = GB_end_transaction(gb_species_data2, error);
269
270    if (error) { // be more verbose :
271        error       = GBS_global_string("%s (internal_name='%s', long_name='%s')", error, internal_name, long_name);
272        GB_export_error(error);
273        gb_species2 = NULL;
274    }
275
276    return gb_species2;
277}
278
279static GB_ERROR create_genelike_entry(const char *internal_name, GBDATA *gb_species_data2, int start_pos, int end_pos, const char *ali_genome, const char *long_name) {
280    GBDATA *gb_genespecies = create_gene_species(gb_species_data2, internal_name, long_name, start_pos, ali_genome+start_pos, end_pos-start_pos+1);
281    return gb_genespecies ? 0 : GB_await_error();
282}
283
284static GB_ERROR create_intergene(GBDATA *gb_species_data2, int start_pos, int end_pos, const char *ali_genome, const char *long_gene_name) {
285    if (start_pos <= end_pos) {
286        char internal_name[128];
287        sprintf(internal_name, "i%x", intergene_counter++);
288        return create_genelike_entry(internal_name, gb_species_data2, start_pos, end_pos, ali_genome, long_gene_name);
289    }
290    return "Illegal inter-gene positions (start behind end)";
291}
292
293static GB_ERROR create_gene(GBDATA *gb_species_data2, int start_pos, int end_pos, const char *ali_genome, const char *long_gene_name) {
294    if (start_pos <= end_pos) {
295        char internal_name[128];
296        sprintf(internal_name, "n%x", gene_counter++);
297        return create_genelike_entry(internal_name, gb_species_data2, start_pos, end_pos, ali_genome, long_gene_name);
298    }
299    return "Illegal gene positions (start behind end)";
300}
301
302static GB_ERROR create_splitted_gene(GBDATA *gb_species_data2, PositionPairList& part_list, const char *ali_genome, const char *long_gene_name) {
303    GB_ERROR                     error    = 0;
304    PositionPairList::iterator list_end = part_list.end();
305
306    int gene_size = 0;
307    for (PositionPairList::iterator part = part_list.begin(); part != list_end; ++part) {
308        int part_size  = part->end-part->begin+1;
309        gp_assert(part_size > 0);
310        gene_size     += part_size;
311    }
312    gp_assert(gene_size > 0);
313    char *gene_sequence = new char[gene_size+1];
314    int   gene_off      = 0;
315
316    char *split_pos_list = 0;   // contains split information: 'gene pos of part2,abs pos of part2;gene pos of part3,abs pos of part3;...'
317
318    for (PositionPairList::iterator part = part_list.begin(); part != list_end;) {
319        int part_size   = part->end-part->begin+1;
320        int genome_pos  = part->begin;
321        memcpy(gene_sequence+gene_off, ali_genome+part->begin, part_size);
322        gene_off       += part_size;
323
324        ++part;
325
326        if (split_pos_list == 0) { // first part
327            split_pos_list = GBS_global_string_copy("%i", gene_off); // gene offset of part 2
328        }
329        else {                  // next parts
330            char *next_split_pos_list;
331            if (part != list_end) { // not last
332                next_split_pos_list = GBS_global_string_copy("%s,%i;%i", split_pos_list, genome_pos, gene_off);
333            }
334            else { // last part
335                next_split_pos_list = GBS_global_string_copy("%s,%i", split_pos_list, genome_pos);
336            }
337            freeset(split_pos_list, next_split_pos_list);
338        }
339    }
340
341    char internal_name[128];
342    sprintf(internal_name, "s%x", splitted_gene_counter++);
343
344    const PositionPair&  first_part  = part_list.front();
345    GBDATA              *gb_species2 = create_gene_species(gb_species_data2, internal_name, long_gene_name, first_part.begin,
346                                                           gene_sequence, first_part.end-first_part.begin+1);
347
348    if (!gb_species2) error = GB_await_error();
349    else {
350#if defined(DEBUG) && 0
351        printf("splitted gene: long_gene_name='%s' internal_name='%s' split_pos_list='%s'\n",
352               long_gene_name, internal_name, split_pos_list);
353#endif // DEBUG
354        error = GBT_write_string(gb_species2, "splitpos", split_pos_list);
355    }
356
357    free(split_pos_list);
358    delete [] gene_sequence;
359
360    return error;
361}
362
363static GB_ERROR scan_gene_positions(GBDATA *gb_gene, PositionPairList& part_list) {
364    GB_ERROR      error    = 0;
365    GEN_position *location = GEN_read_position(gb_gene);
366
367    if (!location) error = GB_await_error();
368    else {
369        GEN_sortAndMergeLocationParts(location);
370        int parts = location->parts;
371        for (int p = 0; p<parts; ++p) {
372            part_list.push_back(PositionPair(location->start_pos[p]-1, location->stop_pos[p]-1));
373        }
374        GEN_free_position(location);
375    }
376    return error;
377}
378
379static GB_ERROR insert_genes_of_organism(GBDATA *gb_organism, GBDATA *gb_species_data2) {
380    // insert all genes of 'gb_organism' as pseudo-species
381    // into new 'species_data' (gb_species_data2)
382
383    GB_ERROR    error            = 0;
384    const char *organism_name = GBT_read_name(gb_organism);
385
386    GenePositionMap geneRanges;
387
388    int gene_counter_old          = gene_counter; // used for statistics only (see end of function)
389    int splitted_gene_counter_old = splitted_gene_counter;
390    int intergene_counter_old     = intergene_counter;
391
392    GBDATA *gb_ali_genom = GBT_read_sequence(gb_organism, GENOM_ALIGNMENT);
393    gp_assert(gb_ali_genom);                                                       // existence has to be checked by caller!
394
395    const char *ali_genom       = GB_read_char_pntr(gb_ali_genom);
396    if (!ali_genom) error       = GB_await_error();
397    PositionPair::genome_length = GB_read_count(gb_ali_genom);     // this affects checks in PositionPair
398
399    for (GBDATA *gb_gene = GEN_first_gene(gb_organism);
400         gb_gene && !error;
401         gb_gene = GEN_next_gene(gb_gene))
402    {
403        const char *gene_name = GBT_read_name(gb_gene);
404
405        PositionPairList part_list;
406        error = scan_gene_positions(gb_gene, part_list);
407
408        if (!error && part_list.empty()) error = "empty position list";
409        if (!error) {
410            int          split_count = part_list.size();
411            PositionPair first_part  = *part_list.begin();
412
413            if (!error) {
414                char *esc_gene_name = GBS_escape_string(gene_name, ";", '\\');
415                char *long_gene_name = GBS_global_string_copy("%s;%s", organism_name, esc_gene_name);
416                if (split_count == 1) { // normal gene
417                    error = create_gene(gb_species_data2, first_part.begin, first_part.end, ali_genom, long_gene_name);
418                    geneRanges.announceGene(first_part);
419                }
420                else {          // splitted gene
421                    error = create_splitted_gene(gb_species_data2, part_list, ali_genom, long_gene_name);
422
423                    for (PositionPairList::iterator p = part_list.begin(); p != part_list.end(); ++p) {
424                        geneRanges.announceGene(*p);
425                    }
426                }
427                free(long_gene_name);
428                free(esc_gene_name);
429            }
430        }
431
432        if (error && gene_name) error = GBS_global_string("in gene '%s': %s", gene_name, error);
433    }
434
435    if (!error) { // add intergenes
436        PositionPairList intergenes;
437        PositionPair     wholeGenome(0, PositionPair::genome_length-1);
438        error = geneRanges.buildIntergeneList(wholeGenome, intergenes);
439
440        for (PositionPairList::iterator i = intergenes.begin(); !error && i != intergenes.end(); ++i) {
441            char *long_intergene_name = GBS_global_string_copy("%s;intergene_%i_%i", organism_name, i->begin, i->end);
442            error                     = create_intergene(gb_species_data2, i->begin, i->end, ali_genom, long_intergene_name);
443            free(long_intergene_name);
444        }
445    }
446
447    if (error && organism_name) error = GBS_global_string("in organism '%s': %s", organism_name, error);
448
449    {
450        int new_genes          = gene_counter-gene_counter_old; // only non-splitted genes
451        int new_splitted_genes = splitted_gene_counter-splitted_gene_counter_old;
452        int new_intergenes     = intergene_counter-intergene_counter_old;
453
454        unsigned long genesSize    = geneRanges.getAllGeneSize();
455        unsigned long overlaps     = geneRanges.getOverlap();
456        double        data_grow    = overlaps/double(PositionPair::genome_length)*100;
457        double        gene_overlap = overlaps/double(genesSize)*100;
458
459        if (new_splitted_genes) {
460
461            printf("  - %s: %i genes (%i splitted), %i intergenes",
462                   organism_name, new_genes+new_splitted_genes, new_splitted_genes, new_intergenes);
463        }
464        else {
465            printf("  - %s: %i genes, %i intergenes",
466                   organism_name, new_genes, new_intergenes);
467        }
468        printf(" (data grow: %5.2f%%, gene overlap: %5.2f%%=%lu bp)\n", data_grow, gene_overlap, overlaps);
469    }
470
471#if defined(DUMP_OVERLAP_CALC)
472    geneRanges.dump();
473#endif // DUMP_OVERLAP_CALC
474
475    return error;
476}
477
478int ARB_main(int argc, char *argv[]) {
479
480    printf("\n"
481           "arb_gene_probe 1.2 -- (C) 2003/2004 Lehrstuhl fuer Mikrobiologie - TU Muenchen\n"
482           "written by Tom Littschwager, Bernd Spanfelner, Conny Wolf, Ralf Westram.\n");
483
484    if (argc != 3) {
485        printf("Usage: arb_gene_probe input_database output_database\n");
486        printf("       Prepares a genome database for Gene-PT-Server\n");
487        return EXIT_FAILURE;
488    }
489
490    const char *inputname  = argv[1];
491    const char *outputname = argv[2];
492
493    // GBK_terminate("test-crash of arb_gene_probe");
494   
495    printf("Converting '%s' -> '%s' ..\n", inputname, outputname);
496
497    GB_ERROR  error   = 0;
498    GB_shell  shell;
499    GBDATA   *gb_main = GB_open(inputname, "rw"); // rootzeiger wird gesetzt
500    if (!gb_main) {
501        error = GBS_global_string("Database '%s' not found", inputname);
502    }
503    else {
504        GB_request_undo_type(gb_main, GB_UNDO_NONE); // disable arbdb builtin undo
505        GB_begin_transaction(gb_main);
506
507        GBDATA *gb_species_data     = GBT_get_species_data(gb_main);
508        GBDATA *gb_species_data_new = GBT_create(gb_main, "species_data", 7); // create a second 'species_data' container
509
510        if (!gb_species_data_new) error = GB_await_error();
511
512        int non_ali_genom_species = 0;
513        int ali_genom_species     = 0;
514
515        for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data);
516             gb_species && !error;
517             gb_species = GBT_next_species(gb_species))
518        {
519            GBDATA *gb_ali_genom = GBT_read_sequence(gb_species, GENOM_ALIGNMENT);
520            if (!gb_ali_genom) {
521                // skip species w/o alignment 'GENOM_ALIGNMENT' (genome DBs often contain pseudo species)
522                ++non_ali_genom_species;
523            }
524            else {
525                error = insert_genes_of_organism(gb_species, gb_species_data_new);
526                ++ali_genom_species;
527            }
528        }
529
530        if (non_ali_genom_species) {
531            printf("%i species had no alignment in '" GENOM_ALIGNMENT "' and have been skipped.\n", non_ali_genom_species);
532        }
533        if (!error && ali_genom_species == 0) {
534            error = "no species with data in alignment '" GENOM_ALIGNMENT "' were found";
535        }
536
537        if (!error) {
538            printf("%i species had data in alignment '" GENOM_ALIGNMENT "'.\n"
539                   "Found %i genes (%i were splitted) and %i intergene regions.\n",
540                   ali_genom_species, gene_counter, splitted_gene_counter, intergene_counter);
541        }
542
543        if (!error) {
544            error = GB_delete(gb_species_data); // delete first (old) 'species_data' container
545        }
546
547        if (!error) {
548            // create map-string
549            char* map_string;
550            {
551                FullNameMap::iterator NameEnd = names.end();
552                FullNameMap::iterator NameIter;
553
554                size_t mapsize = 0;
555                for (NameIter = names.begin(); NameIter != NameEnd; ++NameIter) {
556                    mapsize += strlen(NameIter->first)+NameIter->second.length()+2;
557                }
558
559                map_string  = new char[mapsize+1];
560                size_t moff = 0;
561
562                for (NameIter = names.begin(); NameIter != NameEnd; ++NameIter) {
563                    int len1 = strlen(NameIter->first);
564                    int len2 = NameIter->second.length();
565
566                    memcpy(map_string+moff, NameIter->first, len1);
567                    map_string[moff+len1]  = ';';
568                    moff                  += len1+1;
569
570                    memcpy(map_string+moff, NameIter->second.c_str(), len2);
571                    map_string[moff+len2]  = ';';
572                    moff                  += len2+1;
573                }
574                map_string[moff] = 0;
575
576                gp_assert(moff <= mapsize);
577            }
578
579            GBDATA *gb_gene_map     = GB_create_container(gb_main, "gene_map");
580            if (!gb_gene_map) error = GB_await_error();
581            else    error           = GBT_write_string(gb_gene_map, "map_string", map_string);
582
583            delete [] map_string;
584        }
585
586        if (!error) {
587            // set default alignment for pt_server
588            error = GBT_set_default_alignment(gb_main, "ali_ptgene");
589
590            if (!error) {
591                GBDATA *gb_use     = GB_search(gb_main, "presets/alignment/alignment_name", GB_STRING);
592                if (!gb_use) error = GB_await_error();
593                else {
594                    GB_push_my_security(gb_main);
595                    error = GB_write_string(gb_use, "ali_ptgene");
596                    GB_pop_my_security(gb_main);
597                }
598            }
599        }
600
601        error = GB_end_transaction(gb_main, error);
602
603        if (!error) {
604            printf("Saving '%s' ..\n", outputname);
605            error = GB_save_as(gb_main, outputname, "bfm");
606            if (error) unlink(outputname);
607        }
608
609        GB_close(gb_main);
610    }
611
612    if (error) {
613        printf("Error in arb_gene_probe: %s\n", error);
614        return EXIT_FAILURE;
615    }
616
617    printf("arb_gene_probe done.\n");
618    return EXIT_SUCCESS;
619}
620
Note: See TracBrowser for help on using the repository browser.