+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdarg.h>
-#include <string.h>
-#include <ctype.h>
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "define_header.h"
-#include "dp_lib_header.h"
-
-
-static double mc_delta_matrix ( int ***mat1, int ***mat2, char **alp, int nch);
-static double delta_matrix ( int **mat1,int **mat2, char *alp);
-static double ***pavie_seq2pavie_fmat (Sequence *S,double *gop, double *gep, char **mat, char *idmat, int id_threshold, int sample_size, int nch, char *param );
-static int **pavie_fmat2pavie_logodd_mat (double **fmat, char *alp);
-static double **pavie_aln2fmat(Alignment *A, double **fmat, char *idmat, int id_threshold, int ch, int nch, char *param);
-static int pavie_mat2pavie_id_mat ( int **mat,char *in_name, char *alp, char *ignore, char *force,int T, char *out_name);
-static double paviemat2gep ( int **mat, char *alp);
-static Alignment *align_pavie_sequences (char *seq0,char *seq1,char **mat,double *gop,double *gep,int nch, char *param);
-static int pavie_score (char *s0,int p0, char *s1,int p1,char **mat_file, double *gop, double *gep, int nch, float factor, char *param);
-static char **seq2pavie_alp (Sequence *S, int nch);
-static Sequence * seq2pavie_seq ( Sequence *S, int nch);
-static FILE* output_pavie_aln (Alignment *A, int nch, FILE *fp);
-static char **output_pavie_mat_list ( int ***current_mat, double *gep, char **alp, int nch,char *prefix,int cycle, char **mat_name);
-static float pavie_aln2id ( Alignment *A, int mode);
-static int check_pavie_cl ( char *string);
-float pavie_aln2delta_age ( Alignment *A,int s0, int s1, int a0, int a1);
-
-static float tgep_factor;
-static int id_thres_used_aln;
-static int log_odd_mode;
-Sequence * pavie_seq2noisy_seq ( Sequence *S, int freq, char *alp)
-{
- int a, b, l1, l2;
-
- vsrand(0);
-
- if (alp==NULL)
- {
- char **x;
- x=seq2pavie_alp (S,1);
- alp=x[0];
- }
-
- l2=strlen (alp);
- for (a=0; a< S->nseq; a++)
- {
- l1=strlen (S->seq[a]);
- for ( b=0; b<l1; b++)
- {
- if ( (rand ()%100+1)<freq)
-
- S->seq[a][b]=alp[rand()%l2];
- }
- }
- return S;
-}
-Sequence * pavie_seq2random_seq ( Sequence *S, char *subst)
-{
- int a, b, r, l;
-
-
- vsrand (0);
- r=subst[0]; subst++;
- l=strlen (subst);
- for ( a=0; a< S->nseq; a++)
- for (b=0; b<S->len[a]; b++)
- if ( S->seq[a][b]==r)S->seq[a][b]=subst[rand()%l];
- return S;
-}
-
-double **pavie_seq2pavie_aln(Sequence *S,char *mat, char *mode)
-{
- int a, b,c, nch=0;
- char **mat_list;
- char *buf;
-
- double *gep, *gop;
- Alignment *A;
- char **alp;
- char *pavie_idmat;
- FILE *fp;
- double **dist_mat;
- float score;
-
- check_pavie_cl (mode);
-
- mat_list=declare_char (100, 100);
-
- if ( is_matrix (mat))
- {
- sprintf ( mat_list[nch++], "%s", mat);
- }
- else
- {
- fp=vfopen (mat,"r");
- while ( (c=fgetc(fp))!=EOF)
- {
- ungetc(c, fp);
- fscanf (fp, "%s\n",mat_list[nch++]);
- }
- vfclose (fp);
- }
-
- alp=seq2pavie_alp (S, nch);
- S=seq2pavie_seq (S, nch);
-
- gop=vcalloc (nch, sizeof (double));
- gep=vcalloc (nch, sizeof (double));
-
- for ( a=0; a< nch; a++)
- {
- int **m;
- char *st;
- int v;
- m=read_matrice (mat_list[a]);
- if ((st=vstrstr(mode, "_GEP")))
- {
- sscanf ( st, "_GEP%d_", &v);
- gep[a]=v*-1;
- }
- else if ( m[0][GAP_CODE]==0)
- {
- gep[a]=paviemat2gep(m,alp[a]);
- }
- else
- {
- gep[a]=m[0][GAP_CODE];
- }
- free_int (m, -1);
- }
-
-
- if ( (buf=vstrstr (mode, "_TGEPF")))
- {
-
- sscanf (buf, "_TGEPF%f_", &tgep_factor);
- tgep_factor/=(float)100;
- }
- else
- {
- tgep_factor=0.5;
- }
-
- pavie_idmat=vtmpnam(NULL);
-
- pavie_mat2pavie_id_mat (NULL,"idmat", alp[0],"X","",1,pavie_idmat);
- dist_mat=declare_double ( S->nseq, S->nseq);
-
-
-
- for ( a=0; a< S->nseq-1; a++)
- {
- for ( b=a+1; b< S->nseq; b++)
- {
- int a0, a1;
- float delta_a;
-
-
- if ( ! strstr (mode, "_MSA_"))
- {
- A=align_pavie_sequences (S->seq[a],S->seq[b],mat_list,gop,gep,nch, mode);
- sprintf ( A->name[0], "%s", S->name[a]);
- sprintf ( A->name[1], "%s", S->name[b]);
- }
- else
- {
-
- A=strings2aln ( 2, S->name[a], S->seq[a], S->name[b], S->seq[b]);
- sprintf ( A->seq_al[0], "%s", S->seq[a]);
- sprintf ( A->seq_al[1], "%s", S->seq[b]);
- A->len_aln=strlen (S->seq[a]);
- ungap_aln (A);
- }
-
- if (strm (mode, "_ID01_"))A->score=score=pavie_aln2id (A, 1);
- else if ( vstrstr (mode, "_ID02_"))A->score=score=pavie_aln2id (A, 2);
- else if ( vstrstr (mode, "_ID04_"))A->score=score=pavie_aln2id (A, 4);
- else if ( vstrstr (mode, "_ID05_"))A->score=score=pavie_aln2id (A, 5);
- else if ( vstrstr (mode, "_ID06_"))A->score=score=pavie_aln2id (A, 6);
-
- else A->score=score=pavie_aln2id (A, 1);
-
- a0=S->seq[a][strlen(S->seq[a])+1];
- a1=S->seq[b][strlen(S->seq[b])+1];
-
- delta_a=pavie_aln2delta_age (A, 0, 1, a0, a1);
-
- if ( vstrstr (mode, "_MATDIST_"))
- dist_mat[a][b]=dist_mat[b][a]=(double)(vstrstr (mode, "_ID05") || vstrstr (mode, "_ID06"))?-score*100:(100-score);
- else if ( vstrstr (mode, "_MATSIM_"))
- dist_mat[a][b]=dist_mat[b][a]=(double)(score);
-
-
- if ( !vstrstr (mode, "_MAT") )
- {
- fprintf ( stdout, "#############\nAlignment %s %s: %d %% ID SCORE %d DELTA_AGE %.2f\n", S->name[a], S->name[b], A->score, A->score_aln, delta_a);
- output_pavie_aln (A,nch, stdout);
- }
- free_aln(A);
- }
- }
-
- if ( vstrstr (mode, "_MAT") && !vstrstr ( mode, "_NOPRINT_"))
- {
- if ( vstrstr (mode, "_MFORMAT2"))
- {
- int max, n;
- float *tot,s, bigtot=0;
-
- for (max=0, a=0; a< S->nseq; a++)max=MAX(max,(strlen (S->name[a])));
- tot=vcalloc ( S->nseq, sizeof (float));
- fprintf (stdout, "# TC_DISTANCE_MATRIX_FORMAT_01\n");
- for ( a=0; a<S->nseq; a++)
- fprintf ( stdout, "# SEQ_INDEX %s %d\n",S->name[a],a);
- fprintf ( stdout, "# PW_SEQ_DISTANCES \n");
- for (n=0,a=0;a< S->nseq-1; a++)
- {
- for ( b=a+1; b<S->nseq; b++, n++)
- {
- s=dist_mat[a][b];
-
- fprintf (stdout, "BOT\t %4d %4d\t %5.2f %*s\t %*s\t %5.2f\n", a,b,s,max,S->name[a], max, S->name[b], s);
- fprintf (stdout, "TOP\t %4d %4d\t %5.2f %*s\t %*s\t %5.2f\n", b,a,s,max,S->name[b], max, S->name[a], s);
- tot[a]+=s;
- tot[b]+=s;
- bigtot+=s;
- }
- }
- for ( a=0; a< S->nseq; a++)
- {
- fprintf (stdout, "AVG\t %d\t %*s\t %*s\t %5.2f\n", a,max,S->name[a], max, "*", tot[a]/(S->nseq-1));
- }
- vfree (tot);
- fprintf (stdout, "TOT\t %*s\t %*s\t %5.2f\n", max,"TOT", max, "*", bigtot/n);
- vfclose (stdout);
- }
- else
- {
- for ( a=0; a<S->nseq; a++)
- {
- fprintf ( stdout, "\n%s ", S->name[a]);
- for ( b=0; b< S->nseq; b++)
- fprintf ( stdout, "%6d ", (int)(dist_mat[a][b]*100));
- }
- }
- }
-
- return dist_mat;
-}
-
-float pavie_aln2delta_age ( Alignment *A,int s0, int s1, int a0, int a1)
-{
- int a,r0, r1, g0, g1, n;
- float delta;
- for (n=0, delta=0, a=0; a< A->len_aln; a++)
- {
- r0=A->seq_al[s0][a];
- r1=A->seq_al[s1][a];
-
- g0=!is_gap(r0);
- g1=!is_gap(r1);
-
- a0+=g0;a1+=g1;
- if ( g0 && g1)
- {
- delta+=FABS((a0-a1));
- n++;
- }
- }
- delta/=(float)((n)?n:1);
- return delta;
-}
-
-int **pavie_seq2trained_pavie_mat(Sequence *S, char *param)
-{
- double ***fmat;
- int ***current_mat;
- int ***previous_mat;
- char **alp;
-
- char **mat_file;
- double d,delta_min=10;
- double *gep;
- double *gop;
-
- char ignore[100];
- char force [100];
- char pavie_idmat[100];
- int id_threshold;
- int sample_size;
- char *b;
- int a,n=0,nch=1;
- char *buf;
-
- check_pavie_cl (param);
-
- if ( !param)param=vcalloc (1, sizeof (char));
-
- if ((b=vstrstr(param,"_THR")))sscanf ( b, "_THR%d_", &id_threshold);
- else id_threshold=0;
-
- sample_size=0;
- if ((b=vstrstr(param,"_SAMPLE")))sscanf ( b, "_SAMPLE%d_", &sample_size);
- if ((b=vstrstr(param,"_PARALOGOUS")))
- {
- sscanf ( b, "_PARALOGOUS%d_", &sample_size);
- sample_size*=-1;
- }
-
- if ((b=vstrstr(param,"_CHANNEL")))sscanf ( b, "_CHANNEL%d_", &nch);
- else nch=1;
-
- if ( (buf=vstrstr (param, "_TGEPF")))
- {
- sscanf (buf, "_TGEPF%f_", &tgep_factor);
- tgep_factor/=(float)100;
- }
- else
- {
- tgep_factor=0.5;
- }
- if ( (buf=vstrstr (param, "_PAMLOGODD_")))
- {
- log_odd_mode=1;
- }
- /*Declare Arrays*/
- gep=vcalloc (nch, sizeof (double));
- gop=vcalloc (nch, sizeof (double));
- mat_file=declare_char ( nch, 100);
- current_mat =vcalloc ( nch, sizeof (double**));
- previous_mat=vcalloc ( nch, sizeof (double**));
-
-
- sprintf (ignore, "X");
- force[0]='\0';
- sprintf ( pavie_idmat, "pavie_idmat");
-
-
-
- alp=seq2pavie_alp (S, nch);
-
- S=seq2pavie_seq (S, nch);
-
- pavie_mat2pavie_id_mat (NULL,"idmat", alp[0],ignore,force,1,pavie_idmat);
-
- for ( a=0; a<nch; a++)sprintf (mat_file[a], "idmat");
-
-
- fmat=pavie_seq2pavie_fmat ( S,gop,gep,mat_file,pavie_idmat, id_threshold, sample_size, nch, param);
-
-
- for (a=0; a<nch; a++)
- {
- current_mat[a]=pavie_fmat2pavie_logodd_mat(fmat[a], alp[a]);
- gep[a]=paviemat2gep(current_mat[a], alp[a]);
- }
- free_arrayN((void*)fmat, 3);
-
- mat_file=output_pavie_mat_list ( current_mat,gep, alp, nch,"", n++, mat_file);
-
-
-
- fprintf ( stdout, "\n");
- previous_mat=vcalloc ( nch, sizeof (int**));
- while ((d=mc_delta_matrix (previous_mat, current_mat, alp, nch))>delta_min)
- {
-
- fprintf ( stdout, "\nDelta=%d: ",(int) d);
- for (a=0; a<nch; a++)
- {
- free_int (previous_mat[a], -1);
- previous_mat[a]=current_mat[a];
- }
- fprintf ( stdout, "\n");
-
- fmat=pavie_seq2pavie_fmat (S,gop,gep,mat_file, pavie_idmat, id_threshold, sample_size, nch, param);
-
-
- for (a=0; a< nch; a++)
- {
- current_mat[a]=pavie_fmat2pavie_logodd_mat(fmat[a], alp[a]);
- gep[a]=paviemat2gep(current_mat[a], alp[a]);
- }
-
- mat_file=output_pavie_mat_list ( current_mat,gep, alp, nch,"", n, mat_file);
- free_arrayN((void*)fmat, 3);
- n++;
- }
-
- fprintf ( stdout, "\nDelta=%d Mat: ",(int) d);
- for (a=0; a<nch; a++)
- {
- fprintf ( stdout, "%s ",mat_file[a]);
- }
- fprintf ( stdout, "\n");
-
- return current_mat[0];
-}
-
-double mc_delta_matrix ( int ***mat1, int ***mat2, char **alp, int nch)
-{
- int a;
- double delta=0;
- if ( !mat1 || !mat2) return 100000;
- for ( a=0; a< nch; a++)
- delta+=delta_matrix (mat1[a], mat2[a], alp[a]);
- return delta/nch;
-}
-
-double delta_matrix ( int **mat1,int **mat2, char *alp)
-{
- int ns;
- double delta, v;
- int a, b;
-
- if ( mat1==NULL || mat2==NULL) return 100000;
-
- ns=strlen (alp);
- for (delta=0, a=0; a< ns; a++)
- for ( b=0; b< ns; b++)
- {
- v=mat1[alp[a]-'A'][alp[b]-'A']-mat2[alp[a]-'A'][alp[b]-'A'];
- delta+=v*v;
- }
- delta=sqrt(delta);
-
- return delta;
-}
-
-double ***pavie_seq2pavie_fmat (Sequence *S,double *gop, double *gep, char **mat, char *idmat, int id_threshold, int sample_size, int nch, char *param )
-{
- int a=0, b, chan;
- double ***fmat=NULL;
- Alignment *A;
- int exclude_id=1;
- static int tot=0;
-
-
- id_thres_used_aln=0;
- fmat=vcalloc ( nch, sizeof (double **));
-
- if (sample_size==0)
- {
-
- for (tot=0, a=0; a< S->nseq-exclude_id; a++)
- {
-
- output_completion ( stderr,a+1,S->nseq,1, "");
-
- for ( b=a+exclude_id; b< S->nseq; b++)
- {
- tot++;
-
- A=align_pavie_sequences (S->seq[a],S->seq[b],mat,gop,gep,nch,param);
-
- for ( chan=0; chan< nch; chan++)
- {
-
- fmat[chan]=pavie_aln2fmat (A, fmat[chan], idmat, id_threshold, chan, nch, param);
- }
- free_aln (A);
- }
- }
- }
- else
- {
- int c;
- static int **list;
-
- if ( sample_size>0 && !list)
- {
- if ( exclude_id==0)sample_size*=3;
- if (!list)
- {
- list=declare_int ((sample_size+1), 2);
- vsrand(0);
- tot=0;
- while (tot<sample_size)
- {
- a=rand()%(S->nseq);b=rand()%(S->nseq);
- if ( a!=b)
- {
- list[tot][0]=a;list[tot][1]=b;
- tot++;
- if ( exclude_id==0)
- {
- list[tot][0]=a;list[tot][1]=a;
- tot++;
- list[tot][0]=b;list[tot][1]=b;
- tot++;
- }
- }
- }
- }
- }
- else if ( sample_size<0 && !list)
- {
-
- int **sim;
- int m;
- sim=seq2sim_mat (S, "idmat");
- sample_size*=-1;
- list=declare_int (S->nseq*S->nseq, 2);
-
- m=S->nseq-exclude_id;
- for (a=0; a<m; a++)
- for ( b=a+exclude_id; b<S->nseq; b++)
- {
- if ( sim[a][b]>sample_size)
- {
- list[tot][0]=a;
- list[tot][1]=b;
- tot++;
- fprintf ( stderr, "\n%s %s: %d", S->name[a], S->name[b], sim[a][b]);
- fprintf ( stderr, "\nKeep %s Vs %s : %d%% ID", S->name[a], S->name[b], sim[a][b]);
- }
- }
- free_int(sim, -1);
- }
-
- for (c=0; c<tot; c++)
- {
- a=list[c][0];b=list[c][1];
- A=align_pavie_sequences (S->seq[a],S->seq[b],mat,gop,gep,nch, param);
- for (chan=0; chan< nch; chan++)
- fmat[chan]=pavie_aln2fmat (A, fmat[chan], idmat, id_threshold,chan, nch, param);
-
- free_aln (A);
- output_completion ( stderr,c,tot,1, "");
- }
- }
- fprintf ( stderr, "\n\tSample_size: %d Used alignments: %d\n", tot, id_thres_used_aln);
- return fmat;
-}
-
-
-
-int **pavie_fmat2pavie_logodd_mat (double **fmat, char *alp)
-{
- int s1, s2,S1, S2;
- double r1, r2;
- int **mat;
- int a, b;
- int ns;
- int logodd=0;
-
-
- ns=strlen (alp);
- mat=declare_int (256, 256);
-
- for ( a=0; a<ns; a++)
- {
- s1=tolower(alp[a]);
- fprintf ( stderr, "\n\tSymbol %c Freq %5.2f %%", s1, ((float)fmat[s1][0]*100)/(float)fmat[0][0]);
- }
-
-
-
- for (a=0; a<ns; a++)
- for (b=0; b<ns; b++)
- {
- s1=tolower(alp[a]);S1=toupper(alp[a]);
- s2=tolower(alp[b]);S2=toupper(alp[b]);
-
- if ( log_odd_mode==0)
- {
- r1=(fmat[s1][s2]+1)/(fmat[s1][s1]+1);
- r2=(fmat[s2][s1]+1)/(fmat[s2][s2]+1);
- logodd=(int)(((log(r1)+log(r2))/2)*PAVIE_MAT_FACTOR);
- }
- else if ( log_odd_mode==1)
- {
-
- r1=(fmat[s1][s2]+fmat[s2][s1]+1)/(fmat[1][1]+1);
- r2=((fmat[s1][1]+1)/(fmat[1][1]+1))*((fmat[s2][1]+1)/(fmat[1][1]+1))*2;
- logodd=(int)(log(r1/r2)*PAVIE_MAT_FACTOR);
- }
-
- mat[s1-'A'][s2-'A']=logodd;
- mat[S1-'A'][S2-'A']=logodd;
- mat[S1-'A'][s2-'A']=logodd;
- mat[s1-'A'][S2-'A']=logodd;
-
- }
- return mat;
-}
-
-double **pavie_aln2fmat(Alignment *A, double **fmat, char *idmat, int id_threshold, int ch, int nch, char *param)
-{
- int a;
- int c1, c2;
- int w, id;
- int l,start;
-
- l=(A->len_aln/nch);
- start=l*ch;
- A->len_aln=l;
- A->seq_al[0]+=start;
- A->seq_al[1]+=start;
-
-
-
- if ( fmat==NULL)fmat=declare_double(300, 300);
-
- if ( vstrstr (param, "_TWE00_"))w=100;
- else if ( vstrstr (param, "_TWE01_"))w=pavie_aln2id (A, 1);
- else if ( vstrstr (param, "_TWE02_"))w=pavie_aln2id (A, 2);
- else if ( vstrstr (param, "_TWE03_"))w=pavie_aln2id (A, 3);
- else if ( vstrstr (param, "_TWE04_"))w=pavie_aln2id (A, 4);
- else if ( vstrstr (param, "_TWE05_"))w=pavie_aln2id (A, 5);
- else if ( vstrstr (param, "_TWE06_"))w=pavie_aln2id (A, 6);
-
- else w=pavie_aln2id (A, 3);
-
- id=pavie_aln2id(A, 3);
-
-
- if (id<id_threshold)
- {
- A->len_aln*=nch;
- A->seq_al[0]-=start;A->seq_al[1]-=start;
- return fmat;
- }
- else
- {
- id_thres_used_aln++;
- for ( a=0; a<A->len_aln; a++)
- {
- c1=tolower(A->seq_al[0][a]);
- c2=tolower(A->seq_al[1][a]);
- fmat[c1][c2]+=w;
-
- fmat[c1][0]++;
- fmat[c1][1]+=w;
-
- fmat[c2][0]++;
- fmat[c1][1]+=w;
-
- fmat[0][0]+=2;
- fmat[1][1]+=2*w;
- }
- A->len_aln*=nch;
- A->seq_al[0]-=start;A->seq_al[1]-=start;
-
- return fmat;
- }
-}
-
-int pavie_mat2pavie_id_mat ( int **mat,char *in_name, char *alp, char *ignore, char *force,int T, char *out_name)
-{
- int n1, n2, n3;
- int s1, s2, S1, S2;
- int a, b;
- int **idmat;
-
- if (mat==NULL && in_name==NULL) return 0;
- else if (mat==NULL)
- {
- mat=read_matrice (in_name);
- }
-
-
- idmat=declare_int ( 256, 256);
- n1=strlen (alp);
- n2=strlen (ignore);
- n3=strlen (force);
-
- for (a=0; a< n1; a++)
- for ( b=0; b<n1; b++)
- {
- s1=tolower(alp[a])-'A';S1=toupper(alp[a])-'A';
- s2=tolower(alp[b])-'A';S2=toupper(alp[b])-'A';
- idmat[s1][s2]=idmat[s1][S2]=idmat[S1][S2]=idmat[S1][s2]=(mat[s1][s2]>=T)?PAVIE_MAT_FACTOR:0;
- }
- for (a=0; a<n3; a++)
- for (b=0; b<n1; b++)
- {
- s1=tolower(force[a])-'A';S1=toupper(force[a])-'A';
- s2=tolower(alp[b])-'A';S2=toupper(alp[b])-'A';
- idmat[s1][s2]=idmat[s1][S2]=idmat[S1][S2]=idmat[S1][s2]=PAVIE_MAT_FACTOR;
- }
-
- for (a=0; a<n2; a++)
- for (b=0; b<n1; b++)
- {
- s1=tolower(ignore[a])-'A';S1=toupper(ignore[a])-'A';
- s2=tolower(alp[b])-'A';S2=toupper(alp[b])-'A';
- idmat[s1][s2]=idmat[s1][S2]=idmat[S1][S2]=idmat[S1][s2]=0;
- }
-
-
-
-
- output_pavie_mat (idmat, out_name, 0, alp);
- free_int (idmat, -1);
- return 1;
-}
-double paviemat2gep ( int **mat, char *alp)
-{
- int a, b, l;
- int n=0;
- double gep=0;
- l=strlen ( alp);
-
- for (a=0; a<l-1; a++)
- for ( b=a+1; b< l; b++)
- {
- gep+=mat[alp[a]-'A'][alp[b]-'A'];
- n++;
- }
- gep/=n;
-
- return gep;
-
-}
-
-Alignment *align_pavie_sequences (char *seq0,char *seq1,char **mat,double *gop,double *gep,int nch, char *param)
-{
- double **F;
- int **T;
- Alignment *A;
- int XL, YL, len;
- int i, j, a, b, c;
- double match, gap_inX, gap_inY, MXY=1, GX=2, GY=3;
-
- char *ax, *ay;
- char *bufx, *bufy, *buf;
- char *x,*y;
- float factor;
-
- factor=tgep_factor;
-
-
- /*FActor
- terminal gaps are set to gep=gep*factor
- */
-
- if ( strm (seq0, seq1))
- {
- A=declare_aln2 (2,strlen(seq0)+1);
- A->len_aln=strlen (seq0);
- A->nseq=2;
- A->score=A->score_aln=100;
- sprintf ( A->seq_al[0], "%s", seq1);
- sprintf ( A->seq_al[1], "%s", seq0);
- return A;
- }
-
-
- x=seq0;
- y=seq1;
-
- XL=strlen (x)/nch;
- YL=strlen (y)/nch;
-
-
- ax=vcalloc ( (YL+XL)*nch+1, sizeof (char));
- ay=vcalloc ( (YL+XL)*nch+1, sizeof (char));
- bufx=vcalloc ( (YL+XL)*nch+1, sizeof (char));
- bufy=vcalloc ( (YL+XL)*nch+1, sizeof (char));
-
- F=declare_double (XL+2, YL+2);
- T=declare_int (XL+2, YL+2);
-
-
- /*Fill stage*/
- F[0][0] = 0;
- for(i = 1; i <=XL; i++)
- {
-
- F[i][0] = F[i-1][0]+pavie_score (x,i-1,NULL,GAP_CODE,mat, gop, gep, nch, factor, param) /*CL->M[x[i-1]-'A'][gap]*/;
-
- T[i][0] = GY;
- }
-
- for(j = 1; j <= YL; j++)
- {
-
- F[0][j] = F[0][j-1]+pavie_score (NULL,GAP_CODE,y,j-1,mat, gop, gep, nch, factor, param)/*CL->M[y[j-1]-'A'][gap]*/;
- T[0][j] = GX;
- }
-
-
- for(i = 1; i <= XL; i++)
- {
- for(j = 1; j <= YL; j++)
- {
-
- match = F[i-1][j-1] + /*CL->M[x[i-1]-'A'][y[j-1]-'A']*/pavie_score (x,i-1,y, j-1,mat, gop, gep, nch, 1, param);
- gap_inY= F[i-1][j] + /*CL->M[x[i-1]-'A'][gap]*/ pavie_score (x,i-1, NULL,GAP_CODE,mat, gop, gep, nch, (j==YL)?factor:1, param);
- gap_inX= F[i][j-1] + /*+ CL->M[y[j-1]-'A'][gap]*/ pavie_score (NULL,GAP_CODE,y, j-1,mat, gop, gep, nch, (i==XL)?factor:1, param);
-
- if ( match >= gap_inY && match >=gap_inX){F[i][j]=match; T[i][j]=MXY;}
- else if ( gap_inX>=gap_inY){F[i][j]=gap_inX; T[i][j]=GX;}
- else {F[i][j]=gap_inY; T[i][j]=GY;}
- }
- }
- /*Trace back stage*/
-
-
- i = XL;
- j = YL;
- len=0;
- while(!(i==0 && j==0))
- {
-
- if (T[i][j]==MXY)
- {
- ax[len]=1;i--;
- ay[len]=1;j--;
- }
- else if ( T[i][j]==GY)
- {
- ax[len]=1;i--;
- ay[len]='-';
- }
- else if ( T[i][j]==GX)
- {
- ax[len]='-';
- ay[len]=1;j--;
- }
- len++;
- }
-
- for (a=0; a<len; a++)
- for (b=1; b<nch; b++)
- {
- ax[a+len*b]=ax[len];
- ay[a+len*b]=ay[len];
- }
- len=len*nch;
- ax[len]='\0';
- ay[len]='\0';
-
-
- sprintf ( bufx, "%s", ax);
- sprintf ( bufy, "%s", ay);
-
- for (a=1;a<nch; a++)
- {
- strcat (ax, bufx);
- strcat (ay, bufy);
- }
-
- buf=ax;ax=invert_string (ax);vfree(buf);
- buf=ay;ay=invert_string (ay);vfree(buf);
-
-
- A=declare_aln2 (2,strlen(ax)+1);
-
-
- A->len_aln=strlen (ax);
- A->nseq=2;
- A->score=A->score_aln=F[XL][YL];
-
- for (a=0, b=0, c=0; a<A->len_aln; a++)
- {
- if (ax[a]==1)ax[a]=seq0[b++];
- if (ay[a]==1)ay[a]=seq1[c++];
- }
-
-
-
- sprintf ( A->seq_al[0], "%s", ax);
- sprintf ( A->seq_al[1], "%s", ay);
-
- vfree (ax); vfree(ay);vfree (bufx); vfree (bufy);free_double(F, -1); free_int (T, -1);
- return A;
-}
-
-
-int pavie_score (char *s0,int p0, char *s1,int p1,char **mat_file, double *gop, double *gep, int nch, float factor, char *param)
-
- {
- static char *cmat;
- static int ***mat;
- static int use_age;
- static int mchscore=-1;
-
- int l0, l1, c0, c1;
- int a, score=0;
-
- if ( !use_age)
- {
- strget_param ( param, "_AGECHANNEL", "-1", "%d", &use_age);
-
- }
- if (mchscore==-1)
- {
- strget_param (param, "_MCHSCORE", "0", "%d", &mchscore);
-
- }
-
- if ( !cmat || !mat_file || !strm (cmat, mat_file[0]))
- {
- if ( !cmat)cmat=vcalloc ( 100, sizeof (char));
- sprintf ( cmat, "%s", (mat_file)?mat_file[0]:"idmat");
- if ( !mat)mat=vcalloc ( nch, sizeof (int**));
- for ( a=0; a< nch; a++)
- {
- if ( mat[a])free_int (mat[a], -1);
- mat[a]=read_matrice ((mat_file)?mat_file[a]:"idmat");
-
- }
- }
-
- l0=(s0)?strlen (s0)/nch:0;
- l1=(s1)?strlen (s1)/nch:0;
-
- if (mchscore==0);
- else if (mchscore==1) score=999999;
- else if (mchscore==2)score=-9999999;
- else
- {
- HERE ("Error: mchscore >2 [FATAL]\n");
- exit (EXIT_FAILURE);
- }
- for ( a=0; a< nch; a++)
- {
- int s;
- c0=(s0)?s0[l0*a+p0]-'A':p0;
- c1=(s1)?s1[l1*a+p1]-'A':p1;
- if ( c0==GAP_CODE)s=(gep[a]!=0)?gep[a]:mat[a][c1][GAP_CODE];
- else if ( c1==GAP_CODE)s=(gep[a]!=0)?gep[a]:mat[a][c0][GAP_CODE];
- else s=mat[a][c0][c1];
-
- if (mchscore==0)score+=s;
- else if (mchscore==1)score=MIN(s, score);
- else if (mchscore==2)score=MAX(s, score);
-
-
- }
-
- if ( use_age>0 && s0 && s1)
- {
-
- int a0, a1;
- int s;
-
- a0=s0[strlen(s0)+1];
- a1=s1[strlen(s1)+1];
-
- a0+=p0;
- a1+=p1;
- s=use_age*FABS((a0-a1))*-1;
-
- if (mchscore==0)score+=s;
- else if (mchscore==1)score=MIN(s, score);
- else if (mchscore==2)score=MAX(s, score);
- }
-
-
- score*=factor;
- return score;
- }
-Sequence * seq2pavie_seq ( Sequence *S, int nch)
- {
- char *buf, *p;
- int a, b;
-
- S->nseq/=nch;
-
- for (b=0; b<S->nseq; b++)
- {
-
- buf=vcalloc ((strlen (S->seq[b])*nch)+10, sizeof (char));
- for ( a=0; a< nch; a++)
- {
- strcat (buf, S->seq[b+(S->nseq)*a]);
- vfree ( S->seq[b+(S->nseq)*a]);
- }
- S->seq[b]=buf;
- /*Code Age on the byte just after the string termination*/
-
- if ((p=strstr (S->seq_comment[b], "FIRSTYEAR")))
- {
- sscanf ( p, "FIRSTYEAR%d", (int*)&(S->seq[strlen(buf)+1]));
- }
-
- }
- return S;
- }
-char **seq2pavie_alp (Sequence *S, int nch)
- {
- int a, n;
- char **alp;
-
- n=S->nseq/nch;
- alp=vcalloc (nch, sizeof (char*));
- for ( a=0; a< nch; a++)
- {
- alp[a]=array2alphabet (S->seq+n*a, n, "-.");
- }
- return alp;
- }
-FILE *output_pavie_aln (Alignment *A, int nch, FILE *fp)
-{
- int a, b, c,d, l, start, end;
- Alignment *B;
- Sequence *S;
- B=declare_aln2(A->nseq*nch+nch, A->len_aln);
-
-
-
- l=A->len_aln/nch;
-
- for ( a=0; a< nch; a++)
- {
- for (b=0; b< A->nseq; b++, B->nseq++)
- {
- sprintf (B->name[B->nseq], "%s.c%d", A->name[b], a);
- start=l*a;end=start+l;
- for (d=0,c=start; c<end; c++, d++)B->seq_al[B->nseq][d]=A->seq_al[b][c];
- B->seq_al[B->nseq][d]='\0';
- }
- if ( a!=nch-1)
- {
- B->name[B->nseq][0]='\0';
- for ( b=0; b<l; b++)B->seq_al[B->nseq][b]='^';
- B->nseq++;
- }
- }
-
- B->len_aln=l;
- fp=output_Alignment_without_header (B,fp);
- S=free_aln (B);
- free_sequence (S, S->nseq);
- return fp;
-
-}
-char **output_pavie_mat_list ( int ***current_mat, double *gep, char **alp, int nch,char *prefix,int cycle, char **mat_name)
-{
- int a;
- char mat_list_name[100];
- FILE *fp;
- char latest[1000];
- char current[1000];
- char command[1000];
-
- sprintf ( mat_list_name, "pavie_matrix%s.cycle_%d.mat_list", prefix, cycle+1);
- fp=vfopen ( mat_list_name, "w");
- fprintf ( stderr, "\n\tOutput Pavie Matrix: %s", mat_list_name);
- for ( a=0; a< nch; a++)
- {
- sprintf ( mat_name[a], "pavie_matrix%s.ch_%d.cy_%d.pavie_mat", prefix,a+1, cycle+1);
- sprintf (latest, "pavie_matrix%s.ch_%d.cy_last.pavie_mat",prefix,a+1);
- sprintf ( current, "matrix.ch%d.pavie_mat", a);
- fprintf ( stderr, "\n\t Channel %d Matrix: %s",a+1, mat_name[a]);
- output_pavie_mat (current_mat[a],mat_name[a],gep[a], alp[a]);
- sprintf ( command, "cp %s %s", mat_name[a], latest);
- system (command);
- sprintf ( command, "cp %s %s", latest, current);
- system (command);
- fprintf ( fp, "%s\n", mat_name[a]);
- }
- vfclose (fp);
- return mat_name;
-}
-
-
-int pavie_pair_wise (Alignment *A,int *ns, int **l_s,Constraint_list *CL )
-{
- double **F; int **T;
- char *x,*y;
- char *ax, *ay;
- int XL, YL, len;
- int i, j;
- double match, gap_inX, gap_inY, MXY=1, GX=2, GY=3;
- int gap=GAP_CODE;
- char *ix, *iy;
- float factor=0.5;
-
-
- /*factor:
- decreases terminal gap penalties with a factor X
- factor=1: terminal gap penalties <=> internal gap penalties
- */
-
-
-
- x=A->seq_al[l_s[0][0]];
- y=A->seq_al[l_s[1][0]];
- XL=strlen (x);
- YL=strlen (y);
-
- ax=vcalloc ( YL+XL+1, sizeof (char));
- ay=vcalloc ( YL+XL+1, sizeof (char));
-
-
- F=declare_double (XL+2, YL+2);
- T=declare_int (XL+2, YL+2);
-
-
- /*Fill stage*/
- F[0][0] = 0;
- for(i = 1; i <=XL; i++)
- {
-
- F[i][0] = F[i-1][0]+(CL->M[x[i-1]-'A'][gap]*factor);
- T[i][0] = GY;
- }
-
- for(j = 1; j <= YL; j++)
- {
- F[0][j] = F[0][j-1]+CL->M[y[j-1]-'A'][gap]*factor;
- T[0][j] = GX;
- }
-
-
- for(i = 1; i <= XL; i++)
- {
- for(j = 1; j <= YL; j++)
- {
-
- match = F[i-1][j-1] + CL->M[x[i-1]-'A'][y[j-1]-'A'];
- gap_inY= F[i-1][j] + (CL->M[x[i-1]-'A'][gap]*(j==YL)?factor:1);
- gap_inX= F[i][j-1] + (CL->M[y[j-1]-'A'][gap]*(i==XL)?factor:1);
-
- if ( match >= gap_inY && match >=gap_inX){F[i][j]=match; T[i][j]=MXY;}
- else if ( gap_inX>=gap_inY){F[i][j]=gap_inX; T[i][j]=GX;}
- else {F[i][j]=gap_inY; T[i][j]=GY;}
- }
- }
- /*Trace back stage*/
- A->score=A->score_aln=F[XL][YL];
-
- i = XL;
- j = YL;
- len=0;
- while(!(i==0 && j==0))
- {
-
- if (T[i][j]==MXY)
- {
- ax[len]=x[--i];
- ay[len]=y[--j];
- }
- else if ( T[i][j]==GY)
- {
- ax[len]=x[--i];
- ay[len]='-';
- }
- else if ( T[i][j]==GX)
- {
- ax[len]='-';
- ay[len]=y[--j];
- }
- len++;
- }
- ax[len]='\0';
- ay[len]='\0';
-
- ix=invert_string (ax); iy=invert_string(ay);
- A=realloc_aln (A,len+1);
-
- sprintf ( A->seq_al[l_s[0][0]], "%s", ix);
- sprintf ( A->seq_al[l_s[1][0]], "%s", iy);
- A->nseq=2;
- A->len_aln=len;
-
- vfree (ax); vfree(ay);vfree(ix); vfree(iy); free_double(F, -1); free_int (T, -1);
- return A->score;
-}
-
-float pavie_aln2id ( Alignment *A, int mode)
-{
- int a, id=0, match=0, l1=0, l2=0, r1, r2, is_res1, is_res2;
-
-
-
- for (a=0; a<A->len_aln; a++)
- {
- r1=A->seq_al[0][a];
- r2=A->seq_al[1][a];
-
-
- is_res1=(!is_gap(r1) && r1!='x' && r1!='X')?1:0;
- is_res2=(!is_gap(r2) && r2!='x' && r2!='X')?1:0;
-
- l1+=is_res1;
- l2+=is_res2;
-
-
- if ( is_res1 && is_res2 )
- {
- match++;
- id+=(r1==r2)?1:0;
- }
- }
-
-
- if ( mode==1)return (match==0)?0:((id*100)/match);
- else if (mode ==2) return (A->len_aln==0)?0:((id*100)/A->len_aln);
- else if (mode ==3) return (MIN(l1,l2)==0)?0:((id*100)/(MIN(l1,l2)));
- else if (mode ==4) return (MAX(l1,l2)==0)?0:((id*100)/(MAX(l1,l2)));
- else if (mode ==5)return (A->score_aln * -1)/*/PAVIE_MAT_FACTOR*/;
- else if (mode ==6)return ((MAX(l1,l2)==0)?0:((A->score_aln)/(MAX(l1,l2))))*-1;
- else
- {
- fprintf ( stderr, "\nUnknown Mode [pavie_aln2id:FATAL:%s]", PROGRAM);
- myexit (EXIT_FAILURE);
- return EXIT_FAILURE;
- }
-
-}
-
-int check_pavie_cl ( char *string)
-{
- if ( !string || string[0]=='\0' ) return 1;
- else if (( string[0]!='_') ||string [strlen (string)-1]!='_')
- {
- fprintf ( stderr, "ERROR: parameters must start and finish with an underscore: _parameters_ [FATAL:%s]\n", PROGRAM);
- myexit (EXIT_FAILURE);
- }
- return 1;
-}
-
-Alignment *pavie_seq2pavie_sort ( Sequence *S, char *mat, char *mode)
-{
- int a, b, c=0, avg_c;
- int **avg;
- double **dm;
- Alignment *A=NULL;
- char **new_order;
-
- if ( vstrstr (mode, "_IDSORT_") || vstrstr (mode, "_MASTERSORT"))
- {
- char *buf;
- buf=vcat ( mode, "_MATDIST_NOPRINT_");
- dm=pavie_seq2pavie_aln (S, mat,buf);
- avg=declare_int (S->nseq, 2);
- if ( vstrstr (mode,"_IDSORT_"))
- {
-
- for ( a=0; a< S->nseq; a++)
- {
- avg[a][0]=a;
- for ( b=0; b<S->nseq; b++)
- if ( b!=a)avg[a][1]+=(int) dm[a][b];
- avg[a][1]/=S->nseq-1;
- }
- sort_int ( avg, 2, 1, 0, S->nseq-1);
-
- c=avg[0][0];
- avg_c=avg[0][1];
- fprintf ( stderr, "\nAVG\t %s\t %s\t %d",S->name[c],"avg", avg_c);
- }
- else if ( vstrstr (mode, "_MASTERSORT"))
- {
- char name[100];
- char *s;
- s=vstrstr(mode, "_MASTERSORT");
- mode=substitute ( mode, "_", " ");
- sscanf (s, " MASTERSORT%s", name);
-
- mode=substitute (mode, " ", "_");
- c=name_is_in_list ( name, S->name, S->nseq, 100);
-
- if ( c==-1)
- {
- fprintf ( stderr, "\nERROR: Sequence %s is not in the dataset [FATAL:%s]", name, PROGRAM);
- myexit (EXIT_FAILURE);
- }
- }
-
- for ( a=0; a<S->nseq; a++)
- {
- avg[a][0]=a;
- if ( a!=c)avg[a][1]=dm[c][a];
- else avg[a][1]=-1;
- }
-
- sort_int ( avg, 2, 1, 0, S->nseq-1);
-
- new_order=declare_char ( S->nseq, 100);
- sprintf (new_order[0], "%s", S->name[c]);
- for ( a=1; a<S->nseq; a++)
- {
- sprintf ( new_order[a], "%s", S->name[avg[a][0]]);
-
- fprintf ( stderr, "\nTOP\t %s\t %s\t %d", S->name[c],new_order[a] , avg[a][1]);
- }
-
- fprintf ( stderr, "\n");
- A=seq2aln (S, NULL,RM_GAP);
- A=reorder_aln (A, new_order, A->nseq);
- vfree ( buf);
- free_double(dm, -1);free_int (avg, -1);free_char (new_order, -1);
- }
- else if ( vstrstr ( mode, "_TREESORT_"))
- {
- A=pavie_seq2pavie_msa (S, mat, mode);
- }
- else
- {
- fprintf ( stderr, "\nERROR: pavie_seq2sort <matrix> <_IDSORT_ | _TREESORT_>");
- }
- return A;
-
-}
-NT_node pavie_seq2pavie_tree (Sequence *S, char *mat, char *mode)
-{
- double **dm;
- char *tree_name,*buf;
-
-
- buf=vcat (mode,"_MATDIST_NOPRINT_");
- dm=pavie_seq2pavie_aln (S, mat,buf);
- dist2nj_tree (dm,S->name, S->nseq,tree_name=vtmpnam (NULL));
-
- free_double(dm, -1);vfree (buf);
-
- return main_read_tree (tree_name);
-}
-
-Alignment* pavie_seq2pavie_msa ( Sequence *S, char *mat_in, char *mode)
-{
- Constraint_list *CL;
- char **alp, *s;
- Alignment *A;
- NT_node **FT, T;
- int a;
- char mat[100];
-
-
- A=seq2aln (S, NULL, RM_GAP);
- CL=declare_constraint_list (S, NULL, NULL, 0, NULL, NULL);
- sprintf ( CL->dp_mode, "myers_miller_pair_wise");
- sprintf ( CL->tree_mode, "nj");
- sprintf ( CL->distance_matrix_mode, "idscore");
- CL=choose_extension_mode ("matrix", CL);
-
- if ( !is_matrix (mat_in))
- {
- FILE *fp;
- fp=vfopen ( mat_in, "r");
- fscanf (fp, "%s", mat);
- vfclose (fp);
- add_warning( stderr, "\nWarning: Multiple Channel Not Supported. Used First Channel Only for MSA [Matrix: %s][WARNING:%s]", mat, PROGRAM);
- }
- else
- {
- sprintf ( mat, "%s", mat_in);
- }
-
- CL->M=read_matrice (mat);
- CL->gop=0;
-
- alp=seq2pavie_alp (S, 1);
- CL->gep=paviemat2gep(CL->M, alp[0]);
-
-
- CL->pw_parameters_set=1;
- CL->local_stderr=stderr;
-
- if ( vstrstr (mode, "_QUICKTREE_"))
- {
- FT=make_tree (A, CL, CL->gop, CL->gep,S, NULL,MAXIMISE);
- T=FT[3][0];
- }
- else if ( (s=vstrstr (mode, "_USETREE")))
- {
- char fname[100];
- mode=substitute ( mode, "_", " ");
- sscanf (s, " USETREE%s", fname);
- mode=substitute (mode, " ", "_");
- T=main_read_tree (fname);
- }
- else
- {
- T=pavie_seq2pavie_tree ( S, mat_in, mode);
- }
-
- for ( a=0; a< A->nseq; a++)ungap (A->seq_al[a]);
-
- tree_aln (T->left,T->right,A,(CL->S)->nseq, CL);
- A=reorder_aln ( A,A->tree_order,A->nseq);
-
- return A;
-}
-/*********************************COPYRIGHT NOTICE**********************************/
-/*© Centro de Regulacio Genomica */
-/*and */
-/*Cedric Notredame */
-/*Tue Oct 27 10:12:26 WEST 2009. */
-/*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