+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "fastal_lib_header.h"
-
-
-
-/*!
- * \file diagonal.c
- * \brief Source code for the calculation and preparation of diagonals
- */
-
-
-
-
-
-
-
-int
-diagonal_compare (const void * a, const void * b)
-{
-
- int tmp = ((((Diagonal*)a)->y - ((Diagonal*)a)->x) - (((Diagonal*)b)->y - ((Diagonal*)b)->x));
- if (tmp != 0)
- return tmp;
- return (((Diagonal*)a)->y - ((Diagonal*)b)->y);
-}
-
-int
-max(int a, int b)
-{
- if (a < b)
- return b;
- else
- return a;
-}
-
-
-
-Segment*
-extend_diagonals(Diagonal *diagonals, int *num_diagonal, int l1, int l2)
-{
-// sort diagonal by diagonal number
- int i;
-// printf("INPUT_PRE\n");
-// for (i = 0; i < num_diagonals; ++i)
-// {
-// printf("%i %i %i\n",diagonals[i].x, diagonals[i].y, diagonals[i].length);
-// }
- int num_diagonals = *num_diagonal;
- qsort (diagonals, num_diagonals, sizeof(Diagonal), diagonal_compare);
-
-
-// find nearest diagonal and expand
-// make shure that overlapping segments on the same diagonal are merged
-
- int diff;
-// printf("INPUT\n");
- for (i = 0; i < num_diagonals; ++i)
- {
- printf("%i %i %i\n",diagonals[i].x, diagonals[i].y, diagonals[i].length);
- }
- for (i = 0; i < num_diagonals; ++i)
- {
- diagonals[i].end_exp = 0;
- diagonals[i].front_exp = 0;
- }
- int first = 0;
- int next = 1;
- while (next < num_diagonals)
- {
-
- if (!((diagonals[next].y >= diagonals[first].y) && diagonals[next].x <= diagonals[first].x))
- {
- 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);
- diff = diagonals[next].y - (diagonals[first].y + diagonals[first].length);
- if (diff > 0)
- {
- if ((diagonals[first].y - diagonals[first].x) != ((diagonals[next].y - diagonals[next].x)))
- {
- diagonals[first].end_exp = max(diagonals[first].end_exp, diff);
- diagonals[next].front_exp = max(diagonals[next].front_exp, diff);
- // if (diagonals[next].x - diagonals[next].front_exp < 0)
-
- while (diagonals[++first].x == -1);
- }
- else
- {
- printf("ICH TUWAS\n");
- diagonals[first].length = diagonals[next].y + diagonals[next].length - diagonals[first].y;
- diagonals[first].end_exp = diagonals[next].end_exp;
- diagonals[next].x = -1;
- }
- }
- diff = diagonals[first].x - (diagonals[next].x + diagonals[next].length);
- if (diff > 0)
- {
- // diff = diagonals[next].x + diagonals[next].length - diagonals[first].x;
-
- diagonals[first].front_exp = max(diagonals[first].front_exp, diff);
- diagonals[next].end_exp = max(diagonals[next].end_exp, diff);
- // ++num_segments;
- // if (diagonals[first].x - diagonals[first].front_exp < 0)
-
- while (diagonals[++first].x == -1);
- }
- } else
- ++ first;
- ++next;
-
- }
-
- int num_segments =0;
- Segment *seg = vcalloc(num_diagonals, sizeof(Segment));
- int pos = 0;
-// int diag_num1, diag_num2;
- int *tmp1;
- Diagonal *tmp2;
- for (i = 0; i < num_diagonals; ++i)
- {
- if (diagonals[i].x != -1)
- {
- ++num_segments;
- tmp2 = &diagonals[i];
- seg[pos].segments = vcalloc(6, sizeof(int));
- seg[pos].current_pos = seg[pos].segments;
- tmp1 = seg[pos].segments;
- tmp1[0] = tmp2->x - tmp2->front_exp;
- tmp1[1] = tmp2->y - tmp2->front_exp;
- tmp1[2] = tmp2->front_exp + tmp2->length + tmp2->end_exp;
- tmp1[5] = 1;
- tmp1[4] = l2;
- tmp1[3] = tmp1[0] + (l2 - tmp1[1]);
- ++pos;
- }
- }
-
- for ( i = 0; i < num_segments; ++i )
- {
- 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]);
- }
-// exit(1);
- vrealloc(seg,num_segments*sizeof(Segment));
-
- *num_diagonal = num_segments;
-
-// vfree(diagonals);
- return seg;
-}
-
-
-int **
-segments2int(Segment *diagonals,
- int num_diagonals,
- char *seq1,
- char *seq2,
- Fastal_profile *profile1,
- Fastal_profile *profile2,
- int *num_points,
- Fastal_param *param_set)
-{
-// printf("NUM: %i\n", num_diagonals);
- int l1 = strlen(seq1);
- int l2 = strlen(seq2);
- int gep = param_set->gep;
-
- int dig_length;
- if (seq1 > seq2)
- dig_length = l1;
- else
- dig_length = l2;
-
- int current_size = num_diagonals*dig_length + l1 +l2;
-
- int **list = vcalloc(current_size, sizeof(int*));
-// int *diags = vcalloc(num_diagonals, sizeof(int));
- int i;
-// for (i = 0; i < num_diagonals; ++i)
-// {
-// diags[i] = l1 - diagonals[i*3] + diagonals[i*3+1];
-// }
-
-// qsort (diags, num_diagonals, sizeof(int), fastal_compare);
-
-
-// int *diagx = vcalloc(num_diagonals, sizeof(int));
-// int *diagy = vcalloc(num_diagonals, sizeof(int));
-
-
- //+1 because diagonals start here at position 1, like in "real" dynamic programming
- int a = -1, b = -1;
- for (i = 0; i < num_diagonals; ++i)
- {
- if (diagonals[i].current_pos[1] - diagonals[i].current_pos[0] < 0)
- {
-// diagx[i] = l1 - diags[i];
-// diagy[i] = 0;
- a= i;
- }
- else
- break;
- }
- ++a;
- b=a-1;
- for (; i < num_diagonals; ++i)
- {
-// diagx[i] = 0;
-// diagy[i] = diags[i]-l1;
- diagonals[i].prev_position = diagonals[i].current_pos[1] - diagonals[i].current_pos[0];
- b = i;
- }
-
-// vfree(diags);
- int tmpy_pos;
- int tmpy_value;
- int **M = param_set->M;
- int *last_y = vcalloc(l2+1, sizeof(int));
- int *last_x = vcalloc(l1+1, sizeof(int));
- last_y[0] = 0;
-
- last_x[0] = 0;
- list[0] = vcalloc(6, sizeof(int));
-
- int list_pos = 1;
- int dig_num = l1;
- int tmp_l2 = l2 + 1;
-
- //left border
- for (; list_pos < tmp_l2; ++list_pos)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
- list[list_pos][0] = 0;
- list[list_pos][1] = list_pos;
- last_y[list_pos] = list_pos;
- list[list_pos][2] = list_pos*gep;
- list[list_pos][4] = list_pos-1;
- }
-
- int pos_x = 0;
- int y;
- int tmp_l1 = l1-1;
- int tmpx, tmpy;
- while (pos_x < tmp_l1)
- {
- if (list_pos + num_diagonals+2 > current_size)
- {
- current_size += num_diagonals*1000;
- list = vrealloc(list, current_size * sizeof(int*));
- }
- //upper border
- list[list_pos] = vcalloc(6, sizeof(int));
- list[list_pos][0] = ++pos_x;
- list[list_pos][1] = 0;
- list[list_pos][2] = pos_x * gep;
- list[list_pos][3] = last_y[0];
-// tmpy_value = list_pos;
-// tmpy_pos = 0;
- last_x[pos_x] = list_pos;
- ++list_pos;
-
- //diagonals
- for (i = a; i <= b; ++i)
- {
-
- if (pos_x == diagonals[i].current_pos[0])
- {
- if (diagonals[i].current_pos[2] == 0)
- {
- diagonals[i].current_pos+=3;
- }
- list[list_pos] = vcalloc(6, sizeof(int));
- tmpx = diagonals[i].current_pos[0];
- tmpy = diagonals[i].current_pos[1];
- list[list_pos][0] = tmpx;
- list[list_pos][1] = tmpy;
- list[list_pos][3] = last_y[tmpy];
- list[list_pos][4] = list_pos-1;
- list[list_pos][5] = diagonals[i].prev_position;//last_y[tmpy-1];
- diagonals[i].prev_position = list_pos;
-// printf("A: %i %i\n",tmpx-1, tmpy-1);
- list[list_pos][2] = M[toupper(seq1[tmpx-1])-'A'][toupper(seq2[tmpy-1])-'A'];
- last_y[tmpy] = list_pos;
-// tmpy_value = list_pos;
-// tmpy_pos = tmpy;
- --diagonals[i].current_pos[2];
- ++diagonals[i].current_pos[0];
- ++diagonals[i].current_pos[1];
- ++list_pos;
- }
- }
-// last_y[tmpy_pos] = tmpy_value;
-
-
- //lower border
- if (list[list_pos-1][1] != l2)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
- list[list_pos][0] = pos_x;
- list[list_pos][1] = l2;
- list[list_pos][3] = last_y[l2];
-
- list[list_pos][2] = -1000;
- list[list_pos][4] = list_pos-1;
- if (pos_x > l2)
- list[list_pos][5] = last_x[pos_x-l2];
- else
- list[list_pos][5] = l2-pos_x;
- last_y[l2] = list_pos;
- ++list_pos;
-
- }
-
-
- if ((b >= 0) && (diagonals[b].current_pos[1] > l2))
- --b;
-
- if ((a >0) && (diagonals[a-1].current_pos[0] - diagonals[a-1].current_pos[1] == pos_x))
- {
- --a;
- diagonals[a].prev_position = last_x[pos_x-1];
- }
- }
-
-
- dig_num = -1;
- if (list_pos + l2+2 > current_size)
- {
- current_size += list_pos + l2 + 2;
- list = vrealloc(list, current_size * sizeof(int*));
- }
-
-
-// right border
- list[list_pos] = vcalloc(6, sizeof(int));
- list[list_pos][0] = l1;
- list[list_pos][1] = 0;
- list[list_pos][3] = last_x[l1-1];
- list[list_pos][2] = -1000;
- ++list_pos;
-
-
-
- for (i = 1; i <= l2; ++i)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
- list[list_pos][0] = l1;
- list[list_pos][1] = i;
- list[list_pos][3] = last_y[i];
- list[list_pos][4] = list_pos-1;
- y = last_y[i-1];
- if ((list[y][0] == l1-1) && (list[y][1] == i-1))
- {
- list[list_pos][5] = y;
- list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'];
- }
- else
- {
- if (i <= l1)
- {
- list[list_pos][5] = last_x[l1-i];
- }
- else
- {
- list[list_pos][5] = i-l1;
- }
- list[list_pos][2] = -1000;
- }
- ++list_pos;
- }
-
- list[list_pos - l2][2] = -1000;
-
- *num_points = list_pos;
-// vfree(diagx);
-// vfree(diagy);
-
-
- return list;
-}
-
-
-
-
-int
-seq_pair2blast_diagonal2(char *seq_file_name1,
- char *seq_file_name2,
- Diagonal **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna)
-{
- char *out_file = vtmpnam(NULL);
- char blast_command[200];
- char blast_command2[200];
-
- if (is_dna)
- {
- 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);
- }
- else
- {
- 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);
- }
- system(blast_command);
-
- Diagonal *diags = diagonals[0];
- FILE *diag_f = fopen(out_file,"r");
- char line[300];
- fgets(line, 300, diag_f);
- fgets(line, 300, diag_f);
- fgets(line, 300, diag_f);
-
-
- char delims[] = "\t";
- int length, pos_q, pos_d;
- int current_pos = 0;
- while (fgets(line, 300, diag_f) != NULL)
- {
- strtok(line, delims);
- strtok(NULL, delims);
- strtok(NULL, delims);
- length = atoi(strtok(NULL, delims));
- strtok(NULL, delims);
- strtok(NULL, delims);
- pos_q = atoi(strtok(NULL, delims));
- strtok(NULL, delims);
- pos_d = atoi(strtok(NULL, delims));
-
- if (current_pos >= *dig_length)
- {
- (*dig_length) += 40;
- diags = vrealloc(diags, sizeof(Diagonal)*(*dig_length));
- }
- diags[current_pos].x = pos_q;
- diags[current_pos].y = pos_d;
- diags[current_pos++].length = length;
- }
- fclose(diag_f);
- diagonals[0] = diags;
- return current_pos;
-}
-
-/******************************COPYRIGHT NOTICE*******************************/
-/*© Centro de Regulacio Genomica */
-/*and */
-/*Cedric Notredame */
-/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
-/*All rights reserved.*/
-/*This file is part of T-COFFEE.*/
-/**/
-/* T-COFFEE is free software; you can redistribute it and/or modify*/
-/* it under the terms of the GNU General Public License as published by*/
-/* the Free Software Foundation; either version 2 of the License, or*/
-/* (at your option) any later version.*/
-/**/
-/* T-COFFEE is distributed in the hope that it will be useful,*/
-/* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
-/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
-/* GNU General Public License for more details.*/
-/**/
-/* You should have received a copy of the GNU General Public License*/
-/* along with Foobar; if not, write to the Free Software*/
-/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
-/*............................................... |*/
-/* If you need some more information*/
-/* cedric.notredame@europe.com*/
-/*............................................... |*/
-/**/
-/**/
-/* */
-/******************************COPYRIGHT NOTICE*******************************/
git://source.jalview.org / jabaws.git / blobdiff