binaries/src/tcoffee/t_coffee_source/util_dp_drivers.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "define_header.h"
  10 #include "dp_lib_header.h"
  11
  12 int count_threshold_nodes (Alignment *A, NT_node P, int t);
  13 int set_node_score (Alignment *A, NT_node P, char *mode);
  14 char *split_nodes_nseq (Alignment *A, NT_node P, int max_nseq, char *list);
  15 char *split_nodes_idmax (Alignment *A, NT_node P, int max_id, char *list);
  16 /******************************************************************/
  17 /*                   MAIN DRIVER                                  */
  18 /*                                                                */
  19 /*                                                                */
  20 /******************************************************************/
  21
  22
  23 Constraint_list *profile2list     (Job_TC *job, int nprf)
  24 {
  25   int *seqlist, *cache1, *cache2;
  26
  27   static CLIST_TYPE *entry;
  28   Alignment *A1, *A2, *A;
  29   int a, b, c;
  30   Constraint_list *SCL, *SCL2;
  31   static int *L;
  32   static int max_L_len;
  33   int  max, n_pairs;
  34   int **score, n, s1, s2, si, r1, r2;
  35   char *seqlistb;
  36   int debug=0;
  37   int cons, cons_thres, max_n_pairs, tot_n_pairs;
  38   int *cons_table;
  39   Constraint_list *CL;
  40   char *command;
  41   char *seq;
  42   char *weight;
  43   TC_method *M;
  44   int *iA1, *iA2;
  45   static char *buf1;
  46
  47   if ( !buf1) buf1=vcalloc (1000, sizeof (char));
  48
  49   /*initialize the structure*/
  50   CL=(job->io)->CL;
  51   M=(job->param)->TCM;
  52   command=M->executable;
  53   weight=M->weight;
  54   seq=(job->param)->seq_c;
  55
  56   debug=(getenv("DEBUG_TCOFFEE_profile2list")!=NULL)?1:0;
  57
  58
  59   if ( debug)print_mem_usage (stderr, "IN");
  60   seqlistb=vcalloc (10, sizeof (char));
  61
  62   seqlist=string2num_list (seq)+1;
  63
  64
  65
  66   if (!entry)entry=vcalloc (CL->entry_len, sizeof (int));
  67   entry[SEQ1]=seqlist[1];
  68   entry[SEQ2]=seqlist[2];
  69
  70
  71   A1=seq2profile(CL->S, seqlist[1]);
  72   A2=seq2profile(CL->S, seqlist[2]);
  73
  74   SCL=copy_constraint_list ( CL, SOFT_COPY);
  75
  76   SCL->L=L;SCL->max_L_len=max_L_len;
  77   SCL->fp=NULL;
  78   SCL->M=NULL;
  79   SCL->ne=0;
  80
  81   SCL2=copy_constraint_list ( CL, SOFT_COPY);
  82   SCL2->L=L;SCL2->max_L_len=max_L_len;
  83   SCL2->fp=NULL;
  84   SCL2->M=NULL;
  85   SCL2->ne=0;
  86   /*Merge Sequences*/
  87
  88   SCL->S=merge_seq (A1->S, NULL  );
  89   SCL->S=merge_seq (A2->S, SCL->S);
  90
  91   /*1: Compare the two profiles and identify the N master pairs*/
  92   n=(A1->nseq*A2->nseq);
  93   max=(nprf==0 ||method_uses_structure (M))?n:MIN(nprf,n);
  94
  95   if ( max<n)
  96     {
  97       vsrand(0);
  98       A=align_two_aln (A1, A2, "blosum62mt",-4,-1, "myers_miller_pair_wise");
  99       score=declare_int ( A1->nseq*A2->nseq,3);
 100
 101       for (n=0,a=0; a<A1->nseq; a++)
 102         for ( b=0; b<A2->nseq; b++,n++)
 103           {
 104             score[n][0]=a;
 105             score[n][1]=b;
 106             score[n][2]=(int)get_seq_fsim ( A->seq_al[a], A->seq_al[b+A1->nseq], "-", NOGROUP, NOMATRIX, AVERAGE_POSITIONS);
 107           }
 108       free_aln (A);
 109       sort_int_inv (score,3,2,0, n-1);
 110     }
 111   else
 112     {
 113       score=declare_int (n,3);
 114       for (n=0,a=0; a<A1->nseq; a++)
 115         for ( b=0; b<A2->nseq; b++,n++)
 116           {
 117             score[n][0]=a;
 118             score[n][1]=b;
 119           }
 120     }
 121
 122   iA1=get_name_index (A1->name,A1->nseq, (SCL->S)->name, (SCL->S)->nseq);
 123   iA2=get_name_index (A2->name,A2->nseq, (SCL->S)->name, (SCL->S)->nseq);
 124
 125   /*submit the N pairs*/
 126   for ( n_pairs=1,a=0; a<max; a++)
 127     {
 128       s1=score[a][0];s2=score[a][1];
 129       sprintf ( buf1, "2 %d %d", iA1[s1], iA2[s2]);
 130
 131       job=print_lib_job (job, "param->seq_c=%s",buf1);
 132
 133       /*2: Compute the pairewise library*/
 134       (job->io)->CL=SCL;
 135
 136       SCL=seq2list (job);
 137       /*Unwind the pointer counter:*/job->np--;
 138
 139       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]);
 140
 141
 142       /*3: Update the main library with the pairwise library*/
 143       cache1=seq2inv_pos (A1->seq_al[s1]);
 144       cache2=seq2inv_pos (A2->seq_al[s2]);
 145
 146
 147       if (debug)fprintf ( stderr, " =>%d pairs", SCL->ne);
 148       n_pairs+=(SCL->ne>0)?1:0;
 149       for (c=0; c< SCL->ne; c++)
 150             {
 151               si=vread_clist (SCL, c, SEQ1);
 152               r1=vread_clist(SCL, c, (si==iA1[s1])?R1:R2);
 153               r2=vread_clist(SCL, c, (si==iA1[s1])?R2:R1);
 154
 155               entry[R1]=cache1[r1];
 156               entry[R2]=cache2[r2];
 157
 158               entry[WE]=vread_clist(SCL,c,WE);
 159               entry[CONS]=1;
 160               add_entry2list(entry, SCL2);
 161             }
 162
 163       SCL->ne=0;
 164       vfree ( cache1);vfree ( cache2);
 165       compact_list (SCL2, 0, SCL2->ne, "default");
 166       if (debug)fprintf ( stderr, " =>%d pairs", SCL2->ne);
 167
 168     }
 169
 170   free_sequence (SCL->S,-1);
 171   vfree (iA1); vfree (iA2);
 172
 173   if (debug)fprintf ( stderr, "\nNPairs=%d", n_pairs);
 174
 175
 176   compact_list (SCL2, 0, SCL2->ne, "default");
 177   /*get the concistency distribution*/
 178   cons_table=vcalloc ( 101, sizeof (int));
 179   for (c=0; c< SCL2->ne; c++)
 180     {
 181       entry[R1]=vread_clist(SCL2,c,R1);
 182       entry[R2]=vread_clist(SCL2,c,R2);
 183       entry[WE]=vread_clist(SCL2,c,WE);
 184       entry[CONS]=vread_clist(SCL2,c,CONS);
 185       cons=(entry[CONS]*100)/n_pairs;
 186       cons_table[cons]++;
 187     }
 188
 189   /*Identify the threshold*/
 190   max_n_pairs=(int)((float)(MIN(A1->len_aln, A2->len_aln)*2));
 191
 192   for (cons_thres=0,tot_n_pairs=0,c=100; c>=0; c--)
 193     {
 194
 195       tot_n_pairs+=cons_table[c];
 196       if ( tot_n_pairs>=max_n_pairs){cons_thres=c;c=-1;}
 197     }
 198   vfree (cons_table);
 199
 200   /*Produce the library*/
 201   for (c=0; c< SCL2->ne; c++)
 202     {
 203       entry[R1]=vread_clist(SCL2,c,R1);
 204       entry[R2]=vread_clist(SCL2,c,R2);
 205       entry[WE]=vread_clist(SCL2,c,WE);
 206       entry[CONS]=vread_clist(SCL2,c,CONS);
 207       entry[WE]/=entry[CONS];
 208       cons=(entry[CONS]*100)/n_pairs;
 209
 210
 211       if (cons>=cons_thres)add_entry2list(entry, CL);
 212     }
 213
 214   if ( !seq2R_template_profile (CL->S,seqlist[1]))free_aln (A1);
 215   if ( !seq2R_template_profile (CL->S,seqlist[2]))free_aln (A2);
 216
 217   vfree(seqlistb);
 218   vfree(seqlist-1)  ;
 219   vfree (SCL->L);free_constraint_list ( SCL);
 220
 221
 222   if (debug)print_mem_usage (stderr, "BEF FREE SCL2");
 223   vfree(SCL2->L);free_constraint_list ( SCL2);
 224
 225   if (debug)
 226     {
 227
 228       fprintf ( stderr, "\nCL->ne=%d", CL->ne);
 229       print_mem_usage (stderr, "OUT");
 230     }
 231
 232   (job->io)->CL=CL;
 233   return CL;
 234 }
 235
 236 int method_uses_structure(TC_method *M)
 237 {
 238   if ( strchr (M->seq_type, 'P'))return 1;
 239   else return 0;
 240 }
 241 int method_uses_profile(TC_method *M)
 242 {
 243   if ( strchr (M->seq_type, 'R'))return 1;
 244   else return 0;
 245 }
 246
 247
 248 Constraint_list *seq2list     ( Job_TC *job)
 249     {
 250       char *mode;
 251       Alignment *A=NULL;
 252       Constraint_list *PW_CL;
 253       Constraint_list *RCL=NULL;
 254
 255       int full_prf, nprf;
 256       int *seqlist;
 257
 258
 259
 260
 261       static char *s1, *s2;
 262       Sequence *S, *STL;
 263       Constraint_list *CL;
 264
 265
 266       char *seq;
 267       char *weight;
 268       TC_method *M;
 269
 270
 271       M=(job->param)->TCM;
 272       mode=M->executable;
 273       weight=M->weight;
 274
 275       PW_CL=M->PW_CL;
 276
 277       CL=(job->io)->CL;
 278       seq=(job->param)->seq_c;
 279
 280       S=(CL)?CL->S:NULL;
 281       STL=(CL)?CL->STRUC_LIST:NULL;
 282
 283
 284
 285
 286       seqlist=string2num_list (seq)+1;
 287
 288       if (!s1)s1=vcalloc ( MAXNAMES+1, sizeof (char));
 289       if (!s2)s2=vcalloc ( MAXNAMES+1, sizeof (char));
 290
 291
 292       sprintf (s1, "%s", (CL->S)->name[seqlist[1]]);
 293       sprintf (s2, "%s", (CL->S)->name[seqlist[2]]);
 294
 295 /*Proteins*/
 296
 297
 298       if ( strncmp (CL->profile_comparison, "full", 4)==0)
 299              {
 300                full_prf=1;
 301                if ( CL->profile_comparison[4])nprf=atoi ( CL->profile_comparison+4);
 302                else
 303                  nprf=0;
 304              }
 305       else
 306              {
 307                full_prf=0;
 308              }
 309
 310       if ((method_uses_structure (M)) && profile2P_template_file (CL->S, seqlist[1]) && profile2P_template_file (CL->S, seqlist[2]))
 311         {
 312           RCL=profile2list     (job, nprf);
 313         }
 314       else if ( strm (mode, "ktup_msa"))
 315         {
 316           RCL=hasch2constraint_list (CL->S, CL);
 317         }
 318       else if (     strm (mode, "test_pair") || strm ( mode,"fast_pair")         || strm (mode, "ifast_pair") \
 319                    || strm ( mode, "diag_fast_pair")|| strm (mode, "idiag_fast_pair")\
 320                    || strm ( mode, "blast_pair")    || strm (mode, "lalign_blast_pair") \
 321                    || strm ( mode, "viterbi_pair")  || strm (mode, "slow_pair")      || strm(mode, "glocal_pair") || strm (mode, "biphasic_pair") \
 322                    || strm ( mode, "islow_pair")    || strm (mode, "tm_slow_pair") || strm (mode, "r_slow_pair") \
 323                    || strm ( mode, "lalign_id_pair")|| strm (mode, "tm_lalign_id_pair") || strm (mode , "lalign_len_pair") \
 324                    || strm (mode, "prrp_aln")       || strm ( mode, "test_pair") \
 325                    || strm (mode, "cdna_fast_pair") || strm (mode, "diaa_slow_pair") || strm (mode, "monoaa_slow_pair")\
 326                    || strncmp (mode,"cdna_fast_pair",14)==0     \
 327                    )
 328         {
 329
 330           A=fast_pair (job);
 331           RCL=aln2constraint_list ((A->A)?A->A:A, CL,weight);
 332         }
 333
 334       else if ( strm ( mode, "subop1_pair") || strm ( mode, "subop2_pair") )
 335         {
 336           A=fast_pair (job);
 337           RCL=A->CL;
 338
 339         }
 340       else if ( strm ( mode, "proba_pair") )
 341         {
 342           A=fast_pair (job);
 343           RCL=A->CL;
 344         }
 345       else if ( strm ( mode, "best_pair4prot"))
 346         {
 347           RCL=best_pair4prot (job);
 348         }
 349         else if ( strm ( mode, "best_pair4rna"))
 350         {
 351                 RCL=best_pair4rna (job);
 352         }
 353       else if ( strm ( mode, "exon2_pair"))
 354         {
 355
 356           char weight2[1000];
 357
 358           A=fast_pair (job);
 359           sprintf ( weight2, "%s_subset_objOBJ-",weight);
 360           RCL=aln2constraint_list (A, CL,weight2);
 361         }
 362       else if ( strm ( mode, "exon_pair"))
 363         {
 364
 365           A=fast_pair (job);
 366           RCL=aln2constraint_list (A, CL,weight);
 367
 368         }
 369       else if ( strm ( mode, "exon3_pair"))
 370         {
 371
 372           char weight2[1000];
 373
 374           A=fast_pair (job);
 375           sprintf ( weight2, "%s_subset_objOBJ-",weight);
 376           RCL=aln2constraint_list (A, CL,weight2);
 377         }
 378
 379 /*STRUCTURAL METHODS*/
 380
 381        else if ( strm (mode, "seq_msa"))
 382         {
 383           RCL=seq_msa(M, seq, CL);
 384         }
 385 /*STRUCTURAL METHODS*/
 386       else if (strm (mode, "profile_pair") || strm (mode, "hh_pair"))
 387         {
 388           RCL=profile_pair (M, seq, CL);
 389         }
 390
 391       else if ( strm (mode, "sap_pair"))
 392         {
 393           RCL=sap_pair (seq, weight, CL);
 394         }
 395       else if ( strm (mode, "thread_pair"))
 396         {
 397           RCL=thread_pair (M,seq, CL);
 398         }
 399       else if ( strm (mode, "pdb_pair"))
 400         {
 401           RCL=pdb_pair (M,seq, CL);
 402         }
 403         else if (strm (mode, "rna_pair"))
 404         {
 405                 RCL=rna_pair(M, seq, CL);
 406         }
 407       else if ( strm (mode, "pdbid_pair"))
 408         {
 409           RCL=pdbid_pair (M,seq, CL);
 410         }
 411       else if ( strm (mode, "fugue_pair"))
 412         {
 413           RCL=thread_pair (M,seq, CL);
 414         }
 415       else if ( strm (mode, "lsqman_pair"))
 416         {
 417           RCL=lsqman_pair(seq, CL);
 418         }
 419       else if ( strm ( mode, "align_pdb_pair"))
 420         {
 421           RCL=align_pdb_pair ( seq,"gotoh_pair_wise", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
 422         }
 423       else if ( strm ( mode, "lalign_pdb_pair"))
 424         {
 425           RCL=align_pdb_pair ( seq,"sim_pair_wise_lalign", CL->align_pdb_hasch_mode,CL->align_pdb_param_file,CL, job);
 426         }
 427       else if ( strm ( mode, "align_pdb_pair_2"))
 428         {
 429           RCL=align_pdb_pair_2 ( seq, CL);
 430         }
 431       else
 432         {
 433           fprintf ( CL->local_stderr, "\nERROR: THE FUNCTION %s DOES NOT EXIST [FATAL:%s]\n", mode, PROGRAM);crash("");
 434         }
 435       add_method_output2method_log (NULL,NULL, (A&&A->len_aln)?A:NULL,RCL, NULL);
 436       RCL=(RCL==NULL)?CL:RCL;
 437
 438
 439       vfree ( seqlist-1);
 440       free_aln (A);
 441       return RCL;
 442     }
 443
 444 Constraint_list *method2pw_cl (TC_method *M, Constraint_list *CL)
 445     {
 446       char *mode;
 447       Constraint_list *PW_CL=NULL;
 448       Sequence *S;
 449       char mat[100], *m;
 450       char group_mat[100];
 451
 452
 453
 454       mode=M->executable;
 455       PW_CL=copy_constraint_list ( CL, SOFT_COPY);
 456       PW_CL->pw_parameters_set=1;
 457
 458
 459
 460       S=(PW_CL)?PW_CL->S:NULL;
 461
 462       /*DNA or Protein*/
 463       m=PW_CL->method_matrix;
 464       if (  strm ((PW_CL->S)->type, "PROTEIN"))
 465            {
 466
 467              sprintf ( mat, "%s", (strm(m, "default"))?"blosum62mt":m);
 468              sprintf (group_mat, "vasiliky");
 469              PW_CL->ktup=2;
 470
 471            }
 472       else if (  strm ((PW_CL->S)->type, "DNA") || strm ((PW_CL->S)->type, "RNA") )
 473            {
 474
 475               sprintf(group_mat, "idmat");
 476               sprintf ( mat, "%s", (strm(m, "default"))?"dna_idmat":m);
 477               PW_CL->ktup=5;
 478
 479            }
 480       if ( M->matrix[0])sprintf ( mat, "%s", M->matrix);
 481
 482       PW_CL->M=read_matrice (mat);
 483
 484
 485       if ( M->gop!=UNDEFINED) {PW_CL->gop=M->gop;}
 486       else
 487         {
 488           PW_CL->gop= get_avg_matrix_mm (PW_CL->M, AA_ALPHABET)*10;
 489         }
 490
 491       if ( M->gep!=UNDEFINED)PW_CL->gep=M->gep;
 492       else PW_CL->gep=-1;
 493
 494
 495
 496       if ( strm2 ( mode,"fast_pair", "ifast_pair"))
 497             {
 498               PW_CL->maximise=1;
 499               PW_CL->TG_MODE=1;
 500               PW_CL->use_fragments=0;
 501               if ( !PW_CL->use_fragments)PW_CL->diagonal_threshold=0;
 502               else PW_CL->diagonal_threshold=6;
 503
 504               sprintf (PW_CL->dp_mode, "fasta_pair_wise");
 505               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 506
 507               if ( strm ( mode, "fast_pair"))
 508                 {
 509                   PW_CL->L=NULL;
 510
 511                   PW_CL->get_dp_cost=slow_get_dp_cost;
 512                   PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 513                   PW_CL->extend_jit=0;
 514                 }
 515             }
 516         else if ( strm2 ( mode,"diag_fast_pair","idiag_fast_pair"))
 517             {
 518               PW_CL->L=NULL;
 519               PW_CL->maximise=1;
 520               PW_CL->TG_MODE=1;
 521               PW_CL->S=CL->S;
 522
 523               PW_CL->use_fragments=1;
 524               PW_CL->diagonal_threshold=3;
 525
 526               sprintf (PW_CL->dp_mode, "fasta_pair_wise");
 527               PW_CL->ktup=1;
 528               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 529
 530               PW_CL->get_dp_cost=slow_get_dp_cost;
 531               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 532
 533               PW_CL->extend_jit=0;
 534             }
 535         else if ( strm ( mode,"blast_pair"))
 536             {
 537               PW_CL->L=NULL;
 538               PW_CL->maximise=1;
 539               PW_CL->TG_MODE=1;
 540
 541               PW_CL->use_fragments=0;
 542
 543               PW_CL->pair_wise=gotoh_pair_wise;
 544               PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
 545               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 546               PW_CL->extend_jit=0;
 547             }
 548         else if ( strm ( mode,"lalign_blast_pair"))
 549             {
 550               PW_CL->L=NULL;
 551               PW_CL->maximise=1;
 552               PW_CL->TG_MODE=1;
 553
 554               PW_CL->use_fragments=0;
 555               PW_CL->pair_wise=sim_pair_wise_lalign;
 556               PW_CL->evaluate_residue_pair=evaluate_blast_profile_score;
 557               PW_CL->lalign_n_top=10;
 558
 559               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 560               PW_CL->extend_jit=0;
 561             }
 562         else if ( strm ( mode,"viterbi_pair"))
 563             {
 564               PW_CL->maximise=1;
 565               PW_CL->TG_MODE=1;
 566               PW_CL->use_fragments=0;
 567               sprintf (PW_CL->dp_mode, "viterbi_pair_wise");
 568               PW_CL->L=NULL;
 569               PW_CL->get_dp_cost=slow_get_dp_cost;
 570               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 571               PW_CL->extend_jit=0;
 572             }
 573         else if ( strm ( mode,"glocal_pair"))
 574             {
 575               PW_CL->maximise=1;
 576               PW_CL->TG_MODE=1;
 577               PW_CL->use_fragments=0;
 578               sprintf (PW_CL->dp_mode, "glocal_pair_wise");
 579               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 580
 581               PW_CL->L=NULL;
 582               PW_CL->get_dp_cost=slow_get_dp_cost;
 583               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 584               PW_CL->extend_jit=0;
 585             }
 586         else if ( strm ( mode,"test_pair"))
 587             {
 588               PW_CL->maximise=1;
 589               PW_CL->TG_MODE=1;
 590               PW_CL->use_fragments=0;
 591               sprintf (PW_CL->dp_mode, "test_pair_wise");
 592             }
 593         else if ( strm ( mode,"sticky_pair"))
 594           {
 595             PW_CL->maximise=1;
 596             PW_CL->TG_MODE=1;
 597             PW_CL->use_fragments=0;
 598             PW_CL->L=NULL;
 599             PW_CL->get_dp_cost=cw_profile_get_dp_cost;
 600             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 601             PW_CL->extend_jit=0;
 602             sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp_sticky");
 603             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 604           }
 605
 606         else if ( strm ( mode,"slow_pair")|| strm (mode, "islow_pair" ) )
 607             {
 608
 609               PW_CL->maximise=1;
 610               PW_CL->TG_MODE=1;
 611               PW_CL->use_fragments=0;
 612               sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
 613               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 614
 615               if ( strm ( "islow_pair", mode))
 616                 {
 617                   PW_CL->get_dp_cost=slow_get_dp_cost;
 618                   PW_CL->evaluate_residue_pair=residue_pair_extended_list;
 619                   PW_CL->extend_jit=1;
 620                 }
 621               else if ( strm ("slow_pair", mode) )
 622                 {
 623                   PW_CL->L=NULL;
 624                   PW_CL->get_dp_cost=cw_profile_get_dp_cost;
 625                   PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 626                   PW_CL->extend_jit=0;
 627                 }
 628             }
 629         else if ( strm (mode, "subop1_pair"))
 630             {
 631
 632               PW_CL->maximise=1;
 633               PW_CL->TG_MODE=1;
 634               PW_CL->use_fragments=0;
 635               sprintf (PW_CL->dp_mode, "subop1_pair_wise");
 636               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 637               PW_CL->L=NULL;
 638               PW_CL->get_dp_cost=slow_get_dp_cost;
 639               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 640               PW_CL->extend_jit=0;
 641             }
 642        else if ( strm (mode, "biphasic_pair"))
 643             {
 644               PW_CL->maximise=1;
 645               PW_CL->TG_MODE=1;
 646               PW_CL->use_fragments=0;
 647               sprintf (PW_CL->dp_mode, "biphasic_pair_wise");
 648               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 649               PW_CL->L=NULL;
 650               PW_CL->get_dp_cost=cw_profile_get_dp_cost;
 651               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 652               PW_CL->extend_jit=0;
 653             }
 654       else if ( strm (mode, "proba_pair"))
 655             {
 656
 657               PW_CL->maximise=1;
 658               PW_CL->TG_MODE=1;
 659               PW_CL->use_fragments=0;
 660               sprintf (PW_CL->dp_mode, "proba_pair_wise");
 661               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 662               PW_CL->L=NULL;
 663               PW_CL->get_dp_cost=slow_get_dp_cost;
 664               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 665               PW_CL->extend_jit=0;
 666             }
 667
 668         else if ( strm (mode, "diaa_slow_pair"))
 669             {
 670
 671               PW_CL->maximise=1;
 672               PW_CL->TG_MODE=1;
 673               PW_CL->use_fragments=0;
 674               sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
 675               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 676               PW_CL->L=NULL;
 677               PW_CL->get_dp_cost=slow_get_dp_cost;
 678               PW_CL->evaluate_residue_pair=evaluate_diaa_matrix_score;
 679               PW_CL->extend_jit=0;
 680             }
 681         else if ( strm (mode, "r_slow_pair"))
 682             {
 683               PW_CL->maximise=1;
 684               PW_CL->TG_MODE=1;
 685               PW_CL->use_fragments=0;
 686               sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
 687               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 688               PW_CL->L=NULL;
 689               PW_CL->get_dp_cost=slow_get_dp_cost;
 690               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 691               PW_CL->extend_jit=0;
 692               PW_CL->reverse_seq=1;
 693             }
 694         else if ( strm (mode, "tm_slow_pair"))
 695             {
 696               PW_CL->maximise=1;
 697               PW_CL->TG_MODE=1;
 698               PW_CL->use_fragments=0;
 699               sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
 700               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 701               PW_CL->L=NULL;
 702               PW_CL->get_dp_cost=slow_get_dp_cost;
 703               PW_CL->evaluate_residue_pair=evaluate_tm_matrix_score;
 704               PW_CL->extend_jit=0;
 705             }
 706         else if ( strm (mode, "monoaa_slow_pair"))
 707             {
 708
 709               PW_CL->maximise=1;
 710               PW_CL->TG_MODE=1;
 711               PW_CL->use_fragments=0;
 712               sprintf (PW_CL->dp_mode, "gotoh_pair_wise_lgp");
 713               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 714               PW_CL->L=NULL;
 715               PW_CL->get_dp_cost=slow_get_dp_cost;
 716               PW_CL->evaluate_residue_pair=evaluate_monoaa_matrix_score;
 717               PW_CL->extend_jit=0;
 718             }
 719         else if ( strm (mode, "subop2_pair"))
 720           {
 721
 722             PW_CL->maximise=1;
 723               PW_CL->TG_MODE=1;
 724               PW_CL->use_fragments=0;
 725               sprintf (PW_CL->dp_mode, "subop2_pair_wise");
 726               sprintf (PW_CL->matrix_for_aa_group, group_mat);
 727               PW_CL->L=NULL;
 728               PW_CL->get_dp_cost=slow_get_dp_cost;
 729               PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 730               PW_CL->extend_jit=0;
 731             }
 732
 733         else if (strm ( mode, "exon2_pair"))
 734           {
 735             int a;
 736             PW_CL->maximise=1;
 737             PW_CL->TG_MODE=1;
 738
 739             PW_CL->use_fragments=0;
 740             sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
 741             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 742             PW_CL->L=NULL;
 743
 744             for ( a=0; a<60; a++)
 745               {
 746                 PW_CL->M['x'-'A'][a]=0;
 747                 PW_CL->M[a]['x'-'A']=0;
 748                 PW_CL->M['X'-'A'][a]=0;
 749                 PW_CL->M[a]['X'-'A']=0;
 750               }
 751             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 752             PW_CL->extend_jit=0;
 753           }
 754         else if (strm ( mode, "exon3_pair"))
 755           {
 756             int a;
 757             PW_CL->maximise=1;
 758             PW_CL->TG_MODE=1;
 759
 760             PW_CL->use_fragments=0;
 761             sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
 762             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 763             PW_CL->L=NULL;
 764
 765             for ( a=0; a<60; a++)
 766               {
 767                 PW_CL->M['x'-'A'][a]=0;
 768                 PW_CL->M[a]['x'-'A']=0;
 769                 PW_CL->M['X'-'A'][a]=0;
 770                 PW_CL->M[a]['X'-'A']=0;
 771               }
 772             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 773             PW_CL->extend_jit=0;
 774           }
 775         else if (strm ( mode, "exon_pair"))
 776           {
 777             int a;
 778             PW_CL->maximise=1;
 779             PW_CL->TG_MODE=1;
 780
 781             PW_CL->use_fragments=0;
 782             sprintf (PW_CL->dp_mode, "myers_miller_pair_wise");
 783             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 784             PW_CL->L=NULL;
 785
 786             for ( a=0; a<60; a++)
 787               {
 788                 PW_CL->M['x'-'A'][a]=0;
 789                 PW_CL->M[a]['x'-'A']=0;
 790                 PW_CL->M['X'-'A'][a]=0;
 791                 PW_CL->M[a]['X'-'A']=0;
 792               }
 793             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 794             PW_CL->extend_jit=0;
 795           }
 796         else if ( strm ( mode , "lalign_len_pair"))
 797           {
 798             PW_CL->L=NULL;
 799             PW_CL->maximise=1;
 800             PW_CL->TG_MODE=1;
 801             PW_CL->use_fragments=0;
 802             PW_CL->pair_wise=sim_pair_wise_lalign;
 803             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 804             PW_CL->get_dp_cost=slow_get_dp_cost;
 805             PW_CL->lalign_n_top=CL->lalign_n_top;
 806             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 807             PW_CL->extend_jit=0;
 808           }
 809         else if ( strm ( mode , "lalign_id_pair"))
 810           {
 811             PW_CL->L=NULL;
 812             PW_CL->maximise=1;
 813             PW_CL->TG_MODE=1;
 814             PW_CL->use_fragments=0;
 815             PW_CL->pair_wise=sim_pair_wise_lalign;
 816             PW_CL->evaluate_residue_pair=evaluate_matrix_score;
 817             PW_CL->get_dp_cost=slow_get_dp_cost;
 818             PW_CL->lalign_n_top=CL->lalign_n_top;
 819             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 820             PW_CL->extend_jit=0;
 821           }
 822         else if ( strm ( mode , "tm_lalign_id_pair"))
 823           {
 824             PW_CL->L=NULL;
 825             PW_CL->maximise=1;
 826             PW_CL->TG_MODE=1;
 827             PW_CL->use_fragments=0;
 828             PW_CL->pair_wise=sim_pair_wise_lalign;
 829             PW_CL->evaluate_residue_pair=evaluate_tm_matrix_score;
 830             PW_CL->get_dp_cost=slow_get_dp_cost;
 831             PW_CL->lalign_n_top=CL->lalign_n_top;
 832             sprintf (PW_CL->matrix_for_aa_group, group_mat);
 833             PW_CL->extend_jit=0;
 834           }
 835 /*CDNA*/
 836         else if ( strm ( mode, "cdna_cfast_pair"))
 837             {
 838               PW_CL->L=NULL;
 839               PW_CL->maximise=1;
 840               PW_CL->TG_MODE=1;
 841               PW_CL->S=CL->S;
 842               PW_CL->use_fragments=0;
 843               sprintf (PW_CL->dp_mode, "cfasta_cdna_pair_wise");
 844
 845               PW_CL->M=read_matrice (strcpy ( mat, "blosum62mt"));
 846               PW_CL->extend_jit=0;
 847
 848               PW_CL->f_gop=CL->f_gop;
 849               PW_CL->f_gep=CL->f_gep;
 850               PW_CL->get_dp_cost=get_dp_cost;
 851               PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
 852               PW_CL->ktup=1;
 853             }
 854         else if ( strm ( mode, "cdna_fast_pair") ||  strncmp (mode,"cdna_fast_pair",14)==0)
 855             {
 856
 857               PW_CL->L=NULL;
 858               PW_CL->maximise=1;
 859               PW_CL->TG_MODE=1;
 860               PW_CL->use_fragments=0;
 861               sprintf (PW_CL->dp_mode, "fasta_cdna_pair_wise");
 862
 863               PW_CL->extend_jit=0;
 864               PW_CL->gop=-5;
 865               PW_CL->gep=-1;
 866               PW_CL->f_gop=-15;
 867               PW_CL->f_gep=0;
 868
 869               PW_CL->get_dp_cost=get_dp_cost;
 870               PW_CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
 871               PW_CL->ktup=1;
 872             }
 873         else
 874           {
 875             free_constraint_list (PW_CL);
 876             PW_CL=NULL;
 877           }
 878
 879
 880       if (!strm (CL->method_evaluate_mode, "default"))
 881         {
 882           choose_extension_mode (CL->method_evaluate_mode, PW_CL);
 883         }
 884       return PW_CL;
 885     }
 886 /******************************************************************/
 887 /*                   MULTIPLE ALIGNMENTS                          */
 888 /*                                                                */
 889 /*                                                                */
 890 /******************************************************************/
 891 Alignment * compute_prrp_aln (Alignment *A, Constraint_list *CL)
 892    {
 893        char *tmpseq=NULL;
 894        char *tmpaln=NULL;
 895        char command[10000];
 896        Sequence *S;
 897
 898        tmpseq=vtmpnam(NULL);
 899        tmpaln=vtmpnam(NULL);
 900
 901
 902        A=seq2aln ( CL->S, A, 1);
 903        output_gotoh_seq (tmpseq, A);
 904        sprintf ( command, "prrp -E/dev/null -o%s -F9 %s >/dev/null", tmpaln, tmpseq);
 905        my_system (command);
 906        if (!check_file_exists(tmpaln)){return NULL;}
 907        S=get_fasta_sequence (tmpaln, NULL);
 908
 909        S->contains_gap=0;
 910        A=seq2aln(S, A,0);
 911        free_sequence (S, S->nseq);
 912
 913        return A;
 914    }
 915
 916 Alignment *seq2clustalw_aln  (Sequence *S)
 917 {
 918   return aln2clustalw_aln (seq2aln (S, NULL,RM_GAP), NULL);
 919 }
 920
 921 Alignment * aln2clustalw_aln (Alignment *B, Constraint_list *CL)
 922 {
 923   char *seq=NULL,*aln=NULL, command[1000];
 924
 925   output_fasta_seq (seq=vtmpnam (NULL), B);
 926   sprintf ( command, "clustalw -infile=%s -outorder=input -outfile=%s %s", seq, aln=vtmpnam (NULL), TO_NULL_DEVICE);
 927   my_system (command);
 928
 929   if (!check_file_exists(aln))
 930     return NULL;
 931   else{B->nseq=0;return main_read_aln(aln,B);}
 932 }
 933 Alignment * compute_tcoffee_aln_quick (Alignment *A, Constraint_list *CL)
 934    {
 935        char *tmpseq=NULL;
 936        char *tmpaln=NULL;
 937        char command[10000];
 938
 939
 940        tmpseq=vtmpnam(NULL);
 941        tmpaln=vtmpnam(NULL);
 942
 943        if ( CL)A=seq2aln ( CL->S, A, 1);
 944        output_fasta_seq (tmpseq, A);
 945
 946        sprintf ( command, "t_coffee  -seq %s -very_fast -outfile %s -quiet ",tmpseq,tmpaln);
 947
 948        my_system (command);
 949        if (!check_file_exists(tmpaln))return NULL;
 950        A->nseq=0;
 951        A=main_read_aln(tmpaln,A);
 952
 953        vremove( tmpseq);
 954        vremove (tmpaln);
 955        return A;
 956    }
 957
 958 Alignment * compute_clustalw_aln (Alignment *A, Constraint_list *CL)
 959    {
 960        char *tmpseq=NULL;
 961        char *tmpaln=NULL;
 962        char command[10000];
 963
 964
 965        tmpseq=vtmpnam(NULL);
 966        tmpaln=vtmpnam(NULL);
 967
 968        A=seq2aln ( CL->S, A, 1);
 969        output_fasta_seq (tmpseq, A);
 970
 971        sprintf ( command, "clustalw %sinfile=%s %soutfile=%s %s",CWF, tmpseq,CWF, tmpaln,TO_NULL_DEVICE);
 972
 973        my_system (command);
 974        if (!check_file_exists(tmpaln))return NULL;
 975        A->nseq=0;
 976        A=main_read_aln(tmpaln,A);
 977
 978        vremove( tmpseq);
 979        vremove (tmpaln);
 980        return A;
 981    }
 982
 983 Alignment * realign_block ( Alignment *A, int col1, int col2, char *pg)
 984 {
 985   /*Uses pg: (pg -infile=<infile> -outfile=<outfile> to realign the block [col1 col2[
 986     Only guaranteed if pg can handle empty sequences
 987     set pg to NULL to use the default program
 988   */
 989
 990
 991   Alignment *L, *M, *R;
 992   char *seq_name;
 993   char *aln_name;
 994   char command[1000], script[1000];
 995
 996
 997
 998   seq_name=vtmpnam(NULL);
 999   aln_name=vtmpnam (NULL);
1000
1001   L=copy_aln (A, NULL);
1002   M=copy_aln (A, NULL);
1003   R=copy_aln (A, NULL);
1004
1005   L=extract_aln ( L, 0, col1);
1006   M=extract_aln ( M, col1, col2);
1007   R=extract_aln ( R, col2, A->len_aln);
1008   output_fasta_seq (seq_name, M);
1009
1010   sprintf ( script, "%s", (pg==NULL)?"t_coffee":pg);
1011
1012   sprintf ( command, "%s -infile=%s -outfile=%s %s", script,seq_name, aln_name, TO_NULL_DEVICE);
1013   my_system ( command);
1014   free_aln (M);
1015   M=main_read_aln (aln_name, NULL);
1016
1017
1018   M=reorder_aln (M, L->name,L->nseq);
1019   L=aln_cat (L, M);
1020   L=aln_cat (L, R);
1021   A=copy_aln (L, A);
1022   free_aln (L);free_aln (M); free_aln (R);
1023
1024   return A;
1025 }
1026
1027
1028
1029
1030
1031 /******************************************************************/
1032 /*                  DNA                                           */
1033 /*                                                                */
1034 /*                                                                */
1035 /******************************************************************/
1036
1037
1038 /******************************************************************/
1039 /*                   STRUCTURES                                   */
1040 /*                                                                */
1041 /*                                                                */
1042 /******************************************************************/
1043
1044
1045
1046
1047
1048 Constraint_list * align_pdb_pair_2 (char *seq, Constraint_list *CL)
1049         {
1050             char *tmp_name=NULL;
1051             int s1, s2;
1052
1053
1054             static char *command;
1055             static char *program;
1056
1057             tmp_name=vtmpnam ( NULL);
1058
1059             if ( !program)program=vcalloc ( LONG_STRING, sizeof (char));
1060             if ( !command)command=vcalloc ( LONG_STRING, sizeof (char));
1061
1062 #ifndef     ALIGN_PDB_4_TCOFFEE
1063             if ( getenv ( "ALIGN_PDB_4_TCOFFEE")==NULL)crash ("ALIGN_PDB_4_TCOFFEE IS NOT DEFINED");
1064             else  sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
1065 #else
1066             if ( getenv ( "ALIGN_4_TCOFFEE")==NULL)sprintf (program, "%s", ALIGN_PDB_4_TCOFFEE);
1067             else sprintf ( program, "%s", (getenv ( "ALIGN_PDB_4_TCOFFEE")));
1068 #endif
1069
1070             atoi(strtok (seq,SEPARATORS));
1071             s1=atoi(strtok (NULL,SEPARATORS));
1072             s2=atoi(strtok (NULL,SEPARATORS));
1073
1074             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);
1075
1076             my_system  ( command);
1077             CL=read_constraint_list_file(CL, tmp_name);
1078
1079
1080             vremove ( tmp_name);
1081
1082
1083             return CL;
1084         }
1085
1086 Constraint_list *align_pdb_pair   (char *seq_in, char *dp_mode,char *evaluate_mode, char *file, Constraint_list *CL, Job_TC *job)
1087  {
1088             int s1, s2;
1089             char seq[1000];
1090             char name1[1000];
1091             char name2[1000];
1092
1093
1094             Constraint_list *PWCL;
1095             Alignment *F;
1096
1097             sprintf ( seq, "%s",seq_in);
1098             atoi(strtok (seq,SEPARATORS));
1099             s1=atoi(strtok (NULL,SEPARATORS));
1100             s2=atoi(strtok (NULL,SEPARATORS));
1101
1102
1103
1104             sprintf (name1, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s1], (CL->S)->name[s2], dp_mode);
1105             sprintf (name2, "%s%s_%s.%s.align_pdb", get_cache_dir(),(CL->S)->name[s2], (CL->S)->name[s1], dp_mode);
1106
1107
1108             if ( check_file_exists (name1) &&   is_lib(name1))CL=read_constraint_list_file(CL,name1);
1109             else if ( check_file_exists (name2) &&   is_lib(name2))CL=read_constraint_list_file(CL,name2);
1110             else
1111               {
1112                 PWCL=set_constraint_list4align_pdb ( CL,s1,dp_mode, evaluate_mode, NULL);
1113                 PWCL=set_constraint_list4align_pdb ( CL,s2,dp_mode, evaluate_mode, NULL);
1114                 ((job->param)->TCM)->PW_CL=PWCL;
1115                 F=fast_pair (job);
1116                 output_constraints (name1, "100", F);
1117                 CL=aln2constraint_list (F, CL, "100");
1118                 free_aln (F);
1119               }
1120             return CL;
1121         }
1122
1123 Constraint_list * profile_pair (TC_method *M , char *in_seq, Constraint_list *CL)
1124         {
1125
1126           char seq[1000];
1127           int a, s1, s2;
1128           char *result,*prf1_file, *prf2_file;
1129           Alignment *F=NULL, *A1, *A2;
1130           FILE *fp;
1131           char command[10000];
1132           char *param;
1133
1134           if ( M->executable2[0]=='\0')
1135             fprintf ( stderr, "\nERROR: profile_pair requires a method: thread_pair@EP@executable2@<method> [FATAL:%s]\n", PROGRAM);
1136
1137
1138           sprintf ( seq, "%s", in_seq);
1139           atoi(strtok (seq,SEPARATORS));
1140           s1=atoi(strtok (NULL,SEPARATORS));
1141           s2=atoi(strtok (NULL,SEPARATORS));
1142
1143           A1=seq2R_template_profile(CL->S,s1);
1144           A2=seq2R_template_profile(CL->S,s2);
1145
1146
1147           prf1_file=vtmpnam (NULL);
1148           fp=vfopen (prf1_file, "w");
1149           if ( A1)
1150             {
1151               fprintf (fp, ">%s\n%s\n",(CL->S)->name[s1], aln2cons_seq_mat(A1, "blosum62mt"));
1152               for ( a=0; a< A1->nseq; a++)fprintf (fp, ">prf_seq1_%d\n%s\n", a, A1->seq_al[a]);
1153             }
1154           else
1155             {
1156               fprintf ( fp, ">%s\n%s\n", (CL->S)->name[s1], (CL->S)->seq[s1]);
1157             }
1158           vfclose (fp);
1159
1160           prf2_file=vtmpnam (NULL);
1161           fp=vfopen (prf2_file, "w");
1162           if (A2)
1163             {
1164               fprintf (fp, ">%s\n%s\n",(CL->S)->name[s2], aln2cons_seq_mat(A2, "blosum62mt"));
1165               for ( a=0; a< A2->nseq; a++)fprintf (fp, ">prf_seq2_%d\n%s\n", a, A2->seq_al[a]);
1166             }
1167           else
1168             {
1169               fprintf ( fp, ">%s\n%s\n", (CL->S)->name[s2], (CL->S)->seq[s2]);
1170             }
1171           vfclose (fp);
1172
1173           result=vtmpnam (NULL);
1174           if ( M->param)
1175             {
1176               param=vcalloc(strlen (M->param)+1, sizeof (char));
1177               sprintf ( param, "%s", M->param);
1178               param=substitute ( param, " ", "");
1179               param=substitute ( param, "\n", "");
1180             }
1181
1182           sprintf ( command, "tc_generic_method.pl -mode=profile_pair -method=%s %s%s %s%s %s%s -param=%s -tmpdir=%s", M->executable2,M->in_flag,prf1_file, M->in_flag2,prf2_file,M->out_flag, result, param, get_tmp_4_tcoffee());
1183
1184           my_system ( command);
1185
1186
1187
1188           if ( !check_file_exists (result))
1189             {
1190               fprintf ( stderr, "\n\tprofile_pair/%s failed:\n\t%s\n",M->executable2, command);
1191               myexit (EXIT_FAILURE);
1192             }
1193           else if ( is_lib (result))
1194             {
1195               CL=read_constraint_list_file(CL,result);
1196             }
1197           else if ( is_aln (result))
1198             {
1199               F=main_read_aln (result, NULL);
1200               char *name1, *name2;
1201               name1=(CL->S)->name[s1];
1202               name2=(CL->S)->name[s2];
1203
1204               fp=vfopen (result, "w");
1205               for ( a=0; a< F->nseq; a++)
1206                 if (strm ( F->name[a], name1) || strm (F->name[a], name2))
1207                   fprintf ( fp, ">%s\n%s\n", F->name[a], F->seq_al[a]);
1208               vfclose (fp);
1209               free_aln (F);
1210               F=main_read_aln (result, NULL);
1211               CL=aln2constraint_list (F, CL, "sim");
1212               free_aln (F);
1213             }
1214           return CL;
1215         }
1216 Constraint_list * pdbid_pair (TC_method *M , char *in_seq, Constraint_list *CL)
1217         {
1218
1219           char seq[1000];
1220           char *alternative_method;
1221           char *current_,method;
1222
1223           int s1, s2;
1224           char *result, *pdb1,*pdb1_file, *pdb2, *pdb2_file;
1225           Alignment *F=NULL;
1226           char *command;
1227
1228
1229           if ( M->executable2[0]=='\0')
1230             {
1231             fprintf ( stderr, "\nERROR: pdbid_pair requires a structural alignment method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
1232             myexit (EXIT_FAILURE);
1233             }
1234           sprintf ( seq, "%s", in_seq);
1235
1236
1237           atoi(strtok (seq,SEPARATORS));
1238           s1=atoi(strtok (NULL,SEPARATORS));
1239           s2=atoi(strtok (NULL,SEPARATORS));
1240
1241           pdb1=seq2P_pdb_id(CL->S,s1);
1242           pdb2=seq2P_pdb_id(CL->S,s2);
1243
1244           if (!is_pdb_name (pdb1) || !is_pdb_name(pdb2))
1245             {
1246               return CL;
1247             }
1248
1249
1250           result=vtmpnam (NULL);
1251           command = vcalloc ( 1000, sizeof (char));
1252
1253           sprintf ( command, "tc_generic_method.pl -mode=pdbid_pair -method=%s %s%s %s%s %s%s -email=%s -cache=%s -tmpdir=%s", M->executable2,M->in_flag,pdb1, M->in_flag2,pdb2,M->out_flag, result, Email(ENV,SET),get_cache_dir(), get_tmp_4_tcoffee());
1254           my_system ( command);
1255           vfree (command);
1256           if (file_is_empty (result))return CL;
1257           else
1258             {
1259               F=main_read_aln (result, NULL);
1260
1261               if ( !F)
1262                 {
1263                   fprintf ( stderr, "\n\tpdb_pair/%s failed:\n\t%s\n",M->executable2, command);
1264                 }
1265               else
1266                 {
1267
1268                   sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
1269                   sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
1270                   F=fix_aln_seq (F, CL->S);
1271                   CL=aln2constraint_list (F, CL, "sim");
1272                 }
1273               free_aln (F);
1274             }
1275           return CL;
1276         }
1277
1278 Constraint_list * pdb_pair (TC_method *M , char *in_seq, Constraint_list *CL)
1279         {
1280
1281           char seq[1000];
1282           int s1, s2;
1283           char *result, *pdb1,*pdb1_file, *pdb2, *pdb2_file;
1284           Alignment *F=NULL;
1285
1286
1287           char command[10000];
1288
1289           if ( M->executable2[0]=='\0')
1290             {
1291             fprintf ( stderr, "\nERROR: pdb_pair requires a structural alignment method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
1292             myexit (EXIT_FAILURE);
1293             }
1294
1295           sprintf ( seq, "%s", in_seq);
1296
1297
1298           atoi(strtok (seq,SEPARATORS));
1299           s1=atoi(strtok (NULL,SEPARATORS));
1300           s2=atoi(strtok (NULL,SEPARATORS));
1301
1302           pdb1=seq2P_template_file(CL->S,s1);
1303           pdb2=seq2P_template_file(CL->S,s2);
1304           if ( !pdb1 || !pdb2) return CL;
1305
1306
1307           pdb1_file=vtmpnam (NULL);
1308           pdb2_file=vtmpnam (NULL);
1309
1310
1311
1312           sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST  -nodiagnostic > %s", pdb1, pdb1_file);
1313           my_system (command);
1314
1315
1316
1317           sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST  -nodiagnostic > %s", pdb2, pdb2_file);
1318           my_system (command);
1319
1320
1321           result=vtmpnam (NULL);
1322
1323           sprintf ( command, "tc_generic_method.pl -mode=pdb_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pdb1_file, M->in_flag2,pdb2_file,M->out_flag, result, get_tmp_4_tcoffee());
1324           my_system ( command);
1325
1326           F=main_read_aln (result, NULL);
1327
1328           if ( !F)
1329             {
1330               fprintf ( stderr, "\n\tpdb_pair/%s failed:\n\t%s\n",M->executable2, command);
1331             }
1332           else
1333             {
1334
1335               sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
1336               sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
1337               F=fix_aln_seq (F, CL->S);
1338               CL=aln2constraint_list (F, CL, "sim");
1339             }
1340
1341
1342           free_aln (F);
1343           return CL;
1344         }
1345
1346 Constraint_list * seq_msa (TC_method *M , char *in_seq, Constraint_list *CL)
1347 {
1348   char seq[1000];
1349   char *infile, *outfile;
1350   int a, n, s;
1351   Alignment *F=NULL;
1352   FILE *fp;
1353   char command[1000];
1354
1355
1356   infile=vtmpnam (NULL);
1357   outfile=vtmpnam (NULL);
1358
1359   sprintf ( seq, "%s", in_seq);
1360
1361   n=atoi(strtok (seq,SEPARATORS));
1362
1363   fp=vfopen (infile, "w");
1364   for ( a=0; a<n; a++)
1365     {
1366       s=atoi(strtok (NULL,SEPARATORS));
1367
1368       fprintf (fp, ">%s\n%s\n", (CL->S)->name[s], (CL->S)->seq[s]);
1369     }
1370   vfclose (fp);
1371
1372   sprintf ( command, "t_coffee -other_pg tc_generic_method.pl -mode=seq_msa -method=%s %s%s %s%s -tmpdir=%s", M->executable2, M->in_flag, infile, M->out_flag, outfile, get_tmp_4_tcoffee());
1373   my_system (command);
1374
1375
1376   if ( strm (M->out_mode, "aln") ||  strm (M->out_mode, "A"))
1377     {
1378       F=main_read_aln (outfile, NULL);
1379       if ( !F)
1380         {
1381           fprintf ( stderr, "\n\tseq_msa/%s failed:\n\t%s\n", M->executable2,command);
1382         }
1383       else
1384         {
1385           CL=aln2constraint_list (F, CL, "sim");
1386         }
1387       free_aln (F);
1388     }
1389   else if ( strm (M->out_mode, "fL")|| strm (M->out_mode, "lib"))
1390     {
1391       Constraint_list *NCL;
1392       NCL=read_constraint_list_file(CL,outfile);
1393       if ( !NCL)
1394         {
1395           fprintf ( stderr, "\n\tseq_msa/%s failed:\n\t%s\n", M->executable2,command);
1396         }
1397       else
1398         {
1399           CL=NCL;
1400         }
1401     }
1402   return CL;
1403 }
1404 Constraint_list * thread_pair (TC_method *M , char *in_seq, Constraint_list *CL)
1405         {
1406
1407           char seq[1000];
1408           int s1, s2;
1409
1410           if ( M->executable2[0]=='\0')
1411             {
1412               fprintf ( stderr, "\nERROR: thread_pair requires a threading method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
1413               myexit (EXIT_FAILURE);
1414             }
1415
1416           sprintf ( seq, "%s", in_seq);
1417           atoi(strtok (seq,SEPARATORS));
1418           s1=atoi(strtok (NULL,SEPARATORS));
1419           s2=atoi(strtok (NULL,SEPARATORS));
1420
1421           CL=thread_pair2(M,s1, s2, CL);
1422           CL=thread_pair2(M,s2, s1, CL);
1423
1424           return CL;
1425         }
1426
1427
1428
1429 Constraint_list* thread_pair2 ( TC_method *M, int s1, int s2, Constraint_list *CL)
1430   {
1431     char *result, *pep, *pdb, *pdb1;
1432     Alignment *F=NULL;
1433     Sequence *STL;
1434     FILE *fp;
1435     char command[10000];
1436
1437     STL=(CL)?CL->STRUC_LIST:NULL;
1438
1439
1440     if ( !(CL->S) || !((CL->S)->T[s1]) || !((CL->S)->T[s1])->P || !seq2P_template_file(CL->S,s1))return CL;
1441     else pdb1=seq2P_template_file(CL->S,s1);
1442
1443
1444     pdb=vtmpnam (NULL);
1445     result=vtmpnam (NULL);
1446     pep=vtmpnam (NULL);
1447
1448     sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST  -nodiagnostic > %s", pdb1, pdb);
1449     my_system (command);
1450
1451
1452
1453     fp=vfopen (pep, "w");
1454     fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s2],(CL->S)->seq[s2] );
1455     vfclose (fp);
1456     sprintf ( command, "tc_generic_method.pl -mode=thread_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pep, M->in_flag2,pdb,M->out_flag, result, get_tmp_4_tcoffee());
1457     my_system ( command);
1458     F=main_read_aln (result, NULL);
1459
1460     if ( !F)
1461       {
1462         fprintf ( stderr, "\n\tthread_pair/%s failed:\n\t%s\n", M->executable2,command);
1463       }
1464     else
1465       {
1466         sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
1467         sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
1468         F=fix_aln_seq (F, CL->S);
1469         CL=aln2constraint_list (F, CL, "sim");
1470       }
1471
1472
1473     free_aln (F);
1474     return CL;
1475   }
1476
1477 Constraint_list * lsqman_pair ( char *in_seq, Constraint_list *CL)
1478         {
1479           FILE *fp;
1480           static CLIST_TYPE *entry;
1481           char command[STRING];
1482           char seq[1000];
1483           int s1, s2;
1484           char *seq_file, *lsqman_result, *tmp_name1;
1485           Alignment *F=NULL;
1486           int n_failure=0;
1487
1488           sprintf ( seq, "%s", in_seq);
1489
1490
1491           if ( !entry)entry=vcalloc ( LIST_N_FIELDS, sizeof ( CLIST_TYPE ));
1492
1493           atoi(strtok (seq,SEPARATORS));
1494           s1=atoi(strtok (NULL,SEPARATORS));
1495           s2=atoi(strtok (NULL,SEPARATORS));
1496
1497
1498          tmp_name1=vcalloc (100, sizeof (char));
1499          sprintf ( tmp_name1, "%s_%s.lsqman_aln", (CL->S)->name[s1], (CL->S)->name[s2]);
1500          if ( check_file_exists ( tmp_name1) && (F=main_read_aln(tmp_name1, NULL))!=NULL)
1501               {
1502                 free_aln(F);
1503                 lsqman_result=tmp_name1;
1504               }
1505
1506          else
1507            {
1508              seq_file=vtmpnam (NULL);
1509              lsqman_result=tmp_name1;
1510              fp=vfopen (seq_file, "w");
1511              fprintf ( fp, ">%s\n%s\n",(CL->S)->name[s1],(CL->S)->seq[s2] );
1512              vfclose (fp);
1513              sprintf ( command, "%s -pdb %s -pep %s > %s%s", LSQMAN_4_TCOFFEE, (CL->S)->name[s1], seq_file,get_cache_dir(), lsqman_result);
1514
1515              while (!F)
1516                {
1517                  my_system ( command);
1518                  F=main_read_aln (lsqman_result, NULL);
1519                  if ( !F)
1520                    {
1521                      fprintf ( stderr, "\n\tlsqman failed: will be retried");
1522                      if ( n_failure==0)fprintf ( stderr, "\n\t%s", command);
1523                      n_failure++;
1524                      if ( n_failure==10)
1525                         {
1526                         fprintf ( stderr, "\nCould not run Fugue: will replace it with slow_pair\n");
1527                         vremove (lsqman_result);
1528                         return NULL;
1529                         }
1530                    }
1531                  free_aln (F);
1532                }
1533              vremove( seq_file);
1534
1535            }
1536
1537
1538          F=main_read_aln(lsqman_result, NULL);
1539          /*sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
1540          sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
1541          */
1542          F=fix_aln_seq (F, CL->S);
1543
1544
1545
1546
1547
1548          CL=aln2constraint_list (F, CL, "100");
1549          free_aln (F);
1550          return CL;
1551         }
1552
1553 Constraint_list * sap_pair   (char *seq, char *weight, Constraint_list *CL)
1554  {
1555             register int a;
1556             FILE *fp;
1557             char full_name[FILENAMELEN];
1558             char *tmp_pdb1, *tmp_pdb2;
1559             char *sap_seq1, *sap_seq2;
1560             char *sap_lib,  *tmp_name, *tmp_name1, *tmp_name2;
1561             char *buf=NULL;
1562             int s1, s2, r1=0, r2=0;
1563             int sim=0, tot=0, score=0;
1564
1565             char program[STRING];
1566             char *string1, *string2, *string3, *string4, *string5;
1567             int max_struc_len=10000;
1568             char *template1, *template2;
1569
1570             /*check_program_is_installed ( "sap"                   ,SAP_4_TCOFFEE, "SAP_4_TCOFFEE",MAIL, IS_FATAL);*/
1571
1572
1573
1574
1575             atoi(strtok (seq,SEPARATORS));
1576             s1=atoi(strtok (NULL,SEPARATORS));
1577             s2=atoi(strtok (NULL,SEPARATORS));
1578
1579             template1=seq2T_value(CL->S,s1, "template_name", "_P_");
1580             template2=seq2T_value(CL->S,s2, "template_name", "_P_");
1581
1582
1583             if (!template1 || !template2) return CL;
1584
1585
1586             declare_name (string1);
1587             declare_name (string2);
1588             declare_name (string3);
1589             declare_name (string4);
1590             declare_name (string5);
1591
1592 #ifndef     SAP_4_TCOFFEE
1593             if ( getenv ( "SAP_4_TCOFFEE")==NULL)crash ("SAP_4_TCOFFEE IS NOT DEFINED");
1594             else  sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
1595 #else
1596             if ( getenv ( "SAP_4_TCOFFEE")==NULL)sprintf (program, "%s", SAP_4_TCOFFEE);
1597             else sprintf ( program, "%s", (getenv ( "SAP_4_TCOFFEE")));
1598 #endif
1599
1600
1601
1602             tmp_name1=vcalloc (100, sizeof (char));
1603             sprintf ( tmp_name1, "%s_%s.sap_results",template1,template2);
1604             tmp_name2=vcalloc (100, sizeof (char));
1605             sprintf ( tmp_name2, "%s_%s.sap_results",template2,template1);
1606
1607
1608
1609             if (is_sap_file (tmp_name1))
1610               {
1611                 tmp_name=tmp_name1;
1612               }
1613             else if ( is_sap_file (tmp_name2))
1614               {
1615                 tmp_name=tmp_name2;
1616                 SWAP (s1, s2);
1617               }
1618             else
1619               {
1620                 tmp_name=tmp_name1;
1621
1622                 tmp_pdb1=normalize_pdb_file(seq2P_template_file(CL->S,s1),(CL->S)->seq[s1], vtmpnam (NULL));
1623                 tmp_pdb2=normalize_pdb_file(seq2P_template_file(CL->S,s2),(CL->S)->seq[s2], vtmpnam (NULL));
1624                 sprintf ( full_name, "%s%s", get_cache_dir (), tmp_name);
1625                 printf_system ("%s %s %s >%s",program,tmp_pdb1,tmp_pdb2, full_name);
1626
1627                 if ( !check_file_exists (full_name) || !is_sap_file(full_name))
1628                   {
1629                     add_warning ( stderr, "WARNING: SAP failed to align: %s against %s [%s:WARNING]\n", seq2P_template_file(CL->S,s1),seq2P_template_file(CL->S,s2), PROGRAM);
1630                     if ( check_file_exists (full_name))add2file2remove_list (full_name);
1631                     return CL;
1632                   }
1633                 if ( flag_file2remove_is_on())add2file2remove_list (full_name);
1634                 remove ("super.pdb");
1635               }
1636
1637
1638
1639             sap_seq1=vcalloc (max_struc_len, sizeof (char));
1640             sap_seq2=vcalloc (max_struc_len, sizeof (char));
1641
1642
1643             fp=find_token_in_file ( tmp_name, NULL, "Percent");
1644             fp=find_token_in_file ( tmp_name, fp  , "Percent");
1645             while ( (fgetc (fp))!='\n');
1646             while ((buf=vfgets (buf, fp)))
1647               {
1648
1649                 if ( !strstr (buf, "eighted"))
1650                   {
1651                     remove_charset (buf, "!alnum");
1652                     r1=buf[0];
1653                     r2=buf[strlen(buf)-1];
1654                   }
1655
1656                 sim+=(r1==r2)?1:0;
1657                 if ( tot>max_struc_len)
1658                   {max_struc_len+=max_struc_len;
1659                     sap_seq1=vrealloc ( sap_seq1, sizeof(char)*max_struc_len);
1660                     sap_seq2=vrealloc ( sap_seq2, sizeof(char)*max_struc_len);
1661                   }
1662
1663                 sap_seq1[tot]=r1;
1664                 sap_seq2[tot]=r2;
1665                 tot++;
1666               }
1667             sim=(sim*100)/tot;
1668
1669             if ( is_number (weight))score=atoi(weight);
1670             else if ( strstr ( weight, "OW"))
1671               {
1672                 int ow;
1673                 sscanf ( weight, "OW%d", &ow);
1674                 score=sim*ow;
1675               }
1676             else
1677               {
1678                 score=sim;
1679               }
1680
1681
1682             vfclose (fp);
1683             sap_seq1[tot]=sap_seq2[tot]='\0';
1684
1685
1686             fp=vfopen ( sap_lib=vtmpnam(NULL), "w");
1687             fprintf (fp, "! TC_LIB_FORMAT_01\n");
1688             fprintf (fp, "2\n");
1689             fprintf (fp, "%s %d %s\n", (CL->S)->name[s2],(int)strlen (sap_seq1), sap_seq1);
1690             fprintf (fp, "%s %d %s\n", (CL->S)->name[s1],(int)strlen (sap_seq2), sap_seq2);
1691             fprintf (fp, "#1 2\n");
1692
1693             for ( a=0; a< tot; a++)
1694               {
1695
1696                 fprintf (fp, "%d %d %d 1 0\n", a+1, a+1, score);
1697
1698               }
1699             fprintf (fp, "! CPU 0\n");
1700             fprintf (fp, "! SEQ_1_TO_N\n");
1701             vfclose (fp);
1702             CL=read_constraint_list_file(CL,sap_lib);
1703             vremove (sap_lib);
1704
1705             vfree ( string1);vfree ( string2);vfree ( string3);vfree ( string4);vfree ( string5);
1706             vfree (sap_seq1); vfree(sap_seq2);vfree (tmp_name1); vfree(tmp_name2);
1707             vfree (buf);
1708             return CL;
1709         }
1710
1711
1712
1713 Constraint_list *rna_pair (TC_method *M ,
1714                                                    char *in_seq,
1715                                                    Constraint_list *CL)
1716 {
1717         char seq[1000];
1718         int s1, s2;
1719         char *result, *pdb1, *pdb2, *pdb1_file, *pdb2_file;
1720         Alignment *F=NULL;
1721
1722
1723
1724         char command[10000];
1725
1726         if ( M->executable2[0]=='\0')
1727         {
1728                 fprintf ( stderr, "\nERROR: rna_pair requires a structural alignment method: pdb_pair@EP@EXECUTABLE2@<method> [FATAL:%s]\n", PROGRAM);
1729                 myexit (EXIT_FAILURE);
1730         }
1731
1732         sprintf ( seq, "%s", in_seq);
1733
1734
1735         atoi(strtok (seq,SEPARATORS));
1736         s1=atoi(strtok (NULL,SEPARATORS));
1737         s2=atoi(strtok (NULL,SEPARATORS));
1738
1739         pdb1=seq2P_template_file(CL->S,s1);
1740         pdb1_file=vtmpnam (NULL);
1741         sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST  -nodiagnostic > %s", pdb1, pdb1_file);
1742         my_system (command);
1743         //
1744         pdb2=seq2P_template_file(CL->S,s2);
1745         pdb2_file=vtmpnam (NULL);
1746         sprintf ( command, "extract_from_pdb -infile %s -atom ALL -chain FIRST  -nodiagnostic > %s", pdb2, pdb2_file);
1747         my_system (command);
1748
1749         result=vtmpnam (NULL);
1750
1751         sprintf ( command, "tc_generic_method.pl -mode=rna_pair -method=%s %s%s %s%s %s%s -tmpdir=%s", M->executable2,M->in_flag,pdb1_file, M->in_flag2,pdb2_file,M->out_flag, result, get_tmp_4_tcoffee());
1752
1753         my_system ( command);
1754
1755         F=main_read_aln (result, NULL);
1756
1757
1758         if ( !F)
1759         {
1760                 fprintf ( stderr, "\n\trna_pair/%s failed:\n\t%s\n",M->executable2, command);
1761         }
1762         else
1763         {
1764                 sprintf ( F->name[0],"%s", (CL->S)->name[s1]);
1765                 sprintf ( F->name[1],"%s", (CL->S)->name[s2]);
1766                 F=fix_aln_seq (F, CL->S);
1767                 CL=aln2constraint_list (F, CL, "sim");
1768         }
1769
1770
1771         free_aln (F);
1772         return CL;
1773 }
1774
1775
1776 /******************************************************************/
1777 /*                   GENERIC PAIRWISE METHODS                     */
1778 /*                                                                */
1779 /*                                                                */
1780 /******************************************************************/
1781
1782
1783
1784 Constraint_list * best_pair4rna(Job_TC *job)
1785 {
1786         int n,a;
1787
1788
1789         static char *seq;
1790         Alignment *A;
1791         Constraint_list *PW_CL;
1792         Constraint_list *CL, *RCL;
1793         char *seq_in;
1794         Sequence *S;
1795         TC_method *M, *sara_pairM, *proba_pairM;
1796         int*seqlist;
1797
1798         int id, s1, s2;
1799         Template *T1, *T2;
1800         int ml=0;
1801         struct X_template *r1, *r2, *p1, *p2;
1802         static int **blosum;
1803
1804         if (!seq)seq=vcalloc (100, sizeof (char));
1805
1806         A=(job->io)->A;
1807         M=(job->param)->TCM;
1808         CL=(job->io)->CL;
1809         S=CL->S;
1810         for (a=0; a<S->nseq; a++)ml=MAX(ml, strlen (S->name[a]));
1811
1812
1813         if ( !strm ( retrieve_seq_type(), "RNA") )
1814                 printf_exit (EXIT_FAILURE, stderr, "ERROR: RNA Sequences Only with best4rna_pair  [FATAL:%s]\n",PROGRAM);
1815
1816
1817         seq_in=(job->param)->seq_c;
1818         sprintf (seq, "%s", seq_in);
1819         seqlist=string2num_list (seq);
1820         n=seqlist[1];
1821         if ( n!=2){fprintf ( stderr, "\nERROR: best_pair can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
1822         s1=seqlist[2];
1823         s2=seqlist[3];
1824
1825         T1=S->T[s1];
1826         T2=S->T[s2];
1827         r1=T1->R;
1828         r2=T2->R;
1829         p1=T1->P;
1830         p2=T2->P;
1831
1832         PW_CL=((job->param)->TCM)->PW_CL;
1833         CL=(job->io)->CL;
1834
1835         if (!blosum)blosum=read_matrice ("blosum62mt");
1836
1837 //      id=idscore_pairseq (S->seq[s1], S->seq[s2],-10,-1,blosum, "sim");
1838
1839
1840         proba_pairM=method_file2TC_method(method_name2method_file ("proba_pair"));
1841         proba_pairM->PW_CL=method2pw_cl(proba_pairM, CL);
1842
1843         sara_pairM=method_file2TC_method(method_name2method_file ("sara_pair"));
1844         sara_pairM->PW_CL=method2pw_cl(sara_pairM, CL);
1845
1846         if ( p1 && p2)
1847         {
1848       //Avoid Structural Tem
1849                 T1->R=NULL;
1850                 T2->R=NULL;
1851                 fprintf ( stderr, "\n\t%-*s %-*s: Structure Based Alignment\n", ml,S->name[s1], ml,S->name[s2]);
1852                 (job->param)->TCM=sara_pairM;
1853         }
1854         else
1855         {
1856                 fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment\n", ml,S->name[s1], ml,S->name[s2]);
1857                 (job->param)->TCM=proba_pairM;
1858         }
1859
1860         RCL=seq2list (job);
1861         T1->R=r1;
1862         T2->R=r2;
1863
1864         return RCL;
1865 }
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881 Constraint_list * best_pair4prot      (Job_TC *job)
1882 {
1883   int n,a;
1884
1885
1886   static char *seq;
1887   Alignment *A;
1888   Constraint_list *PW_CL;
1889   Constraint_list *CL, *RCL;
1890   char *seq_in;
1891   Sequence *S;
1892   TC_method *M, *sap_pairM, *proba_pairM;
1893   int*seqlist;
1894
1895   int id, s1, s2;
1896   Template *T1, *T2;
1897   int ml=0;
1898   struct X_template *r1, *r2, *p1, *p2;
1899   static int **blosum;
1900
1901   if (!seq)seq=vcalloc (100, sizeof (char));
1902
1903   A=(job->io)->A;
1904   M=(job->param)->TCM;
1905   CL=(job->io)->CL;
1906   S=CL->S;
1907   for (a=0; a<S->nseq; a++)ml=MAX(ml, strlen (S->name[a]));
1908
1909
1910   if ( strm ( retrieve_seq_type(), "DNA") ||strm ( retrieve_seq_type(), "RNA") )printf_exit (EXIT_FAILURE, stderr, "ERROR: Protein Sequences Only with bestprot_pair  [FATAL:%s]\n",PROGRAM);
1911
1912
1913   seq_in=(job->param)->seq_c;
1914   sprintf (seq, "%s", seq_in);
1915   seqlist=string2num_list (seq);
1916   n=seqlist[1];
1917   if ( n!=2){fprintf ( stderr, "\nERROR: best_pair can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
1918   s1=seqlist[2];
1919   s2=seqlist[3];
1920
1921   T1=S->T[s1];
1922   T2=S->T[s2];
1923   r1=T1->R;
1924   r2=T2->R;
1925   p1=T1->P;
1926   p2=T2->P;
1927
1928   PW_CL=((job->param)->TCM)->PW_CL;
1929   CL=(job->io)->CL;
1930
1931   if (!blosum)blosum=read_matrice ("blosum62mt");
1932
1933   id=idscore_pairseq (S->seq[s1], S->seq[s2],-10,-1,blosum, "sim");
1934
1935
1936   proba_pairM=method_file2TC_method(method_name2method_file ("proba_pair"));
1937   proba_pairM->PW_CL=method2pw_cl(proba_pairM, CL);
1938
1939   sap_pairM=method_file2TC_method(method_name2method_file ("sap_pair"));
1940   sap_pairM->PW_CL=method2pw_cl(sap_pairM, CL);
1941
1942   if ( id>80)
1943     {
1944       //Hide The Template
1945       T1->R=NULL;
1946       T2->R=NULL;
1947       fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment\n", ml,S->name[s1], ml,S->name[s2]);
1948       (job->param)->TCM=proba_pairM;
1949     }
1950   else if ( p1 && p2)
1951     {
1952       //Avoid Structural Tem
1953       T1->R=NULL;
1954       T2->R=NULL;
1955       fprintf ( stderr, "\n\t%-*s %-*s: Structure Based Alignment\n", ml,S->name[s1], ml,S->name[s2]);
1956       (job->param)->TCM=sap_pairM;
1957     }
1958   else if ( r1 || r2)
1959     {
1960       fprintf ( stderr, "\n\tt%-*s %-*s: PSIBLAST Profile Alignment\n", ml,S->name[s1], ml,S->name[s2]);
1961       (job->param)->TCM=proba_pairM;
1962     }
1963   else
1964     {
1965        fprintf ( stderr, "\n\t%-*s %-*s: Direct Sequence Alignment (No Profile)\n", ml,S->name[s1], ml,S->name[s2]);
1966        (job->param)->TCM=proba_pairM;
1967     }
1968
1969   RCL=seq2list (job);
1970   T1->R=r1;
1971   T2->R=r2;
1972
1973   return RCL;
1974 }
1975
1976
1977 Alignment * fast_pair      (Job_TC *job)
1978         {
1979             int s, n,a;
1980             static int **l_s;
1981             static int *ns;
1982             char seq[1000];
1983             Alignment *A;
1984             Constraint_list *PW_CL;
1985             Constraint_list *CL;
1986             char *seq_in;
1987             Sequence *S;
1988             TC_method *M;
1989             int*seqlist;
1990             char **buf;
1991
1992             A=(job->io)->A;
1993
1994             M=(job->param)->TCM;
1995             PW_CL=((job->param)->TCM)->PW_CL;
1996             CL=(job->io)->CL;
1997             seq_in=(job->param)->seq_c;
1998
1999
2000             sprintf (seq, "%s", seq_in);
2001             seqlist=string2num_list (seq);
2002             n=seqlist[1];
2003             if ( n!=2){fprintf ( stderr, "\nERROR: fast_pw_aln can only handle two seq at a time [FATAL]\n");myexit (EXIT_FAILURE);}
2004
2005             S=(CL)->S;
2006
2007             if (!A) {A=declare_aln (CL->S);}
2008             if ( !ns)
2009                 {
2010                 ns=vcalloc ( 2, sizeof (int));
2011                 l_s=declare_int (2,(CL->S)->nseq);
2012                 }
2013             buf=vcalloc ( S->nseq, sizeof (char*));
2014
2015             for ( a=0; a< n; a++)
2016                 {
2017                   s=seqlist[a+2];
2018                   if ( strm (M->seq_type, "G"))
2019                     {
2020                       buf[s]=S->seq[s];
2021                       S->seq[s]=((((S->T[s])->G)->VG)->S)->seq[0];
2022                   }
2023                 else
2024                   buf[s]=S->seq[s];
2025
2026
2027                   sprintf ( A->seq_al[a], "%s",S->seq[s]);
2028                   sprintf ( A->name[a], "%s", (CL->S)->name[s]);
2029                   A->order[a][0]=s;
2030                 }
2031
2032             A->S=CL->S;
2033             PW_CL->S=CL->S;
2034             A->CL=CL;
2035             ns[0]=ns[1]=1;
2036             l_s[0][0]=0;
2037             l_s[1][0]=1;
2038
2039
2040             if (PW_CL->reverse_seq)
2041               {
2042                 invert_string2(A->seq_al[0]);
2043                 invert_string2(A->seq_al[1]);
2044                 invert_string2 ((CL->S)->seq[A->order[0][0]]);
2045                 invert_string2 ((CL->S)->seq[A->order[1][0]]);
2046
2047               }
2048
2049             pair_wise ( A, ns, l_s, PW_CL);
2050
2051             if (PW_CL->reverse_seq)
2052               {
2053
2054                 invert_string2(A->seq_al[0]);
2055                 invert_string2(A->seq_al[1]);
2056                 invert_string2 ((CL->S)->seq[A->order[0][0]]);
2057                 invert_string2 ((CL->S)->seq[A->order[1][0]]);
2058               }
2059             A->nseq=n;
2060
2061             for ( a=0; a<S->nseq; a++)
2062               {
2063                 if ( !buf[a] || buf[a]==S->seq[a]);
2064                 else S->seq[a]=buf[a];
2065               }
2066             vfree (buf);vfree (seqlist);
2067             return A;
2068
2069         }
2070 Alignment * align_two_aln ( Alignment *A1, Alignment  *A2, char *in_matrix, int gop, int gep, char *in_align_mode)
2071         {
2072         Alignment *A=NULL;
2073         Constraint_list *CL;
2074         Sequence *S;
2075         int a;
2076         int *ns;
2077         int **ls;
2078         static char *matrix;
2079         static char *align_mode;
2080
2081         if (!matrix)matrix=vcalloc ( 100, sizeof (char));
2082         if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
2083         sprintf ( matrix, "%s", in_matrix);
2084         sprintf ( align_mode, "%s", in_align_mode);
2085
2086         CL=vcalloc ( 1, sizeof (Constraint_list));
2087         CL->pw_parameters_set=1;
2088         CL->M=read_matrice (matrix);
2089         CL->matrices_list=declare_char (10, 10);
2090
2091         CL->evaluate_residue_pair=evaluate_matrix_score;
2092         CL->get_dp_cost=consensus_get_dp_cost;
2093         CL->normalise=1;
2094
2095         CL->extend_jit=0;
2096         CL->maximise=1;
2097         CL->gop=gop;
2098         CL->gep=gep;
2099         CL->TG_MODE=2;
2100         sprintf (CL->matrix_for_aa_group, "vasiliky");
2101         CL->use_fragments=0;
2102         CL->ktup=5;
2103         if ( !CL->use_fragments)CL->diagonal_threshold=0;
2104         else CL->diagonal_threshold=6;
2105
2106         sprintf (CL->dp_mode, "%s", align_mode);
2107
2108         A=copy_aln (A1, A);
2109         A=stack_aln (A, A2);
2110         CL->S=fill_sequence_struc(A->nseq, A->seq_al,A->name);
2111
2112         ns=vcalloc ( 2, sizeof(int));
2113         ls=declare_int ( 2,A->nseq);
2114         ns[0]=A1->nseq;
2115         ns[1]=A2->nseq;
2116         for ( a=0; a<ns[0]; a++)
2117           ls[0][a]=a;
2118         for ( a=0; a< ns[1]; a++)
2119           ls[1][a]=a+A1->nseq;
2120
2121         A->score_aln=pair_wise (A, ns, ls,CL);
2122
2123         vfree (ns);
2124         free_int (ls, -1);
2125         S=free_constraint_list (CL);
2126         free_sequence (S,-1);
2127         A->S=NULL;
2128         return A;
2129         }
2130
2131 static int align_two_seq_keep_case;
2132 void toggle_case_in_align_two_sequences(int value)
2133 {
2134   align_two_seq_keep_case=value;
2135 }
2136 Alignment * align_two_sequences ( char *seq1, char *seq2, char *in_matrix, int gop, int gep, char *in_align_mode)
2137         {
2138         static Alignment *A;
2139         Constraint_list *CL;
2140         Sequence *S;
2141
2142         int *ns;
2143         int **l_s;
2144
2145         char       **seq_array;
2146         char       **name_array;
2147         static char *matrix;
2148         static int  **M;
2149
2150         static char *align_mode;
2151
2152         if (!matrix)matrix=vcalloc ( 100, sizeof (char));
2153         if (!align_mode)align_mode=vcalloc ( 100, sizeof (char));
2154         sprintf ( align_mode, "%s", in_align_mode);
2155
2156         CL=vcalloc ( 1, sizeof (Constraint_list));
2157
2158         CL->pw_parameters_set=1;
2159
2160         CL->matrices_list=declare_char (10, 10);
2161
2162
2163         if ( !strm (matrix, in_matrix))
2164           {
2165             sprintf ( matrix,"%s", in_matrix);
2166             M=CL->M=read_matrice (matrix);
2167
2168           }
2169         else
2170           {
2171             CL->M=M;
2172           }
2173
2174         if (strstr (in_align_mode, "cdna"))
2175           CL->evaluate_residue_pair=evaluate_cdna_matrix_score;
2176         else
2177           CL->evaluate_residue_pair=evaluate_matrix_score;
2178
2179         CL->get_dp_cost=get_dp_cost;
2180         CL->extend_jit=0;
2181         CL->maximise=1;
2182         CL->gop=gop;
2183         CL->gep=gep;
2184         CL->TG_MODE=2;
2185         sprintf (CL->matrix_for_aa_group, "vasiliky");
2186         CL->use_fragments=0;
2187         CL->ktup=3;
2188         if ( !CL->use_fragments)CL->diagonal_threshold=0;
2189         else CL->diagonal_threshold=6;
2190
2191         sprintf (CL->dp_mode, "%s", align_mode);
2192
2193         seq_array=declare_char ( 2, MAX(strlen(seq1), strlen (seq2))+1);
2194         sprintf (seq_array[0], "%s",seq1);
2195         sprintf (seq_array[1],"%s", seq2);
2196         ungap_array(seq_array,2);
2197         if (align_two_seq_keep_case !=KEEP_CASE)string_array_lower(seq_array,2);
2198
2199         name_array=declare_char (2, STRING);
2200         sprintf ( name_array[0], "A");
2201         sprintf ( name_array[1], "B");
2202
2203
2204         ns=vcalloc ( 2, sizeof(int));
2205         l_s=declare_int ( 2, 1);
2206         ns[0]=ns[1]=1;
2207         l_s[0][0]=0;
2208         l_s[1][0]=1;
2209
2210
2211
2212         CL->S=fill_sequence_struc(2, seq_array, name_array);
2213
2214         A=seq2aln(CL->S, NULL, 1);
2215
2216         ungap (A->seq_al[0]);
2217         ungap (A->seq_al[1]);
2218
2219
2220
2221         A->score_aln=pair_wise (A, ns, l_s,CL);
2222
2223         vfree (ns);
2224         free_int (l_s, -1);
2225         free_char (name_array, -1);free_char ( seq_array,-1);
2226
2227         CL->M=NULL;
2228         S=free_constraint_list (CL);
2229         free_sequence (S,-1);
2230         A->S=NULL;
2231         return A;
2232         }
2233
2234
2235 NT_node make_root_tree ( Alignment *A,Constraint_list *CL,int gop, int gep,Sequence *S,  char *tree_file,int maximise)
2236 {
2237    NT_node **T=NULL;
2238    T=make_tree (A, CL, gop, gep,S,tree_file,maximise);
2239    (T[3][0])->nseq=S->nseq;
2240    return T[3][0];
2241 }
2242 NT_node ** make_tree ( Alignment *A,Constraint_list *CL,int gop, int gep,Sequence *S,  char *tree_file,int maximise)
2243         {
2244           int a, b, ra, rb;
2245         NT_node **T=NULL;
2246         int **distances;
2247         int out_nseq;
2248         char **out_seq_name;
2249         char **out_seq;
2250
2251
2252         if ( !CL || !CL->tree_mode || !CL->tree_mode[0])
2253           {
2254             fprintf ( stderr, "\nERROR: No CL->tree_mode specified (make_tree::util_dp_drivers.c [FATAL:%s]", PROGRAM);
2255             myexit (EXIT_FAILURE);
2256           }
2257         else
2258           fprintf (CL->local_stderr , "\nMAKE GUIDE TREE \n\t[MODE=%s][",CL->tree_mode);
2259
2260
2261         if ( A->nseq==2)
2262           {
2263             int tot_node;
2264             char *tmp;
2265             FILE *fp;
2266             fprintf (CL->local_stderr, "---Two Sequences Only: Make Dummy Pair-Tree ---]");
2267             tmp=vtmpnam (NULL);
2268             fp=vfopen (tmp,"w");
2269             fprintf ( fp, "(%s:0.1, %s:0.1):0.1;\n",S->name[0], S->name[1]);
2270             vfclose (fp);
2271             T=read_tree (tmp, &tot_node, (CL->S)->nseq,(CL->S)->name);
2272
2273             return T;
2274           }
2275         else if ( strm (CL->tree_mode, "cwph"))
2276           {
2277             return  seq2cw_tree ( S, tree_file);
2278           }
2279         else if (strm ( CL->tree_mode, "upgma") || strm ( CL->tree_mode, "nj"))
2280           {
2281             out_nseq=S->nseq;
2282             out_seq_name=S->name;
2283             out_seq=S->seq;
2284
2285             CL->DM=cl2distance_matrix (CL, NOALN,NULL,NULL,0);
2286
2287             if ( CL->S!=S)
2288               {
2289                 /*Shrink the distance matrix so that it only contains the required sequences*/
2290                 distances=declare_int (S->nseq, S->nseq);
2291                 for (a=0; a< S->nseq; a++)
2292                   {
2293                     ra=name_is_in_list ((S)->name[a],(CL->S)->name, (CL->S)->nseq, 100);
2294                     for ( b=0; b< S->nseq; b++)
2295                       {
2296                         rb=name_is_in_list ((S)->name[b],(CL->S)->name, (CL->S)->nseq, 100);
2297                         distances[a][b]=(CL->DM)->score_similarity_matrix[ra][rb];
2298                       }
2299                   }
2300               }
2301             else
2302               {
2303                 distances=duplicate_int ( (CL->DM)->score_similarity_matrix, -1, -1);
2304               }
2305
2306
2307
2308             distances=sim_array2dist_array (distances, MAXID*SCORE_K);
2309             distances=normalize_array (distances, MAXID*SCORE_K, 100);
2310             if ( strm (CL->tree_mode, "order"))
2311               {
2312                 for ( a=0; a< S->nseq; a++)
2313                   for ( b=0; b< S->nseq; b++)
2314                     distances[b][a]=100;
2315                 T=make_nj_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
2316               }
2317             else if ( strm (CL->tree_mode, "nj"))
2318               {
2319                 T=make_nj_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
2320               }
2321             else if ( strm (CL->tree_mode, "upgma"))
2322               T=make_upgma_tree (A,distances,gop,gep,out_seq,out_seq_name,out_nseq, tree_file, CL->tree_mode);
2323             else
2324               {
2325                 printf_exit (EXIT_FAILURE, stderr, "ERROR: %s is an unknown tree computation mode [FATAL:%s]", CL->tree_mode, PROGRAM);
2326               }
2327             free_int (distances, out_nseq);
2328
2329           }
2330
2331         fprintf (CL->local_stderr , "DONE]\n");
2332         return T;
2333         }
2334
2335
2336 Alignment *recompute_local_aln (Alignment *A, Sequence *S,Constraint_list *CL, int scale, int gep)
2337     {
2338     int **coor;
2339     int a;
2340     Alignment *B;
2341
2342     sort_constraint_list (CL, 0, CL->ne);
2343     coor=declare_int (A->nseq, 3);
2344     for ( a=0; a< A->nseq; a++)
2345         {
2346         coor[a][0]=A->order[a][0];
2347         coor[a][1]=A->order[a][1]+1;
2348         coor[a][2]=strlen(S->seq[A->order[a][0]])-coor[a][1];
2349         }
2350     B=stack_progressive_nol_aln_with_seq_coor(CL,0,0,S,coor,A->nseq);
2351     A=copy_aln ( B, A);
2352
2353     CL=compact_list (CL, 0,CL->ne, "shrink");
2354     free_Alignment(B);
2355     return A;
2356     }
2357
2358
2359 Alignment *stack_progressive_nol_aln_with_seq_coor(Constraint_list *CL,int gop, int gep,Sequence *S, int **seq_coor, int nseq)
2360     {
2361
2362     static int ** local_coor1;
2363     static int ** local_coor2;
2364     if ( local_coor1!=NULL)free_int (local_coor1, -1);
2365     if ( local_coor2!=NULL)free_int (local_coor2, -1);
2366
2367     local_coor1=get_nol_seq          ( CL,seq_coor, nseq, S);
2368     local_coor2=minimise_repeat_coor ( local_coor1, nseq, S);
2369
2370     return stack_progressive_aln_with_seq_coor(CL,gop, gep,S, local_coor2,nseq);
2371     }
2372
2373
2374 Alignment *stack_progressive_aln_with_seq_coor (Constraint_list*CL,int gop, int gep, Sequence *S, int **coor, int nseq)
2375     {
2376     Alignment *A=NULL;
2377
2378     A=seq_coor2aln (S,NULL, coor, nseq);
2379
2380     return stack_progressive_aln ( A,CL, gop, gep);
2381     }
2382
2383 Alignment *est_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
2384     {
2385     int a,n;
2386     int**group_list;
2387     int *n_groups;
2388     char *seq;
2389     n_groups=vcalloc ( 2, sizeof (int));
2390     group_list=declare_int ( 2, A->nseq);
2391
2392     n=A->nseq;
2393
2394     n_groups[0]=1;
2395     n_groups[1]=1;
2396     group_list[0][0]=0;
2397     group_list[0][1]=1;
2398
2399     group_list[1][0]=1;
2400     fprintf ( stderr, "\n");
2401     for ( a=1; a<n; a++)
2402         {
2403         sprintf ( A->seq_al[1], "%s", A->seq_al[a]);
2404         fprintf ( stderr, "\t[%30s]->[len=%5d]", A->name[a],(int)strlen ( A->seq_al[0]));
2405         pair_wise ( A,n_groups, group_list, CL);
2406
2407         seq=dna_aln2cons_seq(A);
2408
2409         sprintf ( A->seq_al[0], "%s", seq);
2410         vfree (seq);
2411         fprintf ( stderr, "\n");
2412         }
2413
2414     A->nseq=1;
2415     return A;
2416     }
2417
2418 void analyse_seq ( Alignment *A, int s)
2419    {
2420      int a, b, c;
2421      int r;
2422
2423      int len=0;
2424
2425
2426      int state=0;
2427      int pstate=-1;
2428      float score=0;
2429
2430      for ( a=0; a< A->len_aln; a++)
2431        {
2432          for ( b=0, c=0; b< s; b++)
2433            if ( !is_gap(A->seq_al[b][a])){c=1; break;}
2434
2435          r=!is_gap(A->seq_al[s][a]);
2436
2437          if      (  r &&  c) state=1;
2438          else if ( !r && !c) state=2;
2439          else if ( !r &&  c) state=3;
2440          else if (  r && !c) state=4;
2441
2442          if ( state !=pstate)
2443            {
2444              score+=len*len;
2445              len=0;
2446            }
2447          len+=r;
2448          pstate=state;
2449        }
2450      score=score/(float)(((A->S)->len[s]*(A->S)->len[s]));
2451      fprintf ( stderr, "[%.2f]", score);
2452
2453      return;
2454    }
2455
2456 Alignment *realign_aln ( Alignment*A, Constraint_list *CL)
2457 {
2458   int a, b, c;
2459   int *ns, **ls;
2460   A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
2461
2462   ns=vcalloc (2, sizeof(int));
2463   ls=declare_int ( 2, A->nseq);
2464
2465   for (a=0; a< A->nseq; a++)
2466     {
2467       ns[0]=A->nseq-1;
2468       for (c=0,b=0; b<A->nseq; b++)if (b!=a)ls[0][c++]=b;
2469       ungap_sub_aln ( A, ns[0], ls[0]);
2470
2471       ns[1]=1;
2472       ls[1][0]=a;
2473       ungap_sub_aln ( A, ns[1], ls[1]);
2474       A->score_aln=pair_wise (A, ns, ls,CL);
2475     }
2476
2477   vfree (ns); free_int (ls, -1);
2478   return A;
2479 }
2480
2481 Alignment *realign_aln_random_bipart ( Alignment*A, Constraint_list *CL)
2482 {
2483   int *ns;
2484   int **l_s;
2485
2486   int  a,g;
2487
2488   ns=vcalloc (2, sizeof (int));
2489   l_s=declare_int (2,A->nseq);
2490
2491   for ( a=0; a< A->nseq; a++)
2492     {
2493     g=rand()%2;
2494     l_s[g][ns[g]++]=a;
2495     }
2496
2497   fprintf ( stderr, "\n");
2498   ungap_sub_aln ( A, ns[0], l_s[0]);
2499   ungap_sub_aln ( A, ns[1], l_s[1]);
2500
2501   /* //Display Groups
2502   for (a=0;a<2; a++)
2503     for (b=0; b<ns[a]; b++)
2504       fprintf ( stderr, "[%d-%d]", a, l_s[a][b]);
2505   */
2506   A->score_aln=pair_wise (A, ns, l_s,CL);
2507
2508   vfree(ns);free_int(l_s, -1);
2509   return A;
2510 }
2511 Alignment *realign_aln_random_bipart_n ( Alignment*A, Constraint_list *CL, int n)
2512 {
2513   int *ns;
2514   int **ls;
2515   int *used;
2516
2517   int  a,b,c, p;
2518
2519   if (n>=A->nseq)n=A->nseq/2;
2520   used=vcalloc (A->nseq, sizeof (int));
2521   c=0;
2522   while (c<n)
2523     {
2524       p=rand()%A->nseq;
2525       if (!used[p]){used[p]=1;c++;}
2526     }
2527   ns=vcalloc (2, sizeof (int));
2528   ls=declare_int (2,A->nseq);
2529   for (a=0; a<2; a++)
2530     {
2531       for (b=0; b<A->nseq; b++)
2532         if (used[b]==a)ls[a][ns[a]++]=b;
2533     }
2534   ungap_sub_aln ( A, ns[0], ls[0]);
2535   ungap_sub_aln ( A, ns[1], ls[1]);
2536
2537
2538   A->score_aln=pair_wise (A, ns, ls,CL);
2539   vfree(ns);free_int(ls, -1);vfree (used);
2540   return A;
2541 }
2542 int ** seq2ecl_mat (Constraint_list *CL);
2543 int ** seq2ecl_mat (Constraint_list *CL)
2544
2545 {
2546   int a, b, n;
2547
2548   Alignment *A;
2549   int *ns, **ls;
2550   int **dm;
2551
2552   ns=vcalloc (2, sizeof (int));
2553   ls=declare_int ((CL->S)->nseq, 2);
2554
2555   A=seq2aln (CL->S,NULL, RM_GAP);
2556   n=(CL->S)->nseq;
2557   dm=declare_int (n, n);
2558   for (a=0; a<(CL->S)->nseq-1; a++)
2559     for (b=a+1; b<(CL->S)->nseq; b++)
2560       {
2561         ns[0]=ns[1]=1;
2562         ls[0][0]=a;
2563         ls[1][0]=b;
2564         ungap (A->seq_al[a]);
2565         ungap (A->seq_al[b]);
2566         dm[a][b]=dm[b][a]=linked_pair_wise (A, ns, ls, CL);
2567       }
2568
2569   return dm;
2570 }
2571 Alignment *realign_aln_clust ( Alignment*A, Constraint_list *CL)
2572 {
2573   int *ns;
2574   int **ls;
2575
2576   int a, b, c,n;
2577   static int **rm, **dm, **target;
2578   int score;
2579
2580
2581
2582   if (!A)
2583     {
2584       free_int (dm, -1); free_int (rm, -1);free_int (target, -1);
2585       dm=rm=target=NULL;
2586     }
2587
2588
2589   if (!rm)rm=seq2ecl_mat(CL);
2590   if (!dm)dm=declare_int (A->nseq, A->nseq);
2591   if (!target)target=declare_int (A->nseq*A->nseq, 3);
2592
2593   ns=vcalloc (2, sizeof (int));
2594   ls=declare_int (2,A->nseq);
2595
2596
2597   for (a=0; a<A->nseq-1; a++)
2598     for (b=a+1; b<A->nseq; b++)
2599       {
2600         ns[0]=2;
2601         ls[0][0]=a;
2602         ls[0][1]=b;
2603         score=sub_aln2ecl_raw_score (A, CL, ns[0], ls[0]);
2604         dm[a][b]=dm[b][a]=MAX(0,(rm[a][b]-score));
2605       }
2606   for (n=0,a=0; a<A->nseq; a++)
2607     {
2608       for (b=a; b<A->nseq; b++, n++)
2609         {
2610
2611           target[n][0]=a;
2612           target[n][1]=b;
2613           for ( c=0; c<A->nseq; c++)
2614             {
2615               if (c!=a && c!=b)target[n][2]+=dm[a][c]+dm[b][c];
2616             }
2617         }
2618     }
2619   sort_int_inv (target,3, 2, 0, n-1);
2620
2621   for (a=0; a<A->nseq; a++)
2622     {
2623       if (target[a][0]==target[a][1])
2624         {
2625           ns[0]=1;
2626           ls[0][0]=target[a][0];
2627         }
2628       else
2629         {
2630           ns[0]=2;
2631           ls[0][0]=target[a][0]; ls[0][1]=target[a][1];
2632         }
2633
2634       for (ns[1]=0,b=0; b<A->nseq; b++)
2635         {
2636           if (b!=target[a][0] && b!=target[a][1])ls[1][ns[1]++]=b;
2637         }
2638
2639       ungap_sub_aln (A, ns[0], ls[0]);
2640       ungap_sub_aln (A, ns[1], ls[1]);
2641
2642       A->score_aln=pair_wise (A, ns, ls,CL);
2643       fprintf ( stderr, "\nSEQ: %d %d SCORE=%d\n",target[a][0],target[a][1], aln2ecl_raw_score(A, CL));
2644     }
2645   return A;
2646 }
2647
2648 int get_best_group ( int **used, Constraint_list *CL);
2649 int seq_aln_thr1(Alignment *A, int **used, int threshold, Constraint_list *CL);
2650 int seq_aln_thr2( Alignment*A, int **used, int threshold, int g, Constraint_list *CL);
2651
2652 int get_best_group ( int **used, Constraint_list *CL)
2653 {
2654   int a,b,c,d,n, tot,stot, best_tot, best_seq, nseq;
2655   int ns[2];
2656   int *ls[2];
2657
2658   best_seq=0;
2659   nseq=((CL->S)->nseq);
2660   tot=best_tot=0;
2661   for (a=0; a<nseq; a++)
2662     {
2663       if ( used[a][0]==-1)continue;
2664       for ( tot=0,b=0; b< nseq; b++)
2665         {
2666           if ( a==b) continue;
2667           if ( used[b][0]==-1)continue;
2668           ns[0]=used[a][1];
2669           ls[0]=used[a]+2;
2670           ns[1]=used[b][1];
2671           ls[1]=used[b]+2;
2672           for (stot=0, n=0,c=0; c<ns[0]; c++)
2673             for (d=0; d<ns[1]; d++, n++)
2674               {
2675                 stot+=(CL->DM)->similarity_matrix[ls[0][c]][ls[1][d]];
2676               }
2677           if (n>0)stot/=n;
2678           tot+=stot;
2679         }
2680       if (tot>best_tot)
2681         {
2682           best_tot=tot;
2683           best_seq=a;
2684         }
2685     }
2686   return best_seq;
2687 }
2688
2689
2690 char ** list_file2dpa_list_file (char **list, int *len,int maxnseq, Sequence *S)
2691 {
2692   char **nlist, **profile_list;
2693   int nl, l, a, np, has_lib;
2694   char *seq, *profile;
2695   Alignment *A, *F;
2696
2697
2698   nlist=declare_char (read_array_size_new ((void *)list), read_array_size_new ((void *)list[0]));
2699   nl=0;
2700
2701   profile_list=declare_char (read_array_size_new ((void *)list), read_array_size_new ((void *)list[0]));
2702   np=0;
2703
2704   l=len[0];
2705   for (a=0; a<l; a++)if (list[a][0]!='S')sprintf ( nlist[nl++], "%s", list[a]);
2706
2707   profile=vtmpnam (NULL);
2708   seq=vtmpnam (NULL);
2709
2710
2711   output_fasta_seq (seq, A=seq2aln (S,NULL,RM_GAP));
2712   free_aln (A); A=NULL;
2713   printf_system ( "t_coffee -seq %s -msa_mode groups -outfile=%s -maxnseq=0 -in Xblosum62mt -dp_mode gotoh_pair_wise_lgp -dpa_maxnseq=%d", seq, profile, maxnseq);
2714
2715   if ( !check_file_exists(profile))
2716   {
2717     char command[10000];
2718     sprintf( command,"t_coffee -seq %s -msa_mode groups -outfile=%s -quiet=xxtest -in Xblosum62mt -distance_matrix_mode ktup -dpa_maxnseq=%d -maxnseq=0", seq, profile, maxnseq);
2719     printf_exit ( EXIT_FAILURE,stderr, "Attempt to switch to DPA failed: Could not run\n%s\n[FATAL:%s]\n", command, PROGRAM);
2720   }
2721   F=A=input_conc_aln (profile, NULL);
2722
2723   while (A)
2724     {
2725       char *file;
2726       FILE *fp;
2727
2728       file=vtmpnam (NULL);
2729       fp=vfopen (file, "w");
2730       for (a=0; a<A->nseq; a++)fprintf ( fp, ">%s %s\n%s\n", A->name[a],A->seq_comment[a], A->seq_al[a]);
2731       vfclose (fp);
2732
2733       sprintf (profile_list[np++], "%s", file);
2734       A=A->A;
2735     }
2736   free_aln (F);
2737
2738   for (has_lib=0,a=0; a<l; a++)
2739     {
2740       if ( list[a][0]=='A' || list[a][0]=='L')
2741         has_lib=1;
2742     }
2743
2744   if ( has_lib==1)
2745     nlist=reindex_constraint_list ( profile_list, np, nlist, &nl, S);
2746
2747   for (a=0; a<np; a++)
2748     {
2749       sprintf (nlist[nl++], "R%s", profile_list[a]);
2750     }
2751   for (a=0; a< nl; a++)
2752     {
2753       sprintf (list[a], "%s", nlist[a]);
2754     }
2755   len[0]=nl;
2756   free_char (profile_list, -1);
2757   free_char (nlist, -1);
2758   return list;
2759 }
2760
2761 Alignment * seq2aln_group (Alignment *A, int N, Constraint_list *CL)
2762 {
2763
2764   NT_node P;
2765   int a;
2766
2767   char *list, **list2;
2768   Alignment *F;
2769
2770
2771   fprintf (CL->local_stderr, "\n##### DPA ##### Compute Fast Alignment");
2772   A=iterative_tree_aln (A,1, CL);
2773   fprintf (CL->local_stderr, "\n##### DPA ##### Identify Nodes");
2774   P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
2775   set_node_score (A, P, "idmat_sim");
2776   fprintf (CL->local_stderr, "\n##### DPA ##### Split Nodes");
2777   list=split_nodes_nseq (A,P,N, list=vcalloc (P->nseq*200, sizeof (char)));
2778
2779   list2=string2list (list);
2780   fprintf (CL->local_stderr, "\n##### DPA ##### Save Nodes");
2781
2782   F=A;
2783   for (a=1; a<atoi (list2[0]); a++)
2784     {
2785       A->A=main_read_aln(list2[a], NULL);
2786       A=A->A;
2787     }
2788   fprintf (CL->local_stderr, "\n##### DPA ##### Finished");
2789   vfree (list); free_char (list2, -1);
2790
2791   A=F;
2792   while (A)
2793     {
2794             A=A->A;
2795     }
2796
2797   return F;
2798 }
2799
2800
2801
2802
2803 Alignment * seq_aln ( Alignment*A, int n,Constraint_list *CL)
2804 {
2805
2806   int **used,  a, t,n1, nseq;
2807
2808
2809   n1=nseq=(CL->S)->nseq;
2810   used=declare_int (nseq, nseq+3);
2811
2812
2813   for (a=0; a< nseq; a++)
2814     {
2815       used[a][1]=1;
2816       used[a][2]=a;
2817     }
2818
2819
2820   for (t=50; t>=0 && nseq>1; t-=5)
2821     {
2822       nseq=seq_aln_thr1 (A, used,t, CL);
2823     }
2824
2825   vfree (used);
2826   return A;
2827 }
2828
2829 int seq_aln_thrX(Alignment *A, int **used, int threshold, Constraint_list *CL)
2830 {
2831   int n=0,a;
2832   seq_aln_thr1(A,used,threshold,CL);
2833   for ( a=0; a< (CL->S)->nseq; a++)
2834     n+=(used[a][1]>0)?1:0;
2835
2836   return n;
2837 }
2838 int seq_aln_thr1(Alignment *A, int **used, int threshold, Constraint_list *CL)
2839 {
2840   int a,g, nseq, n_groups;
2841   nseq=(CL->S)->nseq;
2842
2843   g=get_best_group(used, CL);
2844
2845   used[g][0]=1;
2846
2847
2848
2849   while ( seq_aln_thr2 (A, used, threshold,g, CL)!=0)
2850     {
2851       g=get_best_group (used, CL);
2852       used[g][0]=1;
2853     }
2854
2855   for (n_groups=0,a=0; a< nseq; a++)
2856     if ( used[a][1]!=0)
2857       {
2858         n_groups++;
2859         used[a][0]=0;
2860       }
2861   return n_groups;
2862 }
2863
2864
2865 int seq_aln_thr2( Alignment*A, int **used, int threshold, int g, Constraint_list *CL)
2866 {
2867   int a, b,c,d;
2868   int ns[2], *ls[2];
2869   int nseq, n_members;
2870   double sim;
2871
2872   n_members=0;
2873
2874   nseq=((CL->S)->nseq);
2875   used[g][0]=1;
2876   ns[0]=used[g][1];
2877   ls[0]=used[g]+2;
2878
2879   for ( a=0; a< nseq; a++)
2880     {
2881       if (used[a][0]!=0);
2882       else
2883         {
2884           ns[1]=used[a][1];
2885           ls[1]=used[a]+2;
2886
2887           ungap_sub_aln (A, ns[0], ls[0]);
2888           ungap_sub_aln (A, ns[1], ls[1]);
2889
2890           A->score_aln=pair_wise (A, ns, ls,CL);
2891
2892           for (sim=0,b=0; b<ns[0]; b++)
2893             {
2894               for (c=0; c<ns[1]; c++)
2895                 {
2896                   sim+=generic_get_seq_sim (A->seq_al[ls[0][b]], A->seq_al[ls[1][c]], NULL,"idmat_sim2");
2897                 }
2898             }
2899           sim/=(double)(ns[0]*ns[1]);
2900           if (sim>=threshold)
2901             {
2902
2903               used[g][1]+=ns[1];
2904               for (d=0; d<ns[1]; d++)
2905                 ls[0][ns[0]++]=ls[1][d];
2906               used[a][0]=-1;
2907               used[a][1]=0;
2908               b=ns[0];c=ns[1];
2909               n_members++;
2910             }
2911
2912         }
2913     }
2914
2915   if (n_members>0)used[g][0]=-1;
2916   return n_members;
2917 }
2918 /****************************************************************************/
2919 /*                                                                          */
2920 /*                                                                          */
2921 /*            Alignment Methods                                             */
2922 /*                                                                          */
2923 /*                                                                          */
2924 /****************************************************************************/
2925
2926 Alignment * tsp_aln (Alignment *A, Constraint_list *CL, Sequence *S)
2927 {
2928   int a, b   ;
2929   int **  distances;
2930   int *ns, **ls;
2931   int **used;
2932
2933   A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
2934   ns=vcalloc (2, sizeof (int));
2935   ls=declare_int (2, (CL->S)->nseq);
2936   used=declare_int ( A->nseq, 2);
2937
2938
2939   CL->DM=cl2distance_matrix (CL, NOALN,NULL,NULL,0);
2940   distances=declare_int (A->nseq+1, A->nseq+1);
2941   distances=duplicate_int ( (CL->DM)->score_similarity_matrix, -1, -1);
2942
2943   for (a=0; a< A->nseq; a++)
2944     {
2945       used[a][0]=a;
2946       for (b=0; b< A->nseq; b++)
2947         {
2948           used[a][1]+=distances[a][b];
2949         }
2950     }
2951
2952   sort_int_inv (used,2,1,0,(CL->S)->nseq-1);
2953
2954   ls[0][ns[0]++]=used[0][0];
2955   ns[1]=1;
2956
2957   for (a=1; a< S->nseq; a++)
2958     {
2959       fprintf ( stderr, "\n%s %d", (CL->S)->name[used[a][0]], used[a][1]);
2960       ls[1][0]=used[a][0];
2961       pair_wise ( A,ns,ls, CL);
2962       ls[0][ns[0]++]=used[a][0];
2963     }
2964
2965   A->nseq=(CL->S)->nseq;
2966   return A;
2967
2968 }
2969
2970 Alignment *stack_progressive_aln(Alignment *A, Constraint_list *CL, int gop, int gep)
2971     {
2972     int a,n;
2973     int**group_list;
2974     int *n_groups;
2975     char dp_mode[100];
2976
2977
2978     sprintf ( dp_mode, "%s", CL->dp_mode);
2979     sprintf (CL->dp_mode, "gotoh_pair_wise");
2980
2981     n_groups=vcalloc ( 2, sizeof (int));
2982     group_list=declare_int ( 2, A->nseq);
2983
2984     n=A->nseq;
2985
2986     for ( a=0; a<n; a++)ungap(A->seq_al[a]);
2987     for ( a=1; a<n; a++)
2988         {
2989         n_groups[0]=a;
2990         n_groups[1]=1;
2991         group_list[0][a-1]=a-1;
2992         group_list[1][0]  =a;
2993
2994         pair_wise ( A,n_groups, group_list, CL);
2995         fprintf ( stderr, "\n\t[%d]->[%d]", a,(int)strlen ( A->seq_al[0]));
2996         }
2997     fprintf (stderr, "\n");
2998     vfree(n_groups);
2999     free_int ( group_list, -1);
3000     sprintf (CL->dp_mode, "%s",dp_mode);
3001
3002     return A;
3003     }
3004 Alignment *realign_aln_clust ( Alignment*A, Constraint_list *CL);
3005 Alignment *realign_aln_random_bipart_n ( Alignment*A, Constraint_list *CL, int n);
3006 Alignment *iterate_aln ( Alignment*A, int nit, Constraint_list *CL)
3007 {
3008   int it;
3009   int mode=1;
3010   int score, iscore, delta;
3011   fprintf ( CL->local_stderr, "Iterated Refinement: %d cycles START: score= %d\n", nit,iscore=aln2ecl_raw_score (A, CL) );
3012
3013
3014   if ( nit==-1)nit=A->nseq*2;
3015   if ( A->len_aln==0)A=very_fast_aln (A, A->nseq, CL);
3016   A=reorder_aln (A,(CL->S)->name, A->nseq);
3017
3018   for (it=0; it< nit; it++)
3019     {
3020       //CL->local_stderr=output_completion (CL->local_stderr,it, nit,1, "");
3021       if (mode==0)A=realign_aln (A, CL);
3022       else if (mode ==1)A=realign_aln_random_bipart (A, CL);
3023       else if (mode ==2)A=realign_aln_clust (A, CL);
3024       else if (mode ==3)A=realign_aln_random_bipart_n (A, CL,2);
3025
3026
3027       score=aln2ecl_raw_score (A, CL);
3028       delta=iscore-score;
3029       fprintf (CL->local_stderr, "\n\tIteration Cycle: %d Score=%d Improvement= %d", it+1,score, delta);
3030     }
3031   fprintf ( CL->local_stderr, "\nIterated Refinement: Completed Improvement=%d\n", delta);
3032   return A;
3033 }
3034
3035 int get_next_best (int seq, int nseq, int *used, int **dm);
3036 int get_next_best (int seq, int nseq, int *used, int **dm)
3037 {
3038   int a,set, d, bd, bseq;
3039
3040   for (set=0,a=0; a< nseq; a++)
3041     {
3042       if (used[a] || seq==a)continue;
3043       d=dm[seq][a];
3044       if (set==0 || d>bd)
3045         {
3046           bseq=a;
3047           bd=d;
3048           set=1;
3049         }
3050     }
3051   return bseq;
3052 }
3053 Alignment  * full_sorted_aln (Alignment *A, Constraint_list *CL)
3054 {
3055   int a,b;
3056   A=sorted_aln_seq (0, A, CL);
3057   print_aln(A);
3058   for (a=1; a<A->nseq; a++)
3059     {
3060       A=A->A=copy_aln (A, NULL);
3061       for (b=0; b<A->nseq; b++)ungap(A->seq_al[b]);
3062       A=sorted_aln_seq (a, A, CL);
3063       print_aln(A);
3064     }
3065   return A;
3066 }
3067 Alignment * sorted_aln (Alignment *A, Constraint_list *CL)
3068 {
3069   return sorted_aln_seq (-1, A, CL);
3070 }
3071 Alignment * sorted_aln_seq (int new_seq, Alignment *A, Constraint_list *CL)
3072 {
3073   int a, b=0, nseq;
3074   int *ns, **ls, **score, *used, **dm;
3075   int old_seq;
3076
3077   dm=(CL->DM)->score_similarity_matrix;
3078   nseq=(CL->S)->nseq;
3079   score=declare_int (nseq, 3);
3080   used=vcalloc (nseq, sizeof (int));
3081   ls=declare_int (2, nseq);
3082   ns=vcalloc (2, sizeof (int));
3083
3084
3085   if ( new_seq==-1)
3086     {
3087       for (a=0; a<nseq; a++)
3088         {
3089           score[a][0]=a;
3090           score[a][1]=b;
3091           for ( b=0; b<nseq; b++)
3092             score[a][2]+=dm[a][b];
3093         }
3094       sort_int ( score,3, 2, 0, nseq-1);
3095       old_seq=new_seq=score[nseq-1][0];
3096     }
3097
3098   for (a=1; a< nseq; a++)
3099     {
3100       used[new_seq]=1;
3101       ls[0][ns[0]++]=new_seq;
3102       ns[1]=1;
3103       ls[1][0]=get_next_best(new_seq,nseq, used,dm);
3104       old_seq=new_seq;
3105       new_seq=ls[1][0];
3106
3107       A->score_aln=pair_wise (A, ns, ls,CL);
3108
3109     }
3110   return A;
3111 }
3112
3113 Alignment * ungap_aln4tree (Alignment *A);
3114 Alignment * ungap_aln4tree (Alignment *A)
3115 {
3116   int t, n, max_sim, sim;
3117   Alignment *B;
3118
3119
3120
3121   n=35;
3122   max_sim=60;
3123
3124   t=A->len_aln/10;
3125
3126   B=copy_aln (A, NULL);
3127   B=ungap_aln_n(B, n);
3128   return B;
3129
3130   sim=aln2sim (B, "idmat");
3131   while (B->len_aln<t && sim>max_sim && n>0)
3132     {
3133       n-=10;
3134       B=copy_aln (A, B);
3135       B=ungap_aln_n(B, n);
3136       sim=aln2sim (B, "idmat");
3137     }
3138   if ( B->len_aln<t && sim>max_sim)B=copy_aln (A, B);
3139   return B;
3140 }
3141
3142
3143
3144 Alignment * iterative_tree_aln (Alignment *A,int n, Constraint_list *CL)
3145 {
3146   NT_node **T=NULL;
3147   int a;
3148
3149   T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
3150   tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
3151   for ( a=0; a< n; a++)
3152     {
3153
3154       Alignment *B;
3155
3156       B=copy_aln (A, NULL);
3157       B=ungap_aln_n (B, 20);
3158       sprintf ( CL->distance_matrix_mode, "aln");
3159
3160       CL->DM=cl2distance_matrix ( CL,B,NULL,NULL, 1);
3161       free_aln (B);
3162
3163       degap_aln (A);
3164       T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
3165
3166       tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
3167     }
3168   return A;
3169 }
3170
3171 Alignment *profile_aln (Alignment *A, Constraint_list *CL)
3172 {
3173   int a,nseq,nseq2;
3174   int **ls, *ns;
3175   nseq=A->nseq;
3176   nseq2=2*nseq;
3177   ls=declare_int (2, nseq2);
3178   ns=vcalloc (2, sizeof (int));
3179
3180   A=realloc_aln2(A,nseq2, A->len_aln);
3181   for (a=0; a< nseq; a++)
3182     ls[0][ns[0]++]=a;
3183   for ( a=0; a<nseq;a++)
3184     {
3185       sprintf (A->seq_al[a+nseq], "%s", (CL->S)->seq[a]);
3186       sprintf (A->name[a+nseq], "%s", (CL->S)->name[a]);
3187       A->order[a+nseq][0]=a;
3188     }
3189
3190   ns[1]=1;
3191   for (a=0; a<nseq; a++)
3192     {
3193
3194       ls[1][0]=a+nseq;
3195       A->score_aln=pair_wise (A, ns, ls,CL);
3196
3197       ls[0][ns[0]++]=a+nseq;
3198     }
3199   for (a=0; a< nseq; a++)
3200     {
3201       sprintf (A->seq_al[a], "%s", A->seq_al[a+nseq]);
3202     }
3203   A->nseq=nseq;
3204   return A;
3205 }
3206
3207 Alignment * iterative_aln ( Alignment*A, int n,Constraint_list *CL)
3208 {
3209   int *ns,**ls, **score, **dm;
3210   int a,b, nseq, max;
3211   ls=declare_int (2, A->nseq);
3212   ns=vcalloc (2, sizeof (int));
3213   ls[0][ns[0]++]=0;
3214
3215
3216
3217
3218
3219   nseq=(CL->S)->nseq;
3220   score=declare_int (nseq,2);
3221   dm=(CL->DM)->score_similarity_matrix;
3222   for (a=0; a<nseq; a++)
3223     {
3224       score[a][0]=a;
3225       for ( b=0; b<nseq; b++)
3226         score[a][1]+=dm[a][b];
3227       score[a][1]/=nseq;
3228     }
3229   sort_int ( score,2, 1, 0, nseq-1);
3230
3231   max=20;
3232   for (a=0; a<max; a++)
3233     {
3234       ns[0]=nseq-1;
3235       for (ns[0]=0,b=0; b<nseq; b++)
3236         if (b!=score[a][0])ls[0][ns[0]++]=b;
3237
3238       fprintf (stderr, "[%s %s %d]",(CL->S)->name[score[a][0]],A->name[score[a][0]], score[a][1]);
3239       ns[1]=1;
3240       ls[1][0]=score[a][0];
3241       ungap_sub_aln ( A, ns[0], ls[0]);
3242       ungap_sub_aln ( A, ns[1], ls[1]);
3243       A->score_aln=pair_wise (A, ns, ls,CL);
3244       ls[0][ns[0]++]=a;
3245     }
3246
3247   return A;
3248 }
3249 Alignment *simple_progressive_aln (Sequence *S, NT_node **T, Constraint_list *CL, char *mat)
3250 {
3251   int a;
3252   Alignment *A;
3253
3254
3255   A=seq2aln (S, NULL, RM_GAP);
3256
3257   if ( !CL)
3258     {
3259
3260       CL=declare_constraint_list (S, NULL, NULL, 0, NULL, NULL);
3261       sprintf ( CL->dp_mode,   "myers_miller_pair_wise");
3262       sprintf ( CL->tree_mode, "nj");
3263       sprintf ( CL->distance_matrix_mode, "idscore");
3264       CL=choose_extension_mode ("matrix", CL);
3265       CL->gop=-10;
3266       CL->gep=-1;
3267       if (mat)CL->M=read_matrice (mat);
3268       CL->pw_parameters_set=1;
3269       CL->local_stderr=stderr;
3270     }
3271
3272   if ( !T)T=make_tree (A, CL, CL->gop, CL->gep,S, NULL,MAXIMISE);
3273   for ( a=0; a< A->nseq; a++)ungap (A->seq_al[a]);
3274
3275   tree_aln ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
3276   A=reorder_aln ( A,A->tree_order,A->nseq);
3277
3278   return A;
3279 }
3280
3281 Alignment *very_fast_aln ( Alignment*A, int nseq, Constraint_list *CL)
3282 {
3283 char command[10000];
3284 char *tmp_seq;
3285 char *tmp_aln;
3286 FILE *fp;
3287
3288 if ( CL && CL->local_stderr)fp=CL->local_stderr;
3289 else fp=stderr;
3290
3291  fprintf (fp, "\n[Computation of an Approximate MSA...");
3292  tmp_seq= vtmpnam (NULL);
3293  tmp_aln= vtmpnam (NULL);
3294  output_fasta_seq ((tmp_seq=vtmpnam (NULL)), A);
3295  sprintf ( command, "t_coffee -infile=%s -special_mode quickaln -outfile=%s %s -outorder=input", tmp_seq, tmp_aln, TO_NULL_DEVICE);
3296  my_system ( command);
3297  A->nseq=0;
3298  A=main_read_aln (tmp_aln,A);
3299  fprintf (fp, "]\n");
3300  return A;
3301 }
3302
3303 static NT_node* SNL;
3304 NT_node* tree_aln ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3305 {
3306   if ( strm ((CL->TC)->use_seqan, "NO"))return local_tree_aln (LT, RT, A, nseq, CL);
3307   else return seqan_tree_aln (LT, RT, A, nseq, CL);
3308 }
3309
3310 NT_node* seqan_tree_aln ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3311   {
3312
3313
3314     Alignment *B;
3315
3316
3317     char *tree, *lib, *seq, *new_aln;
3318
3319
3320     //Output tree
3321     tree=vtmpnam (NULL);
3322     print_newick_tree (LT->parent, tree);
3323
3324
3325     //Output seq
3326     main_output_fasta_seq (seq=vtmpnam (NULL),B=seq2aln (CL->S,NULL,RM_GAP), NO_HEADER);
3327     free_aln (B);
3328
3329     //Output lib
3330     new_aln=vtmpnam (NULL);
3331     vfclose (save_constraint_list ( CL, 0, CL->ne,lib=vtmpnam(NULL), NULL, "ascii",CL->S));
3332
3333     fprintf (CL->local_stderr, "\n********* USE EXTERNAL ALIGNER: START:\n\tCOMMAND: %s -lib %s -seq %s -usetree %s -outfile %s\n", (CL->TC)->use_seqan,lib, seq, tree, new_aln);
3334     printf_system ( "%s -lib %s -seq %s -usetree %s -outfile %s", (CL->TC)->use_seqan,lib, seq, tree, new_aln);
3335     fprintf (CL->local_stderr, "\n********* USE EXTERNAL ALIGNER: END\n");
3336
3337
3338     main_read_aln (new_aln, A);
3339     return tree2ao (LT,RT, A, A->nseq, CL);
3340
3341
3342
3343   }
3344 NT_node rec_local_tree_aln ( NT_node P, Alignment*A, Constraint_list *CL, int print);
3345 NT_node* local_tree_aln ( NT_node l, NT_node r, Alignment*A,int nseq, Constraint_list *CL)
3346 {
3347   int a;
3348   NT_node P, *NL;
3349   int **min=NULL;
3350
3351   if (!r && !l) return NULL;
3352   else if (!r)P=l;
3353   else if (!l)P=r;
3354   else P=r->parent;
3355
3356   fprintf ( CL->local_stderr, "\nPROGRESSIVE_ALIGNMENT [Tree Based]\n");
3357
3358   //1: make sure the Alignment and the Sequences are labeled the same way
3359   if (CL->translation)vfree (CL->translation);
3360   CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
3361   for ( a=0; a< (CL->S)->nseq; a++)
3362     CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
3363   A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
3364   A->nseq=(CL->S)->nseq;
3365
3366   //2 Make sure the tree is in the same order
3367   recode_tree (P, (CL->S));
3368   index_tree_node(P);
3369
3370   if ( get_nproc()>1 && strstr (CL->multi_thread, "msa"))
3371     {
3372       int max_fork;
3373       max_fork=get_nproc()/2;//number of nodes forked, one node =>two jobs
3374       tree2nnode (P);
3375       NL=tree2node_list (P, NULL);
3376       min=declare_int (P->node+1,3);
3377       for (a=0; a<=P->node; a++)
3378         {
3379           NT_node N;
3380           N=NL[a];
3381           min[a][0]=a;
3382           if (!N);
3383           else if (N && N->nseq==1)min[a][1]=0;
3384           else
3385             {
3386               min[a][1]=MIN(((N->left)->nseq),((N->right)->nseq))*A->nseq+MAX(((N->left)->nseq),((N->right)->nseq));//sort on min and break ties on max
3387               min[a][2]=MIN(((N->left)->nseq),((N->right)->nseq));
3388             }
3389         }
3390       sort_int_inv (min,3, 1, 0, P->node);
3391       for (a=0; a<=P->node && a<max_fork; a++)
3392         {
3393           if (min[a][2]>1)(NL[min[a][0]])->fork=1;
3394         }
3395     }
3396   free_int (min, -1);
3397   rec_local_tree_aln (P, A,CL, 1);
3398   for (a=0; a<P->nseq; a++)sprintf (A->tree_order[a], "%s", (CL->S)->name[P->lseq[a]]);
3399   A->len_aln=strlen (A->seq_al[0]);
3400
3401   fprintf ( CL->local_stderr, "\n\n");
3402
3403   return NULL;
3404 }
3405
3406 NT_node rec_local_tree_aln ( NT_node P, Alignment*A, Constraint_list *CL,int print)
3407 {
3408   NT_node R,L;
3409   int score;
3410
3411
3412   if (!P || P->nseq==1) return NULL;
3413   R=P->right;L=P->left;
3414
3415   if (P->fork )
3416     {
3417       int s, pid1, pid2;
3418       char *tmp1, *tmp2;
3419       tmp1=vtmpnam (NULL);
3420       tmp2=vtmpnam (NULL);
3421
3422       pid1=vfork();
3423       if (pid1==0)
3424         {
3425           if (print==1)
3426             if (L->nseq>R->nseq)print=1;
3427
3428           initiate_vtmpnam (NULL);
3429           rec_local_tree_aln (L, A, CL, print);
3430           dump_msa (tmp1,A, L->nseq, L->lseq);
3431           exit (EXIT_SUCCESS);
3432         }
3433       else
3434         {
3435           pid2=vfork();
3436           if (pid2==0)
3437             {
3438               if (print==1)
3439                 if (L->nseq>R->nseq)print=0;
3440
3441
3442               initiate_vtmpnam (NULL);
3443               rec_local_tree_aln (R, A, CL, print);
3444               dump_msa (tmp2, A, R->nseq, R->lseq);
3445               exit (EXIT_SUCCESS);
3446             }
3447         }
3448       vwaitpid (pid1, &s, 0);
3449       vwaitpid (pid2, &s, 0);
3450
3451
3452       undump_msa (A,tmp1);
3453       undump_msa (A,tmp2);
3454     }
3455   else
3456     {
3457       rec_local_tree_aln (L, A, CL, print);
3458       rec_local_tree_aln (R, A, CL, print);
3459     }
3460
3461   P->score=A->score_aln=score=profile_pair_wise (A,L->nseq, L->lseq,R->nseq,R->lseq,CL);
3462   A->len_aln=strlen (A->seq_al[P->lseq[0]]);
3463   score=node2sub_aln_score (A, CL, CL->evaluate_mode,P);
3464   if (print)fprintf(CL->local_stderr, "\n\tGroup %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->[Score=%4d][Len=%5d][PID:%d]%s",P->index,R->index,R->nseq,L->index,L->nseq,score, A->len_aln,getpid(),(P->fork==1)?"[Forked]":"" );
3465
3466   return P;
3467 }
3468
3469
3470
3471 NT_node* tree2ao ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3472     {
3473     int *n_s;
3474     int ** l_s;
3475     int a, b;
3476
3477     static int n_groups_done, do_split=0;
3478     int  nseq2align=0;
3479     int *translation;
3480
3481
3482     NT_node P=NULL;
3483
3484
3485
3486
3487     if (n_groups_done==0)
3488        {
3489          if (SNL)vfree(SNL);
3490          SNL=vcalloc ( (CL->S)->nseq, sizeof (NT_node));
3491
3492          if (CL->translation)vfree(CL->translation);
3493          CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
3494
3495          for ( a=0; a< (CL->S)->nseq; a++)
3496            CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
3497
3498          n_groups_done=(CL->S)->nseq;
3499          A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
3500          A->nseq=nseq;
3501        }
3502
3503     translation=CL->translation;
3504     n_s=vcalloc (2, sizeof ( int));
3505     l_s=declare_int ( 2, nseq);
3506
3507
3508     if ( RT->parent !=LT->parent)fprintf ( stderr, "Tree Pb [FATAL:%s]", PROGRAM);
3509     else P=RT->parent;
3510
3511     if ( LT->leaf==1 && RT->leaf==0)
3512       tree2ao ( RT->left, RT->right,A, nseq,CL);
3513
3514     else if ( RT->leaf==1 && LT->leaf==0)
3515       tree2ao ( LT->left, LT->right,A,nseq,CL);
3516
3517     else if (RT->leaf==0 && LT->leaf==0)
3518       {
3519         tree2ao ( LT->left, LT->right,A,nseq,CL);
3520         tree2ao ( RT->left, RT->right,A,nseq,CL);
3521       }
3522
3523     if ( LT->leaf==1 && RT->leaf==1)
3524       {
3525         /*1 Identify the two groups of sequences to align*/
3526
3527         nseq2align=LT->nseq+RT->nseq;
3528         n_s[0]=LT->nseq;
3529         for ( a=0; a< LT->nseq; a++)l_s[0][a]=translation[LT->lseq[a]];
3530         if ( LT->nseq==1)LT->group=l_s[0][0];
3531
3532         n_s[1]=RT->nseq;
3533         for ( a=0; a< RT->nseq; a++)l_s[1][a]=translation[RT->lseq[a]];
3534         if ( RT->nseq==1)RT->group=l_s[1][0];
3535
3536
3537         P->group=n_groups_done++;
3538
3539         if (nseq2align==nseq)
3540           {
3541             for (b=0, a=0; a< n_s[0]; a++, b++)sprintf ( A->tree_order[b],"%s", (CL->S)->name[l_s[0][a]]);
3542             for (a=0; a< n_s[1]     ; a++, b++)sprintf ( A->tree_order[b], "%s",(CL->S)->name[l_s[1][a]]);
3543             n_groups_done=0;
3544           }
3545       }
3546     if (P->parent)P->leaf=1;
3547     if ( LT->isseq==0)LT->leaf=0;
3548     if ( RT->isseq==0)RT->leaf=0;
3549
3550     if (RT->isseq){SNL[translation[RT->lseq[0]]]=RT;RT->score=100;}
3551     if (LT->isseq){SNL[translation[LT->lseq[0]]]=LT;LT->score=100;}
3552
3553     do_split=split_condition (nseq2align,A->score_aln,CL);
3554     if (CL->split && do_split)
3555       {
3556
3557         for (a=0; a< P->nseq; a++)SNL[CL->translation[P->lseq[a]]]=NULL;
3558         SNL[CL->translation[RT->lseq[0]]]=P;
3559
3560       }
3561
3562     vfree ( n_s);
3563     free_int ( l_s, 2);
3564     return SNL;
3565
3566     }
3567
3568 NT_node* tree_realn ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3569     {
3570     int *n_s;
3571     int ** l_s;
3572     int a, b;
3573     int score;
3574     static int n_groups_done;
3575     int nseq2align=0;
3576     int *translation;
3577
3578
3579     NT_node P=NULL;
3580
3581
3582
3583
3584     if (n_groups_done==0)
3585        {
3586          if (SNL)vfree(SNL);
3587          SNL=vcalloc ( (CL->S)->nseq, sizeof (NT_node));
3588
3589          if (CL->translation)vfree(CL->translation);
3590          CL->translation=vcalloc ( (CL->S)->nseq, sizeof (int));
3591
3592          for ( a=0; a< (CL->S)->nseq; a++)
3593            CL->translation[a]=name_is_in_list ( (CL->S)->name[a], (CL->S)->name, (CL->S)->nseq, MAXNAMES);
3594          if (nseq>2)fprintf ( CL->local_stderr, "\nPROGRESSIVE_ALIGNMENT [Tree Based]\n");
3595          else fprintf ( CL->local_stderr, "\nPAIRWISE_ALIGNMENT [No Tree]\n");
3596          n_groups_done=(CL->S)->nseq;
3597          A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
3598          A->nseq=nseq;
3599        }
3600
3601     translation=CL->translation;
3602     n_s=vcalloc (2, sizeof ( int));
3603     l_s=declare_int ( 2, nseq);
3604
3605
3606     if ( nseq==2)
3607        {
3608        n_s[0]=n_s[1]=1;
3609        l_s[0][0]=name_is_in_list ((CL->S)->name[0],(CL->S)->name, (CL->S)->nseq, MAXNAMES);
3610        l_s[1][0]=name_is_in_list ((CL->S)->name[1],(CL->S)->name, (CL->S)->nseq, MAXNAMES);
3611        A->score_aln=score=pair_wise (A, n_s, l_s,CL);
3612
3613        vfree ( n_s);
3614        free_int ( l_s, 2);
3615        return SNL;
3616        }
3617     else
3618        {
3619        if ( RT->parent !=LT->parent)fprintf ( stderr, "Tree Pb [FATAL:%s]", PROGRAM);
3620        else P=RT->parent;
3621
3622        if ( LT->leaf==1 && RT->leaf==0)
3623            tree_realn ( RT->left, RT->right,A, nseq,CL);
3624
3625        else if ( RT->leaf==1 && LT->leaf==0)
3626            tree_realn ( LT->left, LT->right,A,nseq,CL);
3627
3628        else if (RT->leaf==0 && LT->leaf==0)
3629           {
3630           tree_realn ( LT->left, LT->right,A,nseq,CL);
3631           tree_realn ( RT->left, RT->right,A,nseq,CL);
3632           }
3633
3634        if ( LT->leaf==1 && RT->leaf==1 && (RT->nseq+LT->nseq)<nseq)
3635           {
3636           /*1 Identify the two groups of sequences to align*/
3637             int *list, s, id1, id2;
3638             list=vcalloc (nseq, sizeof (int));
3639             for (a=0; a<LT->nseq; a++)
3640               {
3641                 s=translation[LT->lseq[a]];
3642                 list[s]=1;
3643               }
3644             for (a=0; a<RT->nseq; a++)
3645               {
3646                 s=translation[RT->lseq[a]];
3647                 list[s]=1;
3648               }
3649             for (a=0; a<nseq; a++)
3650               {
3651                 s=list[a];
3652                 l_s[s][n_s[s]++]=a;
3653               }
3654
3655             vfree (list);
3656
3657             id1=sub_aln2sim (A, n_s, l_s, "idmat_sim");
3658
3659
3660             ungap_sub_aln (A, n_s[0],l_s[0]);
3661             ungap_sub_aln (A, n_s[1],l_s[1]);
3662             P->score=A->score_aln=score=pair_wise (A, n_s, l_s,CL);
3663             id2=sub_aln2sim (A, n_s, l_s, "idmat_sim");
3664
3665
3666
3667
3668             if (nseq2align==nseq)
3669               {
3670                 for (b=0, a=0; a< n_s[0]; a++, b++)sprintf ( A->tree_order[b],"%s", (CL->S)->name[l_s[0][a]]);
3671                 for (a=0; a< n_s[1]     ; a++, b++)sprintf ( A->tree_order[b], "%s",(CL->S)->name[l_s[1][a]]);
3672                 n_groups_done=0;
3673               }
3674           }
3675        if (P->parent)P->leaf=1;
3676        //Recycle the tree
3677        if ( LT->isseq==0)LT->leaf=0;
3678        if ( RT->isseq==0)RT->leaf=0;
3679
3680        if (RT->isseq){SNL[translation[RT->lseq[0]]]=RT;RT->score=100;}
3681        if (LT->isseq){SNL[translation[LT->lseq[0]]]=LT;LT->score=100;}
3682
3683        vfree ( n_s);
3684        free_int ( l_s, 2);
3685        return SNL;
3686        }
3687
3688
3689     }
3690
3691
3692
3693 Alignment* profile_tree_aln ( NT_node P,Alignment*A,Constraint_list *CL, int threshold)
3694 {
3695   int *ns, **ls, a, sim;
3696   NT_node LT, RT, D, UD;
3697   Alignment *F;
3698   static NT_node R;
3699   static int n_groups_done;
3700
3701
3702   //first pass
3703   //Sequences must be in the same order as the tree sequences
3704    if (!P->parent)
3705     {
3706       R=P;
3707       n_groups_done=P->nseq+1;
3708     }
3709
3710   LT=P->left;
3711   RT=P->right;
3712
3713   if (LT->leaf==0)A=delayed_tree_aln1 (LT, A,CL, threshold);
3714   if (RT->leaf==0)A=delayed_tree_aln1 (RT, A,CL, threshold);
3715
3716   ns=vcalloc (2, sizeof (int));
3717   ls=declare_int ( 2,R->nseq);
3718
3719   if ( LT->nseq==1)
3720     {
3721       ls[0][ns[0]++]=LT->lseq[0];
3722       LT->group=ls[0][0]+1;
3723     }
3724   else
3725     node2seq_list (LT,&ns[0], ls[0]);
3726
3727   if ( RT->nseq==1)
3728     {
3729       ls[1][ns[1]++]=RT->lseq[0];
3730       RT->group=ls[1][0]+1;
3731     }
3732   else
3733     node2seq_list (RT,&ns[1], ls[1]);
3734
3735
3736   P->group=++n_groups_done;
3737   fprintf (CL->local_stderr, "\n\tGroup %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,RT->group, ns[1],LT->group, ns[0]);
3738
3739   P->score=A->score_aln=pair_wise (A, ns, ls,CL);
3740   sim=sub_aln2sim(A, ns, ls, "idmat_sim1");
3741
3742   if ( sim<threshold)
3743     {
3744       UD=(ns[0]<=ns[1])?RT:LT;
3745       D= (ns[0]<=ns[1])?LT:RT;
3746
3747       UD->aligned=1;
3748       D->aligned=0;
3749
3750       fprintf (CL->local_stderr,  "[Delayed (Sim=%4d). Kept Group %4d]",sim,UD->group);
3751
3752
3753       ungap_sub_aln (A, ns[0],ls[0]);
3754       ungap_sub_aln (A, ns[1],ls[1]);
3755       A->nseq=MAX(ns[0],ns[1]);
3756
3757       F=A;
3758       while (F->A)F=F->A;
3759       F->A=main_read_aln (output_fasta_sub_aln (NULL, A, ns[(D==LT)?0:1], ls[(D==LT)?0:1]), NULL);
3760       if ( P==R)
3761         {
3762           F=F->A;
3763           F->A=main_read_aln (output_fasta_sub_aln (NULL, A, ns[(D==LT)?1:0], ls[(D==LT)?1:0]), NULL);
3764         }
3765       if (F->A==NULL)
3766         {
3767           printf_exit (EXIT_FAILURE, stderr, "\nError: Empty group");
3768         }
3769     }
3770   else
3771     {
3772       LT->aligned=1; RT->aligned=1;
3773       fprintf (CL->local_stderr, "[Score=%4d][Len=%5d]",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), strlen ( A->seq_al[ls[0][0]]));
3774       A->nseq=ns[0]+ns[1];
3775       if (P==R)
3776         {
3777           F=A;
3778           while (F->A)F=F->A;
3779           F->A=main_read_aln (output_fasta_sub_aln2 (NULL, A, ns, ls), NULL);
3780         }
3781     }
3782   P->nseq=0;
3783   for (a=0; a<LT->nseq;a++)P->lseq[P->nseq++]=LT->lseq[a];
3784   for (a=0; a<RT->nseq;a++)P->lseq[P->nseq++]=RT->lseq[a];
3785
3786   P->aligned=1;
3787
3788   vfree ( ns);
3789   free_int ( ls,-1);
3790   return A;
3791 }////////////////////////////////////////////////////////////////////////////////////////
3792 //
3793 //                               Frame Tree Aln
3794 //
3795 ////////////////////////////////////////////////////////////////////////////////////////
3796
3797 //Alignment *frame_tree_aln (Alignment *A, Constraint_list *CL)
3798 //{
3799
3800
3801 ////////////////////////////////////////////////////////////////////////////////////////
3802 //
3803 //                               Delayed Tree Aln
3804 //
3805 ////////////////////////////////////////////////////////////////////////////////////////
3806 int delayed_pair_wise (Alignment *A, int *ns, int **ls,Constraint_list *CL);
3807 NT_node* delayed_tree_aln_mode1 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL);
3808 NT_node* delayed_tree_aln_mode2 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL);
3809 int paint_nodes2aligned ( NT_node P,char **list, int n);
3810
3811 int reset_visited_nodes ( NT_node P);
3812 int reset_visited_nodes2 ( NT_node P);
3813 Alignment * make_delayed_tree_aln (Alignment *A,int n, Constraint_list *CL)
3814 {
3815   NT_node **T=NULL;
3816   int a;
3817
3818   T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
3819   delayed_tree_aln_mode1 ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
3820
3821   for ( a=0; a< n; a++)
3822     {
3823
3824       sprintf ( CL->distance_matrix_mode, "aln");
3825       CL->DM=cl2distance_matrix ( CL,A,NULL,NULL, 1);
3826       degap_aln (A);
3827       T=make_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
3828       delayed_tree_aln_mode1 ((T[3][0])->left,(T[3][0])->right,A,(CL->S)->nseq, CL);
3829     }
3830
3831   return A;
3832 }
3833 NT_node* delayed_tree_aln_mode1 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3834 {
3835   NT_node P;
3836
3837
3838
3839   P=LT->parent;P->nseq=nseq;
3840   paint_nodes2aligned (P, NULL, 0);
3841
3842   A=delayed_tree_aln1 (P, A, CL,50);
3843   A=delayed_tree_aln2 (P, A, CL, 0);
3844   return NULL;
3845 }
3846
3847 NT_node* delayed_tree_aln_mode2 ( NT_node LT, NT_node RT, Alignment*A, int nseq, Constraint_list *CL)
3848 {
3849   NT_node P;
3850
3851   int thr=50;
3852
3853   P=LT->parent;P->nseq=nseq;
3854
3855   A=delayed_tree_aln1 (P, A, CL,thr);
3856   thr-=10;
3857   while (thr>=0)
3858     {
3859       A=delayed_tree_aln2 (P, A, CL, thr);
3860       thr-=10;
3861     }
3862   return NULL;
3863 }
3864
3865 Alignment* delayed_tree_aln1 ( NT_node P,Alignment*A,Constraint_list *CL, int threshold)
3866 {
3867   int *ns, **ls, a, sim;
3868   NT_node LT, RT, D, UD;
3869
3870   static NT_node R;
3871   static int n_groups_done;
3872
3873
3874   //first pass
3875   //Sequences must be in the same order as the tree sequences
3876    if (!P->parent)
3877     {
3878       R=P;
3879       n_groups_done=P->nseq+1;
3880     }
3881
3882   LT=P->left;
3883   RT=P->right;
3884
3885   if (LT->leaf==0)A=delayed_tree_aln1 (LT, A,CL, threshold);
3886   if (RT->leaf==0)A=delayed_tree_aln1 (RT, A,CL, threshold);
3887
3888   ns=vcalloc (2, sizeof (int));
3889   ls=declare_int ( 2,R->nseq);
3890
3891
3892   node2seq_list (LT,&ns[0], ls[0]);
3893   if ( LT->nseq==1)LT->group=LT->lseq[0]+1;
3894
3895   node2seq_list (RT,&ns[1], ls[1]);
3896   if ( RT->nseq==1)RT->group=RT->lseq[0]+1;
3897
3898
3899   P->group=++n_groups_done;
3900
3901
3902   if ( ns[0]==0 || ns[1]==0)
3903     {
3904       fprintf (CL->local_stderr, "\n\tF-Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->Skipped",P->group,RT->group, ns[1],LT->group, ns[0]);
3905
3906       LT->aligned=(ns[0]==0)?0:1;
3907       RT->aligned=(ns[1]==0)?0:1;
3908     }
3909   else
3910     {
3911       fprintf (CL->local_stderr, "\n\tF-Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,RT->group, ns[1],LT->group, ns[0]);
3912       P->score=A->score_aln=pair_wise (A, ns, ls,CL);
3913       sim=sub_aln2max_sim(A, ns, ls, "idmat_sim1");
3914
3915
3916       if ( sim<threshold)
3917         {
3918           UD=(ns[0]<=ns[1])?RT:LT;
3919           D= (ns[0]<=ns[1])?LT:RT;
3920
3921           UD->aligned=1;
3922           D->aligned=0;
3923
3924           fprintf (CL->local_stderr,  "[Delayed (Sim=%4d). Kept Group %4d]",sim,UD->group);
3925
3926           ungap_sub_aln (A, ns[0],ls[0]);
3927           ungap_sub_aln (A, ns[1],ls[1]);
3928           A->nseq=MAX(ns[0],ns[1]);
3929         }
3930       else
3931         {
3932           LT->aligned=1; RT->aligned=1;
3933           fprintf (CL->local_stderr, "[Score=%4d][Len=%5d]",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), strlen ( A->seq_al[ls[0][0]]));
3934           A->nseq=ns[0]+ns[1];
3935         }
3936       P->nseq=0;
3937       for (a=0; a<LT->nseq;a++)P->lseq[P->nseq++]=LT->lseq[a];
3938       for (a=0; a<RT->nseq;a++)P->lseq[P->nseq++]=RT->lseq[a];
3939
3940       P->aligned=1;
3941     }
3942   vfree ( ns);
3943   free_int ( ls,-1);
3944   return A;
3945 }
3946
3947 Alignment* delayed_tree_aln2 ( NT_node P,Alignment*A,Constraint_list *CL, int thr)
3948 {
3949
3950   NT_node LT, RT, D;
3951
3952   static NT_node R;
3953
3954
3955   LT=P->left;
3956   RT=P->right;
3957   if (!P->parent)
3958     {
3959       R=P;
3960       fprintf (CL->local_stderr, "\n");
3961     }
3962   if (!LT->aligned && !RT->aligned)
3963     {
3964       printf_exit (EXIT_FAILURE, stderr, "ERROR: Unresolved Node On Groups %d  [FATAL:%s]\n", P->group,PROGRAM);
3965     }
3966   else if (!LT->aligned || !RT->aligned)
3967     {
3968       int *ns, **ls, sim;
3969       ns=vcalloc (2, sizeof (int));
3970       ls=declare_int (2, R->nseq);
3971
3972       node2seq_list (R,&ns[0], ls[0]);
3973
3974       D=(!LT->aligned)?LT:RT;
3975       D->aligned=1;
3976       node2seq_list (D,&ns[1], ls[1]);
3977
3978       fprintf (CL->local_stderr, "\tS-Delayed Group %4d: [Group %4d (%4d seq)] with [Group %4d (%4d seq)]-->",P->group,D->group, ns[1],R->group, ns[0]);
3979       P->score=A->score_aln=pair_wise (A, ns, ls,CL);
3980       sim=sub_aln2max_sim(A, ns, ls, "idmat_sim1");
3981       if (sim<thr)
3982         {
3983           fprintf (CL->local_stderr, " [Further Delayed]\n");
3984           ungap_sub_aln (A, ns[0],ls[0]);
3985           ungap_sub_aln (A, ns[1],ls[1]);
3986           D->aligned=0;
3987         }
3988       else
3989         {
3990           fprintf (CL->local_stderr, "[Score=%4d][Len=%5d][thr=%d]\n",sub_aln2sub_aln_score (A, CL, CL->evaluate_mode,ns, ls), (int)strlen ( A->seq_al[ls[0][0]]), thr);
3991           D->aligned=1;
3992         }
3993       vfree (ns);free_int (ls, -1);
3994     }
3995   else
3996     {
3997       ;
3998     }
3999
4000   if (LT->leaf==0)A=delayed_tree_aln2 (LT, A,CL, thr);
4001   if (RT->leaf==0)A=delayed_tree_aln2 (RT, A,CL, thr);
4002
4003   return A;
4004 }
4005
4006 int delayed_pair_wise (Alignment *A, int *ns, int **ls,Constraint_list *CL)
4007 {
4008   int s,s1, s2, a, b;
4009   int **sim;
4010
4011
4012   pair_wise (A, ns, ls, CL);
4013
4014   sim=fast_aln2sim_list (A, "sim3", ns, ls);
4015
4016   sort_int_inv ( sim,3, 2,0, ns[0]*ns[1]-1);
4017
4018   for (a=0; a< 2; a++)
4019     for ( b=0; b< ns[a]; b++)
4020       A->order[ls[a][b]][4]=-1;
4021
4022   for (a=0; a< 10 && sim[a][0]!=-1; a++)
4023     {
4024       s1=sim[a][0];
4025       s2=sim[a][1];
4026       A->order[s1][4]=0;
4027       A->order[s2][4]=0;
4028     }
4029
4030   ungap_sub_aln (A, ns[0],ls[0]);
4031   ungap_sub_aln (A, ns[1],ls[1]);
4032
4033   s=pair_wise (A, ns, ls, CL);
4034
4035   for (a=0; a< 2; a++)
4036     for ( b=0; b< ns[a]; b++)
4037       A->order[ls[a][b]][4]=0;
4038
4039   free_int (sim, -1);
4040   return s;
4041 }
4042
4043 int node2seq_list2 (NT_node P, int *ns, int *ls)
4044 {
4045
4046   if ( !P || P->visited ) return ns[0];
4047   else P->visited=1;
4048
4049   if ( P->isseq)
4050     {
4051       ls[ns[0]++]=P->lseq[0];
4052     }
4053
4054   if (P->left   && (P->left) ->aligned)node2seq_list2 (P->left, ns,ls);
4055   if (P->right  && (P->right)->aligned)node2seq_list2 (P->right,ns,ls);
4056   if (P->aligned && P->parent)node2seq_list2 (P->parent,ns,ls);
4057
4058
4059   return ns[0];
4060 }
4061
4062 int node2seq_list (NT_node P, int *ns, int *ls)
4063 {
4064
4065   if ( P->isseq && P->aligned)
4066     {
4067       ls[ns[0]++]=P->lseq[0];
4068     }
4069   else
4070     {
4071       if (P->left &&  (P->left) ->aligned)node2seq_list (P->left, ns,ls);
4072       if (P->right && (P->right)->aligned)node2seq_list (P->right,ns,ls);
4073     }
4074   return ns[0];
4075 }
4076 int paint_nodes2aligned ( NT_node P,char **list, int n)
4077 {
4078   int r=0;
4079   if ( P->leaf)
4080     {
4081       if ( list==NULL)
4082         P->aligned=1;
4083       else if ( name_is_in_list ( P->name, list, n, 100)!=-1)
4084         P->aligned=1;
4085       else
4086         P->aligned=0;
4087       return P->aligned;
4088     }
4089   else
4090     {
4091       r+=paint_nodes2aligned (P->left, list, n);
4092       r+=paint_nodes2aligned (P->right, list, n);
4093     }
4094   return r;
4095 }
4096
4097 int reset_visited_nodes ( NT_node P)
4098 {
4099   while (P->parent)P=P->parent;
4100   return reset_visited_nodes2 (P);
4101 }
4102 int reset_visited_nodes2 ( NT_node P)
4103 {
4104   int r=0;
4105   if (P->left)r+=reset_visited_nodes2(P->left);
4106   if (P->right)r+=reset_visited_nodes2(P->right);
4107   r+=P->visited;
4108   P->visited=0;
4109   return r;
4110 }
4111
4112 ////////////////////////////////////////////////////////////////////////////////////////
4113 //
4114 //                               DPA_MSA
4115 //
4116 ////////////////////////////////////////////////////////////////////////////////////////
4117
4118 Alignment* dpa_msa2 ( NT_node P,Alignment*A,Constraint_list *CL);
4119 Alignment *dpa_align_node (NT_node P,Alignment*A,Constraint_list *CL);
4120 char *node2profile_list (NT_node P,Alignment*A,Constraint_list *CL, char *list);
4121 char * output_node_aln (NT_node P, Alignment *A, char *name);
4122 int node2nleaf ( NT_node P);
4123
4124 Alignment* dpa_aln (Alignment*A,Constraint_list *CL)
4125 {
4126   NT_node P;
4127
4128   A=iterative_tree_aln (A,1, CL);
4129   P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
4130   degap_aln (A);
4131   while (!P->leaf)
4132     A=dpa_msa2(P, A, CL);
4133   return A;
4134 }
4135
4136 int node2nleaf ( NT_node P)
4137 {
4138   int n=0;
4139   if ( P->leaf) return 1;
4140   else
4141     {
4142       n+=node2nleaf ( P->left);
4143       n+=node2nleaf ( P->right);
4144     }
4145   return n;
4146 }
4147 Alignment* dpa_msa2 ( NT_node P,Alignment*A,Constraint_list *CL)
4148 {
4149   int maxnseq=20;
4150   int n, n_l, n_r;
4151   n=node2nleaf (P);
4152
4153
4154   if ( n>maxnseq)
4155     {
4156         n_l=node2nleaf (P->left);
4157         n_r=node2nleaf (P->right);
4158         if (n_l>n_r)
4159           {
4160             return dpa_msa2 (P->left, A, CL);
4161           }
4162         else
4163           {
4164             return dpa_msa2 (P->right, A, CL);
4165           }
4166     }
4167   A=dpa_align_node (P, A, CL);
4168   P->leaf=1;
4169   return A;
4170 }
4171
4172 Alignment *dpa_align_node (NT_node P,Alignment*A,Constraint_list *CL)
4173 {
4174
4175   char *list, *tmp_aln;
4176   int a, b;
4177   Alignment *B;
4178
4179
4180   list=vcalloc ( P->nseq*100, sizeof (char));
4181   list=node2profile_list (P,A, CL, list);
4182
4183   printf_system ( "t_coffee -profile %s -outfile=%s -dp_mode gotoh_pair_wise_lgp -msa_mode iterative_tree_aln -quiet", list,tmp_aln=vtmpnam (NULL));
4184   B=main_read_aln (tmp_aln, NULL);
4185   A=realloc_aln (A, B->len_aln+1);
4186   for ( a=0; a< B->nseq; a++)
4187     if ( (b=name_is_in_list (B->name[a], A->name, A->nseq, 100))!=-1)
4188       sprintf (A->seq_al[b], "%s", B->seq_al[a]);
4189   A->len_aln=B->len_aln;
4190   free_aln (B);
4191   vfree (list);
4192   return A;
4193 }
4194 char *node2profile_list (NT_node P,Alignment*A,Constraint_list *CL, char *list)
4195 {
4196   if (!P->leaf)
4197     {
4198       list=node2profile_list (P->left, A, CL, list);
4199       list=node2profile_list (P->right, A, CL, list);
4200     }
4201   else
4202     {
4203
4204       list=strcatf (list," %s", output_node_aln (P, A, NULL));
4205       if ( !P->isseq)P->leaf=0;
4206     }
4207   return list;
4208 }
4209 char * output_node_aln (NT_node P, Alignment *A, char *name)
4210 {
4211   FILE *fp;
4212   int a;
4213   if (name==NULL) name=vtmpnam (NULL);
4214   fp=vfopen (name, "w");
4215
4216   for (a=0; a< P->nseq; a++)
4217     fprintf ( fp, ">%s\n%s", A->name[P->lseq[a]], A->seq_al[P->lseq[a]]);
4218   vfclose (fp);
4219   return name;
4220 }
4221 ////////////////////////////////////////////////////////////////////////////////////////
4222 //
4223 //                               NEW_DPA_MSA
4224 //
4225 ////////////////////////////////////////////////////////////////////////////////////////
4226
4227 Alignment * new_dpa_aln (Alignment *A,Constraint_list *CL)
4228 {
4229   NT_node P;
4230
4231   char *tmp_aln;
4232   char *list;
4233
4234   A=make_delayed_tree_aln (A,1, CL);
4235   P=make_root_tree (A, CL, CL->gop, CL->gep,CL->S,NULL, 1);
4236   set_node_score (A, P, "idmat_sim");
4237
4238
4239   list=split_nodes_nseq (A,P,15, list=vcalloc (P->nseq*200, sizeof (char)));
4240   printf_system ( "t_coffee -profile %s -outfile=%s -dp_mode gotoh_pair_wise_lgp -msa_mode iterative_tree_aln", list,tmp_aln=vtmpnam (NULL));
4241   return main_read_aln (tmp_aln, NULL);
4242 }
4243
4244 char *split_nodes_nseq  (Alignment *A, NT_node P, int nseq, char *list)
4245 {
4246   int a,n;
4247
4248   n=P->nseq;
4249   a=100;
4250   while ( n>=nseq)
4251     {
4252       a--;
4253       n=count_threshold_nodes (A, P, a);
4254     }
4255
4256   return split_nodes_idmax (A, P, a,list);
4257 }
4258 char *split_nodes_idmax (Alignment *A, NT_node P, int t, char *list)
4259 {
4260   if (P->isseq || P->score>=t)
4261     {
4262       list=strcatf (list," %s", output_node_aln (P, A, NULL));
4263     }
4264   else if ( P->score<t)
4265     {
4266       list=split_nodes_idmax (A, P->left,t, list);
4267       list=split_nodes_idmax (A, P->right,t,list);
4268     }
4269   return list;
4270 }
4271 int count_threshold_nodes (Alignment *A, NT_node P, int t)
4272 {
4273   int s=0;
4274
4275   if (P->isseq || P->score>=t)
4276     {
4277       s=1;
4278     }
4279   else if ( P->score<t)
4280     {
4281       s+=count_threshold_nodes (A, P->left,t);
4282       s+=count_threshold_nodes (A, P->right,t);
4283     }
4284
4285   return s;
4286 }
4287 int set_node_score (Alignment *A, NT_node P, char *mode)
4288 {
4289   int a;
4290   int ns[2], *ls[2];
4291
4292   if (P->isseq) return 0;
4293   for (a=0; a<2; a++)
4294     {
4295       NT_node N;
4296       N=(a==0)?P->left:P->right;
4297       ns[a]=N->nseq;
4298       ls[a]=N->lseq;
4299     }
4300   P->score=sub_aln2max_sim(A, ns, ls,mode);
4301   set_node_score (A,P->left, mode);
4302   set_node_score (A,P->right, mode);
4303   return 1;
4304 }
4305 /////////////////////////////////////////////////////////////////////////////////////////
4306 int split_condition (int nseq, int score, Constraint_list *CL)
4307 {
4308   int cond1=1, cond2=1;
4309
4310
4311   if ( CL->split_nseq_thres)cond1 =(nseq<=CL->split_nseq_thres)?1:0;
4312   if ( CL->split_score_thres)cond2=(score>=CL->split_score_thres)?1:0;
4313
4314   return (cond1 && cond2);
4315 }
4316 int profile_pair_wise (Alignment *A, int n1, int *l1, int n2, int *l2, Constraint_list *CL)
4317 {
4318   static int *ns;
4319   static int **ls;
4320
4321   if (!ns)
4322     {
4323       ns=vcalloc (2, sizeof (int));
4324       ls=vcalloc (2, sizeof (int*));
4325     }
4326   ns[0]=n1;
4327   ls[0]=l1;
4328   ns[1]=n2;
4329   ls[1]=l2;
4330   return pair_wise (A, ns, ls, CL);
4331 }
4332 int pair_wise (Alignment *A, int*ns, int **l_s,Constraint_list *CL )
4333     {
4334         /*
4335          CL->maximise
4336          CL->gop;
4337          CL->gep
4338          CL->TG_MODE;
4339         */
4340         int score;
4341
4342         int glocal;
4343         Pwfunc function;
4344
4345
4346
4347         /*Make sure evaluation functions update their cache if needed*/
4348         A=update_aln_random_tag (A);
4349
4350         if (! CL->pw_parameters_set)
4351                    {
4352                        fprintf ( stderr, "\nERROR pw_parameters_set must be set in pair_wise [FATAL]\n" );crash("");
4353                    }
4354
4355
4356         function=get_pair_wise_function(CL->pair_wise,  CL->dp_mode,&glocal);
4357         if ( CL->get_dp_cost==NULL)CL->get_dp_cost=get_dp_cost;
4358
4359         if (strlen ( A->seq_al[l_s[0][0]])==0 || strlen ( A->seq_al[l_s[1][0]])==0)
4360           score=empty_pair_wise ( A, ns, l_s, CL, glocal);
4361         else
4362           score=function ( A, ns, l_s, CL);
4363
4364         return score;
4365     }
4366
4367 int empty_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int glocal)
4368 {
4369   int n=0, a, b;
4370   int *l=NULL;
4371   char *string;
4372   int l0, l1, len;
4373
4374   if ( glocal==GLOBAL)
4375     {
4376       l0=strlen (A->seq_al[l_s[0][0]]);
4377       l1=strlen (A->seq_al[l_s[1][0]]);
4378       len=MAX(l1,l0);
4379
4380       if ( len==0)return 0;
4381       else if (l0>l1){n=ns[1];l=l_s[1];}
4382       else if (l0<l1){n=ns[0];l=l_s[0];}
4383       string=generate_null (len);
4384       for ( a=0; a< n; a++)
4385         sprintf ( A->seq_al[l[a]], "%s", string);
4386       A->score=A->score_aln=0;
4387       A->len_aln=len;
4388       vfree ( string);
4389       return 0;
4390     }
4391   else if ( glocal==LALIGN)
4392     {
4393       A->A=declare_aln (A->S);
4394       (A->A)->len_aln=0;
4395       for ( a=0; a< 2; a++)
4396         for ( b=0; b<ns[a]; b++)
4397           A->seq_al[l_s[a][b]][0]='\0';
4398       (A->A)->score_aln=(A->A)->score=0;
4399       return 0;
4400     }
4401   else return 0;
4402 }
4403
4404
4405
4406
4407 Pwfunc get_pair_wise_function (Pwfunc pw,char *dp_mode, int *glocal)
4408   {
4409     /*Returns a function and a mode (Glogal, Local...)*/
4410
4411
4412
4413     int a;
4414     static int npw;
4415     static Pwfunc *pwl;
4416     static char **dpl;
4417     static int *dps;
4418
4419     /*The first time: initialize the list of pairwse functions*/
4420     if ( npw==0)
4421       {
4422         pwl=vcalloc ( 100, sizeof (Pwfunc));
4423         dpl=declare_char (100, 100);
4424         dps=vcalloc ( 100, sizeof (int));
4425
4426         pwl[npw]=fasta_cdna_pair_wise;
4427         sprintf (dpl[npw], "fasta_cdna_pair_wise");
4428         dps[npw]=GLOBAL;
4429         npw++;
4430
4431         pwl[npw]=cfasta_cdna_pair_wise;
4432         sprintf (dpl[npw], "cfasta_cdna_pair_wise");
4433         dps[npw]=GLOBAL;
4434         npw++;
4435
4436         pwl[npw]=idscore_pair_wise;
4437         sprintf (dpl[npw], "idscore_pair_wise");
4438         dps[npw]=GLOBAL;
4439         npw++;
4440
4441         pwl[npw]=gotoh_pair_wise;
4442         sprintf (dpl[npw], "gotoh_pair_wise");
4443         dps[npw]=GLOBAL;
4444         npw++;
4445
4446         pwl[npw]=gotoh_pair_wise_lgp;
4447         sprintf (dpl[npw], "gotoh_pair_wise_lgp");
4448         dps[npw]=GLOBAL;
4449         npw++;
4450
4451         pwl[npw]=gotoh_pair_wise_lgp_sticky;
4452         sprintf (dpl[npw], "gotoh_pair_wise_lgp_sticky");
4453         dps[npw]=GLOBAL;
4454         npw++;
4455
4456         pwl[npw]=proba_pair_wise;
4457         sprintf (dpl[npw], "proba_pair_wise");
4458         dps[npw]=GLOBAL;
4459         npw++;
4460
4461         pwl[npw]=biphasic_pair_wise;
4462         sprintf (dpl[npw], "biphasic_pair_wise");
4463         dps[npw]=GLOBAL;
4464         npw++;
4465
4466         pwl[npw]=subop1_pair_wise;
4467         sprintf (dpl[npw], "subop1_pair_wise");
4468         dps[npw]=GLOBAL;
4469         npw++;
4470
4471         pwl[npw]=subop2_pair_wise;
4472         sprintf (dpl[npw], "subop2_pair_wise");
4473         dps[npw]=GLOBAL;
4474         npw++;
4475
4476         pwl[npw]=myers_miller_pair_wise;
4477         sprintf (dpl[npw], "myers_miller_pair_wise");
4478         dps[npw]=GLOBAL;
4479         npw++;
4480
4481         pwl[npw]=test_pair_wise;
4482         sprintf (dpl[npw], "test_pair_wise");
4483         dps[npw]=GLOBAL;
4484         npw++;
4485
4486         pwl[npw]=fasta_gotoh_pair_wise;
4487         sprintf (dpl[npw], "fasta_pair_wise");
4488         dps[npw]=GLOBAL;
4489         npw++;
4490         pwl[npw]=cfasta_gotoh_pair_wise;
4491         sprintf (dpl[npw], "cfasta_pair_wise");
4492         dps[npw]=GLOBAL;
4493         npw++;
4494
4495         pwl[npw]=very_fast_gotoh_pair_wise;
4496         sprintf (dpl[npw], "very_fast_pair_wise");
4497         dps[npw]=GLOBAL;
4498         npw++;
4499
4500         pwl[npw]=gotoh_pair_wise_sw;
4501         sprintf (dpl[npw], "gotoh_pair_wise_sw");
4502         dps[npw]=LOCAL;
4503         npw++;
4504
4505         pwl[npw]=cfasta_gotoh_pair_wise_sw;
4506         sprintf (dpl[npw], "cfasta_sw_pair_wise");
4507         dps[npw]=LOCAL;
4508         npw++;
4509
4510         pwl[npw]=gotoh_pair_wise_lalign;
4511         sprintf (dpl[npw], "gotoh_pair_wise_lalign");
4512         dps[npw]=LALIGN;
4513         npw++;
4514
4515         pwl[npw]=sim_pair_wise_lalign;
4516         sprintf (dpl[npw], "sim_pair_wise_lalign");
4517         dps[npw]=LALIGN;
4518         npw++;
4519
4520         pwl[npw]=domain_pair_wise;
4521         sprintf (dpl[npw], "domain_pair_wise");
4522         dps[npw]=MOCCA;
4523         npw++;
4524
4525         pwl[npw]=gotoh_pair_wise;
4526         sprintf (dpl[npw], "ssec_pair_wise");
4527         dps[npw]=GLOBAL;
4528         npw++;
4529
4530         pwl[npw]=ktup_pair_wise;
4531         sprintf (dpl[npw], "ktup_pair_wise");
4532         dps[npw]=LOCAL;
4533         npw++;
4534
4535         pwl[npw]=precomputed_pair_wise;
4536         sprintf (dpl[npw], "precomputed_pair_wise");
4537         dps[npw]=GLOBAL;
4538         npw++;
4539
4540         pwl[npw]=myers_miller_pair_wise;
4541         sprintf (dpl[npw], "default");
4542         dps[npw]=GLOBAL;
4543         npw++;
4544
4545         pwl[npw]=viterbi_pair_wise;
4546         sprintf (dpl[npw], "viterbi_pair_wise");
4547         dps[npw]=GLOBAL;
4548         npw++;
4549
4550         pwl[npw]=viterbiL_pair_wise;
4551         sprintf (dpl[npw], "viterbiL_pair_wise");
4552         dps[npw]=GLOBAL;
4553         npw++;
4554
4555         pwl[npw]=viterbiD_pair_wise;
4556         sprintf (dpl[npw], "viterbiD_pair_wise");
4557         dps[npw]=GLOBAL;
4558         npw++;
4559
4560         pwl[npw]=seq_viterbi_pair_wise;
4561         sprintf (dpl[npw], "seq_viterbi_pair_wise");
4562         dps[npw]=GLOBAL;
4563         npw++;
4564
4565         pwl[npw]=pavie_pair_wise;
4566         sprintf (dpl[npw], "pavie_pair_wise");
4567         dps[npw]=GLOBAL;
4568         npw++;
4569
4570         pwl[npw]=glocal_pair_wise;
4571         sprintf (dpl[npw], "glocal_pair_wise");
4572         dps[npw]=GLOBAL;
4573         npw++;
4574
4575         pwl[npw]=linked_pair_wise;
4576         sprintf (dpl[npw], "linked_pair_wise");
4577         dps[npw]=GLOBAL;
4578         npw++;
4579
4580         pwl[npw]=clinked_pair_wise;
4581         sprintf (dpl[npw], "clinked_pair_wise");
4582         dps[npw]=GLOBAL;
4583         npw++;
4584
4585         /*
4586         pwl[npw]=viterbiDGL_pair_wise;
4587         sprintf (dpl[npw], "viterbiDGL_pair_wise");
4588         dps[npw]=GLOBAL;
4589         npw++;
4590         */
4591       }
4592
4593     for ( a=0; a< npw; a++)
4594       {
4595         if ( (dp_mode && strm (dpl[a], dp_mode)) || pwl[a]==pw)
4596              {
4597                pw=pwl[a];
4598                if (dp_mode)sprintf (dp_mode,"%s", dpl[a]);
4599                glocal[0]=dps[a];
4600                return pw;
4601              }
4602       }
4603     fprintf ( stderr, "\n[%s] is an unknown mode for dp_mode[FATAL]\n", dp_mode);
4604     crash ( "\n");
4605     return NULL;
4606   }
4607
4608
4609 /*******************************************************************************/
4610 /*                                                                             */
4611 /*                                                                             */
4612 /*      Util Functions                                                         */
4613 /*                                                                             */
4614 /*                                                                             */
4615 /*******************************************************************************/
4616
4617 char *build_consensus ( char *seq1, char *seq2, char *dp_mode)
4618         {
4619         Alignment *A;
4620         char *buf;
4621         int a;
4622         char c1, c2;
4623         static char *mat;
4624
4625
4626         if ( !mat) mat=vcalloc ( STRING, sizeof (char));
4627
4628
4629         A=align_two_sequences (seq1, seq2, strcpy(mat,"idmat"), 0, 0,dp_mode);
4630         buf=vcalloc ( A->len_aln+1, sizeof (char));
4631
4632         for ( a=0; a< A->len_aln; a++)
4633             {
4634                 c1=A->seq_al[0][a];
4635                 c2=A->seq_al[1][a];
4636                 if (is_gap(c1) && is_gap(c2))buf[a]='-';
4637                 else if (is_gap(c1))buf[a]=c2;
4638                 else if (is_gap(c2))buf[a]=c1;
4639                 else if (c1!=c2){vfree (buf);buf=NULL;free_aln(A);return NULL;}
4640                 else buf[a]=c1;
4641             }
4642         buf[a]='\0';
4643         free_sequence (free_aln (A), -1);
4644         return buf;
4645         }
4646
4647
4648
4649 #ifdef FASTAL
4650
4651 combine_profile () Comobine two profiles into one, using the edit sequence produce by the DP
4652   edit_sequence () insert the gaps using the
4653
4654 int fastal (int argv, char **arg)
4655 {
4656   Sequence *S;
4657   int a, b;
4658   SeqHasch *H=NULL;
4659   int ktup=2;
4660
4661   S=get_fasta_sequence (arg[1], NULL);
4662
4663
4664
4665   for (a=0; a<S->nseq-1; a++)
4666     {
4667       for (b=a+1; b<S->nseq; b++)
4668         {
4669           mat[b][a]=mat[a][b]addrand()%100;
4670         }
4671     }
4672
4673   int_dist2nj_tree (s, S->name, S->nseq, tree_name);
4674   T=main_read_tree (BT);
4675   =fastal_tree_aln (T->L,T->R,S);
4676 }
4677
4678
4679
4680
4681 NT_node fastal_tree_aln ( NT_node P, Sequence *S)
4682 {
4683   int score;
4684
4685
4686   if (!P || P->nseq==1) return NULL;
4687   R=P->right;L=P->left;
4688
4689   fastal_tree_aln (P->left,S);
4690   fastal_tree_aln (P->right,S);
4691   fastal_pair_wise (P);
4692   return P;
4693 }
4694
4695
4696 NT_node fastal_pair_wise (NT_node P)
4697 {
4698   //X- 1
4699   //-X 2
4700   //XX 3
4701   //-- 4
4702
4703   tb=fastal_align_profile ((P->right)->prf, (P->left)->prf);
4704
4705   l=strlen (tb);
4706   for (a=0; a< l; a++)
4707     {
4708       pr1=pr2=0;
4709       if (tb[a]== 1 || tb[a] ==3)pr1=1;
4710       if (tb[a]== 2 || tb[a] ==3)pr2=1;
4711
4712       for (b=0; b<20; b++)
4713         P->prf[a][b]=((pr1==1)?(P->right)->prf[ppr1][b]:0) + ((pr2==1)?(P->left)->prf[ppr2][b]:0);
4714       ppr1+=pr1;
4715       ppr2+=pr2;
4716     }
4717   free_int ((P->left)->prf, -1);
4718   free_int ((P->right)->prf, -1);
4719 }
4720 #endif
4721 /*********************************COPYRIGHT NOTICE**********************************/
4722 /*© Centro de Regulacio Genomica */
4723 /*and */
4724 /*Cedric Notredame */
4725 /*Tue Oct 27 10:12:26 WEST 2009. */
4726 /*All rights reserved.*/
4727 /*This file is part of T-COFFEE.*/
4728 /**/
4729 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
4730 /*    it under the terms of the GNU General Public License as published by*/
4731 /*    the Free Software Foundation; either version 2 of the License, or*/
4732 /*    (at your option) any later version.*/
4733 /**/
4734 /*    T-COFFEE is distributed in the hope that it will be useful,*/
4735 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
4736 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
4737 /*    GNU General Public License for more details.*/
4738 /**/
4739 /*    You should have received a copy of the GNU General Public License*/
4740 /*    along with Foobar; if not, write to the Free Software*/
4741 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
4742 /*...............................................                                                                                      |*/
4743 /*  If you need some more information*/
4744 /*  cedric.notredame@europe.com*/
4745 /*...............................................                                                                                                                                     |*/
4746 /**/
4747 /**/
4748 /*      */
4749 /*********************************COPYRIGHT NOTICE**********************************/