source: tags/svn.1.5.4/PROBE/PT_findEx.cxx

Last change on this file was 8041, checked in by westram, 13 years ago

merge from dev [7990] [7991] [7992] [7993] [7994] [7995] [7996] [7998] [8001] [8003] [8005] [8006] [8007] [8008] [8009] [8011] [8012] [8019]

  • added a faked ecoli to ptserver test-db
  • added unit-tests for gene-ptserver
  • rename ptserver testdb
  • ptserver db-cleanup
    • size of mapfile for SSUREF108 is reduced by 85% (3,4Gb → 519Mb)
  • ptserver
    • refactored
      • prefix tree builders
      • probe_input_data
    • removed ptnd_new_match (dead code)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.3 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PT_findEx.cxx                                     //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "probe.h"
12#include <PT_server_prototypes.h>
13#include "probe_tree.h"
14#include "pt_prototypes.h"
15#include <arb_strbuf.h>
16
17static bool findLeftmostProbe(POS_TREE *node, char *probe, int restlen, int height) {
18    if (restlen==0) return true;
19
20    switch (PT_read_type(node)) {
21        case PT_NT_NODE: {
22            for (int i=PT_A; i<PT_B_MAX; ++i) {
23                POS_TREE *son = PT_read_son(node, PT_BASES(i));
24                if (son) {
25                    probe[0] = PT_BASES(i); // write leftmost probe into result
26                    bool found = findLeftmostProbe(son, probe+1, restlen-1, height+1);
27                    pt_assert(implicated(found, strlen(probe) == (size_t)restlen));
28                    if (found) return true;
29                }
30            }
31            break;
32        }
33        case PT_NT_CHAIN: {
34            pt_assert(0);  // unhandled yet
35            break;
36        }
37        case PT_NT_LEAF: {
38            // here the probe-tree is cut off, because only one species matches
39            DataLoc loc(node);
40            int     pos  = loc.rpos + height;
41            int     name = loc.name;
42
43            if (pos + restlen >= psg.data[name].get_size())
44                break;          // at end-of-sequence -> no probe with wanted length here
45
46            pt_assert(probe[restlen] == 0);
47            const char *seq_data = psg.data[name].get_data();
48            for (int r = 0; r<restlen; ++r) {
49                int data = seq_data[pos+r];
50                if (data == PT_QU || data == PT_N) return false; // ignore probes that contain 'N' or '.'
51                probe[r] = data;
52            }
53            pt_assert(probe[restlen] == 0);
54            pt_assert(strlen(probe) == (size_t)restlen);
55            return true;
56        }
57        default: pt_assert(0); break;  // oops
58    }
59
60    return false;
61}
62//  ---------------------------------------------------------------------
63//      static bool findNextProbe(POS_TREE *node, char *probe, int restlen)
64//  ---------------------------------------------------------------------
65// searches next probe after 'probe' ('probe' itself may not exist)
66// returns: true if next probe was found
67// 'probe' is modified to next probe
68static bool findNextProbe(POS_TREE *node, char *probe, int restlen, int height) {
69    if (restlen==0) return false;  // in this case we found the recent probe
70    // returning false upwards takes the next after
71
72    switch (PT_read_type(node)) {
73        case PT_NT_NODE: {
74            POS_TREE *son   = PT_read_son(node, PT_BASES(probe[0]));
75            bool      found = (son != 0) && findNextProbe(son, probe+1, restlen-1, height+1);
76
77            pt_assert(implicated(found, strlen(probe) == (size_t)restlen));
78
79            if (!found) {
80                for (int i=probe[0]+1; !found && i<PT_B_MAX; ++i) {
81                    son = PT_read_son(node, PT_BASES(i));
82                    if (son) {
83                        probe[0] = PT_BASES(i); // change probe
84                        found = findLeftmostProbe(son, probe+1, restlen-1, height+1);
85                        pt_assert(implicated(found, strlen(probe) == (size_t)restlen));
86                    }
87                }
88            }
89            return found;
90        }
91        case PT_NT_CHAIN:
92        case PT_NT_LEAF: {
93            // species list or single species reached
94            return false;
95        }
96        default: pt_assert(0); break;  // oops
97    }
98
99    pt_assert(0);
100    return false;
101}
102
103int PT_find_exProb(PT_exProb *pep, int) {
104    POS_TREE      *pt      = psg.pt; // start search at root
105    GBS_strstruct *gbs_str = GBS_stropen(pep->numget*(pep->plength+1)+1);
106    bool           first   = true;
107
108    for (int c=0; c<pep->numget; ++c) {
109        bool found = false;
110
111        if (pep->restart) {
112            pep->restart = 0;
113
114            char *probe = (char*)malloc(pep->plength+1);
115            memset(probe, 'N', pep->plength);
116            probe[pep->plength] = 0; // EOS marker
117
118            compress_data(probe);
119
120            pep->next_probe.data = probe;
121            pep->next_probe.size = pep->plength+1;
122
123            found = findLeftmostProbe(pt, pep->next_probe.data, pep->plength, 0);
124
125            pt_assert(pep->next_probe.data[pep->plength] == 0);
126            pt_assert(strlen(pep->next_probe.data) == (size_t)pep->plength);
127        }
128
129        if (!found) {
130            found = findNextProbe(pt, pep->next_probe.data, pep->plength, 0);
131
132            pt_assert(pep->next_probe.data[pep->plength] == 0);
133            pt_assert(strlen(pep->next_probe.data) == (size_t)pep->plength);
134        }
135        if (!found) break;
136
137        // append the probe to the probe list
138
139        if (!first) GBS_strcat(gbs_str, ";");
140        first = false;
141        GBS_strcat(gbs_str, pep->next_probe.data);
142    }
143
144    pep->result = GBS_strclose(gbs_str);
145
146    return 0;
147}
148
Note: See TracBrowser for help on using the repository browser.