# BPP Documentation ## BPP Explained ### Overview Bayesian Phylogenetics & Phylogeography (BPP) is a Bayesian Markov chain Monte Carlo (MCMC) program for analyzing DNA sequence alignments of multiple loci from multiple closely-related species under the multispecies coalescent (MSC) model ([Yang 2002](https://doi.org/10.1093/genetics/162.4.1811); [Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645)). See [Xu 2016](https://doi.org/10.1534/genetics.116.190173) and [Rannala et al 2020](https://inria.hal.science/PGE/hal-02535622) for reviews of the MSC model. The program requires 3 input files: 1. Sequence data file (see [Sequence File](#sequence-file)). 2. Imap file specifying the population source of each sequence (see [Imap file](#imap-file)). 3. Control file specifying other variables needed to run the program (e.g., initial species tree, prior distributions, etc)(see [Control File](#control-file)). This manual will provide details for creating these input files (describing syntax and permissible values for variables). The input files should be in text format and edited using a text format editor (for example, a Unix editor such as emacs, or vi). A basic set of map and control files using sensible default priors based on the sequence data can also be created using the web tool *Minimalist BPP* available at These auto-produced files can then be used as a template for more editing. ### Possible types of analyses The program can be used to conduct four different types of analyses, each specified using two variables in the control file: - **A00** (`speciesdelimitation = 0`, `speciestree = 0`): estimation of species divergence times and population sizes under the MSC, IM and MSC-I models when the species tree is specified by the user ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645); [Flouri et al 2020](https://inria.hal.science/PGE/hal-02536475)) - **A01** (`speciesdelimitation = 0`, `speciestree = 1`): inference of the species tree when sequences are assigned to species by the user ([Rannala and Yang 2017](https://doi.org/10.1093/sysbio/syw119)) - **A10** (`speciesdelimitation = 1`, `speciestree = 0`): species delimitation using a user-specified guide tree ([Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107); [Rannala and Yang 2013](https://doi.org/10.1534/genetics.112.149039)) - **A11** (`speciesdelimitation = 1`, `speciestree = 1`): joint species delimitation and species tree inference using unguided species delimitation [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279) !!! tip "Which analysis should I use?" - **A00**: You know the species tree and want to estimate θ and τ - **A01**: You know the species assignments but want to infer the tree - **A10**: You have a guide tree and want to test species boundaries - **A11**: You want to jointly infer both tree and species boundaries A detailed guide to the use of each of these 4 methods of analysis, as well as example control files, are provided in [Methods of Analysis](#methods-of-analysis). The BPP tutorial also provides examples illustrating the four types of analyses ([Yang 2015](https://doi.org/10.1093/czoolo/61.5.854); [Flouri et al 2020](https://inria.hal.science/PGE/hal-02536475)). ### Model Parameters The basic parameters in the MSC model include the **species divergence times** ( $\tau$ ) and mutation scaled **population sizes** $\theta = 4N\mu$ where $N$ is the effective population size and $\mu$ is the mutation rate per site per generation. Note that $\theta$ specifies the average proportion of sites that have different bases when comparing two sequences sampled at random from the population. Both $\tau$ and $\theta$ are measured in units of expected number of mutations per site ( $\tau$ can be converted to units of years if the substitution rate per site per year is known). For example, if $\tau$ is an estimated divergence time in units of expected substitutions per site (from a BPP analysis) and the substitution rate per site per year is $\mu$, then the estimated divergence time in years is, $\tau'$, is $\tau' = \frac{\tau}{\mu}$ and the diploid effective population size is $N = \frac{\theta}{4\mu}$ For a species tree with $s$ species, there are ($s – 1$) divergence times ($\tau$) and at most ($2s – 1$) population size parameters ($\theta$). Analysis A00 estimates those parameters when the species delimitation and species tree are fixed (specified by the user). Analyses A01, A10, and A11 estimate parameters and model probabilities of the different MSC models. The multispecies-coalescent-with-introgression (MSC-I) model, first implemented in BPP version 4.1 and described in [The MSC-I Model](#the-msc-i-model), adds introgression nodes to the tree, each with a new parameter, the **introgression probability** ($\varphi$) see [Flouri et al 2020](https://doi.org/10.1093/molbev/msz296). The multispecies-coalescent-with--migration (MSC-M) model, first implemented in BPP version 4.6 and described in [The MSC-M Model](#the-msc-m-model), adds a matrix of **instantaneous migration rates** (**M**). For reviews of the MSC model, see [Yang 2014](https://global.oup.com/academic/product/molecular-evolution-9780199602605?lang=es&cc=gb) (Chapter 9), [Xu 2016](https://doi.org/10.1534/genetics.116.190173) and [Rannala et al 2020](https://inria.hal.science/PGE/hal-02535622). ### Model Assumptions The MSC model implemented in BPP makes two basic assumptions about the data: *no recombination between sites of the same locus*, and *free recombination between loci*. The original MSC model ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645)) also assumes complete genetic isolation between populations, but more recent implementations allow gene flow between species (either continuous migration or episodic introgression); the MSC-I model implemented in BPP allows episodic introgression [Flouri et al 2020](https://doi.org/10.1093/molbev/msz296) and the IM model allows continuous migration (see [Introgression and Migration Models](#introgression-and-migration-models) and [Isolation With Migration Models](#isolation-with-migration-models)). The default model assumes a strict molecular clock (constant substitution rates among lineages) and relaxed-clocks (with substitution rates varying among lineages) are implemented through variable clock models (see [Substitution Models](#substitution-models)). ### Nuclear genomic data Ideal data for analysis using the BPP program are loosely-linked short genomic segments. Such data are likely to satisfy the model assumptions because intra-locus recombination is rare for short segments (500 to 1000bp) and distantly spaced segments undergo frequent recombination and thus have nearly independent histories. Genealogical trees are assumed to be a product of neutral evolution (not influenced by natural selection). However, protein-coding gene sequences appear to be useable in BPP analyses since most proteins are performing similar functions in closely related species and the main effect of purifying selection on nonsynonymous mutations is a reduction of the neutral mutation rate ([Shi and Yang 2018](https://doi.org/10.1093/molbev/msx277); [Thawornwattana et al 2018](https://doi.org/10.1093/molbev/msy158)). It is a good idea to separate the noncoding and coding regions of the genome into two datasets. The program allows the option of either a constant rate (strict molecular clock) or variable rates of substitution (relaxed molecular clocks) among lineages. Different substitution models are available, with the most complex (parameter rich) being the General Time-Reversible (GTR) model and the simplest the Jukes-Cantor (JC69) mutation model. All the models correct for multiple hits at individual sites (see [Substitution Models](#substitution-models)). The use of the JC69 model and the strict molecular clock model should be limited to closely related species with sequence divergence not much higher than 10%. Sequences from distantly related species should be analyzed using a GTR substitution model in combination with one of the relaxed clock models (see [Among-species rate variation](#among-species-rate-variation). Models allowing substitution rate variation among sites (see [Among-site rate variation](#among-site-rate-variation)) and among loci (see [Among-locus rate variation](#among-locus-rate-variation))are also implemented in BPP (see [Models of Substitution Rate Variation](#models-of-substitution-rate-variation)). ### Mitochondrial data The MSC model accommodates genealogical heterogeneity across the genome. Thus, different regions of the autosomal genome may have different gene tree topologies and branch lengths (coalescent times). While the mutational process may also vary along the genome, this heterogeneity is believed to be much less important. Thus multiple genes from the mitochondrial genome, which in most species does not undergo recombination, should be treated as one 'locus' in the MSC-based analysis. The mitochondrial genome often has a mutation rate that differs from the nuclear genome, as well as a different effective population size from that of autosomes. While the model of locus-rate variation (option variable `locusrate`) and the heredity scalar (option variable `heredity`) are designed to deal with this, it may be prudent to analyze the loci from the nuclear genome separately from the single mitochondrial locus. ### Identical sequences and phylogenetic signal The model implemented in BPP assumes that the sequences represent random samples from the different species. *Sequences from the same species that are identical should all be used*. It is incorrect to use only the unique haplotypes, which will lead to biased parameter estimates. Similarly, it is incorrect to filter loci based on bootstrap support values and use only those loci with a high "phylogenetic information" signal. ### How to Cite BPP If you publish results obtained using BPP you should cite the version number of the BPP release you used as well as the canonical citation, which as of release v4.3.8 is ``` Tomáš Flouri, Xiyun Jiao, Bruce Rannala, Ziheng Yang (2018) Species Tree Inference with BPP Using Genomic Sequences and the Multispecies Coalescent, Molecular Biology and Evolution 35: 2585--2593. ``` You can also cite the one of the BPP tutorials and the original papers describing the methods you used (see above). Be sure to state the priors that you used since they are necessary for reproducibility. If you conduct a joint analysis of species delimitation and species tree inference, your method description may look like the following (replace the specific values in bold with those you used): > "Joint Bayesian species delimitation and species tree estimation was conducted using the program > BPP (Flouri et al., 2018; Yang, 2015). The method uses the multispecies coalescent model to compare > different models of species delimitation (Yang and Rannala, 2010; Rannala and Yang, 2013) and species > phylogeny (Yang and Rannala, 2014; Rannala and Yang, 2017) in a Bayesian framework, accounting > for incomplete lineage sorting due to ancestral polymorphism and gene tree-species tree discordance. > The population size parameters ($\theta$s) are assigned the **inverse-gamma prior IG(3, 0.002)**, with mean > **0.002/(3–1) = 0.001**. The divergence time at the root of the species tree ($\tau_0$) is assigned the > **inverse-gamma prior IG(3, 0.004)**, with mean **0.002**, while the other divergence time parameters are specified > by the uniform Dirichlet distribution (Yang and Rannala, 2010, eq. 2). Each analysis is run at least twice > to confirm consistency between runs.” ### Features New To BPP 4.0 - Multiple threads - Introgression models (MSC-I) - Migration models (MSC-M) - Local-clock and locus-rates models (clock and locusrate) - Topological constraints (constraint and outgroup) ## Installing and Running BPP This manual applies to BPP versions 4.1 and later. The program is written in C and executables are available for Windows, Mac OSX and Linux. Alternatively, the BPP program can be compiled for Unix (Linux, BSD, etc), Mac OSX, or Windows. If you are using a Windows, Mac OSX or Linux operating system (and you are not a programmer who prefers to compile from source) you should download the latest release executables for your machine (go to [Obtaining BPP](#obtaining-bpp)). If you have not used a command line program before, read [Using the Command Line](#using-the-command-line) then go to [Obtaining BPP](#obtaining-bpp). If you want to compile the program yourself (you must do this if you are using a Unix variant other than Linux or Mac OSX) go to [Compiling the BPP Program](#compiling-the-bpp-program). ### Obtaining BPP Executables and source code for the most recent BPP release are always available at: Download the distribution file containing an executable for your operating system and uncompress the file. Change the current directory to be the root of the subdirectory created by uncompressing the file. For example, on a Linux machine you could type the following to download and install bpp version 4.8.0 (the latest version may be different from this): ``` wget https://github.com/bpp/bpp/releases/download/v4.8.0/bpp-4.8.0-linux-x86_64.tar.gz tar -xzf bpp-4.8.0-linux-x86_64.tar.gz cd bpp-4.8.0-linux-x86_64 ``` You are now ready to proceed to [BPP Trial Run](#bpp-trial-run) ### Using the Command Line BPP is a command-line program, so the preferred way of running it is to open a terminal application and execute the program from the terminal command line, rather than double-clicking on the executable file in your file explorer. If you have not used the command line before, please work through one of the following short tutorials first: Windows: Mac OSX: ### Compiling the BPP Program #### Requirements To compile the program you will need to have a C compiler and the Make program installed on your machine. To test whether this software is installed (on a Unix machine) you can type: ``` cc --version; make --version ``` which should produce output similar to the following if a compiler is installed: ```{ .yaml .no-copy } cc (Ubuntu 9.3.0-10ubuntu2) 9.3.0 Copyright (C) 2019 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. GNU Make 4.2.1 Built for x86_64-pc-linux-gnu Copyright (C) 1988-2016 Free Software Foundation, Inc. License GPLv3+: GNU GPL version 3 or later This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law. ``` If instead you see an error message such as: ```{ .yaml .no-copy } Command 'cc' not found ``` then you (or your system) administrator will need to install the compiler software (using a package manager such as apt on Ubuntu, for example) before you can compile BPP. #### Downloading and compiling the source code If you have the necessary software for compiling you can download the C source code for the latest release from: The link to the source file on github is named `Source code`. You will need to download and uncompress the file, then change to the source code subdirectory before executing the commands outlined below. The program needs to be compiled only once. For example, the following commands use the gcc compiler to compile the program and move the generated executable file (bpp) into the bin/ folder. ``` cd bpp mkdir bin cd src make mv bpp ../bin ``` If you use git you can instead clone the bpp repository and check out the master branch (which contains source code for the latest stable version of bpp) and then compile the program: ``` git clone https://github.com/bpp/bpp.git cd bpp git checkout master cd src make ``` You are now ready to proceed to [BPP Trial Run](#bpp-trial-run) ### BPP Trial Run Run the program from the command line (rather than double-clicking the executable) so that you will see any potential error messages. Change your current directory to the top level of the directory created by uncompressing the bpp distribution file. For example, in Linux the following commands will uncompress the distribution file and move you to the top level of the bpp directory: ``` tar -xvzf bpp-4.8.0-linux-x86_64.tar.gz cd bpp-4.8.0-linux-x86_64 ``` In the bpp/ subdirectory, run the program in Windows by typing the following command within the Terminal application: ``` bin\bpp --cfile examples\frogs\A00.bpp.ctl ``` In Linux or Mac OSX type the following commands within a terminal: ``` cd examples/frogs/ ../../bin/bpp --cfile examples/frogs/A00.bpp.ctl ``` If the program executed successfully, you should see initial output similar to the following: ```{ .yaml .no-copy } bpp v4.8.0_linux_x86_64, 31GB RAM, 12 cores https://github.com/bpp/bpp Auto-selected SIMD ISA: AVX2 Starting timer.. Using seed: -1 Parsing species tree... Done Parsing phylip file... Done ``` Alternatively, you may see one or more error messages. This may be due to a spelling error, or you may not be in the correct subdirectory. If this does not appear to be the case and the problem persists you can ask for help on the BPP discussion group at: When posting on the forum, please specify the exact command you typed and the error message you received (preferably by copying and pasting this information from your terminal or attaching a screen shot). ### BPP Command-line Options Below is a list of command-line options for running BPP. **General options:** | Command-line option | Meaning | |:-------------------------|:--------------------------------------------------------------------------------------------------------| | | With no options, bpp prints out the version number and computer hardware information | | `--help` | Display help information | | `--version` | Display version information | | `--quiet` | Only output warnings and fatal errors to stderr | | `--cfile FILENAME` | Run analysis for the specified control file | | `--simulate FILENAME` | Run simulation for the specified control file | | `--resume FILENAME` | Resume analysis from a specified checkpoint file | | `--msci-create FILENAME` | Construct an MSci graph using a definitions file | | `--summary FILENAME` | Summarize results using specified control file (no MCMC) | **Advanced options:** | Command-line option | Meaning | |:--------------------------|:-------------------------------------------------------------------------------------------------------| | `--arch SIMD` | Force specific vector instruction set (default: auto) | | `--bfdriver FILENAME` | Create control files for marginal likelihood calculation | | `--points INTEGER` | Number of G-L quadrature points (used with --bfdriver) | | `--no-pin` | Do not pin threads to cores | | `--theta_mode INTEGER` | Definition of $\theta$ step lengths (default: 2). 1: one step length for all $\theta$; 2: one step length for tip and one for inner nodes; 3: one step length for each node | | `--theta-prop STRING` | Proposal distribution for $\theta$ gibbs move: 'mg_invg' or 'mg_gamma' | | `--theta-showeps` | Show individual step lengths/pjumps for each $\theta$ | | `--theta-slide-prob FLOAT`| Frequency for $\theta$ sliding window move (default: 0.1) | | `--phi-slide-prob FLOAT` | Frequency for $\phi$ sliding window move (default: 0.1) | | `--wrate_mode INTEGER` | Definition of W step lengths (default: 1) | | `--wrate-slide-prob FLOAT`| Frequency for W sliding window move (default: 0.1) | | `--wrate-showeps` | Show individual step lengths/pjumps for each W | | `--extend INTEGER` | Extend resumed analysis by number of MCMC samples (use with --resume) | | `--keep-labels` | Keep original node labels when summarizing results | ## Input File Formats Two input files are required for every BPP analysis, the sequence file and the control file. If more than one species is being analyzed an Imap file is also required. A few additional input files are required only when specific options are specified in the control file. Here we describe the content and format specifications of the sequence, Imap and control files that will be needed for most BPP analyses. ### Sequence File The sequence alignments must be in the phylip/paml format, with one alignment following the other, all in one file. An alignment is called a locus. Every locus must have at least 2 sequences, different loci can have different numbers of sequences, and some species may have no sequences present at some loci. In the phylip/paml format the first line of the entry for each locus contains two numbers specifying the number of sequences and the number of sites, respectively, that are present at the locus. A sequence name must precede each sequence entry and be separated from it by either a line break or two or more spaces. Each sequence name must contain an individual ID, which is unique among the sequences listed for a given locus, preceded by a caret `^` symbol. For example, the sequence name `GI01234567^Specimen1` has the individual ID `Specimen1`. The label `GI01234567` that precedes the `^` is optional, so `^Specimen1` is also a valid sequence name but `Specimen1` is not. The sequences may be in interleaved format. DNA bases can be specified using either upper or lower case letters, missing data can be specified using a `?` symbol, and IUPAC ambiguity codes may be used. The interpretation of an ambiguity code (as either an uncertain base or a diploid genotype) depends on whether the data are specified as phased in the control file. By default (control file variable `cleandata = 0`, see [Control File](#control-file)), alignment gaps and ambiguity nucleotides are used in the likelihood calculation, with gaps treated as question marks (see [Yang 2014](https://global.oup.com/academic/product/molecular-evolution-9780199602605?lang=es&cc=gb) pp. 111-112). If `cleandata = 1`, all columns with gaps or ambiguity characters are removed before analysis. An example of a sequence data file for 2 loci and 4 sequences per locus is given below (you can also examine the sequence files `frogs.txt` and `yu2001.txt` in the `examples` subdirectory). ``` 4 20 dog^alfy1 attgtgccctctctctctca dog^alfy2 attgtgccctctctctctca cat^smoot attgtgccctctagctctca human^ben attgtgccctctgactctcc 4 10 dog^alfy1 attgtgccct dog^alfy2 attgtgccct fox^sample1 attgtgccct ^ben attgtgaagt ``` Imap File --------- Each sequence name has an individual ID tag after the caret `^` symbol. For example, the sequence name `GI01234567^Specimen1` has the individual ID `Specimen1`. An **Imap file** maps the individual IDs (specimens) to their species/populations. For example, the following Imap file maps `Specimen1` to species A, `Specimen4` to species B, and so on. ``` Specimen1 A Specimen2 A Specimen3 A Specimen4 B Specimen5 C Specimen6 C ``` See also the map files `frogs.Imap.txt` and `Rokas2003-5species-Imap.txt` in the `examples` subdirectory. The MSC model implemented in BPP infers the relationships among species/populations not individuals and so it uses only the population ID for a sequence (such as A, B, ... in the example), and does not use the individual ID. Also when BPP reads the sequence names, it use the individual ID to retrieve the species ID for each sequence but then ignores the sequence name. One motivation for the use of the Imap file and this two-layer design is that one may wish to analyze the same sequence data with different population/species assignments of individual IDs. This can be achieved by editing the small Imap file rather than editing the larger sequence data file. ![imap](https://github.com/bpp/bpp-manual/blob/main/figures/fig-imap.png?raw=true) An example Imap file showing the sequence ID in the Imap file corresponding to that in the sequence file. Note that the caret `^` symbol is not allowed in the sequence ID of the Imap file but is required in the sequence label in the sequence data file. Outgroups and Constraints File ------------------------------ The **outgroups and constraints file** is optional and its name is specified in the control file by setting variable `constraintfile` to the name of the constraint file. Topological constraints on species trees can only be specified for analyses A01 (species tree inference with a fixed species delimitation) and A11 (joint species tree inference and species delimitation). See [Methods of Analysis](#methods-of-analysis) for a detailed description of the A01 and A11 methods of analysis. Note that species trees are always rooted in BPP, whether you use the clock or relaxed clock. The control variable `constraintfile` in the control file has the format ``` constraintfile = constraints.txt ``` Inside the constraint file, three keywords are allowed: `define`, `constraint`, and `outgroup`. The `define` keyword is used to assign a name or alias to a clade, `constraint` defines a clade or subtree, and `outgroup` means that species not on the list (the ingroup species) form a clade. ``` define g1 as (G,H); # identical to constraint (D,E,F,(G,H)); constraint (D,E,F,g1); outgroup A,B,C; constraint (A,B); ``` Either the equal sign (=) or spaces can be used to separate the keyword from the specified value. In the example above, there are eight species in the dataset: A, B, C, D, E, F, G, and H. Line 1 defines an alias `g1` for clade (G,H), which may be used in subsequent `define` or `constraint` statements. Line 3 specifies A, B, C as the outgroups; this is interpreted to mean that species not on the list, (that is, D, E, F, G and H) form a clade. Since this information is already in line 2, line 3 has no effect in this case. The modifier keyword `NOT` (case-insensitive) can be used with `define` to instead specify the complement of the species on the list. Suppose in the data we have 10 species of ratites (flightless birds), and two outgroup species, chicken and ostrich. You can specify the constraint as follows: ``` define ratites as NOT (chicken, ostrich); constraint = (chicken, (ostrich, ratites)); ``` Note that BPP interprets `outgroup` to mean that all species not on the list form a clade. In theory the `outgroup` keyword is unnecessary because you can always achieve the same thing by using `define` (especially with the `NOT` option) and `constraint`. However `outgroup` may be more convenient or explicit in some situations. If there are four species (A,B,C,D), ``` outgroup = B,C,D; ``` achieves nothing since all the 15 rooted species trees are allowed, whereas ``` outgroup = C,D; ``` means (A,B) form a clade so that 3 rooted species trees are allowed. Finally if multiple compatible constraints are specified, BPP merges them into one so that ``` constraint (E,F,G,H); constraint (F,G,H); constraint (G,H); ``` is equivalent to ``` constraint (E,(F,(G,H))); ``` If any two constraints are in conflict, BPP will abort with an error message. For example the following will cause an error. ``` constraint = (E,F,G,H); constraint = (G,H); constraint = (F,G); ``` Date File ------------------------------ The **date file** is optional and its name is specified in the control file by setting variable `datefile` to the name of the date file. When using tip-dating, each sample has an associated sample age. The `datefile` assigns an age to the sample using the individual ID tag. The dates are in units before time present, so larger numbers are older. Any desired time units (e.g. years, thousands of year, etc) can be used so long as the mutation rate prior is on the same time scale (e.g. expected substitutions per year, expected substitutions per thousand years). For example, the following date file assigns `specimen1` an age 500, `specimen2` an age of 10000. The population name can be included before the caret `^` symbol or it can be excluded. ``` A^specimen1 500 specimen2 10000 B^specimen3 45000 B^specimen4 35000 ``` Each line should have one individual followed by the sample age. Each sequence in the dataset must be included in the date file. Also see the date file `mammoth/dates.txt` in the `examples` subdirectory. Control File ------------ The **control file** specifies most aspects of a BPP analysis, including the names of input and output files, choices of models, prior distributions on parameters, etc. The format is line oriented, with each line specifying the value of a particular variable. The default control file name is `bpp.ctl`. Lines beginning with `*` or `#` are comments and are ignored by the program. Most often the order of the lines is unimportant, but there are exceptions (as explained below). Each non-commented line of the control file contains a variable name and an assignment of one or more values to the variable. The general syntax of a line in the control file is: ``` variable = value ``` where `variable` is the name of a valid variable and `value` is a list of one or more permissible values for the variable. Most variable assignments follow this format --- the exceptions are described later. BPP has complicated dependencies among variables, both in terms of whether the variable needs to be defined (or has a default value when undefined) and what its permissible values are. Here we begin by defining all the variable names, the syntax of their arguments and permissible types (integers, floating point numbers, etc), and outline the structure of dependencies among variables | Symbol/Rule | Definition | |-------------|--------------------------------------------------------------| | b | Boolean (0,1) | | +d | Positive integer | | f | floating point number | | s | String (or character) | | t | Tree in Newick format | | x* | List of one or more elements of type x | | y [z] | z is a set of variables whose form is determined by y | | (w,z) | w and z are two different permissible types | | y [(w,z)] | whether w or z is the permissible type depends on value of y | **Table 1.** Symbols used to represent syntax and value types of variables in variable definitions. Detailed definitions of all the control file variables are provided at the end of this section, ordered according to their indexes in Table 2 below. If you are impatient to get started using the program with worked examples, you might want to now skip ahead to [Example control file](#example-control-file) which explains the control file variables in the context of a simple example dataset analysis. You can later consult this section regarding particular control file variables as the need arises. | Index | Variable | Values | Condition | Dependencies | |-------------------------------|-----------------------|---------------------|--------------------|--------------| | [1](#1-seed) | *seed* | (-d,+d) | True | None | | [2](#2-usedata) | *usedata* | b | True | None | | [3](#3-jobname) | **jobname** | s | True | None | | [4](#4-seqfile) | **seqfile** | s | True | <-16 | | [5](#5-finetune) | **finetune** | b [f*] | True | ->7 | | [6](#6-print) | **print** | b b b b b | True | None | | [7](#7-burnin) | `burnin` | +d | True | <-5 | | [8](#8-sampfreq) | **sampfreq** | +d | True | None | | [9](#9-nsample) | **nsample** | +d | True | None | | [10](#10-species&tree) | **species&tree** | +d s* +d* t | True | ->{11,15,24} | | [11](#11-imapfile) | `Imapfile` | s | 10[1]>1 | <-10 | | [12](#12-speciesdelimitation) | *speciesdelimitation* | b [b f*] | True | ->14 | | [13](#13-speciestree) | *speciestree* | b | True | ->13 | | [14](#14-speciesmodelprior) | `speciesmodelprior` | +d | 12[1]=1 OR 13[1]=1 | <-\{12,13\} | | [15](#15-phase) | **phase** | b* | 10[1] | <-10 | | [16](#16-nloci) | *nloci* | +d | True | ->4 <-30 | | [17](#17-model) | *model* | s [s] | True | None | | [18](#18-qrates) | *Qrates* | b +f +f +f +f +f +f | 17=''7'' | <-17 | | [19](#19-basefreqs) | *basefreqs* | b +f +f +f +f | 17=''7'' | <-17 | | [20](#20-alphaprior) | *alphaprior* | f f +d | True | None | | [21](#21-cleandata) | *cleandata* | b | True | None | | [22](#22-thetaprior) | *thetaprior* | s [(f*, f f s)] | True | None | | [23](#23-tauprior) | *tauprior* | s f f | True | None | | [24](#24-phiprior) | `phiprior` | f f | 10[4] | <-10 | | [25](#25-locusrate) | *locusrate* | +d [(f f f s, s, )] | True | <-26 ->26 | | [26](#26-clock) | *clock* | +d [f f f s s] | True | <-25 ->25 | | [27](#27-heredity) | *heredity* | +d [(f f, s)] | True | None | | [28](#28-checkpoint) | *checkpoint* | +d [(+d)] | True | None | | [29](#29-constraintfile) | *constraintfile* | s | True | None | | [30](#30-threads) | `threads` | +d [ +d +d ] | True | ->16 | | [31](#31-datefile) | *datefile* | s | 25[3] | None | **Table 2.** Complete list of BPP control file variables. The column `Variable` contains the control file variable name, the column `Values` specifies (using symbols) the format of the values (symbols are defined in Table 1). The column `Condition` describes the conditions (if any) needed for the variable to be required. An entry of True in this column means that the variable is always required/defined although possibly not explicitly (there may be a default value). The column `Dependencies` lists the variables that either influence the variable (<-) or are influenced by the variable (->) with the variables denoted using their index from column `Index`. Table 2 lists all the variables that can be defined in the BPP control file. The first column lists an index used to identify the variable in other columns of the table and to determine the order in which the variables are defined in subsequent text. The second column lists the variable name. The variable names are highlighted to indicate the structure of their dependencies among variables as follows: - Bold = variable always needs to be defined, has no default value, and does not depend on other variables. - Italic = variable always needs to be defined, has a default value if no value is provided in control file, and does not depend on other variables. - Solid background = whether variable needs to defined (and possibly also its permissible values) depends on the value of another variable. - Regular = variable always needs to be defined and has no default value but its permissible values depend on the value of another variable. Column 3 lists the permissible types of each variable and the syntax for specifying arguments in defining the variable. The variable types are defined in Table 1. Column 4 lists the condition (state of another variable) that either determines whether the variable on that line exists (solid background variables) or determines its permissible values (italic or regular variables). x [y] refers to value of argument y of the variable with index x. Column 5 lists indexes of variables that either are influenced by the variable (the notation ->[c,d] indicates that the variable influences variables c and d) or variables that influence the variable (the notation <-[a] indicates that variable a influences the variable). If you are writing a control file from scratch, these dependencies necessitate that you specify the bold variables first, followed by the italic (only needed if changing defaults) and then the solid background and regular. ## BPP control file variables ### Quick Reference Tables The tables below provide a quick overview of all control file options. See the detailed sections that follow for complete documentation. #### Basic Setup | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `seed` | `seed = -1` or `seed = N` | `-1` | Random seed (-1 = auto-generate) | | `usedata` | `usedata = 0\|1\|2` | `1` | 0: prior only, 1: likelihood+prior, 2: prior with fixed gene tree | | `jobname` | `jobname = prefix` | — | Output file prefix (replaces deprecated `outfile`/`mcmcfile`) | | `seqfile` | `seqfile = path` | — | Path to sequence alignment file | | `Imapfile` | `Imapfile = path` | — | Path to individual-to-species mapping file | | `nloci` | `nloci = N` | — | Number of loci to analyze | #### Analysis Type | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `speciesdelimitation` | `0` or `1 [algorithm options]` | `0` | 0: fixed delimitation, 1: estimate delimitation | | `speciestree` | `0` or `1` | `0` | 0: fixed tree, 1: estimate tree | | `species&tree` | `N name1 name2 ... counts ... (newick);` | — | Species names, sample counts, and tree | **Analysis modes:** A00 (0,0), A01 (0,1), A10 (1,0), A11 (1,1) #### MCMC Settings | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `burnin` | `burnin = N` | — | Number of burn-in iterations | | `sampfreq` | `sampfreq = N` | — | Sample every N iterations | | `nsample` | `nsample = N` | — | Number of samples to collect | | `finetune` | `finetune = 0\|1: params...` | `1` | Auto-tune (1) or manual step sizes | | `checkpoint` | `checkpoint = N N` | — | Checkpoint frequency and number to keep | | `threads` | `threads = N [N N]` | `1` | Number of threads [start pin] | | `loadbalance` | `none` or `zigzag` | `none` | Thread load balancing strategy | #### Priors | Option | Syntax | Example | Description | |--------|--------|---------|-------------| | `thetaprior` | `invgamma a b [e]` | `invgamma 3 0.002` | Prior on θ (population size) | | | `gamma a b` | `gamma 2 100` | | | | `beta p q l u` | `beta 1 1 0 0.01` | | | `tauprior` | `invgamma a b` | `invgamma 3 0.03` | Prior on τ (divergence times) | | | `gamma a b` | `gamma 2 100` | | | `phiprior` | `beta a b` | `beta 1 1` | Prior on φ (introgression probability) | | `wprior` | `wprior = a b` | `wprior = 2 200` | Prior on w (migration rate) — *replaces migprior in v4.8.0* | #### Substitution Model | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `model` | `model = name` | `HKY` | JC69, K80, F81, HKY, T92, TN93, F84, GTR, or amino acid models | | `Qrates` | `Qrates = 0\|1\|2 [a b]` | `0` | Rate matrix: 0=fixed, 1=estimate, 2=prior | | `basefreqs` | `basefreqs = 0\|1\|2 [params]` | `0` | Base frequencies: 0=empirical, 1=estimate, 2=prior | | `alphaprior` | `alphaprior = a b` | — | Gamma prior on α (rate variation among sites) | | `cleandata` | `cleandata = 0\|1` | `0` | Remove ambiguous sites | #### Rate Variation | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `clock` | `clock = 1\|2\|3\|4 [params]` | `1` | 1: strict, 2: independent rates, 3: correlated, 4: simple | | `locusrate` | `locusrate = 0\|1\|2\|3 [params]` | `0` | 0: equal rates, 1-3: variable rates among loci | | `heredity` | `heredity = 0\|1\|2 [params]` | `0` | Heredity multipliers for θ | #### MSC-I Model (Introgression) *(v4.1+)* | Option | Syntax | Description | |--------|--------|-------------| | `hybridization` | `source target as node1 node2 tau=y,y phi=f` | Define introgression event | | `bidirection` | `pop1 pop2 as node1 node2 phi=f,f` | Bidirectional introgression | #### MSC-M Model (Migration) *(v4.6+)* | Option | Syntax | Example | Description | |--------|--------|---------|-------------| | `wprior` | `wprior = a b` | `wprior = 2 200` | Default gamma prior on w *(v4.8.0+, replaces migprior)* | | `migration` | `migration = N` | `migration = 2` | Number of migration connections | | | `source target [a b] [a_w] [pseudo_a pseudo_b]` | `A C 2 100` | Migration specification per connection | | `geneflow` | `geneflow = 0\|1` | `geneflow = 1` | Enable gene flow estimation | #### Output Control | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `print` | `print = N N N N N` | `1 0 0 0 0` | Control output: samples, locusrate, heredity, labels, genetrees | | `printlocus` | `printlocus = N l1 l2 ...` | — | Print gene trees for specific loci | #### Other Options | Option | Syntax | Default | Description | |--------|--------|---------|-------------| | `phase` | `phase = 0\|1 per species` | `0` | Analytical phasing for diploid data | | `speciesmodelprior` | `speciesmodelprior = N` | `1` | Prior on species models (A10, A11) | | `constraintfile` | `constraintfile = path` | — | Topological constraints file | | `datefile` | `datefile = path` | — | Tip dates for dated-tip analysis | | `thetamodel` | `linked-none\|all\|inner\|msci\|mscm` | `linked-none` | Link θ parameters across populations | | `traitfile` | `traitfile = path` | — | Trait data for morphological analysis | --- ### 1 seed ---------------- ``` seed = (-d,+d) ``` **DESCRIPTION** Specifies the seed used for the random number generator. **VALUES** `-d`, use the wall clock (in Windows) or white noise (in Linux) to generate a seed (store seed in file `SeedUsed`). `+d`, use +d as the seed. **DEFAULT** `-1` **COMMENTS** Using the same positive integer seed will produce identical results in different runs, which is useful for debugging, and using different positive integers in different runs produces different results. Using $-d$, for example, $-1$, the computer's entropy generator is used as a source for the seed, and different runs will produce different results. When evaluating MCMC convergence by comparing results across different runs be sure to use -d. The positive integer seed automatically generated when using option -d is stored in a file named SeedUsed and can be used explicitly to replicate a result. !!! tip "Best Practice" Run each analysis at least twice using different seeds (`seed = -1`) to confirm that results are stable across runs. This helps verify MCMC convergence. **EXAMPLES** ``` seed = -1 seed = 278 ``` ### 2 usedata ------------------------------------------------------------------------ ``` usedata = b ``` **DESCRIPTION** Specifies whether data (likelihood+priors) are used in calculating probabilities during MCMC or only priors. **VALUES** `0`, use only the priors to calculate probabilities (likelihood constant). `1`, use likelihood and prior probabilities. `2`, use only the priors (like `0`) but fix the gene tree topology. **DEFAULT** `1` **COMMENTS** The `0` option can be used for debugging, or examining priors. When using option `0` a MCMC run produces samples from the prior for each variable. Option `2` is similar to `0` but additionally fixes the gene tree topology during the MCMC run. **EXAMPLES** ``` usedata = 0 usedata = 1 usedata = 2 ``` ### 3 jobname ------------------------------------------------------------------------ ``` jobname = s ``` **DESCRIPTION** Defines the job name, which serves as a prefix for all output files generated by the analysis, including the main output file and the MCMC sample file. **VALUES** `s`, a string specifying the common prefix for all output files. This string may include a directory path. **EXAMPLES** ``` jobname = /home/mickey/mouse ``` The above will create several files such as: ``` /home/mickey/mouse.txt # main output file /home/mickey/mouse.mcmc.txt # MCMC sample file /home/mickey/mouse.SeedUsed # seed information file ``` Additional files may be generated depending on the type of analysis. **NOTE:** The options `mcmcfile` and `outfile` from previous versions are now obsolete and have been replaced by `jobname`. Using these deprecated options will result in an error. ### 4 seqfile ------------------------------------------------------------------------ ``` seqfile = s ``` **DESCRIPTION** Sets the name of the file containing the sequence data to be the string `s` **VALUES** `s`, a string of characters specifying the directory path and name of file that contains the sequence data **EXAMPLES** ``` seqfile = /home/mickey/mouse_seq.txt seqfile = sequences.txt ``` ### 5 finetune ------------------------------------------------------------------------ ``` finetune = b [key:val ...] ``` **DESCRIPTION** Determines whether step lengths in MCMC proposals are automatically optimized or manually set to fixed values. Step lengths can optionally be specified using a key:value format. !!! warning "Syntax Change in v4.8.1" The finetune syntax changed in BPP v4.8.1. The old positional syntax (`finetune = 1: 5 0.001 0.001 ...`) is no longer supported. Use the new key:value format or simply `finetune = 1` for default auto-tuning. The available proposal keys are: | Key | Parameter Proposal Description | |------|---------------------------------------------------------------| | Gage | Node age on gene tree | | Gspr | Subtree pruning and regrafting move | | tau | Tau (divergence time) | | mix | Mixing step (jointly changing mu, tau and gene tree node ages)| | lrht | Locus rate and heredity | | phis | Introgression probability (phi) | | pi | Stationary frequencies of substitution model | | qmat | Q matrix elements of substitution model | | alfa | Alpha parameter of gamma model for rate variation among sites | | mubr | Mean branch rate (mu-bar) | | nubr | Variance of branch rates (nu-bar) | | mu_i | Locus-specific mu | | nu_i | Locus-specific nu | | brte | Branch rate | **VALUES** `0` Use fixed step lengths (no auto-tuning). Default step lengths are used unless overridden with key:value pairs. `1` Automatically optimize step lengths during burnin. This is the recommended setting. `0 key:val ...` Use fixed step lengths with specified values. `1 key:val ...` Auto-tune with specified initial step lengths. **DEFAULT** `finetune = 1` (auto-tune with default initial values) **DEPENDENCIES** If auto-tuning is enabled (`b = 1`), then `burnin` must be > 200. **COMMENTS** For most analyses, simply using `finetune = 1` is recommended. The program will automatically adjust step lengths during burnin to achieve good mixing. Custom step lengths are only needed for difficult analyses with poor convergence. **EXAMPLES** ``` finetune = 1 # auto-tune with defaults (recommended) finetune = 0 # fixed defaults, no auto-tuning finetune = 1 Gage:5 Gspr:0.001 # auto-tune with custom initial values finetune = 0 Gage:5 Gspr:0.001 mix:0.3 # fixed custom values ``` ### 6 print ------------------------------------------------------------------------ ``` print = b b b b b ``` **DESCRIPTION** Specifies what information is saved to the output files during the MCMC. **VALUES** There are 5 boolean (0,1) variables that specify whether particular variables are saved to file (1) or not saved (0). From left to right: variable 1 specifies MCMC samples; variable 2 specifies locus-specific rate parameters (rate $\mu_i$, variance parameter $\nu_i$ and species-tree branch rates for locus $r_{ij}$); variable 3 specifies locus-specific heredity scalars if those are estimated from the data; variable 4 specifies locus-specific gene trees; and variable 5 specifies locus-specific substitution-rate parameters (qmat for the $Q$ matrix, freqs for base frequencies, and alpha for gamma rates for sites). **EXAMPLES** ``` print = 0 1 1 1 0 print = 1 1 1 1 1 ``` ### 7 burnin ------------------------------------------------------------------------ ``` burnin = +d ``` **DESCRIPTION** Specifies the number of burn-in iterations that are executed in the MCMC run before sampling begins. **VALUES** `+d`, a positive integer specifying the number of burn-in iterations. **DEPENDENCIES** If `finetune = 1` (automatic optimization) then `burnin` has to be $>200$. **COMMENTS** The total number of MCMC iterations is `burnin` + `nsample` $\times$ `sampfreq`. **EXAMPLES** ``` burnin = 100 burnin = 10000 ``` ### 8 sampfreq ------------------------------------------------------------------------ ``` sampfreq = +d ``` **DESCRIPTION** Specifies the interval at which samples from the MCMC are to be written to the MCMC output file specified by `jobname`. **VALUES** `+d`, a positive integer specifying the interval between samples from the MCMC. **COMMENTS** The total number of MCMC iterations is `burnin` + `nsample` $\times$ `sampfreq`. **EXAMPLES** ``` sampfreq = 2 sampfreq = 10 ``` ### 9 nsample ------------------------------------------------------------------------ ``` nsample = +d ``` **DESCRIPTION** Specifies the number of samples from the MCMC that are to be written to the MCMC output file specified by `jobname`. **VALUES** `+d`, a positive integer specifying the number of MCMC samples. **COMMENTS** The total number of MCMC iterations is `burnin` + `nsample` $\times$ `sampfreq`. **EXAMPLES** ``` nsample = 5000 nsample = 25499 ``` ### 10 species&tree ------------------------------------------------------------------------ ``` species&tree = +d s* +d* t ``` **DESCRIPTION** Specifies the number of species, the species names, the maximum number of sequences (at any locus) for each species, and the species tree in Newick format (or species graph in extended Newick format in models with introgression). **VALUES** `+d`, a positive integer specifying the total number of species. `s*`, a list of `+d` strings specifying the name of each species. `+d*`, a list of `+d` positive integers specifying the maximum number of sequences (at any locus) for each species. `t`, a tree (or graph) in Newick format (or Newick extended format) specifying the relationships among species (and possibly introgression events). **DEPENDENCIES** The variable `Imapfile` must be specified only when the number of species `+d` > 1. The length of the boolean array specified for variable `phase` is equal to the number of species `+d`. The variable `phiprior` must be specified only if `t` is an introgression graph. **COMMENTS** See the subsection *Notation for trees and introgression graphs* for details on the format of `t`. The role of `t` in the MCMC analysis depends on the value of the `speciestree` variable. If `speciestree=1` the variable `t` specifies the starting tree for the MCMC analysis, otherwise if `speciestree=0` it is the fixed tree used for estimating parameters or the fixed guide tree for species delimitation. **EXAMPLES** ``` species&tree = 4 H C G O 1 1 1 1 (((H, C), G), O); ``` ### 11 Imapfile ------------------------------------------------------------------------ ``` Imapfile = s ``` **DESCRIPTION** Sets the path/name of the Imap file to be the string `s`. **VALUES** `s`, a string of characters specifying the directory path and name of a file that contains the species/sequence map information. **DEPENDENCIES** The variable `Imapfile` must be specified only when the number of species defined by the `species&tree` variable `+d` > 1. **COMMENTS** See [Imap File](#imap-file) for a complete description of the Imap file format. **EXAMPLES** ``` Imapfile = /home/foo/bar_Imap.txt Imapfile = Imap.txt ``` ### 12 speciesdelimitation ------------------------------------------------------------------------ ``` speciesdelimitation = b [b f*] ``` **DESCRIPTION** Specifies whether species delimitation is performed and if so specifies the parameters of the species delimitation rjMCMC proposal algorithm. **VALUES** `0`, specifies no species delimitation. `1 0 epsilon`, specifies species delimitation with rjMCMC algorithm 0, with parameter $\epsilon$ as specified in equations 3 and 4 of [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107). Reasonable values for $\epsilon$ are 1, 2, 5, etc. `1 1 alpha m`, specifies rjMCMC algorithm 1, with parameters $\alpha$ and $m$ as specified in equations 6 and 7 of [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107). Reasonable values are $\alpha = 1, 1.5, 2$, etc. and $m = 0.5, 1, 2$. **DEFAULT** `0` **DEPENDENCIES** If `speciesdelimitation = 1` then `speciesmodelprior` must be defined. **EXAMPLES** ``` speciesdelimitation = 1 0 2 speciesdelimitation = 1 1 2 1 ``` ### 13 speciestree ------------------------------------------------------------------------ ``` speciestree = b ``` **DESCRIPTION** Specifies whether the species tree is estimated or fixed in the MCMC analysis. **VALUES** `0`, specifies that species tree is fixed. `1`, specifies that species tree is estimated. **DEFAULT** `0` **DEPENDENCIES** If `speciestree = 1` then `speciesmodelprior` must be defined. **COMMENTS** This option invokes the nearest-neighbor interchange (NNI) or subtree-pruning-regrafting (SPR) algorithm to propose changes to the species tree topology during the MCMC run. **EXAMPLES** ``` speciestree = 1 speciestree = 0 ``` ### 14 speciesmodelprior ------------------------------------------------------------------------ ``` speciesmodelprior ``` **DESCRIPTION** Specifies the prior distribution for rooted species trees and species delimitations. **VALUES** `0`, specifies equal probabilities for the labeled histories (rooted trees with the internal nodes ordered by age). `1`, specifies equal probabilities for the rooted species trees. `2`, specifies equal probabilities for the numbers of species ($1/s$ each for $1,2,\ldots,s$ species given $s$ populations) and divides the probability for any specific number of species among the compatible models (of species delimitation and species phylogeny) in proportion to the number of labeled histories. `3`, specifies equal probabilities for the numbers of species ($1/s$ each for $1,2,\ldots,s$ species given $s$ populations) and divides the probability for any specific number of species among the compatible models (of species delimitation and species phylogeny) uniformly. **DEPENDENCIES** The variable `speciesmodelprior` must be defined if either (or both) `speciestree = 1` or `speciesdelimitation = 1`. If `speciestree = 1` and `speciesdelimitation = 0` then only options `speciesmodelprior = 0` or `speciesmodelprior = 1` may be used. **COMMENTS** The priors specified by options 2 and 3 are described in [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279). The prior specified by option 3 may be suitable when there are many populations. **EXAMPLES** ``` speciesmodelprior = 0 speciesmodelprior = 3 ``` ### 15 phase ------------------------------------------------------------------------ ``` phase = b* ``` **DESCRIPTION** Specifies whether the sequence data for each species is phased. **VALUES** `1`, in position $i$ of the list indicates that data is unphased for species $i$ in the list of species in variable `species&tree`, otherwise `0` indicates it is phased. **DEFAULT** `0` **DEPENDENCIES** The number of boolean (0,1) variables in the list must equal the number of species `d+` specified in the `species&tree` variable. **COMMENTS** If at position $i$ the `phase` variable is 1, each sequence from species $i$ is treated as an unphased diploid sequence, with heterozygous sites represented using the ambiguity characters `YRMKSW`. The missing data tokens `N-?` are allowed and are all treated as ?, but other ambiguity characters (`HBVD`) are not allowed. For phasing, each sequence is expanded/resolved into two sequences, and the program averages over possible phase resolutions according to the approach of @Gronau2011. For example, a sequence `R...Y` with two heterozygous sites, `R` (meaning `A/G`) and `Y` (meaning `C/T`), has two possible phase resolutions: 1. `A...T` and `G...C` 2. `A...C` and `G...T` A sequence with no heterozygous sites will be phased/resolved into two identical sequences. The MSC model effectively averages over different phase resolutions of the heterozygous sites in performing the likelihood calculation. The option may increase memory usage and CPU time considerably if many sequences exist with many heterozygote sites at a locus. If the `phase` variable for a species is 0, all sequences from that species are treated as resolved haplotype sequences, and ambiguities are interpreted in the usual way. For example, an ambiguity code `Y` at a position is interpreted to mean that one sequence exists with an uncertain nucleotide at that position that is either a T or a C. The program does not allow some sequences from a species to be phased and other sequences from that same species to be unphased. **EXAMPLES** ``` phase = 0 0 0 0 1 phase = 1 1 1 1 1 ``` ### 16 nloci ------------------------------------------------------------------------ ``` nloci = +d ``` **DESCRIPTION** Specifies the number of loci present in the sequence data file specified by variable `seqfile` that will be analysed. **VALUES** `+d`, a positive integer specifying the number of loci to be analyzed. **DEFAULT** If `nloci` is not specified, all the loci present in the sequence file specified by `seqfile` will be analysed. **DEPENDENCIES** `nloci` must be less than or equal to the number of loci available in the sequence file specified by `seqfile`. **COMMENTS** If `nloci` is less than the number of loci present in the file specified by `seqfile` only the first `nloci` present in that file will be analysed. **EXAMPLES** ``` nloci = 1 nloci = 2000 ``` ### 17 model ------------------------------------------------------------------------ ``` model = s [s] ``` **DESCRIPTION** Specifies the DNA substitution model (or amino acid substitution model) that will be used. For DNA sequences the simplest (default) model that can be chosen is the Jukes-Cantor (JC69) model and the most complex is the General Time-Reversible (GTR) model. A range of models of intermediate complexity are also available. For amino acid sequences, most of the commonly used amino acid substitution matrix models are available. If a model other than `Custom` is specified it will apply to all loci. If `Custom filename` is specified, an additional file `filename` must be provided that contains one or more lines with the format: ``` loci_indices data_type model ``` where `loci_indices` is an index or range of indexes for loci (indexed according to their order of occurrence in the sequence data file), `data_type` is either `DNA` or `AA`, and model specifies the substitution model variable for the indexed loci as described above. For example, the entry ``` 1-10 DNA JC69 11 AA DAYHOFF ``` specifies that loci 1 to 10 comprise DNA sequences that will have the `JC69` model applied to them, while locus 11 is an amino acid sequence that will have the `DAYHOFF` model applied to it. **VALUES** For DNA sequences: `JC69`, `K80`, `F81`, `HKY`, `T92`, `TN93`, `F84`, `GTR`. For amino acid sequences: `DAYHOFF`, `LG`, `DCMUT`, `JTT`, `MTREV`, `WAG`, `RTREV`, `CPREV`, `VT`, `BLOSUM62`, `MTMAM`, `MTART`, `MTZOA`, `PMB`, `HIVB`, `HIVW`, `JTTDCMUT`, `FLU`, and `STMTREV`. Among-locus model variation: `Custom filename`. **DEFAULT** `JC69` **DEPENDENCIES** If `model = GTR` the variables `Qrates` and `basefreqs` can be optionally defined which provide priors on the rate matrix and base frequencies, respectively (default priors will be used otherwise). **EXAMPLES** ``` model = GTR model = CPREV model = Custom models.txt ``` ### 18 Qrates ------------------------------------------------------------------------ ``` Qrates = b +f +f +f +f +f +f ``` **DESCRIPTION** Specifies whether the exchangeability parameters in the GTR model are fixed or variable and either their fixed values or prior relative mean values. **VALUES** `b`, specifies either fixed (1) or variable (0) exchangeability parameters in the GTR model. `+f +f +f +f +f +f`, specify the exchangeability parameters in the order $a, b, c, d, e, f$ for TC, TA, TG, CA, CG, and AG, as described in [Yang 2014](https://global.oup.com/academic/product/molecular-evolution-9780199602605?lang=es&cc=gb). **DEFAULT** `0 1 2 1 1 2 1` **DEPENDENCIES** The variable `Qrates` is defined only if `model = GTR`. **COMMENTS** Two example cases are discussed here with either fixed, or variable rates. The following specifies fixed rates: ``` Qrates = 1 2 1 1 1 1 2 ``` This specifies fixed rates at $a = 2, b = c = d = e = 1$ and $f = 2$, so that the transition rate is twice as high as the transversion rate, and the model corresponds to K80 or HKY. Note that only the relative rates matter, because the rate matrix ($Q$) is scaled so that the average rate (over the base frequencies) is 1 and branch lengths are measured in the expected number of mutations/substitutions per site. Nevertheless the program expects six rate parameters. ``` Qrates = 0 10 5 5 5 5 10 ``` This specifies variable rates with a prior on $a, b, c, d, e$ that is a Dirichlet distribution for every locus. The above specifies Dir(10, 5, 5, 5, 5, 10), with parameters $\alpha_a = 10, \alpha_b = \alpha_c = \alpha_d = \alpha_e = 10$, and $\alpha_f = 10$. The rates generated from the Dirichlet sum to 1, and they are rescaled. Note that larger values for those parameters mean less variance, while the mean for $a$, say, is given by $\alpha_a/\alpha$ with $\alpha = \alpha_a + \alpha_b + \alpha_c + \alpha_d + \alpha_e + \alpha_f$. This specifies an average transition/transversion rate ratio of 2 (if base frequencies are all fixed at $\frac{1}{4}$). **EXAMPLES** ``` Qrates = 1 2 1 1 1 1 2 Qrates = 0 1 1 1 1 1 1 ``` ### 19 basefreqs ------------------------------------------------------------------------ ``` basefreqs = b +f +f +f +f ``` **DESCRIPTION** Specifies whether the base frequency parameters are fixed or variable and either their fixed values or prior relative mean values. **VALUES** `b`, specifies either fixed (1) or variable (0) base frequency parameters in the GTR model. `+f +f +f +f`, specify the base frequency parameters in the order TCAG. **DEFAULT** `0 1 1 1 1` **DEPENDENCIES** The variable `basefreqs` is defined only if `model = GTR`. **COMMENTS** The parameters $\pi_T, \pi_C, \pi_A, \pi_G$ can be either fixed or variable with a prior specified by the Dirichlet distribution with $\alpha$ parameters specifying the relative mean frequencies. **EXAMPLES** ``` basefreqs = 1 0.15 0.35 0.15 0.35 basefreqs = 0 10 10 10 10 ``` ### 20 alphaprior ------------------------------------------------------------------------ ``` alphaprior = f f +d ``` **DESCRIPTION** Specifies the prior probability density for the shape parameter controlling the pattern of among-site rate variation. **VALUES** `f f +d`, specify the $\alpha_p$ and $\beta_p$ parameters of the prior on the Gamma shape parameter and the number of rate categories ncatG, respectively. **DEFAULT** $\alpha=\infty$ (no rate variation among sites) **COMMENTS** Among-site rate variation is assumed to follow a Gamma distribution with a mean of 1 and shape parameter $\alpha$. The variance of rates among sites is $1/\alpha$ and therefore a smaller $\alpha$ specifies more rate variation. A discrete distribution is used to efficiently approximate the Gamma distribution and ncatG is the number of rate categories used in the approximation -- more rate categories provide a better approximation but incur greater computational expense. A value for ncatG of 4 is recommended. The prior on $\alpha$ is also a Gamma distribution with parameters $\alpha_p$ and $\beta_p$. The variance and mean of the prior on $\alpha$ are: $\textrm{Mean}(\alpha) = \frac{\alpha_p}{\beta_p}$, $\textrm{Var}(\alpha) = \frac{\alpha_p}{\beta_p^2}$ **EXAMPLES** ``` alphaprior = 1 1 4 ``` ### 21 cleandata ------------------------------------------------------------------------ ``` cleandata = b ``` **DESCRIPTION** Specifies whether columns in the alignment which have gaps or ambiguity characters will be removed prior to analysis, or will instead be retained for use in the likelihood calculation. **VALUES** `0`, specifies that columns in the alignment which have gaps or ambiguity characters will be treated as missing data and used in the likelihood calculation. `1`, specifies the removal of columns in the alignment which have gaps or ambiguity characters. **DEFAULT** `0` **EXAMPLES** ``` cleandata = 1 cleandata = 0 ``` ### 22 thetaprior ------------------------------------------------------------------------ ``` thetaprior = s [(f*, f f s)] ``` **DESCRIPTION** Specifies the distribution model and parameters of the prior distribution on the contemporary and ancestral population parameters $\theta$ and whether the MCMC integrates over the parameters or the integration is instead performed analytically. **VALUES** `s`, specifies the form of the prior distribution on theta and should be either `invgamma`, `gamma` or `beta` to specify either an inverse gamma, gamma, or beta distribution, respectively. `[(f*, f f s)]`, specifies the parameters of the prior distribution. The form for the parameters depends on the type of distribution specified in the first argument. The permissible combinations of distributions and parameters are: ``` invgamma a b invgamma a b e ``` which specifies an inverse gamma prior with $\alpha$ and $\beta$ parameters a and b, respectively. The first entry specifies that $\theta$s are integrated over analytically and the second entry (with trailing e) specifies that $\theta$s are estimated. ``` gamma a b ``` which specifies a gamma prior with $\alpha$ and $\beta$ parameters a and b, respectively. ``` beta p q l u ``` which specifies a beta prior constrained to the interval $(l,u)$ with parameters p and q, respectively. All the distributions above, except the inverse Gamma distribution (without the e argument) integrate over the $\theta$ parameters numerically as part of the MCMC, estimating the posterior distribution of $\theta$ for each population. **COMMENTS** The inverse-gamma prior `invgamma a b`, where a and b are the parameters $\alpha$ and $\beta$, has mean and variance: $\textrm{Mean}(\theta) = \frac{\beta}{\alpha - 1}$, $\textrm{Var}(\theta) = \frac{\beta^2}{(\alpha - 1)^2 (\alpha - 2)}$. For example, if $\alpha=3$ and $\beta=0.002$ the mean is $0.002/(3 – 1) = 0.001$ (one variable site per kb on average). Note that all $\theta$ parameters in the MSC model (for both modern species and extinct ancestral species) are assigned the same specified prior distribution with identical parameters. The inverse-gamma is a conjugate prior for $\theta$ (Hey and Nielsen, 2007), which means that both the prior and the posterior of $\theta$ will be inverse-gamma. Use of the conjugate prior allows the $\theta$ parameters to be integrated out analytically, and thus the dimension of the parameter space is reduced. This typically leads to improved mixing of the MCMC. !!! note "Trade-off: Analytical Integration vs. Posterior Estimation" Using `invgamma a b` (without `e`) integrates θ analytically for better MCMC mixing, but you won't get posterior distributions for θ. To estimate θ posteriors, use `invgamma a b e` instead. To estimate the $\theta$ parameters when using an inverse gamma prior, add the letter e (or E) on the line, as follows ``` thetaprior = 3 0.002 e ``` Whether $\theta$s are integrated out or estimated, other results (such as the posterior probability of species trees or species-delimitation models) should be identical if the same prior is used. The analytical integration of $\theta$ is only possible with an inverse gamma prior distribution, if any other prior is used the $\theta$s are automatically estimated and the E flag is not needed and should not be used. The gamma prior `gamma a b`, where a and b are the parameters $\alpha$ and $\beta$, has mean and variance: $\textrm{Mean}(\theta) = \frac{\alpha}{\beta}$, $\textrm{Var}(\theta) = \frac{\alpha}{\beta^2}$. For example, if $\alpha=0.001$ and $\beta=1$ the mean is $0.001/1 = 0.001$ (one variable site per kb on average). The beta prior `beta p q l u` has mean and variance: $\textrm{Mean}(\theta) = \frac{p u + q l}{p + q}$, $\textrm{Var}(\theta) = \frac{p q (u - l)^2}{(p + q)^2 (p + q + 1)}$. For example, with `beta 2 18 1e-6 0.1` the mean is: $\frac{2(0.1)+18(10^{-6})}{2+18} = 0.0100009$ which is 10 variable sites per kb on average. **EXAMPLES** ``` thetaprior = invgamma 3 0.002 thetaprior = invgamma 4 0.001 e thetaprior = beta 2 18 1e-6 0.1 thetaprior = beta 2 18 1e-6 0.01 thetaprior = beta 1 10 0.0001 0.2 thetaprior = gamma 0.001 1 ``` ### 23 tauprior ------------------------------------------------------------------------ ``` tauprior = s f f ``` **DESCRIPTION** Specifies the prior probability distribution of the divergence time parameter for the root in the species tree. **VALUES**\ `s f f`, specifies either the gamma or inverse gamma prior distribution and the $\alpha$ and $\beta$ parameters, respectively for $\tau_0$, the divergence time parameter for the root in the species tree. ``` gamma a b invgamma a b ``` specifies a gamma distribution, or an inverse gamma distribution, respectively with a and b to be the parameters $\alpha$ and $\beta$, respectively. **COMMENTS** If the specified prior on $\tau_0$, the age of the root of the tree is an inverse-gamma distribution IG($\alpha, \beta$) with $\alpha > 2$ the mean and variance are: $\textrm{Mean}(\tau_0) = \frac{\beta}{\alpha - 1}$, $\textrm{Var}(\tau_0) = \frac{\beta^2}{(\alpha - 1)^2 (\alpha - 2)}$. If the prior is a gamma distribution the mean and variance are: $\textrm{Mean}(\tau_0) = \frac{\alpha}{\beta}$, $\textrm{Var}(\tau_0) = \frac{\alpha}{\beta^2}$. The remaining species divergence times are generated from the uniform Dirichlet distribution [@Yang2010 eq. 2]. As an example, suppose that `tauprior = 3 0.03` specifies the prior for $\tau_0$, the divergence time parameter for the root in the species tree. The mean of $\tau_0$ is then $0.03/(3 – 1) = 0.015$ (which means 1.5% of sequence divergence between the root of the species tree and the present time). If the mutation rate is $10^{–9}$ mutations/site/year, this distance will translate to an absolute prior mean divergence time of 15 MY. **EXAMPLES** ``` tauprior = invgamma 3 0.03 tauprior = gamma 0.01 1 ``` ### 24 phiprior ------------------------------------------------------------------------ ``` phiprior = f f ``` **DESCRIPTION** Specifies the parameters of the prior distribution on $\varphi$, the introgression probability under the MSC-I model. **VALUES** `f f`, specify the parameters $\alpha$ and $\beta$, respectively, of the beta distribution prior $\textrm{B}(\alpha,\beta)$ for $\varphi$. **DEPENDENCIES** If the `species&tree` variable `t` specifies a tree with introgression nodes then `phiprior` must be defined. **COMMENTS** The introgression probability $\varphi$ takes values on the interval $(0,1)$ and the prior distribution for $\varphi$ is assumed to be beta $\textrm{B}(\alpha, \beta)$ which has mean and variance: $\textrm{Mean}(\varphi) = \frac{\alpha}{\alpha+\beta}$, $\textrm{Var}(\varphi) = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}$. For example, `phiprior = 1 1` specifies the beta prior $\textrm{B}(\alpha,\beta)$ with $\alpha = 1$ and $\beta = 1$ for $\varphi$, which has a mean of $1/2$ and a variance of $1/12$. **EXAMPLES** ``` phiprior = 1 1 phiprior = 1.5 1.5 ``` ### 25 locusrate ------------------------------------------------------------------------ ``` locusrate = d+[(f f f s, s)] ``` **DESCRIPTION** Specifies the model of substitution rate variation among loci, as well as the parameters, and the prior on the parameters of the model. **VALUES** `d+[(f f f s, s)]`, specifies the model and parameters. The variable d+ takes one of 4 possible values with different numbers of additional variables depending on d+. The permissible combinations are: | **d+** | **[args]** | **Description** | |--------|-----------------------------------------------------------------------|------------------------------------| | 0 | | Loci have identical rates | | 1 | f f f s | Locus rates variable and estimated | | | $\alpha_{\bar{\mu}} \,\, \beta_{\bar{\mu}} \,\, \alpha_{\mu_i}$ prior | | | 2 | s | Locus rates variable and specified | | | Filename of file containing locus rates | | | 3 | f f | Loci have identical rates. The rate is estimated using tip-dating | | | $\alpha$ $\beta$ | | where $\alpha_{\bar{\mu}}$ and $\beta_{\bar{\mu}}$ are the parameters of the Gamma distribution prior on the average rate among loci $\bar{\mu}$ and $\alpha_{\mu_i}$ is the parameter of the prior on the rate at the $i$th locus given $\bar{\mu}$. The prior variable must be either `iid` for a conditional iid (hierarchical) prior on $\mu_i$ or `dir` for a Gamma-Dirichlet prior. The prior variable is optional, if it is not specified `iid` is used. **DEFAULT** `0` **DEPENDENCIES** If both the `prior` option for the variable `locusrate` and the `prior` option for the variable `clock` are set, they should take the same value. If the `locusrate = 3 f f` option is used, `datefile` must be specified. **COMMENTS** The setting `locusrate = 0` (default) means that all loci have the same mutation rate. Setting `locusrate = 1 f f f s` specifies a model with rate variation among loci in which rates are estimated from the sequence data. Including the first integer (1) which species the model there are 5 arguments, the last one (prior) being optional. Parameters $\alpha_{\bar{\mu}}$ and $\beta_{\bar{\mu}}$ specify the shape and rate parameters of the gamma distribution for the mean rate across loci ($\bar\mu$), while $\alpha_{\mu_i}$ and prior are used to specify the locus rates ($\mu_i$) given the mean rate ($\bar\mu$). The option prior can take two values dir (for Gamma-Dirichlet rates for loci), and iid (for conditional i.i.d. or hierarchical prior). `locusrate = 2 LocusRateFileName` specifies the fixed-rates model of locus-rate variation ([Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148)). This is the strategy used by [Yang 2002](https://doi.org/10.1093/genetics/162.4.1811), with the relative rates estimated by the distance to an outgroup species. The relative locus rates are listed in the file: there should be as many numbers in the file, separately by spaces or line returns, as the number of loci (nloci). The program re-scales those rates so that the average across all loci is 1 and then use those relative rates as fixed constants. Specifically the mean rate across loci ($\bar{\mu}$) is assigned a gamma prior: $\bar{\mu} = \mathrm{G}(\alpha_{\bar{\mu}},\beta_{\bar{\mu}}),$ with mean and variance: $\textrm{Mean}(\bar{\mu}) = \frac{\alpha_{\bar{\mu}}}{\beta_{\bar{\mu}}}$, $\textrm{Var}(\bar{\mu}) = \frac{\alpha_{\bar{\mu}}}{\beta_{\bar{\mu}}^2}$. When there are no fossil calibrations in the species tree, the rates should all be relative. In this case we suggest fixing $\alpha_{\bar{\mu}}=\beta_{\bar{\mu}}=0$ (with 0 causing the program to fix both parameters at $\infty$), and the program will fix the mean rate across loci at $\bar\mu = 1$. Otherwise one can use equal and large values such as $\alpha_{\bar{\mu}}=100$ and $\beta_{\bar{\mu}}=100$, so that $\bar\mu$ is nearly fixed at 1. Given the mean rate $\bar\mu$, two priors are available to specify the locus rates $\mu_i$, with $i = 1, 2, \cdots, L$, where $L$ is the number of loci (`nloci`). If `prior = dir` (for Gamma-Dirichlet distribution of locus rates), the total rate $L\bar\mu$, given the mean rate ($\bar\mu$), is partitioned into locus rates ($\mu_i$), using the concentration parameter $\alpha_{\mu_i}$ The model and notation follow [Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148) [eq. 4] and [Dos Reis et al 2014](https://doi.org/10.1093/sysbio/syu020) [eqs. 3-5]. If `prior = iid` (for conditional-i.i.d. or hierarchical prior of locus rates), the locus rates ($\mu_i$) are i.i.d. given the mean rate ($\bar\mu$): $\mu_i \sim \mathrm{G}(\alpha_{\mu_i},\alpha_{\mu_i}/\bar{\mu})$. This model is described in [Zhu et al 2015](https://doi.org/10.1093/sysbio/syu109) [eq. 8] and is also implemented in [mcmctree]{.smallcaps}. In both the Gamma-Dirichlet and the conditional i.i.d. models, parameter $\alpha_{\mu_i}$ is inversely related to the extent of rate variation among loci, with a large $\alpha_{\mu_i}$ meaning similar rates among loci. If all loci are noncoding, the rates are probably similar, so $\alpha_{\mu_i} =$ 10 or 20 may be reasonable, while for coding loci or exons, $\alpha_{\mu_i} =$ 2 or 1 may be appropriate. The $\alpha_{\mu_i}$ parameter may affect the estimates of the population size parameter ($\theta$) for the root node on the species tree ([Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148)). The $L$ locus rates ($\mu_i$) are parameters in the model. If $\alpha_{\bar{\mu}} > 0$, the mean rate ($\bar\mu$) is a parameter as well. The setting `locusrate = 3 f f` means the all loci have the same mutation rate, and the rate is estimated using tip-dating. The rate prior is $\Gamma (\alpha, \beta)$. The units of the prior should match the units in the `datefile`. For example, if the dates are specified in years, the rate should be in expected number of substitutions per year. **EXAMPLES** ``` locusrate = 0 locusrate = 2 rates.txt locusrate = 1 0 0 2 locusrate = 1 2 3 2 dir locusrate = 1 1 1 1 iid locusrate = 3 20 1000000 ``` ### 26 clock ------------------------------------------------------------------------ ``` clock = +d [f f f s s] ``` **DESCRIPTION** Specifies whether a strict molecular clock (equal substitution rates among lineages) is used or instead a specific variable clock model (allowing variation in substitution rates among lineages). If a variable clock model is specified, the priors on the parameters of the variable clock model are also specified. **VALUES** `+d`, specifies a strict clock (1), a variable clock with independent rates among branches (2), a variable clock with autocorrelated rates between ancestral and descendent branches (3), or a simple variable clock model (4). `f f f s s`, are required when `clock =` 2 or 3 and specify respectively the parameters $\alpha_{\bar{\nu}}$, $\beta_{\bar{\nu}}$ and $\alpha_{\nu_i}$, the prior distribution for the locus rate ($\mu_i$) which is either `dir` for Dirichlet or `iid` for conditional i.i.d (hierarchical) prior, and the distribution for the branch rates ($r_{ij}$) given the locus rate ($\mu_i$) which is either `G` for Gamma or `LN` for log-normal. Parameters $\alpha_{\bar{\nu}}$ and $\beta_{\bar{\nu}}$ are the parameters of the gamma distribution for the average of variances across loci ($\bar\nu$), while $\alpha_{\nu_i}$ and `prior` are used to specify the variance ($\nu_i$) for locus $i$ given the average ($\bar\nu$). The specification of the distribution of $\nu_i$ given $\bar\nu$ follows the same procedure as the specification of the distribution of $\mu_i$ given the mean rate $\bar\mu$ when modeling among-locus rate variation; see notes above about the `locusrate` variable. `f f`, are required when `clock = 4` and specify the parameters $\alpha_{\bar{\nu}}$ and $\beta_{\bar{\nu}}$ for the simple rates model. **DEFAULT** `1` **DEPENDENCIES** If both the `prior` option for the variable `locusrate` and the `prior` option for the variable `clock` are set, they should take the same value. **COMMENTS** As an example, the specification (`clock = 2 10.0 100.0 5.0 iid G`) specifies a local-clock model in which the mutation rate drifts over branches, independently among loci. A detailed explanation is as follows. First $\bar\nu$ is assigned a gamma distribution $G(10.0, 100.0)$, with mean 0.1. Given $\bar\nu$, the conditional i.i.d. prior means that $\nu_i \sim G(\alpha_{\nu_i}, \alpha_{\nu_i}/\bar\nu)$ with shape parameter $\alpha_{\nu_i} = 5.0$ and mean $\bar\nu$, for $i = 1, 2, \cdots, L$. Note that $\alpha_{\nu_i}$ is inversely related to the variance of $\nu_i$: use small values of $\alpha_{\nu_i}$ (2, 1, or 0.5) if you believe that $\nu_i$ varies among loci (meaning that the clock nearly holds at some loci but is seriously violated at others). As another example, the specification (`clock = 2 10.0 100.0 5.0 dir LN`) means that $\bar\nu$ is assigned a gamma prior $G(10.0, 100.0)$, and then the sum $L\bar\nu$ is partitioned into $\nu_i$ (for $i = 1, 2, \cdots, L$) according to the Dirichlet distribution (`prior = dir`) with concentration parameter $\alpha_{\nu_i} = 5.0$. Again large $\alpha_{\nu_i}$ means the same extent of clock violation at different loci. Given the locus rate $\mu_i$ (specified using the `locusrate` variable) and the variance $\nu_i$ (specified using the `clock` variable) for locus $i$, the different lineages may have different rates at the locus, and those rates are independent among loci. If distribution specified is `G` (for gamma), the rate for species-tree branch $j$ at locus $i$ has the following gamma distribution $r_{ij} | \mu_i, \nu_i \sim \mathrm{G}(\mu_i^2/\nu_i, \mu_i/\nu_i)$. This has mean $\mu_i$ and variance $\nu_i$. Alternatively if distribution = `LN` (for log-normal), the rate for species-tree branch $j$ at locus $i$ has the following the log-normal distribution $r_{ij} | \mu_i, \nu_i \sim \mathrm{LN}(\mu_i, \nu_i)$, where $\mu_i$ is the mean of the LN distribution and $\nu_i$ is the variance parameter of the lognormal. Note that the rate-drift model specifies rates for branches on the species tree (rather than on the gene tree) for each locus, and gene-tree branches residing in the same population or species have the same rate. For example, if all sequences at a locus are from the same species and all coalescent events occur in that species (before reaching an ancestral species), all branches on the gene tree will have the same rate even if the relaxed-clock model allows rates among species. In contrast, if a gene-tree branch passes several species or populations, the different segments of the branch will have different rates. The branch length on the gene tree is calculated by summing up the lengths of those segments (with the length of each segment being the product of the rate and the time duration for the segment). The option `clock = 3` is similar to clock 2 but specifies the autocorrelated-rates model. `clock = 3` with the log-normal distribution (LN) specifies the geometric Brownian motion model of [Rannala and Yang 2007](https://doi.org/10.1080/10635150701420643). This assigns a rate to each species-tree branch, that is, to the mid-point of the branch. Given the rate at the species-tree root ($\mu_i$ at locus $i$), the rates for the two branches around the root are specified. Then given the rate for each ancestral branch, the rates for its two daughter branches are specified, by integrating over the rate at the internal node that is ancestral to the daughter branches. See figure 1 and equations 3-8 in [Rannala and Yang 2007](https://doi.org/10.1080/10635150701420643). The rates for all species-tree branches are thus assigned through a pre-order tree traversal, starting from the root moving to the tips, until all branches are visited. If `clock = 3` is specified with the gamma distribution (option `G`), the model works as follows, using a similar pre-order tree traversal. First the two branches at the species-tree root have the gamma distribution with mean $\mu_i$ and variance $\nu_i$. Then given the rate for each ancestral branch, the rates for its two daughter branches are specified as independent gamma variables with the mean to be the rate of the parental branch and with the variance to be $\nu_i$. Clock 3 is currently implemented for the MSC model only and is unavailable under the MSC-I model. Note that a larger variance $\nu_i$ implies a more serious violation of the molecular clock at locus $i$. Also note that $\nu_i$ will be similar to $\bar\nu$, especially if $\alpha_{\nu_i}$ is large, and $\bar\nu$ has prior mean $\alpha_{\bar{\nu}} / \beta_{\bar{\nu}}$. If $\alpha_{\bar{\nu}} = 10$ and $\beta_{\bar{\nu}} = 100$, the prior mean will be 0.1. Thus, for the log-normal model $\nu = 0.5$ represents a serious violation of the clock while $\nu < 0.1$ represents only a slight violation. **EXAMPLES** ``` clock = 1 clock = 2 10.0 100.0 5.0 dir LN clock = 3 10.0 50.0 3.0 dir G ``` ### 27 heredity ------------------------------------------------------------------------ ``` heredity = +d [(f f, s, )] ``` **DESCRIPTION** Specifies whether proportional differences exist between all $\theta$ values for different loci, and if so whether the proportionality multipliers (called inheritance scalars) are estimated from the data, or set to specified fixed values. If inheritance scalars are estimated, the parameters of the prior also must be specified. **VALUES** `+d`, specifies whether no differences in $\theta$ exist among loci (0), proportional differences exist and inheritance scalars are estimated (1), or proportional differences exist and inheritance scalars are fixed to specified values (2). `f f`, is only specified when `+d` is 1 and specifies the $\alpha$ and $\beta$ parameters, respectively, of the Gamma prior on the scalar $s_i$ for locus $i$. `s`, is only specified when `+d` is 2 and specifies the name of a file containing $L$ (number of loci) heredity scalars (which must be positive numbers) separated by spaces or linebreaks. **DEFAULT** `0` **COMMENTS** `heredity = 0` is the default and constrains $\theta$ to be equal across loci. `heredity = 1` or `2` specifies two models that allow $\theta$ to vary among loci, which may be useful for analyses combining data from autosomal loci, loci on sex chromosomes, and/or mitochondrial loci. Such mixed data, may be expected to have proportional differences in effective population sizes among loci so that inheritance scalars [Hey and Nielsen 2004](https://doi.org/10.1534/genetics.103.024182) should be applied. For example, in diploid sexual organisms with strict maternal inheritance of mtDNA we expect the effective population size of a mitochondrial locus to be $1/4$ that of a nuclear locus. Other factors such as natural selection may also cause $\theta$ to deviate from the neutral expectation even among autosomal loci. BPP implements two options for allowing such variations. The first option (`heredity = 1`) specifies that locus-specific inheritance scalars $s_i$ be estimated, using a Gamma prior with parameters $\alpha$ and $\beta$ specified by the user. The prior mean and variance of $s_i$ are: $\textrm{Mean}(s_i) = \frac{\alpha}{\beta}$, $\textrm{Var}(s_i) = \frac{\alpha}{\beta^2}$. For example, `(heredity = 1 4 4)` specifies a Gamma prior $\textrm{G}(4, 4)$, with mean $4/4 = 1$ and variance $1/4$, for the inheritance scalar at each locus. The MCMC then generates the posterior distribution of $s_i$ for each locus. The second option (`heredity = 2`) allows the user to specify the $s_i$ in a file so they remain fixed constants during the MCMC run. **EXAMPLES** ``` heredity = 0 heredity = 1 4 4 heredity = 2 scalars.txt ``` ### 28 checkpoint ------------------------------------------------------------------------ ``` checkpoint = +d [(+d)] ``` **DESCRIPTION** Specifies whether one or more checkpoint files are created to allow the BPP program to resume using the current state of the MCMC at the iteration during which the checkpoint file was created. **VALUES** `+d[(+d)]`, specifies the iteration at which the first checkpoint file is created (first argument) and how frequently checkpoint files are subsequently created (second optional argument). **DEFAULT** If the checkpoint value is undefined no checkpoint file is created. **COMMENTS** The checkpoint files are named `JOBNAME.Z.chk` where `JOBNAME` is the name specified by the `jobname` variable and `Z` is the number of the checkpoint file (the first checkpoint file is labeled 1, the second 2, and so on). To resume execution of the program starting from a checkpoint use the BPP `–resume` switch followed by the checkpoint file name. For example, `bpp –resume outrun1.1.chk`. **EXAMPLES** ``` checkpoint = 100000 checkpoint = 100000 10000 ``` ### 29 constraintfile ------------------------------------------------------------------------ ``` constraintfile = s ``` **DESCRIPTION** Specifies the name of a file containing information for placing topological constraints on inferred trees and specifying outgroups for tree rooting. **VALUES** `s`, a string specifying the name of the constraint file. **DEFAULT** If `constraintfile` is undefined no constraints or outgroups are used in the analysis. **EXAMPLES** ``` constraintfile = myconstraints.txt ``` ### 30 threads ------------------------------------------------------------------------ ``` threads = +d [+d +d] ``` **DESCRIPTION** Specifies the number of CPU threads to be used and optionally the specific threads to be used. **VALUES** `+d [+d +d\]`, either 1 or 3 positive integers. A single integer argument specifies the number of threads to use. If three integers are given the first specifies the number of threads, the second specifies the index of the first thread used and the third specifies the increment. For example, `12 4 1` specifies that 12 threads should be used, the first with index 4 and incrementing by one so that threads 4-15 are used. **DEFAULT** If `threads` is undefined, the calculations for each locus are placed on different threads when possible. **DEPENDENCIES** The number of specified threads cannot exceed the number of loci specified. **EXAMPLES** ``` threads = 4 threads = 8 4 1 ``` ### 31 datefile ------------------------------------------------------------------------ ``` datefile = s ``` **DESCRIPTION** Sets the path/name of the date file to be the string s **VALUES** `s`, a string specifying the name of the date file. **DEFAULT** If `datefile` is undefined, no tip dates are used in the analysis. **DEPENDENCIES** Model A00 must be used with tip-dating (`speciesdelimitation = 0`, `speciestree = 0`). Checkpointing cannot be used. The tip-dating `locusrate` option must be used, `locusrate = 3 d d`. A global clock must be used. This is the default or can be set with `clock = 1`. **COMMENTS** See [Date File](#date-file) for a complete description of the date file format. **EXAMPLES** ``` datefile = dates.txt datefile = /home/foo/seqDates.txt ``` ### 32 thetamodel ------------------------------------------------------------------------ ``` thetamodel = s ``` **DESCRIPTION** Specifies how theta parameters are linked across populations in the species tree. **VALUES** `linked-none`, no linking of theta parameters (default). `linked-all`, all populations share the same theta. `linked-inner`, inner (ancestral) populations share the same theta. `linked-msci`, linking appropriate for MSC-I models. `linked-mscm`, linking appropriate for MSC-M models. **DEFAULT** `linked-none` **EXAMPLES** ``` thetamodel = linked-none thetamodel = linked-all thetamodel = linked-inner ``` ### 33 printlocus ------------------------------------------------------------------------ ``` printlocus = +d +d +d ... ``` **DESCRIPTION** Specifies which loci should have their gene trees printed to the output. **VALUES** The first integer specifies the number of loci to print, followed by the locus numbers (1-indexed). **EXAMPLES** ``` printlocus = 2 1 5 # print gene trees for loci 1 and 5 printlocus = 3 1 2 3 # print gene trees for loci 1, 2, and 3 ``` ### 34 geneflow ------------------------------------------------------------------------ ``` geneflow = b ``` **DESCRIPTION** Enables or disables gene flow estimation in the model. **VALUES** `0`, disable gene flow estimation. `1`, enable gene flow estimation. **DEFAULT** `0` **EXAMPLES** ``` geneflow = 0 geneflow = 1 ``` ### 35 loadbalance ------------------------------------------------------------------------ ``` loadbalance = s ``` **DESCRIPTION** Specifies the load balancing strategy for distributing computation across threads. **VALUES** `none`, no dynamic load balancing (static distribution). `zigzag`, use zigzag load balancing for better distribution of work across threads. **DEFAULT** `none` **EXAMPLES** ``` loadbalance = none loadbalance = zigzag ``` ### 36 traitfile ------------------------------------------------------------------------ ``` traitfile = s ``` **DESCRIPTION** Specifies the file containing trait data for morphological analysis. **VALUES** `s`, a string specifying the path to the trait data file. **EXAMPLES** ``` traitfile = traits.txt traitfile = /home/user/data/morphology.txt ``` ## Example control file To examine the structure of a typical BPP control file, we consider the example file `A00.bpp.ctl` contained in the BPP distribution subdirectory `examples/frogs`. The goal here is to provide an example of a correctly formatted control file and briefly explain the syntax and meaning of specified variables, as well as proposing some conventions for writing optional comments in control files. An exhaustive description of all control file variables, their effects, and permissible values is provided in the definitions of [BPP Control File Variables](#bpp-control-file-variables). In depth descriptions of specific model and prior choices can be found in [Substitution Models](#substitution-models), [Introgression and Migration Models](#introgression-and-migration-models), and [Isolation with Migration Models](). Control files illustrating the 4 methods of analysis A00, A01, A10 and A11 are found in [Methods of Analysis](#methods-of-analysis). The content of the control file `A00.bpp.ctl` is displayed below: ``` seed = -1 seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out # fixed species tree and delimitation speciesdelimitation = 0 speciestree = 0 species&tree = 4 K C L H 9 7 14 2 (((K, C), L), H); # sequence data are unphased for all 4 populations phase = 1 1 1 1 # 0: no data (prior); 1:seq like usedata = 1 # number of data sets in seqfile nloci = 5 # remove sites with ambiguity data (1:yes, 0:no)? cleandata = 0 # invgamma(a, b) for theta thetaprior = invgamma 3 0.004 E # invgamma(a, b) for root tau & Dirichlet(a) for other tau's tauprior = invgamma 3 0.002 * heredity = 1 4 4 * locusrate = 1 5 # finetune: auto-adjust step lengths during burnin finetune = 1 # MCMC samples, locusrate, heredityscalars, Genetrees print = 1 0 0 0 burnin = 8000 sampfreq = 2 nsample = 100000 ``` **A comment about comments** Optional comments may be added to the control file. If either a '\#' or a '\*' symbol occurs on a line, everything following the symbol up to the end of the line will be ignored by the BPP program. This allows optional comments to be added by the user either to explain choices (or meanings) of variables, or to retain optional alternative variable settings for future use. The control files presented in this documentation use the different comment-initiating symbols for distinct purposes: 1. '\#' initializes a comment that lists and/or explains the possible settings of the variable found immediately below the comment. 2. '\*' initializes a comment that retains deactivated variable settings (for future use). If an alternative to a defined variable it follows that variable. It is good practice to make use of these conventions for comments in your own bpp control files. ### Control file variables explained Here we walk through the variables in the above control file example explaining the meaning of each and the valid values that can be specified. An exhaustive description of all control file variables and permissible values is given in [Control File](#control-file). ``` seed = -1 ``` This line specifies the seed value used to initialize the random number generator. Because BPP uses pseudorandom numbers to drive a stochastic algorithm for inference the results will vary between analyses initiated with different seeds, but will be identical if the same seed is used (because pseudorandom numbers are generated using a deterministic algorithm started with the value of seed). Setting `seed = -1` specifies that BPP use the clock on the computer to generate a seed. With this setting, every run of the program will automatically be started with a different seed. The seed that was used for a run is stored in a file generated by the program called SeedUsed. ``` seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out ``` These lines specify the names of the files used to run the program and summarize the output. The program will assume that all these files are located in the current working directory since no path is specified. The line `seqfile = frogs.txt` specifies that the sequence data are in the file named `frogs.txt`. See [Sequence File](#sequence-file) for a description of the sequence data file format. The line `Imapfile = frogs.Imap.txt` specifies that the map file is named `frogs.Imap.txt`. See [Imap File](#imap-file) for a description of the map file format. The line `jobname = out` specifies that the program output and mcmc samples are printed to files named `out.txt` and `out.mcmc.txt`, respectively. The format of the output file depends on the type of analysis being performed and is described in [Methods of Analysis](#methods-of-analysis). ``` # fixed species tree speciesdelimitation = 0 speciestree = 0 ``` The line `speciesdelimitation = 0 ` specifies that the number of species is fixed (not delimited) and the line `speciestree = 0` specifies that the species phylogeny (topology) is fixed (not estimated). This combination specifies an analysis of type A00, which is parameter estimation under a MSC using multiple-loci and multiple-species data on a fixed species tree ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645); [Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148)). Analysis A00 is described in detail in [A00: Demography and Divergence Time Estimation](#a00-demography-and-divergence-time-estimation). ``` species&tree = 4 K C L H 9 7 14 2 (((K, C), L), H); ``` The line `species&tree = 4 K C L H` specifies the number of species/populations and their names. This example specifies 4 species in the data, which are K, C, L and H. The next line ` 9 7 14 2` specifies the maximum number of sequences at any locus for each species/population (in the same order as the species listing on the previous line). This example specifies a maximum of 9 sequences for K, 7 for C, 14 for L and 2 for H. The program interprets these numbers as follows. First, if the number of sequences for a species is 2 or greater, $\theta$ for that species will be a parameter estimated by the program. Second, the sum of the numbers of sequences over all species specifies the maximum number of sequences at a locus. However, specifying 2 sequences for a species when there are actually 10 sequences for that species at some loci in the data file has no harmful effects. Note that if there are $s$ species in the species tree, the model will minimally contain the following parameters: $s – 1$ species divergence times ($\tau$s) and $s – 1$ ancestral species $\theta$s. A $\theta$ parameter will also exist for each contemporary species with at least two sequences present at some locus. If the data are for a single species/population, the model will contain one parameter only, $\theta$ for that species. The line `(((K, C), L), H);` specifies the (fixed) species tree in Newick format ending with a semicolon (;). This example specifies a tree with no introgression (the MSC model). Specification of the MSC-I model with introgression is different (see [Introgression and Migration Models](#introgression-and-migration-models) for a description of the MSC-I model). The Newick and extended Newick formats for specifying phylogenetic trees are described in [Introgression and Migration Models](#introgression-and-migration-models) and a good overview is found at . ``` # sequence data are unphased for all 4 populations phase = 1 1 1 1 ``` The `phase` variable is specified using a boolean (0/1) indicator variable for each species/population, with 0 (the default) specifying that sequences from that species are fully phased haplotypes and 1 specifying they are unphased diploid sequences. If this line is omitted, all sequences are assumed to be fully phased (that is, value 0 for all species). The program does not allow some sequences from a species to be phased while other sequences from that same species are unphased. If the variable is 1, each sequence from that species is treated as an unphased diploid sequence, with heterozygous sites represented using ambiguity characters `YRMKSW`. The chararcters `N`, `-` and `?` are allowed in the sequence file and are all treated as `?`, but other ambiguities (HBVD) are not allowed. Each sequence is then expanded/resolved into two sequences, and the program averages over all the possible phase resolutions according to the approach of [Gronau et al 2011](https://doi.org/10.1038/ng.937). For example a sequence with two heterozygous sites R and Y is resolved into two possible haplotypes: (i) A\...T and G\...C, and (ii) A\...C and G\...T. Note that this option forces a sequence with no heterozygote sites to be phased/resolved into two identical sequences. The MSC model averages over different phase resolutions of the heterozygote sites in the likelihood calculation. This option may increase the memory usage and CPU time considerably if there are many sequences with many heterozygote sites at one locus. The simulation program MCcoal has been modified to simulate diploid unphased data as well. This allows one to use simulations to examine the difference between the analysis of the full data and the unphased data. If the phase flag for a species is 0, all sequences from that species are treated as resolved haplotype sequences, and ambiguities are interpreted in the usual way (for example, Y means one nucleotide that is either T or C). ``` # 0: no data (prior); 1:seq like usedata = 1 ``` The line `usedata = 1` specifies that the likelihood is used in the MCMC run to generate the posterior distribution of parameters. Setting this to 0 allows one to run the MCMC without sequence data to generate the prior (mostly used for debugging). ``` # number of data sets in seqfile nloci = 5 ``` The line `nloci = 5` specifies the number of loci (alignments) to be 5. There may be more loci in the sequence data file than is specified here. For example, if 200 loci exist in the sequence data file and you specify `nloci = 2` in the control file, BPP will read the first two loci only. ``` # remove sites with ambiguity data (1:yes, 0:no)? cleandata = 0 ``` The line `cleandata = 0` specifies that columns in the alignment which have gaps or ambiguity characters are retained and used in the likelihood calculation. Setting `cleandata = 1` instead specifies that such data are removed prior to analysis. ``` # invgamma(a, b) for theta thetaprior = invgamma 3 0.004 E ``` The line `thetaprior = 3 0.004 E` specifies the inverse-gamma prior IG($\alpha, \beta$) for the $\theta$ parameters, with the mean to be $\beta/(\alpha-1)$. In the example, the mean is $0.004/(3 – 1) = 0.002$ (two substitutions per kb on average) and the option E specifies that the program analytically integrates over $\theta$s. Note that all $\theta$ parameters in the MSC model (for both modern species and extinct ancestral species) are assigned the inverse-gamma prior with the same parameters. **Analytical integration of $\theta$ and the inverse gamma prior** The inverse-gamma is a conjugate prior for $\theta$ ([Hey and Nielsen 2007](https://doi.org/10.1073/pnas.0611164104)) which means that both the prior and the posterior of $\theta$ will be inverse-gamma. Use of the conjugate priors allows the $\theta$ parameters to be integrated out analytically, and thus the dimension of the parameter space is reduced. This typically leads to improved mixing of the MCMC. However, with this option, the posterior of the $\theta$ parameters will not be produced. To estimate the $\theta$ parameters, add the letter e (or E) on the line, as follows ``` thetaprior = invgamma 3 0.002 e ``` Whether $\theta$s are integrated out or estimated, other results (such as the posterior of the $\theta$ parameters or the posterior probability of species trees or species-delimitation models) should be identical. A useful strategy may be to run the A01, A10, and A11 analyses without estimating the $\theta$ parameters, and after the MAP model (the best species tree or the best species delimitation model) is identified, run the A00 analysis with the model fixed to estimate all parameters. Some notes about the inverse-gamma distribution. Since BPP3.4, both the $\theta$ and $\tau$ parameters are assigned the inverse gamma priors rather than the gamma priors in version 3.3 or earlier. One difference is that the gamma is light-tailed while the inverse-gamma is heavy-tailed, so that the inverse-gamma may be less influential than the gamma if your prior mean is much too small. The inverse-gamma distribution IG($\alpha, \beta$) has mean $m = \alpha/(\beta-1)$ if $\alpha > 1$ and variance $s^2 = \beta^2/[(\alpha-1)^2(\alpha-2)]$ if $\alpha > 2$, while the coefficient of variation is $s/m = \sqrt{1/(\alpha-2)}$. If little information is available about the parameters, you can use $\alpha = 3$ for a diffuse prior and then adjust the $\beta$ so that the mean looks reasonable. For example, for the human species, $\theta_H \approx 0.0006$, which means that two random sequences from the human population are different at  0.06% of sites, less than 1 difference per kb. A sensible diffuse prior is then `thetaprior = 3 0.002`, with mean $0.002/(3 – 1) = 0.001$. ``` # invgamma(a, b) for root tau & Dirichlet(a) for other tau's tauprior = invgamma 3 0.002 ``` The line `tauprior = 3 0.002` specifies the inverse-gamma prior IG($\alpha, \beta$) for $\tau_0$, the divergence time parameter for the root in the species tree. Other divergence times are generated from the uniform Dirichlet distribution ([Yang and Rannala 2010, eq.2](https://doi.org/10.1073/pnas.0913022107)). In the example, the mean is $0.002/(3 – 1) = 0.001$ (which means 0.1% average sequence divergence between the root of the species tree and the present time). If the mutation rate is $10^{–9}$ mutations/site/year, this distance will translate to a divergence time of 1 MY. ``` # finetune: auto-adjust step lengths during burnin finetune = 1 ``` The line `finetune = 1` enables automatic optimization of step lengths during burnin. This is the recommended setting. The program will adjust step lengths to achieve good mixing. If auto-tuning is enabled, burnin must be > 200. For custom step lengths, use the key:value format (see section 5). ``` # MCMC samples, locusrate, heredityscalars, Genetrees print = 1 0 0 0 burnin = 8000 sampfreq = 2 nsample = 100000 ``` The line `print = 1 0 0 0` specifies the information that is printed out during the MCMC run. It has several boolean variables (0 for off and 1 for on) that control the printouts during the MCMC. The first variable specifies whether the MCMC samples are written into the main sample file (out.mcmc.txt). The second flag is for locus-specific rate parameters (rate $\mu_i$, variance parameter $\nu_i$ and species-tree branch rates for locus $r_{ij}$). The third flag is for locus-specific heredity scalars if those are estimated from the data. The fourth flag is for printing locus-specific gene trees. The fifth flag is for printing locus-specific substitution-rate parameters (qmat for the $Q$ matrix, freqs for base frequencies, and alpha for gamma rates for sites). MCMC samples are taken after the burnin (8000 iterations), and in this example, are taken every 2 iterations, with a total of 100,000 samples taken. The total number of MCMC iterations is `burnin` + `sampfreq` $\times$ `nsample`. The resulting file can be large. For analysis A00 (`speciesdelimitation = 0, speciestree = 0`), this file can be analysed using the R package or tracer. For other analyses (A01, A10, and A11), the dimension of the sample is changing and cannot be analysed using such traceplot programs. **Combining results from different runs** The variable setting `print = -1` specifies that BPP will not perform an MCMC analysis and will instead read an existing MCMC sample file and summarize the results. Thus, with setting `print = 1`, the `out.mcmc.txt` file receives the MCMC program output, but with `print = -1`, it instead provides the program input. When using this option be careful not to overwrite files that you intended to keep. The `print = -1` option is useful for combining MCMC samples from multiple runs to produce the posterior summary. Suppose you run the same analysis 3 times in different directories -- you can merge the out.mcmc.txt files from the three runs into one file and run the program using the `print = -1` option to summarize the posterior for the combined sample. Note that if you forced one or more of the jobs to terminate (using the Unix kill command for example) the last line of each out.mcmc.txt file may incomplete. If that is the case, you should delete the incomplete line. The header lines of any files being concatenated should also be deleted. Do not combine the MCMC samples from different types of analyses (i.e., using different data files and/or variable settings or priors). ## Substitution Models BPP includes several models that allow more realistic representations of the processes underlying the evolution of DNA and amino acid sequences. Fundamental to all molecular phylogeographic analyses is the model of the nucleotide or amino acid substitution process itself. DNA substitution models can accomodate different rates of substitution between nucleotides due to factors such as transition versus transversion bias and unequal stationary nucleotide frequencies. BPP allows the use of many different DNA base substitution and amino acid substitution models. For DNA sequences, the models range from a simple Jukes-Cantor model (computationally efficient) to a sophisticated General Time-Reversible (GTR) model (computationally expensive but more realistic). For closely related species the Jukes-Cantor model is often adequate, while the GTR often performs better for more distantly related species. The remaining DNA substitution models (K80, F81, HKY, T92, TN93, and F84) are all special cases (submodels) of the GTR model and are of intermediate complexity, usually having one or two free parameters ([see Yang 2014](https://global.oup.com/academic/product/molecular-evolution-9780199602605?lang=es&cc=gb)). Many fixed substitution matrix models (no free parameters) for modeling substitutions among amino acid sequences are also available. Another factor to consider is variation in overall (or average) substitution rates. The overall rate of substitution can vary in multiple ways and, based on empirical evidence, this type of variation is one of the most important factors to accommodate during phylogeographic inference. There are 3 distinct types of substitution rate variation that can be accommodated using the models implemented in BPP: rate variation among sites, rate variation among loci, and rate variation among species (or populations). Substitution rate variation among sites, or among loci, is a common pattern and is likely due to the pervasive effects of negative (purifying) selection. Variation of substitution rates among species (or populations) is more often observed among distantly related species and may reflect the effects of differential selection, evolution of DNA repair mechanisms, or other factors. This Chapter briefly describes the models of DNA substitution, and of substitution rate variation, implemented in BPP, the parameters and priors associated with these models, and how to specify the different models in the BPP control file. ### DNA Substitution Models #### Jukes-Cantor model One of the earliest and simplest models of DNA substitution, the Jukes-Cantor model, assumes that all possible changes between the 4 nucleotides have equal rates and thus the stationary frequencies of bases are $1/4$ each. This model is the default in BPP and will be used if no model is specified in the control file. The control file variable `model` specifies the substitution model and the explicit specification of the JC69 model is: ``` model = JC69 ``` The JC69 model appears exceedingly unrealistic. However, for closely related sequences (with few multiple substitutions) it often performs as well as more sophisticated models while incurring much less computational expense. There are no parameters associated with the JC69 model (other than the substitution rate) so no additional control file variables associated with the JC69 substitution model need to be specified. #### GTR model The General Time-Reversible (GTR) model implemented in BPP allows different rates of substitution between different nucleotide pairs, with the constraint that the process is "reversible" so that, for example, the rate of transition from A to T equals the rate of transition from T to A and so on. The GTR model accommodates factors such as transition-transversion bias and should be used when sequences are more divergent (distantly related species). The control file specification of the GTR model is: ``` model = GTR ``` Note that if the GTR model is used, two other variables may be specified in the control file: `Qrates` and `basefreqs`. The `Qrates` variable specifies whether the exchangeability parameters of the GTR model are fixed or variable as well as either their fixed values, or prior mean values, respectively. The control file specification is: ``` Qrates = B a b c d e f ``` The variable `B` is either 0 (variable exchangeability parameters) or 1 (fixed exchangeability parameters). The variables `a b c d e f` are positive numbers specifying the exchangeability parameters for TC, TA, TG, CA, CG and AG substitutions respectively (see [Yang 1994](https://doi.org/10.1007/BF00178256)). Note that TC is the rate of either T-\>C or C-\>T substitutions. Note that only the relative rates matter because the rate matrix is scaled to have an overall average substitution rate of 1. Thus, for the model of fixed rates if x is any positive number a b c d e f and xa xb xc xd xe xf are equivalent specifications. An example with fixed rates is as follows: ``` Qrates = 1 2 1 1 1 1 2 ``` which specifies $a=2,b=c=d=e=1,f=2$ so that the transition rate is twice as high as the transversion rate, corresponding to a K80 or HKY model. In most cases, users will want to use the GTR model with exchangeability parameters variable. Using this specification, the posterior distribution of the parameters will be estimated from the data. An example with variable rates is as follows: ``` Qrates = 0 10 5 5 5 5 10 ``` The absolute magnitude of the numbers *does matter* in this case as it determines the variance of the prior probability distribution on exchangeability parameters. The variance is inversely proportional to the absolute magnitude of the variables. The prior on the free parameters $a, b, c, d, e$ is a Dirichlet distribution with parameters $\alpha_a = 10, \alpha_b = \alpha_c = \alpha_d = \alpha_e = 5, \texttt{and} \, \alpha_f = 10$ The rates generated from the Dirichlet sum to 1 and are rescaled. If we define $\alpha = \alpha_a + \alpha_b + \alpha_c + \alpha_d + \alpha_e + \alpha_f$ then the mean for $a$ is $\alpha_a/\alpha$ and the variance for $a$ is $\textrm{Var}(a) = \frac{\frac{\alpha_a}{\alpha}\left(1-\frac{\alpha_a}{\alpha}\right)}{\alpha + 1}$. The means and variances of the other parameters are the same (substituting subscripts). Thus, in this example, the mean for $a$ is $\alpha_a/\alpha = 10/40 = 1/4$ and the variance is $(1/4)(3/4)/41 = 4.6 \times 10^{-3}$. The prior specified above has an average transition-transversion rate ratio of 2 (if base frequencies are all fixed at 1/4). The other control file variable that must be specified when using the GTR model is the `basefreqs` variable which specifies whether the stationary nucleotide (base) frequencies are fixed, or variable, and either their fixed values, or relative prior means, respectively. The control file specification is: ``` basefreqs = B T C A G ``` where B is either 0 (variable base frequencies) or 1 (fixed base frequencies) and T C A and G are numbers specifying the frequencies of bases T, C, A and G, respectively. An example with fixed base frequencies is as follows: ``` basefreqs = 1 0.1 0.2 0.2 0.5 ``` In most cases, users will want to use the variable base frequencies setting which specifies that the posterior distribution of base frequencies is estimated from the data. An example with variable base frequencies is as follows: ``` 0 10 10 10 10 ``` Here the 4 variables T C A G specify the alpha parameters of a Dirichlet prior distribution on the nucleotide frequencies. So, the prior mean frequency of base A in this example is $10/40 = 1/4$ and the variance is $(1/4)(3/4)/41 = 4.6 \times 10^{-3}$. As noted above when discussing the Dirichlet prior for exchangeability parameters the absolute magnitude of the variables determines the variance of the Dirichlet prior. If a GTR model is chosen and `Qrates` and `basefreqs` are not specified in the control file default values of `Qrates = 0 1 2 1 1 2 1` and `basefreqs = 0 1 1 1 1` will be used. #### Other DNA substitution models A range of models are available that fall between the JC69 model and the GTR model in terms of complexity and number of parameters. The additional models available include: K80, F81, HKY, T92, TN93, and F84. The priors on the parameters of these models are pre-specified by the program. #### Amino acid substitution models There are many amino acid substitution models available. These models completely specify the substitution matrix so there are no free parameters. Many of the models are described in [Yang 2014](https://global.oup.com/academic/product/molecular-evolution-9780199602605?lang=es&cc=gb). The available models are: DAYHOFF, LG, DCMUT, JTT, MTREV, WAG, RTREV, CPREV, VT, BLOSUM62, MTMAM, MTART, MTZOA, PMB, HIVB, HIV,W, JTTDCMUT, FLU, and STMTREV. ### Custom among-locus model variation If a single substitution model is specified, all loci are assigned that model. The control file option `model = Custom filename` can be used to assign different models among loci, where the models are specified in the file `filename`. This file should contain one or more lines with the format: loci_indices data_type model where `loci_indices` is an index or range of indexes for loci (indexed according to their order of occurrence in the sequence data file), `data_type` is either `DNA` or `AA`, and model specifies the substitution model variable for the indexed loci as described above. For example, the entry 1-10 DNA JC69 11 AA DAYHOFF specifies that loci 1 to 10 comprise DNA sequences that will have the `JC69` model applied to them, while locus 11 is an amino acid sequence that will have the `DAYHOFF` model applied to it. Models of Substitution Rate Variation ------------------------------------- ### Among-site rate variation Substitution rates vary among sites due to many factors, including purifying selection in coding regions. BPP allows among-site rate variation to be accomodated using the canonical Gamma distributed rate variation model ([see Yang 1994](https://doi.org/10.1007/BF00178256)). In this model, the mean rate among sites is constrained to be 1 so that the expected number of substitutions on a branch is proportional to its length (in units of expected substitutions). If the mean of the Gamma distribution is fixed, the only free parameter is the shape parameter $\alpha$. The variance of rates among sites is $1/\alpha$ and therefore a smaller $\alpha$ specifies more rate variation. Gamma distributed rate variation is specified using the control file variable `alphaprior`. If this variable is not set the default is $\alpha = \infty$ so that the variance is zero (rates are constant among sites). The prior on $\alpha$ is also a Gamma distribution with parameters $\alpha_p$ and $\beta_p$. The control file specification is as follows: alphaprior = a b ncatG where $a = \alpha_p$ and $b = \beta_p$ are the parameters of the Gamma distributed prior on the parameter $\alpha$ and `ncatG` is the number of rate categories used to approximate the Gamma distribution. A discrete distribution is used to efficiently approximate the Gamma distribution and ncatG is the number of rate categories used in the approximation -- more rate categories provide a better approximation but incur greater computational expense. A value for ncatG of 4 is recommended. The variance and mean of the prior on $\alpha$ are: $\begin{aligned} \textrm{Mean}(\alpha) & = & \frac{\alpha_p}{\beta_p}, \nonumber \\ \textrm{Var}(\alpha) & = & \frac{\alpha_p}{\beta_p^2} \nonumber\end{aligned}$ As an example, the setting `alphaprior = 1 1 4` specifies a prior on the shape parameter $\alpha$ (with 4 rate categories) that has a mean and variance of 1, thus the expected variance of $\alpha$ is also 1. This is a fairly diffuse (uninformative) prior. ### Among-locus rate variation Empirical evidence suggests that substitution rates can vary substantially among loci, probably reflecting the influences of selection. BPP incorporates a hierarchical prior to model rate variation among loci. The models implemented in BPP allow each locus to have a different substitution rate. For $L$ loci (alignments), there are $L$ substitution rate parameters, where $\mu_i$ is the rate at the $i$th locus. #### Prior distribution of the mean rate among loci The mean rate across loci is $\bar{\mu} = \frac{1}{L} \sum_{i=1}^L \mu_i$. The prior distribution for the mean rate across loci $\bar{\mu}$ is a Gamma distribution with parameters $\alpha_{\bar{\mu}}$ and $\beta_{\bar{\mu}}$ and mean and variance, $\begin{aligned} \textrm{Mean}(\bar{\mu}) & = & \frac{\alpha_{\bar{\mu}}}{\beta_{\bar{\mu}}}, \nonumber \\ \textrm{Var}(\bar{\mu}) & = & \frac{\alpha_{\bar{\mu}}}{\beta_{\bar{\mu}}^2}. \nonumber\end{aligned}$ If there are no fossil calibrations in the species tree, the rates should all be relative. In this case, to avoid overparameterization we suggest fixing the mean substitution rate across loci to be $\bar\mu = 1$. This is done by setting $\alpha_{\bar{\mu}} = \beta_{\bar{\mu}} = \infty$ (see discussion of control file settings below), or setting both these parameters to large values such as $\alpha_{\bar{\mu}} = \beta_{\bar{\mu}} = 100$, in which case $\bar{\mu}$ is nearly fixed. #### Prior distribution for locus-specific rates given mean rate among loci Conditioning on the mean rate $\bar\mu$, two priors are available to specify the probability distribution of locus-specific rates, $\mu_i$. The control file variable argument `prior` (see discussion of control file settings below) can take two values `dir` (for Gamma-Dirichlet rates for loci), and `iid` (for conditional i.i.d. or hierarchical prior). If `prior = dir` (Gamma-Dirichlet distribution of locus rates), the total rate $L\bar\mu$, given the mean rate ($\bar\mu$), is partitioned into locus rate variables , $\mu_i$, that are constrained to sum to $L\bar\mu$ and have a Dirichlet distribution with concentration parameter $\alpha_\mu$. The model and notation follow [Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148) eq.4 and [Dos Reis et al 2014](https://doi.org/10.1093/sysbio/syu020) eqs.3-5. If `prior = iid` (Conditional-i.i.d. or hierarchical prior of locus rates), the locus rates, $\mu_i$, are independent and identically distributed i.i.d. given the mean rate $\bar\mu$, so the prior distribution of $\mu_i$ is: $\mathtt{P}(\mu_i) = \textrm{Gamma}(\alpha_{\mu_i}, \alpha_{\mu_i/\bar{\mu}})$. This model is described in [Zhu et al 2015](https://doi.org/10.1093/sysbio/syu109) eq.8 and implemented in [mcmctree]{.smallcaps}. In both the Gamma-Dirichlet and the Conditional i.i.d. models, the parameter $\alpha_{\mu_i}$ is inversely related to the extent of rate variation among loci, with a large $\alpha_{\mu_i}$ meaning similar rates among loci (less rate variation). If all loci are noncoding, the rates are probably similar, so $\alpha_{\mu_I} =$ 10 or 20 may be reasonable, while for coding loci or exons, $\alpha_{\mu_i} = 2$ or 1 may be appropriate. The $\alpha_{\mu_i}$ parameter may affect the estimates of the population size parameter ($\theta$) for the root node on the species tree ([Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148)). The $L$ locus rates ($\mu_i$) are parameters in the model. If $\alpha_{\bar{\mu}} > 0$, the mean rate, $\bar\mu$, is a parameter as well. #### Control file settings for locusrate The `locusrate` variable takes from 1 and 5 arguments. In the case of 5 arguments the form is shown below, with the last argument (prior) being optional. locusrate = rate_variation a_mubar b_mubar a_mui prior The rate\_variation variable is either 0 (no rate variation, the default), 1 (rate variation with rates inferred) or 2 (rate variation with fixed user-specified rates). If rate\_variation = 0 no other variables are specified, if rate\_variation = 2 then a second argument specifies the name of a file containing $L$ specified fixed rates, one for each locus. If rate\_variation = 1 the additonal variables `a_mubar` and `b_mubar` specify the shape and rate parameters ($\alpha_{\bar{\mu}}$ and $\beta_{\bar{\mu}}$) of the gamma distribution that is the prior for the mean rate across loci (`mubar` or $\bar\mu$), and variable `a_mui` ($\alpha_{\mu_i}$) and `prior` are used to specify parameter of the prior on the locus rates ($\mu_i$) given the mean rate ($\bar\mu$) and the form of that prior (`dir` or `iid`). Some example configuration file settings are: # (0: No rate variation among loci. This is default) locusrate = 0 # (1: estimate locus rates mui) locusrate = 1 10 10 5.0 iid # (1: estimate locus rates mui) locusrate = 1 0 0 5.0 iid # (2: locus rates from file) locusrate = 2 LocusRateFileName Note that `locusrate = 2 LocusRateFileName` specifies the fixed-rates model of locus-rate variation ([Burgess and Yang 2008](https://doi.org/10.1093/molbev/msn148)). This is the strategy used by [Yang 2002](https://doi.org/10.1093/genetics/162.4.1811), with the relative rates estimated by the distance to an outgroup species. The relative locus rates are listed in the file: there should be as many numbers in the file, separately by spaces or line returns, as the number of loci (`nloci`). The program re-scales those rates so that the average across all loci is 1 and then use those relative rates as fixed constants. The following example nearly fixes (first line), or fixes (second line) the mean rate $\bar{\mu}$. This is the recommended practice when no fossil calibrations are available: # locusrate = 1 a_mubar b_mubar a_mui # (1: estimate locus rates mui) locusrate = 1 10 10 5 iid # (1: estimate locus rates mui, mubar=1 fixed) locusrate = 1 0 0 5 iid Note that `a_mubar = 0` and `b_mubar = 0` (with 0 meaning $\infty$) cause the program to fix $\bar{\mu}$, while equal and large values such as `a_mubar = 100` and `b_mubar = 100` cause cause $\bar\mu$ to be nearly fixed at 1. **The locus-rate prior in IMa.** The model of variable rates among loci implemented here has some differences from a similar model implemented in the IMa program [Hey and Nielsen 2004](https://doi.org/10.1534/genetics.103.024182). The biggest difference appears to be the parameterization. BPP defines mutation rate on a per-nucleotide basis, so the prior specifies that the expectation of the mutation rate per site is constant among loci. IMa defines mutation rate on a per-locus basis, so its prior specifies that the expectation of the mutation rate per locus is the same among loci. If locus one has 100 sites and locus two has 1000 sites, then IMa assumes that the per-site rate for locus one is 10 times for locus two, while BPP assumes the same per-site rate. Also IMa constrains the geometric mean of rates across loci to be one, while BPP constrains their arithmetic mean to be one. ### Among-species rate variation The strict molecular clock model of evolution assumes that rates of substitution for a particular gene, or genomic region, are the same in different species. Early applications of relative rates tests suggested that this assumption is often violated in empirical data, probably due to factors such as selection and in more distantly related species evolving DNA repair and proof-reading mechanisms. One can allow rates of substitution to vary in an arbitrary way be inferring unrooted phylogenetic trees but this does not allow species divergence times to be estimated. Relaxed clock models allow the rate of substitution to vary among contemporary and ancestral species but with some constraints which allow species divergence times to be estimated. The BPP program implements two relaxed clock models, as well as a strict molecular clock model. The two relaxed clock models differ in terms of their assumptions about rates in ancestor and descendent species. The autocorrelated rates model assumes that descendents have rates that are correlated with the rate in the ancestor, while independent rates model assumes that the rates on each branch are iid from a common distribution. For the independents rates model, the variance of the distribution of rates among branches determines the degree of departure from a molecular clock - a low variance means that the rates are almost clocklike (vary little among branches). #### Variation in the degree of departure from a molecular clock among loci For a given set of species some loci may have rates that are highly clock-like while other loci may have much greater variation in substitution rates among species (branches of the species tree). To accomodate this, the variance of the substitution rate among species is allowed to vary among loci. The prior distribution of the average variance $\bar{\nu}$ for all loci is $P(\bar{\nu}) = \textrm{Gamma}(\alpha_{\bar{\nu}},\beta_{\bar{\nu}})$, so that the expected average variance is $\alpha_{\bar{\nu}}/\beta_{\bar{\nu}}$. If this value is larger loci depart more from the molecular clock on average. Conditional on the average variance parameter $\bar{\nu}$ the prior for the variance $\nu_i$ at locus $i$ can be either a Gamma-Dirichlet distribution or a Conditional i.i.d distribution (controlfile variables `dir` and `iid`, respectively). Both distributions have a parameter $\alpha_{\nu_i}$ that determines the variance of the variance in rates across loci. The specification of the distribution of $\nu_i$ given the mean $\bar\nu$ follows the same procedure as the specification of the distribution of the locus-specific substitution rate $\mu_i$ given the mean rate $\bar\mu$ in the previous section; see notes above about the locusrate variable. Note that $\alpha$ (control file variable `a_vi`) is inversely related to the variance in $\nu_i$: use small values of $\alpha$ (2, 1, or 0.5) if you believe that $\nu_i$ varies among loci (meaning that the clock nearly holds at some loci but is seriously violated at others). #### Relaxed clock with i.i.d rates among branches The independent rates model of the relaxed molecular clock assumes that for locus $i$ each branch has a rate that is independent and identically distributed (i.i.d). The distribution has mean $\mu_i$ specified by the `locusrate` model and variance $\nu_i$ specified by the prior distribution for rates among loci discussed above. Two distributions are available: a Gamma distribution or a log-normal distribution (control file variables `LN` and `G`, respectively). Given the locus rate $\mu_i$ (specified using the `locusrate` variable) and the variance $\nu_i$ (specified using the `clock` variable) for locus $i$, the different lineages may have different rates at the locus, and those rates are independent among loci. If distribution = G (for gamma), the rate for species-tree branch $j$ at locus $i$ has the following gamma distribution $r_{ij} | \mu_i, \nu_i \sim G(\mu_i^2/\nu_i, \mu_i/\nu_i)$. This has mean $\mu_i$ and variance $\nu_i$. Alternatively if distribution = `LN` (for log-normal), the rate for species-tree branch $j$ at locus $i$ has the following the log-normal distribution $$r_{ij} | \mu_i, \nu_i \sim \mathrm{LN}(\mu_i, \nu_i),$$ where $\mu_i$ is the mean of the LN distribution and $\nu_i$ is the variance parameter of the lognormal. #### Relaxed clock with autocorrelated rates among branches The autocorrelated-rates model (control file option `clock = 3`) with the log-normal distribution (LN) specifies the geometric Brownian motion model of [Rannala and Yang 2007](https://doi.org/10.1080/10635150701420643). This assigns a rate to each species-tree branch, that is, to the mid-point of the branch. Given the rate at the species-tree root ($\mu_i$ at locus $i$), the rates for the two branches around the root are specified. Then given the rate for each ancestral branch, the rates for its two daughter branches are specified, by integrating over the rate at the internal node that is ancestral to the daughter branches. See figure 1 and equations 3-8 in [Rannala and Yang 2007](https://doi.org/10.1080/10635150701420643). The rates for all species-tree branches are thus assigned through a pre-order tree traversal, starting from the root moving to the tips, until all branches are visited. If the gamma distribution (G) is instead used the model works as follows, using a similar pre-order tree traversal. First the two branches at the species-tree root have the gamma distribution with mean $\mu_i$ (`mu_i`) and variance $\nu_i$ (`nu_i`). Then given the rate for each ancestral branch, the rates for its two daughter branches are specified as independent gamma variables with the mean to be the rate of the parental branch and with the variance to be $\nu_i$. The autocorrelated relaxed clock model is only implemented for the MSC model and is not available when using the MSC-I model. **Rates vary among species tree branches**\ The rate-drift model specifies rates for branches on the species tree (rather than on the gene tree) for each locus, and gene-tree branches residing in the same population or species have the same rate. For example, if all sequences at a locus are from the same species and all coalescent events occur in that species (before reaching an ancestral species), all branches on the gene tree will have the same rate even if the relaxed-clock model allows rates to vary among species. In contrast, if a gene-tree branch passes several species or populations, the different segments of the branch will have different rates. In that case, the branch length on the gene tree is the sum of the lengths of those segments (with the length of each segment being the product of the rate and the time duration for the segment). #### Control file settings for clock The `clock` variable takes 6 arguments when using a relaxed clock as follows: ``` clock = clock_type a_vbar b_vbar a_vi prior dist ``` the argument `clock_type` takes on of three values: 1 (strict clock, default), 2 (relaxed clock with i.i.d rates) and 3 (relaxed clock with autocorrelated rates). In the case of the strict clock there are no additional arguments. The next two variables `a_vbar` and `b_vbar` specify the parameters ($\alpha_{\bar{\nu}}$ and $\beta_{\bar{\nu}}$) of the Gamma prior specifying the average rate variance across loci (see above). The variable `a_vi` specifies the variance parameter of the prior on the locus specific variance, $\alpha_{\nu_i}$ (see above). The variable `prior` is either `iid` or `dir` and specifies the prior on $\nu_i$ given $\bar{\nu}$ (see above). The variable `dist` is either `LN` or `G` specifying either a Gamma distribution of a log-normal distribution for the branch rates conditional on the locus rate $\mu_i$ and among-branch variance $\nu_i$. Several example control file entries for the `clock` variables are provided below: ``` # (1: strict clock, default) clock = 1 # (2: independent-rates) clock = 2 10.0 100.0 5.0 iid G # (3: correlated-rates) clock = 3 10.0 100.0 5.0 iid G ``` The specification (`clock = 2 10.0 100.0 5.0 iid G`) means the following. First $\bar\nu$ is assigned a gamma distribution $G(10.0, 100.0)$, with mean 0.1. Given $\bar\nu$, the conditional i.i.d. prior means that $\nu_i \sim G(\alpha, \alpha/\bar\nu)$ with shape parameter $\alpha = 5.0$ (`a_vi`) and mean $\bar\nu$, for $i = 1, 2, \cdots, L$. As another example, the specification (`clock = 2 10.0 100.0 5.0 dir LN`) means that $\bar\nu$ is assigned a gamma prior $G(10.0, 100.0)$, and then the sum $L\bar\nu$ is partitioned into $\nu_i$ (for $i = 1, 2, \cdots, L$) according to the Dirichlet distribution (`prior = dir`) with concentration parameter $\alpha = 5.0$ (`a_vi`). Again large $\alpha$ means the same extent of clock violation at different loci. **Choosing $\alpha_{\bar{\nu}}$ and $\beta_{\bar{\nu}}$** A larger variance $\nu_i$ represents a greater violation of the molecular clock at locus $i$. Note that $\nu_i$ will be similar to $\bar\nu$, especially if $\alpha_{\bar{\nu}}$ (`a_vi`) is large, and $\bar\nu$ has prior mean $\alpha_{\bar{\nu}}/\beta_{\bar{\nu}}$. If (`a_vbar = 10 b_vbar = 100`), the prior mean will be 0.1. For the log-normal model $\nu = 0.5$ suggests a serious violation of the clock while $\nu < 0.1$ represents only a slight violation. ## Introgression and Migration Models Another important factor influencing the results of species tree inference is introgression (gene flow) between species (or populations). BPP also implements a model of episodic introgression between species or populations [Flouri et al 2020](https://doi.org/10.1093/molbev/msz296) as well as a model of ongoing gene flow (MSC-M). Since version 4.1, BPP implements the MSC-I model [Flouri et al 2020](https://doi.org/10.1093/molbev/msz296) and since version 4.4 it implements the MSC-M model. Currently, introgression (MSC-I) and ongoing gene flow (MSC-M) analyses can only be performed using the A00 method of analysis (fixed species tree and delimitation, see [Methods of Analysis](#methods-of-analysis)). In other words, the user has to specify the number of introgression events (or bands of gene flow), their directions, and the populations involved on a fixed species tree. The program will then estimate the parameters of the MSC-I (or MSC-M) models using MCMC. We hope to implement MCMC moves that change the MSC-I (or MSC-M) models and species tree topology in the future. ### The MSC-I Model The MSC-I model involves three types of parameters: the species divergence or hybridization times ($\tau$s), the population size parameters ($\theta$s), and the introgression probabilities ($\varphi$s). The control file must specify an extended Newick format tree with one or more introgression events as well as a prior distribution for the introgression probability $\varphi$. The MSC-I model is specified in the `species&tree` block using the extended Newick notation ([Cardona et al 2008](https://doi.org/10.1186/1471-2105-9-532)) ``` species&tree = 3 A B C 10 10 10 ((A, (C)H[&phi=0.5,&tau-parent=yes])S, (H[&tau-parent=yes], B) T)R; ``` The rules for the extended Newick format are as follows. '(A, B)S' specifies two branches from S to A and from S to B, while '(A)H' specifies one branch from H to A. Every branch in the tree model is represented once. Each tip species occurs once. Internal nodes for speciation nodes may and may not be labeled, but hybridization nodes must be labeled. The extended Newick notation is not unique. Each hybrid node has two parental species. Models A-C are distinguished using a meta-data variable `tau-parent`. The value yes means that the parental population exists and has a separate age parameter $\tau$, with an associated $\varphi$ parameter, while the value no means that such parameters do not exist. For example, in model (A), $\tau_S$, $\tau_T$, $\theta_{Hl}$ (for H-left), and $\theta_{Hr}$ (for H-right) are all parameters; in model (B), $\tau_S$ and $\theta_{Hl}$ are not parameters in the model while $\tau_T$ and $\theta_{Hr}$ are parameters; and finally in model (C), none of those four parameters exists in the model. Model (D) specifies a bidirectional introgression event between species A and B. ![MSCI-models](https://github.com/bpp/bpp-manual/blob/191d107413566b145360597fdfdfef3e68ce9173/figures/fig-msci-models.png?raw=true) **Four different types of MSC-I models implemented in BPP** The models are specified using the extended Newick format, as follows: ``` (A): ((A, (C)H[&phi=0.5,&tau-parent=yes])S, (H[&tau-parent=yes], B) T)R; (B): ((A, (C)H[&phi=0.5,&tau-parent=no])S, (H[&tau-parent=yes], B) T)R; (C): ((A, (C)H[&phi=0.5,&tau-parent=no])S, (H[&tau-parent=no], B) T)R; (D): ((A, (B) Y[&phi=0.3])X, (X[&phi=0.1])Y)R; (D): ((A, (B)Y)X, (X)Y)R; (D): ((A, Y)X, (B, X)Y) R; (D): ((A, y)x, (B, x) y) r; ``` The control file variable `phiprior` specifies the prior probability for the $\varphi$ parameter, which is a beta distribution with parameters $a$ and $b$. The syntax is ``` phiprior = a b ``` where $a$ and $b$ are positive numbers. For example, `phiprior = 1 1` specifies the beta prior beta($a, b$) with $a = 1$ and $b = 1$ for $\varphi$, the introgression probability under the MSC-I model. #### Autogeneration of extended newick graphs Starting with BPP release 4.4.0, we have changed the definition of introgression probability ($\varphi$) so that it is assigned to the horizontal (introgression) branch. This note illustrates the `msci-create` feature in bpp, which generates extended Newick notation for the MSC-I model from a data file which contains a binary species tree with introgression events specified using source and target branches on the tree. The `--msci-create` option has been available since [BPP 4.2.1](https://github.com/bpp/bpp/issues/114). The user prepares a input file (named `msci.txt` in our example), which includes a binary species tree in Newick format and specifies the introgression events by identifying the source and target branches involved in the introgression. Run bpp using the following command: ``` $ bpp --msci-create msci.txt ``` This generates extended Newick notation for the MSC-I model, which can be copied into a bpp control file. The input file `msci.txt` uses four commands: `tree`, `define`, `hybridization`, and `bidirection`. - `tree` defines the binary species tree, in Newick format. - `define` defines an ancestral species or internal node label, as the most recent common ancestor of the tip species. Alternatively you can label the internal nodes on the Newick tree. - `hybridization` defines a hybridization/introgression event by specifying the source and target branches (which represent the source and target populations involved). - `bidirection` defines a bidirectional introgression (BDI) event. One may think of the binary species tree as describing the history of species divergences and add introgression events onto it as new (horizontal) branches. The introgression probability is assigned to the newly created introgression branch. #### Examples using --msci-create A series of examples follow. The introgression graph is displayed in a figure followed by the `msci=create` encoding and the resulting extended Newick format. ![Fig-MSCI2](https://github.com/bpp/bpp-manual/blob/36012225be26f3b82593c9a35eaaa3e7217a5fd0/figures/msci_example2.png?raw=true) ``` # Model A, version 1 tree (A,(B,C)); define T as B,C define R as A,B hybridization R A, T C as S H tau=yes, yes phi=0.10 # Model A, version 2 tree (A,(B,C)T)R; hybridization R A, T C as S H tau=yes, yes phi=0.10 # The generated Newick notation for model A is # ((H[&phi=0.1,tau-parent=yes],A)S, (B,(C)H[&phi=0.9,tau-parent=yes])T)R; # Model B1 tree (A,(B,C)T)R; hybridization R A, T C as S H tau=no, yes phi=0.10 # The generated Newick notation for model B1 is # ((H[&phi=0.1,tau-parent=no],A)S, (B,(C)H[&phi=0.9,tau-parent=yes])T)R; # Model B2 tree ((A,C)S,B)R; hybridization R B, S C as T H tau=no, yes phi=0.10 # The generated Newick notation for model B2 is # ((A,(C)H[&phi=0.9,tau-parent=yes])S, (H[&phi=0.1,tau-parent=no],B)T)R; # Model C tree (A,(B,C)T)R; hybridization R A, T C as S H tau=no, no phi=0.40 # The generated the Newick notation for model C is # ((H[&phi=0.4,tau-parent=no],A)S, (B,(C)H[&phi=0.6,tau-parent=no])T)R; # Model D tree (A,B)R; bidirection  A R, B R as X Y phi=0.1,0.2 ``` The above figure illustrates the specifications of models A, B, C, and D of [Flouri et al 2020](https://doi.org/10.1093/molbev/msz296). In models A, B1, and C, the introgression event is from the source branch RA to the target branch TC, with a new introgression branch SH created. The introgression probability ($\varphi = 0.1$) is assigned to the newly added introgression branch (SH) while $1 – \varphi$ is assigned to the other parent branch (TH). Models A, B1, and C are distinguished by using the keyword `tau` (which means the same as `tau-parent` in the Newick notation of the MSC-I model). Model A assumes that the two parental species S and T of the hybridization node H have distinct ages from H, with $\tau_S > \tau_H$ and $\tau_T > \tau_H$. Thus we have `hybridization R A, T C as S H tau=yes, yes phi=0.10`: the first `yes` means that parental node S on the source branch RA has a distinct tau from node H and that the second `yes` means that parental node T on the target branch TC has a distinct tau from node H. Note that there are often multiple ways of specifying the same MSC-I model and here you can specify model A by starting with the binary tree `((A, C), B)` and adding branch TH as a introgression event. Model A might appear to be nonsensical biologically since only contemporary species can exchange migrants. Nevertheless, model A may be used to represent introgressions from a ghost species not sampled in the sequence dataset. Suppose at time $\tau_S$ a speciation event generated two species SA and SH, and at $\tau_H$, species SH contributed migrants into species THC, but species SH has since become extinct or is otherwise not sampled in the data. This scenario matches model A. Model B1 assumes $\tau_S = \tau_H$ and $\tau_T > \tau_H$ so there is no new tau for parent S of the hybridization node H, and there is a new tau for parent T. Thus we have `tau=no, yes`. Model B2 works similarly. In model C, none of parents S and T has a distinct age from H, so we have `tau=no, no`. One interpretation is that species C is a hybrid species. In model D for the bidirectional introgression (BDI) model, there is no distinction of source and target branches. We specify two phi values, assigned to the introgression (horizontal) branches: $\varphi_X = 0.1$ for node X (or into node X) and $\varphi_Y = 0.2$ for node Y (see above figure). ### The MSC-M Model The multispecies-coalescent-with-migration model (or isolation-with-migration or IM model) is activated in BPP using the keyword migration. A binary species tree is specified using an extended Newick format that includes labels for internal nodes; these are subsequently used in the control file to specify the source and target populations involved in migration. #### Control file specification of MSC-M The basic model is specified in the control file as follows, using a species tree for four species (A, B, C, D) as an example: ``` ((A, B)S, (C, D)T)R; wprior = 2 200 migration = 2 A C S C ``` Here S is the AB common ancestor, T is the CD common ancestor, while R is the ABCD ancestor. Not all internal nodes need to be labeled but those involved in migration have to be. The migration line specifies 2 migration connections: one connecting source population A to target population C, and another connecting source population S to target population C, with migration rates $w_{AC}$ and $w_{SC}$. !!! warning "Version 4.8.0 Breaking Change" The option `migprior` has been replaced by `wprior` with a different parameterization. If you are upgrading from v4.7.x or earlier, you must update your control files. BPP now uses the mutation-scaled migration rate $w$, defined as: $w_{XY} = m_{XY}/\mu$ where $m_{XY}$ is the proportion of immigrants in the recipient population Y from the donor population X every generation. Previous versions (4.7.x and earlier) used the population migration rate $M_{XY} = m_{XY} \cdot N_Y$, which represents the expected number of migrants from X to Y per generation. **Converting from old to new parameterization:** $w_{XY} = 4 M_{XY} / \theta_Y$ and $M_{XY} = \theta_Y \cdot w_{XY} / 4$ See Appendix A for definitions of migration rates used in different programs. #### Prior on migration rates The control file variable `wprior = alpha beta` specifies a default gamma prior density with parameters `alpha` and `beta` for all migration rates. In the example above, both $w_{AC}$ and $w_{SC}$ are assigned the gamma prior $G(2, 200)$, with prior mean $2/200 = 0.01$. A separate gamma prior can instead be defined for each migration rate (where the migration connection is specified using the source and target populations) and this specification takes precedence. For example if we specify ``` wprior = 2 200 migration = 2 A C 2 100 S C ``` the priors are $w_{AC} \sim G(2, 100)$ with the prior mean $0.02$ and $w_{SC} \sim G(2, 200)$ from the default prior specified by `wprior`. #### Variable migration rates among loci In this model, the migration rate $w_i$ at each locus *i* varies according to a gamma distribution $w_i \approx G(\alpha_w,\alpha_w/\overline{w})$, with shape parameter $\alpha_w$ while the mean rate $\overline{w}$ is assigned the gamma prior $\overline{w} \approx G(\alpha,\beta)$. Here $\alpha_w$ is a parameter that characterizes the variation of $w_i$ among loci, with a small $\alpha_w$ (e.g., 0.5 or 1) indicating highly variable rates among loci and a large $\alpha_w$ indicating nearly constant migration rates among loci (when $\alpha_w = \infty$ all loci have the same rate). ``` wprior = 2 200 migration = 2 A C 2 100 5 S C ``` In the above example, $w_{SC}$ is applied to all loci with the prior $G(2, 200)$, but $w_{AC}$ varies among loci according to the shape parameter $\alpha_w = 5$, and the mean rate for all loci $\overline{w}_{AC}$ is assigned the gamma prior $G(2, 100)$. #### Priors and pseudo-priors When the gene flow involves ancestral species, migration may become impossible because of changes to the species divergence times in the MCMC. In the example above, if $\tau_S < \tau_T$ migration from S to C is possible during the time period $(\tau_S, \tau_T)$, but if $\tau_S > \tau_T$ migration from S to C is impossible as the two populations were never contemporary. Thus the MCMC proposal to change species divergence times $\tau_S$ or $\tau_T$ may cause the migration rate parameter $w_{SC}$ to disappear or reappear. When $w_{SC}$ is absent from the model (that is, when $\tau_S > \tau_T$), the MCMC algorithm treats it as a pseudo-parameter (written as $w^*_{SC}$) and assigns it a pseudo-prior to facilitate the trans-dimensional move. The choice of pseudo-priors affects the mixing efficiency of the MCMC, but not the correctness of the algorithm. For good mixing, one should choose the pseudo-prior to be close to the posterior of the parameter $w_{SC}$ (which can be generated by running short chains). ``` wprior = 2 200 migration = 2 A C 2 100 S C 2 200 100 250 ``` In the example above, the prior $w_{SC} \approx G(2, 200)$ is applied when $\tau_S < \tau_T$, and the pseudo-prior $w^*_{SC} \approx G(100, 250)$ with mean 0.4, is applied when $\tau_S > \tau_T$. Note that in this example the pseudo-prior is far more concentrated than the prior (with shape parameter 100 versus 2). In order to summarize the posterior from the MCMC samples, only those samples taken when the migration rate is defined should be used (that is, only samples collected when $\tau_S < \tau_T$ in our example). Samples collected when the migration rate does not exist in the model (that is, when $\tau_S > \tau_T$ in our example) approximate the pseudo-prior. In the case of four species, the following awk commands may be used to split the samples into two files (this assumes that tau_S and tau_T are in column 10 and 11) ``` awk 'BEGIN {cnt=0} NR>1{if ($10 > $11) {cnt=cnt+1;}}} END {print cnt}' out.mcmc.txt # split sample file into ab.txt (samples where t_AB > t_CD) and cd.txt (t_CD > t_AB) head -n 1 out.mcmc.txt > cd.txt; awk 'NR>1 {if ($10<$11) print $0}' out.mcmc.txt >> cd.txt head -n 1 out.mcmc.txt > ab.txt; awk 'NR>1 {if ($10>$11) print $0}' out.mcmc.txt >> ab.txt ``` Use something like this if you understand the syntax. The running mean of the migration rate printed on the monitor during the MCMC run is the posterior mean after the filtering. Note that this problem of the migration rate appearing and disappearing during the MCMC may exist when a migration event involves ancestral populations, and is not limited to the balanced species tree for four species. Later we should automate the summary of the MCMC sample. Under the model of variable $w$ among loci, all migration rates for all loci may appear and disappear if the migration involves ancestral species. The following six ways of specifying the model and priors and pseudo-priors are accepted (with $w_{SC}$ as an example) ``` (a) S C (b) S C a_w (c) S C a b (d) S C a b a_w (e) S C a b pseudo_a pseudo_b (f) S C a b a_w pseudo_a pseudo_b ``` In options (a) and (b), the default prior specified with `wprior` is assigned on $w_{SC}$. ## Combined Analyses of Organelle Genomes and Sex Chromosomes When performing a combined analysis of data from autosomal, mitochondrial, chloroplast or sex chromosomes one needs to account for the differences of effective population sizes between loci from these different sources. This can be done using a *heredity multiplier* (referred to as an inheritance scalar by [Hey and Nielsen 2004](https://doi.org/10.1534/genetics.103.024182)). The `heredity` variable in the BPP control file allows for such differences. Here are some examples: ``` heredity = 0 # (0: No variation) heredity = 1 4 4 # (1: estimate, & a_gamma b_gamma) heredity = 2 heredity.txt # (2: from file) ``` The specification `heredity = 0` is the default and means that $\theta$ is the same for all loci. `heredity = 1` or `2` specifies two models that allow $\theta$ to vary among loci, which may be useful for combined analysis of data from autosomal, mitochondrial, X and Y loci. With such mixed data, the effective population sizes are different among loci, so that a heredity multiplier ([Hey and Nielsen 2004](https://doi.org/10.1534/genetics.103.024182)) should be applied. Other factors such as natural selection may also cause $\theta$ to deviate from the neutral expectation. BPP implements two options for this. The first option (`heredity = 1`) is to estimate the multipliers from the data, using a gamma prior with parameters $\alpha$ and $\beta$ specified by the user. In the example above, a gamma prior $G(4, 4)$, with mean $4/4 = 1$, is specified for the multiplier for each locus. The MCMC should then generate a posterior for the multiplier for each locus. The second option (`heredity = 2`) is for the user to specify the multipliers in a file, and the multipliers will then be used as fixed constants in the MCMC run. The file simply contains as many numerical values as the number of loci, separated by spaces or line breaks. | Genome | Heredity scalar | |:--------------|:----------------| | Autosome | 1 | | X chromosome | 0.75 | | Y chromosome | 0.25 | | Mitochondrial | 0.25 | Table. Examples of Common Heredity Scalars Note.--- The effect of the locus-specific mutation rates and the locus-specific heredity multipliers are different. A locus rate is used to multiply all $\theta$s and $\tau$s for the locus, while a heredity multiplier is used to multiply all $\theta$ parameters for the locus but not the $\tau$s. Nevertheless, those parameters are quite likely to be strongly correlated, especially when the species tree is small. ## Methods of Analysis Here we describe the specifics of the four different types of analyses that BPP can perform and illustrate them using several of the example control files provided in the `examples` subdirectory of the BPP distribution. Briefly, BPP is capable of analyzing sequence data using 4 distinct phylogeographical models. The 4 models are determined by the combinations of two possible settings for each of the two Boolean control file variables `speciesdelimitation` (0 = no species delimitation, 1 = species delimitation) and `speciestree` (0 = fixed species tree, 1 = inferred species tree). Therefore, analysis A00 infers the biogeographical parameters $\theta$ and $\tau$ on a fixed (user specified) tree with a fixed (delimited) number of species/populations as described in ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645). Analysis A10 delimits the number of species and infers biogeographical parameters using a fixed species tree (guide tree) as described in [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107). Analysis A01 infers the species tree and biogeographical parameters using a fixed number of species (delimitation) as described in [Rannala and Yang 2017](https://doi.org/10.1093/sysbio/syw119). Analysis A11 jointly infers the species tree, species delimitation and biogeographical parameters as described in [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279). Under Model A00 it is also possible to include a model of either instantaneous introgression or ongoing gene flow -- examples of such analyses are given in [Introgression and Migration Models](#introgression-and-migration-models). For each of these 4 major types of analyses the control file specifications will differ as well as the form of output printed to the screen and output files. Therefore, we will include two subsections in our discussion of each analysis method: input and output. Here we only consider the control file variables specific to the methods of analysis under consideration. For a general discussion of other control file variables see the example control file discussed in [Control File](#control-file). Note that for most analyses the output printed to screen during the run and that printed to the specified output file will be identical. It is important to note that only analysis A00 applies within-model inference. There is one specified model in that case (the MSC with a fixed species tree) and the parameters are well defined, with the objective being to estimate those parameters. In contrast, analyses A01, A10, and A11 all apply trans-model inference. They move between different models, and the main objective is to calculate the posterior probabilities of the models. Each is an instance of the MSC model, but the species delimitation (the number and nature of the species) and/or the species phylogeny may differ between models. In analyses A01, A10, and A11, the prior specified using control file variable `thetaprior` applies to all $\theta$ parameters in all models. Similarly, all MSC models specifying two or more delimited species have a parameter $\tau_0$ (the divergence time of the root), and these parameters (for different models) are assigned the same prior, specified by the `tauprior` control file variable. ### A00: Demography and Divergence Time Estimation Analysis A00`(speciesdelimitation = 0, speciestree = 0)`, requires a user-specified species tree topology, given by variable `species&tree` in the control file, and infers the $\theta$ and $\tau$ parameters on this fixed tree. The theory underlying estimation of species divergence times and population sizes ($\tau$s and $\theta$s) under the MSC model when the species phylogeny is given is described in ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645). Note that if there are $s$ species in the species tree, the model will minimally contain the following parameters: $s – 1$ species divergence times ($\tau$s) and $s – 1$ ancestral species $\theta$s. A $\theta$ parameter will also exist for each contemporary species with at least two sequences present at some locus. If the data are for a single species/population, the model will contain one parameter only, $\theta$ for that species. If there is only one species, the MSC model becomes the single-population Kingman coalescent. We will provide examples of both a single population analysis (using the control file `yu2001.bpp.ctl`) and a multi-species analysis (using the control file `A00.bpp.ctl`). **How many $\theta$ and $\tau$ parameters exist?** The number of parameters that BPP includes in the MSC or MSC-I/MSC-M models depends on the data configuration. The program always includes $\theta$ and $\tau$ parameters for each internal node (ancestral species) on the species tree, but includes a $\theta$ for an extant species (tip) if and only if that species has at least 2 sequences at some loci. If two or more sequences for a species are specified in the control file but there are at most 0 or 1 sequence per locus for that species in the sequence data file, $\theta$ for that species will not be identifiable or estimable. In that case, the posterior for the parameter will be the prior. Nevertheless, other results (such as the posterior distribution for other parameters or posterior probabilities for species trees and species delimitations) will still be correct. The same applies to other more complex cases of missing data and parameter non-identifiability. As an example, suppose the species tree and the numbers of sequences (for species A, B, C and D from left to right) for two kinds of loci are as follows: ((A,B), (C,D)) locus configuration 1: 1 1 0 0 locus configuration 2: 2 0 0 1 In this case, $\theta_A$, $\theta_{AB}$, $\theta_{ABCD}$, $\tau_{AB}$, and $\tau_{ABCD}$ are identifiable while $\theta_B$, $\theta_C$, $\theta_D$, $\theta_{CD}$ and $\tau_{CD}$ are not. It is impossible to estimate $\theta_D$ , $\theta_{CD}$ and $\tau_{CD}$ as no data are available from species C. #### Input A00 (single population) Our first example is a dataset of human sequence data generated by [Yu et al 2001](https://doi.org/10.1093/oxfordjournals.molbev.a003795). The control file is found in the BPP distribution at location: ``` examples/yu2001/yu2001.bpp.ctl ``` The data comprise 61 phased sequences for a single locus that will be analyzed to estimate the single parameter $\theta$. There is no need to specify an Imap file, or tag the sequence names in the sequence file (`yu2001.txt`); the sequence names are read and then ignored. Multiple loci may be included in the sequence file. The contents of the control file `yu2001.bpp.ctl` are shown below: ``` seed = -1 seqfile = yu2001.txt jobname = out # fixed species delimitation and species tree speciesdelimitation = 0 speciestree = 0 species&tree = 1 H 61 # max number of sequences # 0: no data (prior); 1:seq like usedata = 1 # number of data sets in seqfile nloci = 1 # remove sites with ambiguity data (1:yes, 0:no)? cleandata = 0 # gamma(a, b) for theta thetaprior = gamma 2 2000 # auto (0 or 1): finetune for GBtj, GBspr, theta, tau, mix, locusrate, seqerr finetune = 1 # MCMC samples, locusrate, heredityscalars, genetrees print = 1 0 0 0 burnin = 4000 sampfreq = 2 nsample = 10000 ``` The two lines: ``` speciesdelimitation = 0 speciestree = 0 ``` specify that the species delimitation (number of species) is fixed as well as the species tree topology. In this case, there is no species tree and only a single species exists. There is no need for a species tree, so the block for specifying species names and species tree looks like this: ``` species&tree = 1 H 61 ``` The line: ``` # gamma(a, b) for theta thetaprior = gamma 2 2000 ``` specifies a gamma distribution for the prior on $\theta$ with parameters $\alpha = 2$ and $\beta = 2000$ which has mean $0.001$ and variance $5 \times 10^{-7}$. #### Output A00 (single population) When BPP is run using the above control file a summary of the progress of the MCMC analysis is printed to screen. Here, we briefly explain what this output means. The different analyses presented in this chapter all produce a similar form of output when the program is running, but the number of columns and their specific content varies somewhat. Omitting some mostly irrelevant output related to the adjustment of proposal moves, the output to screen when this control file is run appears as follows: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 thg tau mix | theta1 log-PG log-L ------------------------------------------------------------------------ -15% 0.71 0.18 0.00 0.00 0.00 0.38 0.0004 467.11472 -12721.22305 0:01 (some output omitted here) 5% 0.71 0.29 0.33 0.80 0.00 0.24 0.0004 428.40742 -12721.09413 0:05 10% 0.71 0.29 0.31 0.82 0.00 0.25 0.0003 480.92848 -12720.75185 0:06 15% 0.71 0.29 0.30 0.83 0.00 0.24 0.0003 476.37432 -12720.65199 0:08 20% 0.71 0.29 0.31 0.84 0.00 0.23 0.0003 441.36960 -12720.92326 0:09 25% 0.71 0.29 0.32 0.85 0.00 0.23 0.0003 443.51740 -12720.85802 0:10 30% 0.71 0.29 0.31 0.86 0.00 0.23 0.0003 441.93616 -12720.86033 0:11 35% 0.71 0.29 0.31 0.87 0.00 0.23 0.0003 495.50716 -12720.84565 0:12 40% 0.71 0.29 0.31 0.88 0.00 0.23 0.0003 452.82049 -12720.87810 0:13 45% 0.71 0.29 0.31 0.89 0.00 0.23 0.0003 461.93346 -12720.94582 0:15 50% 0.71 0.29 0.32 0.90 0.00 0.23 0.0004 463.85955 -12721.01005 0:16 55% 0.71 0.29 0.32 0.91 0.00 0.23 0.0004 463.42744 -12720.98006 0:17 60% 0.71 0.29 0.32 0.92 0.00 0.23 0.0004 478.52797 -12721.03353 0:18 65% 0.71 0.29 0.32 0.93 0.00 0.23 0.0004 499.97705 -12720.96168 0:19 70% 0.71 0.29 0.32 0.94 0.00 0.23 0.0004 455.32489 -12720.93611 0:20 75% 0.71 0.29 0.32 0.95 0.00 0.23 0.0004 439.21544 -12720.99190 0:22 80% 0.71 0.29 0.32 0.96 0.00 0.23 0.0004 461.56954 -12721.00180 0:23 85% 0.71 0.29 0.32 0.97 0.00 0.23 0.0004 469.53443 -12720.99823 0:24 90% 0.71 0.29 0.32 0.97 0.00 0.23 0.0004 455.60303 -12720.99213 0:25 95% 0.71 0.29 0.32 0.97 0.00 0.23 0.0004 454.07141 -12721.00745 0:26 100% 0.71 0.29 0.32 0.97 0.00 0.23 0.0004 461.01044 -12720.97549 0:27 0:27 spent in MCMC theta_1H lnL mean 0.000352 -12720.995649 median 0.000337 -12720.726000 S.D 0.000116 2.936769 min 0.000089 -12735.535000 max 0.001293 -12712.774000 2.5% 0.000170 -12727.454000 97.5% 0.000617 -12716.086000 2.5%HPD 0.000145 -12726.920000 97.5%HPD 0.000574 -12715.718000 ESS* 757.838791 1180.434956 Eff* 0.075784 0.118043 ``` Note that in these examples we are using a random seed from the computer clock so your results will differ slightly from those shown. The line: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 thg tau mix | theta1 log-PG log-L ------------------------------------------------------------------------ ``` is a header that explains the content of each column printed during the run. The `Prgs` (progress) column will appear in all analyses and indicates the percentage of the MCMC iterations that have been completed. If this number is negative it indicates that the program is still running the burn-in iterations. For example, `-15%` means that 15 percent of the burn-in iterations remain to be completed. The next 6 columns are the current acceptance proportions for different parameter proposals in the MCMC. The proposals are defined as follows: - `Gage`: proposal to change ages of nodes on gene trees - `Gspr`: proposal to change gene tree topology using SPR move - `th1`: proposal to change $\theta$ using sliding window - `thg`: proposal to change $\theta$ using Gibbs sampler - `tau`: proposal to change $\tau$ parameter - `mix`: mixing step proposal In the output the acceptance proportions are stable. The `tau` column (the $\tau$ proposal acceptance proportion) is zero; this is expected because there are no species divergence times ($\tau$s) in the model when only a single species exists. The column labeled `theta1` gives the current average (mean) value of theta in the MCMC. The final columns are the log probability of the coalescent model (`log-PG`) and the log-likelihood of the sequence data on the gene tree (`log-L`). Check that the acceptance proportions are in a reasonable range (about 20% to 40% for continuous parameters such as $\theta$). If you need to terminate the program before it finishes, you can do so using the key combination `Ctrl-C` (e.g., hold down the Control key and then press the C key). The final screen output, printed after the MCMC is completed, summarizes the posterior parameter estimates for $\theta$ (`theta_1H`) in the second column. The posterior mean in the first row, for example, is $0.000352$ and the lower and upper limits of the highest posterior density (HPD) set of values are in rows 8 and 9, respectively. In this case, the HPD credibility interval for $\theta$ is $(0.000145,0.000574)$. #### Input A00 (multiple populations) Our second example is a dataset of frog sequence data generated by [Zhou et al 2012](https://doi.org/10.1111/j.1365-294X.2011.05411.x). This same dataset will also be used to illustrate analyses A01, A10 and A11. The control file is found in the BPP distribution at location: ``` examples/frogs/A00.bpp.ctl ``` The data comprise 5 loci for 4 populations, with the number of sequences at each locus varying between 21 and 30. We will be estimating $\theta$s for each of the 4 contemporary and 3 ancestral populations as well as the 3 species divergence times ($\tau$s). The contents of the control file `A00.bpp.ctl` are shown below: ``` seed = -1 seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out # fixed number of species/populations speciesdelimitation = 0 # fixed species tree speciestree = 0 species&tree = 4 K C L H 9 7 14 2 (((K, C), L), H); # unphased data for all 4 populations phase = 1 1 1 1 # use sequence likelihood usedata = 1 nloci = 5 # do not remove sites with ambiguity data cleandata = 0 # gamma(a, b) for theta (estimate theta) thetaprior = gamma 2 2000 # gamma(a, b) for root tau & Dirichlet(a) for other tau's tauprior = gamma 2 1000 # finetune: auto-adjust step lengths during burnin finetune = 1 # MCMCsamples, locusrate, heredityscalars, genetrees, substitutionparams print = 1 0 0 0 burnin = 8000 sampfreq = 2 nsample = 100000 ``` As in the previous analysis, the two lines: ``` speciesdelimitation = 0 speciestree = 0 ``` specify that the species delimitation (number of species) is fixed as well as the species tree topology. In this case, a species tree is specified along with population labels. The block for specifying species names and the species tree looks like this: ``` species&tree = 4 K C L H 9 7 14 2 (((K, C), L), H); ``` The first line above lists the number of species (populations) which is 4 in this example. This is followed by the label for each species (these should match the labels in the map file): ``` K C L H ``` The next line specifies the maximum number of sequences per species in the same order as the species labels: ``` 9 7 14 2 ``` the actual number of sequences at any locus must be less than this value. The third line is a rooted phylogenetic tree topology specified in Newick format: ``` (((K, C), L), H); ``` The next line specifies that the loci for each of the 4 populations are unphased: ``` phase = 1 1 1 1 ``` The specified phase variables should include a 0 or 1 (Boolean) entry for each population, where the populations are assumed to be in the same order as specified above (in the `species&tree` block) and the Boolean variable indicates that the sequence data for a population are either unphased (1) or phased (0). Ambiguity characters are used to represent genotypes for unphased sequences (see [BPP control file variables](#bpp-control-file-variables)). In the example data the sequences are unphased for all 4 populations. Currently, it is not possible to have partially phased data (e.g., only some loci unphased in a particular population). The program will infer the phase as part of the analysis. Unlike the earlier single population analysis, in which we estimated a single $\theta$ parameter, in the multipopulation analysis we will be estimating multiple $\theta$s (one for each population) and divergence times ($\tau$s) between populations. Thus, we have prior distributions for both $\theta$ and $\tau_0$. The lines below specify a gamma distribution as the prior for $\theta$ with parameters $\alpha = 2$ and $\beta = 2000$ (the prior mean and variance of $\theta$ are $0.001$ and $0.0000005$): ``` # gamma(a, b) for theta thetaprior = gamma 2 2000 ``` The lines: ``` # gamma(a, b) for root tau & Dirichlet(a) for other tau's tauprior = gamma 1 1000 ``` specify a gamma distribution as a prior for the root age and the remaining $\tau$s have a uniform Dirichlet distribution conditional on the root age (the prior mean and variance of the root age are $0.001$ and $0.000001$,respectively): #### Output A00 (multiple populations) When BPP is run using the above control file a summary of the progress of the MCMC is again printed to screen as follows: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 th2 thg tau mix | theta1 theta2 theta3 tau1 tau2 tau3 log-PG log-L --------------------------------------------------------------------------------------------------------------- -3% 0.65 0.10 0.30 0.30 0.80 0.04 0.12 0.0031 0.0095 0.0068 0.0021 0.0015 0.0014 1204.63468 -4458.38814 (some output omitted here) 5% 0.64 0.26 0.30 0.30 0.85 0.28 0.27 0.0035 0.0096 0.0067 0.0019 0.0018 0.0017 1225.51650 -4443.85539 1:09 10% 0.65 0.27 0.30 0.30 0.86 0.28 0.28 0.0035 0.0095 0.0068 0.0019 0.0018 0.0018 1202.62497 -4445.50020 1:47 15% 0.64 0.27 0.30 0.30 0.87 0.28 0.28 0.0034 0.0095 0.0068 0.0019 0.0018 0.0017 1211.95125 -4446.35694 2:25 20% 0.64 0.26 0.31 0.31 0.88 0.28 0.28 0.0034 0.0095 0.0067 0.0018 0.0018 0.0017 1211.09206 -4445.74292 3:05 25% 0.64 0.26 0.31 0.31 0.89 0.28 0.28 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1222.97603 -4446.12020 3:43 30% 0.64 0.26 0.31 0.31 0.90 0.28 0.28 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1216.60217 -4445.55424 4:21 35% 0.64 0.26 0.31 0.31 0.91 0.28 0.28 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1260.59722 -4445.22113 5:00 40% 0.64 0.26 0.31 0.31 0.92 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1269.45020 -4445.20340 5:38 45% 0.65 0.26 0.31 0.31 0.93 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1243.07668 -4445.00207 6:16 50% 0.65 0.26 0.31 0.31 0.94 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1229.59858 -4445.04036 6:55 55% 0.64 0.26 0.31 0.31 0.95 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1209.31273 -4445.43161 7:33 60% 0.64 0.26 0.31 0.31 0.96 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1206.22162 -4445.53994 8:11 65% 0.64 0.26 0.31 0.31 0.96 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1190.78837 -4445.21618 8:50 70% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1203.82684 -4445.04124 9:30 75% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1186.32176 -4445.00578 10:09 80% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1222.75718 -4444.94315 10:47 85% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1228.39106 -4444.68671 11:26 90% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1228.72385 -4444.54844 12:04 95% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1223.25921 -4444.25455 12:43 100% 0.64 0.26 0.31 0.31 0.97 0.28 0.27 0.0034 0.0096 0.0067 0.0018 0.0018 0.0017 1199.52731 -4444.17857 13:21 13:21 spent in MCMC theta_1K theta_2C theta_3L theta_4H theta_5KCLH theta_6KCL theta_7KC tau_5KCLH tau_6KCL tau_7KC lnL mean 0.003422 0.009588 0.006683 0.003031 0.003720 0.001577 0.001648 0.001823 0.001787 0.001710 -4444.16 median 0.003351 0.009492 0.006599 0.002932 0.003642 0.001412 0.001463 0.001809 0.001773 0.001700 -4443.69 S.D 0.000747 0.001461 0.001118 0.000819 0.000944 0.000899 0.000954 0.000234 0.000230 0.000238 11.615936 min 0.001309 0.004913 0.003258 0.000866 0.001131 0.000029 0.000039 0.001047 0.001043 0.000772 -4496.72 max 0.008147 0.017569 0.012778 0.009724 0.010355 0.008461 0.008987 0.003140 0.003113 0.002847 -4403.99 2.5% 0.002171 0.006999 0.004743 0.001729 0.002100 0.000325 0.000340 0.001404 0.001375 0.001271 -4468.01 97.5% 0.005084 0.012720 0.009104 0.004919 0.005780 0.003751 0.003991 0.002319 0.002280 0.002209 -4422.74 2.5%HPD 0.002071 0.006802 0.004595 0.001566 0.001949 0.000142 0.000164 0.001396 0.001342 0.001252 -4467.06 97.5%HPD 0.004928 0.012469 0.008903 0.004650 0.005573 0.003312 0.003512 0.002308 0.002237 0.002185 -4421.90 ESS* 2721.25 8954.56 5279.29 11762.38 1199.63 4543.27 3299.32 850.33 970.61 932.00 420.40 Eff* 0.027213 0.089546 0.052793 0.117624 0.011996 0.045433 0.032993 0.008503 0.009706 0.009320 0.004204 ``` There are again 5 acceptance proportions for the 5 parameter proposals, but now the current mean values of 6 parameters are printed. The current mean values of $\theta$ for nodes 1, 2 and 3 (`mthet1 mthet2 mthet3`) are printed -- which are $\theta$s of the contemporary species K, C and L, and those of the 3 $\tau$s (`mtau1 mtau2 mtau3`) which are the divergence times for ancestral species KCLH, KCL and KC, respectively. If there are many populations/species only the first few $\tau$s will be printed to screen but all the $\tau$s will be available in the summary created at the end of the run (see below). When the run finishes, a final block is written summarizing the results. The first two rows give the mean and median of the posterior distribution for each parameter. This is followed by summaries of the statistical uncertainty: the standard deviation (S.D.) of the posterior distribution; the minimum (min) and maximum (max) values; the lower and upper bounds of the 95% [Credible Interval](https://en.wikipedia.org/wiki/Credible_interval) (2.5% and 97.5%); the lower and upper bounds of the [Highest Posterior Density (HPD) Interval](https://en.wikipedia.org/wiki/Credible_interval) (2.5%HPD and 97.5%HPD); the Effective Sample Size (ESS\*); and the Efficiency (Eff\*). Larger values for ESS\* and Eff\* indicate better mixing and more reliable estimates. The above information is also printed (along with other technical details of the run) to the output file specified by the `jobname` variable in the control file (in our example, `jobname = out` and thus the output file is called `out.txt`). A tree file in NEXUS/Newick format is also created and placed in the file named FigTree.tre which is formatted for viewing/printing using the [Figtree](http://tree.bio.ed.ac.uk/software/figtree/) program #### MCMC output file A file will be produced (with a name specified by the control file variable `jobname`) that contains the mcmc samples from the run for all the $\theta$ and $\tau$ parameters. When running under analysis method A00 the contents of this file can be analyzed using a program such as [Tracer]() to visually examine the convergence and mixing of the MCMC run. However, the other three analysis methods (A01, A10 and A11) are trans-model analyses in which the number of parameters and or/the meaning of the parameters changes as the chain runs so it is not correct to examine the trace plot without conditioning on a particular model. An example of 5 lines from the `out.mcmc.txt` file produced by the above example are given below: ```{ .yaml .no-copy } Gen theta_1K theta_2C theta_3L theta_4H theta_5KCLH theta_6KCL theta_7KC tau_5KCLH tau_6KCL tau_7KC lnL 2 0.002848 0.010633 0.006000 0.002231 0.003697 0.001594 0.000704 0.001754 0.001717 0.001708 -4466.479 4 0.002848 0.011774 0.004789 0.002231 0.003697 0.002134 0.000798 0.001754 0.001717 0.001708 -4457.097 6 0.002848 0.010660 0.005826 0.002231 0.003697 0.001250 0.000799 0.001754 0.001717 0.001705 -4448.002 8 0.002848 0.012103 0.005826 0.002231 0.002615 0.001250 0.000803 0.001754 0.001717 0.001696 -4455.900 ``` ### A01: Species Tree Estimation Analysis A01 `(speciedelimitation = 0, speciestree = 1)` assumes that the species delimitation is fixed (e.g., species assignments are as specified in the Imap file) but that the species tree is unknown. The program estimated the species tree topology in addition to the $\theta$s and $\tau$s. #### Input A01 Our example uses the dataset of frog sequence data generated by [Zhou et al 2012](https://doi.org/10.1111/j.1365-294X.2011.05411.x) that we considered previously. The contents of the control file `A01.bpp.ctl` are shown below: ``` seed = -1 seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out speciesdelimitation = 0 * fixed species tree speciestree = 1 0.4 0.2 0.1 * speciestree pSlider ExpandRatio ShrinkRatio speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted species&tree = 4 K C L H 9 7 14 2 ((K, C), (L, H)); phase = 1 1 1 1 usedata = 1 * 0: no data (prior); 1:seq like nloci = 5 * number of data sets in seqfile cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)? thetaprior = gamma 2 2000 # gamma(a, b) for theta (estimate theta) tauprior = gamma 2 1000 # gamma(a, b) for root tau & Dirichlet(a) for other tau's finetune = 1 # finetune: auto-adjust step lengths during burnin print = 1 0 0 0 * MCMC samples, locusrate, heredityscalars, Genetrees burnin = 8000 sampfreq = 2 nsample = 100000 ``` The two lines: ``` speciesdelimitation = 0 * fixed species delimitation speciestree = 1 0.4 0.2 0.1 * speciestree pSlider ExpandRatio ShrinkRatio ``` now specify that the species delimitation is fixed and the species tree topology is being estimated `speciestree = 1`. The 3 parameters following the `speciestree` variable specify the probability that the *snakes and ladders* (SNL) move is used, rather than the subtree-pruning-regrafting (SPR) move to modify the species tree and the expand/shrink ratios for the SNL move; these proportions can be adjusted on the interval $(0,1)$ and will affect mixing of the MCMC. For more details, see the description of the control file variable [14 speciestree](#14-speciestree). If no parameters are specified after `speciestree = 1` then the SPR move is used exclusively. See [Rannala and Yang 2017](https://doi.org/10.1093/sysbio/syw119) for a description of the proposals available to change species trees in the MCMC. The line: ``` speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted ``` specifies the prior distribution on the species tree topologies. The specification `speciesmodelprior = 1` uses a prior that gives equal probability to all rooted trees apriori. See the description of the control file variable [15 speciesmodelprior](#15-speciesmodelprior) for more details. The available priors are described in [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279). #### Output A01 When BPP is run using the above control file a summary of the progress of the MCMC is again printed to screen as follows: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 th2 thg tau mix Sspr Ssnl | theta1 tau1 log-PG log-L -------------------------------------------------------------------------------------------------- (some output omitted here) 5% 0.64 0.26 0.31 0.31 0.85 0.28 0.27 0.0980 0.0000 0.0037 0.0019 1213.12189 -4448.22572 1:00 10% 0.64 0.26 0.30 0.30 0.86 0.28 0.27 0.0865 0.0000 0.0037 0.0019 1212.13271 -4447.30801 1:37 15% 0.64 0.26 0.30 0.30 0.87 0.27 0.27 0.0885 0.0000 0.0037 0.0019 1205.76722 -4446.44035 2:20 20% 0.64 0.26 0.30 0.30 0.88 0.27 0.27 0.0867 0.0000 0.0037 0.0019 1209.18617 -4445.20436 3:03 25% 0.64 0.26 0.31 0.31 0.89 0.27 0.27 0.0919 0.0000 0.0038 0.0019 1189.24742 -4444.02665 3:45 30% 0.65 0.26 0.31 0.31 0.90 0.27 0.27 0.0959 0.0000 0.0038 0.0018 1218.30083 -4443.50919 4:28 35% 0.65 0.26 0.31 0.31 0.91 0.27 0.27 0.0908 0.0000 0.0038 0.0019 1271.24840 -4443.30946 5:13 40% 0.65 0.26 0.31 0.31 0.92 0.27 0.27 0.0903 0.0000 0.0038 0.0018 1221.51424 -4442.93455 5:56 45% 0.65 0.26 0.31 0.31 0.93 0.28 0.27 0.0865 0.0000 0.0037 0.0018 1179.64791 -4442.58322 6:39 50% 0.65 0.26 0.31 0.31 0.94 0.28 0.27 0.0883 0.0000 0.0038 0.0018 1220.10587 -4442.33168 7:21 55% 0.65 0.26 0.31 0.31 0.95 0.28 0.27 0.0882 0.0000 0.0037 0.0018 1230.47093 -4442.53042 8:03 60% 0.65 0.26 0.31 0.31 0.96 0.28 0.27 0.0900 0.0000 0.0037 0.0019 1210.69277 -4443.13901 8:46 65% 0.65 0.26 0.31 0.31 0.96 0.28 0.27 0.0906 0.0000 0.0037 0.0019 1259.02426 -4443.42635 9:29 70% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0907 0.0000 0.0037 0.0019 1227.92007 -4443.74417 10:11 75% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0907 0.0000 0.0037 0.0019 1254.82843 -4443.96093 10:55 80% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0912 0.0000 0.0037 0.0019 1211.86344 -4443.93739 11:40 85% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0922 0.0000 0.0037 0.0019 1205.06522 -4443.99758 12:23 90% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0924 0.0000 0.0037 0.0019 1191.95318 -4444.06010 13:05 95% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0918 0.0000 0.0037 0.0019 1213.34484 -4444.05286 13:47 100% 0.65 0.26 0.31 0.31 0.97 0.28 0.27 0.0916 0.0000 0.0037 0.0019 1224.45872 -4443.84148 14:30 14:30 spent in MCMC Species in order: 1. K 2. C 3. L 4. H (A) Best trees in the sample (15 distinct trees in all) 33978 0.33978 0.33978 ((C, (H, L)), K); 17928 0.17928 0.51905 (C, ((H, L), K)); 8961 0.08961 0.60866 ((C, K), (H, L)); 8236 0.08236 0.69102 (C, ((H, K), L)); 6896 0.06896 0.75998 (((C, H), L), K); 4548 0.04548 0.80546 (((C, L), H), K); 4066 0.04066 0.84612 (C, (H, (K, L))); 3037 0.03037 0.87649 ((C, (H, K)), L); 2775 0.02775 0.90424 (((C, K), H), L); 2245 0.02245 0.92669 (((C, H), K), L); 2131 0.02131 0.94800 ((C, L), (H, K)); 1786 0.01786 0.96586 (((C, K), L), H); 1423 0.01423 0.98009 (((C, L), K), H); 1094 0.01094 0.99103 ((C, H), (K, L)); 897 0.00897 1.00000 ((C, (K, L)), H); (B) Best splits in the sample of trees (10 splits in all) 60867 0.608664 0011 45422 0.454215 0111 30230 0.302297 1011 13522 0.135219 1100 13404 0.134039 1001 10235 0.102349 0101 8102 0.081019 0110 8057 0.080569 1101 6057 0.060569 1010 4106 0.041060 1110 (C) Majority-rule consensus tree (K, C, (L, H) #0.608664); (D) Best tree (or trees from the mastertree file) with support values ((C, (H, L) #0.608664) #0.454215, K); [P = 0.339777] ``` The first part of this output are the current acceptance proportions and running means for different parameters as the MCMC was run. Columns 2 to 6 are the five acceptance proportions for the conventional MCMC moves, as discussed above. Columns 7 and 8 are the acceptance proportions for moves that change the species tree topology, both are much lower than the other parameter proposals, this is typical and it is often not possible to improve these acceptance proportions. Columns 9 and 10 are posterior means of $\theta$ and $\tau$ for the root population (-1 is printed in column 9 if $\theta$s are integrated out). The last two numbers are the log MSC gene-tree density ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645)) and the average log sequence likelihood (Felsenstein 1981). Next, there are 4 sections: A, B, C and D summarizing the results. Section (A) lists the species trees in decreasing order of posterior probabilities. From this you can easily identify the 95% or 99% credibility set of species trees. Section (B) lists the splits (or bipartitions) and their posterior probabilities. The splits here take into account the location of the root, and may be different from the splits for unrooted trees. Section (C) prints the majority-rule consensus tree, with posterior probabilities for nodes. #### MCMC output file The MCMC sample of species trees is collected in the file out.mcmc.txt. Below are five lines from that file. The numbers after `:` are the branch lengths ($\tau$s), while those after `#` are the $\theta$s. ```{ .yaml .no-copy } (K #0.002983: 0.001998, ((C #0.009671: 0.001684, H #0.003031: 0.001684) #0.001835: 0.000250, L #0.006266: 0.001934) #0.000560: 0.000065) #0.002770; (K #0.003247: 0.002176, ((C #0.010528: 0.001871, H #0.004726: 0.001871) #0.002902: 0.000147, L #0.006615: 0.002017) #0.000950: 0.000158) #0.003016; (K #0.003234: 0.002167, ((C #0.010486: 0.001843, H #0.003294: 0.001843) #0.002186: 0.000147, L #0.007822: 0.001989) #0.001087: 0.000178) #0.003003; (K #0.002677: 0.001794, ((C #0.010973: 0.001516, H #0.001906: 0.001516) #0.001821: 0.000093, L #0.005827: 0.001609) #0.001261: 0.000185) #0.002486; (K #0.002677: 0.001794, ((C #0.009869: 0.001516, H #0.003071: 0.001516) #0.001821: 0.000109, L #0.005788: 0.001625) #0.002495: 0.000169) #0.002486; ``` ### A10: Species Delimitation Analysis A10 `(speciedelimitation = 1, speciestree = 0)` assumes that the species delimitation is unknown (e.g., population assignments as specified in the Imap file may not represent distinct species) but that the species *guide tree* is fixed (known). The program calculates the posterior probabilities of different possible species delimitations as well as of $\theta$s and $\tau$s. #### Input A10 Our example uses the dataset of frog sequence data generated by [Zhou et al 2012](https://doi.org/10.1111/j.1365-294X.2011.05411.x) that we considered previously. The contents of the control file `A10.bpp.ctl` are shown below: ``` seed = -1 seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out speciesdelimitation = 1 1 2 1 * species delimitation rjMCMC algorithm1 finetune (a m) speciestree = 0 * species tree NNI/SPR speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted species&tree = 4 K C L H 9 7 14 2 ((K, C), (L, H)); phase = 1 1 1 1 usedata = 1 * 0: no data (prior); 1:seq like nloci = 1 * number of data sets in seqfile cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)? thetaprior = gamma 2 2000 # gamma(a, b) for theta (estimate theta) tauprior = gamma 2 1000 # gamma(a, b) for root tau & Dirichlet(a) for other tau's finetune = 1 # finetune: auto-adjust step lengths during burnin print = 1 0 0 0 * MCMC samples, locusrate, heredityscalars, Genetrees burnin = 8000 sampfreq = 2 nsample = 100000 ``` the lines: ``` speciesdelimitation = 1 1 2 1 * species delimitation rjMCMC algorithm1 finetune (a m) speciestree = 0 * species tree NNI/SPR ``` specify that the species tree is fixed `speciestree = 0` (a guide tree is used as specified in the control file) and the rjMCMC algorithm 1, with $\alpha = 2$ and $m = 1$ in equations 6 and 7 of [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107) is used for species delimitation. Another possible specification for delimitation would be: ``` speciesdelimitation = 1 0 2 ``` which would specify rjMCMC algorithm 0 with $\epsilon = 2$ in equations 3 and 4 of [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107). The two algorithms in theory should produce identical results. The line: ``` speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted ``` specifies the prior distribution on topologies (see above) Prior 0 means equal probabilities for labeled histories (which are rooted trees with internal nodes ordered by their age). This is the prior used by [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107). Prior 1 means equal probabilities for rooted trees (now the default). The prior with user specified probabilities for nodes described by [Rannala and Yang 2013](https://doi.org/10.1534/genetics.112.149039) was removed from BPP after version 2.3. We also reduced the number of loci to 1 (`nloci = 1`) in this control file because otherwise species delimitation model 5 (three species) has a posterior probability of 1. #### Output A10 The species delimitation models that can be generated from the fixed guide tree will be listed in the output, together with their prior probabilities calculated by BPP. (As a check, if you use `usedata = 0`, the MCMC should be sampling from this prior distribution.) The species delimitation model (for 4 populations as in this data set) is represented using three 0-1 flags for the three interior nodes 5, 6, 7 in the guide tree, with 0 for 'collapsed' and 1 for 'resolved'. Note that the tips in the guide tree are numbered 1, 2, $\cdots, s$ for $s$ potential species, while the interior (ancestral) nodes are numbered $s + 1, s + 2, \cdots, 2s – 1$, with $s + 1$ to be the root of the guide tree. The numbering is through a tree-traversal algorithm, fixed by the program. This same order is used to specify the divergence time parameters ($\tau$s), so you can work out the order by looking at the list of nodes in the screen output (look at the "population by population table", "# species divergence times in the order:", etc.). When BPP is run using the above control file a summary of the progress of the MCMC is again printed to screen as follows: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 th2 thg tau mix rj | np del mldp theta1 tau1 log-PG log-L ----------------------------------------------------------------------------------------------------------- (some output omitted here) 5% 0.62 0.25 0.30 0.30 0.85 0.33 0.29 0.0061 10 111 P[5]=0.9726 0.0026 0.0009 217.01433 -972.53408 0:08 10% 0.63 0.25 0.30 0.30 0.86 0.34 0.29 0.0037 10 111 P[5]=0.9845 0.0026 0.0009 229.43024 -972.60186 0:12 15% 0.63 0.25 0.30 0.30 0.87 0.34 0.29 0.0033 10 111 P[5]=0.9886 0.0026 0.0009 224.68967 -972.41862 0:16 20% 0.63 0.25 0.30 0.30 0.88 0.34 0.29 0.0037 10 111 P[5]=0.9885 0.0026 0.0009 210.93188 -972.34981 0:21 25% 0.63 0.25 0.30 0.30 0.89 0.33 0.29 0.0038 10 111 P[5]=0.9892 0.0026 0.0009 224.49484 -972.24813 0:26 30% 0.63 0.25 0.30 0.30 0.90 0.33 0.28 0.0035 10 111 P[5]=0.9901 0.0026 0.0009 216.40848 -972.27563 0:30 35% 0.63 0.25 0.30 0.30 0.91 0.33 0.28 0.0033 10 111 P[5]=0.9907 0.0026 0.0009 224.00442 -972.50421 0:35 40% 0.63 0.25 0.30 0.30 0.92 0.33 0.28 0.0031 10 111 P[5]=0.9912 0.0026 0.0009 213.49838 -972.43257 0:40 45% 0.63 0.25 0.30 0.30 0.93 0.33 0.29 0.0034 10 111 P[5]=0.9905 0.0026 0.0009 209.30067 -972.37969 0:45 50% 0.63 0.25 0.30 0.30 0.94 0.33 0.29 0.0033 10 111 P[5]=0.9906 0.0026 0.0009 232.27193 -972.51538 0:50 55% 0.63 0.25 0.30 0.30 0.95 0.33 0.29 0.0035 10 111 P[5]=0.9900 0.0026 0.0009 203.35457 -972.50992 0:55 60% 0.63 0.25 0.30 0.30 0.96 0.33 0.29 0.0037 10 111 P[5]=0.9895 0.0026 0.0009 206.99225 -972.50202 1:00 65% 0.63 0.25 0.30 0.30 0.96 0.33 0.28 0.0038 10 111 P[5]=0.9898 0.0026 0.0009 217.40496 -972.48296 1:05 70% 0.63 0.25 0.30 0.30 0.97 0.33 0.29 0.0038 10 111 P[5]=0.9893 0.0026 0.0009 216.82582 -972.50246 1:10 75% 0.63 0.25 0.30 0.30 0.97 0.33 0.29 0.0037 10 111 P[5]=0.9896 0.0026 0.0009 204.99588 -972.56907 1:15 80% 0.63 0.25 0.30 0.30 0.97 0.33 0.29 0.0036 10 111 P[5]=0.9897 0.0026 0.0009 217.51889 -972.54354 1:21 85% 0.63 0.25 0.30 0.30 0.97 0.33 0.28 0.0037 10 111 P[5]=0.9893 0.0026 0.0009 212.80419 -972.51050 1:27 90% 0.63 0.25 0.30 0.30 0.97 0.33 0.29 0.0037 10 111 P[5]=0.9896 0.0026 0.0009 223.34976 -972.47616 1:33 95% 0.63 0.25 0.30 0.30 0.97 0.33 0.29 0.0036 10 111 P[5]=0.9893 0.0026 0.0009 215.63544 -972.48464 1:39 100% 0.63 0.25 0.30 0.30 0.97 0.33 0.28 0.0037 10 111 P[5]=0.9895 0.0026 0.0009 201.02451 -972.47685 1:45 1:45 spent in MCMC Summarizing the species-delimitation sample in file out.mcmc.txt Number of species-delimitation models = 5 model prior posterior 1 000 0.200000 0.000000 2 100 0.200000 0.000000 3 101 0.200000 0.004630 4 110 0.200000 0.005960 5 111 0.200000 0.989410 Order of ancestral nodes: KCLH KC LH Guide tree with posterior probability for presence of nodes: ((K, C)#0.995370, (L, H)#0.994040)#1.000000;; ``` The starting species delimitation model is generated by choosing at random one of the models defined by the guide tree or starting tree. After the chain has started, columns 2-6 give the acceptance proportions for the conventional MCMC moves discussed above while column 7 (rj) now shows the acceptance proportion for reversible jump proposals (collapsing or expanding species delimitations on the guide tree). In general the larger this proportion, the more efficient the rjMCMC algorithm is. However there is no optimal acceptance proportion for the rjMCMC move, and a value close to 0 may not necessarily mean a problem. If one model has posterior probability close to 1, the acceptance proportion should be near 0 as well. Thus both poor mixing of the rjMCMC algorithm and extreme posterior model probabilities can cause the acceptance proportion for the rjMCMC to be close to 0. It has been noted that if the rjMCMC algorithm is suffering from poor mixing, different starting species trees often lead to different results. Next there are three numbers related to the rjMCMC move. The rest of the line shows the posterior means for $\theta$ and $\tau_0$ of the root ancestral population (which exist in all species-tree or species delimitation models), the current log MSC density and average log sequence likelihood. The three numbers in columns 8-10 that are related to the rjMCMC move, for example, ```{ .yaml .no-copy } 10 111 P[5]=0.9895 ``` represent the number of parameters in the current model (np), in this case 10, the species delimitation, in this case 111, and the current posterior probability of the most probable delimitation P[5]=0.9895. This means that delimitation model 5 (111) is the most favoured model, with the posterior at 0.9895. After the MCMC is finished, the program will summarize the sample. This output is: ```{ .yaml .no-copy } Summarizing the species-delimitation sample in file out.mcmc.txt Number of species-delimitation models = 5 model prior posterior 1 000 0.200000 0.000000 2 100 0.200000 0.000000 3 101 0.200000 0.004630 4 110 0.200000 0.005960 5 111 0.200000 0.989410 Order of ancestral nodes: KCLH KC LH Guide tree with posterior probability for presence of nodes: ((K, C)#0.995370, (L, H)#0.994040)#1.000000;; ``` The five delimitation models are listed, together with their posterior and prior probabilities. #### MCMC output file The MCMC output was saved in the file out.mcmc.txt. As the number of parameters changes when the rjMCMC moves between models, the mcmc sample file may not be very useful, so you can often ignore it. Right now the header line is generated using the starting species tree and should be ignored. After the header line, each line of output has the following numbers, separated by TABs: iteration number, the number of parameters, the tree, the sampled parameter values, and lnL. Here are the first 6 lines of out.mcmc.txt: ```{ .yaml .no-copy } Gen np tree theta_1K theta_2C theta_5KCLH theta_6KC theta_7LH tau_5KCLH tau_6KC lnL 2 10 111 0.001207 0.003438 0.003433 0.000722 0.004404 0.003324 0.000834 0.000736 0.000516 0.000528 -972.035 4 10 111 0.001207 0.003438 0.001733 0.000874 0.004570 0.003914 0.001036 0.000736 0.000596 0.000528 -971.105 6 10 111 0.001207 0.003438 0.003411 0.001107 0.003353 0.001307 0.001944 0.000736 0.000708 0.000386 -973.526 8 10 111 0.001698 0.004529 0.001902 0.000913 0.003032 0.000349 0.001474 0.000607 0.000584 0.000346 -972.899 10 10 111 0.001979 0.005280 0.002217 0.001065 0.004159 0.000407 0.001990 0.000707 0.000680 0.000348 -970.674 ``` If you know the unix command `grep`, you can retrieve the lines for the same tree model to summarize the posterior for parameters in that model, for example, ``` grep "3 111" out.mcmc.txt > result.Tree111.txt ``` In theory this should give you the same posterior as if you run analysis A00 with the species tree fixed at tree 111. In practice it is simpler to edit the Imap file and the control file to run analysis A00. **Running rjMCMC Algorithms.** - The rjMCMC algorithms for species delimitation allow the chain to move from one model to another but can have mixing problems. Check that similar results are produced from multiple runs using Algorithm0 and Algorithm1 regardless of the starting species tree. The starting tree is chosen by the program at random and will vary among runs. Make sure that some of your runs are started with the one-species model, some from the fully resolved tree, and some from other trees in between. If you get consistent results among runs with different starting trees and using the two algorithms, you are unlikely to have a convergence or mixing problem. If you have a computer with multicore, you can run those different combinations or replicates in different folders at the same time. - In medium-sized or large datasets with multiple loci, we have come across cases where the chain is stuck at the one-species model, or have difficulty moving into the one-species model. The problem is less common when there are only 1 or 2 loci. If you have a similar problem, you may start the analysis with 1 locus, then 2 loci, etc. to observe how the results change (you can do this by changing nloci in the control file as there is no need to change the sequence file). - The program used the burnin to adjust the step lengths for the proposals in the MCMC algorithm. If the rjMCMC stays in species model 0, which does not have any $\tau$ parameter, no information is collected during the burnin about the proposals to change $\tau$ and the automatically adjusted step length for $\tau$ can be very poor. - You probably need to evaluate the impact of the priors on $\theta$s and $\tau$. See the note about the gamma distribution later in this document. ### A11: Joint Species Delimitation and Species Tree Estimation Analysis A11 `(speciedelimitation = 1, speciestree = 1)` assumes that the species delimitation is unknown (e.g., population assignments as specified in the Imap file may not represent distinct species) and also that the species tree topology relating the species of any delimitation must be estimated (see [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279)). The program calculates the joint posterior probabilities of different possible species delimitations and species tree topologies as well as of $\theta$s and $\tau$s. #### Input A11 Our example uses the dataset of frog sequence data generated by [Zhou et al 2012](https://doi.org/10.1111/j.1365-294X.2011.05411.x) that we considered previously. The contents of the control file `A11.bpp.ctl` are shown below: ``` seed = -1 seqfile = frogs.txt Imapfile = frogs.Imap.txt jobname = out speciesdelimitation = 1 1 2 1 * species delimitation rjMCMC algorithm1 finetune (a m) speciestree = 1 0.4 0.2 0.1 * speciestree pSlider ExpandRatio ShrinkRatio speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted species&tree = 4 K C L H 9 7 14 2 (((K, C), L), H); phase = 1 1 1 1 usedata = 1 * 0: no data (prior); 1:seq like nloci = 5 * number of data sets in seqfile cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)? thetaprior = gamma 2 2000 # gamma(a, b) for theta (estimate theta) tauprior = gamma 2 1000 # gamma(a, b) for root tau & Dirichlet(a) for other tau's * heredity = 1 4 4 * locusrate = 1 5 finetune = 1 # finetune for GBtj, GBspr, theta, tau, mix, locusrate, seqerr print = 1 0 0 0 * MCMC samples, locusrate, heredityscalars, Genetrees burnin = 8000 sampfreq = 2 nsample = 100000 ``` the lines: ``` speciesdelimitation = 1 1 2 1 * species delimitation rjMCMC algorithm1 finetune (a m) speciestree = 1 0.4 0.2 0.1 * speciestree pSlider ExpandRatio ShrinkRatio ``` specify that the rjMCMC algorithm 1, with $\alpha = 2$ and $m = 1$ in equations 6 and 7 of [Yang and Rannala 2010](https://doi.org/10.1073/pnas.0913022107) is used for species delimitation and the species tree topology is being estimated `speciestree = 1`. As noted above, the 3 parameters following the `speciestree` variable specify the probability that the SNL move is used, rather than the subtree-pruning-regrafting (SPR) move to modify the species tree and the expand/shrink ratios for the SNL move; these proportions can be adjusted on the interval $(0,1)$ and will affect mixing of the MCMC. The line: ``` speciesmodelprior = 1 * 0: uniform LH; 1:uniform rooted trees; 2: uniformSLH; 3: uniformSRooted ``` specifies the prior distribution on delimitations and tree topologies. For analysis A11, `speciesmodelprior` can take the four values 0, 1, 2, 3, which mean Priors 0, 1, 2, 3, respectively. As mentioned above, Prior 1 (which is the default) assigns equal probabilities to the rooted species trees, while Prior 0 means equal probabilities for the labeled histories (rooted trees with the internal nodes ordered by age). Priors 2 and 3 assign equal probabilities for the numbers of species ($1/s$ each for $1,2,\ldots,s$ species given $s$ populations) and then divide up the probability for any specific number of species among the compatible models (of species delimitation and species phylogeny) either uniformly (Prior 3) or in proportion to the labeled histories (Prior 2). Priors 2 and 3 are discussed in [Yang and Rannala 2014](https://doi.org/10.1093/molbev/msu279). Prior 3 may be suitable when there are many populations. #### Output A11 When BPP is run using the above control file a summary of the progress of the MCMC is again printed to screen as follows: ```{ .yaml .no-copy } | Acceptance proportions | Prgs | Gage Gspr th1 th2 thg tau mix Sspr Ssnl rj | sp np mldp theta1 tau1 log-PG log-L ------------------------------------------------------------------------------------------------------------------------ (some output omitted here) 5% 0.65 0.26 0.30 0.30 0.85 0.28 0.30 0.0904 0.0000 0.0000 4 10 P[4]=1.0000 0.0041 0.0017 1239.61701 -4439.20352 0:59 10% 0.65 0.26 0.30 0.30 0.86 0.27 0.30 0.0964 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0017 1221.22458 -4437.63784 1:35 15% 0.65 0.26 0.30 0.30 0.87 0.26 0.30 0.0962 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0017 1238.23420 -4436.36392 2:13 20% 0.65 0.26 0.30 0.30 0.88 0.26 0.30 0.0998 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0016 1207.06023 -4435.09525 2:54 25% 0.65 0.26 0.30 0.30 0.89 0.26 0.30 0.1005 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1211.57636 -4434.31571 3:37 30% 0.65 0.26 0.30 0.30 0.90 0.26 0.30 0.1018 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1199.97256 -4434.32883 4:19 35% 0.65 0.26 0.30 0.30 0.91 0.26 0.30 0.1023 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1200.88102 -4433.94315 4:57 40% 0.65 0.26 0.30 0.30 0.92 0.26 0.30 0.1002 0.0000 0.0000 4 10 P[4]=1.0000 0.0044 0.0016 1191.04407 -4434.06936 5:38 45% 0.65 0.26 0.30 0.30 0.93 0.26 0.30 0.1045 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1193.71849 -4434.12002 6:21 50% 0.65 0.26 0.30 0.30 0.94 0.26 0.30 0.1036 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1247.11644 -4434.33624 7:05 55% 0.65 0.26 0.30 0.30 0.95 0.26 0.30 0.1017 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1218.84106 -4434.33556 7:43 60% 0.65 0.26 0.30 0.30 0.96 0.26 0.30 0.0994 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1233.73293 -4434.00290 8:22 65% 0.65 0.26 0.30 0.30 0.96 0.26 0.30 0.1005 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1235.85799 -4433.85179 9:06 70% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.0995 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1218.90538 -4434.20067 9:50 75% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.1002 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1214.40621 -4434.67870 10:32 80% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.0999 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1215.56021 -4435.19959 11:13 85% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.0980 0.0000 0.0000 4 10 P[4]=1.0000 0.0043 0.0016 1227.95558 -4435.47778 11:52 90% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.0963 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0016 1225.62882 -4435.73626 12:35 95% 0.65 0.26 0.30 0.30 0.97 0.27 0.30 0.0958 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0017 1244.11759 -4435.87327 13:19 100% 0.65 0.26 0.30 0.30 0.97 0.26 0.30 0.0949 0.0000 0.0000 4 10 P[4]=1.0000 0.0042 0.0017 1214.44552 -4436.19194 14:03 14:03 spent in MCMC (A) List of best models (count postP #species SpeciesTree) 24746 0.247460 0.247460 4 (C H K L) ((C, (H, L)), K); 20105 0.201050 0.448510 4 (C H K L) (C, ((H, L), K)); 10918 0.109180 0.557690 4 (C H K L) (C, ((H, K), L)); 9358 0.093580 0.651270 4 (C H K L) ((C, K), (H, L)); 6874 0.068740 0.720010 4 (C H K L) ((C, (H, K)), L); 5122 0.051220 0.771230 4 (C H K L) (((C, H), L), K); 4616 0.046160 0.817390 4 (C H K L) (C, (H, (K, L))); 3746 0.037460 0.854850 4 (C H K L) (((C, K), H), L); 3681 0.036810 0.891660 4 (C H K L) (((C, H), K), L); 3080 0.030800 0.922460 4 (C H K L) (((C, L), H), K); 1741 0.017410 0.939870 4 (C H K L) ((C, L), (H, K)); 1693 0.016930 0.956800 4 (C H K L) (((C, K), L), H); 1691 0.016910 0.973710 4 (C H K L) ((C, H), (K, L)); 1474 0.014740 0.988450 4 (C H K L) ((C, (K, L)), H); 1155 0.011550 1.000000 4 (C H K L) (((C, L), K), H); (B) 1 species delimitations & their posterior probabilities 100000 1.000000 4 (C H K L) (C) 4 delimited species & their posterior probabilities 100000 1.000000 H 100000 1.000000 C 100000 1.000000 L 100000 1.000000 K (D) Posterior probability for # of species P[1] = 0.000000 prior[1] = 0.238095 P[2] = 0.000000 prior[2] = 0.238095 P[3] = 0.000000 prior[3] = 0.285714 P[4] = 1.000000 prior[4] = 0.238095 ``` The arrangement of acceptance proportions, running means etc, are very similar to those of analysis A10, except that there are two additional columns (with headers `Sspr` and `Ssnl`) that provide the acceptance proportions for the SPR and SNL moves proposing changes to the species tree topology. There are four sections of output that summarize the results. Section (A) lists the best models in the decreasing order of the posterior probabilities. Here a model is a full MSC model that species both the species delimitation and species phylogeny. From this section, one can easily construct the 95% or 99% credibility sets of models. This section is further summarized to produce sections B, C, and D. Section (B) gives the posterior probabilities for the top few delimitations. Section (C) lists the delimited species and their posterior probabilities, and section (D) lists the posterior probability for the number of species, together with the prior probabilities calculated by BPP. #### MCMC output file The sampled trees and delimitations from the MCMC are printed to the specified file `out.mcmc.txt`. An example of 5 lines from this file is: ```{ .yaml .no-copy } (C #0.010039: 0.001675, ((L #0.006075: 0.001466, H #0.003617: 0.001466) #0.002155: 0.000205, K #0.003031: 0.001671) #0.000955: 0.000004) #0.004713; 4 (C #0.009000: 0.001462, ((L #0.006482: 0.001228, H #0.001665: 0.001228) #0.001957: 0.000230, K #0.002632: 0.001458) #0.000511: 0.000004) #0.003204; 4 ((C #0.009488: 0.001537, K #0.002525: 0.001537) #0.001082: 0.000004, (L #0.006926: 0.001384, H #0.002947: 0.001384) #0.001931: 0.000157) #0.003378; 4 ((C #0.010464: 0.001806, K #0.002966: 0.001806) #0.001272: 0.000004, (L #0.008138: 0.001487, H #0.002791: 0.001487) #0.002480: 0.000323) #0.004458; 4 ((C #0.009022: 0.001765, K #0.002899: 0.001765) #0.001243: 0.000004, (L #0.007955: 0.001449, H #0.003637: 0.001449) #0.002001: 0.000321) #0.004358; 4 ``` Because the number of species delimited and the tree topology are both potentially changing during the MCMC this file is usually not useful for summarizing the results unless a custom script filtering the MCMC samples is used. This file is used by the program to produce the summaries printed to screen at the end of the run. ### Adjusting step lengths for MCMC moves (finetune) The specified proposal step lengths are the initial values when automatic adjustment is used and are the fixed step lengths when manual specification of step lengths is used. The automatic adjustment of the finetune variable (MCMC step lengths) appears to be reliable and is recommended. However, be sure to check that the resulting acceptance proportions are neither too small nor too large. Here we provide some guidance for manual adjustments of the step lengths. A good strategy is to use automatic adjustments initially to generate good step lengths and then copy the step lengths into the control file, while continuing to use automatic adjustment; that way the next time the control file is run the automatic adjustment procedure will start with near-optimal step lengths and less optimization will therefore be required. #### Adjusting the step lengths manually (finetune = 0) First note the following line found in the control file `yu2001.bpp.ctl`. Here finetune = 1 means that MCMC step lengths will be adjusted automatically, and the specified values are used as initial values. ``` finetune = 1 ``` There are seven finetune steplengths here. They are in a fixed order and always read by the program even if the concerned proposal is not used. The first five of them are $\epsilon_1$, $\epsilon_2$, $\epsilon_3$, $\epsilon_4$, and $\epsilon_5$, described in ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645)). These are the step lengths used in the MCMC proposals that (1) change internal node ages in the gene tree, (2) prune and re-graft nodes in the gene tree (SPR), (3) update $\theta$, (4) update $\tau$s using the rubber-band algorithm, and (5) implements the mixing step. The 6th and 7th are for the proposals that change the locus rates or heredity multipliers and that change the sequencing errors, respectively. If the model assumes the same rate for all loci and does not use heredity multipliers, the 6th proposal step is not used. If the model assumes no sequencing errors, the 7th step is not used. The acceptance proportions for the first five proposals are always printed out on the screen, but those for the 6th and 7th are printed out only if the concerned proposal is used in the model. In the example above, only the first five proposals are used in the algorithm and we will manually change the step lengths so that the acceptance proportions become close to 30%. First, edit the line, changing it to: ``` finetune = 0 ``` The basic strategy for adjusting fine-tune parameters is if the acceptance proportion is too small (e.g., $<0.10$), decrease the corresponding finetune parameter. If the acceptance proportion is too large (e.g., $>0.80$), increase the finetune parameter. Run the program for a small number of iterations and look at the screen output to examine the acceptance proportions: ``` 0% 0.71 0.19 0.34 0.00 0.41 0.0004 463.81031 -12720.98567 0:04 5% 0.72 0.19 0.32 0.00 0.41 0.0003 458.65988 -12720.53426 0:04 10% 0.72 0.19 0.34 0.00 0.42 0.0004 454.99289 -12720.66361 0:05 15% 0.72 0.19 0.34 0.00 0.41 0.0004 481.20520 -12720.78839 0:06 ``` Here the first acceptance proportion, at 0.72, appears too large, which means that the corresponding finetune parameter (2 in the control file) is too small, similarly the second acceptance proportion, at 0.19, appears too small, which means that the corresponding finetune parameter is too large. Terminate the run (Ctrl-C) and edit the finetune values in the control file and then run the program again, checking the acceptance rates. Repeat this process a few times until every acceptance proportion is neither too small nor too large. In this example, changing the second finetune parameter (SPR) from $0.00001$ to $0.000003$ brings the acceptance proportion to about $0.30$, namely: ``` 10% 0.71 0.31 0.34 0.00 0.30 0.0004 447.59035 -12721.04474 0:06 15% 0.71 0.30 0.33 0.00 0.30 0.0004 487.91342 -12720.96697 0:07 20% 0.71 0.30 0.32 0.00 0.30 0.0003 430.67545 -12720.83600 0:08 ``` However, increasing the first finetune parameter (even by an order of magnitude to 20) has little effect. Sometimes it will not be possible to reduce the acceptance proportion for a parameter, as in this case (even using automatic adjustment produces an acceptance proportion of about $0.71$ in this case). The MCMC proposals discussed above are used in all four analyses (A00, A01, A10, A11), so that the description here applies to all of them. Note that the finetune parameters affect the efficiency of the MCMC or how fast one can obtain reliable results. In theory they do not change the results if runs using different finetune parameters (with potentially different acceptance proportions) are all run for long enough to generate reliable results. ## Simulating Gene Trees and Alignments Using BPP The `--simulate` option of BPP4 replaces the earlier MCCOAL program distributed with BPP versions 3.4 and earlier. This option can be used to simulate gene trees and sequence alignments at multiple loci under the MSC ([Rannala and Yang 2003](https://doi.org/10.1093/genetics/164.4.1645)) and MSC-I ([Flouri et al 2020](https://doi.org/10.1093/molbev/msz296)) models. The simulation program allows a GTR+G DNA substitution model (and its special cases) and allows among-loci heterogeneity in the process of sequence evolution. The different loci can have different exchangeability parameters ($a, b, c, d, e, f$) and base compositions in the GTR model, different overall evolutionary rates, differing extents of among-site rate variation (reflected in the alpha parameter for the gamma model), and different among-species rate drift processes (violations of the clock). See @Shi2018 for examples of such simulations. The simulation option also includes a MSC-M continuous migration model, with a user-specified migration rate matrix between species/populations. ### Running the simulation program To run BPP using the simulation option, create a suitable control file (for example, named MCcoal.ctl) and type the following: ``` bpp --simulate MCcoal.ctl ``` Always examine the screen output carefully to confirm that the program has read and executed the control file correctly. ### Simulation under the MSC model #### The control file Here we consider a simple example of a simulation under the MSC model with no migration or introgression. The example control file for simulating sequence alignments and gene trees has the following content: ``` seed = 12345 seqfile = MySeq.txt * comment out this line if you don't want seqs treefile = MyTree.tre * comment out this line if you don't want trees Imapfile = MyImap.txt * concatfile = concatenatedfile.txt * concatenated alignment modelparafile = modelparas.txt * comment out this line if you don't want seqs species&tree = 4 A B C D 3 2 1 1 ((A #0.01, B #0.01) :0.01 #.01, (C, D) :0.011 #0.01) : 0.012 #0.01; phase = 1 1 0 0 * 0: do not phase (fully phased seqs), 1: diploid unphased seqs loci&length = 100 1000 * number of loci & number of sites at each locus * locusrate = 0 (default: same rate for all loci) * locusrate = mu_bar a_mui prior * clock = 1 (default: strict clock) * clock = 2 v_bar a_vi prior dist (independent-rates) * clock = 3 v_bar a_vi prior dist (correlated-rates) model = 7 * model: 0:JC69, 7:REV (GTR) Qrates = 0 10 5 5 5 5 10 * 1: fixed; 0: dirichlet, for TC TA TG CA CG AG basefreqs = 0 10 10 10 10 * 0: random, Dirichlet(aT,aC,aA,aG), for base frequencies alpha_siterate = 0 100 20 5 * G(a, b) for alpha for sites & K for discrete gamma ``` We now go through each control file variable and explain its purpose. *seed* Use `seed = -1` to simulate different datasets each time the program is run (obtaining a random number seed from the computer clock). If you use a positive integer as a seed the data will be exactly the same each time the program is run with that control file; this is usually not what you want to do. *seqfile, treefile, Imapfile, concatfile, modelparafile* These variables specify the names of files to be generated. If you want to simulate gene trees and not sequence files, you can comment out the line for the `seqfile`. The file specified by `concatfile` will contain the sequence alignment concatenated across loci. The sequence files can become very large. The file specified by `modelparafile` contains the parameter values for the GTR model, as well as the locus rate and the rates for nodes on the species tree. There is one line of output per locus. *phase* This variable is used in the same way as in the BPP inference (data analysis) mode. Each species is assigned a switch/flag, with 0 meaning fully phased (haplotype) sequences and 1 for unphased diploid sequences. With the specifications in the example (sampling configuration 3 2 1 1 and phase = 1 1 0 0), the program will simulate 6, 4, 1, 1 sequences for A, B, C, D, respectively, and then combine every two sequences from A, B, C into one single diploid sequence with heterozygote sites represented using ambiguity characters (YRMKSW). If you use the same random number seed, the diploid data generated using the sampling configuration (3 2 1 1) with phase = 1 1 0 0 should match the fully resolved sequences generated using the sampling configuration (6 4 1 1) with phase = 0 0 0 0. The program also prints out the fully resolved sequences in a file with the suffix `_full` (or `MySeq_full.txt` in the example). *loci&length* This variable specifies the number of loci to be simulated and the length (number of sites) for each loci. In this example, 100 loci, each of 1000 sites, will be simulated. *species&tree* This block specifies the number and names of species, the species tree, and the parameters under the multispecies coalescent model, including the species divergence times ($\tau$s) and population size parameters ($\theta = 4N\mu$). For the example above, 100 loci, each of 1000 sites, will be simulated, on the species tree ((A, B), (C, D)), with 3, 2, 1, 1 sequences for A, B, C, D, so that there are 7 sequences at each locus. The divergence time parameters ($\tau$s) are after `:` in the tree, while the population size parameters ($\theta$s) are after `#`. Thus we have 0.01 for all $\theta$s, and $\tau_{ABCD} = 0.012$, $\tau_{AB} = 0.01$, and $\tau_{CD} = 0.011$. We need $\theta_A$ and $\theta_B$ because 2 or more sequences are sampled from species A and B. Parameters $\theta_C$ and $\theta_D$ are unnecessary in this case as only one sequence is sampled from C and D: if you specify them, they will be ignored by the program. Note that both $\theta$s and $\tau$s are measured by the expected number of mutations per site, so that $\theta = 0.01$ means that two sequences sampled from the population are 1% different. For reference, the estimate for the human species is $\theta_H = 0.0006$. *model* This variable can take two possible values: 0 for JC69 and 7 for GTR. The next few lines are relevant for the GTR model. *Qrates* This variable specifies the exchangeability parameters in the GTR model ($a, b, c, d, e, f$ for TC, TA, TG, CA, CG, and AG, in Yang 1994). There are two options, either to have those rates fixed for all loci or to have them sampled at random. The first has the following format (note that the first 0-1 number is a switch): Qrates = 1 2 1 1 1 1 2 * 1: fixed; for TC TA TG CA CG AG The input is one integer value (0 or 1) which acts as a switch followed by six real numbers. The line above fixes the rates at $a = 2, b = c = d = e = 1$ and $f = 2$, so that the transition rate is twice as high as the transversion rate, and the model corresponds to K80 or HKY. Note that only the relative rates matter, because the rate matrix ($Q$) is scaled so that the average rate (over the base frequencies) is 1 and branch lengths are measured in the expected number of mutations/substitutions per site. Nevertheless the program expects six rate parameters in the input here: Qrates = 0 10 5 5 5 5 10 * 0: dirichlet, for TC TA TG CA CG AG The second option, above, is to sample $a, b, c, d, e$ from the Dirichlet distribution for every locus. The above specifies Dir(10, 5, 5, 5, 5, 10), with parameters $\alpha_a = 10, \alpha_b = \alpha_c = \alpha_d = \alpha_e = 10$, and $\alpha_f = 10$. The rates generated from the Dirichlet sum to 1, and they are rescaled. Note that larger values for those parameters mean less variance, while the mean for $a$, say, is given by $\alpha_a/\alpha$ with $\alpha = \alpha_a + \alpha_b + \alpha_c + \alpha_d + \alpha_e + \alpha_f$. The specification here gives on average a transition/transversion rate ratio of 2 (if the base frequencies are fixed at $\frac{1}{4}$ each), but it varies among loci. *basefreqs* The base frequency parameters ($\pi_T, \pi_C, \pi_A, \pi_G$) can also be fixed or sampled from the Dirichlet distribution, like the `Qrates`. The two options are as follows: ``` basefreqs = 1 0.15 0.35 0.15 0.35 * 1: fixed; Base frequencies are for TCAG. basefreqs = 0 10 10 10 10 * 0: random, Dirichlet(aT, aC, aA, aG), for base frequencies ``` *alpha_siterate* The gamma shape parameter ($\alpha$) for rate variation among sites of the same locus can be fixed or sampled from a gamma distribution. The possible options are as follows. When alpha is fixed, the same value is used for all loci. Otherwise different $\alpha$s are used for different loci. Note that given alpha, the rates for sites at the same locus have a G(($\alpha, \alpha$) distribution with mean 1 [@Yang1993; @Yang1994]. The fourth option below will allow the program to sample $\alpha$ from $G(100, 20)$, with mean 5, for each locus and then use $G(\alpha, \alpha)$ to describe the among-site rate variation for the locus: ``` alpha_siterate = 1 0 * 1: alpha fixed at 0(inf), one rate for all sites at the locus alpha_siterate = 1 5.6 0 * 1: alpha fixed at 5.6, K = 0(inf) for continuous gamma alpha_siterate = 1 5.6 5 * 1: alpha fixed at 5.6, K = 5 categories for discrete gamma alpha_siterate = 0 100 20 5 * 0: alpha sampled from G(100, 20), with mean 5, K = 5 for discrete gamma ``` The variables `locusrate` and `clock` are used to specify relaxed-clock models, with rate variation among branches and among loci. These are the same models as implemented in the inference program. Please see [Models of Substitution Rate Variation](#models-of-substitution-rate-variation) for details. *locusrate* This specifies variable rates for different loci. The default is 0, meaning that all loci have the same rate. To specify variable rates for loci, use the following format: ``` locusrate = 0 (default) locusrate = mu_bar a_mui prior locusrate = mu_bar a_mui locusrate = 1.0 5.0 dir ``` Here `mu_bar` is the mean rate across loci. If `prior = dir` (gamma-Dirichlet), the total rate ($L\bar\mu$) is partitioned to generate the locus-rates ($\mu_i$) according to the Dirichlet distribution with concentration parameter $\alpha$ (`a_mui`). Large $\alpha$ (10 or 100, say) means that the rates are similar among loci, while small values (1 or 0.5, say) mean that rates are highly variable among loci. If `prior = iid` (conditional i.i.d.), the locus rates ($\mu_i$) are i.i.d. given the mean rate ($\bar\mu$): ``` mu_i ~ Gamma(a_mui, a_mui/mubar) ``` *clock* This variable is used to specify strict-clock or relaxed-clock models. The default is `clock = 1` for strict clock, while `clock = 2` is the independents-rates model, and `clock = 3` is the correlated-rates model (which is implemented for the MSC model with no introgression and is unavailable for the MSC-I model). ``` clock = 1 (default: strict clock) clock = 2 v_bar a_vi prior dist * independent-rates model clock = 2 0.1 10 dir G * gamma-dirichlet for loci & gamma for branches clock = 2 0.1 10 iid LN * iid for loci & log-normal for branches clock = 3 v_bar a_vi prior dist * correlated-rates model clock = 3 0.1 10 iid G * iid for loci & gamma for branches clock = 3 0.1 10 iid LN * iid for loci & log-normal for branches ``` The specification (`clock = 2 0.1 10 dir G`) means the following. First the average variance parameter is = 0.1. Then the total $L\bar\mu$ is partitioned into $\mu_i$ (for $i = 1, 2, \cdots, L$) according to the Dirichlet distribution with concentration parameter $\alpha = 10$ (`a_vi`). In contrast, the specification (`clock = 2 0.1 10 iid LN`) means that the variance parameter for locus $\nu_i$ is generated as $\nu_i|\bar\nu \sim G(\alpha, \alpha/\bar\nu)$ with shape parameter $\alpha = 10$ (`a_vi`) and mean $\bar\nu = 0.1$. In both prior models (`prior=dir` and `iid`), $\alpha$ (`a_vi`) is inversely related to the variance in $\mu_i$ among loci: use small values of $\alpha$ (2, 1, or 0.5) if the clock nearly holds at some loci but is seriously violated at others, and large values (such as 10 or 100) if the clock is violated in the same extent at different loci. `clock = 2` Given the locus rate $\mu_i$ (specified using the `locusrate` variable) and the variance parameter $\nu_i$ (specified using the `clock` variable) for locus $i$, the rates for species-tree branches are specified as follows. If `dist = G` (for gamma), the rate for branch $j$ at locus $i$ has the following gamma distribution with mean $\mu_i$ and variance $\nu_i$ $r_{ij} | \mu_i, \nu_i \sim G(\mu_i^2/\nu_i, \mu_i/\nu_i).$ If `dist = LN` (for log-normal), the rate has the log-normal distribution $r_{ij} | \mu_i, \nu_i \sim LN(\mu_i, \nu_i).$ `clock = 3` The correlated-rates model is specified similarly, and for the MSC model only. Given the locus rate $\mu_i$ (specified using the `locusrate` variable) and the variance parameter $\nu_i$ (specified using the `clock` variable) for locus $i$, the rates for species-tree branches are specified recursively, starting from the root towards the tips. If `dist = G` (for gamma), the rate for branch $j$ at locus $i$ has the following gamma distribution with shape parameter $\alpha = \mu_i^2/\nu_i$ and with the mean to be the rate for the parental branch of $j$ or $r_{i,p(j)}$: $r_{ij} | \mu_i, \nu_i, r_{i,p(j)} \sim G(\frac{\mu_i^2}{\nu_i}, \frac{\mu_i^2}{\nu_i r_{i,p(j)}}).$ If `dist = LN` (for log-normal), the rate for branch $j$ at locus $i$ has the log-normal distribution with the mean to be the parental rate and variance parameter $\nu_i$: $r_{ij} | \nu_i, r_{i,p(j)} \sim LN(r_{i,p(j)}, \nu_i).$ After the rates for branches (populations) on the species tree are generated for the locus, they are used to calculate the gene-tree branch lengths. Gene-tree branches residing in the same species/populations have the same rate (the rate of that species at the locus), while those residing in different species may have different rates. A branch on a gene tree may have segments residing in different species/populations, and the length of the branch (measured by the expected number of mutations/substitutions per site) is calculated by adding up those segments. #### Program output Running the program on the above control file produces the following screen output, indicating that the simulation has been successful: ```{ .yaml .no-copy } bpp --simulate MCcoal.ctl Site rates: alpha sampled from Gamma(100.000000,20.000000), K = 5 bpp v4.7.0_linux_x86_64, 15GB RAM, 8 cores https://github.com/bpp/bpp Detected CPU features: mmx sse sse2 sse3 ssse3 sse4.1 sse4.2 popcnt avx avx2 Auto-selected SIMD ISA: AVX2 Substitution model: GTR with generated rates from Dirichlet 10.000000 5.000000 5.000000 5.000000 5.000000 10.000000 Base frequencies: generated from Dirichlet a_T=10.000000 a_C=10.000000 a_A=10.000000 a_G=10.000000 4 species: A (3) B (2) C (1) D (1) ((A #0.010000: 0.010000, B #0.010000: 0.010000) #0.010000: 0.002000, (C: 0.011000, D: 0.011000) #0.010000: 0.001000) #0.010000; Map of populations and ancestors (1 in map indicates ancestor): Species 1 2 3 4 5 6 7 1 A 1 0 0 0 1 1 0 tau = 0.000000 theta = 0.010000 2 B 0 1 0 0 1 1 0 tau = 0.000000 theta = 0.010000 3 C 0 0 1 0 1 0 1 tau = 0.000000 theta = 0.000000 4 D 0 0 0 1 1 0 1 tau = 0.000000 theta = 0.000000 5 ABCD 0 0 0 0 1 0 0 tau = 0.012000 theta = 0.010000 6 AB 0 0 0 0 1 1 0 tau = 0.010000 theta = 0.010000 7 CD 0 0 0 0 1 0 1 tau = 0.011000 theta = 0.010000 Sequence data file -> MySeq.txt Trees -> MyTree.tre Model parameters for loci -> modelparas.txt Tags to species mapping (Imap) -> MyImap.txt ``` Several output files are produced. The file `MySeq.txt` contains the simulated sequence alignments in Phylip/BPP format: ```{ .yaml .no-copy } 7 1000 a1^A ACATGCCACT ACACCGGACA TTAGAGCTAC GCGAGCTCTT TGWACTCGTT ACAACCTCAC AATGTAATAG TGTAATTACT TAGCTTAAGT AMT CTTGATC TTTTCTATAR TTACATCGAA CGAATATCTG TACCTGAGCT TTTCACGATC CCTYCCTTCA TCTAGCAATC GTAAACTWTG TAAATTTCCT ACAGTCA TAY TRATAACCCA TATCATGTAA CCAGACCCCC CTWAATTTMT CGAGAATGCC CMTCTTCATT CAATCCCAAA CAGTCAAAGC GTCCACAAAT TAATTGCACA TAACCATCGT AATTGGAGAA ATTCCACCAT TCAATAGTTT ACTTGATTMA AACTAGTGGA TCCCATGTAA GCATAATCCC CCTACYCCAT ATGTATCTTA GAGT MTCGTT TGRGACTTTM CATGCCCCAG CTGTTTTTAA ATGCCAGTAA ACARCGGCGA CTACTATGTC TTACGACAGA ACGAGGTTTC CACTGAGAGT RCAAAGGA TC ACTCGTCTAT AGTTCAAAAW TCATTGCACG ATGTTCTCGA TCTACCTAAA TCTTGCCTGC TTTTTCMCCA TGCAATACAA CACCGTTAAA CAATATCTTA A TATCAATCG TACAGTCGAT GCAGTATAGC CGATACTCCT TAACCTAACT YTTCAATATG ACMTAATACT CAACTGTTGC GRTACCAATT ATCCMTGCTG TACTG ATAAC TTTATTGAAT GCCAATMCTA TWATGATCTA TTTCAGAYTT AATTCACAAT GAATTCAATA CGAATTTTCC TATCGTCAGS ACATCAATTA GATTCACAT G TKGRAAGATC TATTTTGCAT CGACAAATGC ACCACCGCTC CCTCGCCTCT AAATCGCTCA CAATCGAATT AGAAAATATA ACACWTCCAT TGATCAASCT AA AAGCGCAC AGGGAGCGAT TGGTTCGAAA ATGAGTTACA ACGTAAAGTT CTCAGCATAA TGGRGGCTAT ACACCCAAAA a2^A ACATGCCACT ACACCGGACA TTAGAGCTAC GCGAGCTCTT TGTACTCGTT ACAACCTCAC AATGTAATAG TGTAATTACT TAGCTTAAGT AAT CTTGATC TTTTCTATAA TTACATCGAA CGAATATCTG TACCTGAGCT TTTCACGATC CCTCCCTTCA TCTAGCAATC GTAAACTTTG TAAATTTCCT ACAGTCA TAT TAATAACCCA TATCATGTAA CCAGACCCCC CTTAATTTCT CGAGAATGCC CATCTTCATT CAATCCCAAA CAGTCAAAGC GTCCACAAAT TAATTGCACA TAACCATCGT AATTGGAGAA ATTCCACCAT TCAATAGTTT ACTTGATTCA AACTAGTGGA TCCCATGTAA GCATAATCCC CCTACCCCAT ATGTATCTTA GAGT CTCGTT TGGGACTTTC CATGCCCCAG CTGTTTTTAA ATGCCAGTAA ACAGCGGCGA CTACTATGTC TTACGACAGA ACGAGGTTTC CACTGAGAGT GCAAAGGA TC ACTCGTCTAT AGTTCAAAAT TCATTGCACG ATGTTCTCGA TCTACCTAAA TCTTGCCTGC TTTTTCCCCA TGCAATACAA CACCGTTAAA CAATATCTTA A TATCAATCG TACAGTCGAT GCAGTATAGC CGATACTCCT TAACCTAACT CTTCAATATG ACATAATACT CAACTGTTGC GATACCAATT ATCCATGCTG TACTG ATAAC TTTATTGAAT GCCAATCCTA TTATGATCTA TTTCAGATYT AATTCACAAT GAATTCAATA CGAATTTTCC TATCGTCAGC ACATCAATTA GATTCACAT G TGGGAAGATC TATTTTGCAT CGACAAATGC ACCACCGCTC CCTCGCCTCT AAATCGCTCA CAATCGAATT AGAAAATATA ACACTTCCAY TGATCAAGCT AA AAGCGCAC AGGGAGCGAT TGGTTCGAAA ATGAGTTACA ACGTAAAGTT CTCAGCATAA TGGGGGCTAT ACACCCAAAA ... ``` The file MyTree.tre contains the simulated gene tree for each locus in Newick format: ```{ .yaml .no-copy } ((((((a1a^A:0.000217,(a2b^A:0.000182,a3a^A:0.000182):0.000035):0.000362,a2a^A:0.000579):0.000065,a3b^A:0.000644):0 .011379,c1^C:0.012023):0.001408,(a1b^A:0.010443,(((b1a^B:0.001294,b2b^B:0.001294):0.000858,b1b^B:0.002152):0.00144 4,b2a^B:0.003596):0.006847):0.002989):0.010039,d1^D:0.023471):0.000000; [TH=0.023471, TL=0.091502] (((((a1a^A:0.000584,((a2b^A:0.000357,a3b^A:0.000357):0.000057,a3a^A:0.000414):0.000170):0.001213,a1b^A:0.001797):0 .000482,a2a^A:0.002279):0.008361,(b1a^B:0.001291,(b1b^B:0.000558,(b2a^B:0.000345,b2b^B:0.000345):0.000213):0.00073 2):0.009350):0.010901,(d1^D:0.014926,c1^C:0.014926):0.006615):0.000000; [TH=0.021541, TL=0.076273] ((((a1a^A:0.001025,((a1b^A:0.000262,a2a^A:0.000262):0.000003,a2b^A:0.000265):0.000760):0.003850,(a3a^A:0.003915,a3 b^A:0.003915):0.000961):0.005929,((b1a^B:0.000958,(b2a^B:0.000139,b2b^B:0.000139):0.000819):0.000315,b1b^B:0.00127 3):0.009531):0.005512,(d1^D:0.015947,c1^C:0.015947):0.000369):0.000000; [TH=0.016316, TL=0.072093] ... ``` The file `MyImap.txt` contains the Imap data. These files can be used as input to inference programs such as BPP along with a suitable control file that matches the simulated data . ### Simulation with migration under the MSC-M model Here we consider a simple example of a simulation under the MSC-M model with continuous migration. #### The control file The example control file for simulating sequence alignments and gene trees has the following content: ``` seed = -1 treefile = mytree.tre Imapfile = myimap.txt seqfile = mydata.txt *concatfile = concat.txt modelparafile = modelparas.txt #phase = 0 0 0 species&tree = 3 A B C 4 4 4 ((A #0.015, B #0.025) S #0.015 :0.01, C #0.025) :0.02 #0.025; loci&length = 2000 500 model = 0 * JC migration = 2 S C 0.1 A C 0.2 5 ``` There are two major differences from the previous (MSC) example. First, in the species tree the ancestral population to clade (A, B) is given a label S: ``` species&tree = 3 A B C 4 4 4 ((A #0.015, B #0.025) S #0.015 :0.01, C #0.025) :0.02 #0.025; ``` A label is needed if migration is going to be specified for the ancestral population. Second, there is a new `migration` block: ``` migration = 2 S C 0.1 A C 0.2 5 ``` Here, two migration connections are specified with `migration = 2`. The first migration connection is from source population S (ancestor of AB) to target population C and migration occurs at the rate $M_{SC} = 0.1$ migrants per generation. The second migration connection is from source population A to target population C. In this case, the migration rate $M_{AC}$ varies among loci with the mean to be 0.2 and the shape parameter to be 5. The migration rate from populations $i$ to $j$ is defined as $M_{ij} = N_j m_{ij}$ , where $m_{ij}$ is the proportion of individuals in population $j$ that are migrants from population $i$ in each generation. In the example above, $M_{SC} = 0.1$, which means that on average 0.1 individuals are immigrants from population S to population C. We use forward-time for our definition of rates as the rates are then more biologically meaningful. This control file format is used in BPP version 4.4 onward, and is different from the old MCcoal format, in which the labels and order for populations are fixed by the program. The control file included in the release ( `MCcoal.im-3s-saturated.ctl` ) simulates data on a species tree for three species, with eight migration rates. ### Simulation with tip-dating under the MSC Here we consider a simple example of a simulation under the MSC with tip-dating. #### The control file The example control file for simulating sequence alignments and gene trees has the following content: ``` seed = 1 seqfile = simulate.txt treefile=simulate_trees.txt Imapfile = simple.Imap.txt datefile = dates.txt seqDates = seqDates.txt # fixed number of species/populations *speciesdelimitation = 0 # fixed species tree species&tree = 3 A B C 4 4 2 (A #0.001, B #0.001):.007 #0.001, (C #0.001, D #0.001):.004 #.001); phase = 0 0 0 loci&length = 100 1000 clock = 1 locusrate =0 model = 0 ``` The only difference in the control file when using tip-dating is the inclusion of the `datefile` and `seqDates` options. `datefile` specifies the name/path to the file containing the sequence ages The ages are in units of expected number of substitutions, which is not the same as the inference program. The dates are assigned to populations. The number of sample dates must match the number of samples from each population. For example, in the below date file, there are 4 samples from population A with sample ages 0.00005 and 0.00006 substitutions per site. ``` A 0.00006 A 0.00005 A 0.00005 A 0.00006 B 0.00007 B 0 B 0 B 0 C 0 C 0 ``` The `seqDates` option specifies the name/path of the output file for the date file with individual names with sample dates. The `bpp` simulator automatically names the sequences based on the population names. Since the sequences are not named in the input files for the simulator, the simulator must generate a file to match the dates to the sequences. Below is the `seqDates` file that corresponds to the `datefile` shown above. ``` A^a1 0.000050 A^a2 0.000050 A^a3 0.000060 A^a4 0.000060 B^b1 0.000000 B^b2 0.000000 B^b3 0.000000 B^b4 0.000070 C^c1 0.000000 C^c2 0.000000 ``` The `seqDates` and `datefile` options must be used together in the simulator. A strict clock model (the default, which is equivalent to `clock = 1`) is required for simulating with tip dates. The simulator uses dates in units of expected substitutions, but the inference program uses dates in calendar time such as years or thousand years. The dates in the `seqDates` file can be converted to a unit of years by dividing by the per year substition rate. This is required to use the simulator to generate data to use in the inference program. For example, if we assume a substitution rate of $10^{−8}$ per year and use the same example as above, the `datefile` to use for inference would be ``` A^a1 5000 A^a2 5000 A^a3 6000 A^a4 6000 B^b1 0 B^b2 0 B^b3 0 B^b4 7000 C^c1 0 C^c2 0 ``` BPP does not generate a file with the dates in years or similar units. The user must prepare this file after running the simulator. ## Advanced Features of BPP ### Threading and Checkpointing #### Threads BPP4 can use Single-Instruction-Multiple-Data (SIMD) instruction sets (also called vector instructions) available on modern processors. It automatically detects and uses the best instruction set available so you don't have to do anything about this. However, you can manually specify and force an instruction set in the control file using the arch tag. For instance, arch=SSE forces BPP to use the SSE instruction set even if AVX (which has twice the register size of SSE and thus in theory yields twice the speedup) is present on the system. To completely disable vector instruction, one can use the arch=CPU option. The available values for the arch tag are 'CPU', 'SSE', 'AVX', 'AVX2' and 'NEON'. ![Fig-Numa](https://github.com/bpp/bpp-manual/blob/191d107413566b145360597fdfdfef3e68ce9173/figures/fig-numa.png?raw=true) Modern computer architecture using Non-Uniform momery access (NUMA). BPP4 can use multiple cores/threads. The control variable threads has the following syntax. ``` Threads = 8 19 1 ``` This means BPP will use 8 threads, starting from core/thread 19, with increment 1; in other words hardware threads 19-36 will be used. The software threads are pinned onto the hardware threads and they not allowed to migrate. We found that this helps with performance. Modern clusters based on microprocessors use a common shared-memory configuration called NUMA (for non-uniform memory access). The cluster typically consists of two or four microprocessors interconnected on a local bus to a shared memory on a single motherboard. You can use commands like `lscpu` to see the machine architecture. For example, on our server lscpu gives the following output: ``` Architecture: x86_64 CPU op-mode(s): 32-bit, 64-bit Byte Order: Little Endian CPU(s): 144 On-line CPU(s) list: 0-143 Thread(s) per core: 2 Core(s) per socket: 18 Socket(s): 4 NUMA node(s): 4 ... NUMA node0 CPU(s): 0-17,72-89 NUMA node1 CPU(s): 18-35,90-107 NUMA node2 CPU(s): 36-53,108-125 NUMA node3 CPU(s): 54-71,126-143 ``` This means that hardware threads 1-18 are on CPU0, while threads 73-90 are the hyperthreads on CPU0. The specification Threads = 8 19 1 thus uses 8 threads on CPU1. Our recommendation is that you do some tests yourself using threads = 2 or 4 or 8, and use options that work well for your servers. Adjust the step lengths first. Then use burnin=0 and a small nsample so that the running time is about 1-2 minutes. Change threads and record running time. Make sure all threads for the same BPP job are on the same CPU. Use lscpu to see the machine configuration and htop to examine the load on the hardware threads. #### Checkpointing This option instructs BPP to create a checkpoint file (save the progress) after a specified number of MCMC iterations. The MCMC loop in BPP consists of `N = burnin + sampfreq*nsample` iterations, numbered as $1,2,\cdots,N$. There are two formats, where X, Y are whole numbers: ``` checkpoint = X checkpoint = X Y ``` In the first format, a single checkpoint is created after X MCMC iterations (including burnin). In the second format, a checkpoint is created after X steps, and then additional checkpoints are created every Y steps. For example, if X = 10000 and Y = 10000, the first checkpoint is created at MCMC iteration 10000, and further checkpoints at iterations 20000, 30000, etc. The checkpoint files are named `JOBNAME.Z.chk`, where `JOBNAME` is the label specified for the `jobname` option, and `Z` is the number of the checkpoint file, starting from 1 and incrementing with every new checkpoint file. To resume from a checkpoint file, use the resume switch, i.e. ``` bpp --resume checkpoint-file.chk ``` ### Marginal Likelihood Calculations The C program BFdriver generates control files and job subscription scripts for running MCMC (using BPP or MCMCtree) to calculate the marginal likelihood (or the Bayes factor), as described in [Rannala and Yang 2017](https://doi.org/10.1093/sysbio/syw119). The program takes a control file you provide (such as bpp.ctl) and generates $K = 16$ control files with different beta values, which are used to run bpp to sample from the different power posterior distributions. The program also generates job submission scripts and submit the jobs using qsub. All generated control files and output files are in the same current directory. The frogs dataset in the bpp release (([Yang 2015](https://doi.org/10.1093/czoolo/61.5.854)) is used as the example. You need the following: a linux system with SUN grid engine managing job submission (including commonds such as `qsub`, `qstat`, `qdel`, etc.), and a C compiler. If you don't have this job submission system, you can use BFdriver to generate the control files and run the MCMC jobs from the command line. #### Compiling and running BFdriver ``` cc -o BFdriver -O3 BFdriver.c tools.c -lm BFdriver BFdriver A00.ctl 16 tmp.sh ``` You may need to edit the following two lines inside BFdriver.c, and if you do, remember to compile the program after editing. Here bpp is assumed to be on your search path. You can use a full path for the executable prorgam, such as `/home/gooduser/bin/bpp3.3`. Also the second line is for submitting the jobs using `qsub`. Here the limits are set to 4G of RAM and 360 hours of running time. Check those values if necessary (and recompile). ``` fprintf(fcommand, " echo \"bpp %s.b$I.ctl > log.b$I.txt\" > %s\n", ctlf, scriptf); fprintf(fcommand, " qsub -S /bin/bash -l h_vmem=4G -l tmem=4G -l h_rt=360:0:0 -cwd %s\n", scriptf); ``` #### Running the program Create a folder inside `frogs/`, say `bf1`: ``` mkdir frogs/bf1 cd frogs/bf1 ``` Prepare a control file (A00.ctl, say) for the A00 analysis in the folder. Check that it works. This should have a fixed species tree, which is ((K, C), (L, H)). Species delimitation and species tree estimation should be turned off. The control file specifies the priors for theta, tau, and also specifies burnin, nsample, sampfreq etc. Run bpp to confirm that the control file works. Then run BDdriver as follows: ``` BFdriver A00.ctl 16 tree1.sh ``` Here A00.ctl is the control file we have prepared. $K = 16$ is the number of points in the Gauss-Legendre qradrature algorithm for numerical integration. You can use 8 for testing, and 16 or 32 for real calculation. `tree1.sh` is the temporary script file for job submission using `qsub`. The BFdriver command does a few things. First it reads the control file specified (`A00.ctl`) and creates 16 control files with names like `A00.b01.ctl`, \..., `A00.b16.ctl`. Each of those control files has the same content as `A00.ctl` except that one extra line is inserted at the beginning, like the following ``` BayesFactorBeta = 0.122298 * w=0.124629.ctl ``` This specifies the beta value when the control file is used to run BPP. Second `BFdriver` creates a file named `betaweights.txt`, which lists the beta values and Gauss-Legendre weights. I have copied those values into an excel file in `frogs/BFdriver.frogs.xls`. Third `BFdriver` creates a file named commands, which has the bash shell scripts for submitting the 16 jobs using `qsub`. You use the following to submit the jobs. ``` source commands ``` You can look at the content of `tree1.sh` to see the script for the last job: ``` more tree1.sh ``` which should have the content like the following: ``` bpp --cfile bpp.b16.ctl > log.b16.txt ``` You can use `qstat` to check the status of the 16 jobs you have submitted. When the jobs are running, they generate output files in the current folder, such as `mcmc.b01.txt`, `out.b01.txt`, and `log.b01.txt` (which logs the screen output). After all jobs are finished, you can use grep to extract the line with `BFbeta` from screen log (this command is at the bottom of the file commands). ``` grep BFbeta log.b*.txt ``` Then copy the ElnfX values into the excel file, and estimate the logarithm of the marginal likelihood by summing (`weights * ElnfX / 2`) over the 16 points. This gives the log marginal likelihood to be $-3185.93$. #### Exercise and results Duplicate the calculation for the alternative species tree (tree2) in the folder `/bf2`. Change the species tree topology to (((L, H), C), K), and use tree2.sh as the temporary job script file. Everything else should be the same as described above. My runs gave the log marginal likelihood for tree 2 to be $-3186.14$. The ratio of the posterior probabilities for the two trees is then estimated to be $\frac{P_1}{P_2} = \mathrm{e}^{-3185.93 - (-3186.14)} = \mathrm{e}^{0.21} = 1.24$. With BPP4.0 and the inverse-gamma priors on $\theta$ and $\tau$s, the MCMC runs (A01) gave the posterior probabilities for trees 1 and 2 as 0.167 and 0.137, with the ratio 1.22. The MCMC runs and the marginal likelihood calculations seem to agree with each other. (Note that with BPP3 and the gamma priors on $\theta$ and $\tau$s, the posteriors are 0.16 and 0.13 with the ratio 1.2 ([Yang 2015](https://doi.org/10.1093/czoolo/61.5.854), fig.4). #### Common errors and problems Check the bpp control file by running bpp at the command line before submitting the jobs. Make sure that the bpp program is on your path or use a full path. You may have to edit the source file `BFdriver.c` and recompile. ## Experimental Features (Use at Your Own Risk!) *This section is reserved for experimental features that are under development.* --- ## Troubleshooting This section covers common errors, problems, and their solutions. ### Installation Issues #### "Command not found: bpp" **Problem:** The shell cannot find the BPP executable. **Solutions:** - Ensure BPP is compiled: run `make` in the `src/` directory - Add BPP to your PATH: `export PATH=$PATH:/path/to/bpp/src` - Use the full path: `/path/to/bpp/src/bpp --cfile control.ctl` #### Compilation errors on Mac **Problem:** Compilation fails with errors about missing headers or incompatible flags. **Solutions:** - Ensure Xcode command line tools are installed: `xcode-select --install` - For ARM Macs (M1/M2/M3), use the appropriate Makefile: `make -f Makefile.ARM64` ### Control File Errors #### "Option 'migprior' is obsolete" **Problem:** You are using the old `migprior` option which was replaced in v4.8.0. **Solution:** Replace `migprior` with `wprior`. Note the parameterization has changed: - Old: $M_{XY} = m_{XY} \cdot N_Y$ (population migration rate) - New: $w_{XY} = m_{XY} / \mu$ (mutation-scaled migration rate) Convert using: $w_{XY} = 4 M_{XY} / \theta_Y$ #### "Option 'outfile' and 'mcmcfile' are now obsolete" **Problem:** These options have been replaced by `jobname`. **Solution:** Replace both options with a single `jobname` line: ``` # Old (no longer works): outfile = results.txt mcmcfile = results.mcmc.txt # New: jobname = results ``` #### "Option 'diploid' was renamed to 'phase'" **Problem:** The option name changed but the syntax is the same. **Solution:** Simply rename `diploid` to `phase` in your control file. #### "Erroneous format of 'clock'" **Problem:** Invalid clock specification. **Solution:** Valid clock options are: ``` clock = 1 # strict clock (default) clock = 2 a_vbar b_vbar a_vi prior dist # independent rates clock = 3 a_vbar b_vbar a_vi prior dist # correlated rates clock = 4 a_vbar b_vbar # simple rates ``` ### MCMC Issues #### MCMC not converging / poor mixing **Symptoms:** ESS values are low, trace plots show trends or poor mixing. **Solutions:** 1. **Increase run length:** Increase `burnin` and `nsample` 2. **Adjust step sizes:** Set `finetune = 0` initially to let BPP auto-tune, then examine the acceptance rates in the output 3. **Check priors:** Ensure priors are appropriate for your data scale 4. **Run multiple chains:** Use different seeds (`seed = -1`) and compare results 5. **Simplify the model:** Start with a simpler model (e.g., JC69 instead of GTR) #### Acceptance rates too low or too high **Problem:** MCMC proposals are being rejected too often (< 20%) or accepted too often (> 80%). **Solution:** Adjust the `finetune` parameters: - Low acceptance rate → decrease step size - High acceptance rate → increase step size Let BPP auto-tune first (`finetune = 1`), then manually adjust if needed. #### "NaN" or "Inf" in output **Problem:** Numerical overflow or underflow occurred. **Solutions:** - Check that sequence data is valid (no unusual characters) - Ensure prior parameters are reasonable - Try using `cleandata = 1` to remove ambiguous sites - Check for very short or very long branch lengths in starting tree ### Data Issues #### "Number of sequences does not match" **Problem:** Mismatch between Imap file and sequence file. **Solution:** Ensure every sequence in the sequence file has a corresponding entry in the Imap file, and that sample counts in `species&tree` match. #### "Duplicate taxon or node label" **Problem:** The same name appears more than once in the species tree. **Solution:** Ensure all node labels in your species tree are unique. #### Sequences not aligning with species **Problem:** Wrong sequences assigned to species. **Solution:** Check your Imap file carefully. Format is: ``` sequence_name species_name ``` ### Memory and Performance #### Out of memory errors **Solutions:** - Reduce number of threads: `threads = 1` - Analyze fewer loci: reduce `nloci` - Use a machine with more RAM #### Analysis running very slowly **Solutions:** - Increase threads: `threads = N` (where N = number of CPU cores) - Use `loadbalance = zigzag` for better thread utilization - Use simpler substitution model - Reduce number of loci for initial testing ### Getting Help If you cannot resolve your issue: 1. Check the [BPP GitHub issues](https://github.com/bpp/bpp/issues) 2. Post on the [BPP discussion forum](https://groups.google.com/g/bpp-discussion) 3. Include in your report: - BPP version (`bpp --version`) - Operating system - Complete error message - Minimal control file that reproduces the issue --- ## Common Recipes This section provides step-by-step instructions for common analysis tasks. ### Recipe 1: Basic Parameter Estimation (A00) **Goal:** Estimate divergence times and population sizes for a known species tree. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = A00_analysis speciesdelimitation = 0 # fixed delimitation speciestree = 0 # fixed tree species&tree = 3 A B C 5 5 5 ((A, B), C); usedata = 1 nloci = 100 cleandata = 0 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 finetune = 1 print = 1 0 0 0 0 burnin = 10000 sampfreq = 10 nsample = 100000 ``` **Key points:** - Set `speciesdelimitation = 0` and `speciestree = 0` - Provide your species tree in Newick format - Adjust priors based on expected θ and τ values ### Recipe 2: Species Tree Estimation (A01) **Goal:** Infer the species tree when species assignments are known. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = A01_analysis speciesdelimitation = 0 # fixed delimitation speciestree = 1 # estimate tree species&tree = 4 A B C D 5 5 5 5 ((A, B), (C, D)); # starting tree usedata = 1 nloci = 100 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 finetune = 1 burnin = 20000 sampfreq = 10 nsample = 200000 ``` **Key points:** - Set `speciestree = 1` to estimate the tree topology - The tree in `species&tree` is just a starting tree - Run longer chains (more `nsample`) for tree estimation ### Recipe 3: Species Delimitation with Guide Tree (A10) **Goal:** Test species boundaries using a fixed guide tree. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = A10_analysis speciesdelimitation = 1 1 2 1 # rjMCMC algorithm 1 speciestree = 0 # fixed guide tree species&tree = 4 A B C D 5 5 5 5 ((A, B), (C, D)); usedata = 1 nloci = 100 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 speciesmodelprior = 1 # uniform prior on models finetune = 1 burnin = 20000 sampfreq = 10 nsample = 200000 ``` **Key points:** - `speciesdelimitation = 1` enables delimitation - The additional parameters control the rjMCMC algorithm - Output includes posterior probabilities for different delimitation models ### Recipe 4: Adding Migration (MSC-M) **Goal:** Estimate migration rates between populations. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = MSCM_analysis speciesdelimitation = 0 speciestree = 0 species&tree = 3 A B C 10 10 10 ((A, B)AB, C)R; # label internal nodes! usedata = 1 nloci = 100 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 # Migration setup wprior = 2 200 # gamma prior on migration rates migration = 1 # number of migration bands A B # migration from A to B finetune = 1 burnin = 20000 sampfreq = 10 nsample = 200000 ``` **Key points:** - Label internal nodes in the species tree (e.g., `AB`, `R`) - Use `wprior` (not `migprior`) for BPP v4.8.0+ - Specify migration connections with source and target populations ### Recipe 5: Adding Introgression (MSC-I) **Goal:** Model introgression/hybridization events. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = MSCI_analysis speciesdelimitation = 0 speciestree = 0 species&tree = 3 A B C 10 10 10 ((A, B)AB, C)R; usedata = 1 nloci = 100 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 phiprior = beta 1 1 # prior on introgression probability # Introgression: from branch leading to A into branch leading to C hybridization = R A, AB C as S H tau=yes, no phi=0.1 finetune = 1 burnin = 20000 sampfreq = 10 nsample = 200000 ``` **Key points:** - Use `phiprior` to set prior on introgression probability - `hybridization` defines the introgression event - See MSC-I documentation for full syntax options ### Recipe 6: Relaxed Clock Analysis **Goal:** Allow substitution rate variation among lineages. **Control file:** ``` seed = -1 seqfile = sequences.txt Imapfile = species_map.txt jobname = relaxed_clock speciesdelimitation = 0 speciestree = 0 species&tree = 3 A B C 10 10 10 ((A, B), C); usedata = 1 nloci = 100 thetaprior = invgamma 3 0.002 tauprior = invgamma 3 0.03 # Relaxed clock with independent rates among branches clock = 2 10.0 100.0 5.0 iid G finetune = 1 burnin = 20000 sampfreq = 10 nsample = 200000 ``` **Clock options:** - `clock = 1` — strict clock (default) - `clock = 2 ... iid G` — independent rates, gamma distributed - `clock = 3 ... iid G` — autocorrelated rates - `clock = 4 a b` — simple rates model --- ## Glossary ### Parameters | Symbol | Name | Description | |--------|------|-------------| | θ (theta) | Population size parameter | $\theta = 4N\mu$, where N is effective population size and μ is mutation rate | | τ (tau) | Divergence time | Time of species divergence, measured in expected substitutions per site | | φ (phi) | Introgression probability | Probability that a lineage follows the introgression path (MSC-I model) | | w | Migration rate | Mutation-scaled migration rate: $w_{XY} = m_{XY}/\mu$ (MSC-M model, v4.8.0+) | | M | Population migration rate | $M_{XY} = m_{XY} \cdot N_Y$, expected migrants per generation (pre-v4.8.0) | ### Models | Term | Description | |------|-------------| | MSC | Multispecies Coalescent — the basic model relating gene trees to species trees | | MSC-I | MSC with Introgression — allows hybridization/introgression events | | MSC-M | MSC with Migration — allows continuous gene flow between populations | | IM | Isolation with Migration — another name for the MSC-M model | ### Analysis Types | Code | Name | Description | |------|------|-------------| | A00 | Parameter estimation | Fixed species tree and delimitation, estimate θ and τ | | A01 | Species tree estimation | Fixed delimitation, estimate species tree topology | | A10 | Species delimitation | Fixed guide tree, test alternative delimitations | | A11 | Joint estimation | Simultaneously estimate tree and delimitation | ### MCMC Terms | Term | Description | |------|-------------| | Burn-in | Initial MCMC iterations discarded before sampling | | ESS | Effective Sample Size — measure of independent samples in MCMC output | | Mixing | How well the MCMC explores parameter space | | Convergence | When the MCMC has reached the stationary distribution | | Chain | A single MCMC run | | Proposal | A suggested move in parameter space | | Acceptance rate | Fraction of proposals accepted | ### Priors | Prior | Syntax | Use | |-------|--------|-----| | Inverse-gamma | `invgamma a b` | Conjugate prior for θ; allows analytical integration | | Gamma | `gamma a b` | Alternative prior for θ and τ | | Beta | `beta a b` | Prior for bounded parameters like φ | | Dirichlet | `dir` | Prior for rate variation among loci | ### File Types | File | Description | |------|-------------| | Control file | Main configuration file (`.ctl`) specifying all analysis options | | Sequence file | Phylip-format alignment of DNA/protein sequences | | Imap file | Maps individual sequences to species | | MCMC file | Output file containing posterior samples (`.mcmc.txt`) | | Checkpoint file | Saved state for resuming analysis (`.chk`) | ### Substitution Models | Model | Description | |-------|-------------| | JC69 | Jukes-Cantor — equal rates, equal frequencies | | K80 | Kimura 2-parameter — transition/transversion bias | | F81 | Felsenstein 1981 — unequal base frequencies | | HKY | Hasegawa-Kishino-Yano — K80 + F81 | | GTR | General Time Reversible — most general reversible model | ---