+++ /dev/null
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <ctype.h>
-#include <string.h>
-#include <sys/stat.h>
-#include <dirent.h>
-
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "fastal_lib_header.h"
-
-
-
-
-/**
- * \brief Function to calculate the sum of pairs score with quasi affine gap costs.
- *
- * \param alignment_file File with the aligned sequences in fasta_format.
- * \param file_index Positions of the sequences in \a alignment_file.
- * \param number_of_sequences The number of sequences in \a alignment_file.
- * \param score_matrix Matrix containing the scores.
- * \param gop The gap opening costs.
- * \param gep The gep extension costs.
- * \return The score of the alignment.
- */
-double
-calculate_sum_of_pairs_score_affine(char *alignment_file_name,
- int **score_matrix,
- double gop,
- double gep)
-{
-
- const int LINE_LENGTH = 1000;
- int number_of_sequences = 0;
- FILE *alignment_file = fopen(alignment_file_name,"r");
- if (alignment_file == NULL)
- {
- printf("FILE COULD NOT BE OPENED!\n");
- exit(1);
- }
- int alignment_length = 0;
- char line[LINE_LENGTH+1];
- line[0] = '0';
-
- while (line[0] != '>')
- {
- fgets(line, LINE_LENGTH, alignment_file);
- }
- fgets(line, LINE_LENGTH, alignment_file);
- while (line[0] != '>')
- {
- alignment_length += strlen(line);
- fgets(line, LINE_LENGTH, alignment_file);
- }
-
- int alphabet_size = 28;
- int **counting = vcalloc(alphabet_size, sizeof(int*));
-
- int i;
- for (i = 0; i < alphabet_size; ++i)
- {
- counting[i] = vcalloc(alignment_length, sizeof(int));
- }
-
- alphabet_size -= 2;
- char c;
- int was_gap = 0;
- int pos_line, pos_alignment = -1;
- fseek (alignment_file, 0, SEEK_SET);
- pos_line = 0;
-
- int gap_line = 0;
- int current_size = 100;
- int gap_size = 0;
- int gap_pos = 0;
-
- int **gaps = vcalloc(100, sizeof(int*));
-
- while (fgets(line, LINE_LENGTH, alignment_file) != NULL)
- {
- if (line[0] != '>')
- {
- pos_line = -1;
- while (((c = line[++pos_line]) != '\n') && (c != '\0'))
- {
- if (isalpha(c))
- {
- ++counting[toupper(c)-'A'][++pos_alignment];
- if (was_gap)
- {
- gaps[gap_line][gap_pos++] = pos_alignment-1;
- }
- was_gap = 0;
- }
- else
- {
- ++counting[alphabet_size][++pos_alignment];
- if (!was_gap)
- {
- ++counting[alphabet_size+1][pos_alignment];
- was_gap = 1;
- if (gap_pos >= gap_size-4)
- {
- gap_size += 50;
- gaps[gap_line] = vrealloc(gaps[gap_line], gap_size * sizeof(int));
- gaps[gap_line][gap_pos++] = pos_alignment;
- }
- }
- }
- }
- }
- else
- {
- if (number_of_sequences != 0)
- {
- if ((gap_pos % 2) != 0)
- {
- gaps[gap_line][gap_pos++] = pos_alignment-1;
- }
- gaps[gap_line][gap_pos++] = alignment_length+2;
- gaps[gap_line][gap_pos] = alignment_length+2;
- ++gap_line;
- }
- if (current_size == gap_line)
- {
- current_size += 50;
- gaps = vrealloc(gaps, current_size * sizeof(int*));
- }
- gap_size = 2;
- gaps[gap_line] = vcalloc(gap_size, sizeof(int));
- ++number_of_sequences;
- pos_alignment = -1;
- was_gap = 0;
-
- gap_pos = 0;
- }
- }
-
- gaps[gap_line][gap_pos++] = alignment_length+2;
- gaps[gap_line][gap_pos] = alignment_length+2;
-
-
- long double score = 0;
- int j, k, tmp2;
- int non_gap;
- for (i = 0; i < alignment_length; ++i)
- {
- for (j = 0; j < alphabet_size; ++j)
- {
- if ((tmp2 = counting[j][i]) >1)
- {
- while (tmp2 > 1)
- {
- score += score_matrix[j][j] * (--tmp2);
- }
- }
- for (k = j+1; k < alphabet_size; ++k)
- {
- score += score_matrix[j][k] * counting[j][i] * counting[k][i];
- }
-
- }
-
-
- non_gap = number_of_sequences - counting[alphabet_size][i];
- score += counting[alphabet_size][i] * non_gap * gep;
-
- }
-
- ++gap_line;
- alignment_length += 2;
- unsigned long gap_open = 0;
-
- int chunk = 500;
-
-// #pragma omp parallel shared(gaps, chunk, alignment_length, gap_open) private(i, j)
- {
- int alignment_length2 = alignment_length;
-// #pragma omp for schedule(dynamic,chunk) reduction(+:gap_open) nowait
- for (i = 0; i < gap_line; ++i)
- {
- int *gaps1 = gaps[i];
- for (j = i+1; j < gap_line; ++j)
- {
- int *gaps2 = gaps[j];
- int k = 0;
- int l = 0;
- while ((gaps1[k] != alignment_length2) && (gaps2[l] != alignment_length2))
- {
- if (gaps1[k+1] < gaps2[l])
- {
- ++gap_open;
- k+=2;
- continue;
- }
- if (gaps2[l+1] < gaps1[k])
- {
- ++gap_open;
- l+=2;
- continue;
- }
- if (gaps1[k] < gaps2[l])
- {
- if (gaps1[k+1] < gaps2[l+1])
- {
- ++gap_open;
- k+=2;
- continue;
- }
- else
- {
- l+=2;
- continue;
- }
- }
- if (gaps1[k] > gaps2[l])
- {
- if (gaps1[k+1] > gaps2[l+1])
- {
- ++gap_open;
- l+=2;
- continue;
- }
- else
- {
- k+=2;
- continue;
- }
- }
- if (gaps1[k] == gaps2[l])
- {
- if (gaps1[k+1] == gaps2[l+1])
- {
- k+=2;
- l+=2;
- continue;
- }
- else
- {
- if (gaps1[k+1] <gaps2[l+1])
- {
- k+=2;
- }
- else
- {
- l+=2;
- }
- }
- }
-
- }
- while (gaps1[k] != alignment_length2)
- {
- ++gap_open;
- k+=2;
- }
- while (gaps2[l] != alignment_length2)
- {
- ++gap_open;
- l+=2;
- }
- }
- }
- }
-
- score += gop * gap_open;
-
-
- return score / number_of_sequences;
-}
-
-
-int
-compare (const void * a, const void * b)
-{
- return (*(int*)a - *(int*)b);
-}
-
-
-
-int
-get_random_value(int a, int b)
-{
- return rand() % (b - a + 1) + a;
-// j = 1 + (int)( 10.0 * rand() / ( RAND_MAX + 1.0 ) );
-}
-
-
-
-
-
-void
-compute_ref_alignments(char *seq_file_name, char* ref_directory, int num_alignments, int num_seq_in_aln)
-{
-
- printf("%s %s %i %i\n", seq_file_name, ref_directory, num_alignments, num_seq_in_aln);
- //check if target directory exists, if not make it
- char directory[500];
- if (ref_directory[0]== '~')
- {
- sprintf(directory,"%s%s",getenv("HOME"),ref_directory+1);
- }
- else
- {
- sprintf(directory,"%s",ref_directory);
- }
- struct stat st;
- if(stat(directory, &st) != 0)
- {
- mkdir(directory,0700);
- }
-
- printf("made\n");
-
- // get sequences & produce alignments
- //make index
- long *file_positions = NULL;
- long **tmp_pos = &file_positions;
- int number_of_sequences = make_index_of_file(seq_file_name, tmp_pos);
- char *tmp_file_name = vtmpnam(NULL);
- FILE *seq_file = fopen(seq_file_name, "r");
- int *already_taken = vcalloc(num_seq_in_aln, sizeof(int));
- int i, j, k, tmp, already, num_already;
- for (i = 0; i < num_alignments; ++i)
- {
- //choose randomly the sequences
-
- num_already = 0;
- while (num_already < num_seq_in_aln)
- {
- already = 0;
- tmp = get_random_value(0, number_of_sequences-1);
- for (k = 0; k < num_already; ++k)
- {
- if (already_taken[k] == tmp)
- {
- already = 1;
- break;
- }
- }
- if (!already)
- {
- already_taken[num_already] = tmp;
- ++num_already;
- }
- }
- //write sequences into temporary file
- const int LINE_LENGTH = 200;
- char line[LINE_LENGTH];
-
-
- FILE* tmp_file = fopen(tmp_file_name, "w");
- for (k = 0; k < num_already; ++k)
- {
- fseek(seq_file, file_positions[already_taken[k]], SEEK_SET);
- fgets(line, LINE_LENGTH, seq_file);
- fprintf(tmp_file, "%s", line);
- while( (fgets(line, LINE_LENGTH, seq_file)) != NULL)
- {
- if (line[0] == '>')
- break;
- fprintf(tmp_file, "%s", line);
- }
- }
- fclose(tmp_file);
-
- //calculate alignment
- char aln_command[1000];
- sprintf(aln_command,"t_coffee -in %s -outfile %s/%i.fa -output fasta_aln 2>/dev/null >/dev/null", tmp_file_name, directory, i);
- system(aln_command);
- }
-
- vfree(already_taken);
-}
-
-
-
-// void
-// compute agreement_score(char *alignment_file_name, char *alignment_directory)
-// {
-//
-// }
-
-
-
-
-/**
- * This algorithm choses according to a tree a number of sequenes of which it computes a reference algnment.
- * \param seq_file_name The sequence file.
- * \param tree_file_name The given tree.
- * \param ref_aln_name The file where the reference alignment should be written to.
- * \param num_seq_in_ref The number of sequences to choose.
-**/
-void
-make_ref_alignment(char *seq_file_name, char *tree_file_name, char *ref_aln_name, int num_seq_in_ref)
-{
- const int LINE_LENGTH = 200;
- char line[LINE_LENGTH];
-
- FILE *tree_f = fopen(tree_file_name,"r");
-
-
- long *file_positions = NULL;
- long **tmp = &file_positions;
- int num_sequences = make_index_of_file(seq_file_name, tmp);
-
-
- int every_x = num_sequences/(num_seq_in_ref);
-
- int *seq_ids = vcalloc(num_seq_in_ref,sizeof(int));
- char delims[] = " ";
-
- int pos = -1;
- fseek (tree_f, 0, SEEK_SET);
- int node;
- int dist = every_x;
- while(fgets(line, LINE_LENGTH, tree_f)!=NULL)
- {
- node = atoi(strtok(line,delims));
- if (node < num_sequences)
- {
- if (dist == every_x)
- {
-
- seq_ids[++pos] = node;
- dist = 0;
- }
- else
- {
- ++dist;
- }
- }
-
- node = atoi(strtok(NULL,delims));
- if (node < num_sequences)
- {
- if (dist == every_x)
- {
- seq_ids[++pos] = node;
- dist = 0;
- }
- else
- {
- ++dist;
- }
- }
- }
- fclose(tree_f);
-
-
- int i = 0;
- qsort(seq_ids, num_seq_in_ref, sizeof(int), compare);
-
- char *seq_ref_name = vtmpnam(NULL);
- FILE *seq_ref_f = fopen(seq_ref_name, "w");;
- FILE *seq_f = fopen(seq_file_name, "r");
-
- for (i = 0; i < num_seq_in_ref; ++i)
- {
- fseek (seq_f, file_positions[seq_ids[i]], SEEK_SET);
- fgets(line, LINE_LENGTH, seq_f);
- fprintf(seq_ref_f, "%s", line);
- while(fgets(line, LINE_LENGTH, seq_f)!=NULL)
- {
- if (line[0] == '>')
- break;
- fprintf(seq_ref_f, "%s", line);
- }
- }
- fclose(seq_ref_f);
- fclose(seq_f);
-
- char command[500];
- sprintf(command, "mafft --quiet %s > %s", seq_ref_name, ref_aln_name);
- system(command);
-}
-
-
-/**
- * This function calculates the agreement between two alignments in a specified number of sequences.
- * \param ref_file_name The reference alignment.
- * \param aln_file_name The test alignment.
- **/
-double
-agreement_score(char *ref_file_name, char *aln_file_name)
-{
- const int LINE_LENGTH = 200;
- char line[LINE_LENGTH];
-
- int number_of_sequences = 0;
- FILE *ref_f = fopen(ref_file_name,"r");
- while(fgets(line, LINE_LENGTH, ref_f)!=NULL)
- {
- if (line[0] == '>')
- {
- ++number_of_sequences;
- }
- }
- fseek(ref_f, 0, SEEK_SET);
- char **seq_names = vcalloc(number_of_sequences, sizeof(char*));
- int i = 0;
- while(fgets(line, LINE_LENGTH, ref_f)!=NULL)
- {
- if (line[0] == '>')
- {
- seq_names[i] = vcalloc(LINE_LENGTH, sizeof(char));
- sprintf(seq_names[i], "%s",line);
- ++i;
- }
- }
-
-
- fclose(ref_f);
- FILE *aln_f = fopen(aln_file_name,"r");
- char *tmp_name = vtmpnam(NULL);
- FILE *tmp_f = fopen(tmp_name,"w");
- int x = 0;
- long last_pos = -1;
- while(fgets(line, LINE_LENGTH, aln_f)!=NULL)
- {
- if (line[0] == '>')
- {
- for (i = 0; i < number_of_sequences; ++i)
- {
- if (!strcmp(line, seq_names[i]))
- {
- fprintf(tmp_f, "%s", line);
- last_pos = ftell(aln_f);
- while(fgets(line, LINE_LENGTH, aln_f)!=NULL)
- {
- if (line[0] == '>')
- break;
- fprintf(tmp_f, "%s", line);
-
- last_pos = ftell(aln_f);
- }
- ++x;
- fseek(aln_f, last_pos, SEEK_SET);
- break;
- }
- }
- }
- if (x == number_of_sequences)
- break;
- }
- fclose(aln_f);
- fclose(tmp_f);
- char command[500];
-// sprintf(command, "cp %s /users/cn/ckemena/Desktop/1.fa", tmp_name);
-// system(command);
-// fp = popen("ls -l", "r");
- sprintf(command, "t_coffee -other_pg aln_compare -al1 %s -al2 %s ", ref_file_name, tmp_name);
-
-
- FILE *result;
- result = popen(command, "r");
- fgets(line, LINE_LENGTH, result);
- fgets(line, LINE_LENGTH, result);
- fgets(line, LINE_LENGTH, result);
-
-// sprintf(command, "cp %s ~/Destkop/1.fa", tmp_name);
-// system(command);
-
- char delims[] = " ";
- char *tmp_str;
- strtok(line, delims);
- for (i = 0; i < 3; ++i)
- {
- while((tmp_str = strtok(NULL, delims))==NULL);
- }
- return(atof(tmp_str));
-}
-
-
-complete_agreement_score(char *aln_file_name, const char *ref_directory)
-{
- struct dirent *dp;
- char *name = strrchr(aln_file_name,'/')+1;
-// printf("%s ",name);
- // enter existing path to directory below
- char directory[500];
- if (ref_directory[0]== '~')
- {
- sprintf(directory,"%s%s",getenv("HOME"),ref_directory+1);
- }
- else
- {
- sprintf(directory,"%s",ref_directory);
- }
- DIR *dir = opendir(directory);
-
- char ref_file_name[200];
- while ((dp=readdir(dir)) != NULL)
- {
- if ((strcmp(dp->d_name,".")) && (strcmp(dp->d_name,"..")))
- {
-
- sprintf(ref_file_name, "%s/%s",directory, dp->d_name);
- printf("%s %f\n",name, agreement_score(ref_file_name, aln_file_name));
-// printf("%f ", agreement_score(ref_file_name, aln_file_name));
-
- }
- }
- printf("\n");
- closedir(dir);
- return 0;
-
-}
-
-
-// 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));
-/******************************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