Logo Search packages:      
Sourcecode: t-coffee version File versions  Download package

util_dp_gotoh_nw.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>

#include "io_lib_header.h"
#include "util_lib_header.h"
#include "define_header.h"
#include "dp_lib_header.h"

int idscore_pairseq (char *s1, char *s2, int gop, int gep, int **m, char *comp_mode)
{
  int **I, **D, **M, *P;
  int i, j,l1, l2, l,score, id, igop,match;


  l1=strlen (s1); l2=strlen (s2);
  lower_string (s1); lower_string (s2);
  
  I=declare_int (6,l2+1);D=declare_int (6,l2+1);M=declare_int (6,l2+1);
  
  for (j=0; j<=l2; j++)
    {
      D[0][j]=gep*j;M[0][j]=2*gep*j;D[4][j]=0;
    }
 
  for (i=1; i<=l1; i++)
    {

      I[1][0]=i*gep;
      M[1][0]=2*i*gep;
      
      for (j=1; j<=l2; j++)
      {
        score=m[s1[i-1]-'a'][s2[j-1]-'a'];        
        id=(s1[i-1]==s2[j-1])?1:0;
        
        igop=(i==l1 || j==l2)?0:gop;

        if   ((D[0][j]+gep)>(M[0][j]+igop+gep))   {D[1][j]=D[0][j]+gep;      D[3][j]=D[2][j];        D[5][j]=D[4][j];}
        else                                      {D[1][j]=M[0][j]+igop+gep; D[3][j]=M[2][j];        D[5][j]=M[4][j];}
        
        if ( (I[1][j-1]+gep)>(M[1][j-1]+igop+gep)){I[1][j]=I[1][j-1]+gep;      I[3][j]=I[3][j-1];    I[5][j]=I[5][j-1];}
        else                                      {I[1][j]=M[1][j-1]+igop+gep; I[3][j]=M[3][j-1];    I[5][j]=M[5][j-1];}
        
        match=M[0][j-1]+score;
        if (I[1][j]>match && I[1][j]>D[1][j])     {M[1][j]=I[1][j]           ; M[3][j]=I[3][j];      M[5][j]=I[5][j];}
        else if (D[1][j]>match)                   {M[1][j]=D[1][j]           ; M[3][j]=D[3][j];      M[5][j]=D[5][j];}
        else                                      {M[1][j]=match             ; M[3][j]=M[2][j-1]+id; M[5][j]=M[4][j-1]+1;}
      }
      P=I[0]; I[0]=I[1]; I[1]=P;
      P=I[2]; I[2]=I[3]; I[3]=P;
      P=I[4]; I[4]=I[5]; I[5]=P;
      
      P=D[0]; D[0]=D[1]; D[1]=P;
      P=D[2]; D[2]=D[3]; D[3]=P;
      P=D[4]; D[4]=D[5]; D[5]=P;
      
      P=M[0]; M[0]=M[1]; M[1]=P;
      P=M[2]; M[2]=M[3]; M[3]=P;
      P=M[4]; M[4]=M[5]; M[5]=P;
    }
 



  if ( strstr (comp_mode, "sim2"))
    {
       l=MIN(l1,l2);
       score=(l==0)?0:(M[2][l2]*100)/l;
    }
  else if ( strstr (comp_mode, "sim3"))
    {
       l=MAX(l1,l2);
       score=(l==0)?0:(M[2][l2]*100)/l;
    }
  else if ( strstr (comp_mode, "cov"))
    {
      l=MAX(l1,l2);
      score=(l==0)?0:((M[4][l2]*100)/l);
    }
  else
    {
      //default: simple sim
      l=M[4][l2];
      score=(l==0)?0:(M[2][l2]*100)/l;
    }      
  
  free_int (I, -1);
  free_int (D, -1);
  free_int (M, -1);
  
  return score;
}
            
int test_pair_wise (Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int a,l0, l1, n;
  char buf[VERY_LONG_STRING];
  char *gap, *seq;
  
  l0=strlen (A->seq_al[l_s[0][0]]);
  l1=strlen (A->seq_al[l_s[1][0]]);
  
  n=(l0<5)?l0/2:5;
  A=realloc_aln (A,l0+l1+1);
  
  gap=generate_null(l1-n);
  for (a=0;a<ns[0]; a++)
    {
      seq=A->seq_al[l_s[0][a]];
      sprintf (buf, "%s%s",seq, gap);
      sprintf (seq, "%s", buf);
    }
  vfree (gap);
  gap=generate_null(l0-n);
  
  for (a=0;a<ns[1]; a++)
    {
      seq=A->seq_al[l_s[1][a]];
      sprintf (buf, "%s%s",gap, seq);
      sprintf (seq, "%s", buf);
    }
  vfree(gap);
  

  A->len_aln=strlen (A->seq_al[l_s[0][0]]);
  A->score=A->score_aln=100;
  return 100;
}  

int idscore_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
{
  
  A->score_aln=A->score=idscore_pairseq (A->seq_al[l_s[0][0]], A->seq_al[l_s[1][0]], CL->gop, CL->gep,CL->M, "sim3");
  return A->score_aln;
}

int dp_max (int *trace, int n, ...);
int dp_max (int *trace, int n, ...)
{
  va_list ap;
  int a, v, t, best_v=0;
  
  va_start (ap, n);
  for (a=0; a< n; a++)
    {
      t=va_arg (ap, int);
      v=va_arg (ap, int);

      if (a==0)
      {
        best_v=v;
        trace[0]=t;
      }
      else
      {
        if (v>best_v)
          {
            best_v=v;
            trace[0]=t;
          }
      }
    }
  va_end (ap);
  return best_v;
}

void display_mat (int **M, int l1, int l2, char *title);
void display_mat (int **M, int l1, int l2, char *title)
{
  int a, b;
  
  fprintf ( stdout, "\n\nTitle %s\n", title);
  for ( a=0; a<=l1; a++)
    {
      fprintf ( stdout, "\n");
      for ( b=0; b<=l2; b++)
      fprintf ( stdout, "%3d ", M[a][b]);
    }
}
int glocal_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int ***t, ***m;
  int i,j, l1, l2, n, sub, trace,ntrace, a, b, c, score;
  int gop,rgop,tgop, gep, unmatch;
  int M1, M2, I1, D1, LEN;
  char **al, *char_buf, **aln;
  int **pos0;
  

  l1=strlen (A->seq_al[l_s[0][0]]);
  l2=strlen (A->seq_al[l_s[1][0]]);

  n=1;
  M1=n++;D1=n++;I1=n++;M2=n++;
  t=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
  m=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
  
  
  gop=CL->gop*SCORE_K;
  gep=CL->gep*SCORE_K;
  tgop=gop;
  unmatch=gep;
  
  pos0=aln2pos_simple ( A,-1, ns, l_s);
 
  
  for (j=1; j<=l2; j++)
    {
      m[D1][0][j]=gep*j;

      m[M1][0][j]=2*gep*j;
      m[M2][0][j]=4*gep*j;
    }
  
  
  for (i=1; i<=l1; i++)
    {
      m[I1][i][0]=i*gep;
      m[M2][i][0]=4*i*gep;
      m[M1][i][0]=2*i*gep;
                 
      for ( j=1; j<=l2; j++)
      {
        rgop=(i==l1 || j==1)?0:gop;
        rgop=gop;
        sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);  
        m[M1][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1],D1,m[D1][i-1][j-1],M2,m[M2][i-1][j-1])+sub;
        t[M1][i][j]=trace;
        
        m[D1][i][j]=dp_max (&trace,3, M1,m[M1][i][j-1]+rgop,D1, m[D1][i][j-1]+gep, M2, m[M2][i][j-1]);
        t[D1][i][j]=trace;
        
        m[I1][i][j]=dp_max (&trace,3, M1,m[M1][i-1][j]+rgop, I1, m[I1][i-1][j]+gep, M2, m[M2][i-1][j]);
        t[I1][i][j]=trace;

        m[M2][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1]+tgop,I1, m[I1][i-1][j-1]+tgop,D1,m[D1][i-1][j-1]+tgop,M2,m[M2][i-1][j-1])+unmatch;
        t[M2][i][j]=trace;
        
      }
        
    }
  score=dp_max (&trace,4, M1,m[M1][l1][l2],D1,m[D1][l1][l2],I1, m[I1][l1][l2],M2,m[M2][l1][l2]);
  LEN=0;i=l1;j=l2;
  al=declare_char (2, l1+l2);
  

  trace=t[trace][i][j];
  while (!(i==0 &&j==0))
    {
  
      ntrace=t[trace][i][j];
      if (i==0)
      {
        al[0][LEN]=0;
        al[1][LEN]=1;
        j--;
        LEN++;
      }
      else if ( j==0)
      {
        al[0][LEN]=1;
        al[1][LEN]=0;
        i--;
        LEN++;
      }
      else if ( trace==M1)
      {
        al[0][LEN]=1;
        al[1][LEN]=1;
        i--; j--;
        LEN++;
      }
      else if ( trace==M2)
      {
        al[0][LEN]=1;
        al[1][LEN]=0;
        LEN++;

        al[0][LEN]=0;
        al[1][LEN]=1;
        LEN++;

        i--; j--;
        
      }
      else if ( trace==D1)
      {
        al[0][LEN]=0;
        al[1][LEN]=1;
        j--;
        LEN++;
      }
      else if ( trace == I1)
      {
        al[0][LEN]=1;
        al[1][LEN]=0;
        i--;
        LEN++;
      }
      trace=ntrace;     
      
    }
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);   
  if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);   
  
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));     
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
      {           
        int ch=0;
        for ( b=0; b< LEN; b++)
          {          
            if (al[c][b]==1)
            char_buf[b]=aln[l_s[c][a]][ch++];
            else
            char_buf[b]='-';
          }
        char_buf[b]='\0';
        sprintf (aln[l_s[c][a]],"%s", char_buf);
      }
    }
  
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  free_arrayN((void *)m, 3);
  free_arrayN((void *)t, 3);
  vfree (char_buf);
  free_char (al, -1);
  return score;
}
int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
{
  //adapted from gap_count in MAFFT V 5.5
  int p,s,l, c1, c2;
  int gep,gop;
  int open=3, close=4, gap=5;
  
  gop=CL->gop*SCORE_K;
  gep=CL->gep*SCORE_K;
  
  l=strlen (A->seq_al[ls[0]]);
  
  if (!lg)
    {
      lg=declare_int (6, l);
    }
  
  if ( read_array_size_new (lg[0])<l)
    {
      free_int (lg, -1);
      lg=declare_int (6, l);
    }
  
  for( s=0; s<n; s++ ) 
      {
        c1='x';
        for (p=0; p<l; p++)
          {
            c2=A->seq_al[ls[s]][p];
            
            if (c1!='-' && c2=='-')lg[open][p]++;
            if (c1=='-' && c2!='-')lg[close][p]++;
            if ( c1=='-')lg[gap][p]++;
            c1=c2;
          }
      }

  for (p=0; p<l; p++)
    {
      float go, gc, nn;
      nn=n;
      go=lg[open ][p];
      gc=lg[close][p];
     
      
      lg[GOP][p]=0.5*(1-(go/nn))*gop;
      lg[GCP][p]=0.5*(1-(gc/nn))*gop;
      //Checked locacal gep => gives low quality results
      lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
      
      lg[open][p]=lg[close][p]=lg[gap][p]=0;
    }

  return lg;
}
int gotoh_pair_wise_lgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
{
  int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
  int I, J;
  int M1, I1, D1, LEN;
  char **al, *char_buf, **aln;
  int **pos0;
  
  int gop[2], gcp[2], gep[2];
  static int ***gpl, ***t, ***m;
  static int max_li, max_lj;
  
  //gotoh_pair_wise ( A, ns, l_s,CL);
  //ungap_sub_aln (A, ns[0], l_s[0]);
  //ungap_sub_aln (A, ns[1], l_s[1]);
  
  I=0;J=1;

  
  li=strlen (A->seq_al[l_s[I][0]]);
  lj=strlen (A->seq_al[l_s[J][0]]);
  
  if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
  gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
  gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
  
  
  n=1;
  M1=n++;D1=n++;I1=n++;
  
  if ( li>max_li ||lj>max_lj )
    {
      free_arrayN((void**)t, 3);
      free_arrayN((void**)m, 3);
     
      max_li=li;
      max_lj=lj;
      t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
      m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
    }
  pos0=aln2pos_simple ( A,-1, ns, l_s);
 
  
  for (j=1; j<=lj; j++)
    {
      gep[J]=gpl[J][GEP][j-1];
      m[D1][0][j]=gep[J]*j;
      m[I1][0][j]=m[D1][0][j]-1;
      m[M1][0][j]=m[D1][0][j]-1;
    }
  
  //D1: gap in sequence I
  //I1: gap in sequence J
  
  
  for (i=1; i<=li; i++)
    {
      gep[I]=gpl[I][GEP][i-1];
      gop[I]=gpl[I][GOP][i-1];
      gcp[I]=gpl[I][GCP][i-1];
      
      m[I1][i][0]=i*gep[I];
      m[D1][i][0]= m[I1][i][0]-1;
      m[M1][i][0]= m[I1][i][0]-1;
                 
     
      
      gop[I]=(i==1 || i==li )?0:gop[I];
      gcp[I]=(i==1 || i==li )?0:gcp[I];
      
      
      for ( j=1; j<=lj; j++)
      {
        gep[J]=gpl[J][GEP][j-1];
        gop[J]=gpl[J][GOP][j-1];
        gcp[J]=gpl[J][GCP][j-1];
        
        //gep[J]=gep[I]=(gep[J]+gep[I])/2;
        //gop[J]=gop[I]=(gop[J]+gop[I])/2;
        //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
        

        gop[J]=(j==1 || j==lj )?0:gop[J];
        gcp[J]=(j==1 || j==lj )?0:gcp[J];
        
        
        sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);  
        
        m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
        t[M1][i][j]=trace;
        
        m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
        t[D1][i][j]=trace;
        
        m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
        t[I1][i][j]=trace;
      }
        
    }
  score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
  
  LEN=0;i=li;j=lj;
  al=declare_char (2, li+lj);
  

  trace=t[trace][i][j];
  while (!(i==0 &&j==0))
    {
  
      ntrace=t[trace][i][j];
      if (i==0)
      {
        al[0][LEN]=0;
        al[1][LEN]=1;
        j--;
        LEN++;
      }
      else if ( j==0)
      {
        al[0][LEN]=1;
        al[1][LEN]=0;
        i--;
        LEN++;
      }
      else if ( trace==M1)
      {
        al[0][LEN]=1;
        al[1][LEN]=1;
        i--; j--;
        LEN++;
      }
      else if ( trace==D1)
      {
        al[0][LEN]=0;
        al[1][LEN]=1;
        j--;
        LEN++;
      }
      else if ( trace == I1)
      {
        al[0][LEN]=1;
        al[1][LEN]=0;
        i--;
        LEN++;
      }
      trace=ntrace;     
      
    }
  
  invert_list_char ( al[0], LEN);
  invert_list_char ( al[1], LEN);   
  if ( A->declared_len<=LEN)A=realloc_aln  ( A,2*LEN+1);    
  
  aln=A->seq_al;
  char_buf= vcalloc (LEN+1, sizeof (char));     
  for ( c=0; c< 2; c++)
    {
      for ( a=0; a< ns[c]; a++) 
      {           
        int ch=0;
        for ( b=0; b< LEN; b++)
          {          
            if (al[c][b]==1)
            char_buf[b]=aln[l_s[c][a]][ch++];
            else
            char_buf[b]='-';
          }
        char_buf[b]='\0';
        sprintf (aln[l_s[c][a]],"%s", char_buf);
      }
    }
  
  
  A->len_aln=LEN;
  A->nseq=ns[0]+ns[1];
  vfree (char_buf);
  free_char (al, -1);
  free_int (pos0, -1);
  return score;
}
/*******************************************************************************/
/*                GLOCAL 2                                                     */
/*                                                                             */
/*    makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*    for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
int glocal2_pair_wise (Alignment *IN,int*ns, int **ls,Constraint_list *CL)
{
  int a, b, s=0;
  Alignment *A, *R,*L;
  char *seq, *buf;
  
  buf=vcalloc (1000, sizeof (char));
  seq=vcalloc (1000, sizeof (char));
  
  A=copy_aln (IN,NULL);
  L=copy_aln (IN,NULL);
  R=copy_aln (IN,NULL);
  
  gotoh_pair_wise_sw (A, ns, ls, CL);
  
  HERE ("1");
  for (a=0; a<2; a++)
    {
      for (b=0; b<ns[a]; b++)
      {
        s=ls[a][b];
        sprintf ( seq,"%s", IN->seq_al[s]);
        
        seq[A->order[s][2]]='\0';
        sprintf (L->seq_al[s], "%s", seq);
        sprintf (R->seq_al[s], "%s", seq+A->order[s][3]+1);
      }
    }
  HERE ("2");
  print_sub_aln (A, ns, ls);
  gotoh_pair_wise(L, ns, ls, CL);
  print_sub_aln (L, ns, ls);
  gotoh_pair_wise(R, ns, ls, CL);
  print_sub_aln (R, ns, ls);
  
  IN=realloc_aln (IN, A->len_aln+L->len_aln+R->len_aln+1);
  for (a=0; a<2; a++)
    {
      for (b=0; b<ns[a]; b++)
      {
        s=ls[a][b];
        sprintf (IN->seq_al[s], "%s%s%s",L->seq_al[s], A->seq_al[s], R->seq_al[s]);
      }
    }
  IN->len_aln=strlen (IN->seq_al[s]);
            
  print_sub_aln (IN, ns, ls);
  vfree (seq); vfree (buf);
  free_aln (A); free_aln (L);free_aln (R);
  return IN->score_aln;
}


int gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
      {
/*******************************************************************************/
/*                NEEDLEMAN AND WUNSCH (GOTOH)                                 */
/*                                                                             */
/*    makes DP between the the ns[0] sequences and the ns[1]                 */
/*                                                                             */
/*    for MODE, see the function  get_dp_cost                                */
/*******************************************************************************/
        

/*TREATMENT OF THE TERMINAL GAP PENALTIES*/
/*TG_MODE=0---> gop and gep*/
/*TG_MODE=1---> ---     gep*/
/*TG_MODE=2---> ---     ---*/


          int TG_MODE;
          int l_gop, l_gep;
          int gop, gep;
          int maximise;
/*VARIANLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
      int a, b, i, j;

      int *cc;    
      int *dd,*ddg;
      int   e, eg;

      int lenal[2], len;
      int t, c=0,s, ch;
      int sub;
      int fop;
      int score=0;
        int **pos0;
      static char **al;
      char **aln;
      int ala, alb,LEN;
      char *buffer;
      char *char_buf;
/*trace back variables       */
      FILE       *long_trace=NULL;
      TRACE_TYPE *buf_trace=NULL;
      static TRACE_TYPE **trace;
      TRACE_TYPE k;
      TRACE_TYPE *tr;
      int long_trace_flag=0;
      int dim;
/********Prepare penalties*******/
      gop=CL->gop*SCORE_K;
      gep=CL->gep*SCORE_K;
      TG_MODE=CL->TG_MODE;
      maximise=CL->maximise;
      
      
/********************************/  
/*CLEAN UP AFTER USE*/
      if ( A==NULL)
         {
         free_int (trace,-1);
         trace=NULL;
         free_char (al,-1);
         al=NULL;
         return 0;
         }        

/*DO MEMORY ALLOCATION FOR DP*/

      lenal[0]=strlen (A->seq_al[l_s[0][0]]);
      lenal[1]=strlen (A->seq_al[l_s[1][0]]);
      len= MAX(lenal[0],lenal[1])+1;
      buf_trace=vcalloc ( len, sizeof (TRACE_TYPE));  
      buffer=vcalloc ( 2*len, sizeof (char));   
        al=declare_char (2, 2*len);  
      
      char_buf= vcalloc (2*len, sizeof (char)); 
      

      dd = vcalloc (len, sizeof (int));
      

      cc = vcalloc (len, sizeof (int));
      ddg=vcalloc (len, sizeof (int));
      

      
      if ( len>=MAX_LEN_FOR_DP)
          {
          long_trace_flag=1;
          long_trace=vtmpfile();
          }
      else
          {
         
          dim=(trace==NULL)?0:read_size_int ( trace,sizeof (int*));        
          trace    =realloc_int ( trace,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
          }
      
/*END OF MEMORY ALLOCATION*/
      
      
            /*
            0(s)   +(dd)
              \      |
               \     |
                \    |
                 \   |
                  \  |
                   \ |
                    \|
            -(e)----O
            */ 
                   
      pos0=aln2pos_simple ( A,-1, ns, l_s);


      cc[0]=0;          
      tr=(long_trace_flag)?buf_trace:trace[0];
      tr[0]=(TRACE_TYPE)1;
      for ( j=1; j<=lenal[1]; j++)tr[j]=(TRACE_TYPE)-1;
      if (long_trace_flag)fwrite (buf_trace, sizeof ( TRACE_TYPE),lenal[1]+1, long_trace);
      
      
      t=(TG_MODE==0)?gop:0;
      

      for (cc[0]=0,j=1; j<=lenal[1]; j++)
          {
          
          l_gop=(TG_MODE==0)?gop:0;
          l_gep=(TG_MODE==2)?0:gep;

          cc[j]=t=t+l_gep;
          dd[j]=  t+  gop;
          }

      t=(TG_MODE==0)?gop:0;   
      
      for (i=1; i<=lenal[0];i++)
                  {                 
                  tr=(long_trace_flag)?buf_trace:trace[i];
                  s=cc[0];

                  l_gop=(TG_MODE==0)?gop:0;
                  l_gep=(TG_MODE==2)?0:gep;
                  
                  
                  
                  cc[0]=c=t=t+l_gep;
                  e=t+  gop;
                  tr[0]=(TRACE_TYPE)1;

                  

                  for (eg=0,j=1; j<=lenal[1];j++)
                        {                          
                         
                          sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);  
                              
                        /*get the best Insertion*/
                        l_gop=(i==lenal[0] || i==1 )?((TG_MODE==0)?gop:0):gop;
                        l_gep=(i==lenal[0] || i==1)?((TG_MODE==2)?0:gep):gep;
                  

                        if ( a_better_than_b ( e,c+l_gop, maximise))eg++;
                        else eg=1;  
                        e=best_of_a_b (e, c+l_gop, maximise)+l_gep;
                        
                        /*Get the best deletion*/
                        l_gop=(j==lenal[1] || j==1)?((TG_MODE==0)?gop:0):gop;
                        l_gep=(j==lenal[1] || j==1)?((TG_MODE==2)?0:gep):gep;
                        

                        if ( a_better_than_b ( dd[j], cc[j]+l_gop, maximise))ddg[j]++;
                        else ddg[j]=1;
                        dd[j]=best_of_a_b( dd[j], cc[j]+l_gop,maximise)+l_gep;
                        


                        c=best_int(3,maximise,&fop, e, s+sub,dd[j]);
                        /*Chose Substitution for tie breaking*/
                        if ( fop==0 && (s+sub)==e)fop=1;
                        else if ( fop==2 && (s+sub)==dd[j])fop=1;
                        /*Chose Deletion for tie breaking*/
                        else if ( fop==2 && e==dd[j])fop=1;

                        fop-=1;
                        s=cc[j];
                        cc[j]=c;    

      
                        if ( fop<0)
                              {tr[j]=(TRACE_TYPE)fop*eg;
                              }
                        else if ( fop>0)
                                {tr[j]=(TRACE_TYPE)fop*ddg[j];
                              }
                        else if (fop==0)
                              {tr[j]=(TRACE_TYPE)0;   
                              }                             
                        fop= -2;
                        }
                  if (long_trace_flag)
                      {
                      fwrite ( buf_trace, sizeof (TRACE_TYPE), lenal[1]+1, long_trace);
                      }
                  }
      
      score=c;
      
        i=lenal[0];
      j=lenal[1];
      ala=alb=0;
      

      while (i>=0 && j>=0 && ((i+j)!=0))
                  {
                  if ( i==0)
                        k=-1;
                  else if ( j==0)
                        k=1;
                  else if ( j==0 && i==0)
                        k=1;  
                  else
                          {
                        if (long_trace_flag)
                           {
                           fseek ( long_trace, sizeof (TRACE_TYPE)*((lenal[1]+1)*(i)+j),SEEK_SET);
                           fread ( &k, sizeof (TRACE_TYPE), 1, long_trace);
                           }
                        else
                           {
                           
                           k=trace[i][j];
                           }
                        }
                        
                        
                  if (k==0)
                        {
                        
                        al[0][ala++]=1;
                        al[1][alb++]=1;
                        i--;
                        j--;
                        }           
                  else if (k>0)
                        {
                        
                        for ( a=0; a< k; a++)
                              {
                              al[0][ala++]=1;
                              al[1][alb++]=0;
                              i--;
                              }
                        }
                  else if (k<0)
                        {
                        
                        for ( a=0; a>k; a--)
                              {
                              al[0][ala++]=0;
                              al[1][alb++]=1;
                              j--;
                              }
                        }
                  }
      
      LEN=ala;    
      c=LEN-1;  
      
      

      invert_list_char ( al[0], LEN);
      invert_list_char ( al[1], LEN);     
      if ( A->declared_len<=LEN)A=realloc_aln2  ( A,A->max_n_seq, 2*LEN);     
      aln=A->seq_al;

      for ( c=0; c< 2; c++)
          {
          for ( a=0; a< ns[c]; a++) 
            {           
            ch=0;
            for ( b=0; b< LEN; b++)
                {          
                if (al[c][b]==1)
                  char_buf[b]=aln[l_s[c][a]][ch++];
                else
                  char_buf[b]='-';
               }
            char_buf[b]='\0';
            sprintf (aln[l_s[c][a]],"%s", char_buf);
              }
           }
      
      
      A->len_aln=LEN;
      A->nseq=ns[0]+ns[1];
      

      vfree ( cc);
      vfree (dd);       
      vfree (ddg);
      vfree (buffer);
      vfree (char_buf); 
      vfree (buf_trace);
      free_char ( al, -1);
      free_int (pos0, -1);
      if ( long_trace_flag)fclose (long_trace); 



      return score;
      }
     

/*********************************COPYRIGHT NOTICE**********************************/
/* Centre National de la Recherche Scientifique (CNRS) */
/*and */
/*Cedric Notredame */
/*Fri Oct 26 17:03:04     2007. */
/*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**********************************/

Generated by  Doxygen 1.6.0   Back to index