--- /dev/null
+#include "stdio.h"
+#include "stdlib.h"
+#include "math.h"
+#include "ctype.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"
+
+/*!
+ * \file tree.c
+ * \brief Source code for the fastal tree algorithm
+ */
+
+/*
+
+
+int
+main(int argc, char** argv)
+{
+ int alphabet_size = 5;
+ int aa[256];
+ if (alphabet_size == 5)
+ {
+ aa['A'] = 0;
+ aa['B'] = 1;
+ aa['C'] = 1;
+ aa['D'] = 0;
+ aa['G'] = 2;
+ aa['H'] = 0;
+ aa['K'] = 3;
+ aa['M'] = 1;
+ aa['N'] = 0;
+ aa['R'] = 0;
+ aa['S'] = 1;
+ aa['T'] = 3;
+ aa['U'] = 4;
+ aa['W'] = 3;
+ aa['Y'] = 1;
+ }
+ else
+ {
+ aa['A'] = 0;
+ aa['B'] = 1;
+ aa['C'] = 2;
+ aa['D'] = 3;
+ aa['E'] = 4;
+ aa['F'] = 5;
+ aa['G'] = 6;
+ aa['H'] = 7;
+ aa['I'] = 8;
+ aa['J'] = 9;
+ aa['K'] = 10;
+ aa['L'] = 11;
+ aa['M'] = 12;
+ aa['N'] = 13;
+ aa['P'] = 14;
+ aa['Q'] = 15;
+ aa['R'] = 16;
+ aa['S'] = 17;
+ aa['T'] = 18;
+ aa['V'] = 19;
+ aa['W'] = 20;
+ aa['X'] = 21;
+ aa['Y'] = 22;
+ aa['Z'] = 23;
+ }
+ compute_oligomer_distance_tree(argv[1], &aa[0], argv[2], 5, 2, alphabet_size);
+ return 0;
+}*/
+
+
+void
+compute_oligomer_distance_tree(char *seq_file_name, int* char2value, char *tree_file_name, int max_dist, int word_length, int alphabet_size)
+{
+ int i;
+ int *compare = vcalloc(1,sizeof(int));
+ int *num_seq = vcalloc(1,sizeof(int));
+
+ int num_features = (int)pow(alphabet_size,word_length);
+
+ num_features *= num_features * (max_dist+1);
+
+
+ Cluster_info *matrix = feature_extract(seq_file_name, max_dist, word_length, char2value, alphabet_size, compare, num_seq, num_features);
+
+ int *node = vcalloc(1, sizeof(int));
+ *node = *num_seq;
+
+ FILE *tree_f = fopen(tree_file_name,"w");
+ double *mean = vcalloc(num_features, sizeof(double));
+ double *variance = vcalloc(num_features, sizeof(double));
+ cluster_recursive(0, *num_seq-1, matrix, num_features, mean, variance, compare, tree_f, node);
+
+ for (i = 0; i < *num_seq; ++i)
+ {
+// printf("ERG: %i\n", i);
+ vfree(matrix[i].features);
+ }
+ vfree(matrix);
+ fclose (tree_f);
+ vfree(mean);
+ vfree(variance);
+ vfree(compare);
+ vfree(node);
+ vfree(num_seq);
+}
+
+
+
+/**
+ * Recodes a given sequence into a sequence of it's k-mers.
+ *
+ * \param seq The sequence to encode.
+ * \param seq_length The length of \a seq.
+ * \param char2value Array giving each character in seq a unique value.
+ * \param alphabet_size The size of the alphabet.
+ * \param word_length The length of the k-mers to use.
+ *
+ * \return The recodes sequence.
+ */
+int*
+recode_sequence(char * seq, int seq_length, int *char2value, int alphabet_size, int word_length)
+{
+
+ int *recoded = vcalloc(seq_length - word_length + 1, sizeof(int));
+ int i;
+
+ if (word_length == 1)
+ {
+ for (i = 0; i < seq_length; ++i)
+ {
+ recoded[i] = char2value[(short)seq[i]];
+ }
+ }
+ else
+ {
+ int *prod=vcalloc (word_length, sizeof(int));
+ for ( i=0; i<word_length; ++i)
+ {
+ prod[word_length-i-1]=(int)pow(alphabet_size,i);
+ }
+
+ int index = 0;
+ for (i = 0; i < word_length; ++i)
+ {
+ index += char2value[(short)seq[i]] * prod[i];
+ }
+ recoded[0] = index;
+ int z = 0;
+ for (; i < seq_length; ++i)
+ {
+ index -= char2value[(short)seq[z]] * prod[0];
+ index *= alphabet_size;
+ index += char2value[(short)seq[i]];
+ recoded[++z] = index;
+ }
+ }
+
+ return recoded;
+}
+
+
+/**
+ * Extracts the features from a given sequence.
+ *
+ * \param recoded_seq The recoded sequence, consisting of k-mer codes.
+ * \param seq_length The length of \a recoded_seq.
+ * \param k_tup Length of the k-mers.
+ * \param max_coding The maximal number used for coding the alphabet.
+ * \param max_dist The maximal distance to use between to k-mers.
+ *
+ * \return The feature vector of the sequence.
+ */
+double *
+get_features(int *recoded_seq, int seq_length, int k_tup, int max_coding, int max_dist, int num_features)
+{
+
+
+ int vector_length = num_features;//max_coding * max_coding*(max_dist+1);
+
+ double *feature_vector = vcalloc(vector_length, sizeof(double));
+ int i;
+ for (i = 0; i < vector_length; ++i)
+ {
+ feature_vector[i] = 0;
+ }
+ int j, max_length, indJ, indK, indDist;
+ for (i = 0; i < seq_length; ++i)
+ {
+ if (seq_length <= i + max_dist)
+ max_length = seq_length-1;
+ else
+ max_length = i + max_dist;
+ for (j = i; j <= max_length; ++j)
+ {
+
+ indJ = max_coding * (max_dist+1) * recoded_seq[i];
+ indK = (max_dist+1) * recoded_seq[j];
+ indDist = j - i;
+ ++feature_vector[indJ + indK + indDist];
+ }
+ }
+
+ double normalize = 0;
+ for (i = 0; i < vector_length; ++i)
+ {
+ normalize += (feature_vector[i] * feature_vector[i]);
+ }
+ normalize = sqrt(normalize);
+ for (i = 0; i < vector_length; ++i)
+ {
+ feature_vector[i] /= normalize;
+ }
+
+ return feature_vector;
+}
+
+
+/**
+ * Calculates the feature values for all sequences in a file.
+ *
+ * The File has to have fasta format.
+ * \param seq_file_name The name of the file (in fasta format).
+ * \param max_dist The maximal distance to be considered.
+ * \param k_tup Size of the k_tup.
+ * \param alphabet_size The size of the alphabet.
+ *
+ * \return The feature values for every sequence in a Cluster_info object.
+ */
+Cluster_info *
+feature_extract(char *seq_file_name, int max_dist, int k_tup, int *char2value, int alphabet_size, int * elem_2_compare, int *num_seq_p, int num_features)
+{
+ char line[500];
+ FILE *tree_f = fopen(seq_file_name,"r");
+// fgets(line, 500, tree_f);
+ int num_seq = -1;
+ int size = 1000;
+ int matrix_size = 1000;
+ const int STEP = 1000;
+ int seq_pos = -1;
+ char *seq = vcalloc(size, sizeof(char));
+ Cluster_info *matrix = vcalloc(matrix_size, sizeof(Cluster_info));
+ char *c_p;
+ int max_coding = pow(alphabet_size,k_tup);
+ while (fgets(line, 500, tree_f) != NULL)
+ {
+ if (line[0] == '>')
+ {
+ if (num_seq >= 0)
+ {
+ seq[++seq_pos] = '\0';
+ if (num_seq > matrix_size -2)
+ {
+ matrix_size += STEP;
+ matrix = vrealloc(matrix, matrix_size * sizeof(Cluster_info));
+ }
+ int *recoded_seq = recode_sequence(seq, seq_pos, char2value, alphabet_size, k_tup);
+ matrix[num_seq].seq_number = num_seq;
+ matrix[num_seq].elem_2_compare = elem_2_compare;
+ matrix[num_seq].features=get_features(recoded_seq, seq_pos- k_tup +1, k_tup, max_coding, max_dist, num_features);
+ vfree(recoded_seq);
+ seq_pos = -1;
+ }
+ ++num_seq;
+ }
+ else
+ {
+ if (size - seq_pos < 500)
+ {
+ size += STEP;
+ seq = vrealloc(seq, size*sizeof(char));
+ }
+ c_p = &line[0];
+ while ((*c_p != '\n') && (*c_p != '\0'))
+ {
+ seq[++seq_pos] = toupper(*(c_p));
+ ++c_p;
+ }
+ }
+ }
+
+ fclose(tree_f);
+ seq[++seq_pos] = '\0';
+ int *recoded_seq = recode_sequence(seq, seq_pos, char2value, alphabet_size, k_tup);
+ matrix[num_seq].seq_number = num_seq;
+ matrix[num_seq].elem_2_compare = elem_2_compare;
+ matrix[num_seq].features=get_features(recoded_seq, seq_pos- k_tup +1, k_tup, max_coding, max_dist, num_features);
+ *num_seq_p = ++num_seq;
+ vfree(seq);
+ vfree(recoded_seq);
+ return matrix;
+}
+
+
+
+/**
+ * Compare function for qsort given an array of Cluster_info.
+ *
+ * The field used for sorting is done according to the element determined in the Cluster_info object.
+ * \param a void pointer to an object of Cluster_info.
+ * \param b void pointer to an object of Cluster_info.
+ *
+ * \return Value showing wheather \a a is bigger, smaller or equal to \a b.
+ */
+int
+cluster_compare (const void * a, const void * b)
+{
+ Cluster_info *tmp1 = (Cluster_info*)a;
+ Cluster_info *tmp2 = (Cluster_info*)b;
+ int elem = *(tmp1->elem_2_compare);
+ if ((tmp1->features[elem]) > (tmp2->features[elem]))
+ return 1;
+ else if ((tmp1->features[elem]) < (tmp2->features[elem]))
+ return -1;
+ else
+ return 0;
+}
+
+
+
+/**
+ * Recursive function to build a tree according to a given feature matrix.
+ *
+ * \param start The begin of this set in \a matrix.
+ * \param end The end of this set in \a matrix.
+ * \param matrix The feature matrix.
+ * \param dim The number of features.
+ * \param mean Array of size \a dim.
+ * \param variance Array of size \a dim.
+ * \param elem_2_compare Pointer to a single in. This will save the feature according to which shall be sorted (biggest variance).
+ * \param tree_f The file in which the tree will be written.
+ * \param node The current node.
+ *
+ * \return The number of the current node.
+ */
+int
+cluster_recursive(int start,
+ int end,
+ Cluster_info* matrix,
+ int dim,
+ double *mean,
+ double *variance,
+ int *elem_2_compare,
+ FILE* tree_f,
+ int *node)
+{
+
+ //stop recursion
+ if (end-start < 1)
+ {
+ return matrix[start].seq_number;
+ }
+
+ //reset mean/variance
+ int i;
+ int num_seq = end - start +1;
+ for (i = 0; i < dim; ++i)
+ {
+ mean[i] = 0;
+ variance[i] = 0;
+ }
+
+ //calculate mean
+ Cluster_info *start_p = &matrix[start]-1;
+ Cluster_info *end_p = &matrix[end];
+ double *tmp;
+ while (start_p < end_p)
+ {
+ tmp = &(++start_p)->features[0];
+ for (i = 0; i < dim; ++i)
+ {
+ mean[i] += *(tmp++);
+ }
+ }
+ for (i = 0; i < dim; ++i)
+ mean[i] /= num_seq;
+
+
+ //calculate variance
+ start_p = &matrix[start] -1;
+
+ while (start_p < end_p)
+ {
+ tmp = &(++start_p)->features[0];
+ for (i = 0; i < dim; ++i)
+ {
+ variance[i] += ((*tmp) - mean[i]) * ((*tmp) - mean[i]);
+ ++tmp;
+ }
+ }
+
+ //get maximal variance
+ double max = variance[0];
+ int index = 0;
+ for (i = 1; i < dim; ++i)
+ {
+ if (variance[i] > max)
+ {
+ max = variance[i];
+ index = i;
+ }
+ }
+ *elem_2_compare = index;
+
+
+ //devide into to different sets and start recursion on these child sets
+ qsort (&matrix[start], num_seq, sizeof(Cluster_info), cluster_compare);
+ int split = start + (end-start)/2;
+ int split_value = matrix[split].features[index];
+ while ((matrix[split+1].features[index] == split_value) && (split < end-1))
+ ++split;
+
+ int left_child = cluster_recursive(start, split, matrix, dim, mean, variance, elem_2_compare, tree_f, node);
+ int right_child = cluster_recursive(split+1, end, matrix, dim, mean, variance, elem_2_compare, tree_f, node);
+ int node2 = *node;
+ ++(*node);
+ //print node into file
+ fprintf(tree_f, "%i %i %i\n", left_child, right_child, node2);
+
+ return node2;
+}
+
+/******************************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