source: branches/profile/ISLAND_HOPPING/island_hopping.cxx

Last change on this file was 8725, checked in by westram, 7 years ago
  • comments cleanup
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.6 KB
Line 
1// ============================================================= //
2//                                                               //
3//   File      : island_hopping.cpp                              //
4//   Purpose   :                                                 //
5//                                                               //
6//   Institute of Microbiology (Technical University Munich)     //
7//   http://www.arb-home.de/                                     //
8//                                                               //
9// ============================================================= //
10
11#include "island_hopping.h"
12
13#define EXTERN
14#include "i-hopper.h"
15#include "mem.h"
16
17#ifndef ARB_ASSERT_H
18#include <arb_assert.h>
19#endif
20#define ih_assert(bed) arb_assert(bed)
21
22IslandHoppingParameter *IslandHopping::para = 0;
23
24IslandHoppingParameter::IslandHoppingParameter(bool    use_user_freqs_,
25                                               double fT_, double fC_, double fA_, double fG_,
26                                               double rTC_, double rTA_, double rTG_, double rCA_, double rCG_, double rAG_,
27                                               double dist_, double supp_, double gapA_, double gapB_, double gapC_, double thres_)
28{
29    use_user_freqs = use_user_freqs_;
30    fT    = fT_;
31    fC    = fC_;
32    fA    = fA_;
33    fG    = fG_;
34
35    rTC   = rTC_;
36    rTA   = rTA_;
37    rTG   = rTG_;
38    rCA   = rCA_;
39    rCG   = rCG_;
40    rAG   = rAG_;
41
42    dist  = dist_;
43    supp  = supp_;
44    gapA   = gapA_;
45    gapB   = gapB_;
46    gapC   = gapC_;
47    thres = thres_;
48
49    //@@@ init();
50}
51
52IslandHoppingParameter::~IslandHoppingParameter() {
53    // @@ uninit();
54}
55
56
57GB_ERROR IslandHopping::do_align() {
58
59    if (!para) {
60        para = new IslandHoppingParameter(0, 0.25, 0.25, 0.25, 0.25, 0, 4.0, 1.0, 1.0, 1.0, 1.0, 4.0, 0.3, 0.5, 8.0, 4.0, 0.001);
61    }
62
63    int   nX;
64    char *X;
65    int  *secX;
66
67    int nY;
68    char *Y;
69    int *secY;
70
71    char *XX=NULL;
72    char *YY=NULL;
73
74    int i,j,k,J,K,LJ, LK;
75
76    Error = 0;
77
78    nX = 0; nY=0;
79    for(i=0;i<alignment_length;i++) {
80        if(ref_sequence[i]!='-' && ref_sequence[i] != '.') nX++;
81        if(toAlign_sequence[i]!='-' && toAlign_sequence[i] != '.') nY++;
82    }
83
84    X=(char *)malloc((nX+1)*sizeof(char));
85    secX=(int *)malloc((nX)*sizeof(int));
86
87    Y    = (char *)malloc((nY+1)*sizeof(char));
88    secY = (int *)malloc((nY)*sizeof(int));
89
90    // @@@ helix?
91
92    j = 0; J=0; LJ = 0;
93    k = 0; K=0; LK=0;
94
95#if defined(DEBUG)
96    printf("ref_helix     = '%s'\n", ref_helix);
97    printf("toAlign_helix = '%s'\n", toAlign_helix);
98#endif // DEBUG
99
100    for(i=0;i<alignment_length;i++) {
101
102        if(ref_sequence[i]!='-' && ref_sequence[i]!='.') {
103            X[j] = ref_sequence[i];
104            if (ref_helix) {
105                switch(ref_helix[i]) {
106                    case '-': case '.':
107                        if(LJ!=0) J++;
108                        LJ = 0;
109                        break;
110                    case '[': case '<': case '(': case '{':
111                        if(LJ!=1) J++;
112                        LJ = 1;
113                        break;
114                    case ']': case '>': case ')': case '}':
115                        if(LJ!=2) J++;
116                        LJ = 2;
117                        break;
118                    default:
119                        printf("Unknown '%c'\n", ref_helix[i]);
120                        ih_assert(0);
121                        break;
122                }
123            }
124
125            secX[j]=LJ?J:0;
126            j++;
127        }
128        if(toAlign_sequence[i]!='-' && toAlign_sequence[i]!='.') {
129            Y[k] = toAlign_sequence[i];
130            if (toAlign_helix) {
131                switch(toAlign_helix[i]) {
132                    case '-': case '.':
133                        if(LK!=0) K++;
134                        LK=0;
135                        break;
136                    case '[': case '<': case '(': case '{':
137                        if(LK!=1) K++;
138                        LK=1;
139                        break;
140                    case ']': case '>': case ')': case '}':
141                        if(LK!=2) K++;
142                        LK=2;
143                        break;
144                    default:
145                        printf("Unknown '%c'\n", toAlign_helix[i]);
146                        ih_assert(0);
147                        break;
148                }
149            }
150            secY[k]=LK?K:0;
151            k++;
152        }
153    }
154    X[j]='\0'; Y[k]='\0';
155
156    if(output_sequence) {delete output_sequence; output_sequence=0;}
157    if(aligned_ref_sequence) {delete aligned_ref_sequence; aligned_ref_sequence=0;}
158
159    Align(
160          nX,X,secX,&XX,nY,Y,secY,&YY,
161          para->use_user_freqs,para->fT,para->fC,para->fA,para->fG,
162          para->rTC,para->rTA,para->rTG,para->rCA,para->rCG,para->rAG,
163          para->dist,para->supp,para->gapA,para->gapB,para->gapC,para->thres
164          );
165
166    if(!Error) {
167        int nXY                 = strlen(XX);
168        int o;
169        output_alignment_length = nXY;
170
171        {
172            FILE *fp;
173            fp = fopen("alignment.txt","w");
174            for(o=0;o<nXY;o++) fprintf(fp,"%c",XX[o]); fprintf(fp,"\n");
175            for(o=0;o<nXY;o++) fprintf(fp,"%c",YY[o]); fprintf(fp,"\n");
176            fclose(fp);
177        }
178
179        aligned_ref_sequence = new char[nXY+1];
180        output_sequence      = new char[nXY+1];
181
182        for (o = 0;o<nXY;++o) {
183            aligned_ref_sequence[o] = XX[o] == '-' ? '-' : '*';
184            output_sequence[o]      = YY[o] == '-' ? '-' : '*';
185        }
186        aligned_ref_sequence[o] = 0;
187        output_sequence[o]      = 0;
188    }
189
190    free(X);
191    free(secX);
192
193    free(Y);
194    free(secY);
195
196    freeBlock(&XX);
197    freeBlock(&YY);
198
199    return(Error);
200}
Note: See TracBrowser for help on using the repository browser.