+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <ctype.h>
-#include <string.h>
-// #include <unistd.h>
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "define_header.h"
-#include "dp_lib_header.h"
-#include "fastal_lib_header.h"
-#include "fast_tree_header.h"
-
-
-// #include <omp.h>
-// #define CHUNKSIZE 100
-// #define N 1000
-
-
-//Fastal_param *param_set;
-
-
-/*! \mainpage T-Coffee Index Page
- *
- * \section intro_sec Introduction
- *
- * This is the introduction.
- *
- * \section install_sec Installation
- *
- * \subsection step1 Step 1: Opening the box
- *
- * etc...
- * \section fastal_sec Fastal
- *
- * This program is a very fast aligner. It is capable of aligning huge sets of sequences because it keeps as much as necessary on hard disk.
- */
-
-
-
-
-/*!
- * \file fastal.c
- * \brief Source code for the fastal algorithm
- */
-
-
-
-//************************** sparse dynamic aligning **********************************************************
-
-
-void
-fill_arguments_sparse(Sparse_dynamic_param* method_arguments_p)
-{
- method_arguments_p->diagonals = vcalloc(3,sizeof(Diagonal));
- method_arguments_p->dig_length = vcalloc(1,sizeof(int));
- *method_arguments_p->dig_length = 3;
- method_arguments_p->list = NULL;
- method_arguments_p->list_length = vcalloc(1,sizeof(int));
- *method_arguments_p->list_length = 0;
- method_arguments_p->file_name1 = vtmpnam(NULL);
- method_arguments_p->file_name2 = vtmpnam(NULL);
-}
-
-void
-free_sparse(Sparse_dynamic_param* method_arguments_p)
-{
- vfree(method_arguments_p->diagonals);
- vfree(method_arguments_p->dig_length);
- vfree(method_arguments_p->list_length);
-}
-
-
-/**
- * \brief One run of sparse dynamic programming.
- *
- * \param profiles The profiles.
- * \param param_set The fastal parameters.
- * \param method_arguments_p The method arguments.
- * \param is_dna Sequences are DNA (\a is_dna = 1) or protein.
- * \param edit_file The edit file.
- * \param prof_file the profile file.
- * \param number Number of the parent node.
- * \return The length of the alignment.
- */
-int
-sparse_dyn(Fastal_profile **profiles,
- Fastal_param *param_set,
- void *method_arguments_p,
- int is_dna,
- FILE *edit_file,
- FILE *prof_file,
- int number)
-{
-// printf("WHAT THE HELL ARE YOU DOING HERE?\n");
- Sparse_dynamic_param *arguments = (Sparse_dynamic_param*)method_arguments_p;
-// static char *file_name1 = vtmpnam(NULL);
-// static char *file_name2 = vtmpnam(NULL);
- char *file_name1 = arguments->file_name1;
- char *file_name2 = arguments->file_name2;
- char *seq1, *seq2;
- Fastal_profile *tmp1 = profiles[0];
- Fastal_profile *tmp2 = profiles[1];
-
- seq1 = profile2consensus(tmp1, param_set);
- seq2 = profile2consensus(tmp2, param_set);
-
-
- int **diagonals_p = &(arguments->diagonals);
- int num_diagonals = -1;
- if (!strcmp(param_set->diag_method, "blastz"))
- {
- FILE *cons_f = fopen(file_name1,"w");
- fprintf(cons_f, ">%i\n", tmp1->prf_number);
- fprintf(cons_f, "%s", seq1);
- fprintf( cons_f, "\n");
- fclose(cons_f);
- cons_f = fopen(file_name2,"w");
- fprintf(cons_f, ">%i\n", tmp2->prf_number);
- fprintf(cons_f, "%s", seq2);
- fprintf( cons_f, "\n");
- fclose(cons_f);
- num_diagonals = seq_pair2blastz_diagonal(file_name1, file_name2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna);
- }
- else if (!strcmp(param_set->diag_method, "blast"))
- {
- FILE *cons_f = fopen(file_name1,"w");
- fprintf(cons_f, ">%i\n", tmp1->prf_number);
- fprintf(cons_f, "%s", seq1);
- fprintf( cons_f, "\n");
- fclose(cons_f);
- cons_f = fopen(file_name2,"w");
- fprintf(cons_f, ">%i\n", tmp2->prf_number);
- fprintf(cons_f, "%s", seq2);
- fprintf( cons_f, "\n");
- fclose(cons_f);
- int l1 = strlen(seq1);
- int l2 = strlen(seq2);
- num_diagonals = seq_pair2blast_diagonal(file_name1, file_name2, diagonals_p, arguments->dig_length, l1, l2, is_dna);
-// int *num_p = &num_diagonals;
-// Segment* seg = extend_diagonals(*diagonals_p, num_p, l1, l2);
-// printf("A: %i\n", num_diagonals);
- }
- else if (!strcmp(param_set->diag_method, "ktup"))
- {
- num_diagonals = seq_pair2diagonal_own(seq1, seq2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna, 3);
-// num_diagonals = seq_pair2diagonal_swift(seq1, seq2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna, 3);
- }
-
-
-// arguments->diagonals = diagonals_p[0];
-// arguments->list = segments2int(seg, *num_p, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
-
-// t ** segments2int_gap(Segment *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set);
-
-// arguments->list = diagonals2int_dot(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
-// arguments->list = diagonals2int_euclidf(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
- arguments->list = diagonals2int_gap_test(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
-// arguments->list = diagonals2int(arguments->diagonals, num_diagonals, seq1, seq2, arguments->list_length, param_set);
- int alignment_length = list2linked_pair_wise_fastal(profiles[0], profiles[1], param_set, arguments->list, *arguments->list_length, edit_file, prof_file, number);
- int x;
-
- for (x = 0; x < *arguments->list_length; ++x)
- {
- vfree(arguments->list[x]);
- }
- vfree(arguments->list);
- arguments->list = NULL;
- vfree(seq1);
- vfree(seq2);
- return alignment_length;
-}
-
-
-int
-fastal_compare (const void * a, const void * b)
-{
- return (*(int*)a - *(int*)b);
-}
-
-
-/**
- * \brief Makes a sorted list out of diagonals.
- *
- * \param diagonals A list of diagonals to use during dynamic programming.
- * \param num_diagonals Number of diagonals.
- * \param seq1 Sequence 1.
- * \param seq2 Sequence 2.
- * \param num_points Number of points in the list
- * \param param_set Fastal parameters.
- * \return A 2-dim array which contains all points needed for the sparse dynamic programming algorithm.
- */
-int **
-diagonals2int(int *diagonals,
- int num_diagonals,
- char *seq1,
- char *seq2,
- int *num_points,
- Fastal_param *param_set)
-{
-
- 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 (diags[i] < l1)
- {
- 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;
- 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;
- 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)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
-
- list[list_pos][0] = ++diagx[i];
-
- list[list_pos][1] = ++diagy[i];
- list[list_pos][3] = last_y[diagy[i]];
- list[list_pos][4] = list_pos-1;
- list[list_pos][5] = last_y[diagy[i]-1];
- list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'];
- last_y[tmpy_pos] = tmpy_value;
- tmpy_value = list_pos;
- tmpy_pos = diagy[i];
-
- ++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) && (diagy[b] == l2))
- --b;
-
- if ((a >0) && (diagx[a-1] == pos_x))
- --a;
- }
-
-
- 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;
-}
-
-
-/**
- * \brief Makes a sorted list out of diagonals.
- *
- * \param diagonals A list of diagonals to use during dynamic programming.
- * \param num_diagonals Number of diagonals.
- * \param seq1 Sequence 1.
- * \param seq2 Sequence 2.
- * \param num_points Number of points in the list
- * \param param_set Fastal parameters.
- * \return A 2-dim array which contains all points needed for the sparse dynamic programming algorithm.
- */
-int **
-diagonals2int_gap_test(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- int l1 = strlen(seq1);
- int l2 = strlen(seq2);
- int gep = param_set->gep;
-
- int current_size = l2+l1;
-
- 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));
- int *old_pos = 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 (diags[i] < l1)
- {
- 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;
- 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 diags_old = l2;
-
-// int tmp = l1;
- int y;
- int tmp_l1 = l1-1;
- 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)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
-
- list[list_pos][0] = ++diagx[i];
-
- list[list_pos][1] = ++diagy[i];
- list[list_pos][3] = last_y[diagy[i]];
- list[list_pos][4] = list_pos-1;
- list[list_pos][5] = last_y[diagy[i]-1];
-
-
-
-
-//SIMPLEGAP
- int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
- double gap_num = 0;
- int char_c;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
-
- gap_num += profile1->prf[char_c][diagx[i]-1] + profile2->prf[char_c][diagy[i]-1];
- }
-
- gap_num /= num_seq;
-
- list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'] * gap_num;
-
-// CLUSTAL
-// int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
-// double gap_num = 0;
-// int char_c, char_c2;
-// for (char_c = 0; char_c < alphabet_size; ++char_c)
-// for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
-// {
-// gap_num += (profile1->prf[char_c][diagx[i]-1]/profile1->number_of_sequences) * (profile2->prf[char_c][diagy[i]-1]/profile2->number_of_sequences) * M[param_set->pos2char[char_c]-'A'][param_set->pos2char[char_c2]-'A'];
-// }
-// list[list_pos][2] = gap_num;
-
-
- last_y[tmpy_pos] = tmpy_value;
- tmpy_value = list_pos;
- tmpy_pos = diagy[i];
-
- ++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) && (diagy[b] == l2))
- --b;
-
- if ((a >0) && (diagx[a-1] == pos_x))
- --a;
- }
-
-
- 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;
- int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
- double gap_num = 0;
- int char_c;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- gap_num += profile1->prf[char_c][l1-1] + profile2->prf[char_c][i-1];
- }
-
- gap_num /= num_seq;
-
- list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'] * gap_num;
- }
- 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);
- vfree(old_pos);
-
- return list;
-}
-
-
-int **
-diagonals2int_euclidf(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- int l1 = strlen(seq1);
- int l2 = strlen(seq2);
- int gep = param_set->gep;
-
- int current_size = l2+l1;
-
- 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));
- int *old_pos = 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 (diags[i] < l1)
- {
- 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;
- 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 diags_old = l2;
-
-// int tmp = l1;
- int y;
- int tmp_l1 = l1-1;
- 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)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
-
- list[list_pos][0] = ++diagx[i];
-
- list[list_pos][1] = ++diagy[i];
- list[list_pos][3] = last_y[diagy[i]];
- list[list_pos][4] = list_pos-1;
- list[list_pos][5] = last_y[diagy[i]-1];
- int char_c;
- double tmp_score = 0;
- double freq1, freq2;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- freq1 = (double)profile1->prf[char_c][diagx[i]-1] / profile1->number_of_sequences;
-
- freq2 = (double)profile2->prf[char_c][diagy[i]-1] / profile2->number_of_sequences;
-
- tmp_score += ( freq1 - freq2) * (freq1 - freq2);
- }
-
- list[list_pos][2] = 10 - sqrt(tmp_score);
-
- last_y[tmpy_pos] = tmpy_value;
- tmpy_value = list_pos;
- tmpy_pos = diagy[i];
-
- ++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) && (diagy[b] == l2))
- --b;
-
- if ((a >0) && (diagx[a-1] == pos_x))
- --a;
- }
-
-
- 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;
- int char_c;
- int tmp_score = 0;
- double freq1, freq2;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- freq1 = profile1->prf[char_c][l1-1] / profile1->number_of_sequences;
- freq2 = profile2->prf[char_c][i-1] / profile2->number_of_sequences;
- tmp_score += ( freq1 - freq2) * (freq2 - freq1);
- }
- list[list_pos][2] = 10 - sqrt(tmp_score);
-// 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);
- vfree(old_pos);
-
- return list;
-}
-
-int **
-diagonals2int_dot(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- int l1 = strlen(seq1);
- int l2 = strlen(seq2);
- int gep = param_set->gep;
-
- int current_size = l2+l1;
-
- 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));
- int *old_pos = 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 (diags[i] < l1)
- {
- 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;
- 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 diags_old = l2;
-
-// int tmp = l1;
- int y;
- int tmp_l1 = l1-1;
- 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)
- {
- list[list_pos] = vcalloc(6, sizeof(int));
-
- list[list_pos][0] = ++diagx[i];
-
- list[list_pos][1] = ++diagy[i];
- list[list_pos][3] = last_y[diagy[i]];
- list[list_pos][4] = list_pos-1;
- list[list_pos][5] = last_y[diagy[i]-1];
- int char_c;
- double tmp_score = 0;
- double freq1, freq2;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- freq1 = (double)profile1->prf[char_c][diagx[i]-1] / profile1->number_of_sequences;
-
- freq2 = (double)profile2->prf[char_c][diagy[i]-1] / profile2->number_of_sequences;
-
- tmp_score += freq1 * freq2;
- }
-
- list[list_pos][2] = tmp_score * 10;
-
- last_y[tmpy_pos] = tmpy_value;
- tmpy_value = list_pos;
- tmpy_pos = diagy[i];
-
- ++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) && (diagy[b] == l2))
- --b;
-
- if ((a >0) && (diagx[a-1] == pos_x))
- --a;
- }
-
-
- 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;
- int char_c;
- int tmp_score = 0;
- double freq1, freq2;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- freq1 = profile1->prf[char_c][l1-1] / profile1->number_of_sequences;
- freq2 = profile2->prf[char_c][i-1] / profile2->number_of_sequences;
- tmp_score += freq2 * freq1;
- }
- list[list_pos][2] = tmp_score * 10;
-// 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);
- vfree(old_pos);
-
- return list;
-}
-
-
-void
-combine_profiles2file(int **prf1,
- int **prf2,
- int pos1,
- int pos2,
- Fastal_param *param_set,
- FILE *prof_f,
- char state)
-{
- int alphabet_size = param_set->alphabet_size;
- char *pos2aa = &(param_set->pos2char[0]);
- int i;
- int x = 0;
- if (state == 'M')
- {
- for (i = 0; i < alphabet_size; ++i)
- if (prf1[i][pos1] + prf2[i][pos2] > 0)
- {
- if (x)
- fprintf(prof_f," %c%i", pos2aa[i],prf1[i][pos1]+prf2[i][pos2]);
- else
- fprintf(prof_f,"%c%i", pos2aa[i],prf1[i][pos1]+prf2[i][pos2]);
- x = 1;
- }
- fprintf(prof_f,"\n");
- }
- else if (state == 'D')
- {
- for (i = 0; i < alphabet_size; ++i)
- if (prf2[i][pos2] > 0)
- {
- if (x)
- fprintf(prof_f," %c%i", pos2aa[i],prf2[i][pos2]);
- else
- fprintf(prof_f,"%c%i", pos2aa[i],prf2[i][pos2]);
- x = 1;
- }
- fprintf(prof_f,"\n");
- }
- else
- {
- for (i = 0; i < alphabet_size; ++i)
- if (prf1[i][pos1] > 0)
- {
- if (x)
- fprintf(prof_f," %c%i", pos2aa[i],prf1[i][pos1]);
- else
- fprintf(prof_f,"%c%i", pos2aa[i],prf1[i][pos1]);
- x = 1;
- }
- fprintf(prof_f,"\n");
- }
-}
-
-
-
-#define LIN(a,b,c) a[b*5+c]
-/**
- * Calculates a fast and sparse dynamic programming matrix
- *
- * \param prf1 Profile of first sequence.
- * \param prf2 Profile of second sequence.
- * \param param_set The parameter for the alignment.
- * \param list The list of diagonals.
- * \param n number of dots.
- * \param edit_f File to save the edit information.
- * \param prof_f File to save the profile.
- * \param node_number Number of the new profile.
- */
-int
-list2linked_pair_wise_fastal(Fastal_profile *prf1,
- Fastal_profile *prf2,
- Fastal_param *param_set,
- int **list,
- int n,
- FILE *edit_f,
- FILE *prof_f,
- int node_number)
-{
- int a,b, i, j, LEN=0, start_trace = -1;
- int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
-// static int **slist;
- static long *MI, *MJ, *MM,*MT2;
-// static int *sortseq;
- static int max_size;
- int gop, gep, igop, igep;
- int l1, l2, l, ls;
- char **al;
- int ni=0, nj=0;
- long score;
- int nomatch = param_set->nomatch;
-
- l1=prf1->length;
- l2=prf2->length;
-
- al=declare_char (2,l1+l2+1);
-
-
-
- igop=param_set->gop;
- gep=igep=param_set->gep;
- if (n>max_size)
- {
- max_size=n;
-
- vfree (MI);vfree (MJ); vfree (MM);
-
- MI=vcalloc (5*n, sizeof (long));
- MJ=vcalloc (5*n, sizeof (long));
- MM=vcalloc (5*n, sizeof (long));
-
- }
- else
- {
- for (a=0; a<n; a++)
- for (b=0; b<5; b++)
- LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
- }
-
- for (a=0; a<n; a++)
- {
- i=list[a][0];
- j=list[a][1];
-
-
- if (i==l1 || j==l2)gop=0;
- else gop=igop;
-
- if (i==l1 && j==l2)start_trace=a;
- else if ( i==0 || j==0)
- {
- LIN(MM,a,0)=-1000000;
- if (j==0)
- {
- LIN(MJ,a,0)=-10000000;
- LIN(MI,a,0)=gep*i;
- }
- else if (i==0)
- {
- LIN(MI,a,0)=-10000000;
- LIN(MJ,a,0)=gep*j;
- }
-
- LIN(MI,a,1)=LIN(MJ,a,1)=-1;
- LIN(MI,a,2)=LIN(MJ,a,2)=i;
- LIN(MI,a,3)=LIN(MJ,a,3)=j;
- continue;
- }
-
- pi = list[a][3];
- pj = list[a][4];
- ij = list[a][5];
-
- prev_i=list[pi][0];
- prev_j=list[pj][1];
-
- delta_i=list[a][0]-list[pi][0];
- delta_j=list[a][1]-list[pj][1];
-
- /*Linear Notation*/
- LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
- LIN(MI,a,1)=pi;
- LIN(MI,a,2)=delta_i;
- LIN(MI,a,3)=0;
- LIN(MI,a,4)=(LIN(MI,pi,0) >=(LIN(MM,pi,0)+gop))?'i':'m';
-
- LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
- LIN(MJ,a,1)=pj;
- LIN(MJ,a,2)=0;
- LIN(MJ,a,3)=delta_j;
-
- LIN(MJ,a,4)=(LIN(MJ,pj,0) >=LIN(MM,pj,0)+gop)?'j':'m';
-
- if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
- {
- LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*nomatch);
-
- LIN(MM,a,1)=ij;
- LIN(MM,a,2)=ls;
- LIN(MM,a,3)=ls;
- if ( LIN(MM,ij,0) >=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
- else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
- else LIN(MM,a,4)='j';
-
- }
- else
- {
- LIN(MM,a,0)=UNDEFINED;
- LIN(MM,a,1)=-1;
- }
- }
-
- a=start_trace;
- if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
- else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
- else MT2=MJ;
-
- score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
-
- i=l1;
- j=l2;
-
- while (!(i==0 &&j==0))
- {
- int next_a;
- l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
- // HERE ("%c from %c %d %d SCORE=%d [%d %d] [%2d %2d]", T2[a][5],T2[a][4], T2[a][2], T2[a][3], T2[a][0], gop, gep, i, j);
- if (i==0)
- {
- while ( j>0)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--; LEN++;
- }
- }
- else if (j==0)
- {
- while ( i>0)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--; LEN++;
- }
- }
-
-// else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
- else
- {
- for (b=0; b<l; b++, LEN++)
- {
- if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
- else al[0][LEN]=0;
-
- if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
- else al[1][LEN]=0;
- }
-
- next_a=LIN(MT2,a,1);
- if (LIN(MT2,a,4)=='m')MT2=MM;
- else if (LIN(MT2,a,4)=='i')MT2=MI;
- else if (LIN(MT2,a,4)=='j')MT2=MJ;
- a=next_a;
- }
- }
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], LEN);
-
- fprintf(edit_f, "%i\n%i\n%i\n%i\n",prf1->prf_number, prf2->prf_number, prf1->is_leaf, prf2->is_leaf);
- fprintf(prof_f, "%i\n0\n%i\n1\n%i\n", node_number,LEN, prf1->number_of_sequences+prf2->number_of_sequences);
-
- char statec[] = {'M','D','I'};
- int num = 0;
- int state = 0;
- i = 0;
- j = 0;
-
- for ( b=0; b< LEN; b++)
- {
- if ((al[0][b]==1) && (al[1][b]==1))
- {
-
- combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'M');
- ++i;
- ++j;
- if (state != 0)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 0;
- }
- else
- ++num;
- }
- else if (al[0][b]==1)
- {
-// prf1->prf[param_set->alphabet_size-1] += prf2->num_sequences;
- combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'I');
- ++i;
- if (state != 2)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 2;
- }
- else
- ++num;
- }
- else if (al[1][b]==1)
- {
-// prf2->prf[param_set->alphabet_size-1] += prf1->num_sequences;
- combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'D');
- ++j;
- if (state != 1)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 1;
- }
- else
- ++num;
- }
- }
-
-
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 1;
-
-
- fprintf(edit_f,"*\n");
- fprintf(prof_f,"*\n");
- free_char (al, -1);
-// exit(0);
- return LEN;
-}
-
-
-
-
-
-
-/**
- * \brief Turns a profile into a consensus sequence.
- *
- * The character with the highest number of occurences is used as consensus. Gaps are not included. For example: 10 '-' and one 'A' would give 'A' as consensus.
- * \param profile The profile.
- * \param file_name Name of the file to save the consensus sequence in.
- * \param param_set The parameter of the fastal algorithm.
- * \return the sequence
- */
-char*
-profile2consensus(Fastal_profile *profile, Fastal_param *param_set)
-{
-
-// FILE *cons_f = fopen(file_name,"w");
-// fprintf(cons_f, ">%i\n", profile->prf_number);
- char* seq = vcalloc(profile->length+1, sizeof(char));
- int i, j;
- int most_pos = -1, most;
- int alphabet_size = param_set->alphabet_size;
- int **prf = profile->prf;
- char *pos2char = param_set->pos2char;
- for (i = 0; i < profile->length; ++i)
- {
- most = -1;
- for (j = 0; j < alphabet_size; ++j)
- {
- if (prf[j][i] > most)
- {
- most = prf[j][i];
- most_pos = j;
- }
- }
- seq[i] = pos2char[most_pos];
-// fprintf(cons_f, "%c",pos2char[most_pos]);
- }
- return seq;
-}
-
-
-int
-diag_compare (const void * a, const void * b)
-{
- return (((Diagonal_counter*)b)->count - ((Diagonal_counter*)a)->count);
-}
-
-/**
- * \brief Calculates the diagonals between two sequences.
- *
- * Uses to calculate the diagonals.
- * \param seq_file1 File with sequence 1.
- * \param seq_file2 File with sequence 2.
- * \param diagonals An array where the diagonal points will be stored.
- * \param dig_length length of \a diagonals .
- * \param num_points Number of points in all diagonals.
- * \return number of diagonals;
- */
-int
-seq_pair2diagonal_own(char *seq1,
- char *seq2,
- int **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna,
- int word_length)
-{
-// word_length = 7;
- int word_number, i;
- int ng;
- if (is_dna)
- {
- word_number = (int)pow(5, word_length);
- ng = 4;
- }
- else
- {
- word_number = (int)pow(24, word_length);
- ng = 24;
- }
- int **word_index = vcalloc(word_number, sizeof(int*));
- for (i = 0 ; i < word_number; ++i)
- {
- word_index[i] = vcalloc(20, sizeof(int));
- word_index[i][0] = 2;
- word_index[i][1] = 20;
- }
-
-
- //making of k-tup index of seq1
-
- int *prod=vcalloc (word_length, sizeof(int));
- for ( i=0; i<word_length; i++)
- {
- prod[word_length-i-1]=(int)pow(ng,i);
- }
-
- int aa[256];
- if (is_dna)
- {
- aa['A'] = 0;
- aa['C'] = 1;
- aa['G'] = 2;
- aa['T'] = 3;
- aa['U'] = 3;
- }
- else
- {
- aa['A'] = 0;
- aa['B'] = 20;
- aa['C'] = 1;
- aa['D'] = 2;
- aa['E'] = 3;
- aa['F'] = 4;
- aa['G'] = 5;
- aa['H'] = 6;
- aa['I'] = 7;
- aa['J'] = 20;
- aa['K'] = 8;
- aa['L'] = 9;
- aa['M'] = 10;
- aa['N'] = 11;
- aa['P'] = 12;
- aa['Q'] = 13;
- aa['R'] = 14;
- aa['S'] = 15;
- aa['T'] = 16;
- aa['V'] = 17;
- aa['W'] = 18;
- aa['X'] = 20;
- aa['Y'] = 19;
- aa['X'] = 20;
- }
- int index = 0;
- for (i = 0; i < word_length; ++i)
- {
- index += aa[(short)seq1[i]] *prod[i];
- }
- word_index[index][2] = 0;
- word_index[index][0] = 3;
- int z = -1;
- int *tmp;
- for (; i < l1; ++i)
- {
- index -= aa[(short)seq1[++z]] * prod[0];
- index *= ng;
- index += aa[(short)seq1[i]];
- tmp = word_index[index];
- if (tmp[0] == tmp[1])
- {
- tmp[1] += 25;
- tmp = vrealloc(tmp, word_index[index][1] *sizeof(int));
- word_index[index] = tmp;
- }
- tmp[tmp[0]++] = i;
- }
-
-
-
- //counting diagonals
- const int window_length = 14;
-
- Diagonal_counter *diag_index = vcalloc(l1+l2, sizeof(Diagonal_counter));
- int num = l1+l2;
- for (i = 0; i < num; ++i)
- {
- diag_index[i].diagonal = i;
- diag_index[i].count = 0;
- }
- index = 0;
-
- int j;
- for (i = 0; i < word_length; ++i)
- {
- index += aa[(short)seq2[i]] *prod[i];
- for (j = 2; j < word_index[index][0]; ++j)
- {
- ++(diag_index[i - word_index[index][j] + l1].count);
- }
- }
-
- z = -1;
- int i2 = i-1;
- int second_index = index;
- for (; i < window_length; ++i)
- {
- index -= aa[(short)seq2[++z]] * prod[0];
- index *= ng;
- index += aa[(short)seq2[i]];
- tmp = word_index[index];
- for (j = 2; j < tmp[0]; ++j)
- {
- ++(diag_index[i - tmp[j] + l1].count);
- }
- }
- int z2 = -1;
- for (; i < l2; ++i)
- {
- index -= aa[(short)seq2[++z]] * prod[0];
- index *= ng;
- index += aa[(short)seq2[i]];
- second_index -= aa[(short)seq2[++z2]] * prod[0];
- second_index *= ng;
- second_index += aa[(short)seq2[++i2]];
-
- tmp = word_index[index];
- for (j = 2; j < tmp[0]; ++j)
- {
- ++(diag_index[i - tmp[j] + l1].count);
- }
-
-
- tmp = word_index[second_index];
- for (j = 2; j < tmp[0]; ++j)
- {
- if (diag_index[i2 - tmp[j] + l1].count > window_length-3)
- diag_index[i2 - tmp[j] + l1].count = window_length+100;
- else
- --diag_index[i2 - tmp[j] + l1].count;
- }
- }
-
-
- //choose diagonals
- int *diags = diagonals[0];
- int current_pos = 0;
-
-
- qsort (diag_index, num, sizeof(Diagonal_counter*), diag_compare);
-
- i = 0;
- int y, x;
- while (diag_index[i].count > window_length+10)
- {
- if (current_pos > (*dig_length)-3)
- {
- (*dig_length) += 30;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
-
-
- y = diag_index[i].diagonal - l1;
- if (y < 0)
- {
- x = y * (-1);
- y = 0;
- }
- else
- {
- x = 0;
- }
- diags[current_pos++] = x;
- diags[current_pos++] = y;
- diags[current_pos++] = 200;
- ++i;
- }
-
- vfree(diag_index);
- for (i = 0; i < word_number; ++i)
- vfree(word_index[i]);
- vfree(word_index);
- diagonals[0] = diags;
- return current_pos/3;
-}
-
-
-
-int
-seq_pair2diagonal_swift(char *seq1,
- char *seq2,
- int **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna,
- int word_length)
-{
- int word_number, i;
- int ng;
- if (is_dna)
- {
- word_number = (int)pow(5, word_length);
- ng = 5;
- }
- else
- {
- word_number = (int)pow(24, word_length);
- ng = 24;
- }
- int **word_index = vcalloc(word_number, sizeof(int*));
- for (i = 0 ; i < word_number; ++i)
- {
- word_index[i] = vcalloc(20, sizeof(int));
- word_index[i][0] = 2;
- word_index[i][1] = 20;
- }
-
-
- //making of k-tup index of seq1
-
- int *prod=vcalloc (word_length, sizeof(int));
- for ( i=0; i<word_length; i++)
- {
- prod[word_length-i-1]=(int)pow(ng,i);
- }
-
- int aa[256];
- aa['A'] = 0;
- aa['C'] = 1;
- aa['G'] = 2;
- aa['T'] = 3;
- aa['U'] = 4;
- int index = 0;
- for (i = 0; i < word_length; ++i)
- {
- index += aa[(short)seq1[i]] *prod[i];
- }
- word_index[index][2] = 0;
- word_index[index][0] = 3;
- int z = -1;
- int *tmp;
- for (; i < l1; ++i)
- {
- index -= aa[(short)seq1[++z]] * prod[0];
- index *= ng;
- index += aa[(short)seq1[i]];
- tmp = word_index[index];
- if (tmp[0] == tmp[1])
- {
- tmp[1] += 25;
- tmp = vrealloc(tmp, word_index[index][1] *sizeof(int));
- word_index[index] = tmp;
- }
- tmp[tmp[0]++] = i;
- }
-
-
- //counting diagonals
- const int window_length = 14;
- const int threshold = 12;
-
- Swift_diagonal *diag_index = vcalloc(l1+l2, sizeof(Swift_diagonal));
- int num = l1+l2;
- for (i = 0; i < num; ++i)
- {
- diag_index[i].diagonal = i;
- diag_index[i].count = 0;
- diag_index[i].start = -99999;
- diag_index[i].end = -99999;
- }
-
- index = 0;
-
- int j;
- for (i = 0; i < word_length; ++i)
- {
- index += aa[(short)seq2[i]] *prod[i];
- for (j = 2; j < word_index[index][0]; ++j)
- {
- ++(diag_index[i - word_index[index][j] + l1].count);
- }
- }
-
- z = -1;
- int tmp_index;
- for (; i < l2; ++i)
- {
- index -= aa[(short)seq2[++z]] * prod[0];
- index *= ng;
- index += aa[(short)seq2[i]];
- tmp = word_index[index];
- for (j = 2; j < tmp[0]; ++j)
- {
- tmp_index = i - tmp[j] + l1;
- if (i - diag_index[tmp_index].end > window_length)
- {
- if (diag_index[tmp_index].count < threshold)
- {
- diag_index[tmp_index].count = 0;
- diag_index[tmp_index].start = i;
- diag_index[tmp_index].end = i + word_length;
- }
-
- }
- else
- {
- ++(diag_index[tmp_index].count);
- }
- }
-
- }
-
-
-
- // choose diagonals
- int *diags = diagonals[0];
- int current_pos = 0;
- int x, y;
- for (i = 0; i < num; ++i)
- {
- if (diag_index[i].count > threshold)
- {
- if (current_pos > (*dig_length)-3)
- {
- (*dig_length) += 30;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- y = diag_index[i].diagonal - l1;
- if (y < 0)
- {
- x = y * (-1);
- y = 0;
- }
- else
- {
- x = 0;
- }
- diags[current_pos++] = x;
- diags[current_pos++] = y;
- diags[current_pos++] = 200;
- }
- }
-
- vfree(diag_index);
- for (i = 0; i < word_number; ++i)
- vfree(word_index[i]);
- vfree(word_index);
- diagonals[0] = diags;
-
- return current_pos/3;
-}
-
-
-
-/**
- * \brief Calculates the diagonals between two sequences.
- *
- * Uses a k-tup index to choose diagonals.
- * \param seq_file1 File with sequence 1.
- * \param seq_file2 File with sequence 2.
- * \param diagonals An array where the diagonal points will be stored.
- * \param dig_length length of \a diagonals .
- * \param num_points Number of points in all diagonals.
- * \return number of diagonals;
- */
-int
-seq_pair2blast_diagonal(char *seq_file_name1,
- char *seq_file_name2,
- int **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna)
-{
-// static int blast_measure[12]={0,0,0,0,0,0,0,0,0,0,0,0};
- int *diag = vcalloc(l1 + l2, sizeof(int));
- char *out_file = vtmpnam(NULL);
- char blast_command[200];
-// char blast_command2[200];
-// if (x)
-// {
-// int i;
-// printf("BLAST-Types:\n");
-// for (i = 0; i < 11; ++i)
-// {
-// printf("Type %i: %i\n", i, blast_measure[i]);
-// }
-// return 0;
-// }
-// char blast_command2[600];
-// sprintf(blast_command2, "less %s", out_file);
-
- 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);
-
- int *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))-1;
- strtok(NULL, delims);
- pos_d = atoi(strtok(NULL, delims))-1;
-
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- if (diag[l1-(pos_q)+pos_d] == 0)
- {
- diag[l1-(pos_q)+pos_d] =1;
- diags[current_pos++] = pos_q;
- diags[current_pos++] = pos_d;
- diags[current_pos++] = length;
- }
-
- }
- fclose(diag_f);
- int round = 0;
- int e_threshold = 10;
- while ((current_pos == 0) && (round < 10))
- {
- if (is_dna)
- {
- sprintf(blast_command, "bl2seq -p blastn -i %s -j %s -D 1 -g F -o %s -S 1 -F F -W 6 -e %i", seq_file_name1, seq_file_name2, out_file, e_threshold);
- }
- else
- {
- sprintf(blast_command, "bl2seq -p blastp -i %s -j %s -D 1 -g F -o %s -F F -S 1 -e %i", seq_file_name1, seq_file_name2, out_file, e_threshold);
- }
- system(blast_command);
- e_threshold *= 10;
- 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";
- 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))-1;
- strtok(NULL, delims);
- pos_d = atoi(strtok(NULL, delims))-1;
-
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- if (diag[l1-(pos_q)+pos_d] == 0)
- {
- diag[l1-(pos_q)+pos_d] =1;
- diags[current_pos++] = pos_q;
- diags[current_pos++] = pos_d;
- diags[current_pos++] = length;
- }
- }
- fclose(diag_f);
- ++round;
- if (current_pos < 27)
- current_pos = 0;
- }
-// ++blast_measure[round];
-
- if (current_pos == 0)
- {
- printf("BLAST NOT SUCCESFULL\n");
- if (l1 < l2)
- {
- int i;
- int diff = l2 - l1 + 10;
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = 10;
-// printf("A: %i\n", diff);
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
- }
- else
- {
- int i;
- int diff = 10;
-// printf("A: %i\n", diff);
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = l1 - l2 + 10;
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
-
- }
- }
-
-
- vfree(diag);
-
- diagonals[0] = diags;
- return current_pos/3;
-}
-
-
-
-int
-seq_pair2blat_diagonal(char *seq_file_name1,
- char *seq_file_name2,
- int **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna)
-{
- int *diag = vcalloc(l1 + l2, sizeof(int));
- char *out_file = vtmpnam(NULL);
- char blast_command[200];
-// char blast_command2[200];
-// char blast_command2[600];
-// sprintf(blast_command2, "less %s", out_file);
-
- if (is_dna)
- {
- sprintf(blast_command, "blat %s %s %s -out=blast8 -q=dna -t=dna -maxGap=0 >/dev/null 2>/dev/null", seq_file_name2, seq_file_name1, out_file);
- }
- else
- {
- sprintf(blast_command, "blat %s %s %s -out=blast8 -prot -maxGap=0 >/dev/null 2>/dev/null", seq_file_name2, seq_file_name1, out_file);
- }
- system(blast_command);
-
- int *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))-1;
- strtok(NULL, delims);
- pos_d = atoi(strtok(NULL, delims))-1;
-
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- if (diag[l1-(pos_q)+pos_d] == 0)
- {
- diag[l1-(pos_q)+pos_d] =1;
- diags[current_pos++] = pos_q;
- diags[current_pos++] = pos_d;
- diags[current_pos++] = length;
- }
- }
- if (current_pos == 0)
- {
- printf("BLAT NOT SUCCESFULL\n");
- if (l1 < l2)
- {
- int i;
- int diff = l2 - l1 + 10;
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = 10;
-// printf("A: %i\n", diff);
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
- }
- else
- {
- int i;
- int diff = 10;
-// printf("A: %i\n", diff);
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = l1 - l2 + 10;
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
-
- }
- }
-
-// printf("END\n");
- vfree(diag);
- fclose(diag_f);
- diagonals[0] = diags;
- return current_pos/3;
-}
-
-
-
-/**
- * \brief Calculates the diagonals between two sequences.
- *
- * Uses blastz to calculate the diagonals.
- * \param seq_file1 File with sequence 1.
- * \param seq_file2 File with sequence 2.
- * \param diagonals An array where the diagonal points will be stored.
- * \param dig_length length of \a diagonals .
- * \param num_points Number of points in all diagonals.
- * \return number of diagonals;
- */
-int
-seq_pair2blastz_diagonal(char *seq_file_name1,
- char *seq_file_name2,
- int **diagonals,
- int *dig_length,
- int l1,
- int l2,
- int is_dna)
-{
- int *diag = vcalloc(l1 + l2, sizeof(int));
- char *out_file = vtmpnam(NULL);
- char blast_command[200];
-// char blast_command2[200];
-// char blast_command2[600];
-// sprintf(blast_command2, "less %s", out_file);
-
- if (is_dna)
- {
- sprintf(blast_command, "~/Download/blastz-source/blastz %s %s B=0 K=10000> %s", seq_file_name1, seq_file_name2, out_file);
- }
- else
- {
- printf("SORRY - no BLASTZ with amino accid\n");
- exit(0);
- }
- system(blast_command);
-
- int *diags = diagonals[0];
- FILE *diag_f = fopen(out_file,"r");
- char line[300];
- char delims[] = " ";
-// char *result = NULL;
- int length, pos_q, pos_d;
- int current_pos = 0;
- while (fgets(line, 300, diag_f) != NULL)
- {
- if (line[0] == 'a')
- {
- char *line_tmp;
- while (fgets(line, 300, diag_f) != NULL)
- {
- if (line[0] == '}')
- break;
-
- if (line[2] == 'l')
- {
- line_tmp = &line[4];
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- pos_q = atoi(strtok(line_tmp, delims));
- pos_d = atoi(strtok(NULL, delims));
- length = atoi(strtok(NULL, delims) - pos_q);
- if (diag[l1-(pos_q)+pos_d] == 0)
- {
- diag[l1-(pos_q)+pos_d] =1;
- diags[current_pos++] = pos_q;
- diags[current_pos++] = pos_d;
- diags[current_pos++] = length;
- }
- }
- }
- }
- }
-
-
- if (current_pos == 0)
- {
- printf("BLASTZ NOT SUCCESFULL\n");
- if (l1 < l2)
- {
- int i;
- int diff = l2 - l1 + 10;
-
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = 10;
-
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
- }
- else
- {
- int i;
- int diff = 10;
-
- for (i = diff; i > 0; --i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = 0;
- diags[current_pos++] = i;
- diags[current_pos++] = 100;
- }
- diff = l1 - l2 + 10;
- for (i = 0; i < diff; ++i)
- {
- if (current_pos >= *dig_length-20)
- {
- (*dig_length) += 90;
- diags = vrealloc(diags, sizeof(int)*(*dig_length));
- }
- diags[current_pos++] = i;
- diags[current_pos++] = 0;
- diags[current_pos++] = 100;
- }
-
- }
- }
-
- vfree(diag);
- fclose(diag_f);
- diagonals[0] = diags;
- return current_pos/3;
-}
-
-
-
-//************************** needleman-wunsch aligning **********************************************************
-
-
-void
-fill_arguments_nw(Nw_param* method_arguments_p, int alphabet_size)
-{
- method_arguments_p-> dyn_matrix = vcalloc(1,sizeof(double*));
- method_arguments_p->dyn_matrix[0] = vcalloc(1,sizeof(double));
- method_arguments_p->length1 = vcalloc(1,sizeof(int));
- method_arguments_p->length2 = vcalloc(1,sizeof(int));
- *method_arguments_p->length1 = 1;
- *method_arguments_p->length2 = 1;
- method_arguments_p->sumup_prf = vcalloc(alphabet_size+1,sizeof(int*));
- int i;
- for (i = 0; i < alphabet_size+1; ++i)
- method_arguments_p->sumup_prf[i] = vcalloc(1,sizeof(int));
- method_arguments_p->sumup_length = vcalloc(1,sizeof(int));
- *method_arguments_p->sumup_length = 1;
-}
-
-
-void
-free_nw(Nw_param* method_arguments_p, int alphabet_size)
-{
- free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->dyn_matrix);
- int i;
- for (i = 0; i <= alphabet_size; ++i)
- {
- vfree(method_arguments_p->sumup_prf[i]);
- }
- vfree(method_arguments_p->sumup_prf);
- vfree(method_arguments_p->length1);
- vfree(method_arguments_p->length2);
- vfree(method_arguments_p->sumup_length);
-}
-
-
-/**
- * \brief One run of needleman-wunsch dynamic programming.
- *
- * \param profiles The profiles.
- * \param param_set The fastal parameters.
- * \param method_arguments_p The method arguments.
- * \param is_dna Sequences are DNA (\a is_dna = 1) or protein.
- * \param edit_file The edit file.
- * \param prof_file the profile file.
- * \param number Number of the parent node.
- * \return The length of the alignment.
- */
-int
-nw_dyn(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number)
-{
- Nw_param *arguments = (Nw_param*)method_arguments_p;
-// int old_length1 = *arguments->length1;
-// int old_length2 = *arguments->length2;
- arguments->dyn_matrix = resize_dyn_matrix(arguments->dyn_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
- *arguments->length1 = profiles[0]->length+1;
- *arguments->length2 = profiles[1]->length+1;
- int alignment_length = prf_nw(profiles[0], profiles[1], arguments->dyn_matrix, edit_file, arguments->sumup_prf, arguments->sumup_length, param_set);
- write2file(arguments->sumup_prf, alignment_length, prof_file, number,profiles[0]->number_of_sequences + profiles[1]->number_of_sequences, param_set);
- return alignment_length;
-}
-
-
-/**
- * \brief This method takes a profile and turns it into a sumed up version.
- *
- * Required for NW-algorithm.
- * \param profile The profile to sum up.
- * \param sumup A field where the result will be stored.
- * \param param_set Parameters for the fastal algorithm.
- * \return The new \a sumup.
- */
-int**
-sumup_profile(Fastal_profile *profile,
- int **sumup,
- Fastal_param *param_set)
-{
-
- char *pos2aa = &(param_set->pos2char[0]);
- int alphabet_size = param_set->alphabet_size;
- int **M = param_set->M;
- int prof_length = profile->length;
-
- int i,j,k;
-
- for (i = 0; i < prof_length; ++i)
- {
- sumup[alphabet_size][i] = 0;
- for (k = 0; k < alphabet_size; ++k)
- {
- sumup[k][i] = 0;
- sumup[alphabet_size][i] += profile->prf[k][i];
- for (j = 0; j < alphabet_size; ++j)
- {
- sumup[k][i] += profile->weight * profile->prf[j][i] * M[pos2aa[j]-'A'][pos2aa[k]-'A'];
- }
- }
- }
-
- return sumup;
-}
-
-
-/**
- * \brief Turns the dynamic programming matrix into a editfile and calculates the new profile.
- *
- * Required for NW-algorithm.
- * \param prog_matrix The dynamic programming matrix.
- * \param prf1 The first profile.
- * \param prf2 The second profile.
- * \param edit_f A File object (already opened) to write the edit sequence into.
- * \param prf_field A 2-dim array to save the new profile into.
- * \param field_length Length of the new profile.
- * \param param_set Parameters of the Fastal-Algorithm
- */
-int
-nw_matrix2edit_file(double **prog_matrix, //dynamic programming matrix
- Fastal_profile *prf1, //profile of dim1
- Fastal_profile *prf2, //profile of dim2
- FILE *edit_f, //file to safe the edit in
- int **prf_field, //space to safe the new profile
- int *field_length,
- Fastal_param *param_set) //length of prf_field
-{
-// int **M = param_set->M;
- int alphabet_size = param_set->alphabet_size;
- double gap_cost = param_set -> gop;
- fprintf(edit_f, "%i\n%i\n%i\n%i\n",prf1->prf_number, prf2->prf_number, prf1->is_leaf, prf2->is_leaf);
- int sum[] = {0,0,0};
- char sumc[] = {'M','I','D'};
- int last = 0;
- int n = 0;
- int m = 0;
- int field_pos = 0;
- int i;
- int prf1_length = prf1->length;
- int prf2_length = prf2->length;
- while ((n < prf1_length) && (m < prf2_length))
- {
- //if necesarry allocate more memory for result
- if ((*field_length)-alphabet_size < field_pos)
- {
- (*field_length) += ENLARGEMENT_PER_STEP;
-
- for (i = 0; i <alphabet_size+1; ++i)
- {
- prf_field[i] = vrealloc(prf_field[i], (*field_length)*sizeof(int));
- }
- }
-
- if (prog_matrix[n][m] == (prog_matrix[n+1][m] +gap_cost))
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n];
- }
- ++n;
- ++ field_pos;
-
- if (last != 1)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 1;
- ++sum[last];
- }
- else if (prog_matrix[n][m] == (prog_matrix[n][m+1] +gap_cost))
- {
-
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf2->prf[i][m];
- }
- ++m;
- ++ field_pos;
- if (last != 2)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 2;
- ++sum[last];
- }
- else
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n] + prf2->prf[i][m];
- }
- ++n;
- ++m;
- ++ field_pos;
- if (last != 0)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 0;
- ++sum[last];
- }
- }
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
-
- //gaps in prf2
- last = 0;
- while (n < prf1_length)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n];
- }
- ++n;
- ++ field_pos;
- ++last;
- }
- if (last > 0)
- fprintf(edit_f,"I%i\n",last);
-
- //gaps in prf1
- last = 0;
- while (m < prf2_length)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf2->prf[i][m];
- }
- ++m;
- ++ field_pos;
- ++last;
- }
- if (last > 0)
- fprintf(edit_f,"D%i\n",last);
- fprintf(edit_f,"*\n");
- return field_pos;
-}
-
-
-
-/**
- * \brief Pairwise alignments of profile is done here.
- *
- * \param profile1 Profile of sequence 1
- * \param profile2 Profile of sequence 2
- * \param prog_matrix Matrix for dynamic programming
- * \param edit_file_name The edit_file_name
- * \param sumup_prf The sumup version of profile 1, which later contains the aligned profile.
- * \param sumup_length Contains the length of the aligned profile.
- * \return length of the aligned profile
- */
-int
-prf_nw(Fastal_profile *profile1,
- Fastal_profile *profile2,
- double **prog_matrix,
- FILE *edit_file_name,
- int **sumup_prf,
- int *sumup_length,
- Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- double gap_cost = param_set->gop;
-
- int i;
- if (*sumup_length < profile1->length)
- {
- for (i = 0; i < alphabet_size+1; ++i)
- {
- sumup_prf[i] = vrealloc(sumup_prf[i], profile1->length*sizeof(int));
- }
- *sumup_length = profile1->length;
- }
- sumup_prf = sumup_profile(profile1, sumup_prf, param_set);
-
-
-
- int j,k;
- int prof1_length = profile1->length;
- int prof2_length = profile2->length;
-
-// int** M = param_set->M;
- double match_score;
-// int amino_counter;
- int residue_pairs = 0;
-
- for (i = prof2_length; i > 0; --i)
- {
- prog_matrix[prof1_length][i] = gap_cost * (prof2_length-i);
- }
-
- i = prof1_length-1;
- prog_matrix[prof1_length][prof2_length] = 0.0;
- while (i >=0)
- {
- j = prof2_length-1;
-
- prog_matrix[i][prof2_length] = gap_cost*(prof1_length-i);
- while (j >=0)
- {
- match_score = 0.0;
- residue_pairs = 0;
- for (k = 0; k < alphabet_size; ++k)
- {
- residue_pairs += profile2->prf[k][j];
- match_score += (profile2->prf[k][j] * sumup_prf[k][i]);
- }
- match_score /= (residue_pairs * sumup_prf[alphabet_size][i]);
- prog_matrix[i][j] = MAX3(prog_matrix[i+1][j+1]+match_score, prog_matrix[i+1][j]+gap_cost, prog_matrix[i][j+1]+gap_cost);
-
- --j;
- }
- --i;
- }
- return nw_matrix2edit_file(prog_matrix, profile1, profile2, edit_file_name, sumup_prf, sumup_length, param_set);
-}
-
-
-/************** GOTOH ***********************/
-
-
-void
-fill_arguments_gotoh(Gotoh_param* method_arguments_p, int alphabet_size)
-{
- method_arguments_p->m_matrix = vcalloc(1,sizeof(double*));
- method_arguments_p->m_matrix[0] = vcalloc(1,sizeof(double));
- method_arguments_p->d_matrix = vcalloc(1,sizeof(double*));
- method_arguments_p->d_matrix[0] = vcalloc(1,sizeof(double));
- method_arguments_p->i_matrix = vcalloc(1,sizeof(double*));
- method_arguments_p->i_matrix[0] = vcalloc(1,sizeof(double));
- method_arguments_p->length1 = vcalloc(1,sizeof(int));
- method_arguments_p->length2 = vcalloc(1,sizeof(int));
- method_arguments_p->log_saver = vcalloc(alphabet_size+1, sizeof(double*));
- *method_arguments_p->length1 = 1;
- *method_arguments_p->length2 = 1;
- method_arguments_p->sumup_prf = vcalloc(alphabet_size+1,sizeof(int*));
- int i;
- for (i = 0; i < alphabet_size+1; ++i)
- {
- method_arguments_p->sumup_prf[i] = vcalloc(1,sizeof(int));
- method_arguments_p->log_saver[i] = vcalloc(1, sizeof(double));
- }
- method_arguments_p->sumup_length = vcalloc(1,sizeof(int));
- *method_arguments_p->sumup_length = 1;
-}
-
-
-void
-free_gotoh(Gotoh_param* method_arguments_p, int alphabet_size)
-{
- free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->m_matrix);
- free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->i_matrix);
- free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->d_matrix);
-
- int i;
- for (i = 0; i <= alphabet_size; ++i)
- {
- vfree(method_arguments_p->sumup_prf[i]);
- }
- vfree(method_arguments_p->sumup_prf);
- vfree(method_arguments_p->length1);
- vfree(method_arguments_p->length2);
- vfree(method_arguments_p->sumup_length);
-}
-
-
-int
-gotoh_dyn(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number)
-{
- Gotoh_param *arguments = (Gotoh_param*)method_arguments_p;
- arguments->m_matrix = resize_dyn_matrix(arguments->m_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
- arguments->i_matrix = resize_dyn_matrix(arguments->i_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
- arguments->d_matrix = resize_dyn_matrix(arguments->d_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
- int i;
- if (profiles[1]->length > *arguments->length2-1)
- {
- for (i = 0; i < param_set->alphabet_size; ++i)
- {
- arguments->log_saver[i] = vrealloc(arguments->log_saver[i], (profiles[1]->length)*sizeof(double*));
- }
- }
- *arguments->length1 = profiles[0]->length+1;
- *arguments->length2 = profiles[1]->length+1;
- int alignment_length = prf_gotoh(profiles[0], profiles[1], edit_file, arguments, param_set);
- write2file(arguments->sumup_prf, alignment_length, prof_file, number, profiles[0]->number_of_sequences + profiles[1]->number_of_sequences, param_set);
- return alignment_length;
-}
-
-
-int
-gotoh_matrix2edit_file(double **m_matrix, //dynamic programming matrix
- double **v_matrix, //dynamic programming matrix
- double **h_matrix, //dynamic programming matrix
- Fastal_profile *prf1, //profile of dim1
- Fastal_profile *prf2, //profile of dim2
- FILE *edit_f, //file to safe the edit in
- int **prf_field, //space to safe the new profile
- int *field_length,
- Fastal_param *param_set) //length of prf_field
-{
- double comp_num = log((double)prf1->number_of_sequences) + log((double)prf2->number_of_sequences);
- int** M = param_set->M;
- int alphabet_size = param_set->alphabet_size;
- double gep = param_set -> gep;
- fprintf(edit_f, "%i\n%i\n%i\n%i\n",prf1->prf_number, prf2->prf_number, prf1->is_leaf, prf2->is_leaf);
- int sum[] = {0,0,0};
- char sumc[] = {'M','I','D'};
- int last = 0;
- int n = 0;
- int m = 0;
- int field_pos = 0;
- int i;
- int prf1_length = prf1->length;
- int prf2_length = prf2->length;
- int current_mode = 0;
- //determine start mode
- char *pos2char = param_set->pos2char;
- if (h_matrix[n][m] == m_matrix[n][m])
- {
- current_mode = 2;
- }
- else
- {
-
- if (v_matrix[n][m] == m_matrix[n][m])
- {
- current_mode = 1;
- }
- else
- {
- current_mode = 0;
- }
- }
-// printf("%f %f %f - %i\n",h_matrix[n][m],v_matrix[n][m],m_matrix[n][m], current_mode);
- while ((n < prf1_length) && (m < prf2_length))
- {
- //if necesarry allocate more memory for result
- if ((*field_length)-alphabet_size < field_pos)
- {
- (*field_length) += ENLARGEMENT_PER_STEP;
-
- for (i = 0; i <alphabet_size+1; ++i)
- {
- prf_field[i] = vrealloc(prf_field[i], (*field_length)*sizeof(int));
- }
- }
-
-
- if (current_mode == 2)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf2->prf[i][m];
- }
- if (h_matrix[n][m] != (h_matrix[n][m+1]+gep))
- {
- current_mode = 0;
- }
- ++m;
- ++ field_pos;
- if (last != 2)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 2;
- ++sum[last];
- }
- else
- {
- if (current_mode== 1)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n];
- }
- if (v_matrix[n][m] != (v_matrix[n+1][m]+gep))
- {
- current_mode = 0;
- }
- ++n;
- ++ field_pos;
-
- if (last != 1)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 1;
- ++sum[last];
- }
- else
- {
- double match_score = 0.0;
- int char_c, char_c2;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
- {
-
- if ((log(prf1->prf[char_c][n]) != -1) && ( log(prf2->prf[char_c2][m]) != -1))
- {
- match_score += exp(log((double)prf1->prf[char_c][n]) + log((double)prf2->prf[char_c2][m])-comp_num) * M[pos2char[char_c]-'A'][pos2char[char_c2]-'A'];
- }
- }
- }
- if (m_matrix[n+1][m+1] + match_score != m_matrix[n][m])
- {
- if (m_matrix[n][m] == v_matrix[n][m])
- {
- current_mode = 1;
- continue;
- }
- if (m_matrix[n][m] == h_matrix[n][m])
- {
- current_mode = 2;
- continue;
- }
- }
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n] + prf2->prf[i][m];
- }
- ++n;
- ++m;
- ++ field_pos;
- if (last != 0)
- {
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
- sum[last] = 0;
- }
- last = 0;
- ++sum[last];
- }
- }
-
- }
- fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
-
- int needed = MAX(prf1_length -n, prf2_length -m);
-
- if ((*field_length) - needed -10 < field_pos)
- {
- (*field_length) += needed +10;
-
- for (i = 0; i <alphabet_size+1; ++i)
- {
- prf_field[i] = vrealloc(prf_field[i], (*field_length)*sizeof(int));
- }
- }
- //gaps in prf2
- last = 0;
- while (n < prf1_length)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf1->prf[i][n];
- }
- ++n;
- ++ field_pos;
- ++last;
- }
- if (last > 0)
- fprintf(edit_f,"I%i\n",last);
-
- //gaps in prf1
- last = 0;
- while (m < prf2_length)
- {
- for (i = 0; i<alphabet_size; ++i)
- {
- prf_field[i][field_pos] = prf2->prf[i][m];
- }
- ++m;
- ++ field_pos;
- ++last;
- }
- if (last > 0)
- fprintf(edit_f,"D%i\n",last);
-
- fprintf(edit_f,"*\n");
- return field_pos;
-}
-
-
-
-
-/**
- * \brief The gotoh dynamic programming algorithm.
- *
- * \param profile1 The first profile.
- */
-int
-prf_gotoh(Fastal_profile *profile1,
- Fastal_profile *profile2,
- FILE *edit_file_name,
- Gotoh_param *arguments,
- Fastal_param *param_set)
-{
-
-// printf("I AM HERE - again\n");
- int **sumup_prf = arguments->sumup_prf;
- int *sumup_length = arguments->sumup_length;
- int alphabet_size = param_set->alphabet_size;
- double gop = param_set->gop;
- double gep = param_set->gep;
-
- const int INF = -999999;
- int i;
-
- double **m_matrix = arguments->m_matrix;
- double **h_matrix = arguments->i_matrix;
- double **v_matrix = arguments->d_matrix;
-
- int j;
- int prof1_length = profile1->length;
- int prof2_length = profile2->length;
-
- int** M = param_set->M;
- double match_score;
- for (i = prof2_length; i >= 0; --i)
- {
- m_matrix[prof1_length][i] = gop + gep * (prof2_length-i);
- v_matrix[prof1_length][i] = INF;
- h_matrix[prof1_length][i] = m_matrix[prof1_length][i];
- }
-
- m_matrix[prof1_length][prof2_length] = 0.0;
- h_matrix[prof1_length][prof2_length] = INF;
- v_matrix[prof1_length][prof2_length] = INF;
- int l;
- double comp_num = log((double)profile1->number_of_sequences) + log((double)profile2->number_of_sequences);
- static double *log_test = NULL;
- if (!log_test)
- log_test = vcalloc(alphabet_size, sizeof(double));
-// int k;
- int **prf1 = profile1->prf;
- int **prf2 = profile2->prf;
- double **log_test2 = arguments->log_saver;
- for (l = 0; l < alphabet_size; ++l)
- {
- for (i = 0; i < profile2->length; ++i)
- {
- if (prf2[l][i] > 0)
- {
- log_test2[l][i] = log(prf2[l][i]);
- }
- else
- log_test2[l][i] = -1;
- }
- }
-
- char *pos2char = param_set->pos2char;
- i = prof1_length-1;
- while (i >=0)
- {
- j = prof2_length-1;
-
- for (l = 0; l < alphabet_size; ++l)
- {
- if (prf1[l][i] > 0)
- log_test[l] = log((double)prf1[l][i]);
- else
- log_test[l] = -1;
- }
- m_matrix[i][prof2_length] = gop + gep *(prof1_length-i);
- v_matrix[i][prof2_length] = m_matrix[i][prof2_length];
- h_matrix[i][prof2_length] = INF;
- while (j >=0)
- {
-
- match_score = 0.0;
- v_matrix[i][j] = (MAX(v_matrix[i+1][j], m_matrix[i+1][j] + gop) + gep);
- h_matrix[i][j] = (MAX(h_matrix[i][j+1], m_matrix[i][j+1] + gop) + gep);
-
- int char_c, char_c2;
- int num = 0;
- for (char_c = 0; char_c < alphabet_size; ++char_c)
- {
- for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
- {
-
- if ((log_test[char_c] != -1) && (log_test2[char_c2][j] != -1))
- {
- match_score += exp(log_test[char_c] + log_test2[char_c2][j]-comp_num) * M[pos2char[char_c]-'A'][pos2char[char_c2]-'A'];
- }
- }
- }
-
- m_matrix[i][j] = m_matrix[i+1][j+1]+match_score;
-
- if (m_matrix[i][j] < v_matrix[i][j])
- {
- m_matrix[i][j] = v_matrix[i][j];
- }
- if (m_matrix[i][j] < h_matrix[i][j])
- {
- m_matrix[i][j] = h_matrix[i][j];
- }
-
- --j;
- }
- --i;
- }
- return gotoh_matrix2edit_file(m_matrix, v_matrix, h_matrix, profile1, profile2, edit_file_name, sumup_prf, sumup_length, param_set);
-}
-
-
-/************* OTHER STUFF ******************/
-
-/**
- * \brief Writes the alignment into the profile file and the edit file.
- *
- * \param profiles The two profiles to combine.
- * \param alignment The alinment information.
- * \param alignment The length of the alignment.
- * \param edit_f The edit file.
- * \param prof_f The profile file.
- * \param node_number the new node number.
- */
-void
-alignment2files(Fastal_profile **profiles,
- Fastal_param *param_set,
- int **alignment,
- int alignment_length,
- FILE *edit_f,
- FILE *prof_f,
- int node_number)
-{
- fprintf(edit_f, "%i\n%i\n%i\n%i\n",profiles[0]->prf_number, profiles[1]->prf_number, profiles[0]->is_leaf, profiles[1]->is_leaf);
- fprintf(prof_f, "%i\n0\n%i\n1\n", node_number, alignment_length);
-
- int **prf1 = profiles[0]->prf;
- int **prf2 = profiles[1]->prf;
- int i = 0;
- int pos = 0;
- int pos1, pos2;
-
- char statec[] = {'M','D','I'};
- int num = 0;
- int state = 0;
-
- while (i < alignment_length)
- {
-
- pos1 = alignment[0][pos];
- pos2 = alignment[1][pos];
- // match
- if ((pos1 != -1) && (pos2 != -1))
- {
-
- combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'M');
- if (state != 0)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 0;
- }
- else
- ++num;
- ++i;
- }
- // insertion in seq 1
- else if (pos1 != -1)
- {
- combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'I');
- if (state != 2)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 2;
- }
- else
- ++num;
- ++i;
- }
- // deletion in seq 1
- else if (pos2 != -1)
- {
- combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'D');
- if (state != 1)
- {
- fprintf(edit_f, "%c%i\n",statec[state], num);
- num =1;
- state = 1;
- }
- else
- ++num;
- ++i;
- }
- ++pos;
- }
- fprintf(edit_f, "%c%i\n",statec[state], num);
-
- fprintf(edit_f,"*\n");
- fprintf(prof_f,"*\n");
-
-}
-
-
-//******************************* OTHER STUFF ***********************
-
-/**
- * \brief Reads the sequence from a given position in a fasta file and turns it into a profile.
- *
- * \param seq_file The file where the sequence is stored.
- * \param off_set The off_set from the beginning of the file to the position of the sequence name.
- * \param profile The profile where the sequence will be stored into.
- * \param prf_number The number of this profile.
- */
-void
-file_pos2profile(FILE *seq_file, //File with sequences
- long off_set, //offset of sequence from the beginning of file point to the sequence name, not to the sequence itself
- Fastal_profile *profile, //profile to save into
- int prf_number, //number of the profile
- Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- profile->is_leaf = 1;
- profile->number_of_sequences = 1;
- int *aa2pos = &(param_set->char2pos[0]);
- const int LINE_LENGTH = 500;
- char line[LINE_LENGTH];
- profile->num_sequences = 1;
- profile->prf_number = prf_number;
- fseek (seq_file , off_set , SEEK_SET );
-
- fgets (line, LINE_LENGTH , seq_file);
- int seq_length = 0;
- int i, j, x;
-
- while(fgets(line, LINE_LENGTH, seq_file)!=NULL)
- {
- if (line[0] != '>')
- {
- line[LINE_LENGTH-1] = '\n';
- if (seq_length + LINE_LENGTH >= profile->allocated_memory)
- {
- for (i = 0; i < alphabet_size; ++i)
- {
- profile->prf[i] = vrealloc(profile->prf[i], (profile->allocated_memory+PROFILE_ENLARGEMENT)*sizeof(int));
- }
- profile->allocated_memory += PROFILE_ENLARGEMENT;
- }
-
- i = 0;
- x = 0;
- while ((line[i] != '\n') && (line[i] != '\0'))
- {
- if (line[i] != '-')
- {
- for(j = 0; j<alphabet_size; ++j )
- profile->prf[j][seq_length+x] = 0;
- profile->prf[aa2pos[toupper((short)line[i])]][seq_length+x] = 1;
- ++x;
- }
- ++i;
- }
- seq_length += x;
-
- }
- else
- break;
- }
- profile->length = seq_length;
-
-}
-
-
-
-/**
- * \brief Constructs index of fasta_file.
- *
- * The index is of length n (n= number of sequences in the given multi fasta file.). In the order of appearance in the file the position of each sequence in the file is stored.
- * \param file_name The file with the sequences.
- * \param file_positions Array to save the positions in.
- * \return The number of sequences in \a file_name.
- */
-int
-make_index_of_file(char *file_name, //file with sequences
- long **file_positions) //array to save the positions
-{
- const int LINE_LENGTH = 150;
- (*file_positions) = vcalloc(ENLARGEMENT_PER_STEP, sizeof(long));
-
- FILE *file = fopen(file_name,"r");
-
- char line[LINE_LENGTH];
-
- int num_of_sequences = 0;
- int mem_for_pos = ENLARGEMENT_PER_STEP;
-
-
- if (file == NULL)
- {
- printf("FILE NOT FOUND\n");
- exit(1);
- }
- else
- {
- (*file_positions)[num_of_sequences] = ftell(file);
- while(fgets(line, LINE_LENGTH , file)!=NULL)
- {
-// int length = strlen(line);
- if (line[0] == '>')
- {
- ++num_of_sequences;
-
- if (num_of_sequences == mem_for_pos)
- {
- (*file_positions) = vrealloc((*file_positions),(ENLARGEMENT_PER_STEP+mem_for_pos) * sizeof(long));
- mem_for_pos += ENLARGEMENT_PER_STEP;
- }
- }
- (*file_positions)[num_of_sequences] = ftell(file);
- }
- }
-
- fclose(file);
- return num_of_sequences;
-}
-
-
-
-/**
- * \brief Reads a profile from a profile file.
- *
- * \param prof A Fastal_profile object to save the profile in.
- * \param profile_f file where the profile is stored.
- * \param position Position of the profile in \a profile_f.
- * \param param_set The parameter set for Fastal
- */
-
-void
-profile_file2profile(Fastal_profile *prof, //structure to save the profile in
- FILE *profile_f, //file where the profile is stored
- long position, //position in profile_f where the profile is stored
- Fastal_param *param_set)
-{
-
- int alphabet_size = param_set->alphabet_size;
-
- int *aa2pos = &(param_set->char2pos[0]);
-
-
- fseek(profile_f,position,SEEK_SET);
- const int LINE_LENGTH = 500;
- char line[500];
-
- fgets(line, LINE_LENGTH, profile_f);
-
- prof->prf_number = atoi(line);
- fgets(line, LINE_LENGTH, profile_f);
- prof->is_leaf = atoi(line);
-
- fgets(line, LINE_LENGTH, profile_f);
- prof->length = atoi(line);
- fgets(line, LINE_LENGTH, profile_f);
- prof->weight = atoi(line);
- fgets(line, LINE_LENGTH, profile_f);
- prof->number_of_sequences = atoi(line);
- int i,j;
- if (prof->length > prof->allocated_memory)
- for (i = 0;i < alphabet_size; ++i)
- {
- prof->prf[i] = vrealloc(prof->prf[i],prof->length*sizeof(int));
- }
- prof->allocated_memory = prof->length;
- char delims[] = " ";
- char *result = NULL;
- char *result_num = NULL;
-
- int length = prof->length;
-
- for (i = 0; i < length; ++i)
- {
- for(j = 0; j<alphabet_size; ++j )
- prof->prf[j][i] = 0;
- fgets(line, LINE_LENGTH , profile_f);
- result = strtok( line, delims );
-
- while( result != NULL)
- {
- result_num = &result[1];
- prof->prf[aa2pos[(short)result[0]]][i] = atoi(result_num);
- result = strtok( NULL, delims );
- }
- }
-}
-
-
-
-/**
- * \brief Writes a profile into a file
- *
- * \param profile Pointer to the profile which has to be saved.
- * \param file A File object (already opened) to write the profile to.
- * \param param_set The parameters for the fastal algorithm.
- */
-void
-profile2file(Fastal_profile *profile, //the profile to save
- FILE* file, //file to save in
- Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
-
- char *pos2aa = &(param_set->pos2char[0]);
-
- fseek(file,0,SEEK_SET);
-
- fprintf(file,"%i\n", profile->prf_number);
-
-
- fprintf(file,"%i\n", profile->is_leaf);
- fprintf(file,"%i\n", profile->length);
- fprintf(file,"%i\n", profile->weight);
- int i = 0, j = 0;
- int max = profile->length;
- int x= 0;
- --alphabet_size;
- while (i < max)
- {
- for (j = 0; j < alphabet_size; ++j)
- if (profile->prf[j][i] > 0)
- {
- if (x)
- fprintf(file," %c%i", pos2aa[j],profile->prf[j][i]);
- else
- fprintf(file,"%c%i", pos2aa[j],profile->prf[j][i]);
- x = 1;
- }
- if (profile->prf[j][i] > 0)
- {
- if (x)
- fprintf(file," %c%i", pos2aa[j],profile->prf[j][i]);
- else
- fprintf(file,"%c%i", pos2aa[j],profile->prf[j][i]);
- x = 1;
- }
- x = 0;
- fprintf(file,"\n");
- ++i;
- }
- fprintf(file,"*\n");
-}
-
-
-
-/**
-* Reads the profile out of an alignment (NOT IN USE)
-*/
-// void
-// file2profile(FILE* profile_f, //file to read the profile of
-// Fastal_profile *prof, //profile saved in here
-// int prf_number, //number of the profile
-// Fastal_param *param_set)
-// {
-// int alphabet_size = param_set->alphabet_size;
-//
-// int *aa2pos = &(param_set->char2pos[0]);
-//
-//
-// fseek(profile_f,0,SEEK_SET);
-// const int LINE_LENGTH = 500;
-// char line[500];
-//
-// fgets(line, LINE_LENGTH, profile_f);
-// prof->prf_number = atoi(line);
-// fgets(line, LINE_LENGTH, profile_f);
-// prof->is_leaf = atoi(line);
-//
-// fgets(line, LINE_LENGTH, profile_f);
-// prof->length = atoi(line);
-//
-// fgets(line, LINE_LENGTH, profile_f);
-// prof->weight = atoi(line);
-// int i,j;
-// if (prof->length > prof->allocated_memory)
-// for (i = 0;i < alphabet_size; ++i)
-// {
-// prof->prf[i] = vrealloc(prof->prf[i],prof->length*sizeof(int));
-// }
-//
-// char delims[] = " ";
-// char *result = NULL;
-// char *result_num = NULL;
-//
-// int length = prof->length;
-//
-// for (i = 0; i < length; ++i)
-// {
-// for(j = 0; j<alphabet_size; ++j )
-// prof->prf[j][i] = 0;
-// fgets(line, LINE_LENGTH , profile_f);
-// result = strtok( line, delims );
-//
-// while( result != NULL)
-// {
-// result_num = &result[1];
-// prof->prf[aa2pos[(short)result[0]]][i] = atoi(result_num);
-// result = strtok( NULL, delims );
-// }
-// }
-// }
-
-
-
-
-
-
-/**
- * \brief Writes the sequence into the alignment_file.
- *
- * \param aligned_sequence Pattern of aligned sequence.
- * \param sequence_file File with sequences.
- * \param sequence_position Positions of sequences in \a sequence_file.
- * \param alignment_file The file to write the sequence into.
- *
-*/
-void
-edit_seq2aligned_seq(char *aligned_sequence, //pattern for aligned sequence
- FILE *sequence_file, //file with all the sequences
- long sequence_position, //position in sequence file with the correct sequence
- FILE *alignment_file) //file to write the alignment into
-{
- fseek(sequence_file, sequence_position, SEEK_SET);
- const int LINE_LENGTH = 300;
- char line[LINE_LENGTH];
- fgets (line, LINE_LENGTH , sequence_file);
- fprintf(alignment_file,"%s", line); //writing of sequence name
- int pos = 0;
- int i = 0;
- while(fgets(line, LINE_LENGTH, sequence_file)!=NULL)
- {
- if (line[0] != '>')
- {
-
- line[LINE_LENGTH-1] = '\n';
- i = 0;
- while ((line[i] != '\n') && (line[i] != '\0'))
- {
- while (aligned_sequence[pos] == '-')
- {
-
- fprintf(alignment_file,"-");
- ++pos;
- }
- if (line[i] != '-')
- {
-
- fprintf(alignment_file,"%c",line[i]);
- ++pos;
- }
- ++i;
-// ++pos;
- }
- }
- else
- break;
- }
- while (aligned_sequence[pos] != '\n')
- {
- fprintf(alignment_file,"-");
- ++pos;
- }
- fprintf(alignment_file,"\n");
-}
-
-
-
-/**
- * \brief Recursive function to turn the edit_file into the alignment.
- *
- * \param sequence_file File with all sequences.
- * \param sequence_position The array of sequence positions in \a sequence_file
- * \param edit_file File to safe the edit profiles in.
- * \param edit_positions Array saving the coorespondence between edit profile and position in \a edit_file
- * \param node_number The current node.
- * \param number_of_sequences The number of sequences.
- * \param aligned_sequence The sequence that is edited.
- * \param alignment_length The length of the alignment.
- * \param edit_seq_file File that saves the edited_sequences of the internal nodes.
- * \param offset Saves the size of the edited_sequences.
- * \param alignment_file File where the alignment is saved.
- */
-void
-edit2alignment(FILE *sequence_file, //sequence file
- long *seq_positions, //sequence positions
- FILE *edit_file, //file saving the edit profiles
- long *edit_positions, //array saving the correspondence between edit profile and position in edit_file
- int node_number, //the current node
- int number_of_sequences, //number of sequences
- char *aligned_sequence, //the sequence that is edited
- int alignment_length, //length of the alignment - and thus of aligned_sequence
- FILE *edit_seq_file, //file saving the edited_sequences of the internal nodes
- int offset, //saves the size of the edited_sequence
- FILE* alignment_file) //file saving the alignments
-{
- fseek(edit_file, edit_positions[node_number-number_of_sequences], SEEK_SET);
- const int LINE_LENGTH = 50;
- char line[LINE_LENGTH];
- fgets(line, LINE_LENGTH , edit_file);
- int child1 = atoi(line);
- fgets(line, LINE_LENGTH , edit_file);
- int child2 = atoi(line);
- fgets(line, LINE_LENGTH , edit_file);
- int is_leaf1 = atoi(line);
- fgets(line, LINE_LENGTH , edit_file);
- int is_leaf2 = atoi(line);
-
-// static char seq_line[10];
-// printf("SO EINE VERDAMMTE SCHEISE ABER AUCH\n");
- char x;
- int number;
- int pos = 0;
-
- //first child
- while(fgets(line, LINE_LENGTH , edit_file)!=NULL)
- {
-
- x = line[0];
- if (x == '*')
- break;
- number = atoi(&line[1]);
- if (x == 'M')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- --number;
- ++pos;
- }
- }
- else if (x == 'I')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- --number;
- ++pos;
- }
- }
- else if (x == 'D')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- {
- aligned_sequence[pos] = '-';
- --number;
- }
- ++pos;
- }
- }
- }
-
- if (is_leaf1)
- {
-// printf("LEAF\n");
- edit_seq2aligned_seq(aligned_sequence, sequence_file, seq_positions[child1], alignment_file);
- }
- else
- {
- fprintf(edit_seq_file, "%s", aligned_sequence);
- edit2alignment(sequence_file, seq_positions, edit_file, edit_positions, child1, number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);
- }
-
- //second child
- fseek(edit_seq_file, offset, SEEK_CUR);
- fgets(aligned_sequence, alignment_length+3, edit_seq_file);
- fseek(edit_seq_file, offset, SEEK_CUR);
-
- pos = 0;
- fseek(edit_file, edit_positions[node_number-number_of_sequences], SEEK_SET);
- while(fgets(line, LINE_LENGTH , edit_file)!=NULL)
- {
- x = line[0];
- if (x == '*')
- break;
- number = atoi(&line[1]);
- if (x == 'M')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- --number;
- ++pos;
- }
- }
- else if (x == 'I')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- {
- aligned_sequence[pos] = '-';
- --number;
- }
- ++pos;
- }
- }
- else if (x == 'D')
- {
- while (number > 0)
- {
- if (aligned_sequence[pos] == 'X')
- --number;
- ++pos;
- }
- }
- }
-
- if (is_leaf2)
- {
- edit_seq2aligned_seq(aligned_sequence, sequence_file, seq_positions[child2], alignment_file);
- }
- else
- {
- fprintf(edit_seq_file, "%s", aligned_sequence);
- edit2alignment(sequence_file, seq_positions, edit_file, edit_positions, child2, number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);
- }
-}
-
-
-
-
-// * The file has the follwing format (# and text behind are only comments and not included into the file):
-// * 1 # Number of profile.
-// * 1 # is leaf.
-// * 5 # Number of columns in the profile.
-// * 4A 1C # In the first column are 4 'A' and 1 'C'
-// * 3G # In the second column are 3 'G'
-// * 5A # In the third column are 5 'A'
-// * 2A 3C # In the fourth column are 2 'A' and 3 'C'
-// * 5C # In the fifth column are 5 'C'
-// * * # Marks the end of this profile
-
-
-
-/**
- * \brief Writes a profile to a file.
- *
- * \param sumup_prf The profile array, not a real profile.
- * \param length The length of the profile. The format can be seen in ./test.txt
- * \param file The FILE object to write the the profile into.
- * \param is_dna The type of sequence.
- * \param number The number of the profile.
- */
-void
-write2file(int **sumup_prf,
- int length,
- FILE *file,
- int number,
- int num_sequences,
- Fastal_param *param_set)
-{
- char *pos2aa = &(param_set->pos2char[0]);
- fprintf(file,"%i\n0\n%i\n1\n%i\n",number, length, num_sequences );
- int i, j;
- int alphabet_size = param_set->alphabet_size;
-
- i = 0;
- int x = 0;
- while (i < length)
- {
- for (j = 0; j < alphabet_size; ++j)
- if (sumup_prf[j][i] > 0)
- {
- if (x)
- fprintf(file," %c%i", pos2aa[j],sumup_prf[j][i]);
- else
- fprintf(file,"%c%i", pos2aa[j],sumup_prf[j][i]);
- x = 1;
- }
- x = 0;
- fprintf(file,"\n");
- ++i;
- }
- fprintf(file,"*\n");
-}
-
-
-
-//************************************* main function of the fasta algorithm ***********************************************
-
-/**
-* \brief main of the fastal algorithm
-*/
-int
-fastal_main(int argc, //number of arguments
- char **argv) //arguments first = fastal, second = tree
-{
-
- int i;
- //pointer to arguments
- void * method_arguments_p;
- int (*alignment_function)(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number);
-
- struct Fastal_arguments arguments;
-
- arg_parse (argc, argv, &arguments);
-
- Fastal_param *param_set = vcalloc(1,sizeof(Fastal_param));
-
- fill_parameters(arguments.is_dna, param_set, arguments.method, arguments.diag_method, arguments.mat);
- param_set->gep = arguments.gep;
- param_set->gop = arguments.gop;
-
-// printf("%s",arguments.mat);
- if (arguments.evaluate)
- {
- printf("Calculate Sum of pairs Score\n");
- printf("Score: %f\n", calculate_sum_of_pairs_score_affine(arguments.sequence_file, param_set->M, param_set->gop, param_set->gep));
- vfree(param_set);
- exit(0);
- }
-
-
- if (arguments.agreement_score)
- {
- complete_agreement_score(arguments.aln2test, arguments.aln_ref);
- return 0;
- }
-
-
- if (arguments.num_ref_aln)
- {
- compute_ref_alignments(arguments.sequence_file, arguments.aln_ref, arguments.num_ref_aln, arguments.num_seq_in_ref);
- return 0;
- }
-
-
-
- int alphabet_size = param_set->alphabet_size;
-
-
- //sequence file management
-// char **seq_name;
- long *file_positions = NULL;
- long **tmp = &file_positions;
- int number_of_sequences = make_index_of_file(arguments.sequence_file, tmp);
-
-
-
- //edit file management
-
-// long current_edit_pos;
- long *edit_positions = vcalloc(number_of_sequences,sizeof(long));
-
-
- //profile management
- Fastal_profile **profiles = vcalloc(3,sizeof(Fastal_profile*));
- initiate_profiles(profiles, param_set);
- FILE * prof_file = fopen(vtmpnam(NULL),"w+");
- long* profile_positions = vcalloc(4,sizeof(long*));
- int max_prof = 4;
- int saved_prof = 0;
-
-
- printf("METHOD: %s\n",param_set->method);
- if (strcmp(param_set->method, "fast") == 0)
- {
- method_arguments_p = vcalloc(1,sizeof(Sparse_dynamic_param));
- fill_arguments_sparse((Sparse_dynamic_param*)method_arguments_p);
- alignment_function = sparse_dyn;
- }
- else if (strcmp(param_set->method, "nw") == 0)
- {
- method_arguments_p = vcalloc(1,sizeof(Nw_param));
- fill_arguments_nw((Nw_param*)method_arguments_p, alphabet_size);
- alignment_function = nw_dyn;
- }
- else if (strcmp(param_set->method, "gotoh") == 0)
- {
- method_arguments_p = vcalloc(1,sizeof(Gotoh_param));
- fill_arguments_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
- alignment_function = gotoh_dyn;
- }
- else if (strcmp(param_set->method, "udisc") == 0)
- {
- method_arguments_p = vcalloc(1,sizeof(Udisc_param));
- fill_arguments_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
- alignment_function = gotoh_dyn;
-
- }
- else
- {
- printf("ERROR - METHOD");
- exit(1);
- }
-
-
- if (arguments.gap_iterate)
- {
- iterate(param_set, method_arguments_p, arguments.sequence_file, arguments.output_file, arguments.gap_iterate);
- return 0;
- }
-
- if (arguments.tree_file == NULL)
- {
- arguments.tree_file = vtmpnam(NULL);
- printf("CONSTRUCT TREE\n");
- if (strcmp(arguments.tree_method, "parttree")==0)
- {
- make_partTree(arguments.sequence_file, arguments.tree_file, arguments.tree_param1, arguments.tree_param2, arguments.is_dna, 0);
- }
- else if (strcmp(arguments.tree_method, "oligotree") == 0)
- {
- compute_oligomer_distance_tree(arguments.sequence_file, param_set->char2pos, arguments.tree_file, arguments.tree_param2, arguments.tree_param1, param_set->alphabet_size);
- }
-
- if (arguments.tree_only == 1)
- return 0;
- }
-
-
- if (arguments.tree_out == 1)
- {
- char tree_out_file_name[500];
- sprintf(tree_out_file_name, "%s.tree",arguments.output_file);
- char const LINE_LENGTH = 50;
- char line[LINE_LENGTH];
-
- FILE* in = fopen(arguments.tree_file, "r");
- FILE* out = fopen(tree_out_file_name, "w");
- while( (fgets(line, LINE_LENGTH, in)) != NULL)
- fprintf(out, "%s", line);
- fclose(in);
- fclose(out);
- }
-
-
-
-
-
- FILE *seq_file = fopen(arguments.sequence_file,"r");
-// FILE *edit_file = fopen(vtmpnam(NULL),"w+");
- FILE *edit_file = fopen("aha","w+");
-
- printf("CONSTRUCT ALIGNMENT\n");
- FILE *tree_file = fopen(arguments.tree_file,"r");
- const int LINE_LENGTH = 100;
- char line[LINE_LENGTH];
- char delims[] = " ";
- int node[3];
- int alignment_length = -1;
- node[2] = -1;
-
-
- //bottom-up traversal
- while(fgets(line, LINE_LENGTH, tree_file)!=NULL)
- {
- //read profiles
- node[0] = atoi(strtok(line,delims));
- node[1] = atoi(strtok(NULL,delims));
- node[2] = atoi(strtok(NULL,delims));
-
- //getting profile of second child
- if (node[1] < number_of_sequences)
- {
- file_pos2profile(seq_file, file_positions[node[1]], profiles[1], node[1], param_set); //profile to save into
- }
- else
- {
- profile_file2profile(profiles[1], prof_file, profile_positions[--saved_prof], param_set);
- fseek (prof_file , profile_positions[saved_prof] , SEEK_SET);
- }
-
- //getting profile of first child
- if (node[0] < number_of_sequences)
- {
- file_pos2profile(seq_file, file_positions[node[0]], profiles[0], node[0], param_set); //profile to save into
- }
- else
- {
- profile_file2profile(profiles[0], prof_file, profile_positions[--saved_prof], param_set);
- fseek (prof_file , profile_positions[saved_prof] , SEEK_SET);
- }
-
-
- if (saved_prof == max_prof)
- {
- max_prof += 5;
- profile_positions = vrealloc(profile_positions, max_prof*sizeof(long));
- }
-
- edit_positions[node[2]-number_of_sequences] = ftell(edit_file);
- profile_positions[saved_prof] = ftell(prof_file);
- ++saved_prof;
-
- //aligning the sequences
- alignment_length = alignment_function(profiles, param_set, method_arguments_p, arguments.is_dna, edit_file, prof_file, node[2]);
- }
-
-
-// bottom-down traversal (edit_files --> alignment)
-// tmp_out_file_name = vtmpnam(NULL);
-
-// FILE *alignment_file = fopen(tmp_out_file_name, "w");
- FILE *alignment_file = fopen(arguments.output_file, "w");
- FILE *edit_seq_file = fopen(vtmpnam(NULL),"w+");
-
- char *aligned_sequence = vcalloc(alignment_length+3, sizeof(char));
-
-
- long offset = ftell(edit_seq_file);
- for (i = 0; i < alignment_length; ++i)
- {
- fprintf(edit_seq_file, "X");
- aligned_sequence[i] = 'X';
- }
- aligned_sequence[i]= '\n';
- aligned_sequence[i+1]= '\0';
- fprintf(edit_seq_file, "\n");
- offset = (ftell(edit_seq_file) - offset)*-1;
-
-
- edit2alignment(seq_file, file_positions, edit_file, edit_positions, node[2], number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);
- fclose(alignment_file);
- fclose(tree_file);
- fclose(edit_file);
- fclose(seq_file);
- fclose(edit_seq_file);
-
- //set stuff for the next cycle
-// arguments.sequence_file = tmp_out_file_name;
-
-
-// // //DEBUG
-// // char copy_command[500];
-// // sprintf(copy_command, "cp %s %s_%i", tmp_out_file_name, arguments.output_file, cycle);
-// // system(copy_command);
-// ++cycle;
-// }
-
-// printf("HERE_COPY\n");
-// char copy_command[2000];
-// sprintf(copy_command, "mv %s %s", tmp_out_file_name, arguments.output_file);
-// printf("%s\n", copy_command);
-// int error = system(copy_command);
-// printf("ERROR %i\n", error);
-
-
- //free_memory & close files
- fclose(prof_file);
- free_fastal_profile(profiles[0], alphabet_size);
- free_fastal_profile(profiles[1], alphabet_size);
- vfree(profiles);
- vfree(profile_positions);
-
-
-
-// number_of_sequences
-
-
- if (arguments.score)
- {
- printf("Calculate Score\n");
- double aln_score = calculate_sum_of_pairs_score_affine(arguments.output_file, param_set->M, param_set->gop, param_set->gep);
- printf("SCORE: %f\n", aln_score);
- }
-
-
-
-
-
- if (!strcmp(param_set->method, "fast"))
- {
- free_sparse((Sparse_dynamic_param*)method_arguments_p);
- }
- else if (!strcmp(param_set->method, "nw"))
- {
- free_nw((Nw_param*)method_arguments_p, alphabet_size);
- }
- else if (!strcmp(param_set->method, "gotoh"))
- {
- free_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
- }
-
- vfree(param_set);
-
- //free_memory & close files
-
- vfree(edit_positions);
-
-
- return 0;
-}
-
-
-
-
-//****************** toolbox ***************************
-
-
-/**
- * \brief Enlargement of the dynamic programming matrix in case it is to small.
- *
- * \param dyn_matrix The dynamic programming matrix.
- * \param old_length1 Current size of dimension 1.
- * \param old_length2 Current size of dimension 2.
- * \param length1 New size of dimension 1.
- * \param length2 New size of dimension 2.
- * \return Pointer to the new array.
- */
-double**
-resize_dyn_matrix(double **dyn_matrix, //the dynamic programming matrix
- int old_length1, //old length of dimension 1
- int old_length2, //old length of dimension 2
- int length1, //new minimum length of dimension 1
- int length2) //new maximum length of dimension 2
-{
- int i;
- if (old_length1 < length1)
- {
- dyn_matrix = vrealloc(dyn_matrix,length1*sizeof(double*));
- for (i = old_length1; i < length1; ++i)
- dyn_matrix[i] = vcalloc(old_length2,sizeof(double));
- old_length1 = length1;
- }
-
- if (old_length2 < length2)
- {
- for (i = 0;i<old_length1; ++i)
- dyn_matrix[i] = vrealloc(dyn_matrix[i], length2*sizeof(double));
- old_length2 = length2;
- }
- return dyn_matrix;
-}
-
-
-
-/**
- * \brief Frees the memory of a dynamic programming matrix.
- *
- * \param length1 Size of the first dimension of the matrix.
- * \param dyn_matrix The dynamic programming matrix.
- */
-void
-free_dyn_matrix(int length1, //length of first dimension
- double **dyn_matrix) //dynamic matrix
-{
- int i = 0;
- for (; i<length1; ++i)
- vfree(dyn_matrix[i]);
- vfree(dyn_matrix);
-}
-
-
-
-/**
- * \brief Initialises the profiles with basic values.
- *
- * \param profiles Array of 3 profiles.
- * \param param_set The fastal parameters
- */
-void
-initiate_profiles(Fastal_profile **profiles, //profiles pointer
- Fastal_param *param_set)
-{
- int alphabet_size = param_set->alphabet_size;
- int i,j;
- for (i =0; i < 3; ++i)
- {
- profiles[i] = vcalloc(1,sizeof(Fastal_profile));
- profiles[i]->weight = 1;
- profiles[i]->is_leaf = 1;
- profiles[i]->prf = vcalloc(alphabet_size, sizeof(int*));
- for (j = 0; j < alphabet_size; ++j)
- {
- profiles[i]->prf[j] = vcalloc(PROFILE_ENLARGEMENT, sizeof(int));
- }
- profiles[i]->allocated_memory = PROFILE_ENLARGEMENT;
- }
-}
-
-
-
-
-/**
- * \brief frees all memory occupied by the profile
- *
- * \param profile The profile to free.
- * \param alphabet_size The alphabet_size.
- */
-void
-free_fastal_profile(Fastal_profile* profile, int alphabet_size)
-{
- --alphabet_size;
- for (;alphabet_size >= 0; --alphabet_size)
- vfree(profile->prf[alphabet_size]);
- vfree(profile->prf);
-}
-
-
-/**
- * \brief Initialize the Fastal parameter set.
- *
- * \param is_dna 1 when sequences are dna, 0 when not
- * \param param_set The fastal parameter set.
- * \param method The method to use in Fastal.
-*/
-void
-fill_parameters(int is_dna, Fastal_param *param_set, char *method, char *diag_method, char *mat)
-{
- sprintf(param_set->method,"%s",method);
- sprintf(param_set->diag_method,"%s",diag_method);
- int i;
- printf("%s",mat);
- param_set->M = read_matrice(mat);
- if (is_dna)
- {
- param_set->alphabet_size = 4;
- char tmp1[] = {'A','C','G','T'};
-
-// int tmp2[] = { 0, 1, 1, 0, -1, -1, 2, 0, -1, -1, 3, -1, 1, 0, -1, -1, -1, 0, 1, 3, 4, -1, 3, -1, 1, -1};
- for (i = 0; i<param_set->alphabet_size; ++i)
- param_set->pos2char[i] = tmp1[i];
-// for (i = 0; i<26; ++i)
-// param_set->char2pos[i] = tmp2[i];
- param_set->char2pos['A'] = 0;
- param_set->char2pos['B'] = 1;
- param_set->char2pos['C'] = 1;
- param_set->char2pos['D'] = 0;
- param_set->char2pos['G'] = 2;
- param_set->char2pos['H'] = 0;
- param_set->char2pos['K'] = 3;
- param_set->char2pos['M'] = 1;
- param_set->char2pos['N'] = 0;
- param_set->char2pos['R'] = 0;
- param_set->char2pos['S'] = 1;
- param_set->char2pos['T'] = 3;
- param_set->char2pos['U'] = 3;
- param_set->char2pos['W'] = 3;
- param_set->char2pos['Y'] = 1;
-// param_set->M[0][3] = param_set->M[3][0] = -10;
-// param_set->M[1][2] = param_set->M[2][1] = -10;
-// param_set->M[0][1] = param_set->M[0][2] = param_set->M[1][0] = param_set->M[2][0] = -10;
-// param_set->M[3][1] = param_set->M[3][2] = param_set->M[1][3] = param_set->M[2][3] = -10;
- }
- else
- {
- param_set->alphabet_size = 21;
- char tmp1[] = {'A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y','X'};
-// int tmp2[] = { 0, 20, 1, 5, 16, 4, 2, 6, 7, 21, 8, 9, 10, 11, -1, 12, 13, 14, 15, 3, -1, 17, 18, 22, 19, 23};
- for (i = 0; i<param_set->alphabet_size; ++i)
- param_set->pos2char[i] = tmp1[i];
-// for (i = 0; i<26; ++i)
-// param_set->char2pos[i] = tmp2[i];
- param_set->char2pos['A'] = 0;
- param_set->char2pos['B'] = 20;
- param_set->char2pos['C'] = 1;
- param_set->char2pos['D'] = 2;
- param_set->char2pos['E'] = 3;
- param_set->char2pos['F'] = 4;
- param_set->char2pos['G'] = 5;
- param_set->char2pos['H'] = 6;
- param_set->char2pos['I'] = 7;
- param_set->char2pos['J'] = 20;
- param_set->char2pos['K'] = 8;
- param_set->char2pos['L'] = 9;
- param_set->char2pos['M'] = 10;
- param_set->char2pos['N'] = 11;
- param_set->char2pos['P'] = 12;
- param_set->char2pos['Q'] = 13;
- param_set->char2pos['R'] = 14;
- param_set->char2pos['S'] = 15;
- param_set->char2pos['T'] = 16;
- param_set->char2pos['V'] = 17;
- param_set->char2pos['W'] = 18;
- param_set->char2pos['X'] = 20;
- param_set->char2pos['Y'] = 19;
- param_set->char2pos['X'] = 20;
- }
-}
-
-
-
-/******************************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