binaries/src/tcoffee/t_coffee_source/tree_util.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <stdarg.h>
   5 #include <string.h>
   6 #include <ctype.h>
   7 #include "io_lib_header.h"
   8 #include "util_lib_header.h"
   9 #include "dp_lib_header.h"
  10 #include "define_header.h"
  11
  12 #define TOPOLOGY 1
  13 #define WEIGHTED 2
  14 #define LENGTH   3
  15 #define RECODE   4
  16
  17 int distance_tree;
  18 int rooted_tree;
  19 int tot_nseq;
  20 static NT_node compute_fj_tree (NT_node T, Alignment *A, int limit, char *mode);
  21 static NT_node compute_cw_tree (Alignment *A);
  22 static NT_node compute_std_tree (Alignment *A, int n, char **arg);
  23 static NT_node tree2fj_tree (NT_node T);
  24 int tree_contains_duplicates (NT_node T);
  25 int display_tree_duplicates (NT_node T);
  26
  27 static int compare_node1 ( int *b1, int *b2, int n);
  28 static int compare_node2 ( int *b1, int *b2, int n);
  29 static int find_seq_chain (Alignment *A, int **sim,int *used,int seq0,int seq1, int seq2,int chain_length, int limit, int max_chain, int *nseq);
  30 int new_display_tree (NT_node T, int n);
  31 NT_node display_code (NT_node T, int nseq, FILE *fp);
  32 NT_node display_dist (NT_node T, int n, FILE *fp);
  33 /*********************************************************************/
  34 /*                                                                   */
  35 /*                                   dpa_tree_manipulation           */
  36 /*                                                                   */
  37 /*********************************************************************/
  38 static NT_node code_dpa_tree ( NT_node T, int **D);
  39 NT_node collapse_sub_tree ( NT_node T,int nseq, int *list, char *new_name);
  40 NT_node seq2dpa_tree  (Sequence *S, char *mode)
  41 {
  42   Constraint_list *CL;
  43   NT_node **T;
  44   NT_node Tree;
  45   CL=declare_constraint_list_simple (S);
  46   CL->local_stderr=NULL;
  47
  48
  49   CL->DM=cl2distance_matrix (CL,NOALN,(mode==NULL)?"ktup":mode, NULL, 0);
  50
  51   T=int_dist2nj_tree ( (CL->DM)->similarity_matrix, S->name, S->nseq, vtmpnam (NULL));
  52   Tree=T[3][0];
  53
  54   Tree=recode_tree (Tree, S);
  55   Tree=reset_dist_tree (Tree, -1);
  56
  57   Tree=code_dpa_tree (Tree, (CL->DM)->similarity_matrix);
  58   free_distance_matrix (CL->DM);
  59   return Tree;
  60 }
  61
  62 NT_node tree2dpa_tree (NT_node T, Alignment *A, char *mode)
  63 {
  64   /*This Function sets the branches with Length values used by DP*/
  65   /*The tree must be rooted*/
  66   Sequence *S;
  67   int **D;
  68
  69   S=aln2seq (A);
  70   T=recode_tree     (T, S);
  71   T=reset_dist_tree (T, -1);
  72   D=get_sim_aln_array (A,mode);
  73
  74   T=code_dpa_tree (T, D);
  75   return T;
  76 }
  77
  78 NT_node code_dpa_tree ( NT_node T, int **D)
  79 {
  80   if ( !T) return T;
  81   else if ( T->leaf==1)
  82     {
  83       T->dist=100;
  84       return T;
  85     }
  86   else
  87     {
  88       int nl, *ll;
  89       int nr, *lr;
  90       int a, b, min=100;
  91       float tot, n=0;
  92
  93       nl=(T->left)->nseq;ll=(T->left)->lseq;
  94       nr=(T->right)->nseq;lr=(T->right)->lseq;
  95
  96       for (tot=0,n=0, a=0; a< nl; a++)
  97         for ( b=0; b< nr; b++, n++)
  98           {
  99             tot+=D[ll[a]][lr[b]];
 100             min=MIN(min,(D[ll[a]][lr[b]]));
 101           }
 102       /*      T->dist=(mode==AVERAGE)?(tot/n):min:;*/
 103       T->dist=(n>0)?tot/n:0;
 104       T->dist=min;
 105       code_dpa_tree ( T->right, D);
 106       code_dpa_tree ( T->left, D);
 107       return T;
 108     }
 109 }
 110 static int group_number;
 111 char *tree2Ngroup (Alignment *A, NT_node T, int max_n, char *fname, char *mat)
 112 {
 113   double top, bot, mid, pmid;
 114   Sequence *S;
 115   int n;
 116
 117
 118
 119   if (!T)
 120     {
 121       char **list;
 122
 123       list=declare_char ( 2, 100);
 124       sprintf (list[0], "%s",mat);
 125
 126       fprintf ( stderr, "\nCompute Phylogenetic tree [Matrix=%s]", mat);
 127       T=compute_std_tree(A,1, list);
 128       fprintf ( stderr, "\nCompute dpa tree");
 129       T=tree2dpa_tree (T,A, mat);
 130     }
 131
 132   S=tree2seq(T, NULL);
 133
 134   if ( max_n<0)
 135     {
 136       max_n*=-1;
 137       n=tree2group_file (T,S,0, max_n, fname);
 138       fprintf ( stderr, "\n#TrimTC: Split in %d Groups at a minimum of %d%% ID\n",n, (int)max_n);
 139       return fname;
 140
 141     }
 142   else if ( max_n>0)
 143     {
 144       if ( max_n>S->nseq)max_n=S->nseq;
 145
 146       top=100; bot=0;
 147       pmid=0; mid=50;
 148       n=tree2group_file(T, S,0, (int)mid,fname);
 149       mid=dichotomy((double)n, (double)max_n,(pmid=mid), &bot, &top);
 150       while (n!=max_n && (int)pmid!=(int)mid)
 151         {
 152           n=tree2group_file(T, S,0, (int)mid, fname);
 153           mid=dichotomy((double)n, (double)max_n,(pmid=mid), &bot, &top);
 154         }
 155       fprintf ( stderr, "\nDONE2");
 156       fprintf ( stderr, "\n#TrimTC: Split in %d Groups at a minimum of %d%% ID\n",n, (int)mid);
 157       return fname;
 158     }
 159   return NULL;
 160 }
 161 static int group_number;
 162 int tree2group_file ( NT_node T,Sequence *S, int maxnseq, int minsim, char *name)
 163   {
 164     FILE *fp;
 165
 166
 167     fp=vfopen (name, "w");
 168     vfclose (tree2group (T, S,maxnseq,minsim, "tree2ngroup",fp));
 169
 170     return count_n_line_in_file(name);
 171   }
 172
 173
 174 FILE * tree2group ( NT_node T,Sequence *S, int maxnseq, int minsim,char *name, FILE *fp)
 175 {
 176   if ( !T)return fp;
 177   else
 178     {
 179       int m,d;
 180
 181       m=(maxnseq==0)?S->nseq:maxnseq;
 182       d=minsim;
 183
 184
 185
 186       if ( T->nseq<=m && T->dist>=d)
 187         {
 188           int a;
 189           fprintf ( fp, ">%s_%d ", (name)?name:"", ++group_number);
 190           for ( a=0; a< T->nseq; a++)
 191             fprintf ( fp, "%s ", S->name[T->lseq[a]]);
 192           fprintf (fp, "\n");
 193           if (!T->parent)group_number=0;
 194           return fp;
 195         }
 196       else
 197         {
 198           fp=tree2group (T->right, S, maxnseq, minsim, name,fp);
 199           fp=tree2group (T->left, S, maxnseq, minsim, name,fp);
 200           if (!T->parent)group_number=0;
 201           return fp;
 202         }
 203
 204     }
 205 }
 206
 207
 208 NT_node  tree2collapsed_tree (NT_node T, int n, char **string)
 209 {
 210   char ***list;
 211   Sequence *A;
 212   int a, *nlist;
 213
 214
 215   A=tree2seq(T, NULL);
 216   T=recode_tree(T, A);
 217   list=vcalloc (A->nseq, sizeof (char***));
 218   nlist=vcalloc (A->nseq, sizeof (int));
 219   if ( n==0)return T;
 220   else if (n>1)
 221     {
 222       int l;
 223       char *buf;
 224
 225       for (l=0,a=0; a< n; a++)l+=strlen (string[a]);
 226       buf=vcalloc ( 2*n+l+1, sizeof (char));
 227       for (a=0; a< n; a++){buf=strcat (buf,string[a]), buf=strcat ( buf, " ");}
 228       list[0]=string2list (buf);
 229       vfree (buf);
 230     }
 231   else if ( file_exists (NULL,string[0]))
 232     {
 233       list=read_group (string[0]);
 234     }
 235   else
 236     {
 237       fprintf (stderr, "\nERROR: file <%s> does not exist [FATAL:%s]\n",string[0], PROGRAM);
 238       myexit (EXIT_FAILURE);
 239     }
 240
 241
 242   a=0;
 243   while (list[a])
 244     {
 245       int i, b;
 246       n=atoi (list[a][0]);
 247       for (b=0; b<A->nseq; b++)nlist[b]=0;
 248       for (b=2; b<n; b++)
 249         {
 250           i=name_is_in_list (list[a][b], A->name, A->nseq, MAXNAMES);
 251           nlist[i]=1;
 252         }
 253       T=collapse_sub_tree ( T,A->nseq,nlist,list[a][1]);
 254       free_char (list[a], -1);
 255       a++;
 256     }
 257   vfree (list);
 258   return T;
 259 }
 260
 261 NT_node collapse_sub_tree ( NT_node T,int nseq, int *list, char *new_name)
 262 {
 263   if (!T) return T;
 264   else
 265     {
 266       int a=0;
 267
 268
 269       while (a<nseq && list[a]==T->lseq2[a]){a++;}
 270       if (a==nseq)
 271         {
 272           sprintf ( T->name, "%s", new_name);
 273           T->leaf=T->isseq=1;
 274           T->left=T->right=NULL;
 275           return T;
 276         }
 277       else
 278         {
 279          collapse_sub_tree (T->right, nseq, list, new_name);
 280          collapse_sub_tree (T->left, nseq, list, new_name);
 281          return T;
 282         }
 283     }
 284 }
 285
 286 /*********************************************************************/
 287 /*                                                                   */
 288 /*                                   tree pruning                    */
 289 /*                                                                   */
 290 /*                                                                   */
 291 /*********************************************************************/
 292 NT_node remove_leaf ( NT_node T);
 293 NT_node prune_root (NT_node T);
 294 NT_node main_prune_tree ( NT_node T, Sequence *S)
 295 {
 296   T=prune_tree ( T, S);
 297   return T;
 298 }
 299
 300 NT_node prune_tree ( NT_node T, Sequence *S)
 301 {
 302
 303   if (!T ) return T;
 304
 305   if (T->leaf && T->isseq && name_is_in_list (T->name,S->name, S->nseq, 100)==-1)
 306     {
 307       NT_node C, P, PP;
 308
 309       P=T->parent;
 310       if ( !P)
 311         {
 312           int a;
 313           for (a=0; a< S->nseq; a++)
 314             {
 315               HERE ("prune pb ---%s", S->name[a]);
 316             }
 317           myexit (EXIT_FAILURE);
 318         }
 319       C=(P->right==T)?P->left:P->right;
 320       PP=C->parent=P->parent;
 321
 322       if (PP && PP->right==P)PP->right=C;
 323       else if (PP)PP->left=C;
 324       else
 325         {
 326           if (T==P->right)P->right=NULL;
 327           else P->left=NULL;
 328           T=C;
 329
 330         }
 331     }
 332   else
 333     {
 334       prune_tree (T->left, S);
 335       prune_tree (T->right, S);
 336     }
 337   return prune_root(T);
 338 }
 339
 340 NT_node prune_root (NT_node T)
 341 {
 342   //This function prunes the root if needed (and frees it).
 343   if (T->parent)return T;
 344
 345   if (!T->right && T->left)
 346     {
 347       return prune_root (T->left);
 348     }
 349   else if (T->right && !T->left)
 350      {
 351
 352        return prune_root (T->right);
 353     }
 354   else
 355     {
 356       return T;
 357     }
 358 }
 359 /*********************************************************************/
 360 /*                                                                   */
 361 /*                                   tree comparison                 */
 362 /*                                                                   */
 363 /*                                                                   */
 364 /*********************************************************************/
 365 int main_compare_cog_tree (NT_node T1, char *cogfile)
 366 {
 367   char ***array;
 368   int a, nbac, n=0, p, c, b;
 369   Alignment *A;
 370
 371   array=file2list(cogfile, ";\n");
 372   nbac=atoi(array[0][0])-2;
 373
 374   A=declare_aln2 (nbac+1, 10);
 375   for (a=0; a<nbac; a++)
 376     {
 377       sprintf ( A->name[a], "%s", array[0][a+2]);
 378       A->seq_al[a][0]='a';
 379       A->seq_al[a][1]='\0';
 380
 381     }
 382   sprintf ( A->name[nbac], "cons");
 383
 384   A->nseq=nbac+1;
 385   A->len_aln=1;
 386
 387
 388   n=3;
 389   while (array[n]!=NULL)
 390     {
 391       for (b=0; b<nbac; b++)
 392         {
 393           p=atoi (array[1][b+2]);p=(p==0)?'O':'I';
 394           c=atoi (array[n][b+2]);c=(c==0)?'O':'I';
 395           A->seq_al[b][0]=p;
 396           A->seq_al[b][1]=c;
 397           A->seq_al[b][2]='\0';
 398
 399         }
 400       sprintf (A->file[0], "%s", array[n][1]);
 401       A->len_aln=2;
 402       main_compare_aln_tree (T1, A, stdout);
 403       n++;
 404     }
 405   return n;
 406 }
 407
 408
 409 int main_compare_aln_tree (NT_node T1, Alignment *A, FILE *fp)
 410 {
 411   int n=0;
 412
 413   fprintf ( fp, "\nTOT_CLASH COG %s N %d", A->file[0], compare_aln_tree (T1, A, &n, fp));
 414   vfclose (fp);
 415   return n;
 416 }
 417
 418 int compare_aln_tree (NT_node T, Alignment *A, int *n, FILE *fp)
 419 {
 420   if (T->leaf)
 421     {
 422       int i;
 423       i=name_is_in_list (T->name, A->name, A->nseq, 100);
 424       T->seqal=A->seq_al[i];
 425       return 0;
 426     }
 427   else
 428     {
 429       char *seq1, *seq2;
 430       if (!(T->left )->seqal)compare_aln_tree (T->left, A,n, fp);
 431       if (!(T->right)->seqal)compare_aln_tree (T->right, A,n, fp);
 432
 433       seq1=(T->left)->seqal;
 434       seq2=(T->right)->seqal;
 435       (T->left)->seqal=(T->right)->seqal=NULL;
 436       if ( seq1 && seq2)
 437         {
 438           if (strm (seq1, seq2))
 439             {
 440               T->seqal=seq1;
 441
 442             }
 443           else
 444             {
 445
 446               if (seq1[0]!=seq2[0] && seq1[1]!=seq2[1])
 447                 {
 448                   fprintf ( fp, "\nNODE_CLASH: COG %s (%s,%s):(",A->file[0],seq1,seq2 );
 449                   display_leaf_below_node (T->left, fp);
 450                   fprintf ( fp, ");(");
 451                   display_leaf_below_node (T->right, fp);
 452                   fprintf ( fp, ")");
 453                   n[0]++;
 454                 }
 455             }
 456         }
 457     }
 458   return n[0];
 459 }
 460 //**********************************************************************
 461
 462 int compare_split (int *s1, int *s2, int l);
 463 int get_split_size (int *s, int l);
 464
 465 int main_compare_splits ( NT_node T1, NT_node T2, char *mode,FILE *fp)
 466 {
 467   Sequence *S1, *S2, *S;
 468   int  a, b;
 469
 470
 471
 472   int **sl1, n1;
 473   int **sl2, n2;
 474   if ( tree_contains_duplicates (T1))
 475     {
 476       display_tree_duplicates (T1);
 477       printf_exit (EXIT_FAILURE, stderr, "\nFirst Tree Contains Duplicated Sequences [main_compare_trees][FATAL:%s]", PROGRAM);
 478
 479     }
 480   else if ( tree_contains_duplicates (T2))
 481     {
 482       display_tree_duplicates (T2);
 483       printf_exit (EXIT_FAILURE, stderr, "\nSecond Tree Contains Duplicated Sequences [main_compare_trees]");
 484
 485     }
 486
 487   //Identify the commom Sequence Set
 488   S1=tree2seq(T1, NULL);
 489
 490
 491   S2=tree2seq(T2, NULL);
 492
 493
 494   S=trim_seq ( S1, S2);
 495
 496   //Prune the trees and recode the subtree list
 497   T1=prune_tree (T1, S);
 498   T1=recode_tree(T1, S);
 499
 500   T2=prune_tree (T2, S);
 501   T2=recode_tree(T2, S);
 502   HERE ("1");
 503   sl1=declare_int (10*S->nseq, S->nseq);
 504   sl2=declare_int (10*S->nseq, S->nseq);
 505
 506   HERE ("2");
 507   n1=n2=0;
 508   tree2split_list (T1, S->nseq, sl1, &n1);
 509   tree2split_list (T2, S->nseq, sl2, &n2);
 510
 511   for (a=0; a<n1; a++)
 512     {
 513       int n, best, s;
 514
 515       n=get_split_size (sl1[a], S->nseq);
 516       for (best=0,b=0; b<n2; b++)
 517         {
 518           s=compare_split (sl1[a], sl2[b], S->nseq);
 519           best=MAX(s,best);
 520         }
 521       fprintf ( fp, "\n%4d %4d ", MIN(n,(S->nseq)), best);
 522       for (b=0; b<S->nseq; b++)fprintf ( fp, "%d", sl1[a][b]);
 523     }
 524
 525   free_sequence (S, -1);
 526   free_sequence (S1, -1);
 527   free_sequence (S2, -1);
 528   myexit (EXIT_SUCCESS);
 529   return 1;
 530 }
 531 int compare_split (int *s1, int *s2, int l)
 532 {
 533   int n1, n2, score1, score2, a;
 534   n1=get_split_size (s1, l);
 535   n2=get_split_size (s2, l);
 536
 537   for (score1=0,a=0; a< l; a++)
 538     {
 539       score1+=(s1[a]==1 && s2[a]==1)?1:0;
 540     }
 541   score1=(score1*200)/(n1+n2);
 542
 543   for ( score2=0, a=0; a<l; a++)
 544     {
 545       score2+=(s1[a]==0 && s2[a]==1)?1:0;
 546     }
 547   score2=(score2*200)/((l-n1)+n2);
 548   return MAX(score1, score2);
 549 }
 550 int get_split_size (int *s, int l)
 551 {
 552   int a, b;
 553   for (a=b=0; a<l; a++)b+=s[a];
 554   return b;
 555 }
 556 //**********************************************************************
 557   //
 558   //
 559   //                       TREEE COMPARISON FUNCTIONS
 560   //
 561   //
 562   //
 563   //////////////////////////////////////////////////////////////////////
 564
 565 //JM_START
 566 void normalizeScore(float *score, int len)
 567 {
 568         int i;
 569         float SCORE_MIN = FLT_MAX;
 570         float SCORE_MAX = FLT_MIN;
 571         for(i = 0; i < len; i++)
 572         {
 573                 if(score[i] < SCORE_MIN)
 574                         SCORE_MIN = score[i];
 575                 if(score[i] > SCORE_MAX)
 576                         SCORE_MAX = score[i];
 577         }
 578         for(i = 0; i < len; i++)
 579                 score[i] = (9*(score[i]-SCORE_MIN)/(SCORE_MAX-SCORE_MIN));
 580 }
 581
 582 int new_compare_trees ( NT_node T1, NT_node T2, int nseq, Tree_sim *TS);
 583 NT_node new_search_split (NT_node T, NT_node B, int nseq);
 584 int new_compare_split ( int *b1, int *b2, int n);
 585 Tree_sim*  tree_scan_pos (Alignment *A, int start, int end, char *ptree, NT_node RT);
 586 Tree_sim*  tree_scan_pos_woble (Alignment *A, int center, int max, char *ptree, NT_node RT, int *br, int *bl );
 587 Tree_sim*  tree_scan_pair_pos (Alignment *A, int start, int end, int start2, int end2,char *ptree, NT_node RT);
 588 Tree_sim*  tree_scan_multiple_pos (int *poslist, int *wlist,int nl, Alignment *A, char *ptree, NT_node RT);
 589 NT_node aln2std_tree(Alignment *A, int ipara1, int ipara2, char *mode);
 590
 591 NT_node tree_scan (Alignment *A,NT_node RT, char *pscan, char *ptree)
 592 {
 593   int l, a,ax, c, cx, b;
 594   char mode[100];
 595   int start, w;
 596   int nl, *poslist;
 597   char posfile[100];
 598
 599   char *pcFileName = A->file[0];
 600   char prefix[200] ={0};
 601   int len = (strrchr(pcFileName,'.')?strrchr(pcFileName,'.')-pcFileName:strlen(pcFileName));
 602   strncpy(prefix, pcFileName, len);
 603
 604   float *fascore;
 605   char out_format[100];
 606
 607    char *score_csv_file = vcalloc(200, sizeof (char));
 608    char *score_html_file = vcalloc(200, sizeof (char));
 609    char *hit_matrix_file = vcalloc(200, sizeof (char));
 610    char *hit_html_file = vcalloc(200, sizeof (char));
 611    char *tree_file = vcalloc(200, sizeof (char));
 612
 613   sprintf(score_csv_file, "%s%s", prefix, ".score_csv");
 614   sprintf(score_html_file, "%s%s", prefix, ".ts_html");
 615   sprintf(hit_matrix_file, "%s%s", prefix, ".hit_matrix");
 616   sprintf(hit_html_file, "%s%s", prefix, ".hit_html");
 617   sprintf(tree_file, "%s%s", prefix, ".trees_txt");
 618
 619   if ( pscan && strstr ( pscan, "help"))
 620     {
 621       fprintf ( stdout, "\n+tree_scan| _W_     : Window size for the tree computation|STD size in norscan mode");
 622       fprintf ( stdout, "\n+tree_scan| _MODE_  : Mode for the number of windows (single, double, list, scan, pairscan, norscan, hit, norhit)");
 623       fprintf ( stdout, "\n+tree_scan| _MINW_  : Minimum Window size when using the scan mode (4)");
 624       fprintf ( stdout, "\n+tree_scan| _OUTTREE_ : specify the format of outputing tree in every position (default: not ouput)");
 625       myexit (EXIT_SUCCESS);
 626     }
 627
 628   strget_param (pscan, "_W_", "5", "%d",&w);
 629   strget_param (pscan, "_MODE_", "single", "%s",mode);
 630   strget_param (pscan, "_MINW_", "1", "%d",&start);
 631   strget_param (pscan, "_POSFILE_", "NO", "%s", posfile);
 632   strget_param (pscan, "_OUTTREE_", "", "%s", &out_format);
 633
 634         if(strlen(out_format) > 1)
 635                 unlink(tree_file);
 636
 637   l=intlen (A->len_aln);
 638
 639   poslist=vcalloc ( A->len_aln, sizeof (int));
 640   nl=0;
 641   fascore = vcalloc(A->len_aln, sizeof (float));
 642
 643   if ( strm (posfile, "NO"))
 644     {
 645
 646       for ( a=0; a< A->len_aln; a++)poslist[nl++]=a+1;
 647     }
 648   else
 649     {
 650       int *p;
 651       p=file2pos_list (A,posfile);
 652       poslist=pos2list (p, A->len_aln, &nl);
 653       for (a=0; a<nl; a++)poslist[a]++;
 654       vfree (p);
 655     }
 656
 657 //For tree hit
 658   NT_node *TreeArray = vcalloc(nl, sizeof (NT_node));
 659
 660   if ( strm (mode, "woble"))
 661     {
 662       for  (ax=0; ax<nl; ax++)
 663         {
 664           Tree_sim *TS;
 665           int left=0, right=0;
 666           a=poslist[ax];
 667
 668           TS=tree_scan_pos_woble (A,a,w,ptree, RT, &left, &right);
 669           fprintf ( stdout, "P: %*d I: %*d %*d SIM: %6.2f\n", l,a,l,left,l,right,TS->uw);
 670           vfree (TS);
 671         }
 672     }
 673   else if ( strm (mode, "single"))
 674     {
 675       for (b=0,ax=0; ax<nl; ax++)
 676         {
 677
 678           Tree_sim *TS;
 679           int pstart, pend;
 680
 681           a=poslist[ax];
 682
 683           pstart=a-b;
 684           pend=a+b;
 685
 686           if (pstart<1 || pstart>A->len_aln)continue;
 687           if (pend<1 || pend>A->len_aln)continue;
 688           TS=tree_scan_pos (A, pstart,pend, ptree, RT);
 689           fprintf ( stdout, "P: %*d I: %*d %*d SIM: %6.2f L: %2d\n", l,a,l,pstart,l,pend,TS->uw, (w*2)+1);
 690           vfree (TS);
 691         }
 692     }
 693   else if (strm (mode, "scan")||strm (mode, "hit"))
 694     {
 695         FILE *fp_ts;
 696         fp_ts=vfopen (score_csv_file, "w");
 697       fprintf ( fp_ts, "Position,Win_Beg,Win_End,Similarity,Win_Len\n");
 698       for ( ax=0; ax<nl; ax++)
 699         {
 700             float best_score=0;
 701             int best_pos=0, best_w=0, best_start, best_end;
 702
 703             a=poslist[ax];
 704             best_pos = best_start = best_end = a;
 705             for (b=start; b<=w; b++)
 706               {
 707                 Tree_sim *TS;
 708                 int pstart, pend;
 709                 pstart=a-b;
 710                 pend=a+b;
 711
 712                 if (pstart<1 || pstart>A->len_aln)continue;
 713                 if (pend<1 || pend>A->len_aln)continue;
 714                 TS=tree_scan_pos (A, pstart,pend, ptree, RT);
 715                 if (TS->uw>=best_score)
 716                         {best_score=TS->uw;best_w=b;best_start=pstart; best_end=pend;}
 717                 vfree (TS);
 718               }
 719             fprintf (fp_ts, "%*d,%*d,%*d,%6.2f,%2d\n", l,best_pos, l,best_start, l,best_end, best_score,(best_w*2)+1);
 720             fascore[ax]=(float)best_score;
 721                 if(strlen(out_format) > 1)
 722                         vfclose (print_tree (aln2std_tree(A, best_start, best_end, mode), out_format, vfopen (tree_file, "a+")));
 723                 if(strm (mode, "hit"))
 724                         TreeArray[ax] = aln2std_tree(A, best_start, best_end, mode);
 725           }
 726         vfclose(fp_ts);
 727     }
 728 //tree scan by using normal distribution window
 729 //or
 730 //generate hit matrix
 731   else if ( strm (mode, "norscan")||strm (mode, "norhit"))
 732     {
 733       FILE *fp_ts;
 734       ptree=vcalloc(100, sizeof (char));
 735       fp_ts=vfopen (score_csv_file, "w");
 736       fprintf ( fp_ts, "Position,Similarity,STD_Len\n");
 737       for ( ax=0; ax<nl; ax++)
 738           {
 739             float best_score=DBL_MIN;
 740             int best_STD = start;
 741             a=poslist[ax];
 742             for (b=start; b<=w; b++)
 743               {
 744                 Tree_sim *TS;
 745                 sprintf ( ptree, "+aln2tree _COMPARE_nordisaln__STD_%d__CENTER_%d_", b, a); //should be used a or ax
 746                 TS=tree_scan_pos (A, 1, nl, ptree, RT);
 747                 if (TS->uw>=best_score)
 748                         {best_score=TS->uw;best_STD=b;}
 749                 vfree (TS);
 750               }
 751               fascore[ax]=best_score;
 752               fprintf ( fp_ts, "%*d,%6.2f,%d\n", l,a, fascore[ax], best_STD);
 753                 if(strlen(out_format) > 1)
 754                         vfclose (print_tree (aln2std_tree(A, best_STD, a, mode), out_format, vfopen (tree_file, "a+")));
 755                 if(strm (mode, "norhit"))
 756                         TreeArray[ax] = aln2std_tree(A, best_STD, a, mode);
 757           }
 758         vfclose(fp_ts);
 759     }
 760 //generate hit matrix
 761         if (strm (mode, "hit")||strm (mode, "norhit"))
 762     {
 763 //Compute the pair score of tree scan segqtion
 764         fprintf (stdout, "[STRAT] Calculate the hit matrix of the tree scan\n");
 765         float **ffpHitScoreMatrix;
 766
 767         ffpHitScoreMatrix=vcalloc (nl, sizeof (float*));
 768         int i, j;
 769         for(i = 0; i < nl; i++)
 770                 ffpHitScoreMatrix[i]=vcalloc (nl-i, sizeof (float));
 771
 772         fprintf (stdout, "Process positions\n", i);
 773         for(i = 0; i < nl; i++)
 774         {
 775                 fprintf (stdout, "%d, ", i);
 776                 for(j = i; j < nl; j++)
 777                 {
 778                         Tree_sim *TS;
 779                         TS=tree_cmp (TreeArray[i], TreeArray[j]);
 780                         ffpHitScoreMatrix[i][j-i] = TS->uw;
 781                         vfree (TS);
 782                 }
 783         }
 784         vfree(TreeArray);
 785         fprintf (stdout, "\n");
 786         output_hit_matrix(hit_matrix_file, ffpHitScoreMatrix, nl);
 787         fprintf (stdout, "[END]Calculate the hit matrix of the tree scan\n");
 788
 789 //Output Hit Score into color html
 790         output_hit_color_html  (A, ffpHitScoreMatrix, nl, hit_html_file);
 791         vfree(ffpHitScoreMatrix);
 792     }
 793   else if ( strm (mode, "pairscan"))
 794     {
 795       int d, set;
 796
 797       for ( ax=0; ax<nl; ax++)
 798           {
 799             float best_score=0;
 800             int best_pos=0, best_w=0, best_w2=0, best_start, best_end, best_pos2=0, best_start2, best_end2;
 801
 802             Tree_sim *TS;
 803             int pstart, pend, p2start, p2end;
 804             a=poslist[ax];
 805             for ( cx=0; cx<nl; cx++)
 806               {
 807                 set=0;
 808                 c=poslist[cx];
 809                 for (d=start; d<=w; d++)
 810                   {
 811                     for (b=start; b<=w; b++)
 812                       {
 813                         pstart=a-b;
 814                         pend=a+b;
 815                         p2start=c-d;
 816                         p2end=c+d;
 817                         if (pstart<1 || pstart>A->len_aln)continue;
 818                         if (pend<1 || pend>A->len_aln)continue;
 819                         if (p2start<1 || p2start>A->len_aln)continue;
 820                         if (p2end<1 || p2end>A->len_aln)continue;
 821                         if (pstart<=p2start && pend>=p2start) continue;
 822                         if (pstart<=p2end && pend>=p2end) continue;
 823                         TS=tree_scan_pair_pos (A, pstart,pend,p2start, p2end, ptree, RT);
 824
 825                         if (TS->uw>=best_score){best_score=TS->uw; best_pos=a;best_w=b;best_start=pstart; best_end=pend; best_pos2=c, best_w2=d, best_start2=p2start, best_end2=p2end;set=1;}
 826                         vfree (TS);
 827                       }
 828                   }
 829                 if (set)fprintf ( stdout, "P1: %*d  I1: %*d %*d P2: %*d I2: %*d %*d SIM: %6.2f L: %2d\n", l,best_pos, l,best_start, l,best_end, l, best_pos2, l, best_start2, l, best_end2, best_score,(best_w*2)+1 );
 830
 831                 set=0;
 832               }
 833           }
 834     }
 835   else if ( strm (mode, "multiplescan"))
 836     {
 837       int n, **wlist, best_pos;
 838       float best_score;
 839       Tree_sim *TS;
 840       wlist=generate_array_int_list (nl*2,start, w,1, &n, NULL);
 841       HERE ("Scan %d Possibilities", n);
 842
 843       for (best_score=best_pos=0,a=0; a<n; a++)
 844         {
 845           TS=tree_scan_multiple_pos (poslist,wlist[a],nl, A, ptree, RT);
 846           if (TS && TS->uw>best_score)
 847             {
 848               best_score=TS->uw;
 849               fprintf ( stdout, "\n");
 850               for (b=0; b<nl; b++)
 851                 {
 852                   fprintf ( stdout, "[%3d %3d %3d]", poslist[b], wlist[a][b*2], wlist[a][b*2+1]);
 853                 }
 854               fprintf ( stdout, " SCORE: %.2f", best_score);
 855             }
 856           if (TS)vfree (TS);
 857         }
 858     }
 859
 860 //Output Tree Scan core into color html
 861     normalizeScore(fascore, nl);
 862     Alignment *ST;
 863     Sequence *S;
 864     int i, r1;
 865
 866     S=A->S;
 867     ST=copy_aln (A, NULL);
 868     for (a=0; a<ST->nseq; a++)
 869     {
 870         i=name_is_in_list (ST->name[a],S->name, S->nseq, 100);
 871         if ( i!=-1)
 872         {
 873                 for (b=0; b<ST->len_aln; b++)
 874                 {
 875                         r1=ST->seq_al[a][b];
 876                         if ( r1!='-')
 877                                 r1 = (int)fascore[b] + 48;
 878                         ST->seq_al[a][b]=r1;
 879                 }
 880         }
 881     }
 882     output_color_html  ( A, ST, score_html_file);
 883
 884 //free memory
 885         free_aln(ST);
 886         vfree(fascore);
 887   vfree(score_csv_file);
 888   vfree(score_html_file);
 889   vfree(hit_matrix_file);
 890   vfree(hit_html_file);
 891
 892   myexit(EXIT_SUCCESS);return NULL;
 893 }
 894
 895 NT_node aln2std_tree(Alignment *A, int ipara1, int ipara2, char *mode)
 896 {
 897      Alignment *B;
 898      NT_node T;
 899      char *cpSet = vcalloc(100, sizeof (char));
 900
 901         if(strm (mode, "norhit"))
 902         {
 903                 B=extract_aln (A, 1, A->len_aln);
 904                 sprintf ( cpSet, "+aln2tree _COMPARE_nordisaln__STD_%d__CENTER_%d_", ipara1, ipara2);
 905         }
 906         else
 907                 B=extract_aln (A, ipara1, ipara2);
 908
 909         T=compute_std_tree (B, (cpSet)?1:0, (cpSet)?&cpSet:NULL);
 910         free_aln(B);
 911         return T;
 912 }
 913
 914 Tree_sim*  tree_scan_multiple_pos (int *poslist, int *wlist,int nl, Alignment *A, char *ptree, NT_node RT)
 915 {
 916     static Alignment *B;
 917     static int *pos;
 918     int a, b, n, s, p, left, right;
 919     Tree_sim *TS;
 920     NT_node T=NULL;
 921
 922
 923     //poslist positions come [1..n]
 924     vfree(pos);
 925     free_aln (B);
 926
 927     pos=vcalloc ( A->len_aln+1, sizeof (int));
 928     B=copy_aln (A, NULL);
 929
 930     for (a=0; a<nl; a++)
 931       {
 932         p=poslist[a]-1;
 933         left =wlist[2*a];
 934         right=wlist[2*a+1];
 935         for (b=p-left; b<=p+right; b++)
 936           {
 937             if (b<1 ||b>A->len_aln) return NULL;
 938             else pos[b]++;
 939
 940             if (pos[b]>1) return NULL;
 941           }
 942       }
 943
 944     for (s=0; s<A->nseq; s++)
 945       {
 946         for (n=0,a=1; a<=A->len_aln; a++)
 947           {
 948             if (pos[a])B->seq_al[s][n++]=A->seq_al[s][a-1];
 949           }
 950       }
 951
 952     B->len_aln=n;
 953     for (s=0; s<A->nseq; s++)B->seq_al[s][B->len_aln]='\0';
 954
 955     T=compute_std_tree (B, (ptree)?1:0, (ptree)?&ptree:NULL);
 956
 957     TS=tree_cmp (T, RT);
 958
 959     free_tree(T);
 960     return TS;
 961   }
 962
 963 Tree_sim*  tree_scan_pair_pos (Alignment *A, int start, int end, int start2, int end2,char *ptree, NT_node RT)
 964   {
 965     Tree_sim *TS;
 966     Alignment *B,*B1, *B2;
 967     NT_node T=NULL;
 968     int a;
 969
 970
 971     B=copy_aln (A, NULL);
 972      B1=extract_aln (A,start,end);
 973      B2=extract_aln (A,start2, end2);
 974
 975
 976      for ( a=0; a< B->nseq;a++)
 977        sprintf (B->seq_al[a], "%s%s", B1->seq_al[a], B2->seq_al[a]);
 978      B->len_aln=strlen (B->seq_al[0]);
 979
 980      T=compute_std_tree (B, (ptree)?1:0, (ptree)?&ptree:NULL);
 981      TS=tree_cmp (T, RT);
 982
 983      free_tree(T);
 984      free_aln (B);free_aln(B1); free_aln(B2);
 985
 986      return TS;
 987   }
 988
 989 Tree_sim*  tree_scan_pos (Alignment *A, int start, int end, char *ptree, NT_node RT)
 990   {
 991      Tree_sim *TS;
 992      Alignment *B;
 993      NT_node T;
 994
 995      if ( start<1 || start>A->len_aln) return NULL;
 996      if ( end<1 || end>A->len_aln) return NULL;
 997
 998      B=extract_aln (A,start,end);
 999      T=compute_std_tree (B, (ptree)?1:0, (ptree)?&ptree:NULL);
1000      TS=tree_cmp (T, RT);
1001      free_tree(T);free_aln (B);
1002      return TS;
1003   }
1004 Tree_sim*  tree_scan_pos_woble (Alignment *A, int center, int max, char *ptree, NT_node RT, int *br, int *bl )
1005   {
1006     Tree_sim *TS,*BTS;
1007
1008
1009      int left, right;
1010      float best_score=0;
1011      int start, end;
1012
1013      br[0]=bl[0]=0;
1014      BTS=vcalloc (1, sizeof (Tree_sim));
1015
1016      for (left=0; left<max; left++)
1017        for (right=0; right<max; right++)
1018          {
1019            start=center-left;
1020            end=center+right;
1021            TS=tree_scan_pos (A, start, end, ptree, RT);
1022
1023            if (TS && TS->uw >best_score)
1024              {
1025                best_score=TS->uw;
1026                BTS[0]=TS[0];
1027                br[0]=right; bl[0]=left;
1028                vfree(TS);
1029              }
1030          }
1031      return BTS;
1032   }
1033
1034 Tree_sim* tree_cmp( NT_node T1, NT_node T2)
1035 {
1036   Sequence *S1, *S2, *S;
1037   int n;
1038   int a;
1039
1040   Tree_sim *TS1, *TS2;
1041
1042   if ( tree_contains_duplicates (T1))
1043     {
1044       display_tree_duplicates (T1);
1045       printf_exit (EXIT_FAILURE, stderr, "\nFirst Tree Contains Duplicated Sequences [main_compare_trees][FATAL:%s]", PROGRAM);
1046
1047     }
1048   else if ( tree_contains_duplicates (T2))
1049     {
1050       display_tree_duplicates (T2);
1051       printf_exit (EXIT_FAILURE, stderr, "\nSecond Tree Contains Duplicated Sequences [main_compare_trees]");
1052
1053     }
1054
1055   //Identify the commom Sequence Set
1056   S1=tree2seq(T1, NULL);
1057   S2=tree2seq(T2, NULL);
1058
1059   S=trim_seq ( S1, S2);
1060
1061   if ( S->nseq<=2)
1062     {
1063       fprintf ( stderr, "\nERROR: Your two do not have enough common leaf to be compared [FATAL:PROGRAM]");
1064     }
1065
1066   //Prune the trees and recode the subtree list
1067   T1=prune_tree (T1, S);
1068   T1=recode_tree(T1, S);
1069
1070   T2=prune_tree (T2, S);
1071   T2=recode_tree(T2, S);
1072
1073   TS1=vcalloc (1, sizeof (Tree_sim));
1074   TS2=vcalloc (1, sizeof (Tree_sim));
1075
1076   new_compare_trees ( T1, T2, S->nseq, TS1);
1077   new_compare_trees ( T2, T1, S->nseq, TS2);
1078
1079
1080   TS1->n=tree2nnode (T1);
1081   TS1->nseq=S->nseq;
1082
1083   TS2->n=tree2nnode (T2);
1084   /*if (TS1->n !=TS2->n)
1085     printf_exit (EXIT_FAILURE, stderr,"\nERROR: Different number of Nodes in the two provided trees after prunning [FATAL: %s]", PROGRAM);
1086   */
1087
1088   free_sequence (S, -1);
1089   free_sequence (S1, -1);
1090   free_sequence (S2, -1);
1091
1092   TS1->uw=(TS1->uw+TS2->uw)*100/(TS1->max_uw+TS2->max_uw);
1093   TS1->w=(TS1->w+TS2->w)*100/(TS1->max_w+TS2->max_w);
1094   TS1->d=(TS1->d+TS2->d)*100/(TS1->max_d+TS2->max_d);
1095   TS1->rf=(TS1->rf+TS2->rf)/2;
1096   vfree (TS2);
1097   return TS1;
1098 }
1099 int print_node_list (NT_node T, Sequence *RS)
1100 {
1101   NT_node *L;
1102   Sequence *S;
1103   int *nlseq2;
1104
1105   S=tree2seq(T, NULL);
1106   L=tree2node_list (T, NULL);
1107   if (!RS)RS=S;
1108
1109   nlseq2=vcalloc ( RS->nseq, sizeof (int));
1110   while (L[0])
1111     {
1112       int d,b;
1113       d=MIN(((L[0])->nseq), (S->nseq-(L[0])->nseq));
1114       fprintf ( stdout, "Bootstrap: %5.2f Depth: %5d Splits: ", (L[0])->bootstrap, d);
1115
1116       for (b=0; b<RS->nseq; b++)nlseq2[b]='-';
1117       for (b=0; b<S->nseq; b++)
1118         {
1119           int p;
1120           p=name_is_in_list (S->name[b], RS->name, RS->nseq, 100);
1121           if (p!=-1)nlseq2[p]=(L[0])->lseq2[b]+'0';
1122         }
1123       for (b=0; b<RS->nseq; b++) fprintf ( stdout, "%c", nlseq2[b]);
1124       fprintf (stdout, "\n");
1125       L++;
1126     }
1127   return 1;
1128 }
1129 NT_node main_compare_trees_list ( NT_node RT, Sequence *S, FILE *fp)
1130 {
1131   Tree_sim *T;
1132   NT_node *TL;
1133   Sequence *RS;
1134   int a;
1135
1136   T=vcalloc (1, sizeof (Tree_sim));
1137   RS=tree2seq(RT, NULL);
1138
1139   TL=read_tree_list (S);
1140
1141   reset_boot_tree (RT,0.0001);
1142   for (a=0; a<S->nseq; a++)
1143     {
1144       TL[a]=prune_tree(TL[a],RS);
1145       TL[a]=recode_tree (TL[a],RS);
1146       new_compare_trees (RT, TL[a], RS->nseq,T);
1147     }
1148
1149   vfree (T);
1150   return RT;
1151 }
1152 NT_node main_compare_trees ( NT_node T1, NT_node T2, FILE *fp)
1153 {
1154   Tree_sim *T;
1155
1156   T=tree_cmp (T1, T2);
1157   fprintf ( fp, "\n#tree_cmp|T: %.f W: %.2f L: %.2f RF: %d N: %d S: %d", T->uw, T->w, T->d, T->rf, T->n, T->nseq);
1158   fprintf ( fp, "\n#tree_cmp_def|T: ratio of identical nodes");
1159   fprintf ( fp, "\n#tree_cmp_def|W: ratio of identical nodes weighted with the min Nseq below node");
1160   fprintf ( fp, "\n#tree_cmp_def|L: average branch length similarity");
1161   fprintf ( fp, "\n#tree_cmp_def|RF: Robinson and Foulds");
1162   fprintf ( fp, "\n#tree_cmp_def|N: number of Nodes in T1 [unrooted]");
1163   fprintf ( fp, "\n#tree_cmp_def|S: number of Sequences in T1\n");
1164
1165   vfree (T);
1166   return T1;
1167 }
1168
1169 int new_compare_trees ( NT_node T1, NT_node T2, int nseq, Tree_sim *TS)
1170 {
1171   int n=0;
1172   NT_node N;
1173   float t1, t2;
1174
1175   if (!T1 || !T2) return 0;
1176
1177   n+=new_compare_trees (T1->left,  T2, nseq,TS);
1178   n+=new_compare_trees (T1->right, T2, nseq,TS);
1179
1180   //Exclude arbitrary splits (dist==0)
1181   if ((T1->dist==0) && !(T1->parent))return n;
1182
1183   N=new_search_split (T1, T2, nseq);
1184   t1=FABS(T1->dist);
1185   t2=(N)?FABS(N->dist):0;
1186   TS->max_d+=MAX(t1, t2);
1187
1188   if (!N)TS->rf++;
1189   if (T1->nseq>1)
1190     {
1191       int w;
1192       w=MIN((nseq-T1->nseq),T1->nseq);
1193       TS->max_uw++;
1194       TS->max_w+=w;
1195
1196       if (N)
1197         {
1198           TS->uw++;
1199           TS->w+=w;
1200           TS->d+=MIN(t1, t2);
1201           T1->bootstrap++;
1202           //T1->dist=T1->nseq;
1203         }
1204       else
1205         {
1206           //T1->dist=T1->nseq*-1;
1207           ;
1208         }
1209     }
1210   else
1211     {
1212       TS->d+=MIN(t1, t2);
1213       //T1->dist=1;
1214     }
1215   return ++n;
1216 }
1217 NT_node new_search_split (NT_node T, NT_node B, int nseq)
1218 {
1219   NT_node N;
1220   if (!T || !B) return NULL;
1221   else if ( new_compare_split (T->lseq2, B->lseq2, nseq)==1)return B;
1222   else if ( (N=new_search_split (T, B->right, nseq)))return N;
1223   else return new_search_split (T, B->left, nseq);
1224 }
1225 int new_compare_split ( int *b1, int *b2, int n)
1226 {
1227   int a, flag;
1228
1229   for (flag=1, a=0; a<n; a++)
1230     if (b1[a]!=b2[a])flag=0;
1231   if (flag) return flag;
1232
1233   for (flag=1, a=0; a<n; a++)
1234     if (b1[a]==b2[a])flag=0;
1235   return flag;
1236 }
1237
1238 float compare_trees ( NT_node T1, NT_node T2, int nseq,int  mode)
1239 {
1240   /*search each branch of T1 in T2*/
1241   float n=0;
1242
1243
1244   if ( !T1 || !T2)return 0;
1245
1246   if (getenv4debug("DEBUG_TREE_COMPARE"))display_node (T1, "\nNODE ", nseq);
1247
1248   if (T1->parent && T1->nseq>1)n+=search_node ( T1, T2, nseq, mode);
1249
1250   n+=compare_trees ( T1->left, T2, nseq, mode);
1251   n+=compare_trees ( T1->right, T2, nseq, mode);
1252
1253   return n;
1254 }
1255
1256 float search_node ( NT_node B, NT_node T, int nseq, int mode)
1257 {
1258   int n=0;
1259   if ( !B || !T) return -1;
1260   if (getenv4debug("DEBUG_TREE_COMPARE"))display_node ( T, "\n\t", nseq);
1261
1262   n=compare_node ( B->lseq2, T->lseq2, nseq );
1263
1264   if ( n==1)
1265     {
1266       if (getenv4debug("DEBUG_TREE_COMPARE"))fprintf ( stderr, "[1][%d]", (int)evaluate_node_similarity ( B, T, nseq, mode));
1267       if (mode==RECODE)B->dist=B->leaf;
1268       return evaluate_node_similarity ( B, T, nseq, mode);
1269     }
1270   else if ( n==-1)
1271     {
1272       if (getenv4debug("DEBUG_TREE_COMPARE"))fprintf ( stderr, "[-1]");
1273       if (mode==RECODE)B->dist=-B->leaf;
1274       return 0;
1275     }
1276   else
1277     {
1278       if (getenv4debug("DEBUG_TREE_COMPARE"))fprintf ( stderr, "[0]");
1279       n=search_node ( B, T->left, nseq, mode);
1280       if ( n>0) return n;
1281       n=search_node ( B, T->right, nseq, mode);
1282       if ( n>0) return n;
1283       n=search_node ( B, T->bot, nseq, mode);
1284       if ( n>0) return n;
1285     }
1286   return n;
1287 }
1288
1289 float evaluate_node_similarity ( NT_node B, NT_node T, int nseq, int mode)
1290 {
1291 int a, c;
1292
1293  if ( mode==TOPOLOGY || mode ==RECODE)
1294    {
1295      for ( a=0; a< nseq; a++)
1296        if ( B->lseq2[a]!=T->lseq2[a]) return 0;
1297      return 1;
1298    }
1299  else if ( mode == WEIGHTED)
1300    {
1301      for (c=0, a=0; a< nseq; a++)
1302        {
1303          if ( B->lseq2[a]!=T->lseq2[a]) return 0;
1304          else c+=B->lseq2[a];
1305        }
1306      return (float)(MIN(c,nseq));
1307    }
1308  else if ( mode == LENGTH )
1309    {
1310      float d1, d2;
1311
1312      for (c=0, a=0; a< nseq; a++)
1313        {
1314          if ( B->lseq2[a]!=T->lseq2[a]) return 0;
1315        }
1316      d1=FABS((B->dist-T->dist));
1317      d2=MAX(B->dist, T->dist);
1318      return (d2>0)?(d1*100)/d2:0;
1319    }
1320  else
1321    {
1322      return 0;
1323    }
1324 }
1325 int compare_node ( int *b1, int *b2, int nseq)
1326 {
1327   int n1, n2;
1328
1329   n1=compare_node1 ( b1, b2, nseq);
1330   /*fprintf ( stderr, "[%d]", n1);*/
1331   if ( n1==1) return 1;
1332
1333   n2=compare_node2 ( b1, b2, nseq);
1334   /* fprintf ( stderr, "[%d]", n2);*/
1335   if ( n2==1)return 1;
1336   else if ( n2==-1 && n1==-1) return -1;
1337   else return 0;
1338 }
1339 int compare_node1 ( int *b1, int *b2, int n)
1340 {
1341   int a;
1342   int l1, l2;
1343   int r=1;
1344   for ( a=0; a< n; a++)
1345     {
1346       l1=b1[a];
1347       l2=b2[a];
1348       if ( l1==1 && l2==0) return -1;
1349       if ( l1!=l2)r=0;
1350     }
1351   return r;
1352 }
1353 int compare_node2 ( int *b1, int *b2, int n)
1354 {
1355   int a;
1356   int l1, l2;
1357   int r=1;
1358
1359   for ( a=0; a< n; a++)
1360     {
1361       l1=1-b1[a];
1362       l2=b2[a];
1363       if ( l1==1 && l2==0) return -1;
1364       if ( l1!=l2) r=0;
1365     }
1366   return r;
1367 }
1368 void display_node (NT_node N, char *string,int nseq)
1369 {
1370   int a;
1371   fprintf ( stderr, "%s", string);
1372   for (a=0; a< nseq; a++)fprintf ( stderr, "%d", N->lseq2[a]);
1373 }
1374
1375
1376 /*********************************************************************/
1377 /*                                                                   */
1378 /*                                   FJ_tree Computation             */
1379 /*                                                                   */
1380 /*                                                                   */
1381 /*********************************************************************/
1382
1383
1384 NT_node tree_compute ( Alignment *A, int n, char ** arg_list)
1385 {
1386   if (n==0 || strm (arg_list[0], "cw"))
1387     {
1388       return compute_cw_tree (A);
1389     }
1390   else if ( strm (arg_list[0], "fj"))
1391     {
1392       return compute_fj_tree ( NULL, A, (n>=1)?atoi(arg_list[1]):8, (n>=2)?arg_list[2]:NULL);
1393     }
1394
1395   else if ( ( strm (arg_list[0], "nj")))
1396     {
1397       return compute_std_tree (A, n, arg_list);
1398     }
1399   else
1400     return compute_std_tree (A, n, arg_list);
1401 }
1402
1403 NT_node compute_std_tree (Alignment *A, int n, char **arg_list)
1404 {
1405   return compute_std_tree_2 (A, NULL, list2string (arg_list, n));
1406 }
1407 NT_node compute_std_tree_2 (Alignment *A, int **s, char *cl)
1408 {
1409
1410   NT_node T, **BT=NULL;
1411   char *tree_name;
1412
1413   char matrix[100];
1414   char score [100];
1415   char compare[100];
1416   char tmode[100];
1417   int free_s=0;
1418
1419   tree_name =vtmpnam (NULL);
1420
1421   if (strstr (cl, "help"))
1422     {
1423       fprintf ( stdout, "\n+aln2tree| _MATRIX_ : matrix used for the comparison (idmat, sarmat, pam250mt..)\n");
1424       fprintf ( stdout, "\n+aln2tree| _SCORE_  :  score mode used for the distance (sim, raw)\n");
1425       fprintf ( stdout, "\n+aln2tree| _COMPARE_:  comparison mode (aln, ktup, align, nordisaln)\n");
1426       fprintf ( stdout, "\n+aln2tree| _TMODE_  :  tree mode (nj, upgma)\n");
1427       myexit (EXIT_SUCCESS);
1428     }
1429
1430
1431   //matrix: idmat, ktup,sarmat, sarmat2
1432   strget_param (cl, "_MATRIX_", "idmat", "%s",matrix);
1433
1434   //score: sim, raw
1435   strget_param (cl, "_SCORE_", "sim", "%s",score);
1436
1437   //compare: aln, ktup, align
1438   strget_param (cl, "_COMPARE_", "aln", "%s",compare);
1439
1440   //compare: aln, ktup, align
1441   strget_param (cl, "_TMODE_", "nj", "%s",tmode);
1442
1443   int STD, CENTER;
1444   if ( strm (compare, "nordisaln"))
1445   {
1446         strget_param (cl, "_STD_", "1", "%d", &STD);
1447         strget_param (cl, "_CENTER_", "5", "%d", &CENTER);
1448   }
1449   //Use external msa2tree methods
1450   if ( strm (tmode, "cw"))
1451     {
1452       free_int (s, -1);
1453       return compute_cw_tree (A);
1454     }
1455
1456   //compute distance matrix if needed
1457   if ( !s)
1458     {
1459       free_s=1;
1460     if ( strm (compare, "ktup"))
1461         {
1462           ungap_array (A->seq_al, A->nseq);
1463           s=get_sim_aln_array ( A,cl);
1464         }
1465       else if ( strm ( compare, "aln"))
1466         {
1467           if (strm (score, "sim"))
1468             s=get_sim_aln_array(A, matrix);
1469           else if ( strm (score, "raw"))
1470             {
1471               s=get_raw_sim_aln_array (A,matrix);
1472             }
1473         }
1474       else if ( strm ( compare, "nordisaln"))
1475         {
1476           s=get_sim_aln_array_normal_distribution(A, matrix, &STD, &CENTER);
1477         }
1478       s=sim_array2dist_array(s, 100);
1479     }
1480
1481   //Compute the tree
1482   if (strm (tmode, "nj"))
1483     {
1484
1485       BT=int_dist2nj_tree (s, A->name, A->nseq, tree_name);
1486       T=main_read_tree (tree_name);
1487       free_read_tree(BT);
1488     }
1489   else if (strm (tmode, "upgma"))
1490     {
1491       BT=int_dist2upgma_tree (s,A, A->nseq, tree_name);
1492       T=main_read_tree (tree_name);
1493       free_read_tree(BT);
1494     }
1495
1496   if ( strm ( cl, "dpa"))
1497     {
1498       s=dist_array2sim_array(s, 100);
1499       T=code_dpa_tree (T,s);
1500     }
1501
1502   if (free_s)free_int (s, -1);
1503   return T;
1504 }
1505
1506 NT_node similarities_file2tree (char *mat)
1507 {
1508   int **s;
1509   Alignment *A;
1510   char *tree_name;
1511   NT_node T;
1512
1513
1514
1515   tree_name =vtmpnam (NULL);
1516
1517   s=input_similarities (mat,NULL, NULL);
1518
1519
1520   A=similarities_file2aln(mat);
1521   s=sim_array2dist_array(s, 100);
1522
1523
1524   int_dist2nj_tree (s, A->name, A->nseq, tree_name);
1525   T=main_read_tree(tree_name);
1526   free_int (s, -1);
1527   return T;
1528 }
1529
1530 NT_node compute_cw_tree (Alignment *A)
1531 {
1532   char *tmp1, *tmp2, tmp3[1000];
1533
1534   tmp1=vtmpnam (NULL);
1535   tmp2=vtmpnam (NULL);
1536
1537   sprintf ( tmp3, "%s.ph", tmp1);
1538   output_clustal_aln (tmp1, A);
1539   printf_system ("clustalw -infile=%s -tree -newtree=%s %s ", tmp1,tmp3, TO_NULL_DEVICE);
1540   printf_system("mv %s %s", tmp3, tmp2);
1541
1542   return main_read_tree(tmp2);
1543 }
1544
1545 NT_node compute_fj_tree (NT_node T, Alignment *A, int limit, char *mode)
1546 {
1547   static int in_fj_tree;
1548   if (!in_fj_tree)fprintf ( stderr, "\nComputation of an NJ tree using conserved positions\n");
1549
1550   in_fj_tree++;
1551   if (T && T->leaf<=2);
1552   else
1553     {
1554       T=aln2fj_tree(T,A,limit, mode);
1555       T->right=compute_fj_tree ( T->right, A, limit, mode);
1556       T->left=compute_fj_tree  ( T->left, A, limit, mode);
1557     }
1558   in_fj_tree--;
1559   return T;
1560 }
1561
1562
1563 NT_node aln2fj_tree(NT_node T, Alignment *A, int limit_in, char *mode)
1564 {
1565   NT_node NT;
1566   Sequence *S=NULL;
1567   Alignment *subA=NULL;
1568   int fraction_gap;
1569   int l, limit;
1570
1571   if (T)
1572     S=tree2seq (T,NULL);
1573   else
1574     S=aln2seq (A);
1575
1576
1577   l=0;
1578   for ( fraction_gap=100; fraction_gap<=100 && l<1; fraction_gap+=10)
1579     for ( limit=limit_in; limit>0 && l<1; limit--)
1580       {
1581         fprintf ( stderr, "\n%d %d", limit, fraction_gap);
1582         free_aln (subA);
1583         subA=extract_sub_aln2 (A,S->nseq,S->name);
1584         subA=filter_aln4tree (subA,limit,fraction_gap, mode);
1585         l=subA->len_aln;
1586       }
1587
1588   /*  while ( subA->len_aln<1)
1589     {
1590     subA=extract_sub_aln2 (A,S->nseq,S->name);
1591     subA=filter_aln4tree (subA,limit,fraction_gap,mode);
1592     free_aln (subA);
1593     subA=extract_sub_aln2 (A,S->nseq,S->name);
1594     subA=filter_aln4tree (subA,--limit,fraction_gap, mode);
1595     }
1596   */
1597   NT=aln2tree (subA);
1598   NT=tree2fj_tree (NT);
1599
1600   NT=realloc_tree (NT,A->nseq);
1601   fprintf ( stderr, "Limit:%d Gap: %d Columns: %4d Left: %4d Right %4d BL:%4.2f\n",limit,fraction_gap,  subA->len_aln, (NT->right)->leaf,(NT->left)->leaf, (NT->left)->dist+(NT->right)->dist);
1602
1603   if ( T)
1604     {
1605       NT->dist=T->dist;
1606       NT->parent=T->parent;
1607     }
1608   free_tree(T);
1609   free_aln (subA);
1610   free_sequence (S, -1);
1611   return NT;
1612 }
1613
1614 Alignment * filter_aln4tree (Alignment *A, int n,int fraction_gap,char *mode)
1615 {
1616   char *aln_file;
1617   char *ungaped_aln_file;
1618   char *scored_aln_file;
1619   char *filtered_aln_file;
1620
1621   aln_file=vtmpnam(NULL);
1622   ungaped_aln_file=vtmpnam (NULL);
1623   scored_aln_file=vtmpnam (NULL);
1624   scored_aln_file=vtmpnam(NULL);
1625   filtered_aln_file=vtmpnam(NULL);
1626
1627
1628
1629   output_clustal_aln (aln_file, A);
1630   /* 1: remove columns with too many gaps*/
1631   printf_system ("t_coffee -other_pg seq_reformat -in %s -action +rm_gap %d -output clustalw > %s", aln_file,fraction_gap, ungaped_aln_file);
1632
1633   /* 2: evaluate the alignment*/
1634
1635   printf_system ("t_coffee -other_pg seq_reformat -in %s -action +evaluate %s -output clustalw > %s", ungaped_aln_file,(mode)?mode:"categories", scored_aln_file);
1636
1637
1638   /*3 extract the high scoring columns*/
1639   printf_system("t_coffee -other_pg seq_reformat -in %s -struc_in %s -struc_in_f number_aln -action +use_cons +keep '[%d-8]' +rm_gap -output clustalw > %s", ungaped_aln_file, scored_aln_file,n, filtered_aln_file);
1640
1641
1642   free_aln (A);
1643
1644   A=main_read_aln ( filtered_aln_file, NULL);
1645   print_aln (A);
1646
1647   return A;
1648 }
1649
1650 NT_node tree2fj_tree (NT_node T)
1651 {
1652   NT_node L;
1653
1654   return T;
1655
1656   L=find_longest_branch (T, NULL);
1657   T=reroot_tree (T, L);
1658   return T;
1659 }
1660
1661
1662 /*********************************************************************/
1663 /*                                                                   */
1664 /*                                   Tree Filters and MAnipulation   */
1665 /*                                                                   */
1666 /*                                                                   */
1667 /*********************************************************************/
1668 int tree2star_nodes (NT_node R, int n_max)
1669 {
1670   if ( !R) return 0;
1671   else if (!R->left && !R->right)
1672     {
1673       if (n_max>=1)R->dist=0;
1674       return 1;
1675     }
1676   else
1677     {
1678       int n=0;
1679
1680       n+=tree2star_nodes (R->right, n_max);
1681       n+=tree2star_nodes (R->left, n_max);
1682
1683       if (n<n_max)R->dist=0;
1684       return n;
1685     }
1686 }
1687
1688 NT_node aln2tree (Alignment *A)
1689 {
1690   NT_node **T=NULL;
1691
1692
1693   T=make_nj_tree (A, NULL, 0, 0, A->seq_al, A->name, A->nseq, NULL, NULL);
1694   tree2nleaf (T[3][0]);
1695
1696   return T[3][0];
1697 }
1698 NT_node realloc_tree ( NT_node R, int n)
1699 {
1700   if ( !R)return R;
1701   R->right=realloc_tree (R->right,n);
1702   R->left=realloc_tree (R->left,n);
1703   R->bot=realloc_tree (R->bot,n);
1704
1705   R->lseq=vrealloc (R->lseq, n*sizeof (int));
1706   R->lseq2=vrealloc (R->lseq2, n*sizeof (int));
1707   return R;
1708 }
1709
1710 NT_node reset_boot_tree ( NT_node R, int n)
1711 {
1712   if ( !R)return R;
1713   R->right=reset_boot_tree (R->right,n);
1714   R->left=reset_boot_tree (R->left,n);
1715   R->bot=reset_boot_tree (R->bot,n);
1716   R->bootstrap=(float)n;
1717
1718   return R;
1719 }
1720 NT_node tree_dist2normalized_tree_dist ( NT_node R, float max)
1721 {
1722   if (!R)return R;
1723   else
1724     {
1725       tree_dist2normalized_tree_dist ( R->right, max);
1726       tree_dist2normalized_tree_dist ( R->left, max);
1727       R->bootstrap=(int)((R->dist*100)/max);
1728     }
1729   return R;
1730 }
1731 NT_node reset_dist_tree ( NT_node R, float n)
1732 {
1733   if ( !R)return R;
1734   R->right=reset_dist_tree (R->right,n);
1735   R->left=reset_dist_tree (R->left,n);
1736   R->bot=reset_dist_tree (R->bot,n);
1737
1738   if (R->parent && !(R->parent)->parent && !(R->parent)->bot)R->dist=n/2;
1739   else R->dist=n;
1740
1741   return R;
1742 }
1743
1744
1745 NT_node* free_treelist (NT_node *L)
1746 {
1747   int n=0;
1748   while (L[n])free_tree (L[n++]);
1749   vfree (L);
1750   return NULL;
1751 }
1752 NT_node free_tree ( NT_node R)
1753 {
1754   if ( !R)return R;
1755   R->right=free_tree (R->right);
1756   R->left=free_tree (R->left);
1757   R->bot=free_tree (R->bot);
1758   free_tree_node (R);
1759   return R;
1760 }
1761
1762 NT_node free_tree_node ( NT_node R)
1763 {
1764   if (!R)return NULL;
1765
1766   vfree (R->seqal);
1767   vfree (R->idist);
1768   vfree (R->ldist);
1769   vfree (R->file);
1770   vfree ( R->name);
1771   vfree ( R->lseq); vfree ( R->lseq2);
1772   vfree (R);
1773   return NULL;
1774 }
1775
1776 NT_node   rename_seq_in_tree ( NT_node R, char ***list)
1777 {
1778   if ( !R || !list) return R;
1779
1780   if ( R->leaf!=1)
1781     {
1782       R->right=rename_seq_in_tree (R->right, list);
1783       R->left=rename_seq_in_tree (R->left, list);
1784       R->bot=rename_seq_in_tree (R->bot, list);
1785     }
1786   else
1787     {
1788       int n=0;
1789       while ( list[n][0][0])
1790         {
1791           if ( strm (list[n][0], R->name))sprintf (R->name, "%s",list[n][1]);
1792           n++;
1793         }
1794     }
1795   return R;
1796 }
1797 Sequence * tree2seq    (NT_node R, Sequence *S)
1798 {
1799
1800   if ( !R)return S;
1801   if ( !S)
1802     {
1803       S=declare_sequence (10, 10, tree2nseq (R));
1804       S->nseq=0;
1805     }
1806
1807   if (R->leaf==1)
1808     {
1809       sprintf ( S->name[S->nseq++], "%s", R->name);
1810     }
1811   else
1812     {
1813       S=tree2seq (R->left, S);
1814       S=tree2seq (R->right, S);
1815     }
1816   return S;
1817 }
1818
1819 NT_node balance_tree (NT_node T)
1820 {
1821   static int **list;
1822   NT_node NL[3];
1823
1824   if ( !T) return T;
1825   else if ( T->leaf<=2)return T;
1826   else
1827     {
1828       if (!list)list=declare_int (3, 2);
1829
1830       NL[0]=T->left;
1831       NL[1]=T->right;
1832       NL[2]=T->bot;
1833
1834       list[0][0]=(T->left)?(T->left)->leaf:0;
1835       list[0][1]=0;
1836       list[1][0]=(T->right)?(T->right)->leaf:0;
1837       list[1][1]=1;
1838       list[2][0]=(T->bot)?(T->bot)->leaf:0;
1839       list[2][1]=2;
1840
1841       sort_int (list,2,0,0,2);
1842
1843       T->left=NL[list[2][1]];
1844       T->right=NL[list[1][1]];
1845       T->bot=NL[list[0][1]];
1846
1847       T->left=balance_tree (T->left);
1848       T->right=balance_tree (T->right);
1849       T->bot=balance_tree (T->bot);
1850       return T;
1851     }
1852 }
1853 FILE * display_tree (NT_node R, int nseq, FILE *fp)
1854 {
1855   int a;
1856
1857   if ( !R);
1858   else
1859     {
1860       /*
1861         if ( R->nseq==1)fprintf (stderr,"\n[%s] ", R->name);
1862         else fprintf ( stderr, "\n[%d Node] ",R->nseq);
1863         for ( a=0; a< R->nseq; a++) fprintf ( stderr, "[%d]", R->lseq[a]);
1864       */
1865       fprintf (fp, "\n %10s N ", R->name);
1866       for ( a=0; a< nseq; a++)fprintf (fp, "%d", R->lseq2[a]);
1867       fprintf (fp, "\n %10s D ", R->name);
1868       for ( a=0; a< nseq; a++)fprintf (fp, "%d", R->idist[a]);
1869
1870
1871       if (R->leaf==1) fprintf (fp, " %s", R->name);
1872       fprintf (fp, " :%.4f", R->dist);
1873       HERE ("\nGo Left");fp=display_tree (R->left, nseq, fp);
1874       HERE ("\nGo Right");fp=display_tree (R->right, nseq, fp);
1875       HERE ("\nGo Bot");fp=display_tree (R->bot, nseq, fp);
1876     }
1877   return fp;
1878 }
1879 int tree2nnode_unresolved (NT_node R, int *l)
1880 {
1881   if ( !R)return 0;
1882   else if (R->leaf && R->dist==0){return 1;}
1883   else
1884     {
1885       int n=0;
1886       n+=tree2nnode_unresolved (R->right, l);
1887       n+=tree2nnode_unresolved (R->left, l);
1888       if (R->dist==0)
1889         {
1890           return n;
1891         }
1892       else
1893         {
1894           if (n)l[n]++;
1895           return 0;
1896         }
1897     }
1898
1899 }
1900
1901 int tree2nnode ( NT_node R)
1902 {
1903   int n;
1904   if ( !R)n=0;
1905   else if ( R->leaf==1){R->node=1;n=1;}
1906   else
1907     {
1908       n=1;
1909       n+=tree2nnode (R->right);
1910       n+=tree2nnode (R->left);
1911       n+=tree2nnode (R->bot);
1912       R->node=n;
1913     }
1914   return n;
1915 }
1916 int tree2nleaf (NT_node R)
1917 {
1918   if ( !R)return 0;
1919   else if (R->leaf==1){return 1;}
1920   else if (R->right==NULL && R->left==NULL && R->bot==NULL){R->leaf=1; return 1;}
1921   else
1922     {
1923       int n=0;
1924       n+=tree2nleaf (R->right);
1925       n+=tree2nleaf (R->left);
1926       n+=tree2nleaf (R->bot);
1927
1928       R->leaf=n;
1929       return n;
1930     }
1931 }
1932
1933 int tree2nseq ( NT_node R)
1934 {
1935  return tree2nleaf(R);
1936 }
1937
1938 int tree_file2nseq (char *fname)
1939 {
1940   FILE *fp;
1941   char *string;
1942   int p, a, b, c, n;
1943
1944   string=vcalloc (count_n_char_in_file(fname)+1, sizeof (char));
1945
1946   fp=vfopen (fname, "r");
1947   n=0;
1948   while ( (c=fgetc(fp))!=EOF){if (c=='(' || c==')' || c==',' || c==';') string[n++]=c;}
1949   vfclose (fp);string[n]='\0';
1950
1951   for (n=0, p=1; p<strlen (string); p++)
1952     {
1953       a=string[p-1];
1954       b=string[p];
1955
1956       if ( a=='(' && b==',')n++;
1957       else if ( a==',' && b==')')n++;
1958       else if ( a==',' && b==',')n++;
1959     }
1960   if ( getenv4debug("DEBUG_TREE"))fprintf (stderr, "\n[DEBUG_TREE:tree_file2nseq]%s", string);
1961   vfree (string);
1962   return n;
1963 }
1964
1965
1966 void clear_tree ( NT_node R)
1967 {
1968   if (!R) return;
1969
1970   R->visited=0;
1971
1972   if ( R->leaf==1);
1973   else
1974     {
1975       clear_tree ( R->right);
1976       clear_tree ( R->left);
1977       clear_tree ( R->bot);
1978     }
1979 }
1980 int display_leaf_below_node (NT_node T, FILE *fp)
1981 {
1982   int n=0;
1983   if ( !T)return 0;
1984
1985   if ( T->leaf==1)
1986     {
1987       fprintf (fp, " %s", T->name);
1988       return 1;
1989     }
1990   else
1991     {
1992       n+=display_leaf_below_node ( T->right, fp);
1993       n+=display_leaf_below_node ( T->left, fp);
1994       return n;
1995     }
1996 }
1997 int display_leaf ( NT_node T, FILE *fp)
1998 {
1999   int n=0;
2000   if ( !T)return 0;
2001   else if ( T->visited)return 0;
2002   else T->visited=1;
2003
2004   if ( T->leaf==1)
2005     {
2006       fprintf (fp, " %s", T->name);
2007       return 1;
2008     }
2009   else
2010     {
2011       n+=display_leaf ( T->right, fp);
2012       n+=display_leaf ( T->left, fp);
2013       n+=display_leaf ( T->bot, fp);
2014       return n;
2015     }
2016 }
2017
2018
2019
2020
2021 NT_node find_longest_branch ( NT_node T, NT_node L)
2022   {
2023
2024     if ( !L || T->dist>L->dist)
2025       {
2026
2027         L=T;
2028       }
2029
2030     if ( T->leaf==1)return L;
2031     else
2032       {
2033         L=find_longest_branch ( T->right, L);
2034         L=find_longest_branch ( T->left,  L);
2035         return L;
2036       }
2037   }
2038 int node2side (NT_node N);
2039 int test_print (NT_node T);
2040 NT_node straighten_node (NT_node N);
2041 NT_node EMPTY;
2042 NT_node Previous;
2043 NT_node reroot_tree ( NT_node TREE, NT_node Right)
2044 {
2045   /*ReRoots the tree between Node R and its parent*/
2046   NT_node NR;
2047   int n1, n2;
2048
2049   if (!EMPTY)EMPTY=vcalloc (1, sizeof (NT_node));
2050   if ( !Right->parent)return Right;
2051
2052   TREE=unroot_tree (TREE);
2053   if (Right->parent==NULL && Right->bot)
2054     Right=Right->bot;
2055
2056   n1=tree2nleaf (TREE);
2057
2058   NR=declare_tree_node(TREE->maxnseq);
2059
2060   NR->right=Right;
2061   NR->left=Right->parent;
2062   Right->parent=NR;
2063
2064   Right->dist=Right->dist/2;
2065
2066   if ((NR->left)->right==Right)(NR->left)->right=EMPTY;
2067   else if ( (NR->left)->left==Right) (NR->left)->left=EMPTY;
2068
2069   Previous=NULL;
2070
2071
2072   NR->left=straighten_node (NR->left);
2073
2074
2075
2076   (NR->left)->parent=NR;
2077   (NR->left)->dist=Right->dist;
2078
2079
2080
2081   n2=tree2nleaf(NR);
2082
2083   if ( n1!=n2){fprintf ( stderr, "\n%d %d", n1, n2);myexit (EXIT_FAILURE);}
2084   return NR;
2085 }
2086
2087 NT_node straighten_node ( NT_node N)
2088 {
2089   NT_node Child;
2090
2091
2092   if ( N->parent)
2093     {
2094       if (N->right==EMPTY)N->right=N->parent;
2095       else if ( N->left==EMPTY) N->left=N->parent;
2096
2097       Child=N->parent;
2098       if (Child->right==N)
2099         {
2100           Child->right=EMPTY;
2101         }
2102       else if (Child->left==N)
2103         {
2104           Child->left=EMPTY;
2105         }
2106
2107       Previous=N;
2108       Child=straighten_node (Child);
2109       Child->parent=N;
2110       Child->dist=N->dist;
2111       return N;
2112     }
2113   else if ( N->bot && N->bot!=Previous)
2114     {
2115       if ( N->right==EMPTY)N->right=N->bot;
2116       else if ( N->left==EMPTY)N->left=N->bot;
2117
2118       N->bot=NULL;
2119       return N;
2120     }
2121   else
2122     {
2123       N->bot=NULL;
2124       return N;
2125     }
2126 }
2127 int test_print (NT_node T)
2128 {
2129   if ( !T)
2130     {
2131       fprintf ( stderr, "\nEMPTY");
2132     }
2133   else if ( !T->left && !T->right)
2134     {
2135       fprintf ( stderr, "\n%s",T->name);
2136     }
2137   else
2138     {
2139       fprintf ( stderr, "\nGoing Right");
2140       test_print (T->right);
2141       fprintf ( stderr, "\nGoing Left");
2142       test_print (T->left);
2143     }
2144   return 1;
2145 }
2146 int node2side (NT_node C)
2147 {
2148   if ( !C->parent) return UNKNOWN;
2149   else if ( (C->parent)->left==C)return LEFT;
2150   else if ( (C->parent)->right==C)return RIGHT;
2151   else return UNKNOWN;
2152 }
2153 NT_node straighten_tree ( NT_node P, NT_node C, float new_dist)
2154 {
2155   float dist;
2156
2157   if ( C==NULL)return NULL;
2158
2159
2160   dist=C->dist;
2161   C->dist=new_dist;
2162   C->bot=NULL;
2163
2164   if (C->left && C->right)
2165     {
2166       C->parent=P;
2167     }
2168   else if (!C->left)
2169     {
2170       C->left=C->parent;
2171       C->parent=P;
2172     }
2173
2174   if ( C->parent==P);
2175   else if ( C->left==NULL && C->right==NULL)
2176     {
2177       C->parent=P;
2178     }
2179   else if ( C->right==P)
2180     {
2181       C->right=C->parent;
2182       C->parent=P;
2183
2184       C=straighten_tree(C, C->right, dist);
2185     }
2186   else if ( C->left==P)
2187     {
2188       C->left=C->parent;
2189       C->parent=P;
2190       C=straighten_tree (C, C->left, dist);
2191     }
2192   else if ( C->parent==NULL)
2193     {
2194       C->parent=P;
2195     }
2196
2197   return C;
2198 }
2199
2200
2201 NT_node unroot_tree ( NT_node T)
2202 {
2203
2204   if (!T || T->visited) return T;
2205   else T->visited=1;
2206
2207   if (T->parent==NULL)
2208     {
2209
2210       (T->right)->dist=(T->left)->dist=(T->right)->dist+(T->left)->dist;
2211       (T->right)->parent=T->left;
2212       (T->left)->parent=T->right;
2213       T=T->left;
2214       T->leaf=0;
2215       vfree (T->parent);
2216     }
2217   else
2218     {
2219       T->parent=unroot_tree (T->parent);
2220       T->right=unroot_tree (T->right);
2221       T->left=unroot_tree (T->left);
2222     }
2223   T->visited=0;
2224   return T;
2225 }
2226
2227 FILE * print_tree_list ( NT_node *T, char *format,FILE *fp)
2228 {
2229   int a=0;
2230   while ( T[a])
2231     {
2232       fp=print_tree (T[a], format, fp);
2233       a++;
2234     }
2235   return fp;
2236 }
2237 char * tree2string (NT_node T)
2238 {
2239   if (!T) return NULL;
2240   else
2241     {
2242       static char *f;
2243       FILE *fp;
2244
2245       if (!f)f=vtmpnam (NULL);
2246       fp=vfopen (f, "w");
2247       print_tree (T, "newick", fp);
2248       vfclose (fp);
2249       return file2string (f);
2250     }
2251 }
2252 char * tree2file (NT_node T, char *name, char *mode)
2253 {
2254   if (!name)name=vtmpnam (NULL);
2255   string2file (name, mode, tree2string(T));
2256   return name;
2257 }
2258 FILE * print_tree ( NT_node T, char *format,FILE *fp)
2259 {
2260   Sequence *S;
2261
2262   tree2nleaf(T);
2263   S=tree2seq(T, NULL);
2264
2265   recode_tree (T, S);
2266
2267   free_sequence (S, -1);
2268   if ( format && strm (format, "binary"))
2269     fp=display_tree ( T,S->nseq, fp);
2270   else if ( ! format || strm2 (format, "newick_tree","newick"))
2271     {
2272       /*T=balance_tree (T);*/
2273       fp=rec_print_tree (T, fp);
2274       fprintf ( fp, ";\n");
2275     }
2276   else
2277     {
2278       fprintf ( stderr, "\nERROR: %s is an unknown tree format [FATAL:%s]\n", format, PROGRAM);
2279       myexit (EXIT_FAILURE);
2280     }
2281   return fp;
2282 }
2283 int print_newick_tree ( NT_node T, char *name)
2284 {
2285   FILE *fp;
2286   fp=vfopen (name, "w");
2287   fp=rec_print_tree (T,fp);
2288   fprintf (fp, ";\n");
2289   vfclose (fp);
2290   return 1;
2291 }
2292 FILE * rec_print_tree ( NT_node T, FILE *fp)
2293 {
2294
2295
2296
2297   if (!T)return fp;
2298
2299   if ( T->isseq)
2300     {
2301       fprintf ( fp, " %s:%.5f",T->name, T->dist);
2302     }
2303   else
2304     {
2305       if (T->left && T->right)
2306         {
2307           fprintf ( fp, "(");fp=rec_print_tree ( T->left, fp);
2308           fprintf ( fp, ",");fp=rec_print_tree ( T->right, fp);
2309           fprintf ( fp, ")");
2310           if (T->parent || T->dist)
2311             {
2312               if ( T->bootstrap!=0)fprintf (fp, " %d", (int)T->bootstrap);
2313               fprintf (fp, ":%.5f", T->dist);
2314             }
2315         }
2316       else if (T->left)fp=rec_print_tree (T->left, fp);
2317       else if (T->right)fp=rec_print_tree(T->right, fp);
2318     }
2319
2320   return fp;
2321 }
2322
2323
2324
2325
2326
2327 /*********************************************************************/
2328 /*                                                                   */
2329 /*                                  Tree Functions                   */
2330 /*                                                                   */
2331 /*                                                                   */
2332 /*********************************************************************/
2333
2334 int ** make_sub_tree_list ( NT_node **T, int nseq, int n_node)
2335         {
2336
2337
2338 /*This function produces a list of all the sub trees*/
2339
2340
2341 /*        /A */
2342 /*      -* */
2343 /*        \      /B */
2344 /*         \    / */
2345 /*          ---* */
2346 /*                  \ */
2347 /*                   *--C */
2348 /*                    \ */
2349 /*                     \D */
2350
2351 /*      Contains 4 i_nodes */
2352 /*             8   nodes (internal nodes +leaves) */
2353 /*             8 sub trees: */
2354 /*             ABCD */
2355 /*                 1111 */
2356 /*             0111 */
2357 /*             1000 */
2358 /*             0100 */
2359 /*             0011 */
2360 /*             0001 */
2361 /*             0010 */
2362
2363         int **sub_tree_list;
2364         int a, n=0;
2365
2366
2367         if (T)
2368              {
2369                  sub_tree_list=declare_int ( (n_node), nseq);
2370                  make_all_sub_tree_list (T[3][0],sub_tree_list, &n);
2371
2372              }
2373         else
2374              {
2375                  sub_tree_list=declare_int (nseq, nseq);
2376                  for ( a=0; a< nseq; a++)sub_tree_list[a][a]=1;
2377              }
2378
2379         return sub_tree_list;
2380         }
2381
2382 void make_all_sub_tree_list ( NT_node N, int **list, int *n)
2383         {
2384         make_one_sub_tree_list (N, list[n[0]++]);
2385         if (N->leaf!=1)
2386             {
2387             make_all_sub_tree_list (N->left , list, n);
2388             make_all_sub_tree_list (N->right, list, n);
2389             }
2390         return;
2391         }
2392
2393 void make_one_sub_tree_list ( NT_node T,int *list)
2394         {
2395         if (T->leaf==1)
2396             {
2397
2398             list[T->seq]=1;
2399             }
2400         else
2401             {
2402             make_one_sub_tree_list(T->left , list);
2403             make_one_sub_tree_list(T->right, list);
2404             }
2405         return;
2406         }
2407
2408
2409 NT_node old_main_read_tree(char *treefile)
2410 {
2411   /*Reads a tree w/o needing the sequence file*/
2412   NT_node **T;
2413   T=simple_read_tree (treefile);
2414   return T[3][0];
2415 }
2416
2417
2418
2419 NT_node** simple_read_tree(char *treefile)
2420 {
2421   int tot_node=0;
2422   NT_node **T;
2423   T=read_tree ( treefile, &tot_node,tree_file2nseq (treefile),NULL);
2424   return T;
2425 }
2426
2427 void free_read_tree ( NT_node **BT)
2428 {
2429   int a, s;
2430
2431
2432   if (!BT) return;
2433
2434   for (s=0,a=0; a<3; a++)
2435     {
2436       vfree (BT[a]);
2437     }
2438   free_tree (BT[3][0]);
2439   vfree (BT);
2440   return;
2441 }
2442
2443 NT_node** read_tree(char *treefile, int *tot_node,int nseq, char **seq_names)
2444         {
2445
2446             /*The Tree Root is in the TREE[3][0]...*/
2447             /*TREE[0][ntot]--> pointer to each node and leave*/
2448         char ch;
2449         int  a,b;
2450
2451         FILE *fp;
2452         int nseq_read = 0;
2453         int nnodes = 0;/*Number of Internal Nodes*/
2454         int ntotal = 0;/*Number of Internal Nodes + Number of Leaves*/
2455         int flag;
2456         int c_seq;
2457         NT_node **lu_ptr;
2458         NT_node seq_tree, root,p;
2459
2460
2461         tot_nseq=nseq;
2462         rooted_tree=distance_tree=TRUE;
2463
2464         fp = vfopen(treefile, "r");
2465         fp=skip_space(fp);
2466         ch = (char)getc(fp);
2467         if (ch != '(')
2468                 {
2469                 fprintf(stderr, "Error: Wrong format in tree file %s\n", treefile);
2470                 myexit (EXIT_FAILURE);
2471                 }
2472         rewind(fp);
2473
2474
2475         lu_ptr=(NT_node **)vcalloc(4,sizeof(NT_node*));
2476         lu_ptr[0] = (NT_node *)vcalloc(10*nseq,sizeof(NT_node));
2477         lu_ptr[1] = (NT_node *)vcalloc(10*nseq,sizeof(NT_node));
2478         lu_ptr[2] = (NT_node *)vcalloc(10*nseq,sizeof(NT_node));
2479         lu_ptr[3] =(NT_node *) vcalloc(1,sizeof(NT_node));
2480
2481         seq_tree =(NT_node) declare_tree_node(nseq);
2482
2483         set_info(seq_tree, NULL, 0, "  ", 0.0, 0);
2484
2485
2486         fp=create_tree(seq_tree,NULL,&nseq_read, &ntotal, &nnodes, lu_ptr, fp);
2487         fclose (fp);
2488
2489
2490         if (nseq != tot_nseq)
2491                 {
2492                 fprintf(stderr," Error: tree not compatible with alignment (%d sequences in alignment and %d in tree\n", nseq,nseq_read);
2493                 myexit (EXIT_FAILURE);
2494                 }
2495
2496         if (distance_tree == FALSE)
2497                 {
2498                 if (rooted_tree == FALSE)
2499                         {
2500                         fprintf(stderr,"Error: input tree is unrooted and has no distances, cannot align sequences\n");
2501                         myexit (EXIT_FAILURE);
2502                         }
2503                 }
2504
2505         if (rooted_tree == FALSE)
2506                 {
2507                   root = reroot(seq_tree, nseq,ntotal,nnodes, lu_ptr);
2508                   lu_ptr[1][nnodes++]=lu_ptr[0][ntotal++]=root;
2509
2510                 }
2511         else
2512                 {
2513                 root = seq_tree;
2514                 }
2515
2516         lu_ptr[3][0]=root;
2517         tot_node[0]=nnodes;
2518
2519
2520
2521         for ( a=0; a< ntotal; a++)
2522                 {
2523                 (lu_ptr[0][a])->isseq=(lu_ptr[0][a])->leaf;
2524                 (lu_ptr[0][a])->dp=(lu_ptr[0][a])->dist;
2525                 }
2526
2527
2528         for ( a=0; a< nseq; a++)
2529                 {
2530                 if (!seq_names)
2531                   {
2532                     flag=1;
2533                     (lu_ptr[2][a])->order=(lu_ptr[2][a])->seq=a;
2534                   }
2535                 else
2536                   {
2537                     for ( flag=0,b=0; b<nseq; b++)
2538                       {
2539                         if ( strncmp ( (lu_ptr[2][a])->name, seq_names[b], MAXNAMES)==0)
2540                           {
2541                             flag=1;
2542
2543                             (lu_ptr[2][a])->order=(lu_ptr[2][a])->seq=b;
2544                                 /*vfree ( (lu_ptr[2][a])->name);*/
2545                             sprintf ((lu_ptr[2][a])->name, "%s", seq_names[b]);
2546                           }
2547                       }
2548                   }
2549                 /*
2550                 if ( flag==0  && (lu_ptr[0][a])->leaf==1)
2551                       {
2552                         fprintf ( stderr, "\n%s* not in tree",(lu_ptr[2][a])->name);
2553                         for ( a=0; a< ntotal; a++)
2554                           {
2555                             fprintf ( stderr, "\n%d %s",(lu_ptr[2][a])->leaf, (lu_ptr[2][a])->name);
2556                           }
2557                       }
2558                 */
2559                 }
2560
2561         if (seq_names)
2562           {
2563             int tnseq;
2564             char *s;
2565             char **tree_names;
2566             int fail_flag=0;
2567             tnseq=tree_file2nseq(treefile);
2568             tree_names=vcalloc ( tnseq, sizeof (char*));
2569             for (a=0; a<tnseq; a++)
2570               {
2571                 s=(lu_ptr[2][a])->name;
2572                 tree_names[a]=s;
2573                 if ( name_is_in_list(s, seq_names, nseq, MAXNAMES+1)==-1)
2574                   {
2575                     fprintf (stderr, "\nERROR: Sequence %s in the tree [%s] is not in the alignment[FATAL:%s]\n", s, treefile, PROGRAM);
2576                     fail_flag=1;
2577                   }
2578               }
2579             for (a=0; a<nseq; a++)
2580               {
2581                 s=seq_names[a];
2582                 if ( name_is_in_list(s, tree_names, nseq, MAXNAMES+1)==-1)
2583                   {
2584                     fprintf (stderr, "\nERROR: Sequence %s in the sequences is not in the tree [%s][FATAL:%s]\n", s, treefile, PROGRAM);
2585                     fail_flag=1;
2586                   }
2587               }
2588             vfree (tree_names);
2589             if ( fail_flag==1)myexit (EXIT_FAILURE);
2590           }
2591
2592         for ( a=0; a< nseq; a++)
2593                 {
2594                 p=lu_ptr[2][a];
2595                 c_seq=p->seq;
2596
2597                 while ( p!=NULL)
2598                         {
2599                         p->lseq[p->nseq]=c_seq;
2600                         p->nseq++;
2601                         p=p->parent;
2602                         }
2603                 }
2604
2605
2606         return lu_ptr;
2607         }
2608
2609 FILE * create_linear_tree ( char **name, int n, FILE *fp)
2610 {
2611
2612   if (!name || n==0 ||!fp) return NULL;
2613
2614
2615   if (n==2)
2616     fprintf ( fp, "(%s,%s);",name[0],name[1]);
2617   else if ( n==3)
2618     fprintf ( fp, "((%s,%s),%s);",name[0],name[1], name[2]);
2619   else
2620     {
2621       int a;
2622       for (a=0; a<n-2; a++)fprintf (fp, "(");
2623       fprintf (fp, "%s, %s),", name[0], name[1]);
2624       for ( a=2; a<n-2; a++)fprintf ( fp, "%s),",name[a]);
2625       fprintf ( fp, "%s,%s);\n", name[n-2], name[n-1]);
2626     }
2627   return fp;
2628 }
2629 FILE * create_tree(NT_node ptree, NT_node parent,int *nseq,int  *ntotal,int  *nnodes,NT_node **lu, FILE *fp)
2630         {
2631
2632         int i, type;
2633         float dist=0;
2634         float bootstrap=0;
2635         char *name;
2636         int ch;
2637
2638
2639         name=vcalloc ( MAXNAMES+1, sizeof (char));
2640         sprintf ( name, "   ");
2641         fp=skip_space(fp);
2642         ch = (char)getc(fp);
2643
2644         if (ch == '(')
2645                 {
2646                 type = NODE;
2647                 name[0] = '\0';
2648                 lu[0][ntotal[0]] = lu[1][nnodes[0]] = ptree;
2649                 ntotal[0]++;
2650                 nnodes[0]++;
2651                 create_tree_node(ptree, parent);
2652                 fp=create_tree(ptree->left, ptree, nseq,ntotal,nnodes,lu,fp);
2653                 ch = (char)getc(fp);
2654                 if ( ch == ',')
2655                         {
2656                         fp=create_tree(ptree->right, ptree,nseq,ntotal,nnodes,lu,fp);
2657                         ch = (char)getc(fp);
2658                         if ( ch == ',')
2659                                 {
2660
2661                                 ptree = insert_tree_node(ptree);
2662                                 lu[0][ntotal[0]] = lu[1][nnodes[0]] = ptree;
2663                                 ntotal[0]++;
2664                                 nnodes[0]++;
2665                                 fp=create_tree(ptree->right, ptree,nseq,ntotal,nnodes,lu,fp);
2666                                 rooted_tree = FALSE;
2667                                 if ( getenv4debug ( "DEBUG_TREE")){fprintf ( stderr, "\n[DEBUG_TREE:create_tree] Unrooted Tree");}
2668                                 }
2669                         }
2670
2671                 fp=skip_space(fp);
2672                 ch = (char)getc(fp);
2673                 }
2674         else
2675                 {
2676                 type=LEAF;
2677                 lu[0][ntotal[0]] = lu[2][nseq[0]] = ptree;
2678                 ntotal[0]++;
2679                 nseq[0]++;
2680                 name[0] = ch;
2681                 i=1;
2682                 ch = (char)getc(fp);
2683                 if ( name[0]=='\'')
2684                   {
2685                     /*This protects names that are between single quotes*/
2686                     while ( ch!='\'')
2687                       {
2688                         if (i < MAXNAMES) name[i++] = ch;
2689                         ch = (char)getc(fp);
2690                       }
2691                     if (i < MAXNAMES) name[i++] = ch;
2692                     while ((ch != ':') && (ch != ',') && (ch != ')'))ch = (char)getc(fp);
2693                   }
2694                 else
2695                   {
2696                     while ((ch != ':') && (ch != ',') && (ch != ')'))
2697                       {
2698                         if (i < MAXNAMES) name[i++] = ch;
2699                         ch = (char)getc(fp);
2700                       }
2701                   }
2702
2703                 name[i] = '\0';
2704
2705                 if ( i>=(MAXNAMES+1)){fprintf (stderr, "\nName is too long");myexit (EXIT_FAILURE);}
2706                 if (ch != ':' && !isdigit(ch))
2707                         {
2708                           /*distance_tree = FALSE*/;
2709                         }
2710                 }
2711         if (ch == ':')
2712                 {
2713                 fp=skip_space(fp);
2714                 fscanf(fp,"%f",&dist);
2715                 fp=skip_space(fp);
2716                 bootstrap=0;
2717                 }
2718         /*Tree with Bootstrap information*/
2719         else if (isdigit (ch))
2720           {
2721             ungetc(ch,fp);
2722             fscanf(fp,"%f",&bootstrap);
2723             if ( fscanf(fp,":%f",&dist)==1);
2724             else dist=0;
2725             fp=skip_space(fp);
2726           }
2727         else
2728           {
2729             ungetc ( ch, fp);
2730             skip_space(fp);
2731           }
2732
2733         set_info(ptree, parent, type, name, dist, bootstrap);
2734
2735
2736         vfree (name);
2737         return fp;
2738         }
2739
2740 NT_node declare_tree_node (int nseq)
2741         {
2742         NT_node p;
2743
2744         p= (NT_node)vcalloc (1, sizeof ( Treenode));
2745         p->left = NULL;
2746         p->right = NULL;
2747         p->parent = NULL;
2748         p->dist = 0.0;
2749         p->leaf = 0;
2750         p->order = 0;
2751         p->maxnseq=nseq;
2752         p->name=(char*)vcalloc (MAXNAMES+1,sizeof (char));
2753         p->name[0]='\0';
2754         p->lseq=(int*)vcalloc ( nseq, sizeof (int));
2755         return p;
2756
2757         }
2758
2759 void set_info(NT_node p, NT_node parent, int pleaf, char *pname, float pdist, float bootstrap)
2760         {
2761         p->parent = parent;
2762         p->leaf = pleaf;
2763         p->dist = pdist;
2764         p->bootstrap=bootstrap;
2765         p->order = 0;
2766
2767
2768         sprintf (p->name, "%s", pname);
2769
2770         if (pleaf ==1)
2771                 {
2772                 p->left = NULL;
2773                 p->right = NULL;
2774                 }
2775
2776    }
2777 NT_node insert_tree_node(NT_node pptr)
2778         {
2779
2780         NT_node newnode;
2781
2782         newnode = declare_tree_node( pptr->maxnseq);
2783         create_tree_node(newnode, pptr->parent);
2784
2785         newnode->left = pptr;
2786         pptr->parent = newnode;
2787
2788         set_info(newnode, pptr->parent, 0, "", 0.0, 0);
2789
2790         return(newnode);
2791         }
2792
2793 void create_tree_node(NT_node pptr, NT_node parent)
2794         {
2795         pptr->parent = parent;
2796         pptr->left =declare_tree_node(pptr->maxnseq) ;
2797         (pptr->left)->parent=pptr;
2798
2799         pptr->right =declare_tree_node(pptr->maxnseq) ;
2800         (pptr->right)->parent=pptr;
2801         }
2802
2803 FILE * skip_space(FILE *fp)
2804         {
2805         int   c;
2806
2807         do
2808         c = getc(fp);
2809         while(isspace(c));
2810         if ( c==EOF)
2811                 {
2812                 fprintf ( stderr, "\nEOF");
2813                 myexit (EXIT_FAILURE);
2814                 }
2815         ungetc(c, fp);
2816         return fp;
2817         }
2818
2819
2820 NT_node reroot(NT_node ptree, int nseq, int ntotal, int nnodes, NT_node **lu)
2821         {
2822         NT_node p, rootnode, rootptr;
2823         float   diff, mindiff=0, mindepth = 1.0, maxdist;
2824         int   i;
2825         int first = TRUE;
2826
2827
2828
2829         rootptr = ptree;
2830
2831         for (i=0; i<ntotal; i++)
2832                 {
2833                 p = lu[0][i];
2834                 if (p->parent == NULL)
2835                         diff = calc_root_mean(p, &maxdist, nseq, lu);
2836                 else
2837                         diff = calc_mean(p, &maxdist, nseq, lu);
2838
2839                 if ((diff == 0) || ((diff > 0) && (diff < 2 * p->dist)))
2840                          {
2841                         if ((maxdist < mindepth) || (first == TRUE))
2842                                 {
2843                                 first = FALSE;
2844                                 rootptr = p;
2845                                 mindepth = maxdist;
2846                                 mindiff = diff;
2847                                 }
2848                         }
2849
2850                 }
2851         if (rootptr == ptree)
2852                 {
2853                 mindiff = rootptr->left->dist + rootptr->right->dist;
2854                 rootptr = rootptr->right;
2855                 }
2856
2857         rootnode = insert_root(rootptr, mindiff);
2858         diff = calc_root_mean(rootnode, &maxdist, nseq, lu);
2859         return(rootnode);
2860         }
2861
2862
2863 float calc_root_mean(NT_node root, float *maxdist, int nseq, NT_node **lu)
2864         {
2865         float dist , lsum = 0.0, rsum = 0.0, lmean,rmean,diff;
2866         NT_node p;
2867         int i;
2868         int nl, nr;
2869         int direction;
2870
2871
2872         dist = (*maxdist) = 0;
2873         nl = nr = 0;
2874         for (i=0; i< nseq; i++)
2875           {
2876           p = lu[2][i];
2877           dist = 0.0;
2878           while (p->parent != root)
2879                 {
2880                 dist += p->dist;
2881                 p = p->parent;
2882                 }
2883          if (p == root->left) direction = LEFT;
2884          else direction = RIGHT;
2885          dist += p->dist;
2886
2887          if (direction == LEFT)
2888                 {
2889                 lsum += dist;
2890                 nl++;
2891                 }
2892          else
2893                 {
2894                 rsum += dist;
2895                 nr++;
2896                 }
2897
2898         if (dist > (*maxdist)) *maxdist = dist;
2899         }
2900
2901       lmean = lsum / nl;
2902       rmean = rsum / nr;
2903
2904       diff = lmean - rmean;
2905       return(diff);
2906       }
2907
2908 float calc_mean(NT_node nptr, float *maxdist, int nseq,NT_node **lu)
2909         {
2910         float dist , lsum = 0.0, rsum = 0.0, lmean,rmean,diff;
2911         NT_node p, *path2root;
2912         float *dist2node;
2913         int depth = 0, i,j , n;
2914         int nl , nr;
2915         int direction, found;
2916
2917
2918         path2root = (NT_node *)vcalloc(nseq,sizeof(Treenode));
2919         dist2node = (float *)vcalloc(nseq,sizeof(float));
2920
2921         depth = (*maxdist) = dist = 0;
2922         nl = nr = 0;
2923         p = nptr;
2924         while (p != NULL)
2925                 {
2926                 path2root[depth] = p;
2927                 dist += p->dist;
2928                 dist2node[depth] = dist;
2929                 p = p->parent;
2930                 depth++;
2931                 }
2932
2933 /***************************************************************************
2934    *nl = *nr = 0;
2935    for each leaf, determine whether the leaf is left or right of the node.
2936    (RIGHT = descendant, LEFT = not descendant)
2937 ****************************************************************************/
2938         for (i=0; i< nseq; i++)
2939                 {
2940                 p = lu[2][i];
2941                 if (p == nptr)
2942                         {
2943                         direction = RIGHT;
2944                         dist = 0.0;
2945                         }
2946                 else
2947                         {
2948                         direction = LEFT;
2949                         dist = 0.0;
2950
2951                         found = FALSE;
2952                         n = 0;
2953                         while ((found == FALSE) && (p->parent != NULL))
2954                                 {
2955                                 for (j=0; j< depth; j++)
2956                                         if (p->parent == path2root[j])
2957                                                 {
2958                                                 found = TRUE;
2959                                                 n = j;
2960                                                 }
2961                                 dist += p->dist;
2962                                 p = p->parent;
2963                                 }
2964
2965                         if (p == nptr) direction = RIGHT;
2966
2967                         }
2968                 if (direction == LEFT)
2969                         {
2970                         lsum += dist;
2971                         lsum += dist2node[n-1];
2972                         nl++;
2973                         }
2974                 else
2975                         {
2976                         rsum += dist;
2977                         nr++;
2978                         }
2979
2980                 if (dist > (*maxdist)) *maxdist = dist;
2981                 }
2982
2983    vfree(dist2node);
2984    vfree(path2root);
2985
2986
2987
2988    if ( nl==0 || nr==0)
2989      {
2990        myexit (EXIT_FAILURE);
2991      }
2992    lmean = lsum / nl;
2993    rmean = rsum / nr;
2994
2995    diff = lmean - rmean;
2996    return(diff);
2997 }
2998
2999 NT_node insert_root(NT_node p, float diff)
3000 {
3001    NT_node newp, prev, q, t;
3002    float dist, prevdist,td;
3003
3004
3005    newp = declare_tree_node( p->maxnseq);
3006    t = p->parent;
3007
3008
3009    prevdist = t->dist;
3010    p->parent = newp;
3011
3012    dist = p->dist;
3013
3014    p->dist = diff / 2;
3015    if (p->dist < 0.0) p->dist = 0.0;
3016    if (p->dist > dist) p->dist = dist;
3017
3018    t->dist = dist - p->dist;
3019
3020    newp->left = t;
3021    newp->right = p;
3022    newp->parent = NULL;
3023    newp->dist = 0.0;
3024    newp->leaf = NODE;
3025
3026    if (t->left == p) t->left = t->parent;
3027    else t->right = t->parent;
3028
3029    prev = t;
3030    q = t->parent;
3031
3032    t->parent = newp;
3033
3034    while (q != NULL)
3035      {
3036         if (q->left == prev)
3037            {
3038               q->left = q->parent;
3039               q->parent = prev;
3040               td = q->dist;
3041               q->dist = prevdist;
3042               prevdist = td;
3043               prev = q;
3044               q = q->left;
3045            }
3046         else
3047            {
3048               q->right = q->parent;
3049               q->parent = prev;
3050               td = q->dist;
3051               q->dist = prevdist;
3052               prevdist = td;
3053               prev = q;
3054               q = q->right;
3055            }
3056     }
3057
3058 /*
3059    remove the old root node
3060 */
3061    q = prev;
3062    if (q->left == NULL)
3063       {
3064          dist = q->dist;
3065          q = q->right;
3066          q->dist += dist;
3067          q->parent = prev->parent;
3068          if (prev->parent->left == prev)
3069             prev->parent->left = q;
3070          else
3071             prev->parent->right = q;
3072          prev->right = NULL;
3073       }
3074    else
3075       {
3076          dist = q->dist;
3077          q = q->left;
3078          q->dist += dist;
3079          q->parent = prev->parent;
3080          if (prev->parent->left == prev)
3081             prev->parent->left = q;
3082          else
3083             prev->parent->right = q;
3084          prev->left = NULL;
3085       }
3086
3087    return(newp);
3088 }
3089
3090
3091
3092
3093 /*********************************************************************/
3094 /*                                                                   */
3095 /*                                  TrimTC3                          */
3096 /*                                                                   */
3097 /*                                                                   */
3098 /*********************************************************************/
3099
3100 int *aln2seq_chain (Alignment *A, int **sim,int seq1, int seq2, int limit, int max_chain);
3101
3102 Alignment *seq2seq_chain (Alignment *A,Alignment*T, char *arg)
3103 {
3104   int **sim=NULL;
3105   int *buf=NULL, *seq2keep, *list, *tname;
3106   int a, b, c, nl;
3107   int sim_limit;
3108   int min_sim=15;
3109   int max_chain=20;
3110
3111   /*Estimate Similarity within the incoming sequences*/
3112   sim=seq2comp_mat (aln2seq(A), "blosum62mt", "sim2");
3113
3114   /*Read and store the list of sequences to keep*/
3115   seq2keep=vcalloc (A->nseq, sizeof (int));
3116   tname=vcalloc (T->nseq, sizeof (int));
3117   for ( a=0; a< T->nseq; a++)
3118     {
3119       tname[a]=name_is_in_list ( T->name[a], A->name, A->nseq, 100);
3120       if (tname[a]>=0)seq2keep[tname[a]]=1;
3121     }
3122
3123   /*Consider Every Pair of Sequences within the list of sequences to keep*/
3124
3125   fprintf ( stderr, "\n");
3126   for ( a=0; a< T->nseq-1; a++)
3127     {
3128       if (tname[a]<0) continue;
3129       for ( b=a+1;b<T->nseq; b++)
3130         {
3131
3132           if (tname[b]<0) continue;
3133
3134           buf=NULL;sim_limit=90;
3135           while (!buf && sim_limit>min_sim)
3136             {
3137               buf=aln2seq_chain ( A, sim,tname[a],tname[b],sim_limit, max_chain);
3138               sim_limit-=5;
3139             }
3140
3141           if ( buf)
3142             {
3143               for (c=0; c< A->nseq; c++)seq2keep[c]+=buf[c];
3144               vfree (buf);
3145             }
3146           else
3147             {
3148               fprintf ( stderr, "\n#Could Not Find any Intermediate sequence [MAx chain %d MinID %d\n", max_chain, min_sim);
3149             }
3150         }
3151     }
3152
3153   list=vcalloc (A->nseq, sizeof (int));
3154   for ( nl=0,a=0; a< A->nseq; a++)
3155     if ( seq2keep[a])
3156       list[nl++]=a;
3157
3158   A=extract_sub_aln (A, nl, list);
3159
3160   free_int (sim, -1);
3161   vfree (list);
3162   return A;
3163 }
3164 int max_explore=10000000;/*Limits the number of explorations that tends to increase when id is small*/
3165 int n_explore;
3166
3167 int *aln2seq_chain (Alignment *A, int **sim, int seq1, int seq2, int limit, int max_chain)
3168 {
3169   int *used;
3170   int **chain;
3171   char output1[10000];
3172   char output2[10000];
3173   int a;
3174   int *list;
3175   int n, nseq=0;
3176
3177
3178   output1[0]=output2[0]='\0';
3179   used=vcalloc (A->nseq, sizeof(int));
3180   used[seq1]=1;
3181
3182   if (find_seq_chain ( A, sim,used,seq1,seq1, seq2,1,limit, max_chain, &nseq))
3183     {
3184       list=vcalloc (A->nseq, sizeof (int));
3185       chain=declare_int (A->nseq, 2);
3186       for (n=0, a=0; a< A->nseq; a++)
3187         {
3188           if ( used[a])
3189             {
3190               chain[n][0]=used[a];
3191               chain[n][1]=a;
3192               list[used[a]-1]=a;n++;
3193             }
3194         }
3195
3196       sprintf ( output2, "#%s %s N: %d Lower: %d Sim: %d DELTA: %d\n", A->name[list[0]], A->name[list[n-1]],n, limit,sim[list[0]][list[n-1]],limit-sim[list[0]][list[n-1]]);strcat (output1, output2);
3197
3198       sort_int ( chain, 2, 0, 0, n-1);
3199       sprintf ( output2, "#");strcat(output1, output2);
3200
3201       for ( a=0; a< n-1; a++)
3202         {
3203           sprintf (output2, "%s -->%d -->", A->name[chain[a][1]],sim[chain[a][1]][chain[a+1][1]]);strcat ( output1, output2);
3204         }
3205       sprintf ( output2, "%s\n", A->name[chain[n-1][1]]);strcat (output1, output2);
3206
3207       free_int (chain, -1);
3208       vfree (list);
3209     }
3210   else
3211     {
3212       vfree (used);
3213       used=NULL;
3214     }
3215   /*  fprintf ( stdout, "%s", output1);*/
3216   fprintf ( stderr, "%s", output1);
3217   n_explore=0;
3218   return used;
3219 }
3220 static int ***pw_sim;
3221 int find_seq_chain (Alignment *A, int **sim,int *used,int seq0,int seq1, int seq2,int chain_length, int limit, int max_chain, int *nseq)
3222 {
3223   int a,b, seq, seq_sim;
3224
3225   n_explore++;
3226   if ( n_explore>=max_explore)
3227     {
3228       return 0;
3229     }
3230   if (!pw_sim)
3231     {
3232       pw_sim=declare_arrayN(3, sizeof (int), A->nseq, A->nseq, 3);
3233       for ( a=0; a< A->nseq; a++)
3234         {
3235           for ( b=0; b<A->nseq; b++)
3236             {
3237               pw_sim[a][b][0]=b;
3238               pw_sim[a][b][1]=sim[a][b];
3239               pw_sim[a][b][2]=sim[b][seq2];
3240             }
3241           sort_int_inv ( pw_sim[a],3, 1, 0, A->nseq-1);
3242         }
3243     }
3244
3245   if ( chain_length>max_chain)return 0;
3246   else if ( sim[seq1][seq2]>=limit)
3247     {
3248       used[seq2]=chain_length+1;
3249       nseq[0]++;
3250       return 1;
3251     }
3252   else
3253     {
3254       int delta_seq2;
3255       for ( a=0; a< A->nseq; a++)
3256         {
3257           seq=pw_sim[seq1][a][0];
3258           seq_sim=pw_sim[seq1][a][1];
3259           delta_seq2=pw_sim[seq1][a][2]-sim[seq1][seq2];
3260
3261
3262
3263           if ( used[seq])continue;
3264           else if ( seq_sim<limit)continue;
3265           else
3266             {
3267               used[seq]=chain_length+1;
3268               nseq[0]++;
3269               if ( find_seq_chain( A, sim,used,seq0,seq, seq2, chain_length+1,limit, max_chain, nseq))
3270                 {return 1;}
3271               else
3272                 used[seq]=0;
3273             }
3274         }
3275       return 0;
3276     }
3277   return 0;
3278 }
3279
3280
3281 //
3282 //
3283 //
3284 //
3285 /*********************************************************************/
3286 /*                                                                   */
3287 /*                                   New Tree Parsing                */
3288 /*                                                                   */
3289 /*********************************************************************/
3290 NT_node new_insert_node (NT_node T);
3291 int scan_name_and_dist ( FILE *fp, char *name, float *dist);
3292
3293
3294
3295
3296 NT_node check_tree (NT_node T);
3297 NT_node main_read_tree (char *treefile)
3298 {
3299   FILE *fp;
3300   Sequence *S;
3301
3302   NT_node T;
3303
3304
3305
3306
3307   fp=vfopen (remove_charset_from_file (treefile, " \t\n\r"), "r");
3308   T=new_get_node (NULL,fp);
3309   vfclose (fp);
3310
3311   S=tree2seq(T, NULL);
3312
3313   T=recode_tree(T, S);
3314   free_sequence (S,S->nseq);
3315   vfree (T->file);
3316   T->file=vcalloc ( strlen (treefile)+1, sizeof (char));
3317   sprintf ( T->file, "%s", treefile);
3318   return T;
3319 }
3320
3321 //This function codes the tree into lseq and lseq2
3322 //lseq:  list of the N->nseq child sequences of the node
3323 //lsseq2:Array of size Nseq, with lseq[a]=1 if sequence a is child of node N
3324 static int node_index;
3325 NT_node index_tree_node    (NT_node T)
3326 {
3327   if (!T)return T;
3328   if (!T->parent){node_index=tree2nseq (T)+1;}
3329
3330   index_tree_node(T->left);
3331   index_tree_node(T->right);
3332
3333   if (!T->left && !T->right)T->index=T->lseq[0]+1;
3334   else T->index=node_index++;
3335   return T;
3336 }
3337
3338
3339
3340 NT_node simple_recode_tree (NT_node T, int nseq)
3341 {
3342
3343   //recodes atree wher the leafs are already coded
3344   if (!T) return T;
3345
3346
3347
3348
3349   T->nseq=0;
3350
3351   if ( T->isseq)
3352     {
3353
3354       ;
3355
3356     }
3357   else
3358     {
3359       NT_node R,L;
3360       int a;
3361        vfree (T->lseq); T->lseq=vcalloc (nseq, sizeof (int));
3362        vfree (T->lseq2); T->lseq2=vcalloc (nseq, sizeof (int));
3363        vfree (T->idist); T->idist=vcalloc (nseq, sizeof (int));
3364        vfree (T->ldist); T->ldist=vcalloc (nseq, sizeof (int));
3365
3366        R=simple_recode_tree (T->left,nseq);
3367
3368        L=simple_recode_tree (T->right,nseq);
3369
3370        if (R)for (a=0; a<R->nseq; a++)
3371          {
3372            T->lseq2[R->lseq[a]]=1;
3373          }
3374
3375        if (L)for (a=0; a<L->nseq; a++)
3376          {
3377            T->lseq2[L->lseq[a]]=1;
3378          }
3379
3380        for (a=0; a<nseq; a++)
3381          {
3382            if (T->lseq2[a])T->lseq[T->nseq++]=a;
3383            if (T->lseq2[a])T->idist[a]=(!R)?0:R->idist[a]+((!L)?0:L->idist[a])+1;
3384            if (T->lseq2[a])T->ldist[a]=(!R)?0:R->ldist[a]+((!L)?0:L->ldist[a])+(int)(T->dist*10000);
3385          }
3386     }
3387   return T;
3388 }
3389
3390 NT_node recode_tree (NT_node T, Sequence *S)
3391 {
3392
3393
3394   if (!T) return T;
3395
3396
3397   vfree (T->lseq);  T->lseq=vcalloc (S->nseq, sizeof (int));
3398   vfree (T->lseq2); T->lseq2=vcalloc (S->nseq, sizeof (int));
3399   vfree (T->idist); T->idist=vcalloc (S->nseq, sizeof (int));
3400   vfree (T->ldist); T->ldist=vcalloc (S->nseq, sizeof (int));
3401   T->nseq=0;
3402
3403   if ( T->isseq)
3404     {
3405
3406       int i;
3407       i=name_is_in_list (T->name, S->name, S->nseq, -1);
3408
3409       if (i!=-1)
3410         {
3411           T->lseq[T->nseq++]=i;
3412           T->lseq2[i]=1;
3413           T->idist[i]=1;
3414           T->ldist[i]=(int)(T->dist*10000);;
3415         }
3416       else
3417         {
3418           printf_exit ( EXIT_FAILURE, stderr, "\nERROR: Sequence %s is in the Tree but Not in the  Sequence dataset [code_lseq][FATAL:%s]", T->name, PROGRAM);
3419         }
3420
3421     }
3422   else
3423     {
3424       NT_node R,L;
3425       int a;
3426
3427       R=recode_tree (T->left, S);
3428
3429       L=recode_tree (T->right, S);
3430
3431       if (R)
3432         for (a=0; a<R->nseq; a++)
3433           {
3434             T->lseq2[R->lseq[a]]=1;
3435           }
3436
3437       if (L)for (a=0; a<L->nseq; a++)
3438         {
3439           T->lseq2[L->lseq[a]]=1;
3440         }
3441
3442       for (a=0; a<S->nseq; a++)
3443         {
3444           //don't count the root
3445           int d;
3446
3447           if ( !(T->parent) || !(T->parent)->parent)d=0;
3448           else if ( T->dist==0)d=0;
3449           else d=1;
3450
3451           if (T->lseq2[a])T->lseq[T->nseq++]=a;
3452           if (T->lseq2[a])T->idist[a]=(!R)?0:(R->idist[a]+((!L)?0:L->idist[a])+d);
3453           if (T->lseq2[a])T->ldist[a]=(!R)?0:R->ldist[a]+((!L)?0:L->ldist[a])+(int)(T->dist*10000);
3454
3455         }
3456     }
3457   return T;
3458 }
3459 int tree2split_list (NT_node T, int ns,int **sl, int* n)
3460 {
3461   if (!T) return 0;
3462   if (!sl) return 0;
3463
3464   tree2split_list (T->right, ns, sl, n);
3465   tree2split_list (T->left , ns, sl, n);
3466
3467   if (!T->right) return 1;
3468   else if (T->parent && !(T->parent)->parent)return 1;
3469   else if ( T->dist==0)return 1;
3470   else
3471     {
3472       int t=0,t2=0, c=0, a;
3473
3474       for (a=0; a< ns; a++)
3475         {
3476           t2+=(a+1)*T->lseq2[a];
3477           t+=T->lseq2[a];
3478         }
3479
3480       if (t2==0) HERE ("0");
3481       c=(t>(ns-t))?1:0;
3482       sl[n[0]][ns]=t2;//Hash value for quick comparison;
3483
3484       for (a=0; a< ns; a++)sl[n[0]][a]=(c==0)?T->lseq2[a]:(1-T->lseq2[a]);
3485       n[0]++;
3486
3487     }
3488   return 1;
3489 }
3490
3491 NT_node display_splits (NT_node T,Sequence *S, FILE *fp)
3492 {
3493   int a;
3494   if (!T) return T;
3495
3496   if (!S)S=tree2seq (T,NULL);
3497
3498   display_splits (T->right,S, fp);
3499   display_splits (T->left, S, fp);
3500
3501
3502
3503   if (!T->right);
3504   else if (T->parent && !(T->parent)->parent);
3505   else
3506     {
3507       int t=0;
3508       for (a=0; a< S->nseq; a++)
3509         {
3510           fprintf (fp, "%d", T->lseq2[a]);
3511           t+=T->lseq2[a];
3512         }
3513
3514       fprintf ( fp, " %5d \n", MIN(t,((S->nseq)-t)));
3515     }
3516   return T;
3517 }
3518 NT_node display_leaf_nb (NT_node T, int n, FILE *fp, char * name)
3519 {
3520   int a;
3521   if (!T) return T;
3522
3523
3524   display_leaf_nb (T->right, n, fp, name);
3525   display_leaf_nb (T->left, n, fp, name);
3526
3527
3528   if (!T->isseq);
3529   else
3530     {
3531       NT_node P;
3532
3533       P=T->parent;
3534       fprintf (fp, "%s ", T->name);
3535       for (a=0; a< n; a++)fprintf (fp, "%d", P->lseq2[a]);
3536       fprintf ( fp," %s\n", name);
3537     }
3538   return T;
3539 }
3540 static int root4dc;
3541 NT_node display_code (NT_node T, int n, FILE *fp)
3542 {
3543   int a, debug=0, t=0;
3544   if (!T) return T;
3545
3546   if (!T->parent)
3547     root4dc=0;
3548
3549
3550
3551   if (!T->parent && debug) fprintf ( fp, "\nDISPLAY TREE: START");
3552   display_code (T->right, n, fp);
3553   display_code (T->left, n, fp);
3554
3555   fprintf ( fp, "\n");
3556   if (!T->parent) return T;
3557   else if ( !(T->parent)->parent && root4dc==1)return T;
3558   else if ( !(T->parent)->parent && root4dc==0)root4dc=1;
3559
3560   for (a=0; a< n; a++)
3561     t+=T->lseq2[a];
3562   if ( t<=n/2)
3563     for (a=0; a< n; a++)fprintf (fp, "%d", T->lseq2[a]);
3564   else
3565     for (a=0; a< n; a++)fprintf (fp, "%d", 1-T->lseq2[a]);
3566   if (T->isseq && debug)fprintf (fp, "%s", T->name);
3567
3568   if (!T->parent && debug) fprintf (fp, "\nDISPLAY TREE: FINISHED");
3569   return T;
3570 }
3571 NT_node display_dist (NT_node T, int n, FILE *fp)
3572 {
3573   int a;
3574   if (!T) return T;
3575
3576   if (!T->parent)
3577     root4dc=0;
3578
3579   display_dist (T->right, n, fp);
3580   display_dist (T->left, n, fp);
3581
3582   fprintf ( stdout, "\n");
3583   for ( a=0; a< n; a++)
3584     fprintf ( stdout, " %2d ", T->idist[a]);
3585   fprintf ( stdout, "\n");
3586
3587   return T;
3588 }
3589
3590 NT_node check_tree (NT_node T)
3591 {
3592   if (T) HERE("CHECK %s", T->name);
3593   if (!T)
3594     {
3595       HERE ("ERROR: Empty Group");
3596     }
3597
3598   else if (T->isseq)return T;
3599   else
3600     {
3601       HERE ("R");
3602       check_tree (T->right);
3603       HERE ("L");
3604       check_tree (T->left);
3605       return NULL;
3606     }
3607   return 0;}
3608
3609 NT_node new_reroot_tree( NT_node T)
3610 {
3611   T=unroot_tree (T);
3612   return T;
3613 }
3614 NT_node new_get_node (NT_node T, FILE *fp)
3615 {
3616   NT_node NN;
3617   int c;
3618   static int n;
3619
3620   c=fgetc (fp);
3621   if (!T)T=declare_tree_node (100);
3622
3623
3624   if ( c==';')
3625     {
3626
3627       if (!T->right)T=T->left;
3628       else if (!T->left)T=T->right;
3629       vfree (T->parent);T->parent=NULL;
3630       return T;
3631     }
3632   else if ( c==')')
3633     {
3634       --n;
3635       scan_name_and_dist (fp, T->name, &T->dist);
3636       if (T->name && T->name [0])T->bootstrap=atof (T->name);
3637       return new_get_node (T->parent, fp);
3638     }
3639   else if ( c==',')
3640     {
3641       return new_get_node (T, fp);
3642     }
3643   else
3644     {
3645       NN=new_insert_node (T);
3646
3647       if ( c=='(')
3648         {
3649           ++n;
3650           return new_get_node (NN, fp);
3651         }
3652       else
3653         {
3654           ungetc (c, fp);
3655           scan_name_and_dist (fp, NN->name, &NN->dist);
3656
3657           NN->leaf=1;
3658           NN->isseq=1;
3659           return new_get_node (T, fp);
3660         }
3661     }
3662 }
3663 int scan_name_and_dist ( FILE *fp, char *name, float *dist)
3664 {
3665   int a, c;
3666   char number [1000];
3667
3668   a=0;
3669   c=fgetc (fp);ungetc (c, fp);
3670
3671
3672   if ( c==';')return 0;
3673
3674   while ((c=fgetc(fp))!=':' && c!=EOF && c!=')' && c!=';' && c!=',')
3675     {
3676       name[a++]=c;
3677     }
3678   name [a]='\0';
3679
3680   if ( c!=':')
3681     {
3682       ungetc (c, fp);
3683       dist[0]=FLT_MIN;
3684       return 1;
3685     }
3686   a=0;
3687   while (isdigit((c=fgetc(fp))) || c=='.' || c=='-')
3688     {
3689       number[a++]=c;
3690     }
3691
3692   ungetc (c, fp);
3693   number[a]='\0';
3694
3695   dist[0]=atof (number);
3696
3697   return 2;
3698 }
3699 NT_node new_insert_node (NT_node T)
3700 {
3701   NT_node NN;
3702
3703
3704   NN=new_declare_tree_node ();
3705   NN->parent=T;
3706   if (!T)
3707     {
3708       return NN;
3709     }
3710   else if (T->left==NULL)
3711     {
3712       T->left=NN;
3713
3714     }
3715   else if ( T->right==NULL)
3716     {
3717       T->right=NN;
3718     }
3719   else
3720     {
3721       NT_node NN2;
3722       NN2=new_declare_tree_node ();
3723       NN2->left=T->left;
3724       NN2->right=T->right;
3725       NN2->parent=T;
3726
3727       T->left=NN2;
3728       T->right=NN;
3729     }
3730
3731   /*
3732
3733   else
3734     {
3735       NN->right=T->right;
3736       (T->right)->parent=NN;
3737
3738       NN->parent=T;
3739       T->right=NN;
3740       NN->left=new_declare_tree_node ();
3741       (NN->left)->parent=NN;
3742       return NN->left;
3743     }
3744   */
3745   /*
3746     This caused a crash when internal undefined nodes, removed 19/02/08
3747    else
3748     {
3749       NT_node P;
3750       NN->right=T;
3751       P=NN->parent=T->parent;
3752       T->parent=NN;
3753
3754       if (P && P->right==T)P->right=NN;
3755       else if ( P && P->left==T)P->left=NN;
3756
3757       NN->left=new_declare_tree_node ();
3758       (NN->left)->parent=NN;
3759       return NN->left;
3760     }
3761   */
3762   return NN;
3763 }
3764
3765 NT_node new_declare_tree_node ()
3766 {
3767         NT_node p;
3768         static int node_index;
3769         p= (NT_node)vcalloc (1, sizeof ( Treenode));
3770         p->left = NULL;
3771         p->right = NULL;
3772         p->parent = NULL;
3773         p->dist = 0.0;
3774         p->leaf = 0;
3775         p->order = 0;
3776         p->index=++node_index;
3777         p->maxnseq=1000;
3778         p->name=(char*)vcalloc (MAXNAMES+1,sizeof (char));
3779         p->name[0]='\0';
3780         return p;
3781
3782         }
3783 int new_display_tree (NT_node T, int n)
3784 {
3785   int in;
3786
3787   in=n;
3788
3789
3790   if ( T->parent)fprintf (stdout, "\nNode %d: has parents)", in);
3791   else fprintf (stdout, "\nNode %d: NO parents)", in);
3792
3793   if ( T->right)
3794     {
3795       fprintf (stdout, "\nNode %d has Right Child", in);
3796       n=new_display_tree (T->right, n+1);
3797     }
3798   else fprintf ( stdout, "\nNode %d No Right\n", in);
3799
3800   if ( T->left)
3801     {
3802       fprintf (stdout, "\nNode %d has Left Child", in);
3803       n=new_display_tree (T->left, n+1);
3804     }
3805   else fprintf ( stdout, "\nNode %d No Left\n", in);
3806
3807   if ( T->bot)
3808     {
3809       fprintf (stdout, "\nNode %d has Bot Child", in);
3810       n=new_display_tree (T->bot, n+1);
3811     }
3812   else fprintf ( stdout, "\nNode %d No Bot\n", in);
3813
3814
3815   if (T->isseq)
3816     {
3817       fprintf (stdout, "\nNode %d is %s", in, T->name);
3818       return in;
3819     }
3820   else return 0;}
3821 int display_tree_duplicates (NT_node T)
3822 {
3823   static Sequence *S;
3824   static int *dup;
3825   int a, b;
3826
3827   free_sequence (S, -1);
3828   vfree (dup);
3829
3830   S=tree2seq (T, NULL);
3831   dup=vcalloc ( S->nseq, sizeof (int));
3832
3833   for (a=0; a< S->nseq-1; a++)
3834     for ( b=a+1; b<S->nseq; b++)
3835       {
3836         if ( strm (S->name[a], S->name[b]))
3837           {
3838             dup[a]++;
3839           }
3840       }
3841   for (a=0; a< S->nseq-1; a++)
3842     for ( b=a+1; b<S->nseq; b++)
3843       {
3844         if ( strm (S->name[a], S->name[b]) && dup[a])
3845           {
3846             fprintf ( stderr, "\nSequence %s is duplicated %d Times in the tree", S->name[a], dup[a]);
3847             dup[a]=0;
3848           }
3849       }
3850   return 0;
3851 }
3852 int tree_contains_duplicates (NT_node T)
3853 {
3854   static Sequence *S;
3855   int a, b;
3856
3857   free_sequence (S, -1);
3858
3859   S=tree2seq (T, NULL);
3860   for (a=0; a< S->nseq-1; a++)
3861     for ( b=a+1; b<S->nseq; b++)
3862       {
3863         if ( strm (S->name[a], S->name[b]))return 1;
3864       }
3865   return 0;
3866 }
3867
3868 float display_avg_bootstrap ( NT_node T)
3869 {
3870   float tot;
3871   int n;
3872
3873   tot=tree2tot_dist (T, BOOTSTRAP);
3874   n=tree2n_branches (T, BOOTSTRAP);
3875   fprintf ( stdout, "\nAVERAGE BOOTSRAP: %.3f on %d Branches\n", (n>0)?tot/n:0, n);
3876   return (n>0)?tot/n:0;
3877 }
3878
3879
3880 int tree2n_branches(NT_node T, int mode)
3881 {
3882   int n=0;
3883
3884   if (!T) return 0;
3885   if (!T->parent);
3886   else if  ((T->isseq && mode !=BOOTSTRAP) || !T->isseq)
3887     {
3888       n++;
3889     }
3890   n+=tree2n_branches(T->right, mode);
3891   n+=tree2n_branches(T->left, mode);
3892
3893   return n;
3894 }
3895
3896 float tree2tot_dist ( NT_node T, int mode)
3897 {
3898   float t=0;
3899
3900
3901   if ( !T)return 0;
3902
3903   if ( !T->parent);
3904   else if  ((T->isseq && mode !=BOOTSTRAP) || !T->isseq)
3905     {
3906       if ( mode == BOOTSTRAP && T->bootstrap!=0)t+=T->bootstrap;
3907       else t+=T->dist;
3908     }
3909
3910   t+=tree2tot_dist(T->right, mode);
3911   t+=tree2tot_dist(T->left, mode);
3912   return t;
3913 }
3914
3915 //This function displays all the sequences within the tree sorted by node label
3916 int cmp_tree_array ( const void *vp, const void *vq);
3917 int node_sort ( char *name, NT_node T)
3918 {
3919   NT_node N;
3920   int nseq;
3921   int **array, a;
3922   Sequence *S;
3923   while (T->parent)T=T->parent;
3924
3925   nseq=tree2nseq (T);
3926   array=declare_int (nseq, 2);
3927   N=tree2node (name, T);
3928
3929   if (N==NULL)printf_exit (EXIT_FAILURE, stderr, "ERROR: %s is not in the tree [FATAL:%s]\n", name, PROGRAM);
3930   array=display_tree_from_node (N,0,0, array);
3931   qsort ( array, nseq, sizeof (int**), cmp_tree_array);
3932   S=tree2seq(T, NULL);
3933   for (a=0; a<nseq; a++)
3934     fprintf ( stdout, ">%s %d %d\n", S->name[array[a][0]], array[a][1], array[a][2]);
3935   myexit (EXIT_SUCCESS);
3936 }
3937
3938 NT_node tree2root ( NT_node R)
3939 {
3940   if (R)while (R->parent)R=R->parent;
3941   return R;
3942 }
3943
3944 NT_node tree2node (char *name, NT_node T)
3945 {
3946   NT_node T1, T2;
3947   if ( !T) return T;
3948   else if (T->leaf && strm (T->name, name)) return T;
3949   else
3950     {
3951
3952       T1=tree2node ( name, T->right);
3953       T2=tree2node ( name, T->left);
3954       return (T1>T2)?T1:T2;
3955   }
3956
3957 }
3958 static int ni;
3959 NT_node * tree2node_list (NT_node T, NT_node *L)
3960 {
3961
3962
3963   if (!T) return NULL;
3964   if (!L) {ni=0;L=vcalloc (tree2nnode(T)+1, sizeof (NT_node));}
3965   tree2node_list (T->left, L);
3966   tree2node_list (T->right, L);
3967   L[ni++]=T;
3968   return L;
3969 }
3970
3971
3972
3973 int ** display_tree_from_node (NT_node T, int up, int down, int **array)
3974 {
3975
3976   if (!T || T->visited)return array;
3977
3978   T->visited=1;
3979   if (T->isseq)
3980     {
3981       array[T->lseq[0]][0]=T->lseq[0];
3982       array[T->lseq[0]][1]=up;
3983       array[T->lseq[0]][2]=down;
3984
3985     }
3986   else
3987     {
3988       array=display_tree_from_node ( T->left ,up, down+1, array);
3989       array=display_tree_from_node ( T->right,up, down+1, array);
3990     }
3991   array=display_tree_from_node ( T->parent,up+1, 0, array);
3992   T->visited=0;
3993   return array;
3994
3995 }
3996
3997 int cmp_tree_array ( const void *p, const void *q)
3998 {
3999   const int **vp=(const int**)p;
4000   const int **vq=(const int**)q;
4001   if (vp[0][1]>vq[0][1])return 1;
4002   else if ( vp[0][1]<vq[0][1]) return -1;
4003   else if ( vp[0][2]>vq[0][2]) return 1;
4004   else if ( vp[0][2]<vq[0][2]) return -1;
4005   else return 0;
4006 }
4007
4008 NT_node * read_tree_list (Sequence *S)
4009 {
4010   NT_node *T;
4011   int a;
4012
4013   T=vcalloc ( S->nseq+1, sizeof (NT_node));
4014
4015   for ( a=0; a<S->nseq; a++)
4016     {
4017       char *fname;
4018       if (S->seq && S->seq[a] && strlen (S->seq[a])<2)
4019         fname=S->name[a];
4020       else
4021         string2file ((fname=vtmpnam(NULL)), "w", S->seq[a]);
4022
4023       T[a]=main_read_tree (fname);
4024       T[a]->file=vcalloc (strlen (S->name[a])+1, sizeof (char));
4025       sprintf (T[a]->file, "%s", S->name[a]);
4026     }
4027   return T;
4028 }
4029
4030 int treelist2dmat ( Sequence *S)
4031 {
4032   NT_node *T;
4033   int n=0, a, b;
4034   float v;
4035   Sequence *TS;
4036
4037
4038
4039   n=S->nseq;
4040   T=read_tree_list (S);
4041   TS=tree2seq(T[0], NULL);
4042   fprintf (stdout, "\n%d", S->nseq);
4043   for (a=0; a<n; a++)
4044     {
4045       fprintf ( stdout,"\n%-10s ", S->name[a]);
4046       for ( b=0; b<n; b++)
4047         {
4048           v=100-simple_tree_cmp (T[a], T[b], TS, 1);
4049           fprintf ( stdout, "%.4f ", v);
4050         }
4051
4052     }
4053
4054   myexit (EXIT_SUCCESS);
4055   return 0;
4056 }
4057
4058 int treelist2leafgroup ( Sequence *S, Sequence *TS, char *taxon)
4059 {
4060   NT_node *T;
4061   int n=0,nseq, a, c,s;
4062
4063   int *used;
4064
4065   char *split_file, *sorted_split_file;
4066   char *buf=NULL;
4067   FILE *fp;
4068   char *name, *fname, *group, *ref_group, *list;
4069
4070
4071
4072   T=read_tree_list (S);
4073   if (!TS)TS=tree2seq(T[0], NULL);
4074
4075   name=vcalloc (1000, sizeof (char));
4076   fname=vcalloc (1000, sizeof (char));
4077   group=vcalloc (TS->nseq*10, sizeof (char));
4078   ref_group=vcalloc (TS->nseq*10, sizeof (char));
4079   list=vcalloc (100*S->nseq, sizeof (char));
4080   split_file=vtmpnam (NULL);
4081   sorted_split_file =vtmpnam (NULL);
4082
4083   n=S->nseq;
4084   used=vcalloc (n, sizeof (int));
4085
4086   T=read_tree_list (S);
4087   if (!TS)TS=tree2seq(T[0], NULL);
4088   nseq=TS->nseq;
4089   fp=vfopen (split_file, "w");
4090
4091   for ( a=0; a< S->nseq; a++)
4092     {
4093
4094       T[a]=prune_tree  (T[a], TS);
4095       T[a]=recode_tree (T[a], TS);
4096       display_leaf_nb (T[a], TS->nseq,fp, S->name[a]);
4097     }
4098   vfclose (fp);
4099
4100
4101   for (s=0; s< TS->nseq; s++)
4102     {
4103       int i;
4104
4105
4106       if (taxon && !(strm (taxon, TS->name[s]) ))continue;
4107       else
4108         printf_system ( "cat %s | grep %s| sort > %s::IGNORE_FAILURE::", split_file,TS->name[s], sorted_split_file);
4109
4110       vfopen (sorted_split_file, "r");
4111       ref_group[0]=group[0]='\0';
4112
4113       while ( (c=fgetc (fp))!=EOF)
4114         {
4115
4116           ungetc (c, fp);
4117           buf=vfgets (buf, fp);
4118           sscanf (buf, "%s %s %s\n", name, group, fname);
4119
4120           if ( !ref_group[0]|| !strm (group, ref_group))
4121             {
4122               if (ref_group[0])
4123
4124                 {fprintf (stdout, "%s %6.2f %s",name, (((float)n*100)/(float)S->nseq), ref_group);
4125                   for (i=0,a=0; a<nseq; a++)
4126                     if (ref_group[a]=='1' && a!=s)fprintf (stdout, " %s ", TS->name[a]);
4127                   fprintf ( stdout, "\n");
4128                   fprintf (stdout, "\nLIST: %s\n", list);
4129                 }
4130               list[0]='\0';
4131               sprintf ( ref_group, "%s", group);
4132               list=strcatf (list, " %s", fname);
4133               n=1;
4134             }
4135           else
4136             {
4137
4138               list=strcatf (list, " %s", fname);
4139               n++;
4140             }
4141         }
4142
4143       fprintf (stdout, "%s %6.2f %s",name, (((float)n*100)/(float)S->nseq), group);
4144       for (i=0,a=0; a<nseq; a++)
4145         if (group[a]=='1' && a!=s)fprintf (stdout, " %s ", TS->name[a]);
4146       fprintf (stdout, "\nLIST %s\n", list);
4147       fprintf ( stdout, "\n");
4148       vfclose (fp);
4149     }
4150
4151   myexit (0);
4152 }
4153 int count_tree_groups( Sequence *LIST, char *group_file)
4154 {
4155   NT_node *T;
4156   Sequence *S;
4157   int a, b, c,n, w, wo, ng=0;
4158   int  **list, ***rlist, **blist;
4159   char ***l;
4160   int *gs;
4161
4162
4163
4164   T=read_tree_list (LIST);
4165   S=tree2seq(T[0], NULL);
4166   for ( a=0; a< LIST->nseq; a++)
4167     {
4168       T[a]=prune_tree  (T[a], S);
4169       T[a]=recode_tree (T[a], S);
4170     }
4171
4172
4173
4174   gs=vcalloc (2, sizeof (int));
4175   list=declare_int (LIST->nseq*S->nseq*2, S->nseq+1);
4176
4177   blist=declare_int (2, S->nseq+1);
4178   for ( n=0, a=0; a< LIST->nseq; a++)
4179     {
4180       int n2=0;
4181       tree2split_list (T[a], S->nseq, list+n, &n2);
4182       n+=n2;
4183
4184
4185       for (b=0; b<n2; b++)
4186         {
4187           for ( c=0; c<S->nseq; c++)
4188             list[n+b][c]=1-list[n-n2+b][c];
4189         }
4190       n+=n2;
4191     }
4192
4193   if ( group_file)
4194     {
4195       rlist=declare_arrayN(3, sizeof (int), 2,LIST->nseq*S->nseq, S->nseq+1);
4196       l=file2list (group_file, " ");
4197
4198       while (l[ng])
4199         {
4200           int i, b, g;
4201           if (!strstr (l[ng][1], "group")){ng++;continue;}
4202           g=(strm (l[ng][1], "group2"))?0:1;
4203
4204           for (b=2; b<atoi (l[ng][0]); b++)
4205             {
4206               if ((i=name_is_in_list(l[ng][b], S->name, S->nseq, 100))!=-1)rlist[g][gs[g]][i]=1;
4207             }
4208           gs[g]++;
4209           ng++;
4210         }
4211     }
4212   else
4213     {
4214       rlist=vcalloc ( 2, sizeof (int**));
4215       rlist[1]=count_int_strings (list, n, S->nseq);
4216       gs[1]=read_array_size_new (rlist[1]);
4217
4218       rlist[0]=declare_int (S->nseq, S->nseq);
4219       gs[0]=S->nseq;
4220       for ( a=0; a<S->nseq; a++)rlist[0][a][a]=1;
4221     }
4222
4223
4224   for (wo=w=0,a=0; a<gs[0]; a++)
4225     {
4226       for (c=0; c< gs[1]; c++)
4227         {
4228           wo=w=0;
4229           for (b=0; b<S->nseq; b++)
4230             {
4231               blist[0][b]=blist[1][b]=rlist[1][c][b];
4232               blist[0][b]=(rlist[0][a][b]==1)?1:blist[0][b]; //WITH GROUP 1
4233               blist[1][b]=(rlist[0][a][b]==1)?0:blist[1][b]; //wiTHOUT gROUP 1
4234
4235             }
4236           for (b=0; b<n; b++)
4237             {
4238               int x1, x2;
4239               x1=(memcmp (blist[0], list[b], sizeof (int)*S->nseq)==0)?1:0;
4240               w+=x1;
4241               x2=(memcmp (blist[1], list[b], sizeof (int)*S->nseq)==0)?1:0;
4242               wo+=x2;
4243
4244             }
4245           fprintf ( stdout, "\n%d ", MIN(wo, w));
4246           fprintf ( stdout, "(");
4247           for (b=0; b<S->nseq; b++)if (rlist[1][c][b])fprintf ( stdout, "%s ",S->name[b]);
4248           fprintf ( stdout, ") +/- (");
4249           for (b=0; b<S->nseq; b++)if (rlist[0][a][b])fprintf ( stdout, "%s ",S->name[b]);
4250
4251           fprintf (stdout , ") + %d - %d Delta %d", w, wo, FABS((wo-w)));
4252         }
4253     }
4254   myexit (0);
4255 }
4256 Split * print_split ( int n, int **list, Sequence *LIST, Sequence *S, char *buf, char *file);
4257
4258 NT_node split2tree ( NT_node RT,Sequence *LIST, char *param)
4259 {
4260   Split **S;
4261   Alignment *A;
4262   S=count_splits (RT, LIST, param);
4263   A=seq2aln ((S[0])->S,NULL, KEEP_GAP);
4264
4265   return split2upgma_tree (S,A, A->nseq, "no");
4266 }
4267
4268 Split** count_splits( NT_node RT,Sequence *LIST, char *param)
4269 {
4270   NT_node *T, OrderT;
4271   Sequence *S=NULL;
4272   int a, b, c, d, n1, n2;
4273   int  **list1, **list2;
4274   Split **SL;
4275   int nb, tlist;
4276   char *main_buf;
4277   char *in=NULL,*in2=NULL, *out=NULL, order[100], filter[100];
4278   FILE *fp, *fp2;
4279   static char *def_param;
4280   char *cache=NULL;
4281   //+count_splits _NB_x_FILTER_<file>
4282   //_<file is a fasta file containing the list of species to keep>
4283   if (!def_param)def_param=vcalloc ( 10, sizeof (char));
4284
4285
4286
4287   if (!param)param=def_param;
4288
4289
4290   strget_param (param, "_NB_", "0", "%d", &nb);
4291   strget_param (param, "_TLIST_", "0", "%d", &tlist);
4292   strget_param (param, "_ORDER_", "NO", "%s", order);
4293   strget_param (param, "_FILTER_", "NO", "%s", filter);
4294
4295   fprintf ( stderr, "\nREAD TREE LIST [%d Trees...", LIST->nseq);
4296   T=read_tree_list (LIST);
4297   fprintf ( stderr, "..]");
4298
4299   if ( !(strm (order, "NO")))
4300     {
4301       if (is_newick (order))
4302         {
4303           OrderT=main_read_tree (order);
4304         }
4305       else
4306         {
4307           S=main_read_seq (order);
4308         }
4309     }
4310   else
4311     {
4312       OrderT=(RT)?RT:T[0];
4313     }
4314   fprintf ( stderr, "\nTrees Ordered according to: %s", (strm (order, "NO"))?"First Tree":order);
4315
4316
4317   if (!S)S=tree2seq(OrderT, NULL);
4318
4319   for (a=0; a<S->nseq; a++)
4320     {
4321       fprintf ( stdout, "\n#ORDER %15s : %3d", S->name[a], a+1);
4322     }
4323   if ( !strm (filter, "NO"))
4324     {
4325       Sequence *F;
4326       int i;
4327
4328       F=main_read_seq (filter);
4329       cache=vcalloc (S->nseq, sizeof (int));
4330       for ( a=0; a<F->nseq; a++)
4331         {
4332           if ( (i=name_is_in_list (F->name[a], S->name, S->nseq, 100))!=-1)
4333             cache[i]=1;
4334         }
4335       free_sequence (F, -1);
4336     }
4337
4338   main_buf=vcalloc ( S->nseq*(STRING+1), sizeof(int));
4339
4340   list1=declare_int (S->nseq*3, S->nseq+1);
4341   list2=declare_int (S->nseq*3, S->nseq+1);
4342
4343   for ( a=0; a< LIST->nseq; a++)
4344     {
4345       T[a]=prune_tree  (T[a], S);
4346       T[a]=recode_tree (T[a], S);
4347     }
4348
4349
4350
4351   if (!RT)
4352     {
4353       char *buf;
4354       int i,nl;
4355
4356       in=vtmpnam (NULL);in2=vtmpnam(NULL); out=vtmpnam (NULL);
4357
4358       fp=vfopen (in, "w");
4359       fp2=vfopen (in2, "w");
4360       for ( a=0; a< LIST->nseq; a++)
4361         {
4362           n2=0;
4363           tree2split_list (T[a], S->nseq, list2, &n2);
4364           for ( b=0; b<n2; b++)
4365             {
4366               for (c=0; c< S->nseq; c++)
4367                 {fprintf (fp, "%d", list2[b][c]);}
4368               fprintf (fp, "\n");
4369               for (c=0; c< S->nseq; c++)
4370                 {fprintf (fp, "%d", 1-list2[b][c]);}
4371               fprintf (fp, "\n");
4372
4373               for (c=0; c< S->nseq; c++)
4374                 {fprintf (fp2, "%d", list2[b][c]);}
4375               fprintf (fp2, " ");
4376               for (c=0; c< S->nseq; c++)
4377                 {fprintf (fp2, "%d", 1-list2[b][c]);}
4378               fprintf (fp2, " %s\n",LIST->name[a]);
4379             }
4380         }
4381       vfclose (fp2);
4382       vfclose (fp);
4383
4384       count_strings_in_file (in, out);
4385       nl=count_n_line_in_file(out);
4386       list1=declare_int (nl+1, S->nseq+2);
4387
4388       fp=vfopen (out, "r");
4389       n1=0;
4390       buf=vcalloc (measure_longest_line_in_file (out)+1, sizeof (char));
4391       while ( fscanf (fp, "%s %d",buf, &i)==2)
4392         {
4393           for (a=0; a<S->nseq; a++)list1[n1][a]=buf[a]-'0';
4394           list1[n1++][S->nseq+1]=i;
4395         }
4396       vfclose (fp);
4397       vfree (buf);
4398     }
4399   else
4400     {
4401
4402       RT=prune_tree  (RT, S);
4403       RT=recode_tree (RT, S);
4404       n1=0;
4405       tree2split_list (RT, S->nseq, list1,&n1);
4406       for ( a=0; a< LIST->nseq; a++)
4407         {
4408           n2=0;
4409           tree2split_list (T[a], S->nseq, list2, &n2);
4410           for (b=0; b<n1; b++)
4411             {
4412               for ( c=0; c<n2; c++)
4413                 {
4414                   int di=0;
4415                   for (d=0; d<S->nseq; d++)
4416                     {
4417                       if (list1[b][d]!=list2[c][d])di++;
4418                     }
4419                   list1[b][S->nseq+1]+=(di==0 || di== S->nseq)?1:0;
4420                 }
4421
4422             }
4423         }
4424     }
4425   SL=vcalloc ( n1+1, sizeof (Split*));
4426
4427   for (a=0; a<n1; a++)
4428     {
4429       int s1, s2;
4430       int cont=1;
4431       if (nb)fprintf ( stdout, "\nSPLIT: %d and its Neighborhood +/^- %d\n", a+1, nb);
4432
4433       if (cache)
4434         for (b=0; b<S->nseq; b++)if (cache[b]!=list1[a][b])cont=0;
4435       if (!cont) continue;
4436
4437       SL[a]=print_split (a, list1, LIST, S, main_buf, (tlist==1)?in2:NULL);
4438       for (b=0; b<n1; b++)
4439         {
4440           if ( a==b)continue;
4441           else
4442             {
4443
4444               for (d=s1=s2=0,c=0; c<S->nseq; c++)
4445                 {
4446                   s1+=list1[b][c];
4447                   s2+=list1[a][c];
4448                   d+=(list1[a][c]!=list1[b][c])?1:0;
4449                 }
4450
4451             }
4452           if (d<=nb &&((s1==s2)|| ((S->nseq-s1)==s2)))print_split (b, list1, LIST, S, main_buf, (tlist==1)?in2:NULL);
4453         }
4454     }
4455   a=0;
4456   vfree (cache);
4457   return SL;
4458 }
4459 Split * declare_split (int nseq, int ntrees);
4460 Split* print_split ( int a, int **list1, Sequence *LIST, Sequence *S, char *buf, char *split_file)
4461   {
4462     int f1,t,b;
4463     Split *SP=NULL;
4464
4465
4466
4467     SP=declare_split (S->nseq, LIST->nseq);
4468
4469     fprintf ( stdout, "\n>");
4470     for (t=0,b=0; b<S->nseq; b++){fprintf ( stdout, "%d", list1[a][b]);t+=list1[a][b];SP->split[b]='0'+list1[a][b];}
4471     fprintf ( stdout, " NumberSplit %5d SplitSize %5d Score %5.2f %s ", list1[a][S->nseq+1],t, (float)(list1[a][S->nseq+1]*100)/LIST->nseq, (buf)?buf:"");
4472     SP->n= list1[a][S->nseq+1];
4473     SP->score=(float)(list1[a][S->nseq+1]*100)/LIST->nseq;
4474     SP->S=S;
4475
4476     for (f1=1,b=0; b< S->nseq; b++)
4477       {
4478
4479         if (list1[a][b])
4480           {
4481             if (f1==1)fprintf ( stdout, "(");
4482             else fprintf (stdout, ",");
4483             f1=0;
4484             fprintf ( stdout, "%s", S->name [b]);
4485           }
4486       }
4487     fprintf ( stdout, ")");
4488     if (split_file)
4489       {
4490         char *buf=NULL;
4491         FILE *fp;
4492
4493         char c;
4494         fp=vfopen (split_file, "r");
4495         while ( (c=fgetc(fp))!=EOF)
4496           {
4497
4498             c=ungetc (c, fp);
4499             buf=vfgets (buf, fp);
4500             if ( strstr (buf, SP->split))
4501               {
4502                 char **list;
4503                 list=string2list (buf);
4504                 fprintf ( stdout, "\n\t%s %s", SP->split, list[3]);
4505                 free_char (list, -1);
4506               }
4507           }
4508         vfclose (fp);
4509       }
4510
4511     return SP;
4512   }
4513 Split * declare_split (int nseq, int ntrees)
4514 {
4515   Split *S;
4516   S=vcalloc (1, sizeof (Split));
4517   S->split=vcalloc ( nseq+1, sizeof (char));
4518   return S;
4519 }
4520 int treelist2splits( Sequence *S, Sequence *TS)
4521 {
4522   NT_node *T;
4523   int n=0,nseq, a, c;
4524
4525   int *used;
4526
4527   char *split_file, *sorted_split_file;
4528   char *buf=NULL, *ref_buf=NULL;
4529   FILE *fp;
4530
4531   split_file=vtmpnam (NULL);
4532   sorted_split_file =vtmpnam (NULL);
4533
4534   n=S->nseq;
4535   used=vcalloc (n, sizeof (int));
4536
4537   T=read_tree_list (S);
4538   if (!TS)TS=tree2seq(T[0], NULL);
4539   nseq=TS->nseq;
4540   fp=vfopen (split_file, "w");
4541
4542
4543   for ( a=0; a< S->nseq; a++)
4544     {
4545
4546       T[a]=prune_tree  (T[a], TS);
4547       T[a]=recode_tree (T[a], TS);
4548       display_splits (T[a], TS,fp);
4549     }
4550
4551   vfclose (fp);
4552   printf_system ("cp %s split_file::IGNORE_FAILURE::", split_file);
4553   printf_system ( "cat %s | grep 1| sort > %s::IGNORE_FAILURE::", split_file, sorted_split_file);
4554
4555   fp=vfopen (sorted_split_file, "r");
4556   fprintf (stdout, "LEGEND: <#occurences> <coded split> <min group size> <(group1,)> <(group2,>\n");
4557
4558   for ( a=0; a<TS->nseq; a++)fprintf ( stdout, "SEQ_INDEX %d %s\n", a+1, TS->name[a]);
4559   while ( (c=fgetc (fp))!=EOF)
4560     {
4561
4562       ungetc (c, fp);
4563       buf=vfgets (buf, fp);
4564       buf [strlen(buf)-1]='\0';
4565
4566       if ( ref_buf==NULL)
4567         {
4568           ref_buf=vcalloc (strlen (buf)+1, sizeof (char));
4569           sprintf ( ref_buf, "%s", buf);
4570           n=1;
4571         }
4572       else if ( !strm (buf, ref_buf))
4573         {
4574           int i;
4575           fprintf ( stdout, "SPLIT_COUNT %3d %s (", n, ref_buf);
4576           for (i=0,a=0; a<nseq; a++)
4577             if (ref_buf[a]=='1')
4578               {
4579                 if (i==1)fprintf(stdout, ",");
4580                 fprintf (stdout, "%s", TS->name[a]);
4581                 i=1;
4582               }
4583           fprintf ( stdout, "),(");
4584           for (i=0,a=0; a<nseq; a++)
4585             if (ref_buf[a]=='0')
4586               {
4587                 if (i==1) fprintf ( stdout, ",");
4588                 fprintf (stdout, "%s", TS->name[a]);
4589                 i=1;
4590               }
4591
4592           fprintf (stdout, ")\n");
4593           sprintf ( ref_buf, "%s", buf);
4594           n=1;
4595         }
4596       else
4597         {
4598           n++;
4599         }
4600     }
4601   vfclose (fp);
4602
4603
4604   myexit (0);
4605 }
4606
4607 int treelist2splits_old ( Sequence *S, Sequence *TS)
4608 {
4609   NT_node *T;
4610   int n=0,nseq, a,c;
4611
4612   int *used;
4613
4614   char *split_file, *sorted_split_file;
4615   char *buf=NULL, *ref_buf=NULL;
4616   FILE *fp;
4617
4618   split_file=vtmpnam (NULL);
4619   sorted_split_file =vtmpnam (NULL);
4620
4621   n=S->nseq;
4622   used=vcalloc (n, sizeof (int));
4623
4624   T=read_tree_list (S);
4625   if (!TS)TS=tree2seq(T[0], NULL);
4626   nseq=TS->nseq;
4627   fp=vfopen (split_file, "w");
4628
4629   for ( a=0; a< S->nseq; a++)
4630     {
4631
4632       T[a]=prune_tree  (T[a], TS);
4633       T[a]=recode_tree (T[a], TS);
4634       display_leaf_nb (T[a], TS->nseq,fp, S->name[a]);
4635     }
4636   vfclose (fp);
4637   printf_system ("cp %s split_file::IGNORE_FAILURE::", split_file);myexit (0);
4638
4639   printf_system ( "cat %s | grep 1| sort > %s::IGNORE_FAILURE::", split_file, sorted_split_file);
4640
4641   vfopen (sorted_split_file, "r");
4642
4643   while ( (c=fgetc (fp))!=EOF)
4644     {
4645
4646       ungetc (c, fp);
4647       buf=vfgets (buf, fp);
4648       buf [strlen(buf)-1]='\0';
4649
4650       if ( ref_buf==NULL)
4651         {
4652           ref_buf=vcalloc (strlen (buf)+1, sizeof (char));
4653           sprintf ( ref_buf, "%s", buf);
4654           n=1;
4655         }
4656       else if ( !strm (buf, ref_buf))
4657         {
4658           int i;
4659           fprintf ( stdout, "%3d %s(", n, ref_buf);
4660           for (i=0,a=0; a<nseq; a++)
4661             if (ref_buf[a]=='1')
4662               {
4663                 if (i==1)fprintf(stdout, ",");
4664                 fprintf (stdout, "%s", TS->name[a]);
4665                 i=1;
4666               }
4667           fprintf ( stdout, "),(");
4668           for (i=0,a=0; a<nseq; a++)
4669             if (ref_buf[a]=='0')
4670               {
4671                 if (i==1) fprintf ( stdout, ",");
4672                 fprintf (stdout, "%s", TS->name[a]);
4673                 i=1;
4674               }
4675
4676           fprintf (stdout, ")\n");
4677           sprintf ( ref_buf, "%s", buf);
4678           n=1;
4679         }
4680       else
4681         {
4682           n++;
4683         }
4684     }
4685   vfclose (fp);
4686
4687
4688   myexit (0);
4689 }
4690
4691 NT_node *treelist2prune_treelist (Sequence *S, Sequence *TS, FILE *out)
4692 {
4693   NT_node *T;
4694   int a, b, c;
4695
4696   T=read_tree_list (S);
4697   T=vrealloc (T, (S->nseq+1)*sizeof (NT_node));
4698   for (b=0,a=0; a<S->nseq; a++)
4699     {
4700       T[a]=prune_tree  (T[a], TS);
4701       if (tree2nleaf(T[a])<TS->nseq)
4702         {
4703           ;
4704         }
4705       else
4706         {
4707           char *s;
4708           T[b]=T[a];
4709           T[b]=recode_tree (T[b], TS);
4710           sprintf ( S->name[b], "%s", S->name[a]);
4711           s=tree2string (T[a]);
4712           S->seq[b]=vrealloc (S->seq[b], (strlen (s)+1)*sizeof (char));
4713           sprintf (S->seq[b], "%s",s);
4714           sprintf (S->seq_comment[b], " NSPECIES: %d", TS->nseq);
4715           vfree (s);
4716
4717           b++;
4718         }
4719
4720     }
4721
4722   S->nseq=b;
4723   T[S->nseq]=NULL;
4724
4725   if (out)
4726     {
4727       for (a=0; a<S->nseq; a++)
4728         {
4729           print_tree (T[a], "newick", out);
4730         }
4731     }
4732   return T;
4733 }
4734 int** treelist2lti2 ( Sequence *S, Sequence *TS, int ngb, FILE *out);
4735 int treelist2frame (Sequence *S, Sequence *TS)
4736 {
4737   int n, a, b, c,d, **r, **order;
4738   Sequence *temp;
4739
4740   temp=duplicate_sequence (S);
4741   order= treelist2lti (temp, TS,0,stdout);
4742
4743   TS=reorder_seq_2 (TS, order, 0, TS->nseq);
4744   n=TS->nseq;
4745
4746   for (a=3; a<n; a++)
4747     {
4748       NT_node tree;
4749
4750       TS->nseq=a+1;
4751       temp=duplicate_sequence (S);
4752       r=treelist2groups (temp,TS, NULL, NULL);
4753       fprintf ( stdout, "\n>Tree_%d [%d %%]\n ", a+1,r[0][1]);
4754       tree=main_read_tree (temp->name[r[0][0]]);
4755       tree=prune_tree (tree, TS);
4756       print_tree (tree, "newick",stdout);
4757
4758       free_int (r, -1);
4759       free_sequence (temp,-1);
4760     }
4761   myexit (EXIT_SUCCESS);
4762 }
4763 int** treelist2lti2 ( Sequence *S, Sequence *TS, int ngb, FILE *out)
4764 {
4765   NT_node *T;
4766   int a,b, c, d, ****dist, i;
4767   int **score, **order;
4768
4769   score=declare_int (TS->nseq, 3);
4770   order=declare_int (TS->nseq, 2);
4771   vsrand (0);
4772
4773   for (a=0; a<50; a++)
4774     {
4775       Sequence *seq, *trees;
4776       int **r;
4777       trees=duplicate_sequence (S);
4778       seq=duplicate_sequence (TS);
4779       for (b=0; b<TS->nseq; b++){order[b][0]=b;order[b][1]=rand()%10000;}
4780       sort_int (order, 2, 1, 0, TS->nseq-1);
4781       seq=reorder_seq_2(seq, order, 0,5);
4782       r=treelist2groups (trees,seq, NULL, NULL);
4783
4784       for (b=0; b<5; b++)
4785         {
4786           score[order[b][0]][1]+=r[0][1];
4787           score[order[b][0]][2]++;
4788         }
4789       HERE ("Score=%d", r[0][1]);
4790       free_int (r, -1);
4791       free_sequence (seq, -1);
4792       free_sequence (trees, -1);
4793
4794     }
4795
4796   for ( a=0; a< TS->nseq; a++)
4797     {
4798       score[a][0]=a;
4799       HERE ("%s => %d [%d]",TS->name[a], score[a][1]/score[a][2], score[a][2]);
4800       score[a][1]/=(score[a][2])?score[a][2]:1;
4801     }
4802   sort_int_inv (score, 3, 1, 0, TS->nseq-1);
4803
4804   return score;
4805 }
4806
4807
4808 int** treelist2lti ( Sequence *S, Sequence *TS, int ngb, FILE *out)
4809 {
4810   NT_node *T;
4811   int a,b, c, d, ****dist, i;
4812   float score0=0, score1=0;
4813   int **result;
4814
4815
4816   i=S->nseq;
4817   T=treelist2prune_treelist (S, TS,NULL);
4818
4819   if (!ngb)ngb=TS->nseq*2;
4820   dist=vcalloc ( S->nseq, sizeof (int****));
4821   result=declare_int (TS->nseq, 2);
4822   for (a=0; a<TS->nseq; a++)
4823     {
4824       float score_seq=0;
4825       float n_seq=0;
4826       for (b=0; b<TS->nseq;b++)
4827         {
4828           float score_pair=0;
4829           float n_pair=0;
4830           for (c=0; c<S->nseq; c++)
4831             {
4832               if (!dist[c])dist[c]=tree2dist(T[c], TS, NULL);
4833               for (d=0; d<S->nseq; d++)
4834                 {
4835                   float score, d1, d2;
4836
4837                   if (!dist[d])dist[d]=tree2dist(T[d], TS, NULL);
4838                   d1=dist[c][0][a][b];
4839                   d2=dist[d][0][a][b];
4840                   score=FABS((d1-d2));
4841                   if (d1>ngb || d2>ngb);
4842                   else
4843                     {
4844                       score_seq+=score;
4845                       score_pair+=score;
4846                       n_seq++;
4847                       n_pair++;
4848                     }
4849                   // if (d1 && d2) HERE ("%d %d", (int)d1, (int)d2);
4850                 }
4851             }
4852           score_pair=(score_pair*100)/(float)n_pair;
4853           if (out)fprintf ( stdout, "\n>%-20s %-20s LTI: %7.3f [Kept %d Trees Out of %d] ", TS->name[a],TS->name[b], score_pair, S->nseq,i);
4854         }
4855
4856       score_seq=(score_seq*100)/n_seq;
4857       result[a][0]=a;
4858       result[a][1]=(int)(100*score_seq);
4859       if (out)fprintf ( stdout, "\n>%-20s %-20s LTI: %7.3f [Kept %d Trees Out of %d] ", TS->name[a],"*", score_seq, S->nseq, i);
4860     }
4861   sort_int (result,2,1,0, TS->nseq-1);
4862   return result;
4863 }
4864
4865
4866 int ***tree2dist (NT_node T, Sequence *S, int ***d)
4867 {
4868   int *l0, *r0,*l1, *r1, a, b;
4869
4870
4871   if (!T) return d;
4872   if (!S)S=tree2seq(T, NULL);
4873   if (!d)
4874     {
4875       d=declare_arrayN (3, sizeof (float),2, S->nseq, S->nseq);
4876       T=prune_tree(T, S);
4877       T=recode_tree (T, S);
4878     }
4879
4880   if (!T->left)return d;
4881   if (!T->right) return d;
4882
4883   l0=(T->left)->idist;
4884   r0=(T->right)->idist;
4885
4886   l1=(T->left)->ldist;
4887   r1=(T->right)->ldist;
4888
4889
4890
4891   for (a=0; a< S->nseq; a++)
4892     for (b=0; b<S->nseq; b++)
4893       {
4894         if (l0[a]>0 && r0[b]>0)d[0][a][b]=d[0][b][a]=l0[a]+r0[b];
4895         if (l0[a]>0 && r0[b]>0)d[1][a][b]=d[1][b][a]=l1[a]+r1[b];
4896       }
4897
4898   d=tree2dist (T->left, S, d);
4899   d=tree2dist (T->right, S, d);
4900
4901
4902   return d;
4903 }
4904
4905
4906
4907 int **tree2dist_split ( NT_node T, Sequence *S, int **dist)
4908 {
4909
4910   FILE *fp;
4911   int a, b, c, n=0;
4912   char *buf=NULL, **list=NULL, *split_file;
4913
4914
4915   if (!S)S=tree2seq(T, NULL);
4916
4917   T=prune_tree  (T, S);
4918   T=recode_tree (T, S);
4919
4920   split_file=vtmpnam (NULL);
4921   fp=vfopen (split_file, "w");
4922   display_code (T, S->nseq,fp);
4923   vfclose (fp);
4924
4925   list=declare_char (2*S->nseq, S->nseq+1);
4926   fp=vfopen (split_file, "r");
4927
4928   while ((buf=vfgets (buf,fp))!=NULL)
4929     {
4930       if (buf[0]=='1' || buf[0]=='0')sprintf (list[n++], "%s", buf);
4931     }
4932   vfclose (fp);
4933   dist=declare_int ( S->nseq, S->nseq);
4934   for (a=0; a< S->nseq; a++)
4935     for ( b=0; b<S->nseq; b++)
4936       for (c=0; c<n; c++)
4937         if (list[c][a]!=list[c][b])dist[a][b]++;
4938
4939
4940   return dist;
4941 }
4942
4943 int** treelist2groups (Sequence *S, Sequence *TS, char *star_node, FILE *out)
4944    {
4945    NT_node *T;
4946    int a, b, tot, n,i;
4947    int v;
4948    int *used;
4949    int ntop;
4950    int nsn;
4951    int cov=100;
4952    int **results;
4953
4954
4955    i=S->nseq;
4956    T=treelist2prune_treelist (S, TS,NULL);
4957    nsn=(star_node)?atoi(star_node):0;
4958
4959    results=declare_int (S->nseq+1, 2);
4960
4961    if (nsn)
4962      {
4963        for (a=0; a< S->nseq; a++)tree2star_nodes(T[a],nsn);
4964      }
4965
4966    used=vcalloc (S->nseq, sizeof (int));
4967    for (ntop=0,a=0; a<S->nseq; a++)
4968      {
4969
4970        if (used[a]==0)
4971          {
4972            ntop++;
4973            if (out)fprintf ( out, "\nTree %s:",S->name[a]);
4974            used[a]=1;
4975          }
4976        else continue;
4977        tot=1;
4978        for ( b=0; b<S->nseq; b++)
4979          {
4980            v=0;
4981
4982            v=(int)simple_tree_cmp (T[a], T[b], TS, 1);
4983            if ( v==100)
4984              {
4985                used[b]=1;
4986                used[a]++;
4987                if (out)fprintf (stdout," %s ", S->name[b]);
4988                tot++;
4989              }
4990          }
4991
4992        if (out)fprintf ( stdout, "__ N=%d\n", tot-1);
4993      }
4994
4995
4996    for (n=0,a=0; a<S->nseq; a++)
4997      {
4998        if ( used[a]>1)
4999          {
5000            if (out)fprintf ( out, "\n>%-15s %4d %6.2f TOPOLOGY_LIST\n", S->name[a], used[a]-1, (float)(((float)used[a]-1)*100/(float)S->nseq));
5001            if (out)print_tree (T[a], "newick_tree", out);
5002            results[n][0]=a;
5003            results[n][1]=((used[a]-1)*100)/i;
5004            n++;
5005          }
5006      }
5007
5008    for (a=0; a<S->nseq; a++) free_tree(T[a]);
5009    vfree (T);
5010
5011    if (out)fprintf ( stdout, "\nTotal Number of different topologies: %d\n", ntop);
5012    results[n][0]=-1;
5013    sort_int_inv (results,2,1,0, n-1);
5014    for (a=0; a<S->nseq; a++) free_tree(T[a]);
5015    vfree (T);
5016    return results;
5017    }
5018 float simple_tree_cmp (NT_node T1, NT_node T2,Sequence *S, int mode)
5019 {
5020   Tree_sim *TS1, *TS2;
5021   float t, w, l, n;
5022
5023   TS1=vcalloc (1, sizeof (Tree_sim));
5024   TS2=vcalloc (1, sizeof (Tree_sim));
5025
5026
5027   T1=recode_tree(T1, S);
5028   T2=recode_tree(T2, S);
5029
5030   n=new_compare_trees ( T1, T2, S->nseq, TS1);
5031   new_compare_trees ( T2, T1, S->nseq, TS2);
5032
5033
5034
5035   t=(TS1->uw+TS2->uw)*100/(TS1->max_uw+TS2->max_uw);
5036   w=(TS1->w+TS2->w)*100/(TS1->max_w+TS2->max_w);
5037   l=(TS1->d+TS2->d)*100/(TS1->max_d+TS2->max_d);
5038
5039   vfree (TS1); vfree (TS2);
5040   if ( mode ==1)return t;
5041   else if (mode ==2) return w;
5042   else  return l;
5043 }
5044 int treelist2n (NT_node *L)
5045 {
5046   int n=0;
5047   while (L[n])n++;
5048   return n;
5049 }
5050 int **treelist2avg_treecmp (NT_node *L, char *file)
5051 {
5052   int a, b, n;
5053   int **score;
5054
5055   if (file) L=read_tree_list (main_read_seq(file));
5056   n=treelist2n (L);
5057
5058   score=declare_int (n, 2);
5059   for (a=0; a<n; a++)score[a][0]=a;
5060
5061   for (a=0; a<n-1; a++)
5062     {
5063       output_completion (stderr,a,n,1, "Tree Cmp");
5064       for (b=a+1; b<n; b++)
5065         {
5066           Tree_sim *ts;
5067           ts=tree_cmp (L[a],L[b]);
5068           score[a][1]+=ts->uw;
5069           score[b][1]+=ts->uw;
5070           vfree (ts);
5071         }
5072     }
5073   sort_int_inv (score, 2, 1, 0, n-1);
5074   if (file)free_treelist(L);
5075   return score;
5076 }
5077
5078 int treelist_file2consense (char *tree_file, char *outtree, char *outfile)
5079 {
5080   static char *command;
5081   static char *tmp_outtree;
5082   static char *tmp_outfile;
5083   FILE *fp;
5084   int flag1=0;
5085   int flag2=0;
5086
5087   if (!command)
5088     {
5089       command=vtmpnam (NULL);
5090       tmp_outtree=vtmpnam (NULL);
5091       tmp_outfile=vtmpnam (NULL);
5092     }
5093   if (!check_program_is_installed ("consense",NULL,NULL,"www.phylip.com",NO_REPORT))return 0;
5094
5095   fp=vfopen (command, "w");fprintf ( fp, "%s\nY\n", tree_file);fclose (fp);
5096   if ( check_file_exists ("outtree")){flag1=1;printf_system ("mv outtree %s::IGNORE_FAILURE::", tmp_outtree);}
5097   if ( check_file_exists ("outfile")){flag2=1;printf_system ("mv outfile %s::IGNORE_FAILURE::", tmp_outfile);}
5098   printf_system ("consense <%s > /dev/null 2>/dev/null::IGNORE_FAILURE::", command);
5099
5100   if ( outtree)printf_system ("mv outtree %s::IGNORE_FAILURE::", outtree);
5101   remove ("outtree");
5102   if ( outfile)printf_system ("mv outfile %s::IGNORE_FAILURE::", outfile);
5103   remove ("outfile");
5104   if (flag1)printf_system ("mv %s outtree::IGNORE_FAILURE::", tmp_outtree);
5105   if (flag2)printf_system ("mv %s outfile::IGNORE_FAILURE::", tmp_outfile);
5106   return 1;
5107 }
5108
5109
5110 NT_node treelist2filtered_bootstrap ( NT_node *L,char *file, int **score, float t)
5111 {
5112   NT_node BT, *L2;
5113   int n,a;
5114
5115   if (t==1 || t==0 || !score)return treelist2bootstrap (L, file);
5116
5117   if (file)L=read_tree_list (main_read_seq(file));
5118
5119   n=treelist2n(L)*t;
5120
5121   if (n==0) return NULL;
5122
5123   L2=vcalloc ( n+1, sizeof (NT_node));
5124   for (a=0; a<n; a++)
5125     L2[a]=L[score[a][0]];
5126
5127   BT=treelist2bootstrap (L2, NULL);
5128
5129   vfree (L2);
5130   if (file)free_treelist(L);
5131   return BT;
5132 }
5133
5134 NT_node treelist2bootstrap ( NT_node *L, char *file)
5135 {
5136   char *outfile;
5137   NT_node T;
5138   FILE *fp;
5139
5140   if (!file)
5141     {
5142       file=vtmpnam (NULL);
5143       vfclose (print_tree_list (L,"newick", vfopen (file, "w")));
5144     }
5145
5146   outfile=vtmpnam (NULL);
5147
5148   printf_system( "msa2bootstrap.pl -i %s -o %s -input tree >/dev/null 2>/dev/null", file, outfile);
5149
5150   T=main_read_tree (outfile);
5151   T=tree_dist2normalized_tree_dist (T,treelist2n(L));
5152
5153
5154   return T;
5155 }
5156
5157
5158
5159 Sequence * treelist2seq (Sequence *S)
5160 {
5161   int a, b, c, n, i;
5162   char **name;
5163   NT_node *T;
5164   Sequence *TS;
5165   char *fname;
5166   FILE *fp;
5167
5168   name=vcalloc (1, sizeof (char*));
5169   fp=vfopen ((fname=vtmpnam (NULL)), "w");
5170
5171   T=read_tree_list (S);
5172   for (n=0,a=0; a< S->nseq; a++)
5173     {
5174       TS=tree2seq(T[a], NULL);
5175       for (b=0; b<TS->nseq; b++)
5176         {
5177           if ( (i=name_is_in_list (TS->name[b], name, n, 100))==-1)
5178             {
5179               name[n]=vcalloc (100, sizeof (int));
5180               sprintf ( name[n], "%s", TS->name[b]);
5181               n++;
5182               name=vrealloc (name, (n+1)*sizeof (char*));
5183               fprintf ( fp, ">%s\n", TS->name[b]);
5184             }
5185         }
5186       free_sequence(TS, TS->nseq);
5187       free_tree (T[a]);
5188     }
5189
5190   vfclose (fp);
5191   vfree (T);
5192   return get_fasta_sequence (fname, NULL);
5193 }
5194
5195
5196 Sequence * treelist2sub_seq ( Sequence *S, int f)
5197 {
5198   NT_node *T;
5199   int a,b,c, s, i, n, maxnseq, tot;
5200   int **count, **grid;
5201   char *fname;
5202   Sequence *FS, *TS;
5203   FILE *fp;
5204   if (!f)return treelist2seq(S);
5205
5206
5207   //keep as many taxons as possible so that f% of the trees are kept
5208   //1: count the frequency of each taxon
5209
5210   FS=treelist2seq (S);
5211   maxnseq=FS->nseq;
5212
5213   count=declare_int (maxnseq, 3);
5214   grid=declare_int (S->nseq,maxnseq+1);
5215   T=read_tree_list (S);
5216
5217
5218
5219   for (a=0; a<FS->nseq; a++){count[a][0]=a;count[a][2]=1;}
5220   for (n=0,a=0; a< S->nseq; a++)
5221     {
5222       TS=tree2seq(T[a], NULL);
5223       for (b=0; b<TS->nseq; b++)
5224         {
5225           i=name_is_in_list (TS->name[b], FS->name, FS->nseq, 100);
5226           if ( i==-1){myexit (EXIT_FAILURE);}
5227           count[i][1]++;
5228           grid[a][i]=1;
5229         }
5230       free_sequence(TS, TS->nseq);
5231       free_tree (T[a]);
5232     }
5233   vfree (T);
5234   sort_int ( count,3,1, 0, maxnseq-1);
5235
5236   for (a=0; a<maxnseq; a++)
5237     {
5238       count[a][2]=0;
5239       for ( b=0; b< S->nseq; b++)grid[b][maxnseq]=1;//prepare to keep everything
5240       for ( tot=S->nseq, b=0; b< S->nseq; b++)
5241         {
5242           for (c=0; c<maxnseq; c++)
5243             {
5244               s=count[c][0];
5245               if (count[c][2] && !grid[b][s])
5246                 {
5247                   grid[b][maxnseq]=0;
5248                   tot--;
5249                   break;
5250                 }
5251             }
5252         }
5253       tot=(tot*100)/S->nseq;
5254       if ( tot>=f)break;
5255     }
5256   if (tot<f)return NULL;
5257
5258   fname=vtmpnam (NULL);
5259   fp=vfopen (fname, "w");
5260   for (a=0; a<maxnseq; a++)
5261     {
5262       if (count[a][2])
5263         {
5264           fprintf ( fp, ">%s LIMIT: %d %%\n", FS->name[count[a][0]], f);
5265
5266         }
5267     }
5268   vfclose (fp);
5269   free_int (grid, -1); free_int (count, -1);
5270   free_sequence (FS, FS->nseq);
5271
5272   return get_fasta_sequence (fname, NULL);
5273 }
5274 /******************************COPYRIGHT NOTICE*******************************/
5275 /*© Centro de Regulacio Genomica */
5276 /*and */
5277 /*Cedric Notredame */
5278 /*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
5279 /*All rights reserved.*/
5280 /*This file is part of T-COFFEE.*/
5281 /**/
5282 /*    T-COFFEE is free software; you can redistribute it and/or modify*/
5283 /*    it under the terms of the GNU General Public License as published by*/
5284 /*    the Free Software Foundation; either version 2 of the License, or*/
5285 /*    (at your option) any later version.*/
5286 /**/
5287 /*    T-COFFEE is distributed in the hope that it will be useful,*/
5288 /*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
5289 /*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
5290 /*    GNU General Public License for more details.*/
5291 /**/
5292 /*    You should have received a copy of the GNU General Public License*/
5293 /*    along with Foobar; if not, write to the Free Software*/
5294 /*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
5295 /*...............................................                                                                                      |*/
5296 /*  If you need some more information*/
5297 /*  cedric.notredame@europe.com*/
5298 /*...............................................                                                                                                                                     |*/
5299 /**/
5300 /**/
5301 /*      */
5302 /******************************COPYRIGHT NOTICE*******************************/