source: tags/ms_r16q2/PROBE/PT_io.cxx

Last change on this file was 13970, checked in by westram, 9 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 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        if (!seq) {
162            error = GBS_global_string("Could not read data in '%s' for species '%s'\n(Reason: %s)",
163                                      alignment_name, GBT_read_name(gb_species), GB_await_error());
164        }
165        else {
166            size_t seqlen = GB_read_string_count(gb_data);
167            if (seqlen>compressed.get_allowed_size()) {
168                error = GBS_global_string("Sequence too long in '%s' of '%s'\n(Hint: format alignment to fix this problem)",
169                                          alignment_name, GBT_read_name(gb_species));
170            }
171
172            if (!error) {
173                compressed.createFrom(seq, seqlen);
174                {
175                    uint32_t  checksum = GB_checksum(seq, seqlen, 1, ".-");
176                    GBDATA   *gb_cs    = GB_create(gb_species, "cs", GB_INT);
177                    error              = gb_cs
178                        ? GB_write_int(gb_cs, int32_t(checksum))
179                        : GB_await_error();
180                }
181            }
182
183            if (!error) {
184                GBDATA *gb_compr = GB_create(gb_species, "compr", GB_BYTES);
185                error            = gb_compr
186                    ? GB_write_bytes(gb_compr, compressed.get_seq(), compressed.get_size())
187                    : GB_await_error();
188            }
189
190            if (!error) {
191                GBDATA *gb_baseoff = GB_create(gb_species, "baseoff", GB_INTS);
192                error = gb_baseoff
193                    ? GB_write_ints(gb_baseoff, compressed.get_offsets(), compressed.get_size())
194                    : GB_await_error();
195            }
196
197            if (!error) error = GB_delete(gb_ali); // delete original seq data
198        }
199    }
200
201    return error;
202}
203
204GB_ERROR PT_prepare_data(GBDATA *gb_main) {
205    GB_ERROR  error           = GB_begin_transaction(gb_main);
206    GBDATA   *gb_species_data = GBT_get_species_data(gb_main);
207
208    if (!gb_species_data) {
209        error = GB_await_error();
210    }
211    else {
212        int icount = GB_number_of_subentries(gb_species_data);
213        int data_missing = 0;
214
215        char *ali_name = GBT_get_default_alignment(gb_main);
216        long  ali_len  = GBT_get_alignment_len(gb_main, ali_name);
217
218        PT_compressed compressBuffer(ali_len);
219
220        printf("Database contains %i species\n", icount);
221        {
222            {
223                arb_progress progress("Preparing sequence data", icount);
224                for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data);
225                     gb_species && !error;
226                    )
227                {
228                    GBDATA *gb_next = GBT_next_species(gb_species);
229                    bool    no_data;
230
231                    error = PT_prepare_species_sequence(gb_species, ali_name, no_data, compressBuffer);
232                    if (no_data) {
233                        pt_assert(!error);
234                        data_missing++;
235                        error = GB_delete(gb_species);
236                    }
237                    progress.inc();
238                    gb_species = gb_next;
239                }
240                if (error) progress.done();
241            }
242
243            if (!error) {
244                char   *master_data_name  = GBS_global_string_copy("%s/@master_data", GB_SYSTEM_FOLDER);
245                GBDATA *gb_master_data    = GB_search(gb_main, master_data_name, GB_FIND);
246                if (gb_master_data) error = GB_delete(gb_master_data);
247                free(master_data_name);
248            }
249        }
250        if (data_missing) {
251            printf("\n%i species were ignored because of missing data.\n", data_missing);
252        }
253        else {
254            printf("\nAll species contain data in alignment '%s'.\n", ali_name);
255        }
256        fflush_all();
257        free(ali_name);
258    }
259
260    error = GB_end_transaction(gb_main, error);
261    return error;
262}
263
264GB_ERROR PT_init_input_data() {
265    // reads sequence data into psg.data
266
267    GB_begin_transaction(psg.gb_main);
268
269    // read ref SAI (e.g. ecoli)
270    {
271        char   *def_ref     = GBT_get_default_ref(psg.gb_main);
272        GBDATA *gb_sai_data = GBT_get_SAI_data(psg.gb_main);
273        GBDATA *gb_ref      = GBT_find_SAI_rel_SAI_data(gb_sai_data, def_ref);
274
275        psg.ecoli = 0;
276        if (gb_ref) {
277            GBDATA *gb_data = GBT_find_sequence(gb_ref, psg.alignment_name);
278            if (gb_data) {
279                psg.ecoli = GB_read_string(gb_data);
280            }
281        }
282        free(def_ref);
283    }
284
285    GBDATA *gb_species_data = GBT_get_species_data(psg.gb_main);
286    int     icount          = GB_number_of_subentries(gb_species_data);
287
288    psg.data = new probe_input_data[icount];
289    psg.data_count = 0;
290
291    printf("Database contains %i species\n", icount);
292
293    GB_ERROR error = NULL;
294    {
295        arb_progress progress("Checking data", icount);
296        int count = 0;
297
298        for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data);
299             gb_species;
300             gb_species = GBT_next_species(gb_species))
301        {
302            probe_input_data& pid = psg.data[count];
303
304            error = pid.init(gb_species);
305            if (error) break;
306            count++;
307            progress.inc();
308        }
309
310        psg.data_count = count;
311        GB_commit_transaction(psg.gb_main);
312
313        if (error) progress.done();
314    }
315
316    fflush_all();
317    return error;
318}
319
320void PT_build_species_hash() {
321    long i;
322    psg.namehash = GBS_create_hash(psg.data_count, GB_MIND_CASE);
323    for (i=0; i<psg.data_count; i++) {
324        GBS_write_hash(psg.namehash, psg.data[i].get_shortname(), i+1);
325    }
326    unsigned int    max_size;
327    max_size = 0;
328    for (i = 0; i < psg.data_count; i++) {  // get max sequence len
329        max_size = std::max(max_size, (unsigned)(psg.data[i].get_size()));
330        psg.char_count += psg.data[i].get_size();
331    }
332    psg.max_size = max_size;
333
334    if (psg.ecoli) {
335        BI_ecoli_ref *ref = new BI_ecoli_ref;
336        ref->init(psg.ecoli, strlen(psg.ecoli));
337        psg.bi_ecoli = ref;
338    }
339}
340
341
342long PT_abs_2_ecoli_rel(long pos) {
343    if (!psg.ecoli) return pos;
344    return psg.bi_ecoli->abs_2_rel(pos);
345}
346
347// --------------------------------------------------------------------------------
348
349#ifdef UNIT_TESTS
350#ifndef TEST_UNIT_H
351#include <test_unit.h>
352#endif
353
354inline int *intcopy(int i) { int *ip = new int; *ip = i; return ip; }
355
356#define CACHED(p,t) that(p.is_cached()).is_equal_to(t)
357#define CACHED_p123(t1,t2,t3) all().of(CACHED(p1, t1), CACHED(p2, t2), CACHED(p3, t3))
358
359void TEST_CachedPtr() {
360    using namespace cache;
361    {
362        typedef SmartPtr<int> IntPtr;
363
364        CacheHandle<IntPtr> p1;
365        CacheHandle<IntPtr> p2;
366        CacheHandle<IntPtr> p3;
367
368        {
369            Cache<IntPtr> cache(2);
370            TEST_EXPECT_ZERO(cache.entries()); // nothing cached yet
371
372            p1.assign(intcopy(1), cache);
373            TEST_EXPECT_EQUAL(*p1.access(cache), 1);
374            TEST_EXPECT_EQUAL(cache.entries(), 1);
375            TEST_EXPECTATION(CACHED_p123(true, false, false));
376
377            p2.assign(intcopy(2), cache);
378            TEST_EXPECT_EQUAL(*p2.access(cache), 2);
379            TEST_EXPECT_EQUAL(cache.entries(), 2);
380            TEST_EXPECTATION(CACHED_p123(true, true, false));
381
382            p3.assign(intcopy(3), cache);
383            TEST_EXPECT_EQUAL(*p3.access(cache), 3);
384            TEST_EXPECT_EQUAL(cache.entries(), 2);
385            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p1 has been invalidated by caching p3
386
387            p3.assign(intcopy(33), cache); // test re-assignment
388            TEST_EXPECT_EQUAL(*p3.access(cache), 33);
389            TEST_EXPECT_EQUAL(cache.entries(), 2);
390            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p2 still cached
391
392            TEST_EXPECT_EQUAL(*p2.access(cache), 2);          // should make p2 the LRU cache entry
393            TEST_EXPECT_EQUAL(cache.entries(), 2);
394            TEST_EXPECTATION(CACHED_p123(false, true, true)); // p2 still cached
395
396            IntPtr s4;
397            {
398                CacheHandle<IntPtr> p4;
399                p4.assign(intcopy(4), cache);
400                TEST_EXPECT_EQUAL(*p4.access(cache), 4);
401                TEST_EXPECT_EQUAL(cache.entries(), 2);
402
403                s4 = p4.access(cache); // keep data of p4 in s4
404                TEST_EXPECT_EQUAL(s4.references(), 2); // ref'd by s4 and p4
405
406                p4.release(cache); // need to release p4 before destruction (otherwise assertion fails)
407            }
408
409            TEST_EXPECT_EQUAL(*s4, 4);             // check kept value of deleted CacheHandle
410            TEST_EXPECT_EQUAL(s4.references(), 1); // only ref'd by s4
411
412            TEST_EXPECT_EQUAL(cache.entries(), 1);             // contains only p2 (p4 has been released)
413            TEST_EXPECTATION(CACHED_p123(false, true, false)); // p3 has been invalidated by caching p4
414        }
415
416        TEST_EXPECTATION(CACHED_p123(false, false, false)); // Cache was destroyed = > all CacheHandle will be invalid
417        // no need to release p1..p3 (due cache was destroyed)
418    }
419
420    // test cache of SmartCharPtr
421    {
422        Cache<SmartCharPtr> cache(3);
423
424        {
425            const int P = 4;
426            CacheHandle<SmartCharPtr> p[P];
427            const char *word[] = { "apple", "orange", "pie", "juice" };
428
429            for (int i = 0; i<P; ++i) p[i].assign(strdup(word[i]), cache);
430            TEST_REJECT(p[0].is_cached());
431            for (int i = 1; i<P; ++i) TEST_EXPECT_EQUAL(&*p[i].access(cache), word[i]);
432
433            TEST_REJECT(p[0].is_cached());
434            TEST_EXPECT(p[1].is_cached()); // oldest entry
435
436            cache.resize(cache.size()-1);
437            TEST_REJECT(p[1].is_cached()); // invalidated by resize
438
439            for (int i = P-1; i >= 0; --i) p[i].assign(strdup(word[P-1-i]), cache);
440
441            for (int i = 0; i<2; ++i) TEST_EXPECT_EQUAL(&*p[i].access(cache), word[P-1-i]);
442            for (int i = 2; i<P; ++i) TEST_REJECT(p[i].is_cached());
443
444            cache.flush();
445        }
446    }
447}
448
449#endif // UNIT_TESTS
450
451// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.