source: tags/ms_r16q2/GENOM_IMPORT/DBwriter.cxx

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