BAli-Phy User's Guide v4.0-beta13

We typically run BAli-Phy on workstations with at least 8Gb of RAM and 2 cores. More cores will allow you to run more MCMC chains at once, and more RAM will allow you to run larger data sets. However, it is often easier and faster to run BAli-Phy on a (Linux) computing cluster, if you have one available.

If you have previously installed bali-phy, you do not have to remove the old version before installing the new version. Simply follow the installation instructions for the new version. If you are manually adding the new version of bali-phy to your PATH, just make sure that the new version comes before the old version in the PATH, or remove the old version from the PATH.

In order to remove an older version, simply delete the directory bali-phy-oldversion. This will completely uninstall the old version from the system. BAli-Phy does not create hidden files that will remain after you remove its directory.

First check that you have a 64-bit version of the Windows operation system installed. The executables for download will only run on a 64-bit installation of Windows.

% cygpath -w ~/Applications
C:\cygwin64\home\username\Applications\

In order to determine that the software has been correctly installed, and the PATH has been correctly set, run the following commands:

% cp ~/Applications/bali-phy-4.0-beta13/share/doc/bali-phy/examples/sequences/5S-rRNA/25.fasta .
% bali-phy 25.fasta --iter=150
% bali-phy 25.fasta --iter=150
% bp-analyze 25-1 25-2

Furthermore, the directories 25-1 and 25-2 should contain a file called C1.log. You should be able to load these files in Tracer, although the chain will not really have converged yet.

The simplest way to run BAli-Phy is to type all the arguments on the command line:

% bali-phy sequences.fasta

You can run a traditional fixed-alignment Bayesian tree inference by adding -I none:

% bali-phy sequences.fasta -I none

You can also specify a character set for analysis:

% bali-phy sequences.fasta:1-30,90-100

Sequence files can be in FastA or PHYLIP format. FASTA format prefixes sequence names with ">":

>human       this is a comment and is not part of the sequence name
CTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAA
GTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTGGTCTACCCTTGGACCCAGAGGTTCTTT
>tarsier     this is also a comment
CTGACTGCTGAAGAGAAGGCCGCCGTCACTGCCCTGTGGGGCAAGGTAGACGTGGAAGAT
GTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTCGTCTACCCATGGACCCAGAGGTTCTTT
>bushbaby
CTGACTCCTGATGAGAAGAATGCCGTTTGTGCCCTGTGGGGCAAGGTGAATGTGGAAGAA
GTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTTGTCTACCCATGGACCCAGAGGTTCTTT
>hare
CTGTCCGGTGAGGAGAAGTCTGCGGTCACTGCCCTGTGGGGCAAGGTGAATGTGGAAGAA
GTTGGTGGTGAGACCCTGGGCAGGCTGCTGGTTGTCTACCCATGGACCCAGAGGTTCTTC

If the sequence-name line contains a space, BAli-Phy treats everything after the space as a comment.

The sequences in the file do not need to be aligned unless you fix the alignment with -I none.

Sensible defaults are supplied for command line options that are not specified. For example, if sequences.fasta contains DNA sequences, then

% bali-phy sequences.fasta

is equivalent to

% bali-phy sequences.fasta -A DNA -S tn93 -I rs07

Default values that are used will always be displayed on the screen and in the output files so that you do not have to guess. You can specify a more complex substitution model using the -S option. You will generally need to write the substitution model inside single quotes unless it is just a single word.

% bali-phy sequences.fasta -S 'lg08 +> Rates.gamma +> inv'

Every short option like -S has an equivalent long option like --smodel. To see the most frequently-used command-line options, you can run

% bali-phy help

In addition to using the command line, you may also specify options in a file. Option files also use the long form of command line options. Each option is given on its own line using the syntax "option = value" instead of the syntax "--option value". The value can be blank if the option does not take an argument. The align option indicates sequence files. Lines that begin with # are comments, and blank lines are ignored.

For example, consider the following option file:

# sequence data for 3 genes/partitions
align = ITS1.fasta
align = 5.8S.fasta
align = ITS2.fasta

# linked substitution model for 1st and 3rd partition
smodel = 1,3:tn93 +> Rates.free(n=3)

# substitution model for 2nd partition
smodel = 2:tn93

# indel model for second partition
imodel = 2:none

# linked scale for 1st and 3rd partition
scale = 1,3:

# choose a name for output directories
name = ITS-analysis1

Options files are specified with the -c option_file option:

% bali-phy -c analysis1.txt                           # run the analysis
% bali-phy -c analysis1.txt --name ITS-analysis1b     # override the name

Options given on the command line will override values given in the option file.

Running bali-phy on a computing cluster is not necessary, but can speed up the analysis dramatically. This is because a cluster allows you to run several independent MCMC chains simultaneously and pool the resulting samples. You can run multiple chains simultaneously simply by starting several different instances of bali-phy. Each instance of bali-phy runs only one chain and does not require using MPI or special command-line options.

This approach to parallel computation is sometimes more efficient than MCMCMC-based parallelism involving heated chains. It is equivalent to running MCMCMC with no temperature difference between chains, with the exception that it allows results from all chains to be used, instead of just results from the single "cold" chain. Thus, if you run 10 independent chains in parallel, then you may gather samples 10 times faster that a single chain.

BAli-Phy can read in sequences and alignments in both FastA and PHYLIP formats. Filenames for FastA files should end in .fasta, .mpfa, .fna, .fas, .fsa, or .fa. Filenames for PHYLIP files should end in .phy. If one of these extensions is not used, then BAli-Phy will attempt to guess which format is being used.

Large data sets run more slowly than small data sets. We recommend a conservative starting point with few taxa and short sequence lengths. You can then increase the size of your data set until a balance between speed and size is reached. The tool alignment-thin described in Section 12, “Alignment utilities: brief overview” can be used to construct a smaller data set.

The number of MCMC samples that you need depends on whether you are primarily interested in obtaining a point estimate or in obtaining detailed measures of confidence and uncertainty. For detailed measures of confidence and uncertainty you should obtain a minimum of 10,000 samples after the Markov chain converges. For an estimate, you don't need very many samples after convergence. (But you may need many samples to be sure that you've converged!)

See also Section 3.5, “Running on computing clusters”.

BAli-Phy creates a new directory to store its output files each time it is run. By default, the directory name is the name of the sequence file, with a number added on the end to make it unique. BAli-Phy first checks if there is already a directory called file-1/, and then moves on to file-2/, etc. until it finds an unused directory name.

You can specify a different name to use instead of the sequence-file name by using the --name option.

BAli-Phy writes the following output files inside the directory that it creates:

This section is primarily about extracting estimates from output files. See Section 11, “Convergence and Mixing: Is it done yet?” for methods of determining effective sample sizes, and for checking mixing and convergence.

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

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

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

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

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

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

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

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

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

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

% 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

% 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

Instead of manually running each of the steps to analyze the output files, you may instead run the PERL script bp-analyze to execute these commands. The script will create an HTML page Results/index.html that summarizes the posterior distribution.

You may run bp-analyze inside the output directory, like this:

% bp-analyze --skip=iterations

You may also run it with one or more output directories as arguments, like this:

% bp-analyze --skip=iterations directory-1/ directory-2/

In this case, output from multiple runs will be used to assess convergence and mixing, as well as to increase the precision of the estimates.

All the commands that are executed by bp-analyze will be logged to Results/commands.log. You can also see these commands as they are executed by supplying the --verbose option:

% bp-analyze --skip=iterations --verbose

The default substitution model for DNA and RNA is tn93.

gtr(pi={"A":0.1, "C":0.2, "T":0.3, "G":0.4})

gtr(pi=Frequencies.uniform)

The default substitution model for proteins is lg08.

wag +> f({"A":0.047, "R":0.19,...})

wag +> f(pi=Frequencies.uniform)

wag +> fe

The doublets alphabet consists of 16 RNA dinucleotides. It is used to model RNA stems, where two nucleotides matched in the RNA secondary structure are highly correlated.

The default substitution model for doublets is tn93_sym +> x2_sym +> f.

hky85(pi={"A":0.1, "C":0.2, "T":0.3, "G":0.4}) +> x2

hky85_sym +> x2_sym +> f({"AA":0.01, "AC":0.01, "AG":0.01, "AU":0.22, "CA":0.01, "CC":0.01, "CG":0.22, "CU":0.01, "GA":0.01, "GC":0.22, "GG":0.01, "GU":0.01, "UA":0.22, "UC":0.01, "UG":0.01, "UU":0.01})

The triplets alphabet is similar to the codons alphabet, except that stop codons are included. Unlike the codons alphabet, the triplets alphabet has no knowledge of the genetic code.

The default substitution model for triplets is tn93 +> x3.

hky85(pi={"A":0.1, "C":0.2, "T":0.3, "G":0.4}) +> x3

hky85_sym +> x3_sym +> f(pi=f1x4)            // nucleotide frequencies are estimated

hky85_sym +> x3_sym +> fe

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

The default substitution model for codons is gy94.

gy94(pi={"AAA":0.01, "C":0.02,...})
mg94(pi={"A":0.1, "C":0.2, "T":0.3, "G":0.4})

gy94(pi=f1x4)              // nucleotide frequencies are estimated

% bali-phy sequence-file -S gy94 -A Codons
% bali-phy sequence-file -S gy94 -A Codons(RNA)

% bali-phy sequence-file -S gy94 -A Codons(DNA,mt-vert)
% bali-phy sequence-file -S gy94 -A Codons(,mt-vert)

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

Complex substitution models in BAli-Phy are constructed as mixtures of reversible CTMC models that run at different rates (e.g. ${\Gamma }_4+\mathrm{INV}$) or have different parameters (e.g. an M2a codon model).

These models attempt to model the fact that evolutionary rates may change over time within a single column. These models are sometimes called "covarion" models, based on the idea that changes in rate might be caused by changes in an unspecified covarying site.

These models are "Markov modulated" models that create multiple different states for each letter by augmenting each letter with some unobserved hidden state. They attempt to model the fact that substitution processes might not be Markov on the letters, but might become more Markov given the hidden state.

Models, probability distributions, and functions are treated the same in BAli-Phy because all of them have parameters or arguments. Parameters have names in BAli-Phy. Parameter values are specified using square brackets as follows:

hky85(kappa=2)            // model
log(x=2)                  // function
normal(mean=0,sigma=1)    // probability distribution

It is possible to specify parameter values by position instead of by name:

hky85(2)
log(2)
normal(0,1)

It is even possible to mix positional and named arguments, as long as all the positional arguments come before all the named arguments:

normal(0,sigma=1)   // OK
normal(mean=0,1)    // not OK

The order and type of parameters for a function can be found with the help command. For example,

% bali-phy help hky85

A value must be given for each parameter, unless the parameter has a default value (See Section 8.4, “Default values and default priors”).

% bali-phy file.fasta -S 'hky85(kappa=2)'
% bali-phy file.fasta -S 'mixture([tn93,hky85(2)])'

Models in phylogenetics literature are often combined using +. For example, the model WAG + F + G4 + I starts with the WAG amino-acid model, and places several modifiers, like " + G4" on the right.

BAli-Phy follows this convention by treating A +> B as an abbreviation for B(submodel=A). When there are multiple '+>' symbols they associate to the left, so that A +> B +> C is understood to mean (A +> B) +> C. For example:

hky85 + Rates.gamma               // rewritten to Rates.gamma(submodel=hky85)
hky85 +> inv                      // rewritten to inv(submodel=hky85)
wag +> f                          // rewritten to f(submodel=wag)
wag +> f +> Rates.gamma +> inv    // rewritten to inv(submodel=Rates.gamma(submodel=f(submodel=wag)))

This allows a simple method for combining models, when one model is an argument to another model.

Some function arguments have default values. For example, the Rates.gamma parameter n has a default value of 4. Thus the following are equivalent:

hky85 +> Rates.gamma(n=4) +> inv
hky85 +> Rates.gamma +> inv

When the default value is random, then the argument has a default prior. For example, the kappa parameter of hky85 has a default value of ~log_normal(log(2),0.25), so the following are equivalent:

hky85(kappa~log_normal(log(2),0.25))
hky85

The help command can be used to determine the default value for a parameter, if there is one.

Every function has a result type, as well as an argument type for each argument. The argument type specifies what kind of arguments are acceptable, and the result type specifies what kind of result the function produces. Types include Int for integers, Double for double-precision floating point numbers, and String for text strings. Integer arguments are implicitly converted to Double when the argument type is Double.

Some types contain parameters. For example List<Int> indicates a list of integers and List<Double> indicates a list of real numbers. In order to indicate a list of unknown type, we use a type variable a and write List<a>. Type variables always begin with a lower-case letter. They are able to match any specific type, and their value is found by pattern-matching. For example, the function x+y takes two arguments of type a and has a result of type a. Thus:

1 + 2        // arguments are a=Int, so result is of type Int
1.0 + 2.0    // arguments are a=Double, so result is of type Double

(a,b) is a parameterized type that can be specialized to (for example) (String,Double) and (Int,Int).

Types for components of substitution models are often parameterized by type of the alphabet. For example, hky85 has a result type of RevCTMC<a>, where a could be DNA or RNA. The use of alphabet types in substitution models prevents combining substitution models with mismatched alphabets.

You should analyze multiple genes under different evolutionary models by putting each one it its own data partition. Placing different genes in different partitions means that their alignments vary independently. It also prevents sequences in one gene from being aligned against sequences in another gene.

Different partitions share the same tree topology and a common set of unscaled branch lengths. However, branch lengths are scaled by a different factor in each partition, since some genes may evolve faster than others.

To put different genes in different partitions, you can place the sequences from each partition in a different FASTA or Phylip file. The sequence names in files for all partitions should be the same.

% bali-phy gene1.fasta gene2.fasta

You can also select different sites from a single larger file:

% bali-phy sequences.fasta:3-350 sequences.fasta:351-570

By default, each partition will have its own substitution model, insertion/deletion model, and scaled tree length. For example, even if all partitions are assigned a tn93 substitution model, their base frequencies will all be estimated independently. When parameters are estimated separately for two partitions, we say that the parameters for those partitions are "unlinked".

A substitution model or insertion-deletion model that is specified without qualification will apply to every partition. However, each partition will recieve its own copy of each model with unlinked parameter values:

% bali-phy sequence-file1 sequence-file2 -S tn93 -I rs07

You can select partition-specific values for 4 options: -S, -I, -A, and --scale. For example, to specify different substitution models but the same alphabet:

% bali-phy sequence-file1 sequence-file2 -S 1:tn93 -S 2:gtr -A DNA

You can fix the alignment and ignore insertion/deletion information in one partition, while allowing the alignment to vary and using insertion/deletion information in another partition:

% bali-phy sequence-file1 sequence-file2 -I 2:none

Since alignments are estimated by default, the alignment will be estimated in the first partition, but fixed in the second partition.

Specifying specify -I none fixes the alignment in all partitions:

% bali-phy sequence-file1 sequence-file2 -I none

You can also specify that two partitions share a single copy of a single substitution model or indel model. For example, if two partitions both have a tn93 model, linking these models would force the partitions to have the same nucleotide frequencies and substitution rates. Linking partitions reduces the number of parameters that need to be estimated, and also pools information between the partitions:

% bali-phy sequence-file1 sequence-file2 -S 1,2:tn93 -I 1,2:rs07

By default each partition has a separate scale, but you can force groups of partitions to share a scale as follows:

% bali-phy sequence-file1 sequence-file2 --scale 1,2:

The --link command is provided to allow specifying a model for each partition separately, and then afterwards choose which partitions to link.

% bali-phy sequence-file1 sequence-file2 -S 1:tn93 -S 2:tn93 --link=1,2 -t
% bali-phy sequence-file1 sequence-file2 -S tn93             --link=1,2 -t   

If the linked partitions are given different models, BAli-Phy will give an error and refuse to run:

% bali-phy sequence-file1 sequence-file2 -S 1:tn93 --link=1,2 -t
bali-phy: Error! Partitions 1 and 2 cannot be linked because they have differing values 'tn93' and ''

You can also specify which of the 3 attributes "smodel", "imodel", and "scale" are being linked:

% bali-phy sequence-file1 sequence-file2 --link=1,2:smodel,scale -t    // Don't link the indel model

BAli-Phy can reconstruct ancestral sequences for all internal nodes. Unlike some programs, BAli-Phy explicitly infers the presence and absence of characters in ancestral sequences. This means that if the ancestral sequence has no character for a column, the reconstructed ancestors will have a gap there. BAli-Phy reconstructs ancestors for fixed-alignment partitions as well as variable-alignment partitions, but it won't write out fixed alignment samples unless you add the flag --set write-fixed-alignments=true. Additionally, if you have an ambiguous character such as N in an observed sequence bali-phy can impute this character, but will not do so unless you add the flag --set infer-ambiguous-observed=true.

BAli-Phy can now reconstruct ancestral sequences for a given tree topology and (leaf sequence) alignment. This is similar to the ancestor-reconstruction that is usually done for fixed-alignment analyses. However, it is not quite the same, because of uncertainty in the tree and the alignment. When computing the probability of an ancestral residue, this summary averages over uncertainty in the topology, the alignment, and the ancestral state itself.

# Construct the leaf sequence alignment to annotate using posterior decoding
% cut-range dir-1/C1.P1.fastas dir-2/C1.P1.fastas --skip=burn-in | alignment-chop-internal --tree c50.tree | alignment-max > P1-max.fasta
# Construct the tree topology to annotate
% trees-consensus dir-1/C1.trees dir-2/C1.trees | tree-tool - --strip-internal-names --name-all-nodes > c50.tree
# Reconstruct ancestral sequences on the given tree and alignment
% summarize-ancestors P1.max.fasta -A dir-1/C1.P1.fastas -T dir-1/C1.trees -A dir-2/C1.P1.fastas -T dir-2/C1.trees -n c50.tree -g c50.tree > P1.ancestors.fasta

BAli-Phy uses an alignment estimate (here, P1-max.fasta) as a template to construct a summary alignment with ancestral sequences. BAli-Phy doesn't condition on the alignment columns, because (i) many columns occur only once in a posterior sample and (ii) conditioning on the column gives too much weight to the template alignment.

Because the alignment is uncertain, residues in the same column of the template alignment may end up in different columns in an MCMC sample. Therefore, in a given MCMC sample, different leaf residues in the same column may have different ancestors at the same internal node! However, in our ancestral reconstruction, a given column may only display a single ancestral residue.

BAli-Phy addresses this problem by averaging across the different ancestral residues in each column of the template alignment. When identifying the ancestral character to column C from a sampled alignment A, we random select a residue in C and use it to select a column from A. This procedure has the nice property that it will yield the traditional ancestral residue prediction if the alignment column is fixed.

Ancestral sequences are written as part of the alignment matrix in each iteration. Ancestral sequences are given names starting with the letter A. For example, in the following alignment, the sequences A5, A6, and A7 are reconstructed ancestors:

>Halobacterium
-T-TAAGGCGGCCATAGCGGTGGGGTTACTCCCGTAC
>Pyrococcus
GG-TACGGCGGTCATAGCGGGGGGGCCACACCCGGTC
>Sulfolobus
GC-CCACCCGGTCACAGTGAGCGGGCAACACCCGGAC
>Homo
GTCTACGGC---CATACCACCCTGAACGCGCCCGATC
>Escherichia
TG-CCTGGCGGCCGTAGCGCGGTGGTCCCACCTGACC
>A5   
GG-CAAGGCGGCCATAGCGGGGGGGCCACACCCGGCC
>A6   
GT-CAAGGCGGCCATAGCGGGGGGGCTACACCCGGTC
>A7   
GT-CAAGGCGGCCATAGCGGGGGGGCTACACCCGGTC

Sampled alignments for the nth partition are in the file. C1.Pn.fastas.

Ancestral states in these alignments are randomly sampled from their joint posterior and do not represent the most probable ancestral state. The alignment of ancestral sequences is also inferred, so these sequences may contain gaps. The length of ancestral sequences may vary between samples when the length of the ancestral sequence is uncertain.

Each sampled alignment matrix corresponds to a tree in the file C1.trees that is written in the same iteration. This tree specifies the phylogenetic location of each ancestral sequence by labelling the internal nodes of the tree. For example, the tree below shows where the internal nodes A5, A6, and A7 are located on the tree:

(Halobacterium:0.213240,((Escherichia:0.435762,Pyrococcus:0.122678)A5:0.114725,Sulfolobus:0.427210))A6:0.042527,Homo:0.427026))A7;

While the tree is written every iteration, the alignment is only written every 10 iterations (by default) in order to save disk space. One method for extracting the trees that correspond to saved alignments is to extract every 10th tree with the program bali-subsample:

% bali-subsample 10 < C1.trees > C1.10.trees

In Bayesian phylogenetic analyses, the tree is not fixed. Therefore the internal node corresponding to the ancestral sequence you wish to reconstruct may not exist in every iteration. The standard Bayesian approach to tree uncertainty is to reconstruct the ancestor for each node by conditioning on the existence of that node in the tree. This allows the reconstructed ancestor for each node to average over uncertainty about the existence of other nodes.

BAli-Phy additionally allows the researcher to condition on branches, since a branch condition is less restrictive. BAli-Phy does not run a separate MCMC chain with a tree constraint for each node, but instead performs conditioning by selecting samples from a single run that satisfy the condition.

% trees-consensus dir-1/C1.trees dir-2/C1.trees | tree-tool - --strip-internal-names --name-all-nodes > c50.tree
% extract-ancestors -A dir-1/C1.P1.fastas -T dir-1/C1.trees -A dir-2/C1.P1.fastas -T dir-2/C1.trees -n c50.tree -g c50.tree > P1.ancestors.fastas

name = taxon1 taxon2 ... taxonN

Inferences about ancestral sequence are sometimes done by analyzing the summary alignment and its ancestral sequences. A more correct approach is to analyze each sampled alignment separately, and then average the results. (This yield the posterior mean of the analysis result.) When this approach is feasible, it is more statistically rigorous than analyzing the summary alignment, because it operates on joint reconstructions and because it incorporates uncertainty in the phylogeny, alignment, and ancestral sequences.

Analyzing the samples and pooling the results also allows each sampled alignment to be analyzed in combination with the corresponding sampled tree.

Convergence refers to the the tendency of a Markov chain to to "forget" its starting value and become typical of its equilibrium distribution. Note that convergence is a property of the Markov chain itself, not of individual runs of the Markov chain. Ideally a number of individual runs should be examined in order to determine how many initial iterations to discard as "burnin".

In MCMC, each sample is not fully independent of previous samples. In fact, even after a Markov chain has converged, it can get "stuck" in one part of the parameter space for a long time, before jumping to an equally important part. When this happens, each new sample contributes very little new information, and we need to obtain many more samples to get good precision on our parameter estimates. In such a case, we say that the chain isn't "mixing" well.

Potential Scale Reduction Factors check that different runs have similar posterior distributions. Only numerical variables may have a PSRF. To calculate the PSRF for each numerical parameter, you may run:

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

The PSRF is a ratio of the width of the pooled distribution to the average width of each distribution, and should ideally be 1. The PSRF is customarily considered to be small enough if it is less than 1.01.

We compare the PSRF based on the length of 80% credible intervals (Brooks and Gelman 1998) and report the result as PSRF-80%CI. For integer-valued parameters, we avoid excessively large PSRF values by subtracting 1 from the width of the pooled CI.

We also report a new PSRF that is more sensitive for integer distributions. For each individual distribution, we find the 80% credible interval. We divide the probability of that interval (which may be more than 80%) by the probability of the same interval under the pooled distribution. The average of this measure over all distributions gives us a PSRF that we report as PSRF-RCF.

This convergence diagnostic gives a criterion for detecting when a parameter value has stabilized at different values in several independent runs, indicating a lack of convergence. This situation might occur if different runs of the Markov chain were trapped in different modes and failed to adequately mix between modes.

Show basic information about the alignment:

% alignment-info file.fasta
% alignment-info file.fasta file.tree

To select columns from an alignment:

% alignment-cat -c1-10,50-100,600- file.fasta > result.fasta
% alignment-cat -c5-250/3 file.fasta > first_codon_position.fasta
% alignment-cat -c6-250/3 file.fasta > second_codon_position.fasta

To concatenate two or more alignments:

% alignment-cat file1.fasta file2.fasta > all.fasta

Remove columns without a minimum number of letters:

% alignment-thin --min-letters=5 file.fasta > file-thinned.fasta

Remove sequences by name:

% alignment-thin --remove=seq1,seq2 file.fasta > file2.fasta

Remove short sequences:

% alignment-thin --longer-than=250 file.fasta > file-long.fasta

Remove sequences with <= 5 differences from the closest other sequence:

% alignment-thin --cutoff=5 file.fasta > more-than-5-differences.fasta

Like --cutoff, but stop when we have the right number of sequences:

% alignment-thin --down-to=30 file.fasta > file-30taxa.fasta

Protect some sequences from being removed:

% alignment-thin --down-to=30 file.fasta --protect=seq1,seq2 > file-30taxa.fasta

Remove sequences that are missing conserved columns:

% alignment-thin --remove-crazy=10 file.fasta > file2.fasta

Draw an alignment to HTML, optionally coloring residues by AU.

%  alignment-draw file.fasta --show-ruler --color-scheme=DNA+contrast > file.html
%  alignment-draw file.fasta --show-ruler --AU=file-AU.prob --color-scheme=DNA+contrast+fade+fade+fade+fade > file-AU.html

Find the last (or first) FastA alignment in a file.

% alignment-find --first < file.fastas > first.fasta
% alignment-find < file.fastas > last.fasta

Turn columns from a template alignment into alignment constraints:

% alignment-indices template.fasta > constraints.txt
% alignment-indices -c100-110,200,300- template.fasta > constraints.txt

Each line in this file corresponds to one alignment column.

Remove internal-node ancestral sequences from an alignment. (This probably only works for alignments output by bali-phy.)

% alignment-chop-internal file.fasta > file-chopped.fasta

This program analyzes the tree sample contained in file. It reports the MAP topology, the supported taxa partitions (including partial partitions), and the majority consensus topology.

Usage: trees-bootstrap file1 [file2 ... ] --predicates predicate-file [OPTIONS]

This program analyzes the tree samples contained in file1, file2, etc. It gives the support of each tree sample for each predicate in predicate-file, and reports a confidence interval based on the block bootstrap.

Each predicate is the intersection of a set of partitions, and is specified as a list of partitions or (multifurcating) trees, one per line. Predicates are separated by blank lines.

Usage: trees-to-SRQ predicate-file [OPTIONS] trees-file

This program analyzes the tree samples contained in trees-file. It uses them to produce an SRQ plot for each predicate in predicate-file. Plots are produced in gnuplot format, with one point per line and with plots separated by a blank line.

If --mode sum is specified, then a "sum" plot is produced instead of an SRQ plot. In this plot, the slope of the curve corresponds to the posterior probability of the event. If the --invert option is used then the slope of the curve correspond to the probability of the inverse event. This is recommended if the probability of the event is near 1.0, because the sum plot does not distinguish variation in probabilities near 1.0 well.

In order to compile BAli-Phy, you need

We recommend the GNU C++ Compiler (GCC) version 10.0 (or higher) or the Clang compiler version 13 or higher. The Cairo graphics library is optional, but if it is missing, the drawtree tool that is used to draw consensus trees won't be built. See also Section 2.8, “Install programs used for viewing the results”.

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

% sudo apt-get install meson
% dpkg -s meson | grep Version
Version: 1.0.1-5

% conda create -n devel -c conda-forge --strict-channel-priority
% conda activate devel
% conda install meson gxx boost-cpp cmake pkg-config cairo
% export BOOST_ROOT=$CONDA_PREFIX

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

% xcode-select --install

% brew install git meson cairo pandoc

% 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

First check out the code using git:

% git clone https://github.com/bredelings/BAli-Phy.git
% cd BAli-Phy

Then run meson to configure the build process:

% meson setup build --prefix=$HOME/Applications/bali-phy-4.0-beta13/ --buildtype=release

Finally, build and install the software:

% ninja -C build test
% ninja -C build install

The command bali-phy and its associated tools should then be located in ~/Applications/bali-phy-4.0-beta13/bin/. To install to another directory dir, specify --prefix=dir to meson.

You can select the C++ compiler by setting the CXX variable. A useful example of this is to use g++-10 on systems where g++ invokes a compiler that is too old:

% CXX=g++-10 meson setup build --prefix=$HOME/Applications/bali-phy-4.0-beta13 --buildtype=release

You may also set compiler and linker options using the CPPFLAGS, CXXFLAGS, and LDFLAGS variables. For example, you can instruct the compiler to use all the features of your chip, instead of producing generic code that will run anywhere:

% CXXFLAGS="-mtune=native -march=native" meson setup --prefix=$HOME/Applications/bali-phy-4.0-beta13

For example, you can set the CPPFLAGS and LDFLAGS variables to instruct the compiler where to look for libraries, such as cairo:

% CPPFLAGS="-I/usr/local/include" LDFLAGS="-L/usr/local/lib" meson setup build --prefix=$HOME/Applications/bali-phy-4.0-beta13 --buildtype=release

Another useful example of this is to produce an OS X executable on that can run on older versions of OS X:

% CXXFLAGS="-mmacosx-version-min=10.9" LDFLAGS="-mmacosx-version-min=10.9" meson setup build --prefix=$HOME/Applications/bali-phy-4.0-beta13 --buildtype=release