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