source: branches/port5/GENOM_IMPORT/DBwriter.cxx

Last change on this file was 5853, checked in by westram, 16 years ago
  • followup to [5852] (more to follow)
    • GB_get_error() → GB_await_error() where appropriate (GENOM,GENOM_IMPORT)
    • fixed error handling in several functions
File size: 15.6 KB
Line 
1// ================================================================ //
2//                                                                  //
3//   File      : DBwriter.cxx                                       //
4//   Purpose   :                                                    //
5//                                                                  //
6//   Coded by Ralf Westram (coder@reallysoft.de) in November 2006   //
7//   Institute of Microbiology (Technical University Munich)        //
8//   http://www.arb-home.de/                                        //
9//                                                                  //
10// ================================================================ //
11#include "DBwriter.h"
12
13#define AW_RENAME_SKIP_GUI
14#include <algorithm>
15#include <AW_rename.hxx>
16#include <arbdbt.h>
17#include <adGene.h>
18#include <awt_translate.hxx>
19#include <aw_question.hxx>
20#include <GEN.hxx>
21
22#define ARB_GENE_REF "ARB_is_gene"
23
24using namespace std;
25
26typedef SmartCustomPtr(GEN_position, GEN_free_position) GEN_positionPtr;
27
28// --------------------------------------------------------------------------------
29
30void DBerror::init(const string& msg, GB_ERROR gberror) {
31    gi_assert(gberror);
32    if (gberror) err = msg+" (Reason: "+gberror+")";
33    else err = msg; // ndebug!
34}
35
36DBerror::DBerror() { init("", GB_await_error()); }
37DBerror::DBerror(const char *msg) { string errmsg(msg); init(errmsg, GB_await_error()); }
38DBerror::DBerror(const string& msg) { init(msg, GB_await_error()); }
39DBerror::DBerror(const char *msg, GB_ERROR gberror) { string errmsg(msg); init(errmsg, gberror); }
40DBerror::DBerror(const string& msg, GB_ERROR gberror) { init(msg, gberror); }
41
42// --------------------------------------------------------------------------------
43// DB access functions, throwing DBerror on failure
44
45static GBDATA *DB_create_container(GBDATA *parent, const char *name, bool mark) {
46    // create container (optionally set mark flag)
47    GBDATA *gb_container = GB_create_container(parent, name);
48    if (!gb_container) throw DBerror(GBS_global_string("Failed to create container '%s'", name));
49
50    if (mark) {
51        GB_ERROR err = GB_write_flag(gb_container, 1);
52        if (err) throw DBerror(GBS_global_string("Failed to mark %s", name), err);
53    }
54
55    return gb_container;
56}
57
58static GBDATA *DB_create_string_field(GBDATA *parent, const char *field, const char *content) {
59    // create field with content
60
61    gi_assert(content[0]);
62    // do NOT WRITE empty string-fields into ARB DB,
63    // cause ARB DB does not differ between empty content and non-existing fields
64    // (i.e. when writing an empty string, ARB removes the field)
65
66    GBDATA *gb_field = GB_create(parent, field, GB_STRING);
67    if (!gb_field) throw DBerror(GBS_global_string("Failed to create field '%s'", field));
68
69    GB_ERROR err = GB_write_string(gb_field, content);
70    if (err) throw DBerror(GBS_global_string("Failed to write to field '%s'", field), err);
71
72    return gb_field;
73}
74
75static GBDATA *DB_create_byte_field(GBDATA *parent, const char *field, unsigned char content) {
76    // create field with content
77    GBDATA *gb_field = GB_create(parent, field, GB_BYTE);
78    if (!gb_field) throw DBerror(GBS_global_string("Failed to create field '%s'", field));
79
80    GB_ERROR err = GB_write_byte(gb_field, content);
81    if (err) throw DBerror(GBS_global_string("Failed to write to field '%s'", field), err);
82
83    return gb_field;
84}
85
86// --------------------------------------------------------------------------------
87
88void DBwriter::createOrganism(const string& flatfile, const char *importerTag)
89{
90    gi_assert(!gb_organism && !gb_gene_data);
91
92    // create the organism
93    {
94        UniqueNameDetector& UND_species = *session.und_species;
95
96        char *organism_name = AWTC_makeUniqueShortName("genome", UND_species);
97        if (!organism_name) throw DBerror();
98
99        gb_organism = DB_create_container(session.gb_species_data, "species", true);
100        DB_create_string_field(gb_organism, "name", organism_name);
101
102        UND_species.add_name(organism_name);
103        free(organism_name);
104    }
105
106    // store info about source
107    DB_create_string_field(gb_organism, "ARB_imported_from", flatfile.c_str());
108    DB_create_string_field(gb_organism, "ARB_imported_format", importerTag);
109}
110
111class NoCaseCmp {
112    static bool less_nocase(char c1, char c2) { return toupper(c1) < toupper(c2); }
113    static bool nonequal_nocase(char c1, char c2) { return toupper(c1) != toupper(c2); }
114public:
115    NoCaseCmp() {}
116
117    bool operator()(const string& s1, const string& s2) {
118        return lexicographical_compare(s1.begin(), s1.end(),
119                                       s2.begin(), s2.end(),
120                                       less_nocase);
121    }
122
123    static bool has_prefix(const string& s1, const string& s2) {
124        // return true if s1 starts with prefix s2
125        return
126            s1.length() >= s2.length() &&
127            !lexicographical_compare(s1.begin(), s1.begin()+s2.length(),
128                                     s2.begin(), s2.end(),
129                                     nonequal_nocase);
130    }
131};
132
133typedef map<string, string, NoCaseCmp> TranslateMap;
134struct Translator { TranslateMap trans; };
135
136Translator *DBwriter::unreserve = 0;
137
138const string& DBwriter::getUnreservedQualifier(const string& qualifier) {
139    // return a non-reserved qualifier
140    // (some are reserved - e.g. name, pos_start, ... and all starting with 'ARB_')
141    // if a qualifier is reserved, 'ORG_' is prepended.
142    //
143    // (Note: When we'll export data, 'ORG_' shall be removed from qualifiers!)
144
145    static string prefix = "ORG_";
146
147    if (!unreserve) {
148        unreserve = new Translator;
149        const char *reserved[] = {
150            "name", "type",
151            "pos_start", "pos_stop", "pos_complement", "pos_certain", "pos_joined",
152            0
153        };
154
155        for (int i = 0; reserved[i]; ++i) {
156            unreserve->trans[reserved[i]] = prefix+reserved[i];
157        }
158    }
159
160    TranslateMap::const_iterator found = unreserve->trans.find(qualifier);
161    if (found != unreserve->trans.end()) { // qualifier is reserved
162        return found->second;
163    }
164    if (NoCaseCmp::has_prefix(qualifier, prefix) || // qualifier starts with 'ORG_'
165        NoCaseCmp::has_prefix(qualifier, "ARB_")) // or 'ARB_'
166    {
167        unreserve->trans[qualifier] = prefix+qualifier; // add as 'ORG_ORG_' or 'ORG_ARB_' to TranslateMap
168        return unreserve->trans[qualifier];
169    }
170    return qualifier;
171}
172
173void DBwriter::writeFeature(const Feature& feature)
174{
175    gi_assert(gb_organism);
176    if (!gb_gene_data) {
177        gb_gene_data = DB_create_container(gb_organism, "gene_data", false);
178        generatedGenes.clear();
179    }
180
181    // create new gene
182    GBDATA *gb_gene;
183    {
184        string                gene_name = feature.createGeneName();
185        NameCounter::iterator existing  = generatedGenes.find(gene_name);
186        if (existing == generatedGenes.end()) { // first occurrence of that gene name
187            generatedGenes[gene_name] = 1;
188        }
189        else {
190            existing->second++; // increment occurrences
191        }
192
193        gb_gene = DB_create_container(gb_gene_data, "gene", false);
194        DB_create_string_field(gb_gene, "name", gene_name.c_str());
195
196        string type = feature.getType();
197        DB_create_string_field(gb_gene, "type", type.c_str());
198
199        if (type == "source") {
200            DB_create_byte_field(gb_gene, ARB_HIDDEN, 1);
201        }
202    }
203
204    // store location
205    {
206        GEN_positionPtr pos = feature.getLocation().create_GEN_position();
207        GB_ERROR        err = GEN_write_position(gb_gene, &*pos);
208        if (err) throw DBerror("Failed to write location", err);
209    }
210
211    // store qualifiers
212    {
213        const stringMap& qualifiers = feature.getQualifiers();
214        stringMapCIter   e          = qualifiers.end();
215       
216        for (stringMapCIter i = qualifiers.begin(); i != e; ++i) {
217            const string& unreserved = getUnreservedQualifier(i->first);
218            DB_create_string_field(gb_gene, unreserved.c_str(), i->second.c_str());
219        }
220    }
221}
222
223void DBwriter::writeSequence(const SequenceBuffer& seqData)
224{
225    gi_assert(gb_organism);
226    GBDATA *gb_data = GBT_add_data(gb_organism, ali_name, "data", GB_STRING);
227    if (!gb_data) throw DBerror("Failed to create alignment");
228
229    GB_ERROR err = GB_write_string(gb_data, seqData.getSequence());
230    if (err) throw DBerror("Failed to write alignment", err);
231}
232
233void DBwriter::renumberDuplicateGenes() {
234    NameCounter renameCounter; // for those genes which get renumbered, count upwards here
235    {
236        NameCounter::iterator gg_end = generatedGenes.end();
237        for (NameCounter::iterator gg = generatedGenes.begin(); gg != gg_end; ++gg) {
238            if (gg->second > 1) renameCounter[gg->first] = 0;
239        }
240    }
241    NameCounter::iterator rc_end = renameCounter.end();
242
243    for (GBDATA *gb_gene = GB_entry(gb_gene_data, "gene");
244         gb_gene;
245         gb_gene = GB_nextEntry(gb_gene))
246    {
247        GBDATA *gb_name   = GB_entry(gb_gene, "name");
248        string  gene_name = GB_read_char_pntr(gb_name);
249
250        NameCounter::iterator rc = renameCounter.find(gene_name);
251        if (rc != rc_end) { // rename current_gene
252            int maxOccurrences = generatedGenes[gene_name];
253            rc->second++;   // increment occurrence counter
254            gi_assert(rc->second <= maxOccurrences);
255
256            int digits = strlen(GBS_global_string("%i", maxOccurrences));
257            gene_name    += GBS_global_string("_%0*i", digits, rc->second);
258            GB_ERROR err  = GB_write_string(gb_name, gene_name.c_str());
259            if (err) throw DBerror("Failed to write to field 'name' (during gene-renumbering)", err);
260        }
261    }
262}
263
264static void importerWarning(AW_CL cl_importer, const char *message) {
265    Importer *importer = reinterpret_cast<Importer*>(cl_importer);
266    importer->warning(message);
267}
268
269void DBwriter::testAndRemoveTranslations(Importer& importer) {
270    GEN_testAndRemoveTranslations(gb_gene_data, importerWarning, reinterpret_cast<AW_CL>(&importer), session.ok_to_ignore_wrong_start_codon);
271}
272
273// ----------------------------------------------
274//      hide duplicated genes (from genebank)
275
276inline bool operator<(const GEN_positionPtr& A, const GEN_positionPtr& B) {
277    const GEN_position& a = *A;
278    const GEN_position& b = *B;
279
280    int cmp = int(a.start_pos[0]) - int(b.start_pos[0]);
281    if (!cmp) {
282        cmp = int(a.stop_pos[a.parts-1]) - int(b.stop_pos[b.parts-1]);
283        if (!cmp) {
284            cmp = int(a.parts) - int(b.parts); // less parts is <
285            if (!cmp) {
286                for (int p = 0; p<a.parts; ++p) { // compare all parts
287                    cmp = int(b.complement[p]) - int(a.complement[p]); if (cmp) break; // non-complement is <
288                    cmp = int(a.start_pos[p])  - int(b.start_pos[p]);  if (cmp) break;
289                    cmp = int(a.stop_pos[p])   - int(b.stop_pos[p]);   if (cmp) break;
290                }
291            }
292        }
293    }
294
295    return cmp<0;
296}
297
298class PosGene {
299    GBDATA                  *gb_Gene;
300    mutable GEN_positionPtr  pos;
301
302public:
303    PosGene(GBDATA *gb_gene) : gb_Gene(gb_gene) {
304        GEN_position *pp = GEN_read_position(gb_gene);
305        if (!pp) {
306            throw GBS_global_string("Can't read gene position (Reason: %s)", GB_await_error());
307        }
308        pos = pp;
309    }
310
311    // const GEN_position *getPosition() const { return &*pos; }
312    const GEN_positionPtr& getPosition() const { return pos; }
313
314    const char *getName() const {
315        GBDATA *gb_name = GB_entry(gb_Gene, "name");
316        gi_assert(gb_name);
317        return GB_read_char_pntr(gb_name);
318    }
319    const char *getType() const {
320        GBDATA *gb_type = GB_entry(gb_Gene, "type");
321        gi_assert(gb_type);
322        return GB_read_char_pntr(gb_type);
323    }
324    bool hasType(const char *type) const { return strcmp(getType(), type) == 0; }
325
326    void hide() { DB_create_byte_field(gb_Gene, ARB_HIDDEN, 1); }
327    void addRefToGene(const char *name_of_gene) { DB_create_string_field(gb_Gene, ARB_GENE_REF, name_of_gene); }
328
329#if defined(DEBUG)
330    const char *description() const {
331        return GBS_global_string("%zi-%zi (%i parts)", pos->start_pos[0], pos->stop_pos[pos->parts-1], pos->parts);
332    }
333    void dump() {
334        GB_dump(gb_Gene);
335    }
336#endif // DEBUG
337
338};
339
340inline bool operator<(const PosGene& a, const PosGene& b) {
341    return a.getPosition() < b.getPosition();
342}
343inline bool operator == (const PosGene& a, const PosGene& b) {
344    return !(a<b || b<a);
345}
346
347class hasType {
348    const char *type;
349public:
350    hasType(const char *t) : type(t) {}
351    bool operator()(const PosGene& pg) { return pg.hasType(type); }
352};
353
354// typedef SmartPtr<PosGene> PosGenePtr;
355
356void DBwriter::hideUnwantedGenes() {
357    typedef vector<PosGene> Genes;
358    typedef Genes::iterator GeneIter;
359
360    GBDATA *gb_gene;
361    Genes   gps;
362
363    // read all gene positions
364    for (gb_gene = GB_entry(gb_gene_data, "gene"); gb_gene; gb_gene = GB_nextEntry(gb_gene)) {
365        gps.push_back(PosGene(gb_gene));
366    }
367
368    sort(gps.begin(), gps.end()); // sort positions
369
370    // find duplicate geness (with identical location)
371    GeneIter end        = gps.end();
372    GeneIter p          = gps.begin();
373    GeneIter firstEqual = p;
374    GeneIter lastEqual  = p;
375
376    while (firstEqual != end) {
377        ++p;
378        if (p != end && *p == *firstEqual) {
379            lastEqual = p;
380        }
381        else {                  // found different element (or last)
382            int count = lastEqual-firstEqual+1;
383
384            if (count>1) { // we have 2 or more duplicate genes
385                GeneIter equalEnd = lastEqual+1;
386                GeneIter gene     = find_if(firstEqual, equalEnd, hasType("gene")); // locate 'type' == 'gene'
387
388                if (gene != equalEnd) { // found type gene
389                    bool        hideGene  = false;
390                    GeneIter    e         = firstEqual;
391                    const char *gene_name = gene->getName();
392
393                    while (e != equalEnd) {
394                        if (e != gene) { // for all genes that have 'type' != 'gene'
395                            e->addRefToGene(gene_name); // add a reference to the gene
396                            hideGene = true; // and hide the gene
397                        }
398                        ++e;
399                    }
400
401                    if (hideGene) gene->hide();
402                }
403            }
404            firstEqual = lastEqual = p;
405        }
406    }
407}
408
409
410void DBwriter::finalizeOrganism(const MetaInfo& meta, const References& refs, Importer& importer)
411{
412    // write metadata
413    {
414        const stringMap& entries = meta.getEntries();
415        stringMapCIter   e       = entries.end();
416
417        for (stringMapCIter i = entries.begin(); i != e; ++i) {
418            const string& content = i->second;
419            if (!content.empty()) { // skip empty entries
420                DB_create_string_field(gb_organism, i->first.c_str(), content.c_str());
421            }
422        }
423    }
424
425    // write references
426    {
427        stringSet refKeys;
428        refs.getKeys(refKeys);
429
430        stringSetIter e = refKeys.end();
431        for (stringSetIter i = refKeys.begin(); i != e; ++i) {
432            DB_create_string_field(gb_organism, i->c_str(), refs.tagged_content(*i).c_str());
433        }
434    }
435
436    // finalize genes data
437    if (gb_gene_data) {
438        renumberDuplicateGenes();            // renumber genes with equal names
439        testAndRemoveTranslations(importer); // test translations and remove reproducable translations
440        hideUnwantedGenes();
441    }
442    else GB_warning("No genes have been written (missing feature table?)");
443
444    // cleanup
445    generatedGenes.clear();
446    gb_gene_data = 0;
447    gb_organism  = 0;
448}
449
450void DBwriter::deleteStaticData() {
451    if (unreserve) {
452        delete unreserve;
453        unreserve = NULL;
454    }
455}
456
457
458
Note: See TracBrowser for help on using the repository browser.