ExaBayes User's Manual

After you have downloaded the appropriate package, you find all executables in the ./bin folder. Just execute these only using the -h flag and you will be given a help page. If the executables produce error messages (indicating that some library was not found), you either downloaded the wrong package or you have to compile from source.

Download and extract the source archive (sources are also included in all binary distributions).

Let's at first simply run a single chain for some time:

$ ./yggdrasil -f aln.phy -m DNA -n myRun -s $RANDOM

ExaBayes will print output that is separated into various sections:

1. after the header, you find a section (divided by '=') that re-iterates the alignment (number of unique patterns and type of partitions).
2. In the next section, ExaBayes lists the parameters to be integrated. For instance, the tree topology is considered a parameter. It also displays the (default) prior for parameters and initial values.
3. The \(3^{rd}\) section contains the proposals that are instantiated for integrating over the aforementioned parameters.

During the MCMC simulation, ExaBayes prints the log-likelihoods (lnl) of the (sole) chain. You'll find that after some time the lnl will reach a plateau. Stop the run after something like 50,000 generations (using Control-c).

Now examine the output files created by ExaBayes (we neglect the binary alignment file and the checkpoint files):

• ExaBayes_info.myRun
  Contains the same information also printed to the screen.
• ExaBayes_topologies.myRun.0
  Contains all sampled topologies in nexus format.
• ExaBayes_parameters.myRun.0
  Contains values sampled for all non-topological, non-branch length values.
• ExaBayes_diagnostics.myRun
  Contains chain diagnostics (e.g., acceptance ratios for all proposals or topological convergence in form of asdsf).

Now use the post-processing tools to examine the result. First, we create a consensus tree:

$ ./consense -f ExaBayes_topologies.myRun.0 -n myCons

Now, open a tree viewer of your choice (e.g., FigTree, Archaeopteryx or Dendroscope) and have a look at the consensus tree. If you ran just 50,000 generations, you will probably find the confidence in most branches is pretty low.

Let's inspect the \(50\%\) credible set of trees:

$ ./credibleSet -f ExaBayes_topologies.myRun.0 -n cred

The output file ExaBayes_credibleSet.tmp contains all sampled trees ordered by the frequency of their occurrence. You probably will find that no tree occurred more than once.

Finally, let's check, how well the parameters are sampled:

$ ./postProcParam   -f ExaBayes_parameters.myRun.0  -n params

Alternatively, you could also open the parameter file with Tracer and visualize the distributions. If you do not have Tracer installed, have a look at ExaBayes_parameterStatistics.params (spreadsheet tools like Excel are helpful). You'll find summary statistics for each parameter. Specifically, check out the effective sampling size (ESS) value for each parameter. Since samples in a chain are correlated, they are less informative, than if you had drawn the values independently from the original distribution. The ESS of samples indicates the number of samples your samples corresponds to, if they were drawn independently.

You will find that most ESS values are in the range between 30 and 80. This is not bad for the low number of generations, but to assure that each parameter has been sampled adequately, values should be \(>100\) or even better \(>200\). High ESS values indicate that your chain has explored this parameter sufficiently.

Now remove the files from the initial run (otherwise ExaBayes will complain). We will now conduct a proper analysis using several chains on a partitioned alignment.

You have to modify config.nex (remove the square brackets to uncomment). Uncomment numGens and set it to 100000 (or 1e5). Uncomment numRuns and set it to 4. This means that ExaBayes will conduct 4 independent analysis (starting in 4 different random trees). If all runs yield the same split frequencies (i.e., the same confidence in a split), we can be relatively sure that we have sampled the topology parameter sufficiently (while it is still possible that simply all 4 chains got stuck in the same local optimum).

Check out the partitions file (aln.part): it contains 4 partitions. By default, ExaBayes assumes that partitions have distinct branch lengths. Let's link all partitions into a single branch length parameter. To do so, add "brlens = (0-3)" in the params section.

Use this command line to start the analysis:

$ ./yggdrasil -f aln.phy -q aln.part  -n myRun -s $RANDOM -c config.nex

By default, ExaBayes will compute the average standard deviation of split frequencies (asdsf) every 5,000 generations and stops the analysis once it the asdsf is better than \(5\%\) (and once we have at least 100,000 generations for each chain).

The analysis may take a while. If you have a cluster (with for instance 16 cores) at your disposal, consider running the parallel version:

$ mpirun -np 16  ./exabayes -f aln.phy -q aln.part  -n myRun -s $RANDOM -c config.nex -R 4

With -R 4, ExaBayes runs all 4 chains in parallel. In an additional section, ExaBayes will inform you which chains and how many unique site patterns are assigned to each process.

You will notice after 100,000 generations, that the 4 chains converge rather slowly. Thus, enable Metropolis-Coupling to speed up the convergence. Set numCoupledChains to 2 and (for purely computational reasons in this example), reduce the numRuns to 2. Also, set parsimonyStart to true, such that your chains start from a parsimony tree instead of a random tree. Using parsimony trees as initial topologies saves you some time, however theoretically it also increases the probability that all your independent chains become stuck in the same local optimum and you obtain incorrect estimates for posterior probabilities.

If you have many cores, to run the two independent runs as well as the coupled chains in parallel:

$ mpirun -np 16 ./exabayes -f aln.phy -q aln.part -n myRun -s $RANDOM -c config.nex -R 2 -C 2

After \(\approx\) 150,000 generations, the ASDSF should have fallen below \(5\%\), which is acceptable. A look at the consensus tree reveals that this dataset in fact contains several low confidence branches. While in the first run the reason for low confidence branches was due to the low number of generations, we can be more certain now that that low confidence branches are a result of the phylogenetic signal in the alignment. Since you ran 2 independent analyses, you obtain 2 sets of output files (parameter and topologies). You can feed both files into the consensus tree (for both files the initial \(25\%\) of samples will be discarded).

$ ./consense -f ExaBayes_topologies.myRun.* -n myCons

Analyze both parameter files using

$ ./postProcParam -f ExaBayes_parameters.myRun.* -n params

Since we carried out two partitioned analyses, we have a larger number of parameters (where e.g., r{2}(C<->T) is the substitution probability of C to T in the third partition). We find that nearly all parameters have an ESS \(> 100\). The final column now contains the potential scale reduction factor (PSRF). It indicates, whether within-chain variance (of a parameter) is similar to between-chain variance. For this convergence statistic, values should be close to 1, but lower than 1.2 or 1.1 (probably the case after 150,000 generations).

So far, we have neglected the branch length parameter. You may already have noticed, that the consensus tree has been annotated with branch lengths. To extract ESS and PSRF values for each branch length, run:

$ ./extractBips -f ExaBayes_topologies.myRun.* -n bls

Thus, you obtain the following files:

• ExaBayes_bipartitions.tmp
  contains an identifier for each branch/split (that is explicitly printed)
• ExaBayes_fileNames.tmp
  lists the file names used as input (and assigned an id to them)
• ExaBayes_bipartitionBranchLengths.tmp
  contains all raw branch lengths sampled and lists split id and file ids
• ExaBayes_bipartitionStatistics.tmp
  contains statistics for each branch (specifically the ESS and PSRF)

Notice that for low confidence branches you are likely to encounter poor ESS/PSRF values. This means that you substantially have to increase the number of generations in order to obtain accurate estimates of branch lengths distributions.

• -f alignmentFile

  provides an alignment file. If this file is the binary output produced by parser (see Section 7.1), then no further arguments are required. If you provide a plain (un-processed) Phylip file, then either -m (single partition model) or -q (model file) are mandatory.

• -m DNA | PROT

  specifies the data type used, when a Phylip-formatted alignment has been passed via -f. This way, the alignment is parsed as a single partition with either DNA or amino acid (PROT) data.

• -q modelFile

  specifies a raxml-style partitioning/model scheme for the alignment. For this option, a Phylip-formatted alignment must be passed via -f. See Section 9 for a description of the file format.

• -s seed

  provides a random seed. This number makes the run reproducible. The same seed, data set configuration file will result in the exact same result (apart from limitations given in Section 8.5). If you restart from a checkpoint file, this option will be ignored.

• -n id

  provides a run id used for naming output files

• -r runid

  restarts your run from a previous run id. If your previous ExaBayes-run did not finish (because of a manual abort or walltime restrictions), this option can be used for continuing the run. It is essential, that you pass the same configuration and alignment file. Some adaptions to the configuration file are possible (e.g., larger number of generations to be run, lower topological convergence threshold). Furthermore, all files that carry the previous runid in their name must be located in the current folder.

  Example:

  $ mpirun -np 8 ./exabayes -s $RANDOM -n myId -c myConfig -f myBinaryAlnFile.bin 
  $ [runnig....] -> aborted!
  $ mpirun -np 2  ./exabayes -r myId -n myIdContinued -c myConfig -f myBinaryAlnFile.bin -S
  

• -d

  carries out a dry-run. Only checks your config and alignment file and does not compute anything. Very recommendable, before submitting a large run to a cluster.

• -T n

  Executes Yggdrasil with \(n\) threads. We recommend to use the multi-threaded version of yggdrasil only on systems, where no MPI installation is available.

• -c configFile

  passes a configuration file that specifies how the MCMC will be carried out (see ./examples/configFile-all-options.nex and Section 6 for details)

• -w workDir

  specifies a location for output files

• -R num

  (exabayes-only) specifies the number of runs (i.e., independent chains) to be executed in parallel. Large runs should be carried out as separate runs, see Section 8 for further details.

• -C num

  (exabayes-only) specifies the number of chains (i.e., coupled chains per independent run) to be executed in parallel. Employing this option may be less efficient in terms of runtime and memory than data-level parallelism, see Section 8 for further details.

• -S

  try to save memory using the SEV-technique for gap columns on large gappy alignments Please refer to this paper. On very gappy alignments this option yields considerable runtime improvements.

• -M mode

  specifies the memory versus runtime trade-off. <mode> is a value between 0 (fastest, highest memory consumption) and 3 (slowest, least memory consumption). See Section 8.3 for details.

keyword: params

This section allows to declare and link parameters (e.g., branch lengths) across partitions. You should have declared partitions in the partition file (passed via -q). If you provided a partition file to the parser tool, then the binary output file already contains information about partitions. Partition ids start with 0 and refer to the order provided in the partition file.

Currently the following keywords can be used to specify a parameter linking scheme (keywords are case-insensitive):

Note that, by default all parameters are unlinked for all partitions. Specifically regarding branch lengths, most people only want one global branch length parameter. If a partition id is omitted from the scheme, the default behaviour of ExaBayes is to instantiate a new parameter for this partition (i.e., it is unlinked).

You have the following options for specifying linkage (here demonstrated for the branch length parameter):

• use comma to declare partitions as separate parameters
  example: brlens = (0,1,2,3)
  result: v{0}, v{1}, v{2}, v{3}
  
• use plus to link two partitions into one parameter
  example: brlens = (0 + 1 , 2 , 3)
  result: v{0,1}, v{2}, v{3}
  
• use colon to declare a range of unlinked partitions (each one parameter)
  example: brlens = (0:3)
  result: v{0}, v{1}, v{2}, v{3}
  
• use dash to declare a range partitions linked into one parameter
  example: brlens = (0-3)
  result: v{0,1,2,3}
  

For most use cases, you probably will only want to link all branch lengths. However, in case you work with protein partitions, please consider:

• By default ExaBayes creates one aaModel parameter for each of your amino acid partitions. As state frequencies, ExaBayes uses the empirical frequencies provided by the respective amino acid substitution matrix.
• Instead of using the empirical frequencies, you may want to let ExaBayes integrate over these state frequencies. For doing so, you simply have to declare one of the respective partitions when specifying the stateFreq parameter scheme. If you have two AA partitions, then stateFreq = (0) instructs ExaBayes to integrate over the state frequencies of the first amino acid model parameter.
• As an alternative to proposing fixed-rate AA substitution matrices for AA partitions, you can use ExaBayes to integrate over amino acid GTR matrices (189 free parameters). For doing so, declare (and link) the respective AA partitions in the revMat linking scheme (e.g., revMat = (0+1) for 1 shared GTR matrix across 2 AA partitions).

keyword: priors

The prior block let's you declare your prior belief regarding the values of parameters ExaBayes integrates over. This affects parameters implicitly instantiated by ExaBayes or explicitly defined in a params block (see Section 6.1).

By default priors specifications are applied to all matching parameters. You can overwrite these general priors by specifying parameter-specific priors. For doing so, list all at least one partition that is assigned to your target parameter in curly brackets after the prior keyword. For instance:

brlenPr exponential(10)
brlenPr{0,2,10} uniform(1e-6,10)

applies a uniform prior with \([1e-6,10]\) to all branch length parameters that contain the partitions 0,2 or 10 and applies an exponential prior with \(\lambda = 10\) to all remaining branch length parameters.

All of the following options need to be enclosed within a block featuring the keyword run.

ExaBayes allows you to configure proposals that are used to move your chains through the parameter space. For each proposal, a relative weight governs, how often a specific proposal is drawn. You can customize your proposal mixture by modifying these weights. A proposal provides values for a single parameter only, so a change of the relative weight affects all related proposals. Specifically the topological proposals are described in detail in [3].

Using proposal sets does not change how often a proposal is drawn relative to runtime (so no modifications are necessary).

Moreover, the behaviour of the topological proposals can be customized. The eSPR prunes a subtree, follows down a random path (starting with the original pruning position) and chooses the current branch as re-grafting position with a certain stopping probability (keyword: eSprStopProb). In case of the eTBR, the tree is bisected at a branch and the bisected branch traverses the tree on both ends as described for the eSPR (keyword for the stopping probability is eTbrStopProb).

The parsimony-biased SPR (parsSPR) move prunes a subtree and proposes a regraft position proportionally to the parsimony score of the resulting tree. The parsSPR evaluates the parsimony score for regrafting positions that are no more than \(n\) steps (keyword: parsSPRRadius) apart (i.e., it considers branches within a specified radius for re-insertion). Computing the parsimony score is extremely fast and parallelized in ExaBayes. If you are dealing with large trees, consider increasing the radius. It may not increase mixing, but definitely will reduce the burn-in time and the increase in runtime should not be problematic. The default value depends on the logarithm of the number of taxa (a reasonable assumption, if we do not expect comb-like trees).

Similar to MrBayes, parsimony scores are heated (i.e., exponentiated) using the value of parsimonyWarp. If this value is decreased, the probability that trees with low parsimony score are proposed will get higher.

The posterior-guided SPR move is particularly expensive, since after pruning, it evaluates all reattachment locations in a radius (specified by likeSprMaxRadius) and uses a score based on the posterior to propose a SPR move. How heavily ExaBayes relies on this score can be quantified with likeSprWrap (e.g., choose a likeSprWrap = 0.1 or 0.01 to give topological changes that decrease the posterior by -10 resp. -100 log-units a reasonable chance of being proposed).

moveOptMode allows you to enhance topological proposals by simultaneously proposing branch lengths along with topology (turned off by default). Value 1 means that one branch that is remapped by any NNI/SPR will be proposed. For moveOptMode = 2, all branches that are traversed by a moving subtree will be proposed. For moveOptMode = 3, we additionally propose two branches adjacent and for moveOptMode = 4, adjacent subtrees (e.g., also the root branch of the moving subtree) will be proposed. By default, the NR-based Gamma distribution proposal is used, if useMultiplier is true, then a multiplier is used instead (not recommended).

This utility parses an phylip-formatted alignment and creates a binary representation of this alignment. You either have to indicate the data type of a single partition alignment (via -m) or provide a model file via -q (see Section 9).

Parsing large alignment can take a considerable amount of time that is lost manifold when ExaBayes is executed in parallel.

This utility can be used to summarize (similar to sump in MrBayes or the summary statistics in Tracer) all sampled parameters.

This is straight-forward for continuous parameters (such as substitution rates). If you integrate over fixed protein model matrices (e.g, WAG, LG,…), you are integrating over a discrete parameter. The output in the ExaBayes_parameters* will list the respective matrices. In this case, postProcParam will create an extra column that contains the discrete distribution.

You should check, if all ESS values are greater than 100 and (if available) PSRF values are close to 1 (< 1.1 is considered good convergence).

This utility computes deviations of split frequencies (either maximum or average, abbrev. as ASDSF/MSDSF). If you are integrating over topologies (you usually are), ASDSF/MSDSF are an essential convergence criterion. Usually an ASDSF of \(0.5-1\%\) is considered "excellent convergence" and values between \(1-5\%\) are considered to be acceptable.

You will encounter strongest deviations for branches with low posterior probability.

The stand-alone sdsf computes the same result that is also calculated, if you run multiple independent analysis with convergence criterion. If you run an exceptionally large analysis with multiple independent runs and plan on sampling a very large number of trees, it is recommendable to launch each independent run as a distinct ExaBayes session. You could have a master-script that launches the independent runs (to be run for e.g., 2 h), then checks for convergence and restarts the runs from the respective checkpoints, if not converged yet. If an immense number of processes is involved and your cpu-h budget is tight, this saves you sequential overhead.

This utility computes the credible set of topologies (up to a specified percentile) in one or many tree sets. Use it for post-analyses of your tree samples.

This utility extracts bipartitions (AKA splits or edges) from tree sets and the branch lengths associated with these bipartitions. Note that, this utility also examines trivial bipartitions (these correspond to outer branches in a tree).

extractBips produces the following files:

• ExaBayes_bipartitions.* lists the smaller partition of a bipartition (i.e., all taxa omitted are in the complementary partition) and assigns a unique identifier to the bipartition.
• ExaBayes_fileNames.* lists the file names of the input topology files and assigns a for reference in the remaining two files.
• ExaBayes_bipartitionBranchLengths.* contains all unique branch lengths samples associated with a specific bipartition in a specific file. The file id and bipartition id from the previous two files are used for that.
• ExaBayes_bipartitionStatistics.* contains summary statistics for the branch lengths associated with bipartitions (similar to the output of postProcParam). The ESS value indicates, whether you have sufficiently sampled the branch length associated with a branch and the PSRF value can be used to judge, if the samples from different chains converged against the same distribution. You have sufficiently sampled a parameter, if the ESS is > 100 and a PRSF < 1.1 is an indicator of good convergence.

If a bipartition occurs only in one chain, extractBips will produce -nan-values.

This utility allows to build consensus trees from one or more tree sets. If computing the consensus tree (specifically the extended MR consensus) becomes computationally challenging, you may want to give the parallelized consensus tree algorithm in RAxML a try (use -J MRE).

consense will produces two output files (one in Newick format, one in Nexus format). Nodes are annotated with marginal probability (i.e., confidence), median, mean and \(5\%/95\%\) quantile values.

If you ran analyses with unlinked branch lengths (i.e., you had multiple partitions with a branch length parameter each), then you should create one consensus tree for each branch length parameter. ExaBayes writes one topology file per parameter and per analysis. So you can consense all files that have a "tree-x" (where x is the id of the parameter) in their name.

In most cases, if you have \(x\) processes available (e.g., 4 machines with each 12 cores \(\Rightarrow\) 48 processes):

• if you run \(n\) independent analyses, use -R n.
• if you want to run \(m\) coupled chains, check, if ExaBayes becomes faster, if you increase from -C 1 (default) over -C 2 to -C 4 up to -C m (often the optimum is 2 or 4). For each increment (default: 500 generations), ExaBayes prints the time required for the last increment (use this for benchmarking).
• check the load balance output. For good parallel efficiency, each process should at least have ≈ 100 patterns. Otherwise, less processes may be sufficient, but more do not hurt, if you are not under a budget constraint.

ExaBayes implements three levels of parallelism (in descending order of granularity):

• runs-level parallelism,
• chain-level parallelism,
• data parallelism.

For optimal performance, please consider the following example. Assume, you run \(m\) coupled chains and \(n\) independent runs, while you specify that \(m_p\) coupled chains and \(n_p\) independent runs are run in parallel (via -R and -C). For reasons of load balance, \(m\) should be a multiple of \(m_p\) (analog for \(n\)). Assume each of your computing nodes has \(k\) cores and you want to use \(l\) computing nodes for each parallel working unit (working on one coupled chain in an independent run that is executed in parallel). Thus, you will obtain optimal performance, if you execute ExeBayes with a total number of processes of

If the number of cores \(k\) is divisible by 2, \(l = 2^i\) (where \(i < 0\)) works as well. This way several parallel working units fit on a node.

As mentioned earlier, with ExaBayes you can do Bayesian MCMC on alignments of which the size is only limited by the total memory you have available in your computing center.

In addition to that, ExaBayes implements techniques to reduce the overall memory footprint.

The -M x option allows you to trade runtime for reduced memory consumption. The higher x, the slower but less memory-intensive are the likelihood computations. Remember, that for any -M x ExaBayes will yield the exact same results with the limitations described in section 8.5.

This is particularly relevant, if you use chain-level parallelism, since increased parallelism also increases the memory-overhead as explained in section 8.2.2. Recall that for x=0, you need the \((m+m_p)\) sets of likelihood arrays (where m is the number of coupled chains and \(m_p\) the number of coupled chains executed in parallel). This is, because ExaBayes uses an additional set of likelihood arrays to evaluate the likelihood of a new proposal and saves the previous likelihood arrays for the case of rejection of the proposal.

With -M 1, you can instruct ExaBayes to not save likelihood arrays for arrays for inner nodes that are (recursively) computed from two leave (resp. tip) nodes. These nodes are particularly fast to compute, so you will not loose too much runtime. For a balanced binary tree (best case), the memory consumption of the saved likelihood arrays (adding the \(m_p\) to the equation of memory consumption) is reduced by more than \(50\%\). In the worst case (a comb-/caterpillar-like tree), the memory consumption is merely reduced by 1 array.

When run with -M 2, ExaBayes will only save likelihood arrays for the most expensive kind of nodes. These are nodes that have two inner nodes as their descendants. In terms of memory consumption, a balanced binary tree is the worst case (saving only \(> 50\%\) of the additional likelihood arrays). In the best case (here the comb-like tree), ExaBayes will not save any addition likelihood arrays.

For -M 3, ExaBayes by default does not save any likelihood arrays. The run will be executed substantially slower (but still less than a factor of 2), but specifically if you run a lot of chains in parallel, the factor \((m + m_p)\) in the memory consumption formula is reduced to m.

Aside from that, ExaBayes implements a subtree equality vector-technique, that allows you to save memory for dataset that contain many gaps or undetermined characters (see [8]). The amount of memory you save is proportional to the amount of missing data and the runtime penalty should be negligible (resp., there are instances where this actually increases runtime performance).

In a parallel setting, it is less straight-forward to efficiently execute an analysis, if the alignment is highly partitioned. The issue of load distribution in case of data parallelism is discussed in section 8.2.3. The upshot is that you should, whether all processes of a parallel unit have about the same portion of the overall data. In parallel runs, ExaBayes initially lists, how many characters, partitions, (coupled) chains and runs are assigned to each processor. For efficiency it is important that the option proposalSets is by default set to true (default).

For assigning data (possibly a huge number of partitions) to processors, ExaBayes employs the same algorithm as ExaML [9].

ExaBayes comes with a strong guarantee of reproducibility.

Ideally, the same seed, configuration file and alignment file have to result in the exact same outcome (e.g., topology/parameter samples) regardless whether yggdrasil or exabayes were employed. This should hold for any kind of command line parameter governing the specifics of how calculations are to be performed. Furthermore, repeated continuations from a checkpoint file should not influence the output either.

Any change in the configuration file potentially interferes with perfect reproducibility (e.g., increasing the checkpoint frequency).

When parallelism is involved, this guarantee does not hold necessarily. The reason for this is indeterminism in the calculation of the likelihood, when conducted on multiple processors. Compensating for this problem comes at the cost of runtime performance, thus this has not been implement in ExaBayes.

In other words: running ExaBayes with a different number of processes theoretically may yield different results (however, we have not observed this for any MPI implementations yet).

All of the above does not influence the correctness of the results, however it limits the guarantee that the chain is in the exact same state.