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

util_dp_drivers.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"


/******************************************************************/
/*                   MAIN DRIVER                                  */
/*                                                                */
/*                                                                */
/******************************************************************/


Constraint_list *profile2list     (Job_TC *job, int nprf)
{
  int *seqlist, *cache1, *cache2;
   
  static CLIST_TYPE *entry;
  Alignment *A1, *A2, *A;
  int a, b, c;
  Constraint_list *SCL, *SCL2;
  static int **L;
  static int max_L_len;
  int  max, n_pairs;
  int **score, n, s1, s2, si, r1, r2;
  char *seqlistb;
  int debug=1;
  int cons, cons_thres, max_n_pairs, tot_n_pairs;
  int *cons_table;
  Constraint_list *CL;
  char *command;
  char *seq;
  char *weight;
  TC_method *M;
  int *iA1, *iA2;
  static char *buf1;

  if ( !buf1) buf1=vcalloc (1000, sizeof (char));
  
  /*initialize the structure*/
  CL=(job->io)->CL;
  M=(job->param)->TCM;
  command=M->executable;
  weight=M->weight;
  seq=(job->param)->seq_c;
  
  
  debug=(getenv("DEBUG_TCOFFEE_profile2list")!=NULL)?1:0;

  if ( debug)print_mem_usage (stderr, "IN");
  seqlistb=vcalloc (10, sizeof (char));
  
  seqlist=string2num_list (seq)+1;

  

  if (!entry)entry=vcalloc (CL->entry_len, sizeof (int));
  entry[SEQ1]=seqlist[1];
  entry[SEQ2]=seqlist[2];
  
  A1=seq2profile(CL->S, seqlist[1]);
  A2=seq2profile(CL->S, seqlist[2]);
    
  SCL=copy_constraint_list ( CL, SOFT_COPY);

  SCL->L=L;SCL->max_L_len=max_L_len;
  SCL->fp=NULL;
  SCL->M=NULL;
  SCL->ne=0;

  SCL2=copy_constraint_list ( CL, SOFT_COPY);
  SCL2->L=L;SCL2->max_L_len=max_L_len;
  SCL2->fp=NULL;
  SCL2->M=NULL;
  SCL2->ne=0;
  
  /*Merge Sequences*/
  
  SCL->S=merge_seq (A1->S, NULL  );
  SCL->S=merge_seq (A2->S, SCL->S);

  /*1: Compare the two profiles and identify the N master pairs*/
  n=(A1->nseq*A2->nseq);
  max=(nprf==0 ||method_uses_structure (command))?n:MIN(nprf,n);
  if ( max<n)
    {
      vsrand(0);
      A=align_two_aln (A1, A2, "blosum62mt",-4,-1, "myers_miller_pair_wise");
      score=declare_int ( A1->nseq*A2->nseq,3);
    
      for (n=0,a=0; a<A1->nseq; a++)
      for ( b=0; b<A2->nseq; b++,n++)
        {
          score[n][0]=a;
          score[n][1]=b;
          score[n][2]=(int)get_seq_fsim ( A->seq_al[a], A->seq_al[b+A1->nseq], "-", NOGROUP, NOMATRIX, AVERAGE_POSITIONS);
        }
      free_aln (A);
      sort_int_inv (score,3,2,0, n-1);
    }
  else
    {
      score=declare_int (n,3);
      for (n=0,a=0; a<A1->nseq; a++)
      for ( b=0; b<A2->nseq; b++,n++)
        {
          score[n][0]=a;
          score[n][1]=b;
        }
    }

  iA1=get_name_index (A1->name,A1->nseq, (SCL->S)->name, (SCL->S)->nseq);
  iA2=get_name_index (A2->name,A2->nseq, (SCL->S)->name, (SCL->S)->nseq);
  
  /*submit the N pairs*/
  for ( n_pairs=0,a=0; a<max; a++)
    {
      s1=score[a][0];s2=score[a][1];      
      sprintf ( buf1, "2 %d %d", iA1[s1], iA2[s2]);
      job=print_lib_job (job, "param->seq_c=%s",buf1);
      
      /*2: Compute the pairewise library*/
      (job->io)->CL=SCL;
      
      SCL=seq2list (job);
      if (debug)fprintf ( stderr, "\n\tProfile aln %s %s %s [(%d,%d):%d %%id]", (SCL->S)->name[iA1[s1]],(SCL->S)->name[iA2[s2]], command, score[a][0], score[a][1],score[a][2]);
      
      /*3: Update the main library with the pairwise library*/
      cache1=seq2inv_pos (A1->seq_al[s1]);
      cache2=seq2inv_pos (A2->seq_al[s2]);

      
      if (debug)fprintf ( stderr, " =>%d pairs", SCL->ne);
      n_pairs+=(SCL->ne>0)?1:0;
      for (c=0; c< SCL->ne; c++)
          {
            si=vread_clist (SCL, c, SEQ1);
            r1=vread_clist(SCL, c, (si==iA1[s1])?R1:R2);
            r2=vread_clist(SCL, c, (si==iA1[s1])?R2:R1);
            
            entry[R1]=cache1[r1];
            entry[R2]=cache2[r2];
            
            entry[WE]=vread_clist(SCL,c,WE);
            entry[CONS]=1;
            add_entry2list(entry, SCL2);
          }

      SCL->ne=0;
      vfree ( cache1);vfree ( cache2);
      compact_list (SCL2, 0, SCL2->ne, "default");
      if (debug)fprintf ( stderr, " =>%d pairs", SCL2->ne);
      
    }
  
  free_sequence (SCL->S,-1);
  vfree (iA1); vfree (iA2);
  
  if (debug)fprintf ( stderr, "\nNPairs=%d", n_pairs);
  
  
  compact_list (SCL2, 0, SCL2->ne, "default"); 
  /*get the concistency distribution*/    
  cons_table=vcalloc ( 101, sizeof (int));
  for (c=0; c< SCL2->ne; c++)
    {
      entry[R1]=vread_clist(SCL2,c,R1);
      entry[R2]=vread_clist(SCL2,c,R2);
      entry[WE]=vread_clist(SCL2,c,WE);
      entry[CONS]=vread_clist(SCL2,c,CONS);
      cons=(entry[CONS]*100)/n_pairs;
      cons_table[cons]++;
    }
  
  /*Identify the threshold*/
  max_n_pairs=(int)((float)(MIN(A1->len_aln, A2->len_aln)*2));
  
  for (cons_thres=0,tot_n_pairs=0,c=100; c>=0; c--)
    {

      tot_n_pairs+=cons_table[c];
      if ( tot_n_pairs>=max_n_pairs){cons_thres=c;c=-1;}
    }
  vfree (cons_table);

  /*Produce the library*/
  for (c=0; c< SCL2->ne; c++)
    {
      entry[R1]=vread_clist(SCL2,c,R1);
      entry[R2]=vread_clist(SCL2,c,R2);
      entry[WE]=vread_clist(SCL2,c,WE);
      entry[CONS]=vread_clist(SCL2,c,CONS);
      entry[WE]/=entry[CONS];
      cons=(entry[CONS]*100)/n_pairs;
      
      
      if (cons>=cons_thres)add_entry2list(entry, CL);
    }

  if ( !seq2R_template_profile (CL->S,seqlist[1]))free_aln (A1);
  if ( !seq2R_template_profile (CL->S,seqlist[2]))free_aln (A2);
  
  vfree(seqlistb);
  vfree(seqlist-1)  ;
  free_int (SCL->L, -1);free_constraint_list ( SCL);
  
  if (debug)print_mem_usage (stderr, "BEF FREE SCL2");
  free_int (SCL2->L, -1);free_constraint_list ( SCL2);
  
  if (debug)
    {
      
      fprintf ( stderr, "\nCL->ne=%d", CL->ne);
      print_mem_usage (stderr, "OUT");
    }
 
  (job->io)->CL=CL;
  return CL;
}

int method_uses_structure( char *method)
{
  if ( strm ( method, "sap_pair"))return 1;
  else if (  strm ( method, "fugue_pair"))return 1;
  return 0;
}
Constraint_list *seq2list     ( Job_TC *job)
    {
      char *mode;
      Alignment *A=NULL;
      Constraint_list *PW_CL;
      Constraint_list *RCL;
      Constraint_list *fugue_CL;
      int full_prf, nprf;
      int *seqlist;
      int for_sap, for_fugue;
      
      
      int ktup;
      static char *s1, *s2;
      Sequence *S, *STL;
      Constraint_list *CL;
      char *command;
      
      char *seq;
      char *weight;
      TC_method *M;
            
      M=(job->param)->TCM;
      mode=M->executable;
      weight=M->weight;
      
      PW_CL=M->PW_CL;
      
      CL=(job->io)->CL;
      seq=(job->param)->seq_c;
      
      S=(CL)?CL->S:NULL;
      STL=(CL)?CL->STRUC_LIST:NULL;
      
      
            
      seqlist=string2num_list (seq)+1;
      
      if (!s1)s1=vcalloc ( 100, sizeof (char));
      if (!s2)s2=vcalloc ( 100, sizeof (char));
      

      sprintf (s1, "%s", (CL->S)->name[seqlist[1]]);
      sprintf (s2, "%s", (CL->S)->name[seqlist[2]]);

/*Proteins*/      
      
      
      if ( strncmp (CL->profile_comparison, "full", 4)==0)
           {
             full_prf=1;
             if ( CL->profile_comparison[4])nprf=atoi ( CL->profile_comparison+4);
             else
             nprf=0;
           }
      else
           {
             full_prf=0;
           }

      if ((method_uses_structure (mode)||full_prf) && (seq2R_template_profile (CL->S,seqlist[1]) || seq2R_template_profile (CL->S,seqlist[2]) ))
        {
          RCL=profile2list     (job, nprf);
        }
      else if (    strm ( mode,"fast_pair")         || strm (mode, "ifast_pair")\
               || strm ( mode, "diag_fast_pair")|| strm (mode, "idiag_fast_pair")\
               || strm ( mode, "blast_pair")    || strm (mode, "lalign_blast_pair") \
               || strm ( mode, "viterbi_pair")  || strm (mode, "slow_pair") \
               || strm ( mode, "islow_pair")    || strm (mode , "lalign_len_pair")\
               || strm ( mode, "lalign_id_pair")|| strm (mode, "prrp_aln") \
               || strm ( mode, "clustalw_aln")  || strm (mode, "cdna_fast_pair") \
               || strncmp (mode,"cdna_fast_pair",14)==0           \
               )
      {
        
        A=fast_pair (job);
        RCL=aln2constraint_list ((A->A)?A->A:A, CL,weight);
      }
      else if ( strm ( mode, "exon_pair"))
      {
        int a;
        char weight2[1000];
        
        A=fast_pair (job);
        sprintf ( weight2, "%s_subset_objOBJ-",weight);
        RCL=aln2constraint_list (A, CL,weight2);
      }
      else if ( strm ( mode, "exon_pair"))
      {
        int a;
        char xmode[100], xalp[100], ualp[100], lalp[100],mode2[100];
        static char pxalp[100];
        static Sequence *CS;

        if ( strm (mode, "exon_pair"))sprintf (mode2, "pair_obj");
        else sprintf ( mode2, "%s", mode);
        
        sprintf ( xalp, "%s",strchr(mode2, '_')+1); 
        sprintf ( ualp, "%s",strchr(mode2, '_')+1); 
        sprintf ( lalp, "%s",strchr(mode2, '_')+1); 
        
        sprintf ( xmode, "subset_%s%s-", upper_string(ualp), lower_string (lalp));
        
        if (!CS || !strm (xalp,pxalp))
          {
            free_sequence (CS, -1);
            CS=duplicate_sequence (CL->S);
            CS=keep_residues_in_seq (CS,strchr(xmode, '_'), 'X');
            fprintf (stderr, "\n\tALP:%s\n",xmode);
          }
        sprintf ( pxalp, "%s", xalp);
        PW_CL->S=CS;
        A=fast_pair (job);
        RCL=aln2constraint_list (A, CL,xmode);
      }
/*STRUCTURAL METHODS*/
      else if ( strm (mode, "sap_pair"))
        {
          RCL=sap_pair(seq, CL);
        }
      else if ( strm (mode, "thread_pair"))
        {
          RCL=thread_pair (seq, CL);
        }
      else if ( strm (mode, "fugue_pair"))
        {
          RCL=fugue_pair (seq, CL);
        }
      else if ( strm (mode, "lsqman_pair"))
        {
          RCL=lsqman_pair(seq, CL);
        }
      else if ( strm ( mode, "align_pdb_pair"))
          {
            RCL=align_pdb_pair ( seq,"gotoh_pair_wise", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
          }
      else if ( strm ( mode, "lalign_pdb_pair"))
          {
            RCL=align_pdb_pair ( seq,"sim_pair_wise_lalign", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
          }
      else if ( strm ( mode, "align_pdb_pair_2"))
          {
            RCL=align_pdb_pair_2 ( seq, CL);
          }
      else
          {
            fprintf ( CL->local_stderr, "\nERROR: THE FUNCTION %s DOES NOT EXIST [FATAL:%s]\n", mode, PROGRAM);crash("");
          }
      RCL=(RCL==NULL)?CL:RCL;
      
      vfree ( seqlist-1);
      free_aln (A);
      return RCL; 
    }

Constraint_list *method2pw_cl (TC_method *M, Constraint_list *CL)
    {
      char *mode;
      Constraint_list *PW_CL=NULL;
      Sequence *S;
      char mat[100];
      char group_mat[100];
      char dp_mode[100];
      int ktup;
      
      mode=M->executable;
      PW_CL=copy_constraint_list ( CL, SOFT_COPY);
      PW_CL->pw_parameters_set=1;
      
      S=(PW_CL)?PW_CL->S:NULL;
      
      /*DNA or Protein*/
      if (  strm ( (PW_CL->S)->type, "PROTEIN"))
         {
          
           sprintf ( mat, "%s", PW_CL->method_matrix);
           sprintf (group_mat, "vasiliky");
           ktup=2;
         }
      else if (  strm ( (PW_CL->S)->type, "DNA"))
         {
            sprintf(group_mat, "idmat"); 
            sprintf(mat, "idmat");
            ktup=5;
         }  
      
      if ( strm2 ( mode,"fast_pair", "ifast_pair"))
          {
            PW_CL->maximise=1;
              PW_CL->TG_MODE=1;
              PW_CL->use_fragments=0;
              if ( !PW_CL->use_fragments)PW_CL->diagonal_threshold=0;
              else PW_CL->diagonal_threshold=6;
            
              sprintf (PW_CL->dp_mode, "fasta_pair_wise");
              PW_CL->ktup=ktup;
            sprintf (PW_CL->matrix_for_aa_group, group_mat);

              if ( strm ( mode, "fast_pair"))
                {
              PW_CL->L=NULL;
                  PW_CL->M=read_matrice (mat);
                  PW_CL->get_dp_cost=slow_get_dp_cost;
                  PW_CL->evaluate_residue_pair=evaluate_matrix_score;

                  PW_CL->extend_jit=0;
                  PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
                  PW_CL->gep=-1;
                }
          }
      else if ( strm2 ( mode,"diag_fast_pair","idiag_fast_pair"))
          {
            PW_CL->L=NULL;
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            PW_CL->S=CL->S;
            
            PW_CL->use_fragments=1;
            PW_CL->diagonal_threshold=3;
            
            sprintf (PW_CL->dp_mode, "fasta_pair_wise");
            PW_CL->ktup=1;
            sprintf (PW_CL->matrix_for_aa_group, group_mat);
            PW_CL->M=read_matrice (mat);
            PW_CL->get_dp_cost=slow_get_dp_cost;
            PW_CL->evaluate_residue_pair=evaluate_matrix_score;
            
            PW_CL->extend_jit=0;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->gep=-1;             
          }
      else if ( strm ( mode,"blast_pair"))
          {
            PW_CL->L=NULL;
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            
            PW_CL->use_fragments=0;
            
            PW_CL->M=read_matrice (mat);
            PW_CL->pair_wise=gotoh_pair_wise;
            PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->gep=-2;
            PW_CL->ktup=ktup;
            sprintf (PW_CL->matrix_for_aa_group, group_mat);
            PW_CL->extend_jit=0;
          }
      else if ( strm ( mode,"lalign_blast_pair"))
          {
            PW_CL->L=NULL;
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            
            PW_CL->use_fragments=0;
            
            PW_CL->M=read_matrice (mat);
            PW_CL->pair_wise=sim_pair_wise_lalign;
            PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->lalign_n_top=10;
            PW_CL->gep=-2;
            PW_CL->ktup=ktup;
            sprintf (PW_CL->matrix_for_aa_group, group_mat);
            PW_CL->extend_jit=0;
          }
      else if ( strm ( mode,"viterbi_pair"))
          {
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            
            PW_CL->use_fragments=0;
            sprintf (PW_CL->dp_mode, "viterbi_pair_wise");
            
            PW_CL->L=NULL;
            PW_CL->M=read_matrice (mat);
            PW_CL->get_dp_cost=slow_get_dp_cost;
            PW_CL->evaluate_residue_pair=evaluate_matrix_score;
            
            PW_CL->extend_jit=0;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->gep=-1;
          }

      else if ( strm2 ( mode,"slow_pair","islow_pair" ))
          {
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            
            PW_CL->use_fragments=0;
            sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
            PW_CL->ktup=ktup;
            sprintf (PW_CL->matrix_for_aa_group, group_mat);
            
            if ( strm ("slow_pair", mode))
            {
              PW_CL->L=NULL;
              PW_CL->M=read_matrice (mat);
              PW_CL->get_dp_cost=slow_get_dp_cost;
              PW_CL->evaluate_residue_pair=evaluate_matrix_score;
              
              PW_CL->extend_jit=0;
              PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
              PW_CL->gep=-1;
                          
            }
          } 
      
      else if (strm ( mode, "exon_pair"))
        {
          int a;
          char xmode[100];
          char xalp[100];
          char ualp[100];
          char lalp[100];
          char mode2[100];
          
          PW_CL->maximise=1;
          PW_CL->TG_MODE=1;
          
          PW_CL->use_fragments=0;
          sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
          PW_CL->ktup=ktup;
          sprintf (PW_CL->matrix_for_aa_group, group_mat);
          PW_CL->L=NULL;

          PW_CL->M=read_matrice ("exon_mt");
                           
          for ( a=0; a<60; a++)
            {
            PW_CL->M['x'-'A'][a]=0;
            PW_CL->M[a]['x'-'A']=0;
            PW_CL->M['X'-'A'][a]=0;
            PW_CL->M[a]['X'-'A']=0;
            }
          PW_CL->evaluate_residue_pair=evaluate_matrix_score;
          PW_CL->extend_jit=0;
          PW_CL->gop=get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
          PW_CL->gep=-1;
          PW_CL->gop=0;PW_CL->gep=-1;
        }
      else if ( strnm ( mode,"xpair",5) )
        {
          int a;
          char xmode[100];
          char xalp[100];
          char ualp[100];
          char lalp[100];
          char mode2[100];
          
          PW_CL->maximise=1;
          PW_CL->TG_MODE=1;
          
          PW_CL->use_fragments=0;
          sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
          PW_CL->ktup=ktup;
          sprintf (PW_CL->matrix_for_aa_group, group_mat);
          PW_CL->L=NULL;

          PW_CL->M=read_matrice (strm (mode, "exon_pair")?"exon_mt":mat);
                           
          for ( a=0; a<60; a++)
            {
            PW_CL->M['x'-'A'][a]=0;
            PW_CL->M[a]['x'-'A']=0;
            PW_CL->M['X'-'A'][a]=0;
            PW_CL->M[a]['X'-'A']=0;
            }
          PW_CL->evaluate_residue_pair=evaluate_matrix_score;
          PW_CL->extend_jit=0;
          PW_CL->gop=get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
          PW_CL->gep=-1;
          PW_CL->gop=0;PW_CL->gep=-1;
        }
      else if ( strm ( mode , "lalign_len_pair"))
        {   
          PW_CL->L=NULL;
          PW_CL->maximise=1;
          PW_CL->TG_MODE=1;
          PW_CL->use_fragments=0;
          
          PW_CL->M=read_matrice ("blosum50mt");
          PW_CL->pair_wise=sim_pair_wise_lalign;
          PW_CL->evaluate_residue_pair=evaluate_matrix_score;
          PW_CL->get_dp_cost=slow_get_dp_cost;
          PW_CL->lalign_n_top=CL->lalign_n_top;
          PW_CL->gop= -10;
          PW_CL->gep= -4;
          PW_CL->ktup=ktup;
          sprintf (PW_CL->matrix_for_aa_group, group_mat);
          PW_CL->extend_jit=0;
        }
      else if ( strm ( mode , "lalign_id_pair"))
        {   
          PW_CL->L=NULL;
          PW_CL->maximise=1;
          PW_CL->TG_MODE=1;
          PW_CL->use_fragments=0;
          
          PW_CL->M=read_matrice ("blosum50mt");
          PW_CL->pair_wise=sim_pair_wise_lalign;
          
          /*
          PW_CL->get_dp_cost=get_dp_cost_pw_matrix;
          PW_CL->evaluate_residue_pair=NULL;
          */

          PW_CL->evaluate_residue_pair=evaluate_matrix_score;
          PW_CL->get_dp_cost=slow_get_dp_cost;

          
          PW_CL->lalign_n_top=CL->lalign_n_top;
          PW_CL->gop= -10;
          PW_CL->gep= -4;
          PW_CL->ktup=ktup;
          sprintf (PW_CL->matrix_for_aa_group, group_mat);
          PW_CL->extend_jit=0;
        }
/*CDNA*/
      else if ( strm ( mode, "cdna_cfast_pair"))
          { 
            PW_CL->L=NULL;
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            PW_CL->S=CL->S;
            PW_CL->use_fragments=0;
            sprintf (PW_CL->dp_mode, "cfasta_cdna_pair_wise");

            PW_CL->M=read_matrice (strcpy ( mat, "blosum62mt"));
            PW_CL->extend_jit=0;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->gep=-1;
            PW_CL->f_gop=CL->f_gop;
            PW_CL->f_gep=CL->f_gep;
            PW_CL->get_dp_cost=get_dp_cost;
            PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
            PW_CL->ktup=1;
          }
      else if ( strm ( mode, "cdna_fast_pair") ||  strncmp (mode,"cdna_fast_pair",14)==0)
          { 
            PW_CL->L=NULL;
            PW_CL->maximise=1;
            PW_CL->TG_MODE=1;
            PW_CL->use_fragments=0;
            sprintf (PW_CL->dp_mode, "fasta_cdna_pair_wise");
            
            PW_CL->M=read_matrice (strcpy ( mat, "blosum62mt"));
            PW_CL->extend_jit=0;
            PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
            PW_CL->gep=-1;
            
            if ( !strm ( mode, "cdna_fast_pair"))
            {
              sscanf (mode+14, "%d_%d", &PW_CL->f_gop, &PW_CL->f_gep);
              PW_CL->f_gop*=(-1);
              PW_CL->f_gep*=(-1);
            }
            else 
            {
              PW_CL->f_gop=2*PW_CL->gop;
              PW_CL->f_gep=2*PW_CL->gep;
            }
            
            PW_CL->get_dp_cost=get_dp_cost;
            PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
            PW_CL->ktup=1;
          }
      else
        {
          free_constraint_list (PW_CL);
          PW_CL=NULL;
        }
      return PW_CL;
    }
/******************************************************************/
/*                   MULTIPLE ALIGNMENTS                          */
/*                                                                */
/*                                                                */
/******************************************************************/
Alignment * compute_prrp_aln (Alignment *A, Constraint_list *CL)
   {
       char *tmpseq=NULL;
       char *tmpaln=NULL;
       char command[10000];
       Sequence *S;

       tmpseq=vtmpnam(NULL);
       tmpaln=vtmpnam(NULL);
       

       A=seq2aln ( CL->S, A, 1);
       output_gotoh_seq (tmpseq, A);
       sprintf ( command, "prrp -E/dev/null -o%s -F9 %s >/dev/null", tmpaln, tmpseq);
       my_system (command);
       if (!check_file_exists(tmpaln)){return NULL;}
       S=get_fasta_sequence (tmpaln, NULL);

       S->contains_gap=0;
       A=seq2aln(S, A,0);
       free_sequence (S, S->nseq);
       
       return A;
   }
Alignment * aln2clustalw_aln (Alignment *B, Constraint_list *CL)
{
     char *tmpseq=NULL;
     char *tmpaln=NULL;
     char command[10000];
     Alignment *A=NULL;
     FILE *fp;
     int a;

     

     A=copy_aln (B,A); 
     tmpseq=vtmpnam(NULL);
     tmpaln=vtmpnam(NULL);
          
     fp=vfopen (  tmpseq, "w");
     for ( a=0; a< A->nseq; a++)
       {
       ungap ( A->seq_al[a]);
       fprintf ( fp, ">%s\n%s\n", A->name[a], A->seq_al[a]);
       }
     vfclose ( fp);
     
     sprintf ( command, "clustalw %sinfile=%s %soutorder=input %soutfile=%s %s", CWF,tmpseq, CWF, tmpaln, CWF, TO_NULL_DEVICE);
    
     my_system (command);
     if (!check_file_exists(tmpaln))return NULL;
     
    
     
     A->nseq=0;
     A=main_read_aln(tmpaln,A);
     
     for ( a=0; a< B->nseq; a++)
       sprintf (B->seq_al[a], "%s", A->seq_al[a]);
     B->len_aln=A->len_aln;
     
     
     free_aln (A);
     
     vremove( tmpseq);
     vremove (tmpaln);
     return B;
   }  
Alignment * compute_tcoffee_aln_quick (Alignment *A, Constraint_list *CL)
   {
       char *tmpseq=NULL;
       char *tmpaln=NULL;
       char command[10000];
     

       tmpseq=vtmpnam(NULL);
       tmpaln=vtmpnam(NULL);

       if ( CL)A=seq2aln ( CL->S, A, 1);
       output_fasta_seq (tmpseq, A);
       
       sprintf ( command, "t_coffee  -seq %s -very_fast -outfile %s -quiet ",tmpseq,tmpaln);
       
       my_system (command);
       if (!check_file_exists(tmpaln))return NULL;
       A->nseq=0;
       A=main_read_aln(tmpaln,A);
       
       vremove( tmpseq);
       vremove (tmpaln);
       return A;
   }  

Alignment * compute_clustalw_aln (Alignment *A, Constraint_list *CL)
   {
       char *tmpseq=NULL;
       char *tmpaln=NULL;
       char command[10000];
     

       tmpseq=vtmpnam(NULL);
       tmpaln=vtmpnam(NULL);

       A=seq2aln ( CL->S, A, 1);
       output_fasta_seq (tmpseq, A);
       
       sprintf ( command, "clustalw %sinfile=%s %soutfile=%s %s",CWF, tmpseq,CWF, tmpaln,TO_NULL_DEVICE);
       
       my_system (command);
       if (!check_file_exists(tmpaln))return NULL;
       A->nseq=0;
       A=main_read_aln(tmpaln,A);
       
       vremove( tmpseq);
       vremove (tmpaln);
       return A;
   }  

Alignment * realign_block ( Alignment *A, int col1, int col2, char *pg)
{
  /*Uses pg: (pg -infile=<infile> -outfile=<outfile> to realign the block [col1 col2[
    Only guaranteed if pg can handle empty sequences
    set pg to NULL to use the default program
  */
  
 
  Alignment *L, *M, *R;
  char *seq_name;
  char *aln_name;
  char command[1000], script[1000];
 
  

  seq_name=vtmpnam(NULL);
  aln_name=vtmpnam (NULL);

  L=copy_aln (A, NULL);
  M=copy_aln (A, NULL);
  R=copy_aln (A, NULL);

  L=extract_aln ( L, 0, col1);
  M=extract_aln ( M, col1, col2);
  R=extract_aln ( R, col2, A->len_aln);
  output_fasta_seq (seq_name, M);

  sprintf ( script, "%s", (pg==NULL)?"t_coffee":pg);
  
  sprintf ( command, "%s -infile=%s -outfile=%s %s", script,seq_name, aln_name, TO_NULL_DEVICE);
  my_system ( command);
  free_aln (M);
  M=main_read_aln (aln_name, NULL);
 
  
  M=reorder_aln (M, L->name,L->nseq);
  L=aln_cat (L, M);
  L=aln_cat (L, R);
  A=copy_aln (L, A);
  free_aln (L);free_aln (M); free_aln (R);
  
  return A;
}
  
  
  


/******************************************************************/
/*                  DNA                                           */
/*                                                                */
/*                                                                */
/******************************************************************/


/******************************************************************/
/*                   STRUCTURES                                   */
/*                                                                */
/*                                                                */
/******************************************************************/


      


Constraint_list * align_pdb_pair_2 (char *seq, Constraint_list *CL)
        {
          char *tmp_name=NULL;
          int s1, s2;
          
          
          static char *command;
          static char *program;
          
          tmp_name=vtmpnam ( NULL);
          
          if ( !program)program=vcalloc ( LONG_STRING, sizeof (char));
          if ( !command)command=vcalloc ( LONG_STRING, sizeof (char));
          
#ifndef     ALIGN_PDB_4_TCOFFEE
          if ( getenv ( "ALIGN_PDB_4_TCOFFEE")==NULL)crash ("ALIGN_PDB_4_TCOFFEE IS NOT DEFINED");
          else  sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
#else
          if ( getenv ( "ALIGN_4_TCOFFEE")==NULL)sprintf (program, "%s", ALIGN_PDB_4_TCOFFEE);
          else sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
#endif

          atoi(strtok (seq,SEPARATORS));
          s1=atoi(strtok (NULL,SEPARATORS));
          s2=atoi(strtok (NULL,SEPARATORS));

          sprintf ( command , "%s -in P%s P%s -gapopen=-40 -max_delta=2.5 -gapext=0 -scale=0 -hasch_mode=hasch_ca_trace_bubble -maximum_distance=10 -output pdb_constraint_list -outfile stdout> %s%s",program, (CL->S)->file[s1], (CL->S)->file[s2], get_cache_dir(),tmp_name);        
          
          my_system  ( command);
          CL=read_constraint_list_file(CL, tmp_name);
         

          vremove ( tmp_name);


          return CL;
      }

Constraint_list *align_pdb_pair   (char *seq_in, char *dp_mode,char *evaluate_mode, char *file, Constraint_list *CL, Job_TC *job)
 {
          int s1, s2;
          char seq[1000];
          char name1[1000];
          char name2[1000];
          

          Constraint_list *PWCL;
          Alignment *F;
          
          sprintf ( seq, "%s",seq_in); 
          atoi(strtok (seq,SEPARATORS));
          s1=atoi(strtok (NULL,SEPARATORS));
          s2=atoi(strtok (NULL,SEPARATORS));
          


          sprintf (name1, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s1], (CL->S)->name[s2], dp_mode);
          sprintf (name2, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s2], (CL->S)->name[s1], dp_mode);
          
          
          if ( check_file_exists (name1) &&   is_lib(name1))CL=read_constraint_list_file(CL,name1);
          else if ( check_file_exists (name2) &&   is_lib(name2))CL=read_constraint_list_file(CL,name2);
          else
            {
            PWCL=set_constraint_list4align_pdb ( CL,s1,dp_mode, evaluate_mode, NULL);
            PWCL=set_constraint_list4align_pdb ( CL,s2,dp_mode, evaluate_mode, NULL);         
            ((job->param)->TCM)->PW_CL=PWCL;
            F=fast_pair (job);
            output_constraints (name1, "100", F);
            CL=aln2constraint_list (F, CL, "100");
            free_aln (F);
            }
          return CL;
      } 

Constraint_list *align_pdb_pair_old   (char *seq, Constraint_list *CL)
 {
            char *tmp_pdb1, *tmp_pdb2;
          int s1, s2;
          char *tmp_name1, *tmp_name2, *tmp_name;
          char command[STRING];
          Alignment *F;
          
          atoi(strtok (seq,SEPARATORS));
          s1=atoi(strtok (NULL,SEPARATORS));
          s2=atoi(strtok (NULL,SEPARATORS));


          tmp_name1=vcalloc (100, sizeof (char));         
          sprintf ( tmp_name1, "%s_%s.align_pdb_results", (CL->S)->name[s1], (CL->S)->name[s2]);
          tmp_name2=vcalloc (100, sizeof (char));
          sprintf ( tmp_name2, "%s_%s.align_pdb_results", (CL->S)->name[s2], (CL->S)->name[s1]);
          if (is_aln (tmp_name1)){tmp_name=tmp_name1;}
          else if ( is_aln (tmp_name2))tmp_name=tmp_name2;
          else tmp_name=tmp_name1;
          
          

      
      
          tmp_pdb1=normalize_pdb_file(seq2P_template_file(CL->S,s1),(CL->S)->seq[s1], vtmpnam (NULL));
          tmp_pdb2=normalize_pdb_file(seq2P_template_file(CL->S,s2),(CL->S)->seq[s2], vtmpnam (NULL));
      
            
          sprintf ( command , "align_pdb -in P%s P%s -gapopen=-40 -max_delta=2.5 -gapext=0 -scale=0 -hasch_mode=hasch_ca_trace_bubble -pdb_name %s %s -output clustalw_aln -outfile stdout > %s%s",tmp_pdb1,tmp_pdb2,(CL->S)->name[s1], (CL->S)->name[s2],get_cache_dir(),tmp_name);          

          my_system  ( command);
          
          vremove (tmp_pdb1);
          vremove (tmp_pdb2);
 
          F=main_read_aln(tmp_name, NULL);
          F=fix_aln_seq (F, CL->S);
          CL=aln2constraint_list (F, CL, "100");
          free_aln (F);
          
          return CL;
      } 

Constraint_list * fugue_pair ( char *in_seq, Constraint_list *CL)
        {

        char seq[1000];
        int s1, s2;
        
        sprintf ( seq, "%s", in_seq);
            
        check_program_is_installed ( "fugue_client", FUGUE_4_TCOFFEE, "FUGUE_4_TCOFFEE","not public", 1); 
        atoi(strtok (seq,SEPARATORS));
        s1=atoi(strtok (NULL,SEPARATORS));
        s2=atoi(strtok (NULL,SEPARATORS));

        CL=fugue_pair2(s1, s2, CL);
        CL=fugue_pair2(s2, s1, CL);

        return CL;
      }

        

       
Constraint_list* fugue_pair2 ( int s1, int s2, Constraint_list *CL)
  {
    char *fugue_result;
    char *seq_file;
    int is_pdb=1;
    char command[LONG_STRING];
    Alignment *F=NULL;
    Sequence *STL;
    FILE *fp;
    char *pdb_id;
    

    STL=(CL)?CL->STRUC_LIST:NULL;

    if ( !(CL->S) || !((CL->S)->T[s1]) || !((CL->S)->T[s1])->P || seq2T_value ((CL->S), s1, "public_pdb", "_P_"))return CL;
    else pdb_id=seq2T_value(CL->S, s1, "pdb_id", "_P_");
    
    fugue_result=vcalloc (100, sizeof (char));
    sprintf (fugue_result, "%s_%s.fugue_results",pdb_id, (CL->S)->name[s2]);
    

    if ( check_file_exists (fugue_result) &&  (F=main_read_aln(fugue_result, NULL))!=NULL)
      {
      if ( getenv4debug ("DEBUG_FUGUE")){fprintf ( stderr, "\n[DEBUG_FUGUE:fugue_pair2]Fugue_Alignment in Cache: %s", fugue_result);}     
      free_aln (F);
      }
    else
      {
      seq_file=vtmpnam (NULL);
      fp=vfopen (seq_file, "w");
      fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s2],(CL->S)->seq[s2] );
      vfclose (fp);
      sprintf ( command, "%s -pdb %s -pep %s > %s%s", FUGUE_4_TCOFFEE,pdb_id, seq_file,get_cache_dir(), fugue_result);
      
      if ( getenv4debug ("DEBUG_FUGUE")){fprintf ( stderr, "\n[DEBUG_FUGUE:fugue_pair2]command %s",command);}     

      my_system ( command);
      }
    
    F=main_read_aln (fugue_result, NULL);

    if ( !F) 
      {
      fprintf ( stderr, "\n\tFugue failed:\n\t%s\n", command);
      }
    else
      {
    
      sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
      sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
      F=fix_aln_seq (F, CL->S);
      CL=aln2constraint_list (F, CL, "sim");
      }
       
    
    free_aln (F);
    vfree (fugue_result);
    return CL;
  }
Constraint_list * lsqman_pair ( char *in_seq, Constraint_list *CL)
        {
        FILE *fp;
        static CLIST_TYPE *entry;
        char command[STRING];
        char seq[1000];
        int s1, s2;
        char *seq_file, *lsqman_result, *tmp_name1;
        Alignment *F=NULL;
        int n_failure=0;
       
        sprintf ( seq, "%s", in_seq);
            

        if ( !entry)entry=vcalloc ( LIST_N_FIELDS, sizeof ( CLIST_TYPE ));
          
        atoi(strtok (seq,SEPARATORS));
        s1=atoi(strtok (NULL,SEPARATORS));
        s2=atoi(strtok (NULL,SEPARATORS));

             
       tmp_name1=vcalloc (100, sizeof (char));      
       sprintf ( tmp_name1, "%s_%s.lsqman_aln", (CL->S)->name[s1], (CL->S)->name[s2]);
       if ( check_file_exists ( tmp_name1) && (F=main_read_aln(tmp_name1, NULL))!=NULL)
            {
              free_aln(F);
            lsqman_result=tmp_name1;
            }
       
       else
         {
           seq_file=vtmpnam (NULL);
           lsqman_result=tmp_name1;
           fp=vfopen (seq_file, "w");
           fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s1],(CL->S)->seq[s2] );
           vfclose (fp);
           sprintf ( command, "%s -pdb %s -pep %s > %s%s", LSQMAN_4_TCOFFEE, (CL->S)->name[s1], seq_file,get_cache_dir(), lsqman_result);
           
           while (!F)
             {
             my_system ( command);
             F=main_read_aln (lsqman_result, NULL);
             if ( !F)
               {
                 fprintf ( stderr, "\n\tlsqman failed: will be retried");
                 if ( n_failure==0)fprintf ( stderr, "\n\t%s", command);
                 n_failure++;
                 if ( n_failure==10)
                  {
                  fprintf ( stderr, "\nCould not run Fugue: will replace it with slow_pair\n");
                  vremove (lsqman_result);      
                  return NULL;
                  }
               }
             free_aln (F);
             }
           vremove( seq_file);

         }
       

       F=main_read_aln(lsqman_result, NULL);
       /*sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
       sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
       */
       F=fix_aln_seq (F, CL->S);
       




       CL=aln2constraint_list (F, CL, "100");
       free_aln (F);
       return CL;
      }

Constraint_list * sap_pair   (char *seq, Constraint_list *CL)
 {
            register int a;
          FILE *fp;
          char full_name[FILENAMELEN];
          char *tmp_pdb1, *tmp_pdb2;
          char *sap_seq1, *sap_seq2;
          char *sap_lib,  *tmp_name, *tmp_name1, *tmp_name2;
          int s1, s2, r1, r2, l;
          int sim=0, tot=0, score=0;
          char command[STRING];
          char program[STRING];
          char *string1, *string2, *string3, *string4, *string5;
          int max_struc_len=10000;
          char *template1, *template2;
          
          /*check_program_is_installed ( "sap"                   ,SAP_4_TCOFFEE, "SAP_4_TCOFFEE",MAIL, IS_FATAL);*/     
     
          

          atoi(strtok (seq,SEPARATORS));
          s1=atoi(strtok (NULL,SEPARATORS));
          s2=atoi(strtok (NULL,SEPARATORS));
          template1=seq2T_value(CL->S,s1, "template_name", "_P_");
          template2=seq2T_value(CL->S,s2, "template_name", "_P_");
          if (!template1 || !template2) return CL;
          
          declare_name (string1);
          declare_name (string2);
          declare_name (string3);
          declare_name (string4);
          declare_name (string5);
          
#ifndef     SAP_4_TCOFFEE
          if ( getenv ( "SAP_4_TCOFFEE")==NULL)crash ("SAP_4_TCOFFEE IS NOT DEFINED");
          else  sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
#else
          if ( getenv ( "SAP_4_TCOFFEE")==NULL)sprintf (program, "%s", SAP_4_TCOFFEE);
          else sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
#endif

            

          tmp_name1=vcalloc (100, sizeof (char));         
          sprintf ( tmp_name1, "%s_%s.sap_results",template1,template2);
          tmp_name2=vcalloc (100, sizeof (char));
          sprintf ( tmp_name2, "%s_%s.sap_results",template2,template1);
          
          
          if (is_sap_file (tmp_name1))
            {
            tmp_name=tmp_name1;
            }
          else if ( is_sap_file (tmp_name2))
            {
            tmp_name=tmp_name2;
            SWAP (s1, s2);
            }
          else
            {
            tmp_name=tmp_name1;
            


      
            tmp_pdb1=normalize_pdb_file(seq2P_template_file(CL->S,s1),(CL->S)->seq[s1], vtmpnam (NULL));
            tmp_pdb2=normalize_pdb_file(seq2P_template_file(CL->S,s2),(CL->S)->seq[s2], vtmpnam (NULL));
      
            
            sprintf ( command , "%s %s %s > %s%s",program,tmp_pdb1,tmp_pdb2, get_cache_dir(), tmp_name);     

            my_system  ( command);
      
            if ( flag_file2remove_is_on())
              {
                sprintf (full_name, "%s%s",get_cache_dir(), tmp_name);
                add2file2remove_list (full_name);
              }
            
            remove ("super.pdb");
            }

          

          sap_seq1=vcalloc (max_struc_len, sizeof (char));
          sap_seq2=vcalloc (max_struc_len, sizeof (char));
          
          fp=find_token_in_file ( tmp_name, NULL, "Percent");
          fp=find_token_in_file ( tmp_name, fp  , "Percent");
          while ( (fgetc (fp))!='\n');
          while ( fscanf (fp, "%s %s %s %s %s\n", string1, string3, string4, string5, string2)==5 && !strm ( string3, "Weighted") &&  !strm ( string1, "Weighted"))
                {
                
                  l=strlen (string1);
                  r1=tolower(string1[l-1]);
                  r2=tolower(string2[0]);
                  sim+=(r1==r2)?1:0;
                  if ( tot>max_struc_len)
                  {max_struc_len+=max_struc_len;
                       sap_seq1=vrealloc ( sap_seq1, sizeof(char)*max_struc_len);
                       sap_seq2=vrealloc ( sap_seq2, sizeof(char)*max_struc_len);
                  }
                        
                   sap_seq1[tot]=r1;
                   sap_seq2[tot]=r2;                   
                     tot++;
              }
          score=(sim*100)/tot;
          vfclose (fp);
          sap_seq1[tot]=sap_seq2[tot]='\0';
        
           
          fp=vfopen ( sap_lib=vtmpnam(NULL), "w");
          fprintf (fp, "! TC_LIB_FORMAT_01\n");
          fprintf (fp, "2\n");
          fprintf (fp, "%s %d %s\n", (CL->S)->name[s2],strlen (sap_seq1), sap_seq1);
          fprintf (fp, "%s %d %s\n", (CL->S)->name[s1],strlen (sap_seq2), sap_seq2);
          fprintf (fp, "#1 2\n");
          
          score=100;

          for ( a=0; a< tot; a++)
            {
            
            fprintf (fp, "%d %d %d 1 0\n", a+1, a+1, score);
            
            }
          fprintf (fp, "! CPU 0\n");
          fprintf (fp, "! SEQ_1_TO_N\n");
          vfclose (fp);
          CL=read_constraint_list_file(CL,sap_lib);
          vremove (sap_lib); 
          
          vfree ( string1);vfree ( string2);vfree ( string3);vfree ( string4);vfree ( string5);
          vfree (sap_seq1); vfree(sap_seq2);vfree (tmp_name1); vfree(tmp_name2); 
          
          return CL;
      }



Constraint_list * sap_pair_rs (char *seq, Constraint_list *CL)
        {
          FILE *fp;
          char *tmp_pdb1, *tmp_pdb2;
          char *string1, *tmp_name, *string2, *string3, *string4, *string5;
          float score;
          int s1, s2;
          int r1, r2,rs1, l, p1, p2;
          
          static CLIST_TYPE *entry;
          char command[STRING];
          char program[STRING];

          if ( !entry)entry=vcalloc ( LIST_N_FIELDS, sizeof ( CLIST_TYPE ));
          declare_name (string1);
          declare_name (string2);
          declare_name (string3);
          declare_name (string4);
          declare_name (string5);
          
        
          atoi(strtok (seq,SEPARATORS));
          s1=atoi(strtok (NULL,SEPARATORS));
          s2=atoi(strtok (NULL,SEPARATORS));
          

          tmp_name=vtmpnam ( NULL);

#ifndef     SAP_4_TCOFFEE
          if ( getenv ( "SAP_4_TCOFFEE")==NULL)crash ("SAP_4_TCOFFEE IS NOT DEFINED");
          else  sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
#else
          if ( getenv ( "SAP_4_TCOFFEE")==NULL)sprintf (program, "%s", SAP_4_TCOFFEE);
          else sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
#endif


          /*prepare the two files*/
          tmp_pdb1=normalize_pdb_file(seq2P_template_file(CL->S,s1),(CL->S)->seq[s1], vtmpnam (NULL));
          tmp_pdb2=normalize_pdb_file(seq2P_template_file(CL->S,s2),(CL->S)->seq[s2], vtmpnam (NULL));
          

          sprintf ( command , "%s %s %s > %s%s",program,tmp_pdb1,tmp_pdb2,get_cache_dir(), tmp_name);     
          my_system  ( command);
          
          
          rs1=2;
          fp=find_token_in_file ( tmp_name, NULL, "Percent");
          fp=find_token_in_file ( tmp_name, fp  , "Percent");
          while ( (fgetc (fp))!='\n');
          while ( fscanf (fp, "%s %s %s %s %s\n", string1, string3, string4, string5, string2)==5 && !strm ( string3, "Weighted") &&  !strm ( string1, "Weighted"))
                {
                
                  p1=atoi(string3);
                  score=atof(string4);
                  p2=atoi(string5);
                  l=strlen (string1);
                  r1=tolower(string1[l-1]);
                  r2=tolower(string2[0]);
                  
                  
                  if ( p1> (CL->S)->len[s2]){rs1=1; break;}
                  else if ( p2> (CL->S)->len[s1]){rs1=1; break;}
                  else if ( r1!=tolower ( (CL->S)->seq[s2][p1-1])){rs1=1; break;}
                  else if ( r2!=tolower ( (CL->S)->seq[s1][p2-1])){rs1=1; break;}
              }
          vfclose (fp);

          fp=find_token_in_file ( tmp_name, NULL, "Percent");
          fp=find_token_in_file ( tmp_name, fp  , "Percent");
          while ( (fgetc (fp))!='\n');

          
          
          while ( fscanf (fp, "%s %s %s %s %s\n", string1, string3, string4, string5, string2)==5 && !strm ( string3, "Weighted") && !strm ( string1, "Weighted"))
                {
                  

                  p1=(rs1==1)?atoi(string3):atoi(string5);
                  score=atof(string4);
                  p2=(rs1==1)?atoi(string5):atoi (string3);

                              
                  entry[SEQ1]=s1;
                  entry[SEQ2]=s2;
                  entry[R1]=p1;
                  entry[R2]=p2;
                  entry[WE]=(int)score+1;
                  entry[CONS]=1;
                  entry[MISC]=0;
                  CL=add_entry2list(entry, CL);
                  
              }  

          vfclose (fp);

          vfree ( string1);
          vfree ( string2);

          vremove (tmp_pdb1);
          vremove (tmp_pdb2);
          vremove ("super.pdb");

          return CL;
      }

Constraint_list * thread_pair   (char *in_seq, Constraint_list *CL)
 {
   char seq[STRING];
   int s1, s2;
   Alignment* A1=NULL, *A2=NULL, *A3=NULL;
   static char *log_name;

   FILE *fp;
   
   if ( !CL) return NULL;

   if ( CL->method_log[0] && !log_name)
     {
       declare_name(log_name);
       sprintf (log_name, "%s.thread_pair", CL->method_log);
       fp=vfopen ( log_name, "w");
     }
   else if ( CL->method_log[0]);
     {
       fp=vfopen ( log_name, "a");
     }
  
   
   /*1: Identify the two sequences s1 and s2 that are being considered*/
      
   sprintf ( seq, "%s", in_seq);
   atoi(strtok (seq,SEPARATORS));
   s1=atoi(strtok (NULL,SEPARATORS));
   s2=atoi(strtok (NULL,SEPARATORS));
   
   /*2: Replace each sequence with its target structure*/   
   A1=seq2pdb ( (CL->S)->name[s1],(CL->S)->seq[s1],CL->Pdb_Blast);
   A2=seq2pdb ( (CL->S)->name[s2],(CL->S)->seq[s2],CL->Pdb_Blast);
   
      
   /*3: measure the similarity between the two sequences*/
   if (A1 && A2)
     {
       A3=align_two_sequences ((CL->S)->seq[s1],(CL->S)->seq[s2],"blosum62mt",-4,-1, "myers_miller_pair_wise");
       A3->score_aln=aln2sim(A3, "idmat");
     }

   if ( CL->method_log[0])
     {
       fprintf ( fp, "S:  %s\t", (CL->S)->name[s1]);
       fprintf ( fp, "T:  %s\t", (A1)?A1->name[1]:"NO");
       fprintf ( fp, "id: %d\t", (A1)?A1->score_aln:0);
       fprintf ( fp, "co: %d\t", (A1)?aln2coverage(A1,0):0);
       fprintf ( fp, "S:  %s\t", (CL->S)->name[s2]);
       fprintf ( fp, "T:  %s\t", (A2)?A2->name[1]:"NO");
       fprintf ( fp, "id: %d\t", (A2)?A2->score_aln:0);
       fprintf ( fp, "co: %d\t", (A2)?aln2coverage(A2,0):0);
     }
       
   /*4: if one of the structure has no target, or if the targets are less similar than the structures: do a sequence analysis*/
   
   CL=sap_pair     ( in_seq, CL);

#ifdef XXXXX
   if (!A3 || (MIN(A1->score_aln, A2->score_aln)<A3->score_aln))
     
     {
       CL=seq2list     ( CL,"lalign_id_pair", in_seq, "sim");
       CL=seq2list     ( CL,"slow_pair", in_seq, "sim");
       if ( CL->method_log[0])fprintf ( fp, "m: no_sap\n");

     }
   else 
     {
       
       
       CL=seq2list     ( CL,"lalign_id_pair", in_seq, "sim");
       set_file2remove_on();/*Make sure all the sap files will be removed when finished*/
       CL=sap_pair     ( in_seq, CL);
       set_file2remove_off();
        
       if ( CL->method_log[0])fprintf ( fp, "m: sap\n");
     }
#endif

   if ( CL->method_log[0])vfclose(fp);
   free_sequence (free_aln (A1), -1);
   free_sequence (free_aln (A2), -1);
   free_sequence (free_aln (A3), -1);
   
   CL=compact_list (CL, 0, CL->ne, "default");
   return CL;
 }
/******************************************************************/
/*                   DATABASE METHODS                             */
/*                                                                */
/*                                                                */
/******************************************************************/

/******************************************************************/
/*                   GENERIC PAIRWISE METHODS                     */
/*                                                                */
/*                                                                */
/******************************************************************/


Alignment * fast_pair      (Job_TC *job)
        {
          int s, n,a;       
          static int **l_s;
          static int *ns;
          char seq[1000];
          Alignment *A;
          Constraint_list *PW_CL;
          Constraint_list *CL;
          char *seq_in;
          Sequence *S;
          TC_method *M;
          int*seqlist;
          char **buf;
          
          A=(job->io)->A;
          
          M=(job->param)->TCM;
          PW_CL=((job->param)->TCM)->PW_CL;
          CL=(job->io)->CL;
          seq_in=(job->param)->seq_c;
          

          sprintf (seq, "%s", seq_in);
          seqlist=string2num_list (seq);
          n=seqlist[1];
          if ( n!=2){fprintf ( stderr, "\nERROR: fast_pw_aln can only handle two seq at a time [FATAL]\n");exit (EXIT_FAILURE);}
          
          S=(CL)->S;
          
          if (!A) {A=declare_aln (CL->S);}
          if ( !ns)
              {
            ns=vcalloc ( 2, sizeof (int));
            l_s=declare_int (2,(CL->S)->nseq);
            }
          buf=vcalloc ( S->nseq, sizeof (char*));
          
          for ( a=0; a< n; a++)
              {
              s=seqlist[a+2];
              if ( strm (M->seq_type, "G"))
                {
                  buf[s]=S->seq[s];
                  S->seq[s]=((((S->T[s])->G)->VG)->S)->seq[0];
              }
            else
              buf[s]=S->seq[s];
            

              sprintf ( A->seq_al[a], "%s",S->seq[s]);
              sprintf ( A->name[a], "%s", (CL->S)->name[s]);
              A->order[a][0]=s;
            }
        
          A->S=CL->S;
          PW_CL->S=CL->S;
          
          ns[0]=ns[1]=1;
          l_s[0][0]=0;
          l_s[1][0]=1;
       
          pair_wise ( A, ns, l_s, PW_CL);
          A->nseq=n;
          
          for ( a=0; a<S->nseq; a++)
            {
            if ( !buf[a] || buf[a]==S->seq[a]);
            else S->seq[a]=buf[a];
            }
          vfree (buf);vfree (seqlist);
          return A;
          
      }
Alignment * align_two_aln ( Alignment *A1, Alignment  *A2, char *in_matrix, int gop, int gep, char *in_align_mode)
        {
      Alignment *A=NULL;
      Constraint_list *CL;
      Sequence *S;
      int a;
      int *ns;
      int **ls;
      static char *matrix;
      static char *align_mode;
      
      if (!matrix)matrix=vcalloc ( 100, sizeof (char));
      if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
      sprintf ( matrix, "%s", in_matrix);
      sprintf ( align_mode, "%s", in_align_mode);
      
      CL=vcalloc ( 1, sizeof (Constraint_list));
      CL->pw_parameters_set=1;
      CL->M=read_matrice (matrix);
      CL->matrices_list=declare_char (10, 10);

      CL->evaluate_residue_pair=evaluate_matrix_score;
      CL->get_dp_cost=super_fast_get_dp_cost2;
      CL->normalise=1;
      
      CL->extend_jit=0;
      CL->maximise=1;
      CL->gop=gop;
      CL->gep=gep;
      CL->TG_MODE=2;
      sprintf (CL->matrix_for_aa_group, "vasiliky");
      CL->use_fragments=0;
      CL->ktup=5;
      if ( !CL->use_fragments)CL->diagonal_threshold=0;
      else CL->diagonal_threshold=6;

      sprintf (CL->dp_mode, "%s", align_mode);

      A=copy_aln (A1, A);
      A=stack_aln (A, A2);
      CL->S=fill_sequence_struc(A->nseq, A->seq_al,A->name);
      
      ns=vcalloc ( 2, sizeof(int));
      ls=declare_int ( 2,A->nseq);
      ns[0]=A1->nseq;
      ns[1]=A2->nseq;
      for ( a=0; a<ns[0]; a++)
        ls[0][a]=a;
      for ( a=0; a< ns[1]; a++)
        ls[1][a]=a+A1->nseq;
      
      A->score_aln=pair_wise (A, ns, ls,CL);
      
      vfree (ns);
      free_int (ls, -1);
      S=free_constraint_list (CL);
      free_sequence (S,-1);
      A->S=NULL;
      return A;
      }

static int align_two_seq_keep_case;
void toggle_case_in_align_two_sequences(int value)
{
  align_two_seq_keep_case=value;
}
Alignment * align_two_sequences ( char *seq1, char *seq2, char *in_matrix, int gop, int gep, char *in_align_mode)
        {
      static Alignment *A;
      Constraint_list *CL;
      Sequence *S;
      
      int *ns;
      int **l_s;
      
      char       **seq_array;
      char       **name_array;
      static char *matrix;
      static int  **M;
      
      static char *align_mode;
      
      if (!matrix)matrix=vcalloc ( 100, sizeof (char));
      if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
      sprintf ( align_mode, "%s", in_align_mode);
      
      CL=vcalloc ( 1, sizeof (Constraint_list));
      
      CL->pw_parameters_set=1;
      
      CL->matrices_list=declare_char (10, 10);


      if ( !strm (matrix, in_matrix))
        {
          sprintf ( matrix,"%s", in_matrix);
          M=CL->M=read_matrice (matrix);
         
        }
      else
        {
          CL->M=M;
        }
      
      if (strstr (in_align_mode, "cdna")) 
        CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
      else
        CL->evaluate_residue_pair=evaluate_matrix_score;
      
      CL->get_dp_cost=get_dp_cost;
      CL->extend_jit=0;
      CL->maximise=1;
      CL->gop=gop;
      CL->gep=gep;
      CL->TG_MODE=2;
      sprintf (CL->matrix_for_aa_group, "vasiliky");
      CL->use_fragments=0;
      CL->ktup=3;
      if ( !CL->use_fragments)CL->diagonal_threshold=0;
      else CL->diagonal_threshold=6;

      sprintf (CL->dp_mode, "%s", align_mode);

      seq_array=declare_char ( 2, MAX(strlen(seq1), strlen (seq2))+1);
      sprintf (seq_array[0], "%s",seq1);
      sprintf (seq_array[1],"%s", seq2);
      ungap_array(seq_array,2);
      if (align_two_seq_keep_case !=KEEP_CASE)string_array_lower(seq_array,2);

      name_array=declare_char (2, STRING);
      sprintf ( name_array[0], "A");
      sprintf ( name_array[1], "B");
        
      
      ns=vcalloc ( 2, sizeof(int));
      l_s=declare_int ( 2, 1);
      ns[0]=ns[1]=1;
      l_s[0][0]=0;
      l_s[1][0]=1;
      
      
      
      CL->S=fill_sequence_struc(2, seq_array, name_array);

      A=seq2aln(CL->S, NULL, 1);

      ungap (A->seq_al[0]);
      ungap (A->seq_al[1]);



      A->score_aln=pair_wise (A, ns, l_s,CL);

      vfree (ns);
      free_int (l_s, -1);
      free_char (name_array, -1);free_char ( seq_array,-1);
      
      CL->M=NULL;
        S=free_constraint_list (CL);
      free_sequence (S,-1);
      A->S=NULL;
      return A;
      }

NT_node ** make_tree ( Alignment *A,Constraint_list *CL,int gop, int gep,Sequence *S,  char *tree_file, char *tree_mode, int maximise)
      {
      int a, b, ra, rb;
      NT_node **T;
      int **distances;
      int out_nseq;
      char **out_seq_name;
      char **out_seq;
      
      
      out_nseq=S->nseq;
      out_seq_name=S->name;
      out_seq=S->seq;
      
      
      get_constraint_list_similarity_matrix (CL);

      if ( CL->S!=S)
        {
          /*Shrink the distance matrix so that it only contains the required sequences*/
          distances=declare_int (S->nseq, S->nseq);
          for (a=0; a< S->nseq; a++)
            {
            ra=name_is_in_list ((S)->name[a],(CL->S)->name, (CL->S)->nseq, 100);
            for ( b=0; b< S->nseq; b++)
            {
              rb=name_is_in_list ((S)->name[b],(CL->S)->name, (CL->S)->nseq, 100);
              distances[a][b]=CL->score_similarity_matrix[ra][rb];
            }
            }
        }
      else
        {
          distances=duplicate_int ( CL->score_similarity_matrix, -1, -1);
        }

      

      distances=sim_array2dist_array (distances, 100);
            
      if ( strm (tree_mode, "order"))
         {
           for ( a=0; a< S->nseq; a++)
             for ( b=0; b< S->nseq; b++)
             distances[b][a]=100;
         }
      
                  
      
      fprintf (CL->local_stderr , "\nMAKE NEIGHBOR JOINING DENDROGRAM\n\t[MODE=%s][",tree_mode);
      
      
      T=make_nj_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, tree_mode);
      fprintf (CL->local_stderr , "DONE]\n");
      free_int (distances, out_nseq);
      return T;
      }


Alignment *recompute_local_aln (Alignment *A, Sequence *S,Constraint_list *CL, int scale, int gep)
    {
    int **coor;
    int a;
    Alignment *B;
   
    sort_constraint_list (CL, 0, CL->ne);
    coor=declare_int (A->nseq, 3);
    for ( a=0; a< A->nseq; a++)
        {
        coor[a][0]=A->order[a][0];
      coor[a][1]=A->order[a][1]+1;
      coor[a][2]=strlen(S->seq[A->order[a][0]])-coor[a][1];
      }
    B=build_progressive_nol_aln_with_seq_coor(CL,0,0,S,coor,A->nseq);
    A=copy_aln ( B, A);
 
    CL=compact_list (CL, 0,CL->ne, "shrink");
    free_Alignment(B);
    return A;
    }
    

Alignment *build_progressive_nol_aln_with_seq_coor(Constraint_list *CL,int gop, int gep,Sequence *S, int **seq_coor, int nseq)
    {

    static int ** local_coor1;
    static int ** local_coor2;
    if ( local_coor1!=NULL)free_int (local_coor1, -1);
    if ( local_coor2!=NULL)free_int (local_coor2, -1);

    local_coor1=get_nol_seq          ( CL,seq_coor, nseq, S);
    local_coor2=minimise_repeat_coor ( local_coor1, nseq, S);

    return build_progressive_aln_with_seq_coor(CL,gop, gep,S, local_coor2,nseq);
    }


Alignment *build_progressive_aln_with_seq_coor (Constraint_list*CL,int gop, int gep, Sequence *S, int **coor, int nseq)
    {
    Alignment *A=NULL;

    A=seq_coor2aln (S,NULL, coor, nseq);
   
    return build_progressive_aln ( A,CL, gop, gep);
    }

Alignment *est_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
    {
    int a,n;
    int**group_list; 
    int *n_groups;
    char *seq;
    n_groups=vcalloc ( 2, sizeof (int));
    group_list=declare_int ( 2, A->nseq);
    
    n=A->nseq;

    n_groups[0]=1;
    n_groups[1]=1; 
    group_list[0][0]=0;
    group_list[0][1]=1;
    
    group_list[1][0]=1;
    fprintf ( stderr, "\n");  
    for ( a=1; a<n; a++)
        {
      sprintf ( A->seq_al[1], "%s", A->seq_al[a]);
      fprintf ( stderr, "\t[%30s]->[len=%5d]", A->name[a],strlen ( A->seq_al[0]));  
      pair_wise ( A,n_groups, group_list, CL);
      
      seq=dna_aln2cons_seq(A);
      
      sprintf ( A->seq_al[0], "%s", seq);
      vfree (seq);
      fprintf ( stderr, "\n");  
      }
   
    A->nseq=1;
    return A;
    }

void analyse_seq ( Alignment *A, int s)
   {
     int a, b, c;
     int r;

     int len=0;


     int state;
     int pstate=-1;
     float score=0;
     
     for ( a=0; a< A->len_aln; a++)
       {
       for ( b=0, c=0; b< s; b++)
         if ( !is_gap(A->seq_al[b][a])){c=1; break;}
       
       r=!is_gap(A->seq_al[s][a]);
       
       if      (  r &&  c) state=1;
       else if ( !r && !c) state=2;
       else if ( !r &&  c) state=3;
       else if (  r && !c) state=4;
       
       if ( state !=pstate)
         {
           score+=len*len;
           len=0;
         }
       len+=r;
       pstate=state;
       }
     score=score/(float)(((A->S)->len[s]*(A->S)->len[s]));
     fprintf ( stderr, "[%.2f]", score);
     
     return;
   }
Alignment * tsp_aln (Alignment *A, Constraint_list *CL, Sequence *S)
{
  int a, b   ;
  int **  distances;
  int *ns, **ls;
  int **used;
      
  A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
  ns=vcalloc (2, sizeof (int));
  ls=declare_int (2, (CL->S)->nseq);
  used=declare_int ( A->nseq, 2);
  
  
  get_constraint_list_similarity_matrix (CL);  
  distances=declare_int (A->nseq+1, A->nseq+1);
  distances=duplicate_int ( CL->score_similarity_matrix, -1, -1);

  for (a=0; a< A->nseq; a++)
    {
      used[a][0]=a;
      for (b=0; b< A->nseq; b++)
      {
        used[a][1]+=distances[a][b];
      }
    }

  sort_int_inv (used,2,1,0,(CL->S)->nseq-1);

  ls[0][ns[0]++]=used[0][0];
  ns[1]=1;
  
  for (a=1; a< S->nseq; a++)
    {
      fprintf ( stderr, "\n%s %d", (CL->S)->name[used[a][0]], used[a][1]);
      ls[1][0]=used[a][0];
      pair_wise ( A,ns,ls, CL);
      ls[0][ns[0]++]=used[a][0];
    }

  A->nseq=(CL->S)->nseq;
  return A;

}

Alignment *build_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
    {
    int a,n;
    int**group_list; 
    int *n_groups;
    char dp_mode[100];


    sprintf ( dp_mode, "%s", CL->dp_mode);
    sprintf (CL->dp_mode, "gotoh_pair_wise");

    n_groups=vcalloc ( 2, sizeof (int));
    group_list=declare_int ( 2, A->nseq);
    
    n=A->nseq;
    
    for ( a=0; a<n; a++)ungap(A->seq_al[a]);
    for ( a=1; a<n; a++)
        {
      n_groups[0]=a;
      n_groups[1]=1;
      group_list[0][a-1]=a-1;
      group_list[1][0]  =a;
      
      pair_wise ( A,n_groups, group_list, CL);
      fprintf ( stderr, "\n\t[%d]->[%d]", a,strlen ( A->seq_al[0]));
      }
    fprintf (stderr, "\n");
    vfree(n_groups);
    free_int ( group_list, -1);
    sprintf (CL->dp_mode, "%s",dp_mode);

    return A;
    }

Alignment *iterate_aln ( Alignment*A, int nit, Constraint_list *CL)
{
  int it;
  Constraint_list *B;

  B=vcalloc ( 1, sizeof (Constraint_list));
  B->evaluate_residue_pair=CL->evaluate_residue_pair;
  B->get_dp_cost= CL->get_dp_cost;
  B->pair_wise=CL->pair_wise;
  
  fprintf ( CL->local_stderr, "Iterated Refinement: %d cycles\n", nit);
  choose_extension_mode ( "slow_triplet",CL);

  
  if ( nit==-1)nit=A->nseq*2;
  if ( A->len_aln==0)A=very_fast_aln (A, A->nseq, CL);
  for (it=0; it< nit; it++)
    {
      CL->local_stderr=output_completion (CL->local_stderr,it, nit,1, "");
      A=realign_aln (A, CL);
    }
  fprintf ( CL->local_stderr, "\nIterated Refinement: Completed\n");
  CL->evaluate_residue_pair=B->evaluate_residue_pair;
  CL->get_dp_cost=B->get_dp_cost;
  CL->pair_wise=B->pair_wise;
  
  vfree (B);

  return A;
}

Alignment *realign_aln ( Alignment*A, Constraint_list *CL)
{
  int *ns;
  int **l_s;
  int **s, **t;
  int n;
  int  a,b;
  
  s=get_sim_aln_array ( A, "idmat");
  t=declare_int (A->nseq, 2);
  for (a=0; a< A->nseq; a++)
    {
      t[a][0]=a;
      for ( b=0; b< A->nseq; b++)
      t[a][1]+=s[a][b];
    }
  sort_int_inv (t,2,1,0, A->nseq-1);
  for (a=0; a< A->nseq; a++)
    fprintf ( stderr, "\n%3d %3d", t[a][0], t[a][1]/A->nseq);
  
  ns=vcalloc (2, sizeof (int));
  l_s=declare_int (2,A->nseq);
  
  n=A->nseq;
  for ( a=n/2; a<n; a++)
    {
      
      ns[0]=ns[1]=0;
      for (b=0; b<a      ; b++)l_s[0][ns[0]++]=t[b][0];
      for (b=a; b<=a     ; b++)if ( b!=a)l_s[1][ns[1]++]=t[b][0];
      
      fprintf ( stderr, "\n%d %s",t[a][0], A->seq_al[t[a][0]]);
      /*ungap_sub_aln ( A, ns[0], l_s[0]);*/
      ungap_sub_aln ( A, ns[1], l_s[1]);
      
      A->score_aln=pair_wise (A, ns, l_s,CL);
      fprintf ( stderr, "\n%d %s",t[a][0], A->seq_al[t[a][0]]);
      
    }
  
  /*
  for ( a=0; a< A->nseq; a++)
    {
      
      ns[0]=ns[1]=0;
      for (b=a; b<=a      ; b++)l_s[0][ns[0]++]=t[b][0];
      for (b=0; b<A->nseq; b++)if ( b!=a)l_s[1][ns[1]++]=t[b][0];
      
      fprintf ( stderr, "\n%s", A->seq_al[t[a][0]]);
      ungap_sub_aln ( A, ns[0], l_s[0]);
      ungap_sub_aln ( A, ns[1], l_s[1]);
      
      A->score_aln=pair_wise (A, ns, l_s,CL);
      fprintf ( stderr, "\n%s", A->seq_al[t[a][0]]);
      
    }
  */
  
  vfree(ns);free_int(l_s, -1);
  return A;
}
Alignment *realign_aln2 ( Alignment*A, Constraint_list *CL)
{
  int *ns;
  int **l_s;
  
  int  a,g;
  
  ns=vcalloc (2, sizeof (int));
  l_s=declare_int (2,A->nseq);
  
  for ( a=0; a< A->nseq; a++)
    {
    g=rand()%2;
    l_s[g][ns[g]++]=a;
    }
 
  for ( a=0; a<2; a++)
    {
      int b;
      fprintf ( stderr, "\nG%d:", a);
      for ( b=0; b< ns[a]; b++)
      fprintf ( stderr, " [%d]", l_s[a][b]);
    }
  fprintf ( stderr, "\n");
  ungap_sub_aln ( A, ns[0], l_s[0]);
  ungap_sub_aln ( A, ns[1], l_s[1]);
  
  
  A->score_aln=pair_wise (A, ns, l_s,CL);
  
  vfree(ns);free_int(l_s, -1);
  return A;
}

      
      
    
Alignment *very_fast_aln ( Alignment*A, int nseq, Constraint_list *CL)
{
char command[10000];
char *tmp_seq;
char *tmp_aln;
FILE *fp;

if ( CL && CL->local_stderr)fp=CL->local_stderr;
else fp=stderr; 
 
 fprintf (fp, "\n[Computation of an Approximate MSA...");
 tmp_seq= vtmpnam (NULL);
 tmp_aln= vtmpnam (NULL);
 output_fasta_seq ((tmp_seq=vtmpnam (NULL)), A);
 sprintf ( command, "t_coffee -infile=%s -very_fast -outfile=%s %s -outorder=input", tmp_seq, tmp_aln, TO_NULL_DEVICE);
 my_system ( command);
 A->nseq=0;
 A=main_read_aln (tmp_aln,A);
 fprintf (fp, "]\n");
 return A;
} 

static NT_node* SNL;
NT_node* tree_aln ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
    {
    int *n_s;
    int ** l_s;
    int a, b, c;
    int score;
    static int n_groups_done;
    int *list_seq, nseq2align;
    int *translation;
    int do_split;
    
    NT_node P;
    Alignment *SUB1=NULL, *SUB2=NULL, *SUB3=NULL;
    
      
        

    if (n_groups_done==0)
       {
       if (SNL)vfree(SNL);
       SNL=vcalloc ( (CL->S)->nseq, sizeof (NT_node));
       
       if (CL->translation)vfree(CL->translation);
       CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
       
       for ( a=0; a< (CL->S)->nseq; a++)
         CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, 100);
       if (nseq>2)fprintf ( CL->local_stderr, "\nPROGRESSIVE_ALIGNMENT [Tree Based]\n");
       else fprintf ( CL->local_stderr, "\nPAIRWISE_ALIGNMENT [No Tree]\n");
       n_groups_done=(CL->S)->nseq;
       A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
       A->nseq=nseq;
       }

    translation=CL->translation;    
    n_s=vcalloc (2, sizeof ( int));
    l_s=declare_int ( 2, nseq);

    
    if ( nseq==2)
       {
       n_s[0]=n_s[1]=1;
       l_s[0][0]=name_is_in_list ((CL->S)->name[0],(CL->S)->name, (CL->S)->nseq, 100);
       l_s[1][0]=name_is_in_list ((CL->S)->name[1],(CL->S)->name, (CL->S)->nseq, 100);
       A->score_aln=score=pair_wise (A, n_s, l_s,CL);
       
       vfree ( n_s);
       free_int ( l_s, 2);
       return SNL;
       }
    else
       {
       if ( RT->parent !=LT->parent)fprintf ( stderr, "Tree Pb [FATAL:%s]", PROGRAM);
       else P=RT->parent;
       
       if ( LT->leaf==1 && RT->leaf==0)
         tree_aln ( RT->left, RT->right,A, nseq,CL);

       else if ( RT->leaf==1 && LT->leaf==0)
         tree_aln ( LT->left, LT->right,A,nseq,CL);
 
       else if (RT->leaf==0 && LT->leaf==0)
          {
        tree_aln ( LT->left, LT->right,A,nseq,CL);
        tree_aln ( RT->left, RT->right,A,nseq,CL);
        }

       if ( LT->leaf==1 && RT->leaf==1)
          {   
        /*1 Identify the two groups of sequences to align*/
        
        nseq2align=LT->nseq+RT->nseq;
        n_s[0]=LT->nseq;
        for ( a=0; a< LT->nseq; a++)l_s[0][a]=translation[LT->lseq[a]];
        if ( LT->nseq==1)LT->group=l_s[0][0];
       
        n_s[1]=RT->nseq;
        for ( a=0; a< RT->nseq; a++)l_s[1][a]=translation[RT->lseq[a]];
        if ( RT->nseq==1)RT->group=l_s[1][0];
        

        list_seq=vcalloc (nseq2align, sizeof (int));
        for (b=0, a=0; a< 2; a++){for (c=0;c< n_s[a]; c++)list_seq[b++]=l_s[a][c];}
        
        
        P->group=n_groups_done++;
        fprintf (CL->local_stderr, "\n\tGroup %3d: [Group %3d (%3d seq)] with [Group %3d (%3d seq)]-->",P->group+1,RT->group+1, n_s[1],LT->group+1, n_s[0]);
        
        P->score=A->score_aln=score=pair_wise (A, n_s, l_s,CL);
        if (getenv4debug ("DEBUG_MALN")){fprintf (stderr, "\nDEBUG_MALN:tree_aln]");print_aln (A);}
        if (getenv4debug ("DEBUG_MALN")){fprintf (stderr, "\nEvaluation mode: %s", CL->evaluate_mode);}

        if (!strm (CL->evaluate_mode, "no"))
            { 
            SUB1=copy_aln (A, NULL);        
            SUB2=extract_sub_aln (SUB1,n_s[0]+n_s[1],list_seq);
            
            SUB3=main_coffee_evaluate_output (SUB2,CL,CL->evaluate_mode);
            score=SUB3->score_aln;
          
            free_aln (SUB1);SUB1=NULL;
            free_aln (SUB2);SUB2=NULL;
            free_aln (SUB3);SUB3=NULL;
            }
        
        fprintf (CL->local_stderr, "[Score=%4d][Len=%5d]",score, strlen ( A->seq_al[l_s[0][0]]));
        
        if (nseq2align==nseq) 
            {
            fprintf (CL->local_stderr, "\n\n");
            
            for (b=0, a=0; a< n_s[0]; a++, b++)sprintf ( A->tree_order[b],"%s", (CL->S)->name[l_s[0][a]]);
            for (a=0; a< n_s[1]     ; a++, b++)sprintf ( A->tree_order[b], "%s",(CL->S)->name[l_s[1][a]]);
            n_groups_done=0;
            }
        }
       if (P->parent)P->leaf=1;
       if ( LT->isseq==0)LT->leaf=0;
       if ( RT->isseq==0)RT->leaf=0;
              
       if (RT->isseq){SNL[translation[RT->lseq[0]]]=RT;RT->score=100;}
       if (LT->isseq){SNL[translation[LT->lseq[0]]]=LT;LT->score=100;}
       
       do_split=split_condition (nseq2align,A->score_aln,CL);
       if (CL->split && do_split)
       {
           
         for (a=0; a< P->nseq; a++)SNL[CL->translation[P->lseq[a]]]=NULL;
         SNL[CL->translation[RT->lseq[0]]]=P;
         
       }
       
       
       vfree ( n_s);
       free_int ( l_s, 2);
       return SNL;
       }
    
    }
int split_condition (int nseq, int score, Constraint_list *CL)
{
  int cond1=1, cond2=1;


  if ( CL->split_nseq_thres)cond1 =(nseq<=CL->split_nseq_thres)?1:0;
  if ( CL->split_score_thres)cond2=(score>=CL->split_score_thres)?1:0;

  return (cond1 && cond2);
}
int pair_wise (Alignment *A, int*ns, int **l_s,Constraint_list *CL )
    {
      /*
       CL->maximise
       CL->gop;
       CL->gep
       CL->TG_MODE;
      */
      int score;
      static int tag;
      int glocal;
      Pwfunc function;
      
      A->random_tag=++tag;
      
      if (! CL->pw_parameters_set)
               {
                   fprintf ( stderr, "\nERROR pw_parameters_set must be set in pair_wise [FATAL]\n" );crash("");
               }


      function=get_pair_wise_function(CL->pair_wise,  CL->dp_mode,&glocal); 
      if ( CL->get_dp_cost==NULL)CL->get_dp_cost=get_dp_cost;     
      
      if (strlen ( A->seq_al[l_s[0][0]])==0 || strlen ( A->seq_al[l_s[1][0]])==0)
        score=empty_pair_wise ( A, ns, l_s, CL, glocal);
      else
        score=function ( A, ns, l_s, CL);
      return score;
    }
      
int empty_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int glocal)
{
  int n, a, b;
  int *l;
  char *string;
  int l0, l1, len;
  
  if ( glocal==GLOBAL)
    {
      l0=strlen (A->seq_al[l_s[0][0]]);
      l1=strlen (A->seq_al[l_s[1][0]]);
      len=MAX(l1,l0);
      
      if ( len==0)return 0;
      else if (l0>l1){n=ns[1];l=l_s[1];}
      else if (l0<l1){n=ns[0];l=l_s[0];}
      string=generate_null (len);
      for ( a=0; a< n; a++)
      sprintf ( A->seq_al[l[a]], "%s", string);
      A->score=A->score_aln=0;
      A->len_aln=len;
      vfree ( string);
      return 0;
    }
  else if ( glocal==LALIGN)
    {
      A->A=declare_aln (A->S);
      (A->A)->len_aln=0;
      for ( a=0; a< 2; a++)
      for ( b=0; b<ns[a]; b++)
        A->seq_al[l_s[a][b]][0]='\0';
      (A->A)->score_aln=(A->A)->score=0;
      return 0;
    }
  else return 0;
}
  
      
      
      
Pwfunc get_pair_wise_function (Pwfunc pw,char *dp_mode, int *glocal)
  {
    /*Returns a function and a mode (Glogal, Local...)*/
      


    int a;
    static int npw;
    static Pwfunc *pwl;
    static char **dpl;
    static int *dps;

    /*The first time: initialize the list of pairwse functions*/
    if ( npw==0)
      {
      pwl=vcalloc ( 100, sizeof (Pwfunc));
      dpl=declare_char (100, 100);
      dps=vcalloc ( 100, sizeof (int));
      
      pwl[npw]=fasta_cdna_pair_wise;
      sprintf (dpl[npw], "fasta_cdna_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=cfasta_cdna_pair_wise;
      sprintf (dpl[npw], "cfasta_cdna_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=gotoh_pair_wise;
      sprintf (dpl[npw], "gotoh_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=myers_miller_pair_wise;
      sprintf (dpl[npw], "myers_miller_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=fasta_gotoh_pair_wise;
      sprintf (dpl[npw], "fasta_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      pwl[npw]=cfasta_gotoh_pair_wise;
      sprintf (dpl[npw], "cfasta_pair_wise");
      dps[npw]=GLOBAL;
      npw++;

      pwl[npw]=very_fast_gotoh_pair_wise;
      sprintf (dpl[npw], "very_fast_pair_wise");
      dps[npw]=GLOBAL;
      npw++;

      pwl[npw]=gotoh_pair_wise_sw;
      sprintf (dpl[npw], "gotoh_pair_wise_sw");
      dps[npw]=LOCAL;
      npw++;
      
      pwl[npw]=cfasta_gotoh_pair_wise_sw;
      sprintf (dpl[npw], "cfasta_sw_pair_wise");
      dps[npw]=LOCAL;
      npw++;
      
      pwl[npw]=gotoh_pair_wise_lalign;
      sprintf (dpl[npw], "gotoh_pair_wise_lalign");
      dps[npw]=LALIGN;
      npw++;
      
      pwl[npw]=sim_pair_wise_lalign;
      sprintf (dpl[npw], "sim_pair_wise_lalign");
      dps[npw]=LALIGN;
      npw++;
      
      pwl[npw]=domain_pair_wise;
      sprintf (dpl[npw], "domain_pair_wise");
      dps[npw]=MOCCA;
      npw++;
    
      pwl[npw]=gotoh_pair_wise;
      sprintf (dpl[npw], "ssec_pair_wise");
      dps[npw]=GLOBAL;
      npw++;

      pwl[npw]=ktup_pair_wise;
      sprintf (dpl[npw], "ktup_pair_wise");
      dps[npw]=LOCAL;
      npw++;
      
      pwl[npw]=precomputed_pair_wise;
      sprintf (dpl[npw], "precomputed_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=myers_miller_pair_wise;
      sprintf (dpl[npw], "default");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=viterbi_pair_wise;
      sprintf (dpl[npw], "viterbi_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=viterbiL_pair_wise;
      sprintf (dpl[npw], "viterbiL_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=viterbiD_pair_wise;
      sprintf (dpl[npw], "viterbiD_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=seq_viterbi_pair_wise;
      sprintf (dpl[npw], "seq_viterbi_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      
      pwl[npw]=pavie_pair_wise;
      sprintf (dpl[npw], "pavie_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      

      
      /*
      pwl[npw]=viterbiDGL_pair_wise;
      sprintf (dpl[npw], "viterbiDGL_pair_wise");
      dps[npw]=GLOBAL;
      npw++;
      */
      }

    for ( a=0; a< npw; a++)
      {
      if ( (dp_mode && strm (dpl[a], dp_mode)) || pwl[a]==pw)
           {
             pw=pwl[a];
             if (dp_mode)sprintf (dp_mode,"%s", dpl[a]);
             glocal[0]=dps[a];
             return pw;
           }
      }
    fprintf ( stderr, "\n[%s] is an unknown mode for pair_wise[FATAL]\n", dp_mode);
    crash ( "\n");
    return NULL;
  }
 

/*******************************************************************************/
/*                                                                             */
/*                                                                             */
/*    Util Functions                                                         */
/*                                                                             */
/*                                                                           */
/*******************************************************************************/

char *build_consensus ( char *seq1, char *seq2, char *dp_mode)
        {
      Alignment *A;
      char *buf;
      int a;
      char c1, c2;
      static char *mat;

      
      if ( !mat) mat=vcalloc ( STRING, sizeof (char));


      A=align_two_sequences (seq1, seq2, strcpy(mat,"idmat"), 0, 0,dp_mode);
      buf=vcalloc ( A->len_aln+1, sizeof (char));

      for ( a=0; a< A->len_aln; a++)
          {
            c1=A->seq_al[0][a];
            c2=A->seq_al[1][a];
            if (is_gap(c1) && is_gap(c2))buf[a]='-';
            else if (is_gap(c1))buf[a]=c2;
            else if (is_gap(c2))buf[a]=c1;
            else if (c1!=c2){vfree (buf);buf=NULL;free_aln(A);return NULL;}
            else buf[a]=c1;
          }
      buf[a]='\0';
      free_sequence (free_aln (A), -1);
      return buf;
      }
       


/*********************************COPYRIGHT NOTICE**********************************/
/* Centre National de la Recherche Scientifique (CNRS) */
/*and */
/*Cedric Notredame */
/*Tue May 10 12:08:44     2005. */
/*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