+++ /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