--- /dev/null
+/******************************COPYRIGHT NOTICE*******************************/
+/* (c) Centro de Regulacio Genomica */
+/* and */
+/* Cedric Notredame */
+/* 12 Aug 2014 - 22:07. */
+/*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*******************************/
+
+// #include "km_coffee.h"
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <limits.h>
+#include <float.h>
+#include <math.h>
+#include <stdarg.h>
+#include <fcntl.h>
+#include <sys/file.h>
+#include <sys/types.h>
+#include <unistd.h>
+
+#include "io_lib_header.h"
+#include "util_lib_header.h"
+#include "define_header.h"
+
+#include "km_coffee_header.h"
+
+/**
+ * \brief Writes the sequences to the file
+ * \param node The node containing the data to write
+ * \param vecs The vector set
+ * \param seq_set The sequence Set
+ * \param name_f The file to write the sequence into.
+ */
+void write_files(KM_node *node, int *vecs, const SeqSet *seq_set, char *name_f)
+{
+ size_t end = node->end;
+ FILE * tmp_F = fopen(name_f, "w");
+ size_t i;
+ for (i = node->start; i<end; ++i)
+ {
+ fprintf(tmp_F, ">%s\n%s\n", seq_set->seqs[vecs[i]]->name, seq_set->seqs[vecs[i]]->seq);
+ }
+ fclose(tmp_F);
+}
+
+/**
+ * \brief Frees all memory occupied by this structure
+ * \param root The root of the tree to delete.
+ */
+void
+del_tree(KM_node* root)
+{
+ Stack *to_do = Stack_init();
+
+ Node_pair *tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=root;
+ tmp->id = 0;
+ push(to_do, tmp);
+
+ KM_node* current;
+ size_t child;
+ while (to_do->size != 0)
+ {
+ current = ((Node_pair*)to_do->last)->node;
+ child = ((Node_pair*)to_do->last)->id;
+ if (current->n_children==0)
+ {
+ free(pop(to_do));
+ delKM_node(current);
+ }
+ else
+ {
+ if (child == current->n_children)
+ {
+ delKM_node(current);
+ free(pop(to_do));
+ }
+ else
+ {
+ ++((Node_pair*)to_do->last)->id;
+ tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=current->children[child];
+ tmp->id = 0;
+ push(to_do, tmp);
+ }
+ }
+ }
+ delStack(to_do);
+}
+
+
+/**
+ * \brief Traverses the tree and calls t_coffee to produce the alignments
+ * \param root The root of the tree.
+ * \param vecs The vector set.
+ * \param seq_set The sequence set.
+ * \param out_f The output file.
+ * \param n_core The number of cores to use (OPEN-MP for kmeans/fork in T-Coffee)
+ * \param gapopen The gapopening costs
+ * \param gapext The gapextension costs
+ * \param method The method to use in the alignment
+ */
+void
+traverse_km_tree(KM_node* root, int *vecs, const SeqSet *seq_set, char *out_f, int n_cores, int gapopen, int gapext, char *method)
+{
+ Stack *to_do =Stack_init();
+ Node_pair *tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=root;
+ tmp->id = 0;
+ push(to_do, tmp);
+
+ KM_node* current;
+ size_t child;
+ size_t j;
+// size_t pos;
+ char command[1000];
+ while (to_do->size != 0)
+ {
+
+ current = ((Node_pair*)to_do->last)->node;
+ child = ((Node_pair*)to_do->last)->id;
+ if (current->n_children==0)
+ {
+ if(current->end-current->start > 2)
+ {
+ sprintf(command, "%li",current->id);
+ write_files(current, vecs, seq_set, command);
+// sprintf(command,"clustalo --guidetree-out co_tree.dnd -i %li --force >/dev/null 2>/dev/null", current->id);
+// if (system(command))
+// {
+// printf("%s\n",command);
+// exit(1);
+// }
+ if (current->id!=0)
+ sprintf(command,"t_coffee -in %li -output fasta_aln -outfile %li.fa -n_core %i -gapopen %i -gapext %i -method %s -quiet >/dev/null 2>/dev/null", current->id, current->id, n_cores, gapopen, gapext, method);
+ else
+ sprintf(command,"t_coffee -in %li -output fasta_aln -outfile %s -n_core %i -gapopen %i -gapext %i -method %s -quiet >/dev/null 2>/dev/null ", current->id, out_f, n_cores, gapopen, gapext, method);
+ if (system(command))
+ {
+ printf("ERROR when running: %s\n",command);
+ exit(1);
+ }
+ }
+ else
+ if(current->end-current->start > 1)
+ {
+ sprintf(command, "%li",current->id);
+ write_files(current, vecs, seq_set, command);
+ if (current->id!=0)
+ sprintf(command,"t_coffee -in %li -output fasta_aln -outfile %li.fa -n_core %i -gapopen %i -gapext %i -method %s -quiet >/dev/null 2>/dev/null", current->id, current->id, n_cores, gapopen, gapext, method);
+ else
+ sprintf(command,"t_coffee -in %li -output fasta_aln -outfile %s -n_core %i -gapopen %i -gapext %i -method %s -quiet >/dev/null 2>/dev/null ", current->id, out_f, n_cores, gapopen, gapext, method);
+
+ printf("%s\n", command);
+ if (system(command))
+ {
+ printf("ERROR when running: %s\n",command);
+ exit(1);
+ }
+ }
+ else
+ {
+ sprintf(command, "%li.fa",current->id);
+ write_files(current, vecs, seq_set, command);
+ }
+ tmp =(Node_pair*) pop(to_do);
+ }
+ else
+ {
+ if (child == current->n_children)
+ {
+ if (current->id!=0)
+ sprintf(command, "t_coffee -output fasta_aln -method %s -quiet -outfile %li.fa -n_core %i -gapopen %i -profile FILE::prf.fa ", method, current->id,n_cores, gapopen);
+ else
+ sprintf(command, "t_coffee -output fasta_aln -method %s -quiet -outfile %s -n_core %i -gapopen %i -profile FILE::prf.fa ", method, out_f, n_cores, gapopen);
+ FILE *prf_F = my_fopen("prf.fa", "w");
+ for(j=0; j<current->n_children;++j)
+ {
+ fprintf(prf_F, "%li.fa\n", current->children[j]->id);
+ }
+ fclose(prf_F);
+
+ printf("%s\n", command);
+ if (system(command))
+ {
+
+ printf("ERROR when running: %s\n",command);
+ exit(1);
+ }
+ pop(to_do);
+ }
+ else
+ {
+ ++((Node_pair*)to_do->last)->id;
+ tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=current->children[child];
+ tmp->id = 0;
+ push(to_do, tmp);
+ }
+ }
+ }
+ exit(1);
+}
+
+
+int
+my_seq_sort (const void *i, const void *j)
+{
+ return strcmp((*(Seq**)i)->seq,(*(Seq**)j)->seq);
+}
+
+
+
+
+void
+print_km_tree(KM_node *root, int *vecs, const SeqSet *seq_set, char *out_f)
+{
+ FILE *out_F = fopen(out_f, "w");
+ Stack *to_do =Stack_init();
+ Node_pair *tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=root;
+ tmp->id = 0;
+ push(to_do, tmp);
+
+ KM_node* current;
+ size_t child;
+ size_t k;
+ // size_t pos;
+// char command[1000];
+ while (to_do->size != 0)
+ {
+
+ current = ((Node_pair*)to_do->last)->node;
+ child = ((Node_pair*)to_do->last)->id;
+ if (current->n_children==0)
+ {
+ if (current->end-current->start >1)
+ fprintf(out_F, "(");
+ for (k=current->start; k<current->end; ++k)
+ {
+ fprintf(out_F, "%s", seq_set->seqs[vecs[k]]->name);
+ if (k<current->end-1)
+ fprintf(out_F, ",");
+ }
+ tmp =(Node_pair*) pop(to_do);
+ if (current->end-current->start >1)
+ fprintf(out_F, ")");
+ }
+ else
+ {
+ if (child == 0)
+ {
+ fprintf(out_F, "(");
+ }
+ if (child == current->n_children)
+ {
+ fprintf(out_F, ")");
+ pop(to_do);
+ }
+ else
+ {
+ if (child != 0)
+ fprintf(out_F, ",");
+ ++((Node_pair*)to_do->last)->id;
+ tmp = (Node_pair*)my_malloc(sizeof(Node_pair));
+ tmp->node=current->children[child];
+ tmp->id = 0;
+ push(to_do, tmp);
+ }
+ }
+ }
+ fprintf(out_F, ";");
+}
+
+
+typedef struct
+{
+ size_t id;
+ float sim;
+} sorter;
+
+
+
+
+int
+my_val_compare (const void *i, const void *j)
+{
+ return (*(Vector**)j)->val < (*(Vector**)j)->val;
+}
+
+/*
+KM_node*
+simple_clust(VectorSet *vecSet, unsigned int k)
+{
+ size_t n_vecs = vecSet->n_vecs;
+ size_t i, max_id=0, max_val=0;
+ Vector **vecs = vecSet->vecs;
+ size_t dim=vecSet->dim;
+
+ // determine longest_seq
+ for (i=0; i<n_vecs; ++i)
+ {
+ if (vecs[i]->seq_len > max_val)
+ {
+ max_val=vecs[i]->seq_len;
+ max_id=i;
+ }
+ }
+
+ size_t num_nodes= ceil(n_vecs*1.0/k);
+ KM_node **nodes = (KM_node**)malloc(sizeof(KM_node*)*num_nodes);
+ size_t j;
+
+ //determine leaves
+ size_t node_id=0;
+ for (j=0; j<num_nodes; ++j)
+ {
+ for (i=0; i<n_vecs; ++i)
+ vecs[i]->val=km_sq_dist(vecs[max_id], vecs[i],dim);
+ qsort (vecs, n_vecs, sizeof(Vector*), my_val_compare);
+ nodes[node_id] = malloc(sizeof(KM_node));
+ nodes[node_id]->children = NULL;
+ nodes[node_id]->start = node_id*k;
+ nodes[node_id]->end = (n_vecs>k )? (node_id+1)*k: (node_id*k)+n_vecs;
+ nodes[node_id]->n_children=0;
+ nodes[node_id]->id=++node_id;
+ n_vecs -= k;
+ max_id =0;
+ vecs+=k;
+ }
+
+ // generate inner nodes
+ j=num_nodes;
+ KM_node *tmp;
+ size_t pos=0;
+ size_t overall_pos=0;
+// printf("1: %li %li\n", num_nodes, node_id);
+ while (1)
+ {
+ for(i=0; i<num_nodes; ++i)
+ {
+// printf("RUN %li\n", num_nodes);
+ if ((i%k)==0)
+ {
+ tmp=malloc(sizeof(KM_node));
+ tmp->children=malloc(sizeof(KM_node*)*num_nodes);
+ tmp->children[0]=nodes[i];
+ tmp->n_children=1;
+ tmp->id=++node_id;
+ nodes[overall_pos]=tmp;
+ pos=0;
+ ++overall_pos;
+ }
+ else
+ {
+ tmp->children[++pos]=nodes[i];
+ ++tmp->n_children;
+ }
+ }
+ overall_pos=0;
+ if (num_nodes==1)
+ break;
+ num_nodes=ceil(num_nodes*1.0/k);
+
+ }
+// printf("%li\n", node_id);
+ free(nodes);
+ tmp->id=0;
+ return tmp;
+}*/
+
+/*
+typedef struct KM_node{
+ size_t start;
+ size_t end;
+ struct KM_node **children;
+ size_t id;
+ size_t n_children;
+ } KM_node;*/
+
+
+
+int
+km_coffee_align3(char *seq_f, int k, int k_leaf, char *method, char *aln_f, int n_cores, int gapopen, int gapext, char *init)
+{
+ char *use_as_temp = get_tmp_4_tcoffee();
+
+ #ifdef _OPENMP
+ omp_set_num_threads(n_cores);
+ #endif
+
+ SeqSet *seq_set = read_fasta(seq_f);
+ qsort(seq_set->seqs, seq_set->n_seqs, sizeof(Seq*), my_seq_sort);
+ srand(time(0));
+
+
+ short j = -1;
+ short i;
+ /****************************************************
+ Sequences to vector using k-mers
+ *****************************************************/
+ short alphabet[256];
+
+ // standard alphabet
+ for (i = 65; i < 91; ++i)
+ if ((i==66) || (i==74) || (i==79) || (i==88) || (i==90))
+ alphabet[i] = 0;
+ else
+ alphabet[i] = ++j;
+ j=-1;
+ for (i = 97; i < 123; ++i)
+ if ((i==98) || (i==106) || (i==111) || (i==120) || (i==122))
+ alphabet[i] = 0;
+ else
+ alphabet[i] = ++j;
+
+ // shrinked alphabet
+// for (i = 0; i < 256; ++i)
+// alphabet[i] = 0;
+
+// char *groups[]={"LlVvIiMmCcAaGgSsTtPpFfYyWw","EeDdNnQqKkRrHh"};
+// char *groups[]={"LlVvIiMmCc","AaGgSsTtPp","FfYyWw","EeDdNnQqKkRrHh"};
+// size_t n_groups = 4;
+// size_t len;
+// char *group;
+// for (i=0; i<n_groups; ++i)
+// {
+// group=groups[i];
+// len=strlen(group);
+// for (j=0; j<len; ++j)
+// alphabet[group[j]]=i+1;
+// }
+
+ VectorSet *vec_set = seqset2vecs_kmer(seq_set, 2, 21, alphabet);
+
+
+ /****************************************************
+ Sequences to vector using distances
+ *****************************************************/
+ // char *groups[]={"LVIMC","AGSTP","FYW","EDNQKRH"};
+// size_t n_groups = 4;
+// char *groups[]={"LlVvIiMmCc","AaGgSsTtPp","FfYyWw","EeDdNnQqKkRrHh"};
+// char *groups[]={"LlVvIiMmCcAaGgSsTtPpFfYyWw","EeDdNnQqKkRrHh"};
+// size_t n_groups = 2;
+// VectorSet *vec_set = seqset2vecs_whatever(seq_set, groups, n_groups);
+
+
+
+
+
+
+// char vec_file[500];
+// sprintf(vec_file, "%s_2_8_%li_%li.txt", strrchr(seq_f, '/')+1, vec_set->n_vecs, vec_set->dim);
+
+// print_vecs(vec_set, &vec_file[0]);
+// read_vecs(vec_set, "matrix_59");
+// exit(1);
+// normalize(vec_set);
+ KM_node *root = hierarchical_kmeans(vec_set, k, k_leaf, init, 0.001);
+// KM_node *root = simple_clust(vec_set, k);
+
+
+ char templatee[400];
+ sprintf(templatee, "%s/km_coffee_tmp_XXXXXX", use_as_temp);
+ char tmp_str[FILENAME_MAX];
+ km_cwd = getcwd(tmp_str, FILENAME_MAX);
+
+ km_tmp_dir = my_make_temp_dir(templatee, "main");
+ chdir(km_tmp_dir);
+ char out_f[500];
+ if (aln_f[0] != '/')
+ sprintf(out_f, "%s/%s", km_cwd, aln_f);
+ else
+ sprintf(out_f, "%s", aln_f);
+
+
+
+ size_t n_vecs = seq_set->n_seqs;
+ int *assignment = (int*)malloc(n_vecs*sizeof(int));
+ size_t l;
+ for (l = 0; l< n_vecs; ++l)
+ assignment[l]=vec_set->vecs[l]->id;
+
+// printf("TRAVERSE\n");
+ delVecSet(vec_set);
+ traverse_km_tree(root, assignment, seq_set, out_f, n_cores, gapopen, gapext, method);
+ free( assignment);
+ del_tree(root);
+ delSeqSet(seq_set);
+
+ free(km_tmp_dir);
+
+
+
+
+ return EXIT_SUCCESS;
+}
git://source.jalview.org / jabaws.git / blobdiff