Advanced Topics 3.2
From MrBayes
| Table of contents |
Advanced Topics
Compiling MrBayes
Compiling the MrBayes executable from the source code can be done on several different compilers targeting all the common operating systems: Macintosh, Windows, and Unix. The easiest way to build MrBayes is to use the included configure script and with a make tool for the Makefile this script generates. One can also compile MrBayes with the Metroworks Codewarrior and Microsoft Visual Studio suites.
Some performance gain can be achieved by using different compiler optimization flags. Keep in mind, that not all options are safe to use. A number of optimization methods might give a slowdown instead of a speedup and some options relax numerical assumptions, leading to a program that might produce incorrect results. When you want to file a bug report, make sure that you compiled the program with the default options.
Compiling with GNU Make and configure
The configure script that is distributed with the MrBayes source code examines your build environment for appropriate compiler flags and program settings. In most cases you can run this script (./configure) without any arguments and it will produce a Makefile which can be used to compile the MrBayes executable with the command make. There are a number of options you can give the configure script to alter it's behavior. These are:
--enable-mpi. Use this option if you want to compile the parallel (MPI) version of MrBayes. This option is off by default, meaning that a serial version of MrBayes will be built. See below for more information on how to compile and run the MPI version of MrBayes.
--enable-fastlog. With this option a fast approximation algorithm is used instead of the normal log math functions. Since this approximation algorithm can actually slow down the program on some computer architectures, this option is turned off by default.
--enable-debug. If you want to compile a debug version of MrBayes, you can use this option. This adds the appropriate flag for the GNU gdb debugger.
Before running the configure script, one can set a number of shell-environment variables to influence the compiler and it's flags.
CC. This variable defines which compiler to use. For example, gcc for the GNU compiler or icc for the Intel C compiler. The default setting is the GNU compiler.
CFLAGS. Sets the optimization flags for the compiler. The default is set to -O3, which yields good results for every platform. It is, however, possible to perform some tuning with this variable. We give a few possibilities below for some common processor types, assuming you are using gcc version 3. See the gcc manual for further information on optimization flags.
Intel x86/AMD/AMD64
Some compiler flags for gcc under unix and for gcc/cygwin under windows:
-march=X, with X one of pentium4, athlon-xp or opteron.
If you have one of these processors this will generate instructions specifically tailored for that processor.
-mfpmath=sse attempts to use the SSE extension for numerical calculations.This flag is only effective in combination with the above mentioned -march flag. This flag can provide a big performance gain. However, using this flag in combination with other optimization flags might yield numerically incorrect code. For example, one can set -mfpmath=sse,386
- but this flag leads to incorrect results
-fomit-frame-pointer saves some function overhead. -O2 instead of -O3 turns on a smaller number of optimization flags. However, a number of optimizations turned on by -O3 might give a slowdown instead of producing faster code, especially in combination with the -mfpmath=sse option.
Mac G4 and G5
Some compiler flags for gcc for OS X:
-fast. This flag is specific for the gcc version delivered by Apple.
It turns on a set of compiler flags which tries to optimize for maximum performance. This is the recommended setting if you have a G5 processor and this version of gcc. Code compiled with this flag will not run on a G4 processor.
-O2 or -O3. -mcpu=X, with X one of G4 or G5.
Setting -mcpu or -fast on the Mac results in gcc enabling a number of different flags. Read the gcc manual carefully if you want to experiment with other flags.
Compilation MrBayes version 3.2 will look like,
|
~ $ tar zxf mrbayes-3.2.tar.gz |
Compiling with Code Warrior or Visual Studio
We provide MrBayes project files for both Metrowerks Code Warrior and Microsoft Visual Studio in the source code package. All the relevant flags are set in these files, so you should be able to compile the code without any further modifications.
Compiling and Running the Parallel Version of MrBayes
Metropolis coupling or heating is well suited for parallelization. MrBayes 3 takes advantage of this and uses MPI to distribute heated and cold chains among available processors (Altekar et al., 2004). The optimal number of processors is therefore equal to the total number of chains (the number of runs times the number of chains per run, the default is 2 times 4 equals 8).
As of version 3.2, POOCH is no longer supported by the MrBayes team. If you wish to use POOCH for your Macintosh cluster and need support by the MrBayes team, you can buy a recent version of the POOCH application and donate it to Paul van der Mark.
The MPI version for Unix clusters, including Xserve clusters, has to be compiled before you can run it. To tell the compiler that you want the MPI version, you need to run the configure script with the '--enable-mpi option. The script and the resulting Makefile assume that you have the mpi compiler mpicc in your path and that your system is set up correctly so it can find the relevant libraries.
A typical make session would look as follows:
|
~/mrbayes-3.2 $ ./configure --enable-mpi |
This produces an MPI-enabled version of MrBayes called mb. Make sure that the mpicc compiler is invoked. It is perfectly normal if the build process stops for a few minutes on the mcmc.c file; this is the largest source file and it takes the compiler some time to optimize the code. How you run the resulting executable depends on the MPI implementation on your cluster. A simple approach would use LAM/MPI. First, the LAM virtual machine is set up as usual with lamboot. Then the parallel MrBayes job is started with a line such as
$mpirun –np 4 mb batch.nex > log.txt &
to have MrBayes process the file batch.nex and run all analyses on four processors (-np 4), saving screen output to the file log.txt. If you keep both a serial and a parallel version of MrBayes on your system, make sure you are using the parallel version with your mpirun command.
If your analysis takes a lot of time, we advise you to turn on checkpointing by setting the checkpoint option (mcmcp checkpoint=yes). In case of a crash or power failure, this allows you to restart the MrBayes application from the last checkpoint.
Working with the Source Code
MrBayes 3 is written entirely in ANSI C. If you are interested in investigating or working with the source code, you can download the latest (bleeding edge) version from the MrBayes CVS repository at SourceForge. You can access the CVS repository from the MrBayes home page at SourceForge (http://sourceforge.net/projects/mrbayes) (http://sourceforge.net/projects/mrbayes). SourceForge gives detailed instructions for anonymous access to the CVS repository on their documentation pages.
If you are interested in contributing code with bug fixes, the best way is to send a diff with respect to the most recent file versions in the CVS repository to Paul van der Mark (paulvdm[at]scs.fsu.edu), and we will include your fixes in the main development branch as soon as possible. If you would like to add functionality to MrBayes or improve some of the algorithms, please contact Paul for directions before you start any extensive work on your project to make sure your additions will be compatible with other ongoing development activities. You should also consider whether you want to work with version 3 or version 4 of the program. We are currently shifting our focus to the development of MrBayes 4. Unlike version 3, which is written in C, this version will be written in C++ and our goal is to provide a cleaner, faster, and more extensively documented implementation of Bayesian MCMC phylogenetic analysis. This means, among other things, that the code will be better organized, and all important sections will be documented using Doxygen (http://www.doxygen.org) for easy access to other developers. You are welcome to examine this project as it develops by downloading the source code, doxygen documentation, or programming style directives from the MrBayes CVS repository at SourceForge.
Advanced Options
LSet UseGibbs Option
As described in the Gamma-distributed rates section, MrBayes can accommodate rate heterogeneity across sites using a discrete approximation to the Gamma-distribution. The discrete approximation to the Gamma-distribution is a form of hidden Markov model, in which there are ngammacat rate categories that any site can belong to. There is a "hidden" state that identifies the appropriate rate category for each site. Because we cannot see the rate category for a site, we must treat it as an unknown variable. The LSet UseGibbs option controls how this variable is handled.
When LSet UseGibbs=no is in effect, the likelihoods calculated by MrBayes will be comparable to the likelihood from other software (and versions of MrBayes earlier than 3.2). In this mode, a weighted sum over all rate categories is performed to calculate each site's likelihood. The weights are the prior probabilities that a site would belong to that category. Because the gamma is discretized by breaking it up into ngammacat equal size chunks, the probability that any site i, would be in rate category c is simply 1/ngammacat (and this quantity is independent of the site of the rate). The likelihood for site i is a sum over all categories c of the quantity
(1/ngammacat)L(i | c)
where L(i | c) is the likelihood of site i conditional upon it being in category c. Typically (at least for datasets with a large number of taxa), one of the rate categories contributes the vast majority of the likelihood to this sum, because its conditional likelihood is much larger than the conditional likelihoods from the other rate categories. Thus the likelihoods under this model are close to 1/ngammacat times the conditional likelihood of the best fitting rate category for each site.
The fact that the site likelihoods are dominated by one term implies that it is wasteful to calculate the conditional likelihoods over all ngammacat rate categories --most of the calculations will be lost in the rounding error of the final summation over rate categories.
The LSet UseGibbs=yes option addresses this concern. Instead of calculating a site's likelihood by summing over all of the hidden states (the rate categories), an indicator variable that corresponds to the rate category for a site is sampled during the MCMC. The "likelihood" reported for each site is actually the conditional likelihood - the likelihood of the site conditional upon it belonging to rate category. Each site has its own value for the hidden state. The fact that we are uncertain of each site's "true" rate category is accommodated by sampling with MCMC instead of summation.
Conveniently, when you calculate the likelihood by summing over all rate categories, then you are also in position to calculate the posterior probability that each site belongs to a particular category. To update the hidden states for each site, you have to:
- Calculate the conditional likelihood for each rate category.
- Calculate the full site-likelihood by summing over all rate categories (weighted by their prior probabilities).
- Calculate the posterior probability of each of the hidden state assignments, by dividing the conditional likelihood for each rate category by the site likelihood (so that the probabilities sum to one).
- Select a new hidden state for each site by choosing a randomly -- giving each rate category a probability of being chosen that is equal to the posterior probability that the site belongs to that category
This is a Gibbs Sampling (http://en.wikipedia.org/wiki/Gibbs_sampling) update of the hidden state. The LSet Gibbsfreq option controls how frequently the hidden states are updated. Thus the likelihood for all ngammacat categories only has to be calculated in the iteration in which the hidden state is being updated In all other iterations the likelihood is calculated from only the rate category to which that site is currently assigned.
Note that in the LSet UseGibbs=no mode, the priors associated with each rate category (the 1/ngammacat terms) are used in calculating the site likelihood. In the LSet UseGibbs=yes mode these priors are dealt with when rate categories are sampled. Thus, when using Gibbs Sampling the likelihood reported is not decreased by weighting it by a prior probability (thus the likelihoods that are reported will be higher).
