Before you can use BAli-Phy on Windows, you need to install a Unix command-line environment. We recommand installing Cygwin. You may then access the Unix command line environment by running the Cygwin shell (not the normal windows command line). The Cygwin shell mounts the C: drive on /cygdrive/c/, so you can access the directory C:/Users/ as /cygdrive/c/Users/ from within the Cygwin shell, for example.

While running the Cygwin installer setup-x86_64.exe, you will be given an opportunity to select additional packages.

You can re-run the installer to add packages that you did not add during the initial install.

You can optionally use MSYS2 instead of Cygwin. Both MSYS2 and Cygwin can be installed at the same time. After installing MSYS2, You may access the Unix command line environment by running the MSYS2 shell (not the normal windows command line). The MSYS2 shell mounts the C: drive on /c/, so you can access the directory C:/Users/ as /c/Users/ from within the MSYS2 shell, for example.

After installing MSYS2 you will need to install a few packages before you proceed. Run the MSYS2 shell, and enter the command:

% pacman -S perl tar

After installing a Unix command line, use it to download and extract the bali-phy executables from the website:

% mkdir -p ~/Applications
% cd ~/Applications
% wget http://www.bali-phy.org/files/bali-phy-3.4-win64.tar.gz
% tar -zxf bali-phy-3.4-win64.tar.gz

Second, check that the bali-phy executable runs:

% ~/Applications/bali-phy-3.4/bin/bali-phy --version

You still need to add it to your PATH as described in Section 2.6, “Add BAli-Phy to your PATH”.

First install the XCode (version 6 or higher) command line tools:

% xcode-select --install

Then install homebrew and use homebrew to compile and install bali-phy:

% brew tap brewsci/bio
% brew install bali-phy

Check that the executable runs:

% bali-phy --version

If you install with homebrew, you don't need to do anything extra to put bali-phy in your PATH.

Open a windows in the Terminal app to access the UNIX command line. Then download and extract the executables:

% mkdir -p ~/Applications
% cd ~/Applications
% curl -O http://www.bali-phy.org/files/bali-phy-3.4-mac64.tar.gz
% tar -zxf bali-phy-3.4-mac64.tar.gz

Check that the executable runs:

% ~/Applications/bali-phy-3.4/bin/bali-phy --version

You still need to add it to your PATH as described in Section 2.6, “Add BAli-Phy to your PATH”.

You can install gnuplot via homebrew:

%  brew install gnuplot

You can install R via homebrew:

%  brew tap caskroom/cask
%  brew cask install xquartz
%  brew install r

However, note that this might conflict with R installed from other places, such as MRAN.

You can install Figtree with homebrew:

% brew tap caskroom/cask
% brew cask install figtree

However, Seaview and Tracer don't have homebrew packages at the moment.

% sudo apt-get install bali-phy

% bali-phy --version

First install wget. If you have Debian or Ubuntu Linux, type:

% sudo apt-get install wget

Then download and extract the executables:

% mkdir -p ~/Applications
% cd ~/Applications
% wget http://www.bali-phy.org/files/bali-phy-3.4-linux64.tar.gz
% tar -zxf bali-phy-3.4-linux64.tar.gz

Second, check that the executable runs:

% ~/Applications/bali-phy-3.4/bin/bali-phy --version

You still need to add it to your PATH as described in Section 2.6, “Add BAli-Phy to your PATH”.

If you have Debian or Ubuntu Linux, you can install other recommended programs by typing:

% sudo apt-get install gnuplot
% sudo apt-get install r-base

% sudo apt-get install seaview
% sudo apt-get install figtree

However, there isn't a Debian or Ubuntu package for Tracer at the moment.

First check if the executable is in your PATH.

% bali-phy --version

If this shows version info, then bali-phy is already in your PATH and you can skip this section. This should be true if you installed bali-phy using a package manager such as homebrew or apt, or if you've already added it to your PATH.

If bali-phy is not in your path, then you should see:

% bali-phy --version
bali-phy: command not found.

If bali-phy is not in your PATH, then continue with this section.

Add bali-phy to your PATH, so that the shell knows where to find it. This command only affects the terminal in which it is typed, and will not affect new terminals:

% export PATH=~/Applications/bali-phy-3.4/bin:$PATH

To set the PATH automatically for new terminals, type:

% test -r ~/.bash_profile && echo 'export PATH=~/Applications/bali-phy-3.4/bin:$PATH' >> ~/.bash_profile
% echo 'export PATH=~/Applications/bali-phy-3.4/bin:$PATH' >> ~/.profile

This will affect new terminals only after you log out and log back in though.

Now check that the executable runs:

% bali-phy --version

If it does, then your PATH is set up correctly, and you can probably skip the rest of this section.

If you installed BAli-Phy to the directory ~/Applications, then you can run bali-phy by typing ~/Applications/bali-phy-3.4/bin/bali-phy. However, it would be much nicer to simply type bali-phy and let the computer find the executable for you. This can be achieved by putting the directory that contains the BAli-Phy executables into your "path". The "path" is a colon-separated list of directories that is searched to find program names that you type. It is stored in an environment variable called PATH.

Setting your PATH is also a pre-requisite for running the bp-analyze script to summarize your MCMC runs.

You can examine the current value of this environment variable by typing:

% echo $PATH

We will assume that you extracted the bali-phy archive in ~/Applications and so you want to add $HOME/Applications/bali-phy-3.4/bin to your PATH. (If you installed to another directory, replace $HOME/Applications/bali-phy-3.4/ with that directory.)

The commands for doing this depend on what "shell" you are using. Type echo $SHELL to find out. If your shell is sh or bash then the command looks like this:

% PATH=$HOME/Applications/bali-phy-3.4/bin:$PATH

If your shell is csh or tcsh, then the command looks like this:

% setenv PATH $HOME/Applications/bali-phy-3.4/bin:$PATH

Note that these commands will only affect the window you are typing in, and will vanish when you reboot.

To make this change survives when you logout or reboot, open your shell configuration file in a text editor, and add the command on a line by itself. This will ensure that it is run every time you log in.

To find the right configuration file, look in your $HOME directory for .profile (for the Bourne shell sh), .bash_profile (for BASH), or .login (for tcsh). You may have to create the file if it is not present. On Cygwin, you should put the change in the file .bashrc.

If you do not know which directory is your home directory, you can find its full name by typing:

% echo $HOME

BAli-Phy is quite CPU intensive, and so we recommend using 150 or fewer taxa in order to limit the time required to accumulate enough MCMC samples.

When designing an MCMC analysis, I recommend performing an initial analysis with a much smaller number of sequences. This smaller analysis will run much faster, and allow discovering mistakes much more quickly. Then, after you are sure that you are running the program correctly and have chosen the best model, you can ramp up the number of sequences towards your desired number.

Aligning just a pair of sequences takes $O(L^2)$ time and memory, where $L$ represents the sequence length. Therefore sequences longer than (say) 1000 letters become increasingly impractical. However, you might try to see how long you can make your sequences before you run out of memory, or the program becomes too slow.

For multi-gene analyses, two separate data partitions (i.e. genes) of 500 letters will be twice as fast to align as one data partition of 1000 letters. So, it may be possible to analyze several genes as long as each gene individually is not too long.

Also, note that you can sometimes speed up the analysis of protein sequences by coding them as amino acids or codons, rather than nucleotides. This is because it decreases the sequence length.

This section explains the meaning of the various field names in the file C1.log.

The "prior" field includes the probability of the alignment, since the alignment is not observed.

The likelihood is the probabilistic analogue to summed mismatch penalties.

The prior_A is the probabilistic analogue to summed gap penalties.

The prefixes "S$n$/" and "I$n$/" will be dropped if not necessary to disambiguate parameters with the same name in different sub-models.

To compute the majority consensus tree, do the following. (The program FigTree allows you to view the resulting tree file graphically.)

% trees-consensus dir-1/C1.trees dir-2/C1.trees > c50.PP.tree

By default, the first 10% of tree samples are skipped as burn-in (--skip=10% or -s 10%) and every generation is analyzed (--subsample=1 or -x 1). To discard the first 1000 tree samples and analyze every 10th sample:

% trees-consensus -s 1000 -x 10 dir-1/C1.trees dir-2/C1.trees > c50.PP.tree

By default, splits are included in the consensus tree if they have a PP greater than 0.5. You can specify a more stringent level (e.g. 0.66) by adding the option --consensus-PP=0.66 as follows:

% trees-consensus -s20% -x10 --consensus-PP=0.66 dir-1/C1.trees dir-2/C1.trees > c66.PP.tree

You may also make the program write directly to the output file (e.g. c66.PP.tree) by using the more general form --consensus-PP=0.66:c66.PP.tree. Leaving off the ":c66.PP.tree" part (as we did above) or specifying ":-" sends the output to the standard output (e.g. the terminal, if not redirected).

% trees-consensus -s20% -x10 dir-1/C1.trees dir-2/C1.trees --consensus-PP=0.66:c66.PP.tree

You can supply multiple levels and filenames separated by commas. This is faster than running the program multiple times with different consensus levels.

% trees-consensus -s20% -x10 dir-1/C1.trees dir-2/C1.trees --consensus-PP=0.5:c50.PP.tree,0.66:c66.PP.tree

Finally, you may use the option --consensus= instead of the option --consensus-PP= if you do not wish the resulting tree to contain embedded posterior probabilities on branches, as well as branch lengths.

% trees-consensus -s20% -x10 dir-1/C1.trees dir-2/C1.trees --consensus=0.5:c50.PP.tree,0.66:c66.PP.tree

Both the --consensus= and --consensus-PP= options may be given simultaneously.

See trees-consensus --help for a complete list of options.

The greedy consensus tree may be used instead of a majority-consensus tree when a fully resolved (e.g. bifurcating) tree is required. When the topology has many tips and each topology may be sampled only once, the greedy consensus should be higher quality than the estimate of the MAP topology. To obtained a fully resolved tree, the greedy consensus strategy starts with the majority consensus and then adds the highest-supported split that does not conflict.

To compute the greedy consensus tree do:

% trees-consensus --skip=burnin dir-1/C1.trees dir-2/C1.trees --greedy-consensus=greedy.tree

To compute the maximum a posteriori tree do:

% trees-consensus --skip=burnin dir-1/C1.trees dir-2/C1.trees --map-tree=MAP.tree

When the tree has many tips, each topology may be sampled only once, leading to low quality estimates of the MAP topology. As a result, when you need a bifurcating tree you should probably use the greedy consensus instead.

% trees-bootstrap dir-1/C1.trees dir-2/C1.trees

This command computes the effective sample size for the posterior probability of each split. It also computes the Average Standard Deviation of Split Frequencies (ASDSF) between two or more independent runs.

See Section 10, “Convergence and Mixing: Is it done yet?” for more information.

This command gives a median and confidence interval, ESS, and a stabilization time:

% statreport dir-1/C1.log dir-2/C1.log > Report 

When multiple runs are analyzed, this command gives PSRF and joint ESS values. The program Tracer allows you to view the same summaries graphically.

See Section 10, “Convergence and Mixing: Is it done yet?” for more information.

To construct an alignment estimate via posterior decoding, select any tree file tree that corresponds to your alignment. It does not need to be fully resolved.

% cut-range dir-1/C1.Pp.fastas dir-2/C1.Pp.fastas --skip=burn-in | alignment-chop-internal --tree tree | alignment-max > Pp-max.fasta

You can optionally replace --tree tree with -N n_sequences, where n_sequences is the number of non-ancestral sequences in your alignment.

You can use the program SeaView to view the alignment graphically.

To annotate a specific alignment alignment.fasta, choose a fully resolved tree estimate tree:

% cut-range dir-1/C1.Pp.fastas dir-2/C1.Pp.fastas --skip=burn-in | alignment-chop-internal --tree tree  | alignment-gild alignment.fasta tree  > alignment-AU.prob
% alignment-draw alignment.fasta --AU alignment-AU.prob > alignment-AU.html

The majority consensus tree is usually not fully resolved, so we recommend using the greedy consensus instead.

The Results/ directory will contain the following useful files:

The following files will be generated to summarize alignment uncertainty, unless the analysis uses a fixed alignment.

The following files describe convergence and mixing:

The SRQ plots can be viewed by typing "plot 'file' with lines" in gnuplot.

This file reports the quality of estimates of support for each partition in terms of the posterior probability (PP) and log-10 odds (LOD). It also reports the auto-correlation time (ACT), the effective sample size (Ne), the number of samples that support (1) or do not support (0) the partition, and the number of regenerations. Only partitions with PP > 0.1 are shown by default.

All the DNA models are special cases of the GTR model.

Frequencies are estimated by default. Frequencies can be fixed by setting the pi parameter to a constant value, if the model allows unequal frequencies.

Constant frequencies are specified as a list of pairs that associates each letter with its frequency:

gtr[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]

Frequencies can also be specified using functions:

gtr[pi=Frequencies.uniform]

Frequencies are estimated by default. Frequencies can be fixed by setting the pi parameter to a constant value, if the model allows unequal frequencies.

Constant frequencies are specified as a list of pairs that associates each letter with its frequency:

wag+f[pi=List[Pair["A",0.047],Pair["R",0.19],...]]

Frequencies can also be specified using functions:

wag+f[pi=Frequencies.uniform]]

The +fe model is shorthand for +f[pi=Frequencies.uniform]:

wag+fe

As of version 3.4, BAli-Phy does not yet allow specifying which nucleotides are paired either with a string like ((.)) or with a "pairs" file. Instead you must manually extract the paired nucleotides and put them in their own partition (for stems), and then manually extract each loop and put it in its own partition.

The stems should be arranged so that paired nucleotides are adjacent. For example, suppose the sequence AGGCT was paired according to ((.)). Then the input file for the stems should contain a sequence of doublets that looks like ATGC, where AT is the first pair, and GC is the second pair. Later versions of the software should allow extracting stems and loops from nucleotide sequences using parenthesis notation or a "pairs" file.

Frequencies are estimated by default. Frequencies can be fixed by setting the pi parameter to a constant value, if the model allows unequal frequencies.

Constant frequencies are specified as a list of pairs that associates each letter with its frequency.

hky85[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]+x2

hky85_sym+x2_sym+f[pi=List[Pair["AA",0.01],Pair["AC",0.01],Pair["AG",0.01],Pair["AU",0.22],Pair["CA",0.01],Pair["CC",0.01],Pair["CG",0.22],Pair["CU",0.01],Pair["GA",0.01],Pair["GC",0.22],Pair["GG",0.01],Pair["GU",0.01],Pair["UA",0.22],Pair["UC",0.01],Pair["UG",0.01],Pair["UU",0.01]]]

Frequencies can also be specified using functions:

BAli-Phy interprets branch lengths for doublet models as 1/2 the number of substitutions per doublet. Thus, they should be comparable to branch lengths under DNA/RNA nucleotide models.

Frequencies are estimated by default. Frequencies can be fixed by setting the pi parameter to a constant value, if the model allows unequal frequencies.

Constant frequencies are specified as a list of pairs that associates each letter with its frequency.

hky85[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]+x3

Frequencies can also be specified using functions:

hky85_sym+x3_sym+f[pi=f1x4]            // nucleotide frequencies are estimated

The +fe model is shorthand for +f[pi=Frequencies.uniform]:

hky85_sym+x3_sym+fe

BAli-Phy interprets branch lengths for codon models as 1/3 of the number of substitutions per codon. Thus, they should be comparable to branch lengths under DNA/RNA models.

The x3, x3_sym, x3x3, dNdS, and mut_sel models can be used to build up codon models piecewise:

Frequencies are estimated by default. Frequencies can be fixed by setting the pi parameter to a constant value, if the model allows unequal frequencies.

Constant frequencies are specified as a list of pairs that associates each letter with its frequency.

gy94[pi=List[Pair["AAA",0.01],Pair["C",0.02],...]]
mg94[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]

Frequencies can also be specified using functions:

gy94[pi=f1x4]              // nucleotide frequencies are estimated

When using a codon-based substitution model like gy94, you may select the genetic code by specifying -A Codons[,genetic-code]. Available genetic codes are standard, mt-vert, mt-invert, mt-yeast, mt-protozoan.

If the genetic code is not specified, then the standard code is used:

% bali-phy sequence-file -S gy94 -A Codons
% bali-phy sequence-file -S gy94 -A Codons[RNA]

These examples specify the vertebrate mitochondrial code:

% bali-phy sequence-file -S gy94 -A Codons[DNA,mt-vert]
% bali-phy sequence-file -S gy94 -A Codons[,mt-vert]

In order to use the branch-site substitution model, the user needs to specify an unrooted tree topology and disable topology changes in order to keep the topology fixed:

% bali-phy alignment.fasta -S branch_site -T tree.tree --disable=topology

The tree file should be in Newick format, with foreground branches labelled using NHX attributes. The NHX attribute must be applied to the branch, not the node, so it must occur after a colon.

The posterior probability of positive selection is the posterior mean of the posSelection parameter. This may be computed using the statreport program with the --mean option. In case this probability is extremely close to 1 or 0, you may wish to add the option --Rao-Blackwellize branch_site:posSelection. This will report the log-probability of positive selection each iteration. The user may exponentiate the reported values and then average them (using R, for example) in order to compute a more accurate estimate of the posterior probability of positive selection.

Priors on model parameters are specified by giving a random value. Random values can be obtained from distributions using the function sample. For example, this places a log-normal prior on the parameter kappa of the hky85 model:

hky85[kappa=sample[log_normal[1,1]]]

You can write ~Dist as a shorthand for sample[Dist]:

hky85[kappa=~log_normal[1,1]]

The =~ can be further shortened to just ~:

hky85[kappa~log_normal[1,1]]

It also is possible to use random values as inputs to other functions. For example:

add[1,~exponential[10]]

In such cases the parameter value should be specified with =, as in the following example:

rs07[mean_length=add[1,~exponential[10]]]]

Random values and distributions have different types. For example, the following is of type Distribution[Double]:

uniform[0,1]

In contrast, the following are both of type Double:

sample[uniform[0,1]]
~uniform[0,1]

This is important when passing distributions as arguments to other distributions and functions. For example, the distribution iid is used to generate a specific number of samples from another distribution. Thus, it needs to receive a distribution as an argument:

~iid[4,normal[0,1]]      // OK    : 4 samples from the normal[0,1] distribution
~iid[4,~normal[0,1]]     // not OK: 4 samples from ... a random number?

(See Section 8.5, “Argument and result types”.)

To calculate the ASDSF and MSDSF run:

% trees-bootstrap dir-1/C1.trees dir-2/C1.trees ... dir-n/C1.trees > partitions.bs

For each split, the SDSF value is just the standard deviation across runs of the Posterior Probabilities for that split. By averaging the resulting SDSF values across splits, we may obtain the ASDSF value (Huelsenbeck and Ronquist 2001). This is commonly considered acceptable if it is < 0.01.

However, it is also useful to consider the maximum of the SDSF values (MSDSF). This represents the range of variation in PP across the runs for the split with the most variation.

To generate the split-frequency comparison plot, you must have R installed. Locate the script compare-runs.R. Then run:

% trees-bootstrap dir-1/C1.trees dir-2/C1.trees ... dir-n/C1.trees --LOD-table=LOD-table > partitions.bs 
% R --slave --vanilla --args LOD-table compare-SF.pdf < compare-runs.R

Following Beiko et al (2006), this displays the variation in estimates of split frequencies across runs. Splits are arranged on the x-axis in increasing order of Posterior Probability (PP), which is obtained by averaging over runs. We then plot a vertical bar from the minimum PP to the maximum PP.

To calculate the split ESS values, run:

% statreport dir-1/C1.log dir-2/C1.log ... dir-n/C1.log > Report 

We calculate effective sample sizes based on integrated autocorrelation times. This method has the nice property that simply duplicating every sample does not increase the ESS.

The program Tracer also computes ESS values.

To calculate the split ESS values, run:

% trees-bootstrap dir-1/C1.trees dir-2/C1.trees ... dir-n/C1.trees > partitions.bs

To compute the ESS for a split, we consider the presence or absence of a split in each iteration as a series of binary values. We compute the integrated autocorrelation time for this binary sequence, which leads to an ESS. This approach is similar to dividing the iterations into blocks and computing the ESS on the PP estimates in the blocks. It is also similar to estimating the variance reduction under a block bootstrap.

To obtain estimates of the stabilization time for each numerical parameter, you may run:

% statreport C1.log > Report 

Each series of values is counted as having stabilized after the series crosses its upper and then lower 95% confidence bounds twice (if the initial value is below the median) or crosses its lower and then upper confidence bounds twice (if the initial value is above the median). The confidence bounds are those based on its equilibrium distribution as calculated from the last third of the values in the sequence.

In addition to examining convergence diagnostics for continuous parameters, it is important to examine convergence diagnostics for the topology as well (Beiko et al 2006). In theory, we recommend the web tool Are We There Yet (AWTY) (Wilgenbush et al, 2004). However, AWTY gives incorrect results if you upload plain NEWICK tree samples -- which is what BAli-Phy outputs. Therefore, if you wish to use AWTY, you must convert the tree samples files to NEXUS before you upload them to AWTY in order to get correct results.

It is also be possible to assess stabilization of tree topologies using tools distributed with bali-phy by using commands like the following. Here, sub-sampling and burnin does not apply to the equilibrium tree files. Also, note that you need to manually construct the equilibrium samples, which we recommend to contain at least 500 trees; you might do this by sub-sampling using the BAli-Phy tool sub-sample.

Note that the running time is the product of the number of trees in the two files. Therefore, comparing two complete tree samples without sub-sampling will take too long.

On Debian and Ubuntu, you can type:

% sudo apt-get install g++ git libcairo2-dev pandoc

% sudo apt-get install meson
% dpkg -s meson | grep Version
Version: 0.45.1-2

% sudo apt-get install python3 python3-pip ninja
% python3 -m venv meson
% source meson/bin/activate
% pip3 install meson

On Mac OS X, the simplest way to get a compiler is to install XCode (version 6 or newer) command line tools, which come with clang.

% xcode-select --install

To get the other tools, first install homebrew, and then type:

% brew install git meson cairo pandoc

The MSYS2 project provides an MINGW64 compiler that can create native windows executables. MSYS2 itself is actually non-native (it is derived from cygwin), and therefore the MSYS2 shell refers to drives as /c/ instead of C:/.

% pacman --needed --noconfirm -Sy pacman-mirrors
% pacman -Sy
% pacman -S mingw-w64-x86_64-ninja
% pacman -S mingw-w64-x86_64-toolchain
% pacman -S mingw-w64-python3-pip
% PATH=/c/msys64/mingw64/bin:$PATH # Put the mingw64 executables into your path
% pip3 install meson

Keep in mind that MSYS2 keeps its (non-native) executables in C:/msys64/usr/bin, while it keeps the (native) MINGW executables in C:/msys64/mingw64/bin. If you want to use the native MINGW executables, you need to make sure that /c/msys64/mingw64/bin/ is in your PATH. If you forget to put the MINGW executables in the path, some of the installed MINGW programs (such as pip3 above) will show up as missing when you try to run them.