Makefile

@@ -21,7 +21,7 @@
 CONFIG ?= fltk-config
 PDF_PS_FLAGS ?= -DNO_PDF
 
-EXTRALIBS = `$(CONFIG) --use-images --ldflags`
+EXTRALIBS = `$(CONFIG) --use-images --ldflags` -ldl -lpthread -lz
 
 CC = gcc
 CXX = g++
@@ -38,7 +38,7 @@ OBJECTS = seaview.o custom.o use_mase_files.o regions.o load_seq.o align.o xfmat
 comlines.o resource.o nexus.o \
 viewasprots.o racnuc_fetch.o concatenate.o statistics.o \
 trees.o treedraw.o addbootstrap.o least_squares_brl.o \
-pseudoterminal.o unrooted.o pdf_or_ps.o svg.o threads.o Fl_Native_File_Chooser_FLTK.o
+pseudoterminal.o unrooted.o pdf_or_ps.o svg.o threads.o tbe.o Fl_Native_File_Chooser_FLTK.o
  
 COBJECTS = raa_acnuc.o parser.o md5.o zsockr.o misc_acnuc.o dnapars.o protpars.o seq.o phylip.o lwl.o bionj.o phyml_util.o
 

README

+           HOW TO COMPILE SEAVIEW ON UNIX/LINUX SYSTEMS
+
+The seaview program requires the FLTK library for its graphical user interface.
+Version 1.3.0 or later of FLTK is necessary.
+
+The FLTK library is available as packages for a number of Linux distributions.
+If this package is installed on your system, go to step II.
+If no such package is available for your system, or if version 1.3.0 is not available,
+or if you don't want to use such package, you can compile the FLTK library as follows.
+
+I COMPILING THE FLTK LIBRARY FROM SOURCE
+========================================
+
+The source code of the FLTK library is available at the "Download" item of http://www.fltk.org/.
+- Get the stable release source archive, currently fltk-1.3.2-source.tar.gz
+- Compile the FLTK library on your system:
+   cd fltk-1.3.2
+   ./configure
+   make
+It's not necessary to do "make install".
+
+II COMPILING THE SEAVIEW PROGRAM
+================================
+
+First unpack the seaview source archive (seaview.tar.gz).
+
+IIa Preparing the seaview makefile
+
+If the FLTK library is installed on your system (i.e., include files in /usr/include/FL,
+and libraries at /usr/lib/libfltk*), go to step IIb.
+Otherwise, edit Makefile from the seaview source directory which contains these 5 lines:
+FLTK = /usr/include
+#uncomment the next 3 lines to use custom installed FLTK library
+#FLTK = $(HOME)/fltk-1.3.2
+#IFLTK = -I$(FLTK)
+#LFLTK = -L$(FLTK)/lib
+
+Transform these lines into:
+#FLTK = /usr/include
+#uncomment the next 3 lines to use custom installed FLTK library
+FLTK = $(HOME)/fltk-1.3.2
+IFLTK = -I$(FLTK)
+LFLTK = -L$(FLTK)/lib
+and also change the FLTK = $(HOME)/fltk-1.3.2 line so it gives the location in your system
+of the top of the FLTK source tree.
+
+IIb Compile/link the seaview program
+
+In the seaview source directory, run
+make
+
+III Install seaview external programs
+=====================================
+
+Seaview drives 4 external programs (clustalo, muscle, Gblocks, phyml). Linux executables
+of these are available from the linux items of http://pbil.univ-lyon1.fr/software/seaview.html
+They are also available as packages of most Linux distributions.
+
+These programs should be put in your PATH or in the same directory as the seaview
+executable for seaview to drive them.
+
+By default seaview searches for the phyml executable under the names PhyML-3.1_linux32
+or PhyML-3.1_linux64. If they are installed under other names on your system, seaview 
+will ask you for their name and location which are remembered for subsequent runs.
+

debian/changelog

+seaview (1:4.7-1) UNRELEASED; urgency=medium
+
+  * New upstream version
+  * Drop patch since upstream has fixed Makefile
+  * Re-add example which was lost by accident
+
+ -- Andreas Tille <tille@debian.org>  Wed, 09 May 2018 14:09:47 +0200
+
 seaview (1:4.6.5-1) unstable; urgency=medium
 
   * New upstream version

debian/patches/re-add_missing_ldlibs.patch

-Description: Add ldflags used in Makefile of last version to build successfully
-Author: Andreas Tille <tille@debian.org>
-Last-Update: Fri, 27 Apr 2018 17:20:34 +0200
-
---- a/Makefile
-+++ b/Makefile
-@@ -21,7 +21,7 @@
- CONFIG ?= fltk-config
- PDF_PS_FLAGS ?= -DNO_PDF
- 
--EXTRALIBS = `$(CONFIG) --use-images --ldflags`
-+EXTRALIBS = `$(CONFIG) --use-images --ldflags` -ldl -lpthread -lz
- 
- CC = gcc
- CXX = g++

debian/patches/series

-re-add_missing_ldlibs.patch

debian/seaview.examples

 debian/example_files/*
-# example.nxs
+example.nxs

seaview.appdata.xml

+<?xml version="1.0" encoding="UTF-8"?>
+<!-- Copyright 2015 Manolo Gouy <manolo.gouy_at_univ-lyon1.fr> -->
+<component type="desktop">
+ <id>seaview.desktop</id>
+ <metadata_license>CC0-1.0</metadata_license>
+ <project_license>GPL-3.0</project_license>
+ <name>SeaView</name>
+ <summary>GUI for multiple sequence alignment and molecular phylogeny</summary>
+ <description>
+<p>SeaView is a multiplatform, graphical user interface for multiple sequence alignment and molecular phylogeny.</p>
+<ul>
+	<li>SeaView reads and writes various file formats (NEXUS, MSF, CLUSTAL, FASTA, PHYLIP,
+	MASE, Newick) of DNA and protein sequences and of phylogenetic trees.</li>
+	<li>SeaView drives programs muscle or Clustal-Omega for multiple sequence alignment, 
+	and also allows to use any external alignment algorithm able to read and write FASTA-formatted files.</li>
+	<li>Seaview drives the Gblocks program to select blocks of evolutionarily conserved sites.</li>
+	<li>SeaView computes phylogenetic trees by:</li>
+	<li>- parsimony, using PHYLIP's dnapars/protpars algorithm,</li>
+	<li>- distance, with NJ or BioNJ algorithms on a variety of evolutionary distances,</li>
+	<li>- maximum likelihood, driving program PhyML 3.1.</li>
+	<li>SeaView prints and draws phylogenetic trees on screen, SVG, PDF or PostScript files.</li>
+	<li>SeaView allows to download sequences from EMBL/GenBank/UniProt using the Internet.</li>
+</ul>
+ </description>
+ <screenshots>
+  <screenshot type="default">
+   <image>http://doua.prabi.fr/binaries/seaview4.png</image>
+   <caption>The main alignment window used in the color-by-codon mode</caption>
+  </screenshot>
+  <screenshot>
+   <image>http://doua.prabi.fr/binaries/seaview-tree.png</image>
+   <caption>The tree window with some zoom-in applied</caption>
+  </screenshot>
+  <screenshot>
+   <image>  http://doua.prabi.fr/software/seaview_data/PhyML-dialog.png</image>
+   <caption>Dialog window to perform Maximum-Likelihood tree-building</caption>
+  </screenshot>
+ </screenshots>
+ <url type="homepage">http://doua.prabi.fr/software/seaview</url>
+ <updatecontact>manolo.gouy_at_univ-lyon1.fr</updatecontact>
+ 
+ <developer_name>Manolo Gouy</developer_name>
+
+  <keywords>
+  <keyword>molecular phylogeny</keyword>
+  <keyword>phylogenetic tree</keyword>
+  <keyword>sequence alignment</keyword>
+  <keyword>parsimony</keyword>
+  <keyword>neighbor-joining/keyword>
+  <keyword>maximum likelihood</keyword>
+  </keywords>
+
+</component>

seaview.html

@@ -44,7 +44,7 @@ Molecular Biology and Evolution 27(2):221-224.
 SEAVIEW and PHYLO_WIN: two graphic tools for sequence alignment and molecular phylogeny.</a> 
 Comput. Appl. Biosci., 12:543-548.
 <p>
-Version 4.6.5
+Version 4.7
 
 <p>Binaries and full source code available from <a href=http://doua.prabi.fr/software/seaview>http://doua.prabi.fr/software/seaview</a>
 
@@ -56,16 +56,17 @@ Version 4.6.5
 <p>Seaview drives the Muscle, Clustal Omega, Gblocks, and PhyML programs and uses code from the PHYLIP
 package for parsimony. Please quote:
 <ul><li>Edgar R.C. (2004) <a href=http://dx.doi.org/10.1093/nar/gkh340>MUSCLE: multiple sequence alignment with high accuracy and 
-    high throughput.</a> Nucleic Acids Res. 32(5):1792-1797.
+    high throughput.</a> <i>Nucleic Acids Res.</i> 32(5):1792-1797.
 <li>Sievers F. et al. (2011) <a href=http://dx.doi.org/10.1038/msb.2011.75>Fast, scalable generation of high-quality protein 
 multiple sequence alignments using Clustal Omega.</a>
-    Molecular Systems Biology 7:539.
+    <i>Molecular Systems Biology</i> 7:539.
 <li>Castresana J. (2000) <a href=http://mbe.oxfordjournals.org/content/17/4/540.full>Selection of conserved blocks from multiple alignments for their 
-  use in phylogenetic analysis.</a> Molecular Biology and Evolution 17(4):540-552. 
+  use in phylogenetic analysis.</a> <i>Molecular Biology and Evolution</i> 17(4):540-552.
 <li>Guindon S., Dufayard J.F., Lefort V., Anisimova M., Hordijk W., Gascuel O. (2010) 
   <a href=http://dx.doi.org/10.1093/sysbio/syq010>New Algorithms and Methods to Estimate Maximum-Likelihood 
-  Phylogenies: Assessing the Performance of PhyML 3.0</a>. Systematic Biology 59(3):307-321.
+  Phylogenies: Assessing the Performance of PhyML 3.0</a>. <i>Systematic Biology</i> 59(3):307-321.
 <li>Felsenstein J. (2013) <a href=http://evolution.genetics.washington.edu/phylip.html>PHYLIP</a> version 3.696.
+<li>Lemoine F., Domelevo Entfellner J.-B., Wilkinson E., Correia D., Dávila Felipe M., De Oliveira T. & Gascuel O. (2018) <a href=http://dx.doi.org/10.1038/s41586-018-0043-0>Renewing Felsenstein’s phylogenetic bootstrap in the era of big data.</a> <i>Nature</i> 556:452-456.
 </ul>
 Seaview uses the <a href=http://www.fltk.org>FLTK</a> project for its user interface.
 
@@ -592,6 +593,7 @@ The weights of all equally parsimonious trees add up to 1 for each bootstrap rep
   <br><u> ignore all gap sites:</u> if on, all gap-containing sites are excluded from analysis;
       if off, not all sequence pairs use the same set of sites for computation of distances.
   <br><u> Bootstrap:</u> performs bootstrap evaluation of clade statistical support (can be interrupted).
+  Bootstrap supports can be optionally computed with the <a href=https://doi.org/10.1038/s41586-018-0043-0>Transfer Bootstrap Expectation</a> method.
   Furthermore, if "Show bootstrap trees" is checked,
   individual bootstrap trees will appear in a second tree window, in addition to the main distance tree window.
   <br><u> User tree:</u> computes least squares branch lengths for selected user tree topology.
@@ -606,8 +608,17 @@ Phylogenies: Assessing the Performance of PhyML 3.0</a>.
 Systematic Biology 59(3):307-321. 
  <br><u>   Model:</u> select one among a variety of evolutionary models.
  <br><u>   Branch support:</u> can be omitted (None) or estimated either by the approximate likelihood 
-       ratio test approach (aLRT) or by bootstrap. In the last case, if "Show bootstrap trees" is checked,
-       individual bootstrap trees will appear in a second tree window, in addition to the main PhyML tree window.
+ ratio test approach (aLRT) or by bootstrap. The bootstrap branch support can be computed by
+ two methods:
+ <ul><li> The 'standard' method (Felsenstein's) counts the fraction of replicates containing each internal
+   branch of the reference tree.
+   <li> The 'Transfer Bootstrap Expectation' method (TBE) proposed by Lemoine <i>et al.</i>:
+     <a href=https://doi.org/10.1038/s41586-018-0043-0>Renewing Felsenstein’s phylogenetic bootstrap in the era of big data</a>.
+     </ul>
+ If "Transfer Bootstrap Expectation method" is checked, both methods for computing
+ bootstrap supports are used and displayed in two tree windows.
+ With both bootstrap methods, if "Show bootstrap trees" is checked,
+ individual bootstrap trees will appear in yet another tree window.
  <br><u> Nucleotide/Amino acid equilibrium frequencies:</u> Residue frequencies can be estimated by counting
  their frequencies in the sequence alignment ("Empirical" option). Alternatively, they can be optimized by 
  maximum likelihood ("Optimized" option for DNA data) or fixed to a set of model-given values ("Model-given"

tbe.cxx

+#include "tbe.h"
+#include <assert.h>
+
+extern const char *preptree(TBE_Tree *fd_nj_plot);
+extern char *ecrit_arbre_parenth_unrooted(TBE_Tree *fd_nj_plot, TBE_Node *root);
+
+char* TBE_Compute(char *reftree,char* boottrees){
+  TBE_Branch *boot_branches;
+  TBE_Branch *ref_branches;
+
+  // Parsing reference tree
+  TBE_Tree *fd_reftree = (TBE_Tree *)calloc(1, sizeof(TBE_Tree));
+  fd_reftree->trees = reftree;
+  fd_reftree->current_tree = reftree;
+  const char *mess = preptree(fd_reftree);
+  //TBE_PrintTree(fd_reftree);
+
+  // Parsing bootstrap trees
+  TBE_Tree *fd_boot = (TBE_Tree *)calloc(1, sizeof(TBE_Tree));
+  fd_boot->trees = boottrees;
+  fd_boot->current_tree = boottrees;
+  fd_boot->rank = 0;
+
+  ref_branches = TBE_Assign_Branch_Ids(fd_reftree);
+  
+  // For each bootstrap tree
+  char *p = fd_boot->trees;
+  int nboot = 0;
+  int i;
+  TBE_Branch *cur_edge;
+  short unsigned** i_matrix;
+  short unsigned** c_matrix;
+  short unsigned** hamming;
+  short unsigned* min_dist;
+  short unsigned* min_dist_edge;
+  int* cluster_sizes;
+  char*  newtree;
+
+  while (p) {
+    // We compare the reference tree with the current bootstrap tree
+    char *oldtrees = fd_boot->trees;
+    fd_boot->trees = NULL;
+    char *oldname = fd_boot->tree_name;
+    fd_boot->tree_name = NULL;
+    free_tree(fd_boot);
+    fd_boot->trees = oldtrees;
+    fd_boot->tree_name = oldname;
+    mess = preptree(fd_boot);
+    //TBE_PrintTree(fd_boot);
+
+    if (mess) {
+      fl_alert("%s", mess);
+      return NULL;
+    }
+    
+    boot_branches = TBE_Assign_Branch_Ids(fd_boot);
+    
+    Alloc_TBE_Matrices(fd_reftree->notu, &i_matrix, &c_matrix, &hamming, &min_dist, &min_dist_edge, &cluster_sizes);
+    
+    Update_All_IC_Ref_Tree(fd_reftree, fd_boot,
+			   i_matrix, c_matrix,
+			   cluster_sizes,
+			   ref_branches, boot_branches);
+    Update_All_IC_Boot_Tree(fd_reftree, fd_boot,
+			    i_matrix, c_matrix,
+			    hamming, min_dist, min_dist_edge,
+			    cluster_sizes,
+			    ref_branches, boot_branches);
+    
+    for(i=0; i<2*fd_reftree->notu-2; i++){
+      cur_edge = &(ref_branches[i]);
+      cur_edge->score+=min_dist[cur_edge->id];
+      cur_edge->p=cluster_sizes[cur_edge->id];
+      if(cur_edge->p>(fd_reftree->notu/2)){
+	cur_edge->p = fd_reftree->notu - cur_edge->p;
+      }
+    }
+    Free_TBE_Matrices(fd_reftree->notu, &i_matrix, &c_matrix, &hamming, &min_dist, &min_dist_edge, &cluster_sizes);
+    free(boot_branches);
+    nboot++;
+    p = get_next_tree(fd_boot);
+  }
+  free(ref_branches);
+  newtree = (char *)malloc(500000);
+  TBE_PrintTree(newtree,fd_reftree,ref_branches, nboot);
+  free(reftree);
+  return newtree;
+}
+
+
+//////
+
+/* Matrices used in transfer bootstrap computation (tbe.c) */
+void Alloc_TBE_Matrices(int n_otu, short unsigned*** i_matrix,
+			short unsigned*** c_matrix, short unsigned*** hamming,
+			short unsigned** min_dist, short unsigned**  min_dist_edge,
+			int** cluster_sizes){
+  int i;
+  int nb_edges = 2*n_otu-2;
+  (*min_dist) = (short unsigned*) malloc(nb_edges*sizeof(short unsigned)); /* array of min Hamming distances */
+  (*min_dist_edge) = (short unsigned*) malloc(nb_edges*sizeof(short unsigned)); /* array of edge ids corresponding to min Hamming distances */
+  (*cluster_sizes) = (int*) malloc(nb_edges*sizeof(int)); /* array of sizes of clusters associated to each branch (in the post order traversal) */
+  (*c_matrix) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of cardinals of complements */
+  (*i_matrix) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of cardinals of intersections */
+  (*hamming) = (short unsigned**) malloc(nb_edges*sizeof(short unsigned*)); /* matrix of Hamming distances */
+  for (i=0; i<nb_edges; i++){
+    (*c_matrix)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
+    (*i_matrix)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
+    (*hamming)[i] = (short unsigned*) malloc(nb_edges*sizeof(short unsigned));
+    (*min_dist)[i] = n_otu; /* initialization to the nb of taxa */
+    (*cluster_sizes)[i] = 0;
+  }
+}
+
+void Free_TBE_Matrices(int n_otu,  short unsigned*** i_matrix, short unsigned*** c_matrix,
+		       short unsigned*** hamming, short unsigned** min_dist,
+		       short unsigned**  min_dist_edge, int** cluster_sizes){
+  int i;
+  int nb_edges = 2*n_otu-2;
+  for (i=0; i<nb_edges; i++) {
+    free((*c_matrix)[i]);
+    free((*i_matrix)[i]);
+    free((*hamming)[i]);
+  }
+  free((*c_matrix));
+  free((*i_matrix));
+  free((*hamming));
+  free((*min_dist));
+  free((*min_dist_edge));
+  free((*cluster_sizes));
+}
+
+///// Distance computation
+
+
+/* UNION AND INTERSECT CALCULATIONS (FOR THE TRANSFER METHOD) */
+void Update_IC_Ref_Tree(TBE_Tree *ref_tree, TBE_Node * orig, TBE_Node* target, TBE_Branch *my_br, TBE_Tree *boot_tree,
+			short unsigned** i_matrix, short unsigned** c_matrix, int* cluster_sizes,
+			TBE_Branch *ref_branches,
+			TBE_Branch *boot_branches){
+  /* this function does the post-order traversal (recursive from the pseudoroot to the leaves, updating knowledge for the subtrees)
+     of the reference tree, examining only leaves (terminal edges) of the bootstrap tree.
+     It sends a probe from the orig node to the target node (nodes in ref_tree), calculating I_ij and C_ij
+     (see Brehelin, Gascuel, Martin 2008). */
+  int j;
+  int edge_id; /* its id */
+  int next_edge_id;
+  TBE_Node *tip;
+  int boot_edge_id;
+  edge_id = my_br->id; /* all this is in ref_tree */
+
+  // A tip
+  if(TBE_Is_Taxon(target)) {
+    cluster_sizes[edge_id] = 1;
+    for (j=0; j < 2*boot_tree->notu-2; j++) { /* for all the terminal edges of boot_tree */
+      boot_edge_id = boot_branches[j].id;
+      // If not tip, continue
+      tip = NULL;
+      if(TBE_Is_Taxon(boot_branches[j].bouta))
+	tip = boot_branches[j].bouta;
+      if(TBE_Is_Taxon(boot_branches[j].boutb))
+	tip = boot_branches[j].boutb;
+      if(tip == NULL) continue;
+      /* we only want to scan terminal edges of boot_tree, where the right son is a leaf */
+      /* else we update all the I_ij and C_ij with i = edge_id */
+      if (strcmp(ref_tree->labels[target->rank],boot_tree->labels[tip->rank])) {
+	/* here the taxa are different */
+	i_matrix[edge_id][boot_edge_id] = 0;
+	c_matrix[edge_id][boot_edge_id] = 1;
+      } else {
+	/* same taxa here in T_ref and T_boot */
+	i_matrix[edge_id][boot_edge_id] = 1;
+	c_matrix[edge_id][boot_edge_id] = 0;
+      }
+    } /* end for on all edges of T_boot, for my_br being terminal */
+  } else {
+    cluster_sizes[edge_id] = 0;
+    /* now the case where my_br is not a terminal edge */
+    /* first initialise (zero) the cells we are going to update */
+    for (j=0; j < 2*boot_tree->notu-2; j++){
+      /* We initialize the i and c matrices for the edge edge_id with :
+       * 0 for i : because afterwards we do i[edge_id] = i[edge_id] || i[next_edge_id]
+       * 1 for c : because afterwards we do c[edge_id] = c[edge_id] && c[next_edge_id]
+       */
+      if(TBE_Is_Taxon(boot_branches[j].bouta) || TBE_Is_Taxon(boot_branches[j].boutb)){
+	boot_edge_id = boot_branches[j].id;
+ 	i_matrix[edge_id][boot_edge_id] = 0;
+	c_matrix[edge_id][boot_edge_id] = 1;
+      }
+    }
+    if(target->v1!=orig){
+      Update_IC_Ref_Tree(ref_tree, target, target->v1, &(ref_branches[target->br_id1]), boot_tree, i_matrix, c_matrix, cluster_sizes, ref_branches, boot_branches);
+      next_edge_id = target->br_id1;
+      cluster_sizes[edge_id] += cluster_sizes[next_edge_id];
+      for (j=0; j < 2*boot_tree->notu-2; j++) { /* for all the terminal edges of boot_tree */
+	boot_edge_id = boot_branches[j].id;
+	// If not a tip, continue
+	if(!TBE_Is_Taxon(boot_branches[j].bouta) && !TBE_Is_Taxon(boot_branches[j].boutb)) continue;
+	i_matrix[edge_id][boot_edge_id] = i_matrix[edge_id][boot_edge_id] || i_matrix[next_edge_id][boot_edge_id];
+	/* above is an OR between two integers, result is 0 or 1 */
+	c_matrix[edge_id][boot_edge_id] = c_matrix[edge_id][boot_edge_id] && c_matrix[next_edge_id][boot_edge_id];
+	/* above is an AND between two integers, result is 0 or 1 */
+      } /* end for j */
+    }
+    if(target->v2!=orig){
+      Update_IC_Ref_Tree(ref_tree, target, target->v2, &(ref_branches[target->br_id2]), boot_tree, i_matrix, c_matrix, cluster_sizes, ref_branches, boot_branches);
+      next_edge_id = target->br_id2;
+      cluster_sizes[edge_id] += cluster_sizes[next_edge_id];
+      for (j=0; j < 2*boot_tree->notu-2; j++) { /* for all the terminal edges of boot_tree */
+	boot_edge_id = boot_branches[j].id;
+	// If not a tip, continue
+	if(!TBE_Is_Taxon(boot_branches[j].bouta) && !TBE_Is_Taxon(boot_branches[j].boutb)) continue;
+	i_matrix[edge_id][boot_edge_id] = i_matrix[edge_id][boot_edge_id] || i_matrix[next_edge_id][boot_edge_id];
+	/* above is an OR between two integers, result is 0 or 1 */
+	c_matrix[edge_id][boot_edge_id] = c_matrix[edge_id][boot_edge_id] && c_matrix[next_edge_id][boot_edge_id];
+	/* above is an AND between two integers, result is 0 or 1 */
+      } /* end for j */
+    }
+    
+    if(target->v3!=orig){
+      Update_IC_Ref_Tree(ref_tree, target, target->v3, &(ref_branches[target->br_id3]), boot_tree, i_matrix, c_matrix, cluster_sizes, ref_branches, boot_branches);
+      next_edge_id = target->br_id3;
+      cluster_sizes[edge_id] += cluster_sizes[next_edge_id];
+      for (j=0; j < 2*boot_tree->notu-2; j++) { /* for all the terminal edges of boot_tree */
+	boot_edge_id = boot_branches[j].id;
+	// If not a tip, continue
+	if(!TBE_Is_Taxon(boot_branches[j].bouta) && !TBE_Is_Taxon(boot_branches[j].boutb)) continue;
+	i_matrix[edge_id][boot_edge_id] = i_matrix[edge_id][boot_edge_id] || i_matrix[next_edge_id][boot_edge_id];
+	/* above is an OR between two integers, result is 0 or 1 */
+	c_matrix[edge_id][boot_edge_id] = c_matrix[edge_id][boot_edge_id] && c_matrix[next_edge_id][boot_edge_id];
+	/* above is an AND between two integers, result is 0 or 1 */
+      } /* end for j */
+    }
+  } /* ending the case where my_br is an internal edge */
+} /* end update_i_c_post_order_ref_tree */
+
+void Update_All_IC_Ref_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree,
+			    short unsigned** i_matrix, short unsigned** c_matrix, int* cluster_sizes,
+			    TBE_Branch *ref_branches,
+			    TBE_Branch *boot_branches){
+  /* this function is the first step of the union and intersection calculations */
+  TBE_Node *root;
+  
+  root=ref_tree->racine;
+  Update_IC_Ref_Tree(ref_tree, root, root->v1, &(ref_branches[root->br_id1]), boot_tree, i_matrix, c_matrix, cluster_sizes, ref_branches, boot_branches);
+  Update_IC_Ref_Tree(ref_tree, root, root->v2, &(ref_branches[root->br_id2]), boot_tree, i_matrix, c_matrix, cluster_sizes, ref_branches, boot_branches);
+} /* end update_all_i_c_post_order_ref_tree */
+
+void Update_IC_Boot_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree, TBE_Node* orig, TBE_Node* target,
+			 TBE_Branch *my_br, short unsigned** i_matrix, short unsigned** c_matrix,
+			 short unsigned** hamming, short unsigned* min_dist,
+			 short unsigned* min_dist_edge, int* cluster_sizes,
+			 TBE_Branch *ref_branches,
+			 TBE_Branch *boot_branches){
+  /* here we implement the second part of the Brehelin/Gascuel/Martin algorithm:
+     post-order traversal of the bootstrap tree, and numerical recurrence. 
+     in this function, orig and target are nodes of boot_tree (aka T_boot).
+     min_dist is an array whose size is equal to the number of edges in T_ref.
+     It gives for each edge of T_ref its min distance to a split in T_boot. */
+  int i;
+  int boot_edge_id /* its id */, next_boot_edge_id /* id of descending branches. */;
+  int N = ref_tree->notu;
+  /* we first have to determine which is the direction of the edge (orig -> target and target -> orig) */
+  boot_edge_id = my_br->id; /* here this is an edge_id corresponding to T_boot */
+
+  // Not a taxon
+  if(!TBE_Is_Taxon(target)){
+    /* because nothing to do in the case where target is a leaf: intersection and union already ok. */
+    /* otherwise, keep on posttraversing in all other directions */
+    
+    /* first initialise (zero) the cells we are going to update */
+    for (i=0; i < 2*ref_tree->notu-2; i++) i_matrix[i][boot_edge_id] = c_matrix[i][boot_edge_id] = 0;
+    
+    if(target->v1 != orig){
+      next_boot_edge_id = target->br_id1;
+      Update_IC_Boot_Tree(ref_tree, boot_tree, target, target->v1, &(boot_branches[next_boot_edge_id]),
+			  i_matrix, c_matrix, hamming, min_dist, min_dist_edge, cluster_sizes, ref_branches, boot_branches);
+      for (i=0; i < 2*ref_tree->notu-2; i++) { /* for all the edges of ref_tree */ 
+	i_matrix[i][boot_edge_id] += i_matrix[i][next_boot_edge_id];
+	c_matrix[i][boot_edge_id] += c_matrix[i][next_boot_edge_id];
+      } /* end for i */
+    }
+    if(target->v2 != orig){
+      next_boot_edge_id = target->br_id2;
+      Update_IC_Boot_Tree(ref_tree, boot_tree, target, target->v2, &(boot_branches[next_boot_edge_id]),
+			  i_matrix, c_matrix, hamming, min_dist, min_dist_edge, cluster_sizes, ref_branches, boot_branches);
+      for (i=0; i < 2*ref_tree->notu-2; i++) { /* for all the edges of ref_tree */ 
+	i_matrix[i][boot_edge_id] += i_matrix[i][next_boot_edge_id];
+	c_matrix[i][boot_edge_id] += c_matrix[i][next_boot_edge_id];
+      } /* end for i */
+    }
+    if(target->v3 != orig){
+      next_boot_edge_id = target->br_id3;
+      Update_IC_Boot_Tree(ref_tree, boot_tree, target, target->v3, &(boot_branches[next_boot_edge_id]),
+			  i_matrix, c_matrix, hamming, min_dist, min_dist_edge, cluster_sizes, ref_branches, boot_branches);
+      for (i=0; i < 2*ref_tree->notu-2; i++) { /* for all the edges of ref_tree */ 
+	i_matrix[i][boot_edge_id] += i_matrix[i][next_boot_edge_id];
+	c_matrix[i][boot_edge_id] += c_matrix[i][next_boot_edge_id];
+      } /* end for i */
+    }
+  } /* end if target is not a leaf: the following loop is performed in all cases */
+
+  for (i=0; i< 2*ref_tree->notu-2; i++) { /* for all the edges of ref_tree */ 
+    /* at this point we can calculate in all cases (internal branch or not) the Hamming distance at [i][boot_edge_id], */
+    /* card of union minus card of intersection */
+    hamming[i][boot_edge_id] =
+      cluster_sizes[i] /* #taxa in the cluster i of T_ref */
+      + c_matrix[i][boot_edge_id] /* #taxa in cluster edge_id of T_boot BUT NOT in cluster i of T_ref */
+      - i_matrix[i][boot_edge_id]; /* #taxa in the intersection of the two clusters */
+    
+    /* Let's immediately calculate the right ditance, taking into account the fact that the true disance is min (dist, N-dist) */
+    if (hamming[i][boot_edge_id] > N/2 /* floor value */) hamming[i][boot_edge_id] = N - hamming[i][boot_edge_id];
+    /*   and update the min of all Hamming (TRANSFER) distances hamming[i][j] over all j */
+    if (hamming[i][boot_edge_id] < min_dist[i]){
+      min_dist[i] = hamming[i][boot_edge_id];
+      min_dist_edge[i] = boot_edge_id;
+    }
+  } /* end for on all edges of T_ref */
+} /* end update_i_c_post_order_boot_tree */
+
+
+void Update_All_IC_Boot_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree,
+			     short unsigned** i_matrix, short unsigned** c_matrix,
+			     short unsigned** hamming, short unsigned* min_dist,
+			     short unsigned* min_dist_edge, int* cluster_sizes,
+			     TBE_Branch *ref_branches,
+			     TBE_Branch *boot_branches){
+  /* this function is the second step of the union and intersection calculations */
+  int i;
+  TBE_Node *root;
+  root=boot_tree->racine;
+  Update_IC_Boot_Tree(ref_tree, boot_tree, root, root->v1, &(boot_branches[root->br_id1]), i_matrix, c_matrix, hamming, min_dist, min_dist_edge, cluster_sizes, ref_branches, boot_branches);
+  Update_IC_Boot_Tree(ref_tree, boot_tree, root, root->v2, &(boot_branches[root->br_id2]), i_matrix, c_matrix, hamming, min_dist, min_dist_edge, cluster_sizes, ref_branches, boot_branches);
+  
+  /* and then some checks to make sure everything went ok */
+  for(i=0; i<2*ref_tree->notu-2; i++) {
+    assert(min_dist[ref_branches[i].id] >= 0);
+    if(TBE_Is_Taxon(ref_branches[i].bouta) || TBE_Is_Taxon(ref_branches[i].boutb)){
+      assert(min_dist[ref_branches[i].id] == 0); /* any terminal edge should have an exact match in any bootstrap tree */
+    }
+  }
+} /* end update_all_i_c_post_order_boot_tree */
+
+
+// Assigns identifiers to all branches of the tree
+// And sets branch ids to each nodes br_id1 br_id2 and br_id3
+TBE_Branch* TBE_Assign_Branch_Ids(TBE_Tree *tree){
+  int id = -1;
+  TBE_Branch * branches = (TBE_Branch*) malloc((2*tree->notu-2)*sizeof(TBE_Branch));
+  for(int i=0;i<2*tree->notu-2;i++){
+    branches[i].score = 0;
+    branches[i].p = 0;
+  }
+  TBE_Init_Branches(tree->racine->v1, tree->racine, branches, &id);
+  TBE_Init_Branches(tree->racine->v2, tree->racine, branches, &id);
+  return branches;
+}
+
+
+void TBE_Init_Branches(TBE_Node *cur, TBE_Node *prev, TBE_Branch *branches, int *id){
+  (*id)++;
+  branches[*id].id = *id;
+  branches[*id].bouta = prev;
+  branches[*id].boutb = cur;
+
+  if(cur ==  prev->v1){
+    prev->br_id1 = *id;
+  }else if(cur == prev->v2){
+    prev->br_id2 = *id;
+  }else if(cur == prev->v3){
+    prev->br_id3 = *id;
+  }
+
+  if(prev ==  cur->v1){
+    cur->br_id1 = *id;
+  }else if(prev == cur->v2){
+    cur->br_id2 = *id;
+  }else if(prev == cur->v3){
+    cur->br_id3 = *id;
+  }
+  
+  if(cur->v1!=NULL && cur->v1!=prev){
+    TBE_Init_Branches(cur->v1, cur, branches,id);
+  }
+  if(cur->v2!=NULL && cur->v2!=prev){
+    TBE_Init_Branches(cur->v2, cur, branches,id);
+  }
+  if(cur->v3!=NULL && cur->v3!=prev){
+    TBE_Init_Branches(cur->v3, cur, branches,id);
+  }
+}
+
+int TBE_Is_Taxon(TBE_Node *n){
+  int nneigh = 0;
+  nneigh+=(n->v1!=NULL);
+  nneigh+=(n->v2!=NULL);
+  nneigh+=(n->v3!=NULL);
+  return nneigh==1;
+}
+
+void TBE_PrintTree(char* nwstring, TBE_Tree *t, TBE_Branch *branches, int nboot){
+  int pos=0;
+  nwstring[pos++]='(';
+  TBE_PrintTree_Recur(&pos,nwstring,t, t->racine->v1,t->racine, t->racine->br_id1, t->racine->l1, branches, nboot);
+  nwstring[pos++]=',';
+  TBE_PrintTree_Recur(&pos,nwstring,t, t->racine->v2,t->racine, t->racine->br_id2,t->racine->l2, branches, nboot);
+  nwstring[pos++]=')';
+  nwstring[pos++]=';';
+  nwstring[pos++]='\0';
+}
+
+void TBE_PrintTree_Recur(int *pos, char *nwstring, TBE_Tree *t, TBE_Node *cur,TBE_Node *orig, int cur_edge, double length, TBE_Branch *branches, int nboot){
+  if(TBE_Is_Taxon(cur)){
+    (*pos)+=snprintf(nwstring+(*pos),500000-(*pos),"%s",t->labels[cur->rank]);
+  }else{
+    nwstring[(*pos)++]='(';
+    int nb = 0;
+    if(cur->v1!=orig){
+      TBE_PrintTree_Recur(pos, nwstring, t, cur->v1, cur, cur->br_id1, cur->l1, branches, nboot);
+      nb++;
+    }
+    if(cur->v2!=orig){
+      if(nb>0)
+	nwstring[(*pos)++]=',';
+      TBE_PrintTree_Recur(pos, nwstring, t, cur->v2, cur, cur->br_id2, cur->l2, branches, nboot);
+      nb++;
+    }
+    if(cur->v3!=orig){
+      nwstring[(*pos)++]=',';
+      TBE_PrintTree_Recur(pos,nwstring, t, cur->v3, cur, cur->br_id2, cur->l3, branches, nboot);
+    }
+    nwstring[(*pos)++]=')';
+
+    // Compute and write support
+    if(branches[cur_edge].p>1){
+      double support = 100.0*(1.0 - (branches[cur_edge].score * 1.0 / nboot) / (branches[cur_edge].p-1));
+      (*pos)+=snprintf(nwstring+(*pos),500000-(*pos),"%d",int(support+0.5));
+    }
+  }
+  // Write length
+  (*pos)+=snprintf(nwstring+(*pos),500000-(*pos),":%f",length);
+}

tbe.h

+#ifndef TBE_H
+#define TBE_H
+
+#include "treedraw.h"
+/* Main function for computing TBE support
+   - reftree: 1 tree in newick string
+   - boottrees: several trees in a single newick string
+   - returns: newick string of the tree with supports
+*/
+char* TBE_Compute(char* reftree, char* boottrees);
+
+typedef struct { /* une branche definie par ses deux extremites */
+	struct noeud *bouta;
+	struct noeud *boutb;
+	char *br_label;
+	int id;
+	double score; // For TBE bootstrap support
+	int p; // Number of taxa on the light side og the bipartition (only defined after TBE_Bootstrap)
+} TBE_Branch;
+
+
+typedef struct FD_nj_plot TBE_Tree;
+typedef struct noeud TBE_Node;
+
+/* Matrices used in transfer bootstrap computation (tbe.c) */
+void Alloc_TBE_Matrices(int n_otu, short unsigned*** i_matrix,
+			short unsigned*** c_matrix,short unsigned*** hamming,
+			short unsigned** min_dist, short unsigned**  min_dist_edge,
+			int** cluster_sizes);
+
+void Free_TBE_Matrices(int n_otu,  short unsigned*** i_matrix, short unsigned*** c_matrix,
+		       short unsigned*** hamming, short unsigned** min_dist,
+		       short unsigned**  min_dist_edge, int** cluster_sizes);
+
+
+/* Transfer distance computation */
+void Update_All_IC_Ref_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree,
+			    short unsigned** i_matrix, short unsigned** c_matrix, int* cluster_sizes,
+			     TBE_Branch *ref_branches,
+			     TBE_Branch *boot_branches);
+void Update_IC_Ref_Tree(TBE_Tree *ref_tree, TBE_Node * orig, TBE_Node* target, TBE_Branch *my_br, TBE_Tree *boot_tree,
+			short unsigned** i_matrix, short unsigned** c_matrix, int* cluster_sizes,
+			     TBE_Branch *ref_branches,
+			     TBE_Branch *boot_branches);
+void Update_IC_Boot_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree, TBE_Node* orig, TBE_Node* target,
+			 TBE_Branch *my_br, short unsigned** i_matrix, short unsigned** c_matrix,
+			 short unsigned** hamming, short unsigned* min_dist,
+			 short unsigned* min_dist_edge, int* cluster_sizes,
+			     TBE_Branch *ref_branches,
+			     TBE_Branch *boot_branches);
+void Update_All_IC_Boot_Tree(TBE_Tree* ref_tree, TBE_Tree* boot_tree,
+			     short unsigned** i_matrix, short unsigned** c_matrix,
+			     short unsigned** hamming, short unsigned* min_dist,
+			     short unsigned* min_dist_edge, int* cluster_sizes,
+			     TBE_Branch *ref_branches,
+			     TBE_Branch *boot_branches);
+
+
+// Assigns identifiers to all TBE_branchs of the tree
+// And sets branch ids to each nodes br_id1 br_id2 and br_id3
+TBE_Branch* TBE_Assign_Branch_Ids(TBE_Tree *tree);
+void TBE_Init_Branches(TBE_Node *cur, TBE_Node *prev, TBE_Branch *branches, int *id);
+int TBE_Is_Taxon(TBE_Node *n);
+void TBE_PrintTree(char *nwstring, TBE_Tree *t, TBE_Branch *branches, int nboot);
+void TBE_PrintTree_Recur(int *pos, char *nwstring, TBE_Tree *t, TBE_Node *cur,TBE_Node *orig, int cur_edge, double length, TBE_Branch *branches, int nboot);
+#endif

treedraw.h

@@ -16,6 +16,9 @@ struct noeud {
 	struct noeud *v1,*v2,*v3;
 	char *nom;
 	int rank;
+	// Branches may be associated with id
+	// And this may be used in some cases
+	int br_id1, br_id2, br_id3;
 };
 
 struct nom {
@@ -130,5 +133,7 @@ struct FD_nj_plot {
 };
 
 extern Fl_Window *treedraw(char *tree, SEA_VIEW *view, const char *name, int from_tree_menu, int count = 1);
+char *get_next_tree(FD_nj_plot *fd_nj_plot);
+void free_tree(FD_nj_plot *fd_nj_plot);
 
 #endif // TREEDRAW_H

trees.cxx

@@ -3,6 +3,7 @@
 #endif
 #include "seaview.h"
 #include "treedraw.h"
+#include "tbe.h"
 #include <FL/Fl_Round_Button.H>
 #include <FL/Fl_Check_Button.H>
 #include <FL/Fl_Select_Browser.H>
@@ -43,7 +44,7 @@ void distance_method_dialog(SEA_VIEW *view);
 matrix *run_calc_dist_matrix(allseq *seqs, int distkind, int protein, char **lwlseqs, int in_bootstrap, char **p_err_mess, int no_gui);
 char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int use_bionj, int use_bootstrap,
 			 int replicates, Fl_Box *box, const char *distance_name, int using_kaks,
-			 const char* distancefname, const char *user_tree, int no_gui, void (*alert)(const char*, ...), int keep_b_trees );
+			  const char* distancefname, const char *user_tree, int no_gui, void (*alert)(const char*, ...), int keep_b_trees, int compute_tbe );
 char *bootstrap_reformat(char *tree, int replicates);
 matrix *Obs_Dist(allseq *data, model *mod);
 matrix *Kimura_p_Dist(allseq *data);
@@ -221,12 +222,17 @@ static void distance_callback(Fl_Widget *wgt, void *data)
 
 static void dist_b_cb(Fl_Widget *wgt, void *data)
 {
+  Fl_Button **buttons = (Fl_Button **)data;
   Fl_Button *me = (Fl_Button*)wgt;
-  Fl_Button *target = (Fl_Button*)data;
-  if (me->value()) target->activate();
-  else {
-    target->deactivate();
-    target->value(0);
+  for (int i=0; i<2; i++){
+    Fl_Button *target = buttons[i];
+    if (target){
+      if (me->value()) target->activate();
+      else {
+        target->deactivate();
+        target->value(0);
+      }
+    }
   }
 }
 
@@ -249,10 +255,11 @@ void distance_method_dialog(SEA_VIEW *view)
 	static Fl_Button *interrupt;
 	static Fl_Return_Button *go;
 	static int started, maxchoice;
-  static Fl_Check_Button *keep_b_trees;
+	static Fl_Check_Button *keep_b_trees;
+	static Fl_Check_Button *compute_tbe;
 	if(first) {
 		first = FALSE;
-		w = new Fl_Window(250, 240);
+		w = new Fl_Window(270, 265);
 		w->label("Distance analysis");
 		w->set_modal();
 		nj = new Fl_Round_Button(3, 3, 60, 20, "NJ");
@@ -264,7 +271,7 @@ void distance_method_dialog(SEA_VIEW *view)
 		distance = new Fl_Choice(90, nj->y() + nj->h() + 15, 100, 20, "Distance");
 		distance->align(FL_ALIGN_LEFT);
 		ignore_all_gaps = new Fl_Check_Button(160, distance->y() + distance->h() + 5, 50, 20, 
-						      "ignore all gap sites");
+						      "Ignore all gap sites");
 		ignore_all_gaps->value(1);
 		ignore_all_gaps->align(FL_ALIGN_LEFT);
 		int y = ignore_all_gaps->y() + ignore_all_gaps->h() + 25;
@@ -274,10 +281,19 @@ void distance_method_dialog(SEA_VIEW *view)
 		char tmp[10];
 		sprintf(tmp, "%d", def_replicates);
 		replicates->value(tmp);
-    keep_b_trees = new Fl_Check_Button(3, bootstrap->y() + bootstrap->h() + 5, 180, 20, "Show bootstrap trees");
-    keep_b_trees->deactivate();
-    bootstrap->callback(dist_b_cb, keep_b_trees);
-    bootstrap->when(FL_WHEN_CHANGED);
+
+		compute_tbe = new Fl_Check_Button(3, bootstrap->y() + bootstrap->h() + 5, 280, 20, "Transfer Bootstrap Expectation method");
+		compute_tbe->deactivate();
+
+		keep_b_trees = new Fl_Check_Button(3, compute_tbe->y() + compute_tbe->h() + 5, 180, 20, "Show bootstrap trees");
+		keep_b_trees->deactivate();
+
+    static Fl_Check_Button *buttons[2];
+		buttons[0]=compute_tbe;
+		buttons[1]=keep_b_trees;
+		bootstrap->callback(dist_b_cb, buttons);
+		bootstrap->when(FL_WHEN_CHANGED);
+		
 		box = new Fl_Box(3, keep_b_trees->y() + keep_b_trees->h() + 5, 150, 20, "");
 		usertree = new Fl_Round_Button(nj->x(), box->y() + box->h() + 15, 85, nj->h(), "User tree:");
 		usertree->type(FL_RADIO_BUTTON);
@@ -362,7 +378,7 @@ void distance_method_dialog(SEA_VIEW *view)
 	  char *tree = run_distance_method(view, distance_choice, ignore_all_gaps->value(), bionj->value(), 
 						  bootstrap->value(), def_replicates, box, distance->text(), using_kaks,
 			      savetofile->value()?"":NULL, usertree->value() ? view->trees[treechoice->value()] : NULL,
-			      false, fl_alert, keep_b_trees->value());
+					   false, fl_alert, keep_b_trees->value(), compute_tbe->value());
 	  const char *choice = bionj->value()? "BioNJ_tree" : "NJ_tree";
 	  const char *n = extract_filename(view->masename);
 	  char *title = NULL;
@@ -438,7 +454,7 @@ matrix *run_calc_dist_matrix(allseq *seqs, int distkind, int protein, char **lwl
 
 char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int use_bionj, int use_bootstrap,
 	int replicates, Fl_Box *box, const char *distance_name, int using_kaks, const char* distancefname,
-	const char *user_tree, int no_gui, void (*alert)(const char*, ...), int keep_b_trees )
+			  const char *user_tree, int no_gui, void (*alert)(const char*, ...), int keep_b_trees, int compute_tbe )
 {
 	char *display_tree, *bootstrap_trees = NULL, *err_mess, tree_label[100]="";
 	int total = 0, i, j, lbtrees;
@@ -511,7 +527,7 @@ char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int
 	phyml_mat->tree = Make_Tree_From_Scratch(phyml_seqs->n_otu, phyml_seqs);
 	Bionj(phyml_mat);
 	sprintf(tree_label, "%s %d sites %s", use_bionj ? "BioNJ" : "NJ",
-			phyml_seqs->clean_len, distance_name );
+		phyml_seqs->clean_len, distance_name);
 	if (use_bootstrap) {
 		char *full_tree = Write_Tree(phyml_mat->tree);
 		Free_Tree(phyml_mat->tree);
@@ -532,7 +548,7 @@ char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int
 				phyml_mat->tree = Make_Tree_From_Scratch(phyml_seqs->n_otu, phyml_seqs);
 				Bionj(phyml_mat);
 				char *one_replicate_tree = Write_Tree(phyml_mat->tree);
-        if (keep_b_trees) {
+        if (keep_b_trees || compute_tbe) {
           if (r == 0) {
             lbtrees = replicates * (strlen(one_replicate_tree) + 1);
             bootstrap_trees = (char *)malloc(lbtrees + 1);
@@ -564,7 +580,7 @@ char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int
 		}
 		if(total > 0) {
 			display_tree = finish_add_bootstrap(total);
-			sprintf(tree_label + strlen(tree_label), " %d repl.", total);
+			sprintf(tree_label + strlen(tree_label), " %d repl.%s", total, compute_tbe ? " (TBE)" : "" );
 			}
 		else display_tree = NULL;
 	}
@@ -581,12 +597,15 @@ char *run_distance_method(SEA_VIEW *view, int distkind, int remove_all_gaps, int
 	  if( !using_kaks) alert("Can't compute distances because sequences are too divergent.\nAre they aligned?");
 	  return NULL;
 	}
+	if (bootstrap_trees && compute_tbe) {
+	  display_tree=TBE_Compute(display_tree,bootstrap_trees);
+	}
 	//add label to tree
 	char *tree = (char *)malloc(strlen(tree_label) + strlen(display_tree) + 4);
 	sprintf(tree, "[%s] %s", tree_label, display_tree);
 	free(display_tree);
 	if (box) box->window()->hide();//necessary under X11 so tree is in foreground
-  if (bootstrap_trees) {
+  if (bootstrap_trees && keep_b_trees) {
     treedraw(bootstrap_trees, view, "bootstrap trees", FALSE);
   }
   return tree;
@@ -733,17 +752,21 @@ Fl_Round_Button *g_no; Fl_Choice *g_cats; Fl_Round_Button *a_est; Fl_Float_Input
 Fl_Round_Button *nni_b, *spr_b, *nni_spr_b;
 Fl_Round_Button *u_bionj, *u_menutree, *u_random;
 Fl_Choice *u_choice;
-Fl_Check_Button *u_optimize, *u_quiet, *b_keep_trees;
+Fl_Check_Button *u_optimize, *u_quiet, *b_keep_trees, *b_tbe;
 Fl_Int_Input *u_random_count;
 char *phyml_path;
 
 static void support_cb(Fl_Widget *wid, void *data)
 {
   Fl_Button *b = (Fl_Button*)data;
-  if (b->value()) b_keep_trees->activate();
+  if (b->value()) {
+    b_keep_trees->activate();
+    b_tbe->activate();
+  }
   else {
     b_keep_trees->deactivate();
     b_keep_trees->value(0);
+    b_tbe->deactivate();
   }
 }
 
@@ -831,7 +854,7 @@ void phyml_dialog(SEA_VIEW *view)
 		}
 	y += mymodel->h() + 25;
 	//bootstrap
-	Fl_Group *b_group = new Fl_Group(0, y, w->w(), 75, "Branch Support");
+	Fl_Group *b_group = new Fl_Group(0, y, w->w(), 95, "Branch Support");
 	b_group->box(FL_ROUNDED_BOX);
 	b_group->align(FL_ALIGN_TOP|FL_ALIGN_CENTER);
 	y += 5;
@@ -846,9 +869,14 @@ void phyml_dialog(SEA_VIEW *view)
 	b_count = new Fl_Int_Input(b_yes->x() + b_yes->w() + 5, y, 50, 20, "replicates");
 	b_count->align(FL_ALIGN_RIGHT);
 	b_count->static_value("100");
+  
   y += b_count->h();
+  b_tbe = new Fl_Check_Button(b_yes->x(), y, 265, 20, "Transfer Bootstrap Expectation method");
+  b_tbe->deactivate();
+  y += b_tbe->h();
   b_keep_trees = new Fl_Check_Button(b_yes->x(), y, 180, 20, "Show bootstrap trees");
   b_keep_trees->deactivate();
+
   b_alrt->callback(support_cb, b_yes);
   b_no->callback(support_cb, b_yes);
   b_yes->callback(support_cb, b_yes);
@@ -1219,7 +1247,7 @@ void run_phyml_callback(Fl_Widget *ob, void *data)
         fclose(in);
       }
     }
-    if (replicates > 0 && b_keep_trees->value()) {
+    if (replicates > 0 && (b_keep_trees->value() || b_tbe->value())) {
       sprintf(input, "%s.phy_phyml_boot_trees.txt", base_fname);
       in = fopen(input, "r");
       if (!in) {
@@ -1271,18 +1299,47 @@ void run_phyml_callback(Fl_Widget *ob, void *data)
       }
     }
 	  //clean tree internal labels
-	  if(replicates > 0) display_tree = bootstrap_reformat(display_tree, replicates);
+    char *display_tree_tbe = NULL;
+    if(replicates > 0) {
+      if (b_tbe->value()){
+        char *tmp = strdup(display_tree);
+        display_tree_tbe = TBE_Compute(tmp, bootstrap_trees);
+      }
+      display_tree = bootstrap_reformat(display_tree, replicates);
+      if(bootstrap_trees && !b_keep_trees->value()) {
+        free(bootstrap_trees);
+        bootstrap_trees = NULL;
+      }
+    }
 	  else if(replicates < 0) display_tree = alrt_reformat(display_tree);
+	  
 	  //assemble tree description
 	  char tree_label[100];
 	  sprintf(tree_label, "PhyML ln(L)=%.1f %d sites %s", logL, seqlen, modelname);
-	  if(replicates > 0) sprintf(tree_label + strlen(tree_label), " %d replic.", replicates);
+	  if(replicates > 0) {
+      sprintf(tree_label + strlen(tree_label), " %d replic.", replicates);
+	  }
 	  if(cats > 1) sprintf(tree_label + strlen(tree_label), " %d rate classes", cats);
 	  if(user_tree) strcat(tree_label, " User-tree");
 	  //add label to tree
 	  char *tree = (char *)malloc(strlen(tree_label) + strlen(display_tree) + 4);
 	  sprintf(tree, "[%s] %s", tree_label, display_tree);
 	  free(display_tree);
+    
+    if (replicates > 0 && b_tbe->value()) {
+      //assemble TBE tree description
+      char tree_label_tbe[100];
+      sprintf(tree_label_tbe, "PhyML ln(L)=%.1f %d sites %s", logL, seqlen, modelname);
+      sprintf(tree_label_tbe + strlen(tree_label_tbe), " %d replic. (TBE)", replicates);
+      if(cats > 1) sprintf(tree_label_tbe + strlen(tree_label_tbe), " %d rate classes", cats);
+      if(user_tree) strcat(tree_label_tbe, " User-tree");
+      //add label to tree
+      char *tree = (char *)malloc(strlen(tree_label_tbe) + strlen(display_tree_tbe) + 4);
+      sprintf(tree, "[%s] %s", tree_label_tbe, display_tree_tbe);
+      free(display_tree_tbe);
+      display_tree_tbe = tree;
+    }
+    
 	  if (is_view_valid(view)) {
 	    const char *n = extract_filename(view->masename);
 	    char *title = NULL;
@@ -1293,13 +1350,19 @@ void run_phyml_callback(Fl_Widget *ob, void *data)
 	      sprintf(title, "%s-%s", n, "PhyML_tree");
 	      if (q) *q = '.';
 	    }
-      Fl_Window *btreew, *treew;
-      if (bootstrap_trees) {
+      Fl_Window *btreew=0, *treew=0;
+      if (replicates > 0 && b_keep_trees->value()) {
         btreew = treedraw(bootstrap_trees, view, "bootstrap trees", FALSE);
+      }
+      if (replicates > 0 && b_tbe->value()) {
+        const char *p = n ? title:"PhyML_tree";
+        char *q = (char*)malloc(strlen(p) + 5);
+        sprintf(q, "%s_TBE", p);
+        treedraw(display_tree_tbe, view, q, FALSE);
       }
 	    treew = treedraw(tree, view, n ? title:"PhyML_tree", FALSE);
 #if !defined(__APPLE__)
-      if (bootstrap_trees) treew->position(btreew->x() + 30, btreew->y() + 40);
+      if (treew && btreew) treew->position(btreew->x() + 30, btreew->y() + 40);
 #endif
 	    delete[] title;
 	  }
@@ -2133,7 +2196,8 @@ void parsimony_dialog(SEA_VIEW *view)
 	w_btrepl->align(FL_ALIGN_RIGHT);
   Fl_Check_Button *keep_b_trees = new Fl_Check_Button(2, w_dobootstrap->y() + w_dobootstrap->h() + 5, 180, 20, "Show bootstrap trees");
   keep_b_trees->deactivate();
-  w_dobootstrap->callback(dist_b_cb, keep_b_trees);
+  static Fl_Check_Button *buttons[2] = {keep_b_trees, NULL};
+  w_dobootstrap->callback(dist_b_cb, buttons);
   w_dobootstrap->when(FL_WHEN_CHANGED);
 
 	Fl_Box *w_count = new Fl_Box(0, keep_b_trees->y() + keep_b_trees->h() + 5, w->w(), 20, NULL);
@@ -2744,7 +2808,7 @@ void command_line_phylogeny(int argc, char **argv)
     
     view->tot_trees = 0;
     tree = run_distance_method(view, distkind, isarg(argc, argv, "-nogaps"), !isarg(argc, argv, "-NJ"), replicates > 1, replicates, NULL,
-			       distname, distkind == Ka || distkind == Ks, argname(argc, argv, "-distance_matrix"), utree, true, Fl::fatal, false);
+			       distname, distkind == Ka || distkind == Ks, argname(argc, argv, "-distance_matrix"), utree, true, Fl::fatal, false, false);
     if (isarg(argc, argv, "-distance_matrix")) exit(0);
   }