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