IA32 and x64 tcoffee binaries are added
[jabaws.git] / binaries / src / tcoffee / t_coffee_source / diagonal.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4
5 #include "io_lib_header.h"
6 #include "util_lib_header.h"
7 #include "fastal_lib_header.h"
8
9
10
11 /*!
12  *      \file diagonal.c
13  *      \brief Source code for the calculation and preparation of diagonals
14  */
15
16
17
18
19
20
21
22 int 
23 diagonal_compare (const void * a, const void * b)
24 {
25         
26         int tmp = ((((Diagonal*)a)->y - ((Diagonal*)a)->x) - (((Diagonal*)b)->y - ((Diagonal*)b)->x));
27         if (tmp != 0)
28                 return tmp;
29         return (((Diagonal*)a)->y - ((Diagonal*)b)->y);
30 }
31
32 int
33 max(int a, int b)
34 {
35         if (a < b)
36                 return b;
37         else
38                 return a;
39 }
40
41
42
43 Segment*
44 extend_diagonals(Diagonal *diagonals, int *num_diagonal, int l1, int l2)
45 {
46 //      sort diagonal by diagonal number
47         int i;
48 //      printf("INPUT_PRE\n");
49 //      for (i = 0; i < num_diagonals; ++i)
50 //      {
51 //              printf("%i %i %i\n",diagonals[i].x, diagonals[i].y, diagonals[i].length);
52 //      }
53         int num_diagonals = *num_diagonal;
54         qsort (diagonals, num_diagonals, sizeof(Diagonal), diagonal_compare);
55
56
57 //      find nearest diagonal and expand
58 //      make shure that overlapping segments on the same diagonal are merged
59
60         int diff;
61 //      printf("INPUT\n");
62         for (i = 0; i < num_diagonals; ++i)
63         {
64                 printf("%i %i %i\n",diagonals[i].x, diagonals[i].y, diagonals[i].length);
65         }
66         for (i = 0; i < num_diagonals; ++i)
67         {
68                 diagonals[i].end_exp = 0;
69                 diagonals[i].front_exp = 0;
70         }
71         int first = 0;
72         int next = 1;
73         while (next < num_diagonals)
74         {
75
76                 if (!((diagonals[next].y >= diagonals[first].y) && diagonals[next].x <= diagonals[first].x))
77                 {
78                         printf("%i %i %i %i %i %i\n",diagonals[first].x, diagonals[first].y, diagonals[next].x, diagonals[next].y, diagonals[first].y - diagonals[first].x, diagonals[next].y - diagonals[next].x);
79                         diff = diagonals[next].y - (diagonals[first].y + diagonals[first].length);
80                         if (diff > 0)
81                         {
82                                 if ((diagonals[first].y - diagonals[first].x) != ((diagonals[next].y - diagonals[next].x)))
83                                 {
84                                         diagonals[first].end_exp = max(diagonals[first].end_exp, diff);
85                                         diagonals[next].front_exp = max(diagonals[next].front_exp, diff);
86         //                              if (diagonals[next].x - diagonals[next].front_exp < 0)
87         
88                                         while (diagonals[++first].x == -1);
89                                 }
90                                 else
91                                 {
92                                         printf("ICH TUWAS\n");
93                                         diagonals[first].length = diagonals[next].y + diagonals[next].length - diagonals[first].y;
94                                         diagonals[first].end_exp = diagonals[next].end_exp;
95                                         diagonals[next].x = -1;
96                                 }
97                         }
98                         diff = diagonals[first].x - (diagonals[next].x + diagonals[next].length);
99                         if (diff > 0)
100                         {
101         //                      diff = diagonals[next].x + diagonals[next].length - diagonals[first].x;
102         
103                                 diagonals[first].front_exp =  max(diagonals[first].front_exp, diff);
104                                 diagonals[next].end_exp = max(diagonals[next].end_exp, diff);
105         //                      ++num_segments;
106         //                      if (diagonals[first].x - diagonals[first].front_exp < 0)
107         
108                                 while (diagonals[++first].x == -1);
109                         }
110                 } else
111                         ++ first;
112                 ++next;
113                 
114         }
115
116         int num_segments =0;
117         Segment *seg = vcalloc(num_diagonals, sizeof(Segment));
118         int pos = 0;
119 //      int diag_num1, diag_num2;
120         int *tmp1;
121         Diagonal *tmp2;
122         for (i = 0; i < num_diagonals; ++i)
123         {
124                 if (diagonals[i].x != -1)
125                 {
126                         ++num_segments;
127                         tmp2 = &diagonals[i];
128                         seg[pos].segments = vcalloc(6, sizeof(int));
129                         seg[pos].current_pos = seg[pos].segments;
130                         tmp1 = seg[pos].segments;
131                         tmp1[0] = tmp2->x - tmp2->front_exp;
132                         tmp1[1] = tmp2->y - tmp2->front_exp;
133                         tmp1[2] = tmp2->front_exp + tmp2->length + tmp2->end_exp;
134                         tmp1[5] = 1;
135                         tmp1[4] = l2;
136                         tmp1[3] = tmp1[0] + (l2 -  tmp1[1]);
137                         ++pos;
138                 }
139         }
140
141         for ( i = 0; i < num_segments; ++i )
142         {
143                 printf("%i %i %i || %i %i %i\n", seg[i].current_pos[0], seg[i].current_pos[1], seg[i].current_pos[2], seg[i].current_pos[3], seg[i].current_pos[4], seg[i].current_pos[5]);
144         }
145 //      exit(1);
146         vrealloc(seg,num_segments*sizeof(Segment));
147         
148         *num_diagonal = num_segments;
149         
150 //      vfree(diagonals);
151         return seg;
152 }
153
154
155 int **
156 segments2int(Segment *diagonals,
157                         int num_diagonals,
158                         char *seq1,
159                         char *seq2,
160                         Fastal_profile *profile1,
161                         Fastal_profile *profile2,
162                         int *num_points,
163                         Fastal_param *param_set)
164 {
165 //      printf("NUM: %i\n", num_diagonals);
166         int l1 = strlen(seq1);
167         int l2 = strlen(seq2);
168         int gep = param_set->gep;
169
170         int dig_length;
171         if (seq1 > seq2)
172                 dig_length = l1;
173         else
174                 dig_length = l2;
175
176         int current_size = num_diagonals*dig_length + l1 +l2;
177
178         int **list = vcalloc(current_size, sizeof(int*));
179 //      int *diags = vcalloc(num_diagonals, sizeof(int));
180         int i;
181 //      for (i = 0; i < num_diagonals; ++i)
182 //      {
183 //              diags[i] = l1 - diagonals[i*3] + diagonals[i*3+1];
184 //      }
185
186 //      qsort (diags, num_diagonals, sizeof(int), fastal_compare);
187
188
189 //      int *diagx = vcalloc(num_diagonals, sizeof(int));
190 //      int *diagy = vcalloc(num_diagonals, sizeof(int));
191
192
193         //+1 because diagonals start here at position 1, like in "real" dynamic programming
194         int a = -1, b = -1;
195         for (i = 0; i < num_diagonals; ++i)
196         {
197                 if (diagonals[i].current_pos[1] - diagonals[i].current_pos[0] < 0)
198                 {
199 //                      diagx[i] = l1 - diags[i];
200 //                      diagy[i] = 0;
201                         a= i;
202                 }
203                 else
204                         break;
205         }
206         ++a;
207         b=a-1;
208         for (; i < num_diagonals; ++i)
209         {
210 //              diagx[i] = 0;
211 //              diagy[i] = diags[i]-l1;
212                 diagonals[i].prev_position = diagonals[i].current_pos[1] - diagonals[i].current_pos[0];
213                 b = i;
214         }
215         
216 //      vfree(diags);
217         int tmpy_pos;
218         int tmpy_value;
219         int **M = param_set->M;
220         int *last_y = vcalloc(l2+1, sizeof(int));
221         int *last_x = vcalloc(l1+1, sizeof(int));
222         last_y[0] = 0;
223
224         last_x[0] = 0;
225         list[0] = vcalloc(6, sizeof(int));
226
227         int list_pos = 1;
228         int dig_num = l1;
229         int tmp_l2 = l2 + 1;
230
231         //left border
232         for (; list_pos < tmp_l2; ++list_pos)
233         {
234                 list[list_pos] = vcalloc(6, sizeof(int));
235                 list[list_pos][0] = 0;
236                 list[list_pos][1] = list_pos;
237                 last_y[list_pos] = list_pos;
238                 list[list_pos][2] = list_pos*gep;
239                 list[list_pos][4] = list_pos-1;
240         }
241
242         int pos_x = 0;
243         int y;
244         int tmp_l1 = l1-1;
245         int tmpx, tmpy;
246         while (pos_x < tmp_l1)
247         {
248                 if (list_pos + num_diagonals+2 > current_size)
249                 {
250                         current_size += num_diagonals*1000;
251                         list = vrealloc(list, current_size * sizeof(int*));
252                 }
253                 //upper border
254                 list[list_pos] = vcalloc(6, sizeof(int));
255                 list[list_pos][0] = ++pos_x;
256                 list[list_pos][1] = 0;
257                 list[list_pos][2] = pos_x * gep;
258                 list[list_pos][3] = last_y[0];
259 //              tmpy_value = list_pos;
260 //              tmpy_pos = 0;
261                 last_x[pos_x] = list_pos;
262                 ++list_pos;
263
264                 //diagonals
265                 for (i = a; i <= b; ++i)
266                 {
267                         
268                         if (pos_x == diagonals[i].current_pos[0])
269                         {
270                                 if (diagonals[i].current_pos[2] == 0)
271                                 {
272                                         diagonals[i].current_pos+=3;
273                                 }
274                                 list[list_pos] = vcalloc(6, sizeof(int));
275                                 tmpx = diagonals[i].current_pos[0];
276                                 tmpy = diagonals[i].current_pos[1];
277                                 list[list_pos][0] = tmpx;
278                                 list[list_pos][1] = tmpy;
279                                 list[list_pos][3] = last_y[tmpy];
280                                 list[list_pos][4] = list_pos-1;
281                                 list[list_pos][5] = diagonals[i].prev_position;//last_y[tmpy-1];
282                                 diagonals[i].prev_position = list_pos;
283 //                              printf("A: %i %i\n",tmpx-1, tmpy-1);
284                                 list[list_pos][2] = M[toupper(seq1[tmpx-1])-'A'][toupper(seq2[tmpy-1])-'A'];
285                                 last_y[tmpy] = list_pos;
286 //                              tmpy_value = list_pos;
287 //                              tmpy_pos = tmpy;
288                                 --diagonals[i].current_pos[2];
289                                 ++diagonals[i].current_pos[0];
290                                 ++diagonals[i].current_pos[1];
291                                 ++list_pos;
292                         }
293                 }
294 //              last_y[tmpy_pos] = tmpy_value;
295
296                 
297                 //lower border
298                 if (list[list_pos-1][1] != l2)
299                 {
300                         list[list_pos] = vcalloc(6, sizeof(int));
301                         list[list_pos][0] = pos_x;
302                         list[list_pos][1] = l2;
303                         list[list_pos][3] = last_y[l2];
304                         
305                         list[list_pos][2] = -1000;
306                         list[list_pos][4] = list_pos-1;
307                         if (pos_x > l2)
308                                 list[list_pos][5] = last_x[pos_x-l2];
309                         else
310                                 list[list_pos][5] = l2-pos_x;
311                         last_y[l2] = list_pos;
312                         ++list_pos;
313                         
314                 }
315
316
317                 if ((b >= 0) && (diagonals[b].current_pos[1] > l2))
318                         --b;
319                 
320                 if ((a >0) && (diagonals[a-1].current_pos[0] - diagonals[a-1].current_pos[1] == pos_x))
321                 {
322                         --a;
323                         diagonals[a].prev_position = last_x[pos_x-1];
324                 }
325         }
326
327
328         dig_num = -1;
329         if (list_pos + l2+2 > current_size)
330         {
331                 current_size += list_pos + l2 + 2;
332                 list = vrealloc(list, current_size * sizeof(int*));
333         }
334
335
336 //      right border    
337         list[list_pos] = vcalloc(6, sizeof(int));
338         list[list_pos][0] = l1;
339         list[list_pos][1] = 0;
340         list[list_pos][3] = last_x[l1-1];
341         list[list_pos][2] = -1000;
342         ++list_pos;
343         
344         
345
346         for (i = 1; i <= l2; ++i)
347         {
348                 list[list_pos] = vcalloc(6, sizeof(int));
349                 list[list_pos][0] = l1;
350                 list[list_pos][1] = i;
351                 list[list_pos][3] = last_y[i];
352                 list[list_pos][4] = list_pos-1;
353                 y = last_y[i-1];
354                 if ((list[y][0] == l1-1) && (list[y][1] == i-1))
355                 {
356                         list[list_pos][5] = y;
357                         list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'];
358                 }
359                 else
360                 {
361                         if (i <= l1)
362                         {
363                                 list[list_pos][5] = last_x[l1-i];
364                         }
365                         else
366                         {
367                                 list[list_pos][5] = i-l1;
368                         }
369                         list[list_pos][2] = -1000;
370                 }
371                 ++list_pos;
372         }
373         
374         list[list_pos - l2][2] = -1000;
375
376         *num_points = list_pos;
377 //      vfree(diagx);
378 //      vfree(diagy);
379
380
381         return list;
382 }
383
384
385
386
387 int
388 seq_pair2blast_diagonal2(char *seq_file_name1,
389                                                  char *seq_file_name2,
390                                                  Diagonal **diagonals,
391                                                  int *dig_length,
392                                                  int l1,
393                                                  int l2,
394                                                  int is_dna)
395 {
396         char *out_file = vtmpnam(NULL);
397         char blast_command[200];
398         char blast_command2[200];
399
400         if (is_dna)
401         {
402                 sprintf(blast_command, "bl2seq -p blastn -i %s -j %s -D 1 -g F -o %s -S 1 -F F", seq_file_name1, seq_file_name2, out_file);
403         }
404         else
405         {
406                 sprintf(blast_command, "bl2seq -p blastp -i %s -j %s -D 1 -g F -o %s -F F -S 1", seq_file_name1, seq_file_name2, out_file);
407         }
408         system(blast_command);
409
410         Diagonal *diags = diagonals[0];
411         FILE *diag_f = fopen(out_file,"r");
412         char line[300];
413         fgets(line, 300, diag_f);
414         fgets(line, 300, diag_f);
415         fgets(line, 300, diag_f);
416         
417         
418         char delims[] = "\t";
419         int length, pos_q, pos_d;
420         int current_pos = 0;
421         while (fgets(line, 300, diag_f) != NULL)
422         {
423                 strtok(line, delims);
424                 strtok(NULL, delims);
425                 strtok(NULL, delims);
426                 length =  atoi(strtok(NULL, delims));
427                 strtok(NULL, delims);
428                 strtok(NULL, delims);
429                 pos_q = atoi(strtok(NULL, delims));
430                 strtok(NULL, delims);
431                 pos_d = atoi(strtok(NULL, delims));
432
433                 if (current_pos >= *dig_length)
434                 {
435                         (*dig_length) += 40;
436                         diags = vrealloc(diags, sizeof(Diagonal)*(*dig_length));
437                 }
438                 diags[current_pos].x = pos_q;
439                 diags[current_pos].y = pos_d;
440                 diags[current_pos++].length = length;
441         }
442         fclose(diag_f);
443         diagonals[0] = diags;
444         return current_pos;
445 }
446
447 /******************************COPYRIGHT NOTICE*******************************/
448 /*© Centro de Regulacio Genomica */
449 /*and */
450 /*Cedric Notredame */
451 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
452 /*All rights reserved.*/
453 /*This file is part of T-COFFEE.*/
454 /**/
455 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
456 /*    it under the terms of the GNU General Public License as published by*/
457 /*    the Free Software Foundation; either version 2 of the License, or*/
458 /*    (at your option) any later version.*/
459 /**/
460 /*    T-COFFEE is distributed in the hope that it will be useful,*/
461 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
462 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
463 /*    GNU General Public License for more details.*/
464 /**/
465 /*    You should have received a copy of the GNU General Public License*/
466 /*    along with Foobar; if not, write to the Free Software*/
467 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
468 /*...............................................                                                                                      |*/
469 /*  If you need some more information*/
470 /*  cedric.notredame@europe.com*/
471 /*...............................................                                                                                                                                     |*/
472 /**/
473 /**/
474 /*      */
475 /******************************COPYRIGHT NOTICE*******************************/