--- /dev/null
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+// #include <ctype.h>
+// #include <string.h>
+
+#include "io_lib_header.h"
+#include "util_lib_header.h"
+#include "fastal_lib_header.h"
+
+
+/**
+ * Adds a single sequence to a profile.
+ * \param seq The sequence.
+ * \param prf The profile.
+ */
+void
+seq2profile2(char *seq, Fastal_profile *prf, int *char2pos)
+{
+// printf("SEDQ: %s\n", seq);
+ int **prf_in = prf->prf;
+ int i = -1;
+ --seq;
+ while (*(++seq) !='\0')
+ {
+ if (*seq != '-')
+ ++(prf_in[char2pos[*seq]][++i]);
+ else
+ ++i;
+// ++seq;
+ }
+}
+
+
+
+
+// typedef struct
+// {
+// /// Number of sequences in this profile
+// int num_sequences;
+// /// number of the profile
+// int prf_number;
+// ///0 = combination of two profiles, 1 = profile of a single sequence -> name1 = seq_name
+// int is_leaf;
+// ///length of the profile
+// int length;
+// ///weight of the sequence
+// int weight;
+// ///saves the amount of allocated memory
+// int allocated_memory;
+// ///the profile itself [alphabet_size][profile_length]
+// int **prf;
+// ///number_of_sequences
+// int number_of_sequences;
+// }
+// Fastal_profile;
+
+
+
+
+
+/**
+ * Deletes all gap columns from the profile.
+ *
+ * \param prf The profile.
+ * \param alphabet_size The size of the alphabet.
+ * \param gap_list The gap_list.
+ * \param gap_list_length The length of the gap list.
+ * \param num_gaps The number of gaps.
+ * \return The new gap list.
+ */
+int *
+del_gap_from_profile(Fastal_profile *prf, int alphabet_size, int *gap_list, int *gap_list_length, int *num_gaps)
+{
+ int i;
+ int pos = -1;
+ int not_gap = 0;
+ int gap_pos = 0;
+ int write_pos = 0;
+// int gap_counter = 0;
+ int **prf_in = prf->prf;
+ int prf_length = prf->length;
+
+ while (gap_pos < prf_length)
+ {
+// printf("PRF1: %i %i %i %i\n", prf_in[0][gap_pos], prf_in[1][gap_pos], prf_in[2][gap_pos], prf_in[3][gap_pos]);
+
+ not_gap = 0;
+ for (i = 0; i < alphabet_size; ++i)
+ {
+ if (prf_in[i][gap_pos])
+ {
+ not_gap = 1;
+// printf("ARG: %i\n", i);
+ break;
+ }
+ }
+ if (not_gap)
+ {
+
+ if (gap_pos != write_pos)
+ {
+ for (i = 0; i < alphabet_size; ++i)
+ {
+ prf_in[i][write_pos] = prf_in[i][gap_pos];
+ }
+ }
+ ++write_pos;
+ }
+ else
+ {
+
+ if (++pos >= *gap_list_length)
+ {
+ *gap_list_length += 10;
+ gap_list = vrealloc(gap_list, (*gap_list_length)*sizeof(int));
+ }
+ gap_list[pos] = gap_pos;
+// ++gap_counter;
+ }
+ ++gap_pos;
+ }
+// printf("DEEEEEEEEEEEEEEEEEL %i %i\n", gap_counter, pos+1);
+
+ *num_gaps = pos+1;
+ prf->length -= *num_gaps;
+ return gap_list;
+}
+
+
+/**
+ * Splits a given alignment according to a given position.
+ * \param alignment_file The alignment file.
+ * \param gap_prf Saves the profile with the gap sequences.
+ * \param no_gap_prf Saves the profile of the non gap sequences.
+ * \param seq A temporary place to save the sequence.
+ * \param index The index where to look for the gaps.
+ * \param char2pos The character to position translation array.
+ */
+void
+split_set(FILE *alignment_file, Fastal_profile *gap_prf, Fastal_profile *no_gap_prf, char *seq, int index, int *char2pos, char* split_file_name)
+{
+ FILE *split_set_file1 = fopen("GAP", "w");
+ FILE *split_set_file2 = fopen("NOGAP", "w");
+// FILE *shorted = fopen("SHORTT.fa", "w");
+ printf("INDEX: %i\n",index);
+ fseek (alignment_file , 0 , SEEK_SET);
+ FILE *split_file = fopen(split_file_name, "w");
+ int i,j;
+ const int LINE_LENGTH = 200;
+ char line[LINE_LENGTH];
+ int pos = -1;
+ int overall_pos = -1;
+ fgets(line, LINE_LENGTH , alignment_file);
+// fprintf(shorted, "%s", line);
+ while(fgets(line, LINE_LENGTH , alignment_file)!=NULL)
+ {
+ pos = -1;
+ if (line[0] == '>')
+ {
+// fprintf(shorted, "%s", line);
+ if (seq[index] == '-')
+ {
+
+ seq[++overall_pos] = '\0';
+ fprintf(split_file, "g\n");
+ fprintf(split_set_file1,"%s\n",seq);
+ seq2profile2(seq, gap_prf, char2pos);
+ ++(gap_prf->number_of_sequences);
+ }
+ else
+ {
+ seq[++overall_pos] = '\0';
+ fprintf(split_file, "n\n");
+ fprintf(split_set_file2,"%s\n",seq);
+ seq2profile2(seq, no_gap_prf, char2pos);
+ ++(no_gap_prf->number_of_sequences);
+ }
+ overall_pos = -1;
+
+ }
+ else
+ {
+// for (j = 35; j < 96; ++j)
+// {
+// fprintf(shorted, "%c", line[j]);
+// }
+// fprintf(shorted, "\n");
+ while ((line[++pos] != '\n') && (line[pos] != '\0'))
+ {
+ seq[++overall_pos] = line[pos];
+ }
+ }
+ }
+
+ if (seq[index] == '-')
+ {
+ seq[++overall_pos] = '\0';
+ fprintf(split_file, "g\n");
+ seq2profile2(seq, gap_prf, char2pos);
+ ++(gap_prf->number_of_sequences);
+ }
+ else
+ {
+ seq[++overall_pos] = '\0';
+ fprintf(split_file, "n\n");
+ seq2profile2(seq, no_gap_prf, char2pos);
+ ++(no_gap_prf->number_of_sequences);
+ }
+ fclose(split_file);
+}
+
+
+Fastal_profile *
+enlarge_prof(Fastal_profile *prof, int new_length, int alphabet_size)
+{
+ if (new_length > prof->allocated_memory)
+ {
+ int i;
+ for (i = 0;i < alphabet_size; ++i)
+ {
+ prof->prf[i] = vrealloc(prof->prf[i],new_length*sizeof(int));
+ }
+ prof->allocated_memory = new_length;
+ }
+ return prof;
+}
+
+
+
+
+void
+iterate(Fastal_param *param, void *method_arguments_p, char *aln_file_name, char *out_file_name_end, int iteration_number)
+{
+// calculate alignment length
+
+//count gap
+ int it_coutner_2 = 0;
+ const int LINE_LENGTH = 200;
+ char line[LINE_LENGTH];
+ char *seq1 = vcalloc(1,sizeof(char));
+ char *seq2 = vcalloc(1,sizeof(char));
+ Fastal_profile **profiles = vcalloc(3,sizeof(Fastal_profile*));
+ initiate_profiles(profiles, param);
+ Fastal_profile *gap_prf = profiles[0];
+ Fastal_profile *no_gap_prf = profiles[1];
+ int alphabet_size = param->alphabet_size;
+
+ int *gap_list_1 = vcalloc(1, sizeof(int));
+ int *gap_list_1_length = vcalloc(1, sizeof(int));
+ *gap_list_1_length = 1;
+ int num_gaps_1 = 0;
+ int *gap_list_2 = vcalloc(1, sizeof(int));
+ int *gap_list_2_length = vcalloc(1, sizeof(int));
+ *gap_list_2_length = 1;
+ int num_gaps_2 = 0;
+
+ int alignment_length = 1;
+ //from here repeat!
+ int it_counter = 0;
+ char *out_file_name = aln_file_name;
+ int *gap_profile = vcalloc(alignment_length, sizeof(int));
+
+// while (it_counter < alignment_length)
+// {
+
+
+ FILE *alignment_file = fopen(aln_file_name,"r");
+ if (alignment_file == NULL)
+ {
+ printf("Could not open alignment file %s\n", aln_file_name);
+ exit(1);
+ }
+
+
+
+ alignment_length = 0;
+ int tmp_len = -1;
+ fgets(line, LINE_LENGTH , alignment_file);
+ while(fgets(line, LINE_LENGTH , alignment_file)!=NULL)
+ {
+ tmp_len = -1;
+ if (line[0] == '>')
+ {
+ break;
+ }
+ while ((line[++tmp_len] != '\n') && (line[tmp_len] != '\0'));
+ alignment_length += tmp_len;
+ }
+ // printf("ALN_LENGTH %i\n", alignment_length);
+ seq1 =vrealloc(seq1, (1+alignment_length)*sizeof(char));
+
+ gap_profile = vrealloc(gap_profile, alignment_length * sizeof(int));
+ int i;
+ for (i = 0; i < alignment_length; ++i)
+ {
+ gap_profile[i] = 0;
+ }
+
+ int number_of_sequences = 0;
+ int pos = -1;
+ fseek (alignment_file , 0 , SEEK_SET);
+ while(fgets(line, LINE_LENGTH , alignment_file)!=NULL)
+ {
+ if (line[0] == '>')
+ {
+ ++number_of_sequences;
+ pos = -1;
+ }
+ else
+ {
+ i = -1;
+ while ((line[++i] != '\n') && (line[i] != '\0'))
+ {
+ ++pos;
+ if (line[i] == '-')
+ {
+ ++gap_profile[pos];
+ }
+ }
+ }
+ }
+
+ double max_entrop = 0;
+ int column_index = 0;
+ double entrop = 0;
+ double last = 0;
+ double p;
+
+
+// for (i = it_counter; i<=it_counter; ++i)
+// {
+// p = gap_profile[i]/(double)number_of_sequences;
+// if (!p)
+// {
+// entrop = 0;
+// }
+// else
+// {
+// entrop = (-1)*(p*log(p) + (1-p)*log(1-p) ) ;
+// }
+// if (entrop > max_entrop)
+// {
+// column_index = i;
+// max_entrop = entrop;
+// }
+// last = entrop;
+// }
+// ++it_counter;
+// if (max_entrop < 0.6)
+// {
+// printf("CONTINUE %f\n",entrop);
+// continue;
+// }
+
+// column_index = 18;//it_counter-1;
+// if (column_index == 19)
+ column_index = 58;
+ out_file_name = vtmpnam(NULL);
+
+ char *edit_file_name = "EDIT";//vtmpnam(NULL);
+ FILE *edit_file = fopen(edit_file_name,"w");
+ char *profile_file_name = vtmpnam(NULL);
+ FILE *profile_file = fopen(profile_file_name,"w");
+ char *split_file_name = "AHA";//vtmpnam(NULL);
+ ++it_coutner_2;
+
+ gap_prf = enlarge_prof(gap_prf, alignment_length, alphabet_size);
+ no_gap_prf = enlarge_prof(no_gap_prf, alignment_length, alphabet_size );
+ no_gap_prf->number_of_sequences = 0;
+ gap_prf->number_of_sequences = 0;
+
+ split_set(alignment_file, gap_prf, no_gap_prf, seq1, column_index, param->char2pos, split_file_name);
+ gap_prf -> length = alignment_length;
+ no_gap_prf -> length = alignment_length;
+
+
+ gap_list_1 = del_gap_from_profile(gap_prf, alphabet_size, gap_list_1, gap_list_1_length, &num_gaps_1 );
+ gap_list_2 = del_gap_from_profile(no_gap_prf, alphabet_size, gap_list_2, gap_list_2_length, &num_gaps_2 );
+
+
+ fclose(alignment_file);
+ profiles[0] = gap_prf;
+ profiles[1] = no_gap_prf;
+
+ alignment_length = gotoh_dyn(profiles, param, method_arguments_p, 0, edit_file, profile_file, 0);
+ seq1 =vrealloc(seq1, (1+alignment_length)*sizeof(char));
+ seq2 =vrealloc(seq2, (1+alignment_length)*sizeof(char));
+
+ fclose(edit_file);
+ edit_file = fopen(edit_file_name,"r");
+ edit2seq_pattern(edit_file, seq1, seq2);
+ fclose(edit_file);
+ fclose(profile_file);
+
+ write_iterated_aln(aln_file_name, out_file_name, split_file_name, seq1, gap_list_1, num_gaps_1, seq2, gap_list_2, num_gaps_2);
+ aln_file_name = out_file_name;
+
+// if (it_coutner_2 == 2)
+// break;
+// }
+ char command[1000];
+ sprintf(command, "mv %s %s",out_file_name, out_file_name_end);
+ system(command);
+
+ vfree(seq1);
+ vfree(seq2);
+ vfree(gap_list_1);
+ vfree(gap_list_2);
+ if (!strcmp(param->method, "fast"))
+ {
+ free_sparse((Sparse_dynamic_param*)method_arguments_p);
+ }
+ else if (!strcmp(param->method, "nw"))
+ {
+ free_nw((Nw_param*)method_arguments_p, alphabet_size);
+ }
+ else if (!strcmp(param->method, "gotoh"))
+ {
+ free_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
+ }
+}
+
+
+
+
+void
+write_iterated_aln(char* old_aln_file_name, char* new_aln_file_name, char *gap_file_name, char *seq1, int *gap_list1, int num_gap1, char *seq2, int *gap_list2, int num_gap2)
+{
+ int i;
+// for (i = 0; i < num_gap1; ++i)
+// {
+// printf("%i ", gap_list1[i]);
+// }
+// printf("\n-----\n");
+// for (i = 0; i < num_gap2; ++i)
+// {
+// printf("%i ", gap_list2[i]);
+// }
+// printf("\n");
+//
+// printf("%s\n%s\n",seq1, seq2);
+// printf("EX: %i\n", gap_list1[0]);
+ FILE *gap_file = fopen(gap_file_name, "r");
+ FILE *new_aln_file = fopen(new_aln_file_name, "w");
+// if (new_aln_file == NULL)
+// printf("AHA\n");
+ FILE *old_aln_file = fopen(old_aln_file_name, "r");
+
+ const int GAPLINE_LENGTH = 5;
+ char gapline[GAPLINE_LENGTH];
+ char *seq;
+ int *gap_list;
+ int num_gap;
+
+ const int ALNLINE_LENGTH = 200;
+ char alnline[ALNLINE_LENGTH];
+
+ int alnline_pos = -1;
+ int overall_aln_pos = -1;
+ int pattern_pos = -1;
+ int gap_pos = 0;
+
+
+// fprintf(new_aln_file, "%s", alnline);
+ while(fgets(gapline, GAPLINE_LENGTH , gap_file)!=NULL)
+ {
+
+ //this is the sequence name
+ fgets(alnline, ALNLINE_LENGTH , old_aln_file);
+ fprintf(new_aln_file,"%s",alnline);
+ //getting (part if) the sequence
+ fgets(alnline, ALNLINE_LENGTH , old_aln_file);
+ if (gapline[0] == 'g')
+ {
+ seq = seq1;
+ gap_list = gap_list1;
+ num_gap = num_gap1;
+ }
+ else
+ {
+ seq = seq2;
+ gap_list = gap_list2;
+ num_gap = num_gap2;
+ }
+ overall_aln_pos = 0;
+ pattern_pos = 0;
+ gap_pos = 0;
+ alnline_pos = -1;
+
+ //go along the pattern
+ while (seq[pattern_pos] != '\0')
+ {
+ if (seq[pattern_pos] == '-')
+ {
+ ++pattern_pos;
+ fprintf(new_aln_file,"-");
+ continue;
+ }
+ //get next part of the sequence if you are at the end
+ if ((alnline[++alnline_pos] == '\n') || (alnline[alnline_pos] == '\0'))
+ {
+ alnline_pos = 0;
+ fgets(alnline, ALNLINE_LENGTH, old_aln_file);
+ }
+ if ((gap_pos < num_gap) && (overall_aln_pos == gap_list[gap_pos]))
+ {
+ ++gap_pos;
+ ++overall_aln_pos;
+ continue;
+ }
+ fprintf(new_aln_file,"%c",alnline[alnline_pos]);
+ ++overall_aln_pos;
+ ++pattern_pos;
+ }
+ fprintf(new_aln_file,"\n");
+ }
+
+ fclose(new_aln_file);
+ fclose(old_aln_file);
+ fclose(gap_file);
+}
+
+
+void
+edit2seq_pattern(FILE *edit_file, char *seq1, char *seq2)
+{
+
+ fseek(edit_file, 0, 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);
+
+ char x;
+ int number;
+ int pos = -1;
+ while(fgets(line, LINE_LENGTH , edit_file)!=NULL)
+ {
+// printf("%s\n",line);
+ x = line[0];
+ if (x == '*')
+ break;
+ number = atoi(&line[1]);
+ if (x == 'M')
+ {
+ while (--number >= 0)
+ {
+ seq1[++pos] = 'X';
+ seq2[pos] = 'X';
+ }
+ }
+ else if (x == 'I')
+ {
+ while (--number >= 0)
+ {
+ seq1[++pos] = 'X';
+ seq2[pos] = '-';
+ }
+ }
+ else if (x == 'D')
+ {
+ while (--number >= 0)
+ {
+// printf("POS: %i\n", pos);
+ seq1[++pos] = '-';
+ seq2[pos] = 'X';
+ }
+ }
+ }
+ seq1[++pos] = '\0';
+ seq2[pos] = '\0';
+// printf("%s\n%s\n",seq1, seq2);
+}
+/******************************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