ARB
island_hopping.cxx
Go to the documentation of this file.
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 
22 IslandHoppingParameter *IslandHopping::para = 0;
23 
24 IslandHoppingParameter::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 
52 IslandHoppingParameter::~IslandHoppingParameter() {
53  // @@ uninit();
54 }
55 
56 
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=NULp;
72  char *YY=NULp;
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 }
const char * GB_ERROR
Definition: arb_core.h:25
static void do_align(int &score, long act_seq_length)
Definition: ClustalV.cxx:807
static int J
Definition: align.cxx:489
#define freeBlock(v)
Definition: mem.h:25
#define NULp
Definition: cxxforward.h:114
#define ih_assert(bed)
void Align(int nX, char X[], int secX[], char **XX, int nY, char Y[], int secY[], char **YY, int freqs, double fT, double fC, double fA, double fG, double rTC, double rTA, double rTG, double rCA, double rCG, double rAG, double dist, double supp, double gapA, double gapB, double gapC, double thres)
Definition: align.cxx:799
static int K
Definition: align.cxx:489