source: tags/ms_r16q2/SL/FAST_ALIGNER/seq_search.cxx

Last change on this file was 9225, checked in by westram, 11 years ago
  • rename unit test macros
    • TEST_ASSERT.. → TEST_EXPECT..
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.1 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : seq_search.cxx                                    //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Coded by Ralf Westram (coder@reallysoft.de)                   //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "seq_search.hxx"
13
14#include <arbdbt.h>
15#include <cstdarg>
16#include <climits>
17#include <cctype>
18
19#define MESSAGE_BUFFERSIZE 300
20
21#if defined(DEBUG)
22inline int check_equal(int o, int n) {
23    if (o != n) {
24        fprintf(stderr, "o=%i\nn=%i\n", o, n);
25        fflush(stderr);
26    }
27    fa_assert(o == n);
28    return n;
29}
30#endif
31
32__ATTR__FORMAT(1) static void messagef(const char *format, ...) {
33    va_list argp;
34    char buffer[MESSAGE_BUFFERSIZE];
35
36    va_start(argp, format);
37
38    IF_ASSERTION_USED(int chars =) vsprintf(buffer, format, argp);
39    fa_assert(chars<MESSAGE_BUFFERSIZE);
40
41    va_end(argp);
42    GB_warning(buffer);
43}
44
45void Dots::append(Dots *neu) {
46    if (Next) Next->append(neu);
47    else Next = neu;
48}
49
50static CharPredicate pred_is_ali_gap(is_ali_gap);
51
52CompactedSequence::CompactedSequence(const char *Text, int Length, const char *name, int start_offset)
53    : basepos(Text, Length, pred_is_ali_gap),
54      myName(strdup(name)),
55      myStartOffset(0),  // corrected at end of ctor (otherwise calculations below go wrong)
56      dots(NULL)
57      // referred(1)
58{
59    fa_assert(Length>0);
60
61    {
62        long firstBase = basepos.first_base_abspos();
63        long lastBase  = basepos.last_base_abspos();
64        int  dotOffset = firstBase + strcspn(Text+firstBase, "."); // skip non-dots
65
66        while (dotOffset <= lastBase) {
67            fa_assert(Text[dotOffset] == '.');
68            storeDots(basepos.abs_2_rel(dotOffset));
69            dotOffset += strspn(Text+dotOffset, ".");  // skip dots
70            dotOffset += strcspn(Text+dotOffset, "."); // skip non-dots
71        }
72    }
73
74    int cLen     = length();
75    myText       = new char[cLen+1];
76    myText[cLen] = 0;
77
78    gapsBeforePosition = new int[cLen+1];
79
80    for (int cPos = 0; cPos<cLen; ++cPos) {
81        int xPos                 = basepos.rel_2_abs(cPos);
82        myText[cPos]             = toupper(Text[xPos]);
83        gapsBeforePosition[cPos] = no_of_gaps_before(cPos);
84    }
85    if (cLen>0) {
86        gapsBeforePosition[cLen] = no_of_gaps_before(cLen);    // gaps before end of sequence
87    }
88
89    myStartOffset += start_offset;
90}
91
92CompactedSequence::~CompactedSequence() {
93    delete [] myText;
94    delete[] gapsBeforePosition;
95    free(myName);
96    delete dots;
97}
98
99FastAlignInsertion::~FastAlignInsertion() {
100    if (myNext) delete myNext;
101}
102
103FastSearchSequence::FastSearchSequence(const CompactedSubSequence& seq) {
104    memset((char*)myOffset, 0, MAX_TRIPLES*sizeof(*myOffset));
105    SequencePosition triple(seq);
106
107    mySequence = &seq;
108
109    while (triple.rightOf()>=3) { // enough text for triple?
110        int tidx = triple_index(triple.text());
111        TripleOffset *top = new TripleOffset(triple.leftOf(), myOffset[tidx]);
112
113        myOffset[tidx] = top;
114        ++triple;
115    }
116}
117
118FastSearchSequence::~FastSearchSequence() {
119    for (int tidx = 0; tidx<MAX_TRIPLES; ++tidx) {
120        delete myOffset[tidx];
121    }
122}
123
124void AlignBuffer::correctUnalignedPositions() {
125    long off  = 0;
126    long rest = used;
127
128    while (rest) {
129        const char *qual  = quality();
130        char       *found = (char*)memchr(qual+off, '?', rest); // search for next '?' in myQuality
131       
132        if (!found || (found-qual) >= rest) break;
133
134        long cnt;
135        for (cnt=0; found[cnt]=='?'; cnt++) found[cnt] = '+';   // count # of unaligned positions and change them to '+'
136
137        long from  = found-myQuality;
138        long b_off = from-1;   // position before unaligned positions
139        long a_off = from+cnt; // positions after unaligned positions
140
141        long before, after;
142        for (before=0; b_off>=0 && isGlobalGap(b_off); before++, b_off--) ;             // count free positions before unaligned positions
143        for (after=0; a_off<used && isGlobalGap(a_off); after++, a_off++) ;             // count free positions after unaligned positions
144
145        if (b_off<=0) before = LONG_MAX;
146        if (a_off>=used) after = LONG_MAX;
147
148        if (before<after) { // move to the left
149            if (before) {
150                long to = from-before;
151                while (cnt--) moveUnaligned(from++, to++);
152            }
153        }
154        else if (after!=LONG_MAX) { // move to the right
155            if (after) {
156                from += cnt-1;          // copy from right to left!
157                long to = from+after;
158
159                while (cnt--) moveUnaligned(from--, to--);
160            }
161        }
162
163        rest -= (a_off-off);
164        off = a_off;
165    }
166}
167
168void AlignBuffer::restoreDots(CompactedSubSequence& slaveSequence) {
169    long rest = used;
170    long off = 0;               // real position in sequence
171    long count = 0;             // # of bases left of this position
172    long nextDot = slaveSequence.firstDotPosition();
173
174    while (nextDot!=-1 && rest) {
175        unsigned char gap = myBuffer[off];
176        if (gap=='-' || gap=='.') {
177            if (nextDot==count) {
178                myBuffer[off] = '.';
179            }
180        }
181        else {
182            count++;
183            if (nextDot==count) {
184                messagef("Couldn't restore dots at offset %li of '%s' (gap removed by aligner)",
185                         off, slaveSequence.name());
186            }
187            if (nextDot<count) {
188                nextDot = slaveSequence.nextDotPosition();
189            }
190        }
191
192        off++;
193        rest--;
194    }
195}
196
197void AlignBuffer::setDotsAtEOSequence() {
198    int i; 
199    for (i=0; (myBuffer[i]=='.' || myBuffer[i]=='-') && i<used; i++) {
200        myBuffer[i] = '.';
201    }
202
203    for (i=used-1; (myBuffer[i]=='.' || myBuffer[i]=='-') && i>=0; i--) {
204        myBuffer[i] = '.';
205    }
206}
207
208// --------------------------------------------------------------------------------
209
210#ifdef UNIT_TESTS
211#include <test_unit.h>
212#include <test_helpers.h>
213
214static int get_dotpos(const CompactedSubSequence& css, int i) {
215    int pos = css.firstDotPosition();
216    while (i) {
217        --i;
218        pos = css.nextDotPosition();
219    }
220    return pos;
221}
222static int count_dotpos(const CompactedSubSequence& css) {
223    int pos   = css.firstDotPosition();
224    int count = 0;
225    while (pos != -1) {
226        ++count;
227        pos = css.nextDotPosition();
228    }
229    return count;
230}
231
232struct bind_css {
233    CompactedSubSequence& css;
234    int test_mode;
235    bind_css(CompactedSubSequence& css_) : css(css_), test_mode(-1) {}
236
237    int operator()(int i) const {
238        switch (test_mode) {
239            case 0: return css.compPosition(i);
240            case 1: return css.expdPosition(i);
241            case 2: return css[i];
242            case 3: return css.no_of_gaps_before(i);
243            case 4: return css.no_of_gaps_after(i);
244            case 5: return get_dotpos(css, i);
245        }
246        fa_assert(0);
247        return -666;
248    }
249};
250
251#define TEST_EXPECT_CSS_SELF_REFLEXIVE(css) do {        \
252        for (int b = 0; b<css.length(); ++b) {          \
253            int x = css.expdPosition(b);                \
254            int c = css.compPosition(x);                \
255            TEST_EXPECT_EQUAL(c,b);                     \
256        }                                               \
257    } while(0)
258
259#define CSS_COMMON(in,offset)                                           \
260    int len  = strlen(in);                                              \
261    fprintf(stderr, "in='%s'\n", in);                                   \
262    CompactedSubSequence css(in, len, "noname", offset);                \
263    TEST_EXPECT_CSS_SELF_REFLEXIVE(css);                                \
264    bind_css bound_css(css)
265
266#define GEN_COMP_EXPD()                                                 \
267    bound_css.test_mode = 0;                                            \
268    char *comp = collectIntFunResults(bound_css, 0, css.expdLength()-1, 3, 0, 1); \
269    bound_css.test_mode = 1;                                            \
270    char *expd = collectIntFunResults(bound_css, 0, css.length(), 3, 0, 0)
271
272#define GEN_TEXT(in)                                                    \
273    bound_css.test_mode = 2;                                            \
274    char *text = collectIntFunResults(bound_css, 0, css.length()-1, 3, 0, 0)
275   
276#define GEN_GAPS()                                                      \
277    bound_css.test_mode = 3;                                            \
278    char *gaps_before = collectIntFunResults(bound_css, 0, css.length(), 3, 0, 0); \
279    bound_css.test_mode = 4;                                            \
280    char *gaps_after = collectIntFunResults(bound_css, 0, css.length()-1, 3, 0, 1)
281
282#define GEN_DOTS(in)                                                    \
283    bound_css.test_mode = 5;                                            \
284    char *dots = collectIntFunResults(bound_css, 0, count_dotpos(css)-1, 3, 0, 1)
285   
286#define FREE_COMP_EXPD()                                                \
287    free(expd);                                                         \
288    free(comp)
289   
290#define FREE_GAPS()                                                     \
291    free(gaps_before);                                                  \
292    free(gaps_after)
293   
294#define COMP_EXPD_CHECK(exp_comp,exp_expd)                              \
295    GEN_COMP_EXPD();                                                    \
296    TEST_EXPECT_EQUAL(comp, exp_comp);                                  \
297    TEST_EXPECT_EQUAL(expd, exp_expd);                                  \
298    FREE_COMP_EXPD()
299
300#define GAPS_CHECK(exp_before,exp_after)                                \
301    GEN_GAPS();                                                         \
302    TEST_EXPECT_EQUAL(gaps_before, exp_before);                         \
303    TEST_EXPECT_EQUAL(gaps_after, exp_after);                           \
304    FREE_GAPS()
305
306#define DOTS_CHECK(exp_dots)                                            \
307    GEN_DOTS();                                                         \
308    TEST_EXPECT_EQUAL(dots, exp_dots);                                  \
309    free(dots)
310
311// ------------------------------------------------------------
312       
313#define TEST_CS_EQUALS(in,exp_comp,exp_expd) do {                       \
314        CSS_COMMON(in, 0);                                              \
315        COMP_EXPD_CHECK(exp_comp,exp_expd);                             \
316    } while(0)
317
318#define TEST_CS_EQUALS_OFFSET(in,offset,exp_comp,exp_expd) do {         \
319        CSS_COMMON(in, offset);                                         \
320        COMP_EXPD_CHECK(exp_comp,exp_expd);                             \
321    } while(0)
322       
323#define TEST_GAPS_EQUALS_OFFSET(in,offset,exp_before,exp_after) do {    \
324        CSS_COMMON(in, offset);                                         \
325        GAPS_CHECK(exp_before,exp_after);                               \
326    } while(0)
327
328#define TEST_DOTS_EQUALS_OFFSET(in,offset,exp_dots) do {                \
329        CSS_COMMON(in, offset);                                         \
330        DOTS_CHECK(exp_dots);                                           \
331    } while(0)
332
333#define TEST_CS_TEXT(in,exp_text) do {                                  \
334        CSS_COMMON(in, 0);                                              \
335        GEN_TEXT(in);                                                   \
336        TEST_EXPECT_EQUAL(text, exp_text);                              \
337        free(text);                                                     \
338    } while(0)
339
340#define TEST_CS_CBROKN(in,exp_comp,exp_expd) do {                       \
341        CSS_COMMON(in, 0);                                              \
342        GEN_COMP_EXPD();                                                \
343        TEST_EXPECT_EQUAL__BROKEN(comp, exp_comp);                      \
344        TEST_EXPECT_EQUAL(expd, exp_expd);                              \
345        FREE_COMP_EXPD();                                               \
346    } while(0)
347
348void TEST_EARLY_CompactedSequence() {
349    // reproduce a bug in compPosition
350
351    // no base
352    TEST_CS_EQUALS("-",          "  0  [0]",       "  1");
353    TEST_CS_EQUALS("--",         "  0  0  [0]",    "  2");
354    TEST_CS_EQUALS("---",        "  0  0  0  [0]", "  3");
355
356    TEST_CS_TEXT("---", "");
357    TEST_CS_TEXT("-?~", "");
358    TEST_CS_TEXT("-A-", " 65");    // "A"
359    TEST_CS_TEXT("A-C", " 65 67"); // "AC"
360
361    TEST_CS_EQUALS("----------", "  0  0  0  0  0  0  0  0  0  0  [0]", " 10");
362
363    // one base
364    TEST_CS_EQUALS("A---------", "  0  1  1  1  1  1  1  1  1  1  [1]", "  0 10");
365    TEST_CS_EQUALS("-A--------", "  0  0  1  1  1  1  1  1  1  1  [1]", "  1 10"); 
366    TEST_CS_EQUALS("---A------", "  0  0  0  0  1  1  1  1  1  1  [1]", "  3 10"); 
367    TEST_CS_EQUALS("-----A----", "  0  0  0  0  0  0  1  1  1  1  [1]", "  5 10"); 
368    TEST_CS_EQUALS("-------A--", "  0  0  0  0  0  0  0  0  1  1  [1]", "  7 10"); 
369    TEST_CS_EQUALS("---------A", "  0  0  0  0  0  0  0  0  0  0  [1]", "  9 10"); 
370
371    // two bases
372    TEST_CS_EQUALS("AC--------", "  0  1  2  2  2  2  2  2  2  2  [2]", "  0  1 10");
373
374    TEST_CS_EQUALS("A-C-------", "  0  1  1  2  2  2  2  2  2  2  [2]", "  0  2 10");
375    TEST_CS_EQUALS("A--------C", "  0  1  1  1  1  1  1  1  1  1  [2]", "  0  9 10");
376    TEST_CS_EQUALS("-AC-------", "  0  0  1  2  2  2  2  2  2  2  [2]", "  1  2 10");
377    TEST_CS_EQUALS("-A------C-", "  0  0  1  1  1  1  1  1  1  2  [2]", "  1  8 10");
378    TEST_CS_EQUALS("-A-------C", "  0  0  1  1  1  1  1  1  1  1  [2]", "  1  9 10");
379    TEST_CS_EQUALS("-------A-C", "  0  0  0  0  0  0  0  0  1  1  [2]", "  7  9 10");
380    TEST_CS_EQUALS("--------AC", "  0  0  0  0  0  0  0  0  0  1  [2]", "  8  9 10");
381
382    // 3 bases
383    TEST_CS_EQUALS("ACG-------", "  0  1  2  3  3  3  3  3  3  3  [3]", "  0  1  2 10");
384    TEST_CS_EQUALS("AC---G----", "  0  1  2  2  2  2  3  3  3  3  [3]", "  0  1  5 10"); 
385    TEST_CS_EQUALS("A-C--G----", "  0  1  1  2  2  2  3  3  3  3  [3]", "  0  2  5 10"); 
386    TEST_CS_EQUALS("A-C--G----", "  0  1  1  2  2  2  3  3  3  3  [3]", "  0  2  5 10"); 
387    TEST_CS_EQUALS("A--C-G----", "  0  1  1  1  2  2  3  3  3  3  [3]", "  0  3  5 10"); 
388    TEST_CS_EQUALS("A---CG----", "  0  1  1  1  1  2  3  3  3  3  [3]", "  0  4  5 10");
389
390    TEST_CS_EQUALS("A---C---G-", "  0  1  1  1  1  2  2  2  2  3  [3]", "  0  4  8 10");
391    TEST_CS_EQUALS("A---C----G", "  0  1  1  1  1  2  2  2  2  2  [3]", "  0  4  9 10");
392    TEST_CS_EQUALS("A~~~C????G", "  0  1  1  1  1  2  2  2  2  2  [3]", "  0  4  9 10");
393
394    // 4 bases
395    TEST_CS_EQUALS("-AC-G--T--", "  0  0  1  2  2  3  3  3  4  4  [4]", "  1  2  4  7 10"); 
396
397    // 9 bases
398    TEST_CS_EQUALS("-CGTACGTAC", "  0  0  1  2  3  4  5  6  7  8  [9]", "  1  2  3  4  5  6  7  8  9 10");
399    TEST_CS_EQUALS("A-GTACGTAC", "  0  1  1  2  3  4  5  6  7  8  [9]", "  0  2  3  4  5  6  7  8  9 10");
400    TEST_CS_EQUALS("ACGTA-GTAC", "  0  1  2  3  4  5  5  6  7  8  [9]", "  0  1  2  3  4  6  7  8  9 10");
401    TEST_CS_EQUALS("ACGTACGT-C", "  0  1  2  3  4  5  6  7  8  8  [9]", "  0  1  2  3  4  5  6  7  9 10");
402    TEST_CS_EQUALS("ACGTACGTA-", "  0  1  2  3  4  5  6  7  8  9  [9]", "  0  1  2  3  4  5  6  7  8 10");
403
404    // all 10 bases
405    TEST_CS_EQUALS("ACGTACGTAC", "  0  1  2  3  4  5  6  7  8  9 [10]", "  0  1  2  3  4  5  6  7  8  9 10");
406
407    TEST_CS_EQUALS_OFFSET("A--C-G-", 0, "  0  1  1  1  2  2  3  [3]",             "  0  3  5  7");
408    TEST_CS_EQUALS_OFFSET("A--C-G-", 2, "  0  0  0  1  1  1  2  2  3  [3]",       "  2  5  7  9");
409    TEST_CS_EQUALS_OFFSET("A--C-G-", 3, "  0  0  0  0  1  1  1  2  2  3  [3]",    "  3  6  8 10");
410    TEST_CS_EQUALS_OFFSET("A--C-G-", 4, "  0  0  0  0  0  1  1  1  2  2  3  [3]", "  4  7  9 11");
411
412    // test no_of_gaps_before() and no_of_gaps_after()
413    TEST_GAPS_EQUALS_OFFSET("-AC---G",     0, "  1  0  3  0", "  0  3  0 [-1]");
414    TEST_GAPS_EQUALS_OFFSET(".AC..-G",     0, "  1  0  3  0", "  0  3  0 [-1]");
415    TEST_GAPS_EQUALS_OFFSET("~AC?\?-G",     0, "  1  0  3  0", "  0  3  0 [-1]");
416    TEST_GAPS_EQUALS_OFFSET("A--C-G-",     0, "  0  2  1  1", "  2  1  1 [-1]");
417    TEST_GAPS_EQUALS_OFFSET("A--C-G-",  1000, "  0  2  1  1", "  2  1  1 [-1]"); // is independent from offset
418
419    // test dots
420    TEST_DOTS_EQUALS_OFFSET("----ACG--T--A-C--GT----", 0, " [-1]");    // no dots
421    TEST_DOTS_EQUALS_OFFSET("....ACG--T--A-C--GT....", 0, " [-1]");    // no extraordinary dots
422    TEST_DOTS_EQUALS_OFFSET("....ACG--T..A-C--GT....", 0, "  4 [-1]"); // i.e. "before base number 4"
423    TEST_DOTS_EQUALS_OFFSET("....ACG..T--A.C--GT....", 0, "  3  5 [-1]");
424    TEST_DOTS_EQUALS_OFFSET("....ACG..T~~A.C??GT....", 0, "  3  5 [-1]");
425
426    TEST_DOTS_EQUALS_OFFSET("AC-----GTA-----CGT", 0, " [-1]");
427    // just care THAT dots occur in a gap, dont care how many
428    TEST_DOTS_EQUALS_OFFSET("A-.-G-...-C...T", 0, "  1  2  3 [-1]");
429    TEST_DOTS_EQUALS_OFFSET("A.--G...--C...T", 0, "  1  2  3 [-1]");
430    TEST_DOTS_EQUALS_OFFSET("A--.G--...C...T", 0, "  1  2  3 [-1]");
431}
432
433#endif // UNIT_TESTS
434
Note: See TracBrowser for help on using the repository browser.