Detecting evolutionary change-points with branch-specific substitution models and shrinkage priors
Data files
Nov 07, 2026 version files 376.51 MB
-
apobec_ESS.pdf
5.96 KB
-
apobec_HMC_save
82.57 KB
-
apobec_hmc_timer.txt
25 B
-
apobec_MLE_analytic_printScreen.txt
204.20 KB
-
apobec_MLE_numeric_printScreen.txt
250.46 KB
-
apobec_Serial_save
76.73 KB
-
apobec_serial_timer.txt
25 B
-
apobec.shrinkage_HMC.xml
27.29 MB
-
apobec.shrinkage_MLE_analytic.xml
27.29 MB
-
apobec.shrinkage_MLE_numeric.xml
27.29 MB
-
apobec.shrinkage_Serial.xml
27.29 MB
-
apobec.shrinkage_simulation_0.4.xml
27.30 MB
-
apobec.shrinkage_simulation_0.6.xml
27.30 MB
-
apobec.shrinkage_simulation_0.8.xml
27.30 MB
-
apobec.shrinkage_simulation_1.0.xml
27.30 MB
-
apobec.shrinkage.xml
27.31 MB
-
apobecMMM.HKY.shrinkage.large.parameters.hmc.log
49.90 MB
-
apobecMMM.HKY.shrinkage.large.parameters.serial.log
66.51 MB
-
BRCA1_coding_shrinkage_fixedRates.large.parameters.HMC.log
6.10 MB
-
BRCA1_coding_shrinkage_fixedRates.large.parameters.Serial.log
7.31 MB
-
BRCA1_ESS.pdf
5.10 KB
-
BRCA1_hmc_save
11.32 KB
-
BRCA1_hmc_timer.txt
25 B
-
BRCA1_MLE_analytic_printScreen.txt
36.52 KB
-
BRCA1_MLE_numeric_printScreen.txt
20.50 KB
-
BRCA1_mse.py
15.57 KB
-
BRCA1_serial_save
11.03 KB
-
BRCA1_serial_timer.txt
26 B
-
BRCA1_shrinkage_fixedRates_HMC.xml
37.18 KB
-
BRCA1_shrinkage_fixedRates_MLE_analytic.xml
37.09 KB
-
BRCA1_shrinkage_fixedRates_MLE_numeric.xml
37.31 KB
-
BRCA1_shrinkage_fixedRates_MLE.xml
36.98 KB
-
BRCA1_shrinkage_fixedRates_Serial.xml
36.80 KB
-
BRCA1_shrinkage_fixedRates.xml
40.51 KB
-
ESS.Rmd
8.10 KB
-
mse.py
5.03 KB
-
README.md
22.74 KB
Abstract
Branch-specific substitution models are popular for detecting evolutionary change-points, such as shifts in selective pressure. However, applying such models typically requires prior knowledge of change-point locations on the phylogeny or faces scalability issues with large data sets. To address both limitations, we integrate branch-specific substitution models with shrinkage priors to automatically identify change-points without prior knowledge, while simultaneously estimating distinct substitution parameters for each branch. To enable tractable inference under this high-dimensional model, we develop an analytical gradient algorithm for the branch-specific substitution parameters where the computational time is linear in the number of parameters. We apply this gradient algorithm to infer selection pressure dynamics in the evolution of the BRCA1 gene in primates and mutational dynamics in viral sequences from the recent mpox epidemic. Our novel algorithm enhances inference efficiency, achieving up to a 126-fold speedup per iteration in maximum likelihood optimization when compared to central difference numerical gradient method and up to a 2026-fold improvement in computational performance within a Bayesian framework using Hamiltonian Monte Carlo sampler compared to conventional univariate random walk sampler.
This dataset contains the instructions and files needed to reproduce the analyses in the paper "Detecting Evolutionary Change-Points with Branch-Specific Substitution Models and Shrinkage Priors" by Ji et al.
Description of the data and file structure
The analyses use two empirical data examples:
- BRCA1: a codon alignment of the BRCA1 gene from 8 primate and other mammal taxa.
- MPXV: an alignment of monkeypox virus (MPXV) genomes, used to study APOBEC3-type mutations. File names for this example start with
apobec.
Every analysis is a BEAST XML file run with the hmc-clock branch of BEAST (see "Software setup" below). The *_save files store a saved MCMC state. BEAST loads them with -load_state, so the benchmarking and optimization runs start from the same converged state.
The repository is organized as follows:
docker/ Dockerfile used to build the BEAST + BEAGLE image
xmls/BRCA1/ BEAST XMLs and saved states for the BRCA1 example
xmls/MPXV/ BEAST XMLs and saved states for the MPXV example
xmls/Simulation/ Simulation and bootstrap XMLs and Python 3 scripts
logFiles/ Output from the benchmarking and optimization runs, the R script that summarizes them, and the resulting figures
Files and variables
BEAST input files (xmls/BRCA1/ and xmls/MPXV/)
| File | Description |
|---|---|
BRCA1_shrinkage_fixedRates.xml |
Full Bayesian analysis of BRCA1 with branch-specific substitution parameters under a shrinkage prior. |
BRCA1_shrinkage_fixedRates_HMC.xml |
BRCA1 benchmarking run that samples the branch-specific parameters with Hamiltonian Monte Carlo (HMC). |
BRCA1_shrinkage_fixedRates_Serial.xml |
BRCA1 benchmarking run that uses univariable (serial) random-walk operators. |
BRCA1_shrinkage_fixedRates_MLE_analytic.xml |
Maximum-likelihood estimation (MLE) for BRCA1 with analytic gradients. |
BRCA1_shrinkage_fixedRates_MLE_numeric.xml |
MLE for BRCA1 with numerical (finite-difference) gradients. |
BRCA1_hmc_save, BRCA1_serial_save |
Saved BEAST states loaded with -load_state by the BRCA1 benchmarking and MLE runs. |
apobec.shrinkage.xml |
Full Bayesian analysis of the MPXV example. |
apobec.shrinkage_HMC.xml |
MPXV benchmarking run with the HMC sampler. |
apobec.shrinkage_Serial.xml |
MPXV benchmarking run with univariable (serial) operators. |
apobec.shrinkage_MLE_analytic.xml |
MLE for MPXV with analytic gradients. |
apobec.shrinkage_MLE_numeric.xml |
MLE for MPXV with numerical gradients. |
apobec_HMC_save, apobec_Serial_save |
Saved BEAST states for the MPXV benchmarking and MLE runs. |
Simulation and bootstrap files (xmls/Simulation/)
| File | Description |
|---|---|
apobec.shrinkage_simulation_0.4.xml, apobec.shrinkage_simulation_0.6.xml, apobec.shrinkage_simulation_0.8.xml, apobec.shrinkage_simulation_1.0.xml |
MPXV simulation XMLs. Each one simulates sequences on the MPXV tree under a branch-specific model and then estimates the branch-specific rates by MLE. The suffix is the true APOBEC rate on the branches with elevated rates (the other branches have a rate of 0.1). In the repository these files are named apobec.shrinkage_simulation.xml and sit in subfolders 0.4/, 0.6/, 0.8/ and 1.0/, which is where mse.py expects them. |
mse.py |
Python 3 script for the MPXV simulation study. For each true rate (0.4, 0.6, 0.8, 1.0) it runs BEAST on the matching simulation XML with seeds 1 to 100. It parses the MLE rates and writes all_rates.txt (estimated rates, one row per replicate) and failed.txt (replicates that failed) into each rate subfolder. These outputs are used to compute the mean squared error (MSE) of the estimates. Requires numpy. |
BRCA1_shrinkage_fixedRates_MLE.xml |
Template BRCA1 MLE XML used by BRCA1_mse.py. It contains the original BRCA1 alignment. |
BRCA1_mse.py |
Python 3 script for the BRCA1 bootstrap study. For seeds 1 to 100 it resamples codon sites from the alignment in BRCA1_shrinkage_fixedRates_MLE.xml with replacement. It writes one analytic-gradient and one numeric-gradient MLE XML per replicate to ./xmls/ and runs both in BEAST (loading BRCA1_hmc_save). The log-likelihoods and screen output go to ./logs/. |
BRCA1_hmc_save |
Copy of the saved BRCA1 state, used by BRCA1_mse.py. |
Both Python scripts call BEAST through a hard-coded path, /Users/.../beast-mcmc/build/dist/beast.jar. Change beast_cmd (in mse.py) or beast_analytic_cmd/beast_numeric_cmd (in BRCA1_mse.py) to point to your own beast.jar before running them.
Benchmarking output (logFiles/)
These files compare the HMC sampler with univariable (serial) operators by effective sample size (ESS) per minute.
| File | Description |
|---|---|
BRCA1_coding_shrinkage_fixedRates.large.parameters.HMC.log |
BEAST log from BRCA1_shrinkage_fixedRates_HMC.xml. Tab-delimited, with 3 comment lines starting with # and then a header row. Each row is one MCMC sample (state = iteration number) and holds the joint posterior, prior and likelihood plus the sampled branch-specific parameter values. Columns 1–18 are summary and global quantities, and the remaining columns are the branch-specific parameters used in the ESS comparison. |
BRCA1_coding_shrinkage_fixedRates.large.parameters.Serial.log |
Same as above for BRCA1_shrinkage_fixedRates_Serial.xml. |
apobecMMM.HKY.shrinkage.large.parameters.hmc.log |
BEAST log from apobec.shrinkage_HMC.xml. Same format as above. The first 278 columns are summary and global quantities, and the remaining columns are the branch-specific parameters used in the ESS comparison. |
apobecMMM.HKY.shrinkage.large.parameters.serial.log |
Same as above for apobec.shrinkage_Serial.xml. |
BRCA1_hmc_timer.txt, BRCA1_serial_timer.txt |
Total wall-clock run time of the BRCA1 HMC and serial runs. Each file has one line with a number and a unit (minutes or hours), e.g. 4.14 hours. |
apobec_hmc_timer.txt, apobec_serial_timer.txt |
Total wall-clock run time of the MPXV HMC and serial runs, in the same format. |
ESS.Rmd |
R Markdown script that reads the four .log files and four *_timer.txt files. It drops the first 10% of samples as burn-in and computes the ESS of each branch-specific parameter with the coda package. It divides by 90% of the run time to get ESS per minute, prints minimum/median/maximum ESS per minute and the HMC speed-up, and draws the ESS histograms. Requires the R packages coda, ggplot2, scales and stringr. |
BRCA1_ESS.pdf |
Histogram of log2 ESS per minute for the branch-specific parameters, HMC vs. serial, BRCA1 example (made by ESS.Rmd). |
apobec_ESS.pdf |
Same histogram for the MPXV example. |
Optimization output (logFiles/)
| File | Description |
|---|---|
BRCA1_MLE_analytic_printScreen.txt |
Screen output from running BRCA1_shrinkage_fixedRates_MLE_analytic.xml. It shows the optimizer's progress at each iteration (fx = negative log-likelihood, xnorm, gnorm = gradient norm, step). At the end it lists the final estimates on the original scale (Untransformed X) and the transformed scale (X), plus the final gradient, the gradient type, the final objective value (Fx), the run time (Time) and the number of likelihood/gradient evaluations (Count). |
BRCA1_MLE_numeric_printScreen.txt |
Same output from BRCA1_shrinkage_fixedRates_MLE_numeric.xml (numerical gradients). |
apobec_MLE_analytic_printScreen.txt |
Same output from apobec.shrinkage_MLE_analytic.xml. |
apobec_MLE_numeric_printScreen.txt |
Same output from apobec.shrinkage_MLE_numeric.xml. |
Comparing Time and Count between the analytic and numeric files gives the speed-up from analytic gradients reported in the manuscript.
Sharing/Access information
The source code for the software used is publicly available:
- BEAST (
hmc-clockbranch, commit9899f203d8d43fab1b4242420c0f6f508a088aec): https://github.com/beast-dev/beast-mcmc - BEAGLE (
hmc-clockbranch, commitedfb106eb12efae798945ed3dd0a2f918e8c1a28): https://github.com/beagle-dev/beagle-lib - Pre-built Docker images:
beast1/beastx:x86_64andbeast1/beastx:arm_64on Docker Hub
Software setup
We provide two ways to set up the software:
- Docker image. See the
Docker imagesection if you want to run BEAST inside a local Docker container that won't affect your other BEAST installations. - Global installation. See the
Global installationsection if you want to install BEAST and BEAGLE system-wide. This could interfere with your other BEAST installations.
Each installation method has its own Reproducing the analyses section.
Docker image
We provide a Docker image to make software setup easier. You can either use the pre-built image or build it from the Dockerfile.
Pull the pre-built Docker image
For x86_64 systems (e.g., Intel/AMD chips), pull the x86 image:
docker pull beast1/beastx:x86_64
and start an interactive bash session:
cd where_this_repository_is_stored
docker run --mount type=bind,src="$(pwd)",target=/tmp -it beast1/beastx:x86_64 bash
For arm_64 systems (e.g., Apple silicon chips), pull the arm image (e.g., in a Docker Desktop command line environment):
docker pull beast1/beastx:arm_64
and start an interactive bash session:
cd where_this_repository_is_stored
docker run --mount type=bind,src="$(pwd)",target=/tmp -it beast1/beastx:arm_64 bash
You can now skip ahead to "Reproducing the analyses (Docker)" (the next section explains how to build the images yourself).
Build the image from the Dockerfile
To build the image for x86_64 systems (e.g., Intel/AMD chips), run (e.g., in a terminal window):
cd where_this_repository_is_stored
docker build --platform x86_64 -t beast docker
To build the image for arm_64 systems (e.g., Apple silicon chips), run (e.g., in a Docker Desktop command line environment):
cd where_this_repository_is_stored
docker build --platform linux/arm64 -t beast docker
Then open an interactive bash session in the container:
docker run --mount type=bind,src="$(pwd)",target=/tmp -it beast bash
Reproducing the analyses (Docker)
Run each of the commands below in the interactive bash shell to reproduce the two data examples in the manuscript.
BRCA1
-
Benchmarking
-
HMC (writes
BRCA1_coding_shrinkage_fixedRates.large.parameters.HMC.log)beast -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_HMC.xml -
Univariable (writes
BRCA1_coding_shrinkage_fixedRates.large.parameters.Serial.log)beast -load_state xmls/BRCA1/BRCA1_serial_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_Serial.xml
-
-
Optimization
-
Analytic (save the screen output as
BRCA1_MLE_analytic_printScreen.txt)beast -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_MLE_analytic.xml > BRCA1_MLE_analytic_printScreen.txt -
Numeric (save the screen output as
BRCA1_MLE_numeric_printScreen.txt)beast -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_MLE_numeric.xml > BRCA1_MLE_numeric_printScreen.txt
-
-
Full analysis
beast -seed 666 -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates.xml
MPXV
-
Benchmarking
-
HMC (writes
apobecMMM.HKY.shrinkage.large.parameters.hmc.log)beast -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_HMC.xml -
Univariable (writes
apobecMMM.HKY.shrinkage.large.parameters.serial.log)beast -load_state xmls/MPXV/apobec_Serial_save -overwrite xmls/MPXV/apobec.shrinkage_Serial.xml
-
-
Optimization
-
Analytic
beast -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_MLE_analytic.xml > apobec_MLE_analytic_printScreen.txt -
Numeric
beast -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_MLE_numeric.xml > apobec_MLE_numeric_printScreen.txt
-
-
Full analysis
beast -seed 666 -overwrite xmls/MPXV/apobec.shrinkage.xml
Global installation
Setting up BEAGLE
Follow the BEAGLE installation instructions, but check out the hmc-clock branch.
For Mac users, the following commands compile the CPU version of BEAGLE. See the instructions if you need to install any other required software.
xcode-select --install
brew install libtool autoconf automake
git clone -b hmc-clock https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib
git checkout edfb106eb12efae798945ed3dd0a2f918e8c1a28
mkdir build
cd build
cmake -DBUILD_CUDA=OFF -DBUILD_OPENCL=OFF ..
sudo make install
For Linux users, the commands are similar:
sudo apt-get install build-essential autoconf cmake automake libtool git pkg-config openjdk-9-jdk
git clone -b hmc-clock https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib
git checkout edfb106eb12efae798945ed3dd0a2f918e8c1a28
mkdir build
cd build
cmake -DBUILD_CUDA=OFF -DBUILD_OPENCL=OFF ..
sudo make install
The libraries are installed in /usr/local/lib. You can list them with ls /usr/local/lib/*beagle*.
Setting up BEAST
The following commands compile the hmc-clock branch of BEAST:
git clone -b hmc-clock https://github.com/beast-dev/beast-mcmc.git
cd beast-mcmc
git checkout 9899f203d8d43fab1b4242420c0f6f508a088aec
ant
Mac users may need to install ant with brew install ant through Homebrew.
Linux users can install ant with sudo apt-get install ant.
This compiles the jar files into beast-mcmc/build/dist/, where you will find beast.jar, beauti.jar and trace.jar.
Reproducing the analyses (global installation)
Use the following commands for each of the data examples in the manuscript.
First, change your working directory to where you want to save the output log files:
cd where_you_want_to_save_results
In the commands below, replace xmls/... with the path to the xmls folder of this dataset if it is not inside your working directory.
BRCA1
-
Benchmarking
-
HMC
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_HMC.xml -
Univariable
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/BRCA1/BRCA1_serial_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_Serial.xml
-
-
Optimization
-
Analytic
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_MLE_analytic.xml > BRCA1_MLE_analytic_printScreen.txt -
Numeric
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/BRCA1/BRCA1_hmc_save -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates_MLE_numeric.xml > BRCA1_MLE_numeric_printScreen.txt
-
-
Full analysis
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -seed 666 -overwrite xmls/BRCA1/BRCA1_shrinkage_fixedRates.xml
MPXV
-
Benchmarking
-
HMC
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_HMC.xml -
Univariable
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/MPXV/apobec_Serial_save -overwrite xmls/MPXV/apobec.shrinkage_Serial.xml
-
-
Optimization
-
Analytic
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_MLE_analytic.xml > apobec_MLE_analytic_printScreen.txt -
Numeric
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -force_resume -load_state xmls/MPXV/apobec_HMC_save -overwrite xmls/MPXV/apobec.shrinkage_MLE_numeric.xml > apobec_MLE_numeric_printScreen.txt
-
-
Full analysis
java -jar -Djava.library.path=/usr/local/lib where_beast_is_git_cloned/beast-mcmc/build/dist/beast.jar -seed 666 -overwrite xmls/MPXV/apobec.shrinkage.xml
ESS comparison (HMC vs. univariable)
After the four benchmarking runs finish, put the four .log files and the four *_timer.txt files in the same folder as ESS.Rmd (all of them are already in logFiles/). Then knit the script in R:
cd logFiles
Rscript -e 'rmarkdown::render("ESS.Rmd")'
This prints the ESS-per-minute summaries and speed-ups, and writes BRCA1_ESS.pdf and apobec_ESS.pdf.
Simulation and bootstrap studies
These scripts call BEAST directly, so run them with the global installation (or inside the Docker container after you change the beast.jar path in the scripts). Edit the hard-coded beast.jar path in each script first.
MPXV simulation (MSE of the estimated APOBEC rates):
cd xmls/Simulation
python3 mse.py
This expects 0.4/, 0.6/, 0.8/ and 1.0/ subfolders, each holding the matching simulation XML named apobec.shrinkage_simulation.xml. If you downloaded the flat files from Dryad, set them up first:
for r in 0.4 0.6 0.8 1.0; do mkdir -p $r; cp apobec.shrinkage_simulation_$r.xml $r/apobec.shrinkage_simulation.xml; done
BRCA1 bootstrap (analytic vs. numeric MLE on 100 resampled alignments):
cd xmls/Simulation
python3 BRCA1_mse.py
This needs BRCA1_shrinkage_fixedRates_MLE.xml and BRCA1_hmc_save in the working directory.
