source: branches/lib/GDE/MUSCLE/src/estring.cpp

Last change on this file was 10390, checked in by aboeckma, 12 years ago

added muscle sourcles amd makefile
File size: 13.1 KB
Line 
1#include "muscle.h"
2#include "pwpath.h"
3#include "estring.h"
4#include "seq.h"
5#include "msa.h"
6
7/***
8An "estring" is an edit string that operates on a sequence.
9An estring is represented as a vector of integers.
10It is interpreted in order of increasing suffix.
11A positive value n means copy n letters.
12A negative value -n means insert n indels.
13Zero marks the end of the vector.
14Consecutive entries must have opposite sign, i.e. the
15shortest possible representation must be used.
16
17A "tpair" is a traceback path for a pairwise alignment
18represented as two estrings, one for each sequence.
19***/
20
21#define c2(c,d) (((unsigned char) c) << 8 | (unsigned char) d)
22
23unsigned LengthEstring(const short es[])
24        {
25        unsigned i = 0;
26        while (*es++ != 0)
27                ++i;
28        return i;
29        }
30
31short *EstringNewCopy(const short es[])
32        {
33        unsigned n = LengthEstring(es) + 1;
34        short *esNew = new short[n];
35        memcpy(esNew, es, n*sizeof(short));
36        return esNew;
37        }
38
39void LogEstring(const short es[])
40        {
41        Log("<");
42        for (unsigned i = 0; es[i] != 0; ++i)
43                {
44                if (i > 0)
45                        Log(" ");
46                Log("%d", es[i]);
47                }
48        Log(">");
49        }
50
51static bool EstringsEq(const short es1[], const short es2[])
52        {
53        for (;;)
54                {
55                if (*es1 != *es2)
56                        return false;
57                if (0 == *es1)
58                        break;
59                ++es1;
60                ++es2;
61                }
62        return true;
63        }
64
65static void EstringCounts(const short es[], unsigned *ptruSymbols,
66  unsigned *ptruIndels)
67        {
68        unsigned uSymbols = 0;
69        unsigned uIndels = 0;
70        for (unsigned i = 0; es[i] != 0; ++i)
71                {
72                short n = es[i];
73                if (n > 0)
74                        uSymbols += n;
75                else if (n < 0)
76                        uIndels += -n;
77                }
78        *ptruSymbols = uSymbols;
79        *ptruIndels = uIndels;
80        }
81
82static char *EstringOp(const short es[], const char s[])
83        {
84        unsigned uSymbols;
85        unsigned uIndels;
86        EstringCounts(es, &uSymbols, &uIndels);
87        assert((unsigned) strlen(s) == uSymbols);
88        char *sout = new char[uSymbols + uIndels + 1];
89        char *psout = sout;
90        for (;;)
91                {
92                int n = *es++;
93                if (0 == n)
94                        break;
95                if (n > 0)
96                        for (int i = 0; i < n; ++i)
97                                *psout++ = *s++;
98                else
99                        for (int i = 0; i < -n; ++i)
100                                *psout++ = '-';
101                }
102        assert(0 == *s);
103        *psout = 0;
104        return sout;
105        }
106
107void EstringOp(const short es[], const Seq &sIn, Seq &sOut)
108        {
109#if     DEBUG
110        unsigned uSymbols;
111        unsigned uIndels;
112        EstringCounts(es, &uSymbols, &uIndels);
113        assert(sIn.Length() == uSymbols);
114#endif
115        sOut.Clear();
116        sOut.SetName(sIn.GetName());
117        int p = 0;
118        for (;;)
119                {
120                int n = *es++;
121                if (0 == n)
122                        break;
123                if (n > 0)
124                        for (int i = 0; i < n; ++i)
125                                {
126                                char c = sIn[p++];
127                                sOut.push_back(c);
128                                }
129                else
130                        for (int i = 0; i < -n; ++i)
131                                sOut.push_back('-');
132                }
133        }
134
135unsigned EstringOp(const short es[], const Seq &sIn, MSA &a)
136        {
137        unsigned uSymbols;
138        unsigned uIndels;
139        EstringCounts(es, &uSymbols, &uIndels);
140        assert(sIn.Length() == uSymbols);
141
142        unsigned uColCount = uSymbols + uIndels;
143
144        a.Clear();
145        a.SetSize(1, uColCount);
146
147        a.SetSeqName(0, sIn.GetName());
148        a.SetSeqId(0, sIn.GetId());
149
150        unsigned p = 0;
151        unsigned uColIndex = 0;
152        for (;;)
153                {
154                int n = *es++;
155                if (0 == n)
156                        break;
157                if (n > 0)
158                        for (int i = 0; i < n; ++i)
159                                {
160                                char c = sIn[p++];
161                                a.SetChar(0, uColIndex++, c);
162                                }
163                else
164                        for (int i = 0; i < -n; ++i)
165                                a.SetChar(0, uColIndex++, '-');
166                }
167        assert(uColIndex == uColCount);
168        return uColCount;
169        }
170
171void PathToEstrings(const PWPath &Path, short **ptresA, short **ptresB)
172        {
173// First pass to determine size of estrings esA and esB
174        const unsigned uEdgeCount = Path.GetEdgeCount();
175        if (0 == uEdgeCount)
176                {
177                short *esA = new short[1];
178                short *esB = new short[1];
179                esA[0] = 0;
180                esB[0] = 0;
181                *ptresA = esA;
182                *ptresB = esB;
183                return;
184                }
185
186        unsigned iLengthA = 1;
187        unsigned iLengthB = 1;
188        const char cFirstEdgeType = Path.GetEdge(0).cType;
189        char cPrevEdgeType = cFirstEdgeType;
190        for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
191                {
192                const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
193                char cEdgeType = Edge.cType;
194
195                switch (c2(cPrevEdgeType, cEdgeType))
196                        {
197                case c2('M', 'M'):
198                case c2('D', 'D'):
199                case c2('I', 'I'):
200                        break;
201
202                case c2('D', 'M'):
203                case c2('M', 'D'):
204                        ++iLengthB;
205                        break;
206
207                case c2('I', 'M'):
208                case c2('M', 'I'):
209                        ++iLengthA;
210                        break;
211
212                case c2('I', 'D'):
213                case c2('D', 'I'):
214                        ++iLengthB;
215                        ++iLengthA;
216                        break;
217
218                default:
219                        assert(false);
220                        }
221                cPrevEdgeType = cEdgeType;
222                }
223
224// Pass2 for seq A
225        {
226        short *esA = new short[iLengthA+1];
227        unsigned iA = 0;
228        switch (Path.GetEdge(0).cType)
229                {
230        case 'M':
231        case 'D':
232                esA[0] = 1;
233                break;
234
235        case 'I':
236                esA[0] = -1;
237                break;
238
239        default:
240                assert(false);
241                }
242
243        char cPrevEdgeType = cFirstEdgeType;
244        for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
245                {
246                const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
247                char cEdgeType = Edge.cType;
248
249                switch (c2(cPrevEdgeType, cEdgeType))
250                        {
251                case c2('M', 'M'):
252                case c2('D', 'D'):
253                case c2('D', 'M'):
254                case c2('M', 'D'):
255                        ++(esA[iA]);
256                        break;
257
258                case c2('I', 'D'):
259                case c2('I', 'M'):
260                        ++iA;
261                        esA[iA] = 1;
262                        break;
263
264                case c2('M', 'I'):
265                case c2('D', 'I'):
266                        ++iA;
267                        esA[iA] = -1;
268                        break;
269
270                case c2('I', 'I'):
271                        --(esA[iA]);
272                        break;
273
274                default:
275                        assert(false);
276                        }
277
278                cPrevEdgeType = cEdgeType;
279                }
280        assert(iA == iLengthA - 1);
281        esA[iLengthA] = 0;
282        *ptresA = esA;
283        }
284
285        {
286// Pass2 for seq B
287        short *esB = new short[iLengthB+1];
288        unsigned iB = 0;
289        switch (Path.GetEdge(0).cType)
290                {
291        case 'M':
292        case 'I':
293                esB[0] = 1;
294                break;
295
296        case 'D':
297                esB[0] = -1;
298                break;
299
300        default:
301                assert(false);
302                }
303
304        char cPrevEdgeType = cFirstEdgeType;
305        for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
306                {
307                const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
308                char cEdgeType = Edge.cType;
309
310                switch (c2(cPrevEdgeType, cEdgeType))
311                        {
312                case c2('M', 'M'):
313                case c2('I', 'I'):
314                case c2('I', 'M'):
315                case c2('M', 'I'):
316                        ++(esB[iB]);
317                        break;
318
319                case c2('D', 'I'):
320                case c2('D', 'M'):
321                        ++iB;
322                        esB[iB] = 1;
323                        break;
324
325                case c2('M', 'D'):
326                case c2('I', 'D'):
327                        ++iB;
328                        esB[iB] = -1;
329                        break;
330
331                case c2('D', 'D'):
332                        --(esB[iB]);
333                        break;
334
335                default:
336                        assert(false);
337                        }
338
339                cPrevEdgeType = cEdgeType;
340                }
341        assert(iB == iLengthB - 1);
342        esB[iLengthB] = 0;
343        *ptresB = esB;
344        }
345
346#if     DEBUG
347        {
348        const PWEdge &LastEdge = Path.GetEdge(uEdgeCount - 1);
349        unsigned uSymbols;
350        unsigned uIndels;
351        EstringCounts(*ptresA, &uSymbols, &uIndels);
352        assert(uSymbols == LastEdge.uPrefixLengthA);
353        assert(uSymbols + uIndels == uEdgeCount);
354
355        EstringCounts(*ptresB, &uSymbols, &uIndels);
356        assert(uSymbols == LastEdge.uPrefixLengthB);
357        assert(uSymbols + uIndels == uEdgeCount);
358
359        PWPath TmpPath;
360        EstringsToPath(*ptresA, *ptresB, TmpPath);
361        TmpPath.AssertEqual(Path);
362        }
363#endif
364        }
365
366void EstringsToPath(const short esA[], const short esB[], PWPath &Path)
367        {
368        Path.Clear();
369        unsigned iA = 0;
370        unsigned iB = 0;
371        int nA = esA[iA++];
372        int nB = esB[iB++];
373        unsigned uPrefixLengthA = 0;
374        unsigned uPrefixLengthB = 0;
375        for (;;)
376                {
377                char cType;
378                if (nA > 0)
379                        {
380                        if (nB > 0)
381                                {
382                                cType = 'M';
383                                --nA;
384                                --nB;
385                                }
386                        else if (nB < 0)
387                                {
388                                cType = 'D';
389                                --nA;
390                                ++nB;
391                                }
392                        else
393                                assert(false);
394                        }
395                else if (nA < 0)
396                        {
397                        if (nB > 0)
398                                {
399                                cType = 'I';
400                                ++nA;
401                                --nB;
402                                }
403                        else
404                                assert(false);
405                        }
406                else
407                        assert(false);
408
409                switch (cType)
410                        {
411                case 'M':
412                        ++uPrefixLengthA;
413                        ++uPrefixLengthB;
414                        break;
415                case 'D':
416                        ++uPrefixLengthA;
417                        break;
418                case 'I':
419                        ++uPrefixLengthB;
420                        break;
421                        }
422
423                PWEdge Edge;
424                Edge.cType = cType;
425                Edge.uPrefixLengthA = uPrefixLengthA;
426                Edge.uPrefixLengthB = uPrefixLengthB;
427                Path.AppendEdge(Edge);
428
429                if (nA == 0)
430                        {
431                        if (0 == esA[iA])
432                                {
433                                assert(0 == esB[iB]);
434                                break;
435                                }
436                        nA = esA[iA++];
437                        }
438                if (nB == 0)
439                        nB = esB[iB++];
440                }
441        }
442
443/***
444Multiply two estrings to make a third estring.
445The product of two estrings e1*e2 is defined to be
446the estring that produces the same result as applying
447e1 then e2. Multiplication is not commutative. In fact,
448the reversed order is undefined unless both estrings
449consist of a single, identical, positive entry.
450A primary motivation for using estrings is that
451multiplication is very fast, reducing the time
452needed to construct the root alignment.
453
454Example
455
456        <-1,3>(XXX)     = -XXX
457        <2,-1,2>(-XXX) = -X-XX
458
459Therefore,
460
461        <-1,3>*<2,-1,2> = <-1,1,-1,2>
462***/
463
464static bool CanMultiplyEstrings(const short es1[], const short es2[])
465        {
466        unsigned uSymbols1;
467        unsigned uSymbols2;
468        unsigned uIndels1;
469        unsigned uIndels2;
470        EstringCounts(es1, &uSymbols1, &uIndels1);
471        EstringCounts(es2, &uSymbols2, &uIndels2);
472        return uSymbols1 + uIndels1 == uSymbols2;
473        }
474
475static inline void AppendGaps(short esp[], int &ip, int n)
476        {
477        if (-1 == ip)
478                esp[++ip] = n;
479        else if (esp[ip] < 0)
480                esp[ip] += n;
481        else
482                esp[++ip] = n;
483        }
484
485static inline void AppendSymbols(short esp[], int &ip, int n)
486        {
487        if (-1 == ip)
488                esp[++ip] = n;
489        else if (esp[ip] > 0)
490                esp[ip] += n;
491        else
492                esp[++ip] = n;
493        }
494
495void MulEstrings(const short es1[], const short es2[], short esp[])
496        {
497        assert(CanMultiplyEstrings(es1, es2));
498
499        unsigned i1 = 0;
500        int ip = -1;
501        int n1 = es1[i1++];
502        for (unsigned i2 = 0; ; ++i2)
503                {
504                int n2 = es2[i2];
505                if (0 == n2)
506                        break;
507                if (n2 > 0)
508                        {
509                        for (;;)
510                                {
511                                if (n1 < 0)
512                                        {
513                                        if (n2 > -n1)
514                                                {
515                                                AppendGaps(esp, ip, n1);
516                                                n2 += n1;
517                                                n1 = es1[i1++];
518                                                }
519                                        else if (n2 == -n1)
520                                                {
521                                                AppendGaps(esp, ip, n1);
522                                                n1 = es1[i1++];
523                                                break;
524                                                }
525                                        else
526                                                {
527                                                assert(n2 < -n1);
528                                                AppendGaps(esp, ip, -n2);
529                                                n1 += n2;
530                                                break;
531                                                }
532                                        }
533                                else
534                                        {
535                                        assert(n1 > 0);
536                                        if (n2 > n1)
537                                                {
538                                                AppendSymbols(esp, ip, n1);
539                                                n2 -= n1;
540                                                n1 = es1[i1++];
541                                                }
542                                        else if (n2 == n1)
543                                                {
544                                                AppendSymbols(esp, ip, n1);
545                                                n1 = es1[i1++];
546                                                break;
547                                                }
548                                        else
549                                                {
550                                                assert(n2 < n1);
551                                                AppendSymbols(esp, ip, n2);
552                                                n1 -= n2;
553                                                break;
554                                                }
555                                        }
556                                }
557                        }
558                else
559                        {
560                        assert(n2 < 0);
561                        AppendGaps(esp, ip, n2);
562                        }
563                }
564        esp[++ip] = 0;
565
566#if     DEBUG
567        {
568        int MaxLen = (int) (LengthEstring(es1) + LengthEstring(es2) + 1);
569        assert(ip < MaxLen);
570        if (ip >= 2)
571                for (int i = 0; i < ip - 2; ++i)
572                        {
573                        if (!(esp[i] > 0 && esp[i+1] < 0 || esp[i] < 0 && esp[i+1] > 0))
574                                {
575                                Log("Bad result of MulEstring: ");
576                                LogEstring(esp);
577                                Quit("Assert failed (alternating signs)");
578                                }
579                        }
580        unsigned uSymbols1;
581        unsigned uSymbols2;
582        unsigned uSymbolsp;
583        unsigned uIndels1;
584        unsigned uIndels2;
585        unsigned uIndelsp;
586        EstringCounts(es1, &uSymbols1, &uIndels1);
587        EstringCounts(es2, &uSymbols2, &uIndels2);
588        EstringCounts(esp, &uSymbolsp, &uIndelsp);
589        if (uSymbols1 + uIndels1 != uSymbols2)
590                {
591                Log("Bad result of MulEstring: ");
592                LogEstring(esp);
593                Quit("Assert failed (counts1 %u %u %u)",
594                  uSymbols1, uIndels1, uSymbols2);
595                }
596        }
597#endif
598        }
599
600static void test(const short es1[], const short es2[], const short esa[])
601        {
602        unsigned uSymbols1;
603        unsigned uSymbols2;
604        unsigned uIndels1;
605        unsigned uIndels2;
606        EstringCounts(es1, &uSymbols1, &uIndels1);
607        EstringCounts(es2, &uSymbols2, &uIndels2);
608
609        char s[4096];
610        memset(s, 'X', sizeof(s));
611        s[uSymbols1] = 0;
612
613        char *s1 = EstringOp(es1, s);
614        char *s12 = EstringOp(es2, s1);
615
616        memset(s, 'X', sizeof(s));
617        s[uSymbols2] = 0;
618        char *s2 = EstringOp(es2, s);
619
620        Log("%s * %s = %s\n", s1, s2, s12);
621
622        LogEstring(es1);
623        Log(" * ");
624        LogEstring(es2);
625        Log(" = ");
626        LogEstring(esa);
627        Log("\n");
628
629        short esp[4096];
630        MulEstrings(es1, es2, esp);
631        LogEstring(esp);
632        if (!EstringsEq(esp, esa))
633                Log(" *ERROR* ");
634        Log("\n");
635
636        memset(s, 'X', sizeof(s));
637        s[uSymbols1] = 0;
638        char *sp = EstringOp(esp, s);
639        Log("%s\n", sp);
640        Log("\n==========\n\n");
641        }
642
643void TestEstrings()
644        {
645        SetListFileName("c:\\tmp\\muscle.log", false);
646        //{
647        //short es1[] = { -1, 1, -1, 0 };
648        //short es2[] = { 1, -1, 2, 0 };
649        //short esa[] = { -2, 1, -1, 0 };
650        //test(es1, es2, esa);
651        //}
652        //{
653        //short es1[] = { 2, -1, 2, 0 };
654        //short es2[] = { 1, -1, 3, -1, 1, 0 };
655        //short esa[] = { 1, -1, 1, -1, 1, -1, 1, 0 };
656        //test(es1, es2, esa);
657        //}
658        //{
659        //short es1[] = { -1, 3, 0 };
660        //short es2[] = { 2, -1, 2, 0 };
661        //short esa[] = { -1, 1, -1, 2, 0 };
662        //test(es1, es2, esa);
663        //}
664        //{
665        //short es1[] = { -1, 1, -1, 1, 0};
666        //short es2[] = { 4, 0 };
667        //short esa[] = { -1, 1, -1, 1, 0};
668        //test(es1, es2, esa);
669        //}
670        //{
671        //short es1[] = { 1, -1, 1, -1, 0};
672        //short es2[] = { 4, 0 };
673        //short esa[] = { 1, -1, 1, -1, 0};
674        //test(es1, es2, esa);
675        //}
676        //{
677        //short es1[] = { 1, -1, 1, -1, 0};
678        //short es2[] = { -1, 4, -1, 0 };
679        //short esa[] = { -1, 1, -1, 1, -2, 0};
680        //test(es1, es2, esa);
681        //}
682        {
683        short es1[] = { 106, -77, 56, -2, 155, -3, 123, -2, 0};
684        short es2[] = { 50, -36, 34, -3, 12, -6, 1, -6, 18, -17, 60, -5, 349, -56, 0 };
685        short esa[] = { 0 };
686        test(es1, es2, esa);
687        }
688        exit(0);
689        }
Note: See TracBrowser for help on using the repository browser.