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

Last change on this file was 10946, checked in by westram, 11 years ago
  • remove outdated cppcheck-suppressions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.2 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : ClustalV.cxx                                      //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Coded based on ClustalV by Ralf Westram (coder@reallysoft.de) //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "ClustalV.hxx"
13#include "seq_search.hxx"
14
15#include <cctype>
16#include <climits>
17
18
19// ----------------------------------------------------------------
20
21#define MASTER_GAP_OPEN                 50
22#define CHEAP_GAP_OPEN                  20      // penalty for creating a gap (if we have gaps in master)
23
24#define DEFAULT_GAP_OPEN                30      // penalty for creating a gap
25
26#define MASTER_GAP_EXTEND               18
27#define CHEAP_GAP_EXTEND                5       // penalty for extending a gap (if we have gaps in master)
28
29#define DEFAULT_GAP_EXTEND              10      // penalty for extending a gap
30
31#define DEFAULT_IMPROBABLY_MUTATION     10      // penalty for mutations
32#define DEFAULT_PROBABLY_MUTATION       4       // penalty for special mutations (A<->G,C<->U/T)        (only if 'weighted')
33
34#define DYNAMIC_PENALTIES                       // otherwise you get FIXED PENALTIES    (=no cheap penalties)
35
36// ----------------------------------------------------------------
37
38#define MAX_GAP_OPEN_DISCOUNT           (DEFAULT_GAP_OPEN-CHEAP_GAP_OPEN)       // maximum subtracted from DEFAULT_GAP_OPEN
39#define MAX_GAP_EXTEND_DISCOUNT         (DEFAULT_GAP_EXTEND-CHEAP_GAP_EXTEND)   // maximum subtracted from DEFAULT_GAP_EXTEND
40
41#define MAXN    2       // Maximum number of sequences (both groups)
42
43static bool module_initialized = false;
44
45typedef int Boolean;
46
47static Boolean dnaflag;
48static Boolean is_weight;
49
50#define MAX_BASETYPES 21
51
52static int xover;
53static int little_pam;
54static int big_pam;
55static int pamo[(MAX_BASETYPES-1)*MAX_BASETYPES/2];
56static int pam[MAX_BASETYPES][MAX_BASETYPES];
57
58static int             pos1;
59static int             pos2;
60static int           **naa1;                 // naa1[basetype][position]     counts bases for each position of all sequences in group1
61static int           **naa2;                 // naa2[basetype][position]     same for group2
62static int           **naas;                 //
63static int             seqlen_array[MAXN+1]; // length of all sequences
64static unsigned char  *seq_array[MAXN+1];    // the sequences
65static int             group[MAXN+1];        // group of sequence
66static int             alist[MAXN+1];        // indices of sequences to be aligned
67static int             fst_list[MAXN+1];
68static int             snd_list[MAXN+1];
69
70static int nseqs;               // # of sequences
71static int weights[MAX_BASETYPES][MAX_BASETYPES];     // weights[b1][b2] : penalty for mutation from base 'b1' to base 'b2'
72
73#if defined(ASSERTION_USED)
74static size_t displ_size = 0;
75#endif // ASSERTION_USED
76
77static int *displ;                                  // displ == 0 -> base in both , displ<0 -> displ gaps in slave, displ>0 -> displ gaps in master
78static int *zza;                                    // column (left->right) of align matrix (minimum of all paths to this matrix element)
79static int *zzb;                                    // -------------- " ------------------- (minimum of all paths, where gap inserted into slave)
80static int *zzc;                                    // column (left<-right) of align matrix (minimum of all paths to this matrix element)
81static int *zzd;                                    // -------------- " ------------------- (minimum of all paths, where gap inserted into slave)
82static int  print_ptr;
83static int  last_print;
84
85static const int *gaps_before_position;
86
87#if defined(DEBUG)
88// #define MATRIX_DUMP
89// #define DISPLAY_DIFF
90#endif // DEBUG
91
92#ifdef MATRIX_DUMP
93# define IF_MATRIX_DUMP(xxx) xxx
94# define DISPLAY_MATRIX_SIZE 3000
95
96static int vertical             [DISPLAY_MATRIX_SIZE+2][DISPLAY_MATRIX_SIZE+2];
97static int verticalOpen         [DISPLAY_MATRIX_SIZE+2][DISPLAY_MATRIX_SIZE+2];
98static int diagonal             [DISPLAY_MATRIX_SIZE+2][DISPLAY_MATRIX_SIZE+2];
99static int horizontal           [DISPLAY_MATRIX_SIZE+2][DISPLAY_MATRIX_SIZE+2];
100static int horizontalOpen       [DISPLAY_MATRIX_SIZE+2][DISPLAY_MATRIX_SIZE+2];
101
102#else
103# define IF_MATRIX_DUMP(xxx)
104#endif
105
106inline int master_gap_open(int beforePosition) {
107#ifdef DYNAMIC_PENALTIES
108    long gaps = gaps_before_position[beforePosition-1];
109    return gaps ? MASTER_GAP_OPEN - MAX_GAP_OPEN_DISCOUNT : MASTER_GAP_OPEN;
110
111    /*    return
112          gaps >= MAX_GAP_OPEN_DISCOUNT
113          ? DEFAULT_GAP_OPEN-MAX_GAP_OPEN_DISCOUNT
114          : DEFAULT_GAP_OPEN-gaps; */
115#else
116    return DEFAULT_GAP_OPEN;
117#endif
118}
119inline int master_gap_extend(int beforePosition) {
120#ifdef DYNAMIC_PENALTIES
121    long gaps = gaps_before_position[beforePosition-1];
122    return gaps ? MASTER_GAP_EXTEND - MAX_GAP_EXTEND_DISCOUNT : MASTER_GAP_EXTEND;
123    /*    return
124          gaps >= MAX_GAP_EXTEND_DISCOUNT
125          ? DEFAULT_GAP_EXTEND-MAX_GAP_EXTEND_DISCOUNT
126          : DEFAULT_GAP_EXTEND-gaps; */
127#else
128    return DEFAULT_GAP_EXTEND;
129#endif
130}
131
132inline int master_gapAtWithOpenPenalty(int atPosition, int length, int penalty) {
133    if (length<=0) return 0;
134
135    int beforePosition = atPosition;
136    int afterPosition  = atPosition-1;
137
138    while (length--) {
139        int p1 = master_gap_extend(beforePosition);
140        int p2 = master_gap_extend(afterPosition+1);
141
142        if (p1<p2 && beforePosition>1) {
143            penalty += p1;
144            beforePosition--;
145        }
146        else {
147            penalty += p2;
148            afterPosition++;
149        }
150    }
151
152    return penalty;
153}
154
155inline int master_gapAt(int atPosition, int length) {
156    return master_gapAtWithOpenPenalty(atPosition, length, master_gap_open(atPosition));
157}
158
159inline int slave_gap_open(int /* beforePosition */) {
160    return DEFAULT_GAP_OPEN;
161}
162
163inline int slave_gap_extend(int /* beforePosition */) {
164    return DEFAULT_GAP_EXTEND;
165}
166
167inline int slave_gapAtWithOpenPenalty(int atPosition, int length, int penalty) {
168    return length<=0 ? 0 : penalty + length*slave_gap_extend(atPosition);
169}
170
171inline int slave_gapAt(int atPosition, int length) {
172    return slave_gapAtWithOpenPenalty(atPosition, length, slave_gap_open(atPosition));
173}
174
175#define UNKNOWN_ACID 255
176static const char *amino_acid_order   = "XCSTPAGNDEQHRKMILVFYW";
177
178#define NUCLEIDS 16
179static const char *nucleic_acid_order = "-ACGTUMRWSYKVHDBN";
180static const char *nucleic_maybe_A    = "-A----AAA---AAA-A";
181static const char *nucleic_maybe_C    = "--C---C--CC-CC-CC";
182static const char *nucleic_maybe_G    = "---G---G-G-GG-GGG";
183static const char *nucleic_maybe_T    = "----T---T-TT-TTTT";
184static const char *nucleic_maybe_U    = "-----U--U-UU-UUUU";
185static const char *nucleic_maybe[6]   = { NULL, nucleic_maybe_A, nucleic_maybe_C, nucleic_maybe_G, nucleic_maybe_T, nucleic_maybe_U };
186
187/*
188 *      M = A or C              S = G or C              V = A or G or C         N = A or C or G or T
189 *      R = A or G              Y = C or T              H = A or C or T
190 *      W = A or T              K = G or T              D = A or G or T
191 */
192
193#define cheap_if(cond)  ((cond) ? 1 : 2)
194static int baseCmp(unsigned char c1, unsigned char c2) {
195    // c1,c2 == 1=A,2=C (==index of character in nucleic_acid_order[])
196    // returns  0 for equal
197    //          1 for probably mutations
198    //          2 for improbably mutations
199#define COMPARABLE_BASES 5
200
201    if (c1==c2) return 0;
202
203    if (c2<c1) { // swap
204        int c3 = c1;
205
206        c1 = c2;
207        c2 = c3;
208    }
209
210    if (c2<=COMPARABLE_BASES) {
211        fa_assert(c1<=COMPARABLE_BASES);
212
213        switch (c1) {
214            case 0: return 2;
215            case 1: return cheap_if(c2==3);                     // A->G
216            case 2: return cheap_if(c2==4 || c2==5);            // C->T/U
217            case 3: return cheap_if(c2==1);                     // G->A
218            case 4: if (c2==5) return 0;                        // T->U
219            case 5: if (c2==4) return 0;                        // U->T
220                return cheap_if(c2==2);                         // T/U->C
221            default: fa_assert(0);
222        }
223    }
224
225    int bestMatch = 3;
226    if (c1<=COMPARABLE_BASES) {
227        for (int i=1; i<=COMPARABLE_BASES; i++) {
228            if (isalpha(nucleic_maybe[i][c2])) { // 'c2' maybe a 'i'
229                int match = baseCmp(c1, i);
230                if (match<bestMatch) bestMatch = match;
231            }
232        }
233    }
234    else {
235        for (int i=1; i<=COMPARABLE_BASES; i++) {
236            if (isalpha(nucleic_maybe[i][c1])) { // 'c1' maybe a 'i'
237                int match = baseCmp(i, c2);
238                if (match<bestMatch) bestMatch = match;
239            }
240        }
241    }
242
243    fa_assert(bestMatch>=0 && bestMatch<=2);
244    return bestMatch;
245}
246#undef cheap_if
247
248int baseMatch(char c1, char c2) {
249    // c1,c2 == ACGUTRS...
250    // returns  0 for equal
251    //          1 for probably mutations
252    //          2 for improbably mutations
253    //          -1 if one char is illegal
254
255    const char *p1 = strchr(nucleic_acid_order, c1);
256    const char *p2 = strchr(nucleic_acid_order, c2);
257
258    fa_assert(c1);
259    fa_assert(c2);
260
261    if (p1 && p2) return baseCmp(p1-nucleic_acid_order, p2-nucleic_acid_order);
262    return -1;
263}
264
265static int getPenalty(char c1, char c2) {
266    // c1,c2 = A=0,C=1,... s.o.
267    switch (baseCmp(c1, c2)) {
268        case 1: return DEFAULT_PROBABLY_MUTATION;
269        case 2: return DEFAULT_IMPROBABLY_MUTATION;
270        default: break;
271    }
272
273    return 0;
274}
275
276static char *result[3];         // result buffers
277
278static char pam250mt[] = {
279    12,
280     0, 2,
281    -2, 1, 3,
282    -3, 1, 0, 6,
283    -2, 1, 1, 1, 2,
284    -3, 1, 0, -1, 1, 5,
285    -4, 1, 0, -1, 0, 0, 2,
286    -5, 0, 0, -1, 0, 1, 2, 4,
287    -5, 0, 0, -1, 0, 0, 1, 3, 4,
288    -5, -1, -1, 0, 0, -1, 1, 2, 2, 4,
289    -3, -1, -1, 0, -1, -2, 2, 1, 1, 3, 6,
290    -4, 0, -1, 0, -2, -3, 0, -1, -1, 1, 2, 6,
291    -5, 0, 0, -1, -1, -2, 1, 0, 0, 1, 0, 3, 5,
292    -5, -2, -1, -2, -1, -3, -2, -3, -2, -1, -2, 0, 0, 6,
293    -2, -1, 0, -2, -1, -3, -2, -2, -2, -2, -2, -2, -2, 2, 5,
294    -6, -3, -2, -3, -2, -4, -3, -4, -3, -2, -2, -3, -3, 4, 2, 6,
295    -2, -1, 0, -1, 0, -1, -2, -2, -2, -2, -2, -2, -2, 2, 4, 2, 4,
296    -4, -3, -3, -5, -4, -5, -4, -6, -5, -5, -2, -4, -5, 0, 1, 2, -1, 9,
297     0, -3, -3, -5, -3, -5, -2, -4, -4, -4, 0, -4, -4, -2, -1, -1, -2, 7, 10,
298    -8, -2, -5, -6, -6, -7, -4, -7, -7, -5, -3, 2, -3, -4, -5, -2, -6, 0, 0, 17
299};
300
301static char *matptr = pam250mt;
302
303#if (defined(DISPLAY_DIFF) || defined(MATRIX_DUMP))
304static void p_decode(const unsigned char *naseq, unsigned char *seq, int l) {
305    int len = strlen(amino_acid_order);
306
307    for (int i=1; i<=l && naseq[i]; i++) {
308        fa_assert(naseq[i]<len);
309        seq[i] = amino_acid_order[naseq[i]];
310    }
311}
312
313static void n_decode(const unsigned char *naseq, unsigned char *seq, int l) {
314    int len = strlen(nucleic_acid_order);
315
316    for (int i=1; i<=l && naseq[i]; i++) {
317        fa_assert(naseq[i]<len);
318        seq[i] = nucleic_acid_order[naseq[i]];
319    }
320}
321#endif
322
323inline ARB_ERROR MAXLENtooSmall() {
324    return "ClustalV-aligner: MAXLEN is dimensioned to small for this sequence";
325}
326
327inline void *ckalloc(size_t bytes, ARB_ERROR& error) {
328    if (error) return NULL;
329
330    void *ret = malloc(bytes);
331
332    if (!ret) error = "out of memory";
333    return ret;
334}
335
336static ARB_ERROR init_myers(long max_seq_length) {
337    ARB_ERROR error;
338
339    naa1 = (int **)ckalloc(MAX_BASETYPES * sizeof (int *), error);
340    naa2 = (int **)ckalloc(MAX_BASETYPES * sizeof (int *), error);
341    naas = (int **)ckalloc(MAX_BASETYPES * sizeof (int *), error);
342
343    for (int i=0; i<MAX_BASETYPES && !error; i++) {
344        naa1[i] = (int *)ckalloc((max_seq_length+1)*sizeof(int), error);
345        naa2[i] = (int *)ckalloc((max_seq_length+1)*sizeof(int), error);
346        naas[i] = (int *)ckalloc((max_seq_length+1)*sizeof(int), error);
347    }
348    return error;
349}
350
351static void make_pamo(int nv) {
352    little_pam=big_pam=matptr[0];
353    for (int i=0; i<210; ++i) {
354        int c=matptr[i];
355        little_pam=(little_pam<c) ? little_pam : c;
356        big_pam=(big_pam>c) ? big_pam : c;
357    }
358    for (int i=0; i<210; ++i) {
359        pamo[i] = matptr[i]-little_pam;
360    }
361    nv -= little_pam;
362    big_pam -= little_pam;
363    xover = big_pam - nv;
364    /*
365      fprintf(stdout,"\n\nxover= %d, big_pam = %d, little_pam=%d, nv = %d\n\n"
366      ,xover,big_pam,little_pam,nv);
367    */
368}
369
370static void fill_pam() {
371    int pos = 0;
372
373    for (int i=0; i<20; ++i) {
374        for (int j=0; j<=i; ++j) {
375            pam[i][j]=pamo[pos++];
376        }
377    }
378
379    for (int i=0; i<20; ++i) {
380        for (int j=0; j<=i; ++j) {
381            pam[j][i]=pam[i][j];
382        }
383    }
384
385    if (dnaflag) {
386        xover=4;
387        big_pam=8;
388        for (int i=0; i<=NUCLEIDS; ++i) {
389            for (int j=0; j<=NUCLEIDS; ++j) {
390                weights[i][j] = getPenalty(i, j);
391            }
392        }
393    }
394    else {
395        for (int i=1; i<MAX_BASETYPES; ++i) {
396            for (int j=1; j<MAX_BASETYPES; ++j) {
397                weights[i][j] = big_pam - pam[i-1][j-1];
398            }
399        }
400        for (int i=0; i<MAX_BASETYPES; ++i) {
401            weights[0][i] = xover;
402            weights[i][0] = xover;
403        }
404    }
405}
406
407static void exit_myers() {
408    for (int i=0; i<MAX_BASETYPES; i++) {
409        free(naas[i]);
410        free(naa2[i]);
411        free(naa1[i]);
412    }
413    free(naas);
414    free(naa2);
415    free(naa1);
416}
417
418static ARB_ERROR init_show_pair(long max_seq_length) {
419    ARB_ERROR error;
420
421#if defined(ASSERTION_USED)
422    displ_size = (2*max_seq_length + 1);
423#endif // ASSERTION_USED
424
425    displ      = (int *) ckalloc((2*max_seq_length + 1) * sizeof (int), error);
426    last_print = 0;
427
428    zza = (int *)ckalloc((max_seq_length+1) * sizeof (int), error);
429    zzb = (int *)ckalloc((max_seq_length+1) * sizeof (int), error);
430
431    zzc = (int *)ckalloc((max_seq_length+1) * sizeof (int), error);
432    zzd = (int *)ckalloc((max_seq_length+1) * sizeof (int), error);
433
434    return error;
435}
436
437static void exit_show_pair() {
438    freenull(zzd);
439    freenull(zzc);
440    freenull(zzb);
441    freenull(zza);
442    freenull(displ);
443#if defined(ASSERTION_USED)
444    displ_size = 0;
445#endif // ASSERTION_USED
446}
447
448inline int set_displ(int offset, int value) {
449    fa_assert(offset >= 0 && offset<int(displ_size));
450    displ[offset] = value;
451    return displ[offset];
452}
453inline int decr_displ(int offset, int value) {
454    fa_assert(offset >= 0 && offset<int(displ_size));
455    displ[offset] -= value;
456    return displ[offset];
457}
458
459
460inline void add(int v) {
461    // insert 'v' gaps into master ???
462    if (last_print<0 && print_ptr>0) {
463        set_displ(print_ptr-1, v);
464        set_displ(print_ptr++, last_print);
465    }
466    else {
467        last_print = set_displ(print_ptr++, v);
468    }
469}
470
471inline int calc_weight(int iat, int jat, int v1, int v2) {
472    fa_assert(pos1==1 && pos2==1);
473    unsigned char j = seq_array[alist[1]][v2+jat-1];
474    return j<128 ? naas[j][v1+iat-1] : 0;
475}
476
477#ifdef MATRIX_DUMP
478
479inline const unsigned char *lstr(const unsigned char *s, int len) {
480    static unsigned char *lstr_ss = 0;
481
482    freeset(lstr_ss, (unsigned char*)strndup((const char*)s, len));
483    return lstr_ss;
484}
485
486inline const unsigned char *nstr(unsigned char *cp, int length) {
487    const unsigned char *s = lstr(cp, length);
488    (dnaflag ? n_decode : p_decode)(s-1, const_cast<unsigned char *>(s-1), length);
489    return s;
490}
491
492inline void dumpMatrix(int x0, int y0, int breite, int hoehe, int mitte_vert) {
493    char *sl = (char*)malloc(hoehe+3);
494    char *ma = (char*)malloc(breite+3);
495
496    sprintf(ma, "-%s-", nstr(seq_array[1]+x0, breite));
497    sprintf(sl, "-%s-", nstr(seq_array[2]+y0, hoehe));
498
499    printf("                 ");
500
501    {
502        int b;
503        for (b=0; b<=mitte_vert; b++) {
504            printf("%5c", ma[b]);
505        }
506        printf("  MID");
507        for (b++; b<=breite+1+1; b++) {
508            printf("%5c", ma[b-1]);
509        }
510    }
511
512    for (int h=0; h<=hoehe+1; h++) {
513        printf("\n%c vertical:      ", sl[h]);
514        for (int b=0; b<=breite+1+1; b++) printf("%5i", vertical[b][h]);
515        printf("\n  verticalOpen:  ");
516        for (int b=0; b<=breite+1+1; b++) printf("%5i", verticalOpen[b][h]);
517        printf("\n  diagonal:      ");
518        for (int b=0; b<=breite+1+1; b++) printf("%5i", diagonal[b][h]);
519        printf("\n  horizontal:    ");
520        for (int b=0; b<=breite+1+1; b++) printf("%5i", horizontal[b][h]);
521        printf("\n  horizontalOpen:");
522        for (int b=0; b<=breite+1+1; b++) printf("%5i", horizontalOpen[b][h]);
523        printf("\n");
524    }
525
526    printf("--------------------\n");
527
528    free(ma);
529    free(sl);
530}
531#endif
532
533static int diff(int v1, int v2, int v3, int v4, int st, int en) {
534    /* v1,v3    master sequence (v1=offset, v3=length)
535     * v2,v4    slave sequence (v2=offset, v4=length)
536     * st,en    gap_open-penalties for start and end of the sequence
537     *
538     * returns  costs for inserted gaps
539     */
540#ifdef DEBUG
541# ifdef MATRIX_DUMP
542    int display_matrix = 0;
543
544    fa_assert(v3<=(DISPLAY_MATRIX_SIZE+2));         // width
545    fa_assert(v4<=(DISPLAY_MATRIX_SIZE));           // height
546
547# define DISPLAY_MATRIX_ELEMENTS ((DISPLAY_MATRIX_SIZE+2)*(DISPLAY_MATRIX_SIZE+2))
548
549    memset(vertical,       -1, DISPLAY_MATRIX_ELEMENTS*sizeof(int));
550    memset(verticalOpen,   -1, DISPLAY_MATRIX_ELEMENTS*sizeof(int));
551    memset(diagonal,       -1, DISPLAY_MATRIX_ELEMENTS*sizeof(int));
552    memset(horizontal,     -1, DISPLAY_MATRIX_ELEMENTS*sizeof(int));
553    memset(horizontalOpen, -1, DISPLAY_MATRIX_ELEMENTS*sizeof(int));
554
555# endif
556#endif
557    static int deep;
558    deep++;
559
560#if (defined (DEBUG) && 0)
561    {
562        char *d = lstr(seq_array[1]+v1, v3);
563
564        (dnaflag ? n_decode : p_decode)(d-1, d-1, v3);
565
566        for (int cnt=0; cnt<deep; cnt++) putchar(' ');
567        printf("master = '%s' (penalties left=%i right=%i)\n", d, st, en);
568
569        d = lstr(seq_array[2]+v2, v4);
570        (dnaflag ? n_decode : p_decode)(d-1, d-1, v4);
571
572        for (cnt=0; cnt<deep; cnt++) putchar(' ');
573        printf("slave  = '%s'\n", d);
574    }
575#endif
576
577    if (v4<=0) {                                                // if slave sequence is empty
578        if (v3>0) {
579            if (last_print<0 && print_ptr>0) {
580                last_print = decr_displ(print_ptr-1, v3);         // add ..
581            }
582            else {
583                last_print = set_displ(print_ptr++, -(v3));       // .. or insert gap of length 'v3' into slave
584            }
585        }
586
587        deep--;
588        return master_gapAt(v1, v3);
589    }
590
591    int ctrj = 0;
592    if (v3<=1) {
593        if (v3<=0) {                                            // if master sequence is empty
594            add(v4);                                            // ??? insert gap length 'v4' into master ???
595            deep--;
596            return slave_gapAt(v2, v4);
597        }
598
599        fa_assert(v3==1);
600        // if master length == 1
601
602        if (v4==1) {
603            if (st>en) st=en;
604
605            //!*************if(!v4)*********BUG*******************************
606
607            int ctrc = slave_gapAtWithOpenPenalty(v2, v4, st);
608            ctrj = 0;
609
610            for (int j=1; j<=v4; ++j) {
611                int k = slave_gapAt(v2, j-1) + calc_weight(1, j, v1, v2) + slave_gapAt(v2+j, v4-j);
612                if (k<ctrc) { ctrc = k; ctrj = j; }
613            }
614
615            if (!ctrj) {
616                add(v4);
617                if (last_print<0 && print_ptr>0) {
618                    last_print = decr_displ(print_ptr-1, 1);
619                }
620                else {
621                    last_print = set_displ(print_ptr++, -1);
622                }
623            }
624            else {
625                if (ctrj>1) add(ctrj-1);
626                set_displ(print_ptr++, last_print = 0);
627                if (ctrj<v4) add(v4-ctrj);
628            }
629
630            deep--;
631            return ctrc;
632        }
633    }
634
635    fa_assert(v3>=1 && v4>=1);
636
637    // slave length  >= 1
638    // master length >= 1
639
640# define AF 0
641
642    // first column of matrix (slave side):
643    IF_MATRIX_DUMP(vertical[0][0]=)
644        zza[0] = 0;
645    int p = master_gap_open(v1);
646    for (int j=1; j<=v4; j++) {
647        p += master_gap_extend(v1);
648        IF_MATRIX_DUMP(vertical[0][j]=)
649            zza[j] = p;
650        IF_MATRIX_DUMP(verticalOpen[0][j]=)
651            zzb[j] = p + master_gap_open(v1);
652    }
653
654    // left half of the matrix
655    p = st;
656    int ctri = v3 / 2;
657    for (int i=1; i<=ctri; i++) {
658        int n = zza[0];
659        p += master_gap_extend(v1+i+AF);
660        int k = p;
661        IF_MATRIX_DUMP(vertical[i][0]=)
662            zza[0] = k;
663        int l = p + master_gap_open(v1+i+AF);
664
665        for (int j=1; j<=v4; j++) {
666            // from above (gap in master (behind position i))
667            IF_MATRIX_DUMP(verticalOpen[i][j]=)         k += master_gap_open(v1+i+AF)+master_gap_extend(v1+i+AF);       // (1)
668            IF_MATRIX_DUMP(vertical[i][j]=)             l += master_gap_extend(v1+i+AF);                                // (2)
669            if (k<l) l                                     = k;                                                         // l=min((1),(2))
670
671            // from left (gap in slave (behind position j))
672            IF_MATRIX_DUMP(horizontalOpen[i][j]=)       k = zza[j] + slave_gap_open(v2+j+AF)+slave_gap_extend(v2+j+AF);         // (3)
673            int m;
674            IF_MATRIX_DUMP(horizontal[i][j]=)           m = zzb[j] + slave_gap_extend(v2+j+AF);                                 // (4)
675            if (k<m) m = k;     // m=min((3),(4))
676
677            // diagonal (no gaps)
678            IF_MATRIX_DUMP(diagonal[i][j]=)             k = n + calc_weight(i, j, v1, v2);                              // (5)
679            if (l<k) k = l;
680            if (m<k) k = m;     // k = minimum of all paths
681
682            n      = zza[j];    // minimum of same row; one column to the left
683            zza[j] = k;         // minimum of all paths to this matrix position
684            zzb[j] = m;         // minimum of those two paths, where gap was inserted into slave
685        }
686    }
687
688# define MHO 1  // X-Offset for second half of matrix-arrays (only used to insert MID-column into dumpMatrix())
689# define BO  1  // Offset for second half of matrix (cause now the indices i and j are smaller as above)
690
691    // last column of matrix:
692
693    zzb[0]=zza[0];
694    IF_MATRIX_DUMP(vertical[v3+1+MHO][v4+1]=)
695        zzc[v4]=0;
696    p = master_gap_open(v1+v3);
697    for (int j=v4-1; j>-1; j--) {
698        p += master_gap_extend(v1+v3);
699        IF_MATRIX_DUMP(vertical[v3+1+MHO][j+BO]=)
700            zzc[j] = p;
701        IF_MATRIX_DUMP(verticalOpen[v3+1+MHO][j+BO]=)
702            zzd[j] = p+master_gap_open(v1+v3);
703    }
704
705    // right half of matrix (backwards):
706    p = en;
707    for (int i=v3-1; i>=ctri; i--) {
708        int n = zzc[v4];
709        p += master_gap_extend(v1+i);
710        int k = p;
711        IF_MATRIX_DUMP(vertical[i+BO+MHO][v4+1]=)
712            zzc[v4] = k;
713        int l = p+master_gap_open(v1+i);
714
715        for (int j=v4-1; j>=0; j--) {
716            // from below (gap in master (in front of position (i+BO)))
717            IF_MATRIX_DUMP(verticalOpen[i+BO+MHO][j+BO]=)       k += master_gap_open(v1+i)+master_gap_extend(v1+i);     // (1)
718            IF_MATRIX_DUMP(vertical[i+BO+MHO][j+BO]=)           l += master_gap_extend(v1+i);                           // (2)
719            if (k<l) l                                             = k;                                                 // l=min((1),(2))
720
721            // from right (gap in slave (in front of position (j+BO)))
722            IF_MATRIX_DUMP(horizontalOpen[i+BO+MHO][j+BO]=)     k = zzc[j] + slave_gap_open(v2+j) + slave_gap_extend(v2+j);     // (3)
723            int m;
724            IF_MATRIX_DUMP(horizontal[i+BO+MHO][j+BO]=)         m = zzd[j] + slave_gap_extend(v2+j);                            // (4)
725            if (k<m) m = k;     // m=min((3),(4))
726
727            // diagonal (no gaps)
728            IF_MATRIX_DUMP(diagonal[i+BO+MHO][j+BO]=)           k = n + calc_weight(i+BO, j+BO, v1, v2);                // (5)
729            if (l<k) k = l;
730            if (m<k) k = m;     // k = minimum of all paths
731
732            n      = zzc[j];    // minimum of same row; one column to the right
733            zzc[j] = k;         // minimum of all paths to this matrix position
734            zzd[j] = m;         // minimum of those two paths, where gap was inserted into slave
735        }
736    }
737
738#undef BO
739
740    // connect rightmost column of first half (column ctri)
741    // with leftmost column of second half (column ctri+1)
742
743    zzd[v4] = zzc[v4];
744
745    int ctrc = INT_MAX;
746    int flag = 0;
747
748    for (int j=(ctri ? 0 : 1); j<=v4; j++) {
749        int k;
750        IF_MATRIX_DUMP(vertical[ctri+MHO][j]=)
751            k = zza[j] + zzc[j];                // sum up both calculations (=diagonal=no gaps)
752
753        if (k<ctrc || ((k==ctrc) && (zza[j]!=zzb[j]) && (zzc[j]==zzd[j]))) {
754            ctrc = k;
755            ctrj = j;
756            flag = 1;
757        }
758    }
759
760    for (int j=v4; j>=(ctri ? 0 : 1); j--) {
761        int k;
762        IF_MATRIX_DUMP(verticalOpen[ctri+MHO][j]=)
763            k = zzb[j] + zzd[j]                 // paths where gaps were inserted into slave (left and right side!)
764            - slave_gap_open(j);                // subtract gap_open-penalty which was added twice (at left and right end of gap)
765
766        if (k<ctrc) {
767            ctrc = k;
768            ctrj = j;
769            flag = 2;
770        }
771    }
772
773    fa_assert(flag>=1 && flag<=2);
774
775#undef MHO
776
777#ifdef MATRIX_DUMP
778    if (display_matrix) {
779        dumpMatrix(v1, v2, v3, v4, ctri);
780    }
781#endif
782
783    // Conquer recursively around midpoint
784
785    if (flag==1) { // Type 1 gaps
786        diff(v1, v2, ctri, ctrj, st, master_gap_open(v1+ctri));                 // includes midpoint ctri and ctrj
787        diff(v1+ctri, v2+ctrj, v3-ctri, v4-ctrj, master_gap_open(v1+ctri), en);
788    }
789    else {
790        diff(v1, v2, ctri-1, ctrj, st, 0);                                      // includes midpoint ctrj
791
792        if (last_print<0 && print_ptr>0) { // Delete 2
793            last_print = decr_displ(print_ptr-1, 2);
794        }
795        else {
796            last_print = set_displ(print_ptr++, -2);
797        }
798
799        diff(v1+ctri+1, v2+ctrj, v3-ctri-1, v4-ctrj, 0, en);
800    }
801
802    deep--;
803    return ctrc;       // Return the score of the best alignment
804}
805
806static void do_align(int& score, long act_seq_length) {
807    pos1 = pos2 = 0;
808
809    // clear statistics
810
811    for (int i=1; i<=act_seq_length; ++i) {
812        for (int j=0; j<MAX_BASETYPES; ++j) {
813            naa1[j][i] = naa2[j][i]=0;
814            naas[j][i]=0;
815        }
816    }
817
818    // create position statistics for each group
819    // [here every group contains only one seq]
820
821    int l1 = 0;
822    int l2 = 0;
823    for (int i=1; i<=nseqs; ++i) {
824        if (group[i]==1) {
825            fst_list[++pos1]=i;
826            for (int j=1; j<=seqlen_array[i]; ++j) {
827                unsigned char b = seq_array[i][j];
828                if (b<128) {
829                    ++naa1[b][j];
830                    ++naa1[0][j];
831                }
832            }
833            if (seqlen_array[i]>l1) l1=seqlen_array[i];
834        }
835        else if (group[i]==2) {
836            snd_list[++pos2]=i;
837            for (int j=1; j<=seqlen_array[i]; ++j) {
838                unsigned char b = seq_array[i][j];
839                if (b<128) {
840                    ++naa2[b][j];
841                    ++naa2[0][j];
842                }
843            }
844            if (seqlen_array[i]>l2) l2=seqlen_array[i];
845        }
846    }
847
848    if (pos1>=pos2) {
849        int t_arr[MAX_BASETYPES];
850
851        for (int i=1; i<=pos2; ++i) alist[i]=snd_list[i];
852        for (int n=1; n<=l1; ++n) {
853            for (int i=1; i<MAX_BASETYPES; ++i) t_arr[i]=0;
854            for (int i=1; i<MAX_BASETYPES; ++i) {
855                if (naa1[i][n]>0) {
856                    for (int j=1; j<MAX_BASETYPES; ++j) {
857                        t_arr[j] += (weights[i][j]*naa1[i][n]);
858                    }
859                }
860            }
861            int k = naa1[0][n];
862            if (k>0) {
863                for (int i=1; i<MAX_BASETYPES; ++i) {
864                    naas[i][n]=t_arr[i]/k;
865                }
866            }
867        }
868    }
869    else {
870        fa_assert(0);       // should never occur if MAXN==2
871    }
872
873    score = diff(1, 1, l1, l2, master_gap_open(1), master_gap_open(l1+1)); // Myers and Miller alignment now
874}
875
876static int add_ggaps(long /* max_seq_length */) {
877    int pos   = 1;
878    int to_do = print_ptr;
879
880    for (int i=0; i<to_do; ++i) { // was: 1 .. <=to_do
881        if (displ[i]==0) {
882            result[1][pos]=result[2][pos]='*';
883            ++pos;
884        }
885        else {
886            int k = displ[i];
887            if (k>0) {
888                for (int j=0; j<=k-1; ++j) {
889                    result[2][pos+j]='*';
890                    result[1][pos+j]='-';
891                }
892                pos += k;
893            }
894            else {
895                k = (displ[i]<0) ? displ[i] * -1 : displ[i];
896                for (int j=0; j<=k-1; ++j) {
897                    result[1][pos+j]='*';
898                    result[2][pos+j]='-';
899                }
900                pos += k;
901            }
902        }
903    }
904
905#ifdef DEBUG
906    result[1][pos] = 0;
907    result[2][pos] = 0;
908#endif
909
910    return pos-1;
911}
912
913
914static int res_index(const char *t, char c) {
915    if (t[0]==c) return UNKNOWN_ACID;
916    for (int i=1; t[i]; i++) {
917        if (t[i]==c) return i;
918    }
919    return 0;
920}
921
922static void p_encode(const unsigned char *seq, unsigned char *naseq, int l) {
923    // code seq as ints .. use -2 for gap
924    bool warned = false;
925
926    for (int i=1; i<=l; i++) {
927        int c = res_index(amino_acid_order, seq[i]);
928
929        if (!c) {
930            if (seq[i] == '*') {
931                c = -2;
932            }
933            else {
934                if (!warned && seq[i] != '?') { // consensus contains ? and it means X
935                    char buf[100];
936                    sprintf(buf, "Illegal character '%c' in sequence data", seq[i]);
937                    aw_message(buf);
938                    warned = true;
939                }
940                c = res_index(amino_acid_order, 'X');
941            }
942        }
943
944        fa_assert(c>0 || c == -2);
945        naseq[i] = c;
946    }
947}
948
949static void n_encode(const unsigned char *seq, unsigned char *naseq, int l) {
950    // code seq as ints .. use -2 for gap
951    bool warned = false;
952    for (int i=1; i<=l; i++) {
953        int c = res_index(nucleic_acid_order, seq[i]);
954
955        if (!c) {
956            if (!warned && seq[i] != '?') {
957                char buf[100];
958                sprintf(buf, "Illegal character '%c' in sequence data", seq[i]);
959                aw_message(buf);
960                warned = true;
961            }
962            c = res_index(nucleic_acid_order, 'N');
963        }
964
965        fa_assert(c>0);
966        naseq[i] = c;
967    }
968}
969
970ARB_ERROR ClustalV_align(int           is_dna,
971                         int           weighted,
972                         const char   *seq1,
973                         int           length1,
974                         const char   *seq2,
975                         int           length2,
976                         const int    *gapsBefore1, // size of array = length1+1
977                         int           max_seq_length,
978                         // result params:
979                         const char*&  res1,
980                         const char*&  res2,
981                         int&          reslen,
982                         int&          score)
983{
984    ARB_ERROR error;
985    gaps_before_position = gapsBefore1;
986
987    if (!module_initialized) { // initialize only once
988        dnaflag = is_dna;
989        is_weight = weighted;
990
991        error             = init_myers(max_seq_length);
992        if (!error) error = init_show_pair(max_seq_length);
993        make_pamo(0);
994        fill_pam();
995
996        for (int i=1; i<=2 && !error; i++) {
997            seq_array[i] = (unsigned char*)ckalloc((max_seq_length+2)*sizeof(char), error);
998            result[i] = (char*)ckalloc((max_seq_length+2)*sizeof(char), error);
999            group[i] = i;
1000        }
1001
1002        if (!error) module_initialized = true;
1003    }
1004    else {
1005        if (dnaflag!=is_dna || is_weight!=weighted) {
1006            error = "Please call ClustalV_exit() between calls that differ in\n"
1007                "one of the following parameters:\n"
1008                "       is_dna, weighted";
1009        }
1010    }
1011
1012    if (!error) {
1013        nseqs = 2;
1014        print_ptr = 0;
1015
1016#if defined(DEBUG) && 1
1017        memset(&seq_array[1][1], 0, max_seq_length*sizeof(seq_array[1][1]));
1018        memset(&seq_array[2][1], 0, max_seq_length*sizeof(seq_array[2][1]));
1019        memset(&displ[1], 0xff, max_seq_length*sizeof(displ[1]));
1020        seq_array[1][0] = '_';
1021        seq_array[2][0] = '_';
1022#endif
1023
1024        {
1025            typedef void (*encoder)(const unsigned char*, unsigned char*, int);
1026            encoder encode = dnaflag ? n_encode : p_encode;
1027
1028            encode((const unsigned char*)(seq1-1), seq_array[1], length1);
1029            seqlen_array[1] = length1;
1030            encode((const unsigned char*)(seq2-1), seq_array[2], length2);
1031            seqlen_array[2] = length2;
1032
1033            do_align(score, max(length1, length2));
1034            int alignedLength = add_ggaps(max_seq_length);
1035
1036            result[1][alignedLength+1] = 0;
1037            result[2][alignedLength+1] = 0;
1038
1039            res1   = result[1]+1;
1040            res2   = result[2]+1;
1041            reslen = alignedLength;
1042        }
1043    }
1044
1045    return error;
1046}
1047
1048void ClustalV_exit() {
1049    if (module_initialized) {
1050        module_initialized = false;
1051
1052        for (int i=1; i<=2; i++) {
1053            free(result[i]);
1054            free(seq_array[i]);
1055        }
1056        exit_show_pair();
1057        exit_myers();
1058    }
1059}
1060
1061// --------------------------------------------------------------------------------
1062
1063#ifdef UNIT_TESTS
1064#ifndef TEST_UNIT_H
1065#include <test_unit.h>
1066#endif
1067
1068static arb_test::match_expectation clustal_aligns(const char *i1, const char *i2, const char *expected1, const char *expected2, int expected_score) {
1069    size_t l1 = strlen(i1);
1070    size_t l2 = strlen(i2);
1071
1072    const char *result1 = NULL;
1073    const char *result2 = NULL;
1074    int         result_len; // test result value?
1075    int         score;
1076
1077    int gaps_size = l1+1;
1078    int no_gaps_before1[gaps_size];
1079    memset(no_gaps_before1, 0, gaps_size*sizeof(*no_gaps_before1));
1080
1081    ARB_ERROR error = ClustalV_align(0, 0,
1082                                     i1, l1,
1083                                     i2, l2,
1084                                     no_gaps_before1,
1085                                     max(l1, l2),
1086                                     result1,
1087                                     result2,
1088                                     result_len,
1089                                     score);
1090
1091    using namespace   arb_test;
1092    expectation_group expected;
1093
1094    expected.add(doesnt_report_error(error.deliver()));
1095    expected.add(that(result1).is_equal_to(expected1));
1096    expected.add(that(result2).is_equal_to(expected2));
1097    expected.add(that(score).is_equal_to(expected_score));
1098
1099    return all().ofgroup(expected);
1100}
1101
1102#define TEST_CLUSTAL_ALIGNS(i1,i2,o1,o2,score) TEST_EXPECTATION(clustal_aligns(i1,i2,o1,o2,score))
1103
1104void TEST_clustalV() {
1105    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1106                        "ACGTTGCAACGT",
1107                        "************",
1108                        "************",
1109                        138);
1110
1111    TEST_CLUSTAL_ALIGNS("ACGTTAACGT",
1112                        "ACGTTGCAACGT",
1113                        "*****--*****",
1114                        "************",
1115                        207);
1116   
1117    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1118                        "ACCAACGT",
1119                        "************",
1120                        "*----*******",
1121                        156);
1122
1123    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1124                        "AGCTGTCACAGT",
1125                        "************",
1126                        "************",
1127                        187);
1128
1129    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1130                        "ACGTACGTACGT",
1131                        "************",
1132                        "************",
1133                        164);
1134
1135    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1136                        "ACGTACGT",
1137                        "************",
1138                        "****----****",
1139                        162);
1140
1141    TEST_CLUSTAL_ALIGNS("ACGTTGCAACGT",
1142                        "AACGGTTC",
1143                        "************",
1144                     // "ACGTTGCAACGT",
1145                     // "A----ACGGTTC",
1146                        "*----*******",
1147                        193);
1148}
1149
1150#endif // UNIT_TESTS
1151
1152// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.