source: tags/arb-6.0.4/PROBE/PT_io.cxx

Last change on this file was 11086, checked in by westram, 11 years ago
  • ptserver failed to build when species w/o data in the selected alignment existed
    • PT_prepare_species_sequence no longer reports an error for missing data (that error was ignored anyway)
    • delete species w/o data during '-build_clean' phase
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.3 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PT_io.cxx                                         //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11
12#include "probe.h"
13#include "pt_prototypes.h"
14#include "PT_compress.h"
15
16#include <arbdbt.h>
17#include <BI_basepos.hxx>
18#include <arb_progress.h>
19#include <arb_file.h>
20
21int compress_data(char *probestring) {
22    //! change a sequence with normal bases the PT_? format and delete all other signs
23    char    c;
24    char    *src,
25        *dest;
26    dest = src = probestring;
27
28    while ((c=*(src++)))
29    {
30        switch (c) {
31            case 'A':
32            case 'a': *(dest++) = PT_A; break;
33            case 'C':
34            case 'c': *(dest++) = PT_C; break;
35            case 'G':
36            case 'g': *(dest++) = PT_G; break;
37            case 'U':
38            case 'u':
39            case 'T':
40            case 't': *(dest++) = PT_T; break;
41            case 'N':
42            case 'n': *(dest++) = PT_N; break;
43            default: break;
44        }
45
46    }
47    *dest = PT_QU;
48    return 0;
49}
50
51ARB_ERROR probe_read_data_base(const char *name, bool readOnly) { // goes to header: __ATTR__USERESULT
52    ARB_ERROR error;
53    GB_set_verbose();
54
55    psg.gb_shell = new GB_shell;
56
57    if (!readOnly && !GB_is_writeablefile(name)) {
58        error = GBS_global_string("Database '%s' is write-protected - aborting", name);
59    }
60    if (!error) {
61        GBDATA *gb_main     = GB_open(name, readOnly ? "r" : "rw");
62        if (!gb_main) error = GB_await_error();
63        else {
64            error = GB_begin_transaction(gb_main);
65            if (!error) psg.gb_main = gb_main;
66            error = GB_end_transaction(gb_main, error);
67        }
68    }
69    return error;
70}
71
72uchar PT_compressed::translate[256]; 
73bool  PT_compressed::translation_initialized = false;
74
75#if defined(COUNT_COMPRESSES_BASES)
76BaseCounter PT_compressed::base_counter;
77#endif
78
79size_t probe_compress_sequence(char *seq, size_t seqsize) {
80    // translates a readable sequence into PT_base
81    // (see also: probe_2_readable)
82
83    PT_compressed compressed(seqsize);
84
85    compressed.createFrom(reinterpret_cast<unsigned char*>(seq), seqsize);
86    pt_assert(compressed.get_size() <= (seqsize+1));
87
88    memcpy(seq, compressed.get_seq(), compressed.get_size());
89    return compressed.get_size();
90}
91
92char *readable_probe(const char *compressed_probe, size_t len, char T_or_U) {
93    static SmartMallocPtr(uchar) smart_tab;
94    uchar *tab = NULL;
95
96    if (smart_tab.isNull()) {
97        tab = (uchar *) malloc(256);
98        memset(tab, '?', 256);
99
100        tab[PT_A]  = 'A';
101        tab[PT_C]  = 'C';
102        tab[PT_G]  = 'G';
103        tab[PT_QU] = '.';
104        tab[PT_N]  = 'N';
105
106        tab[PT_B_UNDEF] = '!';
107
108        smart_tab = tab;
109    }
110
111    tab = &*smart_tab;
112    tab[PT_T] = T_or_U;
113   
114    char *result = (char*)malloc(len+1);
115    for (size_t i = 0; i<len; ++i) {
116        result[i] = tab[safeCharIndex(compressed_probe[i])];
117    }
118    result[len] = 0;
119    return result;
120}
121
122inline GBDATA *expect_entry(GBDATA *gb_species, const char *entry_name) {
123    GBDATA *gb_entry = GB_entry(gb_species, entry_name);
124    if (!gb_entry) {
125        GB_export_errorf("Expected entry '%s' is missing for species '%s'",
126                         entry_name, GBT_read_name(gb_species));
127    }
128    return gb_entry;
129}
130
131cache::Cache<SmartCharPtr>                  probe_input_data::seq_cache(1);     // resized later
132cache::Cache<probe_input_data::SmartIntPtr> probe_input_data::rel2abs_cache(1); // resized later
133
134GB_ERROR probe_input_data::init(GBDATA *gb_species_) {
135    GBDATA *gb_cs      = expect_entry(gb_species_, "cs");
136    GBDATA *gb_compr   = expect_entry(gb_species_, "compr");
137    GBDATA *gb_baseoff = expect_entry(gb_species_, "baseoff");
138
139    GB_ERROR error = NULL;
140    if (!gb_cs || !gb_compr || !gb_baseoff) error = GB_await_error();
141    else {
142        gb_species = gb_species_;
143        size       = GB_read_count(gb_compr);
144    }
145
146    return error;
147}
148
149inline GB_ERROR PT_prepare_species_sequence(GBDATA *gb_species, const char *alignment_name, bool& data_missing, PT_compressed& compressed) {
150    GB_ERROR  error   = NULL;
151    GBDATA   *gb_ali  = GB_entry(gb_species, alignment_name);
152    GBDATA   *gb_data = gb_ali ? GB_entry(gb_ali, "data") : NULL;
153
154    data_missing = false;
155
156    if (!gb_data) {
157        data_missing = true;
158    }
159    else {
160        const char *seq = GB_read_char_pntr(gb_data);
161
162        if (!seq) {
163            error = GBS_global_string("Could not read data in '%s' for species '%s'\n(Reason: %s)",
164                                      psg.alignment_name, GBT_read_name(gb_species), GB_await_error());
165        }
166        else {
167            size_t seqlen = GB_read_string_count(gb_data);
168            compressed.createFrom(seq, seqlen);
169
170            {
171                uint32_t  checksum = GB_checksum(seq, seqlen, 1, ".-");
172                GBDATA   *gb_cs    = GB_create(gb_species, "cs", GB_INT);
173                error              = gb_cs
174                    ? GB_write_int(gb_cs, int32_t(checksum))
175                    : GB_await_error();
176            }
177
178            if (!error) {
179                GBDATA *gb_compr = GB_create(gb_species, "compr", GB_BYTES);
180                error            = gb_compr
181                    ? GB_write_bytes(gb_compr, compressed.get_seq(), compressed.get_size())
182                    : GB_await_error();
183            }
184
185            if (!error) {
186                GBDATA *gb_baseoff = GB_create(gb_species, "baseoff", GB_INTS);
187                error = gb_baseoff
188                    ? GB_write_ints(gb_baseoff, compressed.get_offsets(), compressed.get_size())
189                    : GB_await_error();
190            }
191
192            if (!error) error = GB_delete(gb_ali); // delete original seq data
193        }
194    }
195
196    return error;
197}
198
199GB_ERROR PT_prepare_data(GBDATA *gb_main) {
200    GB_ERROR  error           = GB_begin_transaction(gb_main);
201    GBDATA   *gb_species_data = GBT_get_species_data(gb_main);
202
203    if (!gb_species_data) {
204        error = GB_await_error();
205    }
206    else {
207        int icount = GB_number_of_subentries(gb_species_data);
208        int data_missing = 0;
209
210        char *ali_name = GBT_get_default_alignment(gb_main);
211        long  ali_len  = GBT_get_alignment_len(gb_main, ali_name);
212
213        PT_compressed compressBuffer(ali_len);
214
215        printf("Database contains %i species\n", icount);
216        {
217            {
218                arb_progress progress("Preparing sequence data", icount);
219                for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data);
220                     gb_species && !error;
221                    )
222                {
223                    GBDATA *gb_next = GBT_next_species(gb_species);
224                    bool    no_data;
225
226                    error = PT_prepare_species_sequence(gb_species, ali_name, no_data, compressBuffer);
227                    if (no_data) {
228                        pt_assert(!error);
229                        data_missing++;
230                        error = GB_delete(gb_species);
231                    }
232                    progress.inc();
233                    gb_species = gb_next;
234                }
235                if (error) progress.done();
236            }
237
238            if (!error) {
239                char   *master_data_name  = GBS_global_string_copy("%s/@master_data", GB_SYSTEM_FOLDER);
240                GBDATA *gb_master_data    = GB_search(gb_main, master_data_name, GB_FIND);
241                if (gb_master_data) error = GB_delete(gb_master_data);
242                free(master_data_name);
243            }
244        }
245        if (data_missing) {
246            printf("\n%i species were ignored because of missing data.\n", data_missing);
247        }
248        else {
249            printf("\nAll species contain data in alignment '%s'.\n", ali_name);
250        }
251        fflush_all();
252        free(ali_name);
253    }
254
255    error = GB_end_transaction(gb_main, error);
256    return error;
257}
258
259GB_ERROR PT_init_input_data() {
260    // reads sequence data into psg.data
261
262    GB_begin_transaction(psg.gb_main);
263
264    // read ref SAI (e.g. ecoli)
265    {
266        char   *def_ref     = GBT_get_default_ref(psg.gb_main);
267        GBDATA *gb_sai_data = GBT_get_SAI_data(psg.gb_main);
268        GBDATA *gb_ref      = GBT_find_SAI_rel_SAI_data(gb_sai_data, def_ref);
269
270        psg.ecoli = 0;
271        if (gb_ref) {
272            GBDATA *gb_data = GBT_read_sequence(gb_ref, psg.alignment_name);
273            if (gb_data) {
274                psg.ecoli = GB_read_string(gb_data);
275            }
276        }
277        free(def_ref);
278    }
279
280    GBDATA *gb_species_data = GBT_get_species_data(psg.gb_main);
281    int     icount          = GB_number_of_subentries(gb_species_data);
282
283    psg.data = new probe_input_data[icount];
284    psg.data_count = 0;
285
286    printf("Database contains %i species\n", icount);
287
288    GB_ERROR error = NULL;
289    {
290        arb_progress progress("Checking data", icount);
291        int count = 0;
292
293        for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data);
294             gb_species;
295             gb_species = GBT_next_species(gb_species))
296        {
297            probe_input_data& pid = psg.data[count];
298
299            error = pid.init(gb_species);
300            if (error) break;
301            count++;
302            progress.inc();
303        }
304
305        psg.data_count = count;
306        GB_commit_transaction(psg.gb_main);
307
308        if (error) progress.done();
309    }
310
311    fflush_all();
312    return error;
313}
314
315void PT_build_species_hash() {
316    long i;
317    psg.namehash = GBS_create_hash(psg.data_count, GB_MIND_CASE);
318    for (i=0; i<psg.data_count; i++) {
319        GBS_write_hash(psg.namehash, psg.data[i].get_shortname(), i+1);
320    }
321    unsigned int    max_size;
322    max_size = 0;
323    for (i = 0; i < psg.data_count; i++) {  // get max sequence len
324        max_size = std::max(max_size, (unsigned)(psg.data[i].get_size()));
325        psg.char_count += psg.data[i].get_size();
326    }
327    psg.max_size = max_size;
328
329    if (psg.ecoli) {
330        BI_ecoli_ref *ref = new BI_ecoli_ref;
331        ref->init(psg.ecoli, strlen(psg.ecoli));
332        psg.bi_ecoli = ref;
333    }
334}
335
336
337long PT_abs_2_ecoli_rel(long pos) {
338    if (!psg.ecoli) return pos;
339    return psg.bi_ecoli->abs_2_rel(pos);
340}
341
342// --------------------------------------------------------------------------------
343
344#ifdef UNIT_TESTS
345#ifndef TEST_UNIT_H
346#include <test_unit.h>
347#endif
348
349inline int *intcopy(int i) { int *ip = new int; *ip = i; return ip; }
350
351#define CACHED(p,t) that(p.is_cached()).is_equal_to(t)
352#define CACHED_p123(t1,t2,t3) all().of(CACHED(p1, t1), CACHED(p2, t2), CACHED(p3, t3))
353
354void TEST_CachedPtr() {
355    using namespace cache;
356    {
357        typedef SmartPtr<int> IntPtr;
358
359        CacheHandle<IntPtr> p1;
360        CacheHandle<IntPtr> p2;
361        CacheHandle<IntPtr> p3;
362
363        {
364            Cache<IntPtr> cache(2);
365            TEST_EXPECT_ZERO(cache.entries()); // nothing cached yet
366
367            p1.assign(intcopy(1), cache);
368            TEST_EXPECT_EQUAL(*p1.access(cache), 1);
369            TEST_EXPECT_EQUAL(cache.entries(), 1);
370            TEST_EXPECTATION(CACHED_p123(true, false, false));
371
372            p2.assign(intcopy(2), cache);
373            TEST_EXPECT_EQUAL(*p2.access(cache), 2);
374            TEST_EXPECT_EQUAL(cache.entries(), 2);
375            TEST_EXPECTATION(CACHED_p123(true, true, false));
376
377            p3.assign(intcopy(3), cache);
378            TEST_EXPECT_EQUAL(*p3.access(cache), 3);
379            TEST_EXPECT_EQUAL(cache.entries(), 2);
380            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p1 has been invalidated by caching p3
381
382            p3.assign(intcopy(33), cache); // test re-assignment
383            TEST_EXPECT_EQUAL(*p3.access(cache), 33);
384            TEST_EXPECT_EQUAL(cache.entries(), 2);
385            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p2 still cached
386
387            TEST_EXPECT_EQUAL(*p2.access(cache), 2);          // should make p2 the LRU cache entry
388            TEST_EXPECT_EQUAL(cache.entries(), 2);
389            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p2 still cached
390
391            IntPtr s4;
392            {
393                CacheHandle<IntPtr> p4;
394                p4.assign(intcopy(4), cache);
395                TEST_EXPECT_EQUAL(*p4.access(cache), 4);
396                TEST_EXPECT_EQUAL(cache.entries(), 2);
397
398                s4 = p4.access(cache); // keep data of p4 in s4
399                TEST_EXPECT_EQUAL(s4.references(), 2); // ref'd by s4 and p4
400
401                p4.release(cache); // need to release p4 before destruction (otherwise assertion fails)
402            }
403
404            TEST_EXPECT_EQUAL(*s4, 4);             // check kept value of deleted CacheHandle
405            TEST_EXPECT_EQUAL(s4.references(), 1); // only ref'd by s4
406
407            TEST_EXPECT_EQUAL(cache.entries(), 1);             // contains only p2 (p4 has been released)
408            TEST_EXPECTATION(CACHED_p123(false, true, false)); // p3 has been invalidated by caching p4
409        }
410
411        TEST_EXPECTATION(CACHED_p123(false, false, false)); // Cache was destroyed = > all CacheHandle will be invalid
412        // no need to release p1..p3 (due cache was destroyed)
413    }
414
415    // test cache of SmartCharPtr
416    {
417        Cache<SmartCharPtr> cache(3);
418
419        {
420            const int P = 4;
421            CacheHandle<SmartCharPtr> p[P];
422            const char *word[] = { "apple", "orange", "pie", "juice" };
423
424            for (int i = 0; i<P; ++i) p[i].assign(strdup(word[i]), cache);
425            TEST_REJECT(p[0].is_cached());
426            for (int i = 1; i<P; ++i) TEST_EXPECT_EQUAL(&*p[i].access(cache), word[i]);
427
428            TEST_REJECT(p[0].is_cached());
429            TEST_EXPECT(p[1].is_cached()); // oldest entry
430
431            cache.resize(cache.size()-1);
432            TEST_REJECT(p[1].is_cached()); // invalidated by resize
433
434            for (int i = P-1; i >= 0; --i) p[i].assign(strdup(word[P-1-i]), cache);
435
436            for (int i = 0; i<2; ++i) TEST_EXPECT_EQUAL(&*p[i].access(cache), word[P-1-i]);
437            for (int i = 2; i<P; ++i) TEST_REJECT(p[i].is_cached());
438
439            cache.flush();
440        }
441    }
442}
443
444#endif // UNIT_TESTS
445
446// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.