BDBV 2026 outbreak: phylogenetic dating analysis
Vinicius Franceschi and Erik Volz
2026-07-16
Summary
This report analyses genomic data generated by other researchers and originally described in the Virological.org post. There is no intention to publish these analyses. Instead, this report provides complementary analyses to help understand the ongoing Bundibugyo ebolavirus (BDBV) outbreak.
Data download and curation
The sequence alignment and associated metadata were downloaded from Pathoplexus on 9 July 2026. The original post describes 139 sequences sampled between 2 May 2026 and 23 June 2026; however, 149 sequences were available on Pathoplexus at the time of download.
Before analysis, we performed several curation steps, largely following the recommendations in the original post:
- Removed sequence PP_00764QW because of evidence of ADAR editing.
- Trimmed the alignment to 18,900 nt to remove sequencing artefacts at the 3’ end (introduced by sequence PP_006XHL9.3, which was retained).
- Masked the first 65 nt and last 432 nt, as these positions were covered in only a small number of sequences.
To make the analyses as comparable as possible with the original report, we downloaded the corresponding Pathoplexus SeqSet and matched sequence identifiers. This yielded 127 matching sequences, indicating that the datasets are unfortunately not identical. In addition, sequence PP_00711T3 failed the IQ-TREE composition test and was removed prior to downstream analyses.
Alignment and tree reconstruction
Following curation, the Pathoplexus alignment was aligned to the
NC_014373 2007 BDBV reference genome using MAFFT with
the --keeplength and --add options. A
maximum-likelihood (ML) phylogeny was then reconstructed using IQ-TREE 2
with ModelFinder, which selected the K3Pu+F
substitution model (-m MFP -B 1000 --bnni -alrt 1000).
Quality-control checks using Raccoon
(aln-qc and tree-qc) did not identify any
issues requiring further action.
All subsequent analyses were based on this ML tree and included:
- data loading and tree rooting;
- molecular clock dating using multiple models in treedater and BactDating;
- identification and removal of temporal outliers using treedater and DiagnoDating;
- molecular clock model comparison and selection;
- root-to-tip regression;
- estimation of uncertainty in evolutionary rates and node dates;
- estimation of effective population size through time; and
- reconstruction of geographic transitions using maximum parsimony.
1. Load libraries
2. Load tree & metadata
2.1 Tree & alignment settings
set.seed(9876)
TR_DIR <- "results/2026_07_09_5pm"
tr <- read.tree(glue("{TR_DIR}/bdbv.treefile"))
tr##
## Phylogenetic tree with 127 tips and 125 internal nodes.
##
## Tip labels:
## NC_014373, PP_006XHKB.2, PP_00764EH.1, PP_00765UM.1, PP_0076617.1, PP_00765M1.1, ...
## Node labels:
## , 0/72, 75/73, 0/21, 0/4, 0/7, ...
##
## Unrooted; includes branch length(s).2.2 Metadata
DATA_DIR <- "data/2026_07_09_5pm/"
F_NAME <- "ebola-bdbv_metadata_2026-07-09T1550.tsv"
md <- read.csv(glue("{DATA_DIR}/{F_NAME}"), sep = "\t")
#colnames(md)
table(md$geoLocAdmin2)##
## ARU Arua BAMBU Bunia BUNIA Bunya Congo
## 1 1 3 2 14 36 2
## DRC Gety Hoima Kampala Katwa KILO Komanda
## 2 2 1 5 1 2 1
## LITA MAMBASA MANGALA Mongbwalu MONGBWALU MONGWALU MUNGWALU
## 4 1 1 2 1 2 4
## NIZI Not provided Ntoroko NYAKUNDE Nyankunde NYANKUNDE RIMBA
## 9 2 1 2 2 1 1
## Rwampara RWAMPARA SOTA unknown Wakiso
## 11 27 1 1 32.3 Match metadata to tree tips
## [1] 126md_match <- md %>%
filter(accessionVersion %in% tr$tip.label) %>%
dplyr::slice(match(tr$tip.label, accessionVersion))
md_match$sampleCollectionDate <- as.Date(md_match$sampleCollectionDate)
md_match$num_date <- decimal_date(md_match$sampleCollectionDate)
# Named vector of sample times
sts <- md_match$num_date
names(sts) <- md_match$accessionVersion3. Root tree & drop outgroup
4. Initial clock dating (strict)
td0 <- dater(
tr, sts, s = aln_len,
clock = "strict",
meanRateLimits = c(0.0005, 0.002),
numStartConditions = 20, quiet = FALSE,
ncpu = 4
)## Note: Initial guess of substitution rate not provided. Will attempt to guess starting conditions. Provide initial guesses of the rate using *omega0* parameter.
## Note: Minimum temporal branch length (*minblen*) set to 0.000113068058273516. Increase *minblen* in the event of convergence failures.
## Tree is rooted. Not estimating root position.
## Initial guesses of substitution rate: 8.78437382887111e-06,2.10974324671581e-05,3.48291496089972e-05,5.95540891219653e-05,6.35892509370373e-05,8.32221642627345e-05,0.000100560829622948,0.000106040694262712,0.000128178033814052,0.000133670741449473,0.000149774752289877,0.000163939494996624,0.000170951430309609,0.000191814317753028,0.00019204430936788,0.000212459617011739,0.000218460878004584,0.000232976884136016,0.000243214384353623,0.000243538769524812,0.000253451884571231,0.000267544709784504,0.000273969151695508,0.000290688977526408,0.000294614450954219,0.000313142330729711,0.000315477338397638,0.00033504734541521,0.00033665401641737,0.000356526323507719,0.000358250734893195,0.000377687031619072,0.000380388074444535,0.000398627045443617,0.000403206604444512,0.000419437199456908,0.000426874679585282,0.000429820841437541,0.000440204483418174,0.00045159961909825,0.000461014637431465,0.000477644394878376,0.000481954651256011,0.000503115359367364,0.000505353655634688,0.000524594337459873,0.000535196399789737,0.000546499352145372,0.00056784031510974,0.000568952705348674,0.000592096973090578,0.000604291547302476,0.00061610291335027,0.000641180804870498,0.000646185326156106,0.000667597373507202,0.000695702187878458,0.000696481952247756,0.000725970941425609,0.000759080853252135,0.000761537642182704,0.000796052431938045,0.000838544250407924,0.000860276515502149,0.000889558867948174,0.000955543237054781,0.00105569162175096,0.00122945913758186,0.0014301443779574## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0005538747 0 2026.183 -250.3063
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.00061125 0 2026.203 -249.9776
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.0007188471 0 2026.222 -250.2246
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0005945648 0 2026.184 -250.266
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0006781569 0 2026.217 -249.9833
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.0007753669 0 2026.24 -251.465
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0005945648 0 2026.192 -250.0323
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0006781569 0 2026.217 -249.9833
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.0007753669 0 2026.24 -251.465
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0005897992 0 2026.196 -250.0536
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0006781569 0 2026.215 -249.9699
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.0007753669 0 2026.24 -251.465
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0006352549 0 2026.201 -249.9637
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0007188471 0 2026.229 -250.4141
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0006323524 0 2026.208 -249.9291
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0007188471 0 2026.228 -250.3708
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0006339887 0 2026.209 -249.9276
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0007188471 0 2026.228 -250.3945
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0006781569 0 2026.216 -249.9745
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0007753669 0 2026.24 -251.465
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007188471 0 2026.22 -250.2419
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008134118 0 2026.249 -251.9879
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007188471 0 2026.223 -250.2319
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008134118 0 2026.249 -251.9879
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007188471 0 2026.23 -250.5201
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008134118 0 2026.249 -251.9879
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007343893 0 2026.232 -250.6092
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008134118 0 2026.252 -252.4934
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007659935 0 2026.237 -251.1446
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.000854102 0 2026.258 -255.3773
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007917293 0 2026.244 -252.0332
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008947921 0 2026.263 -267.8198
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0007924749 0 2026.244 -252.0618
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008947921 0 2026.263 -267.7991
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0008134118 0 2026.25 -252.1248
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0008947921 0 2026.266 -267.2933
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.000854102 0 2026.257 -254.3846
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.000947133 0 2026.272 -266.4799
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.000854102 0 2026.257 -254.701
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.000947133 0 2026.272 -266.4799
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0008947921 0 2026.263 -267.7014
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.0009984486 0 2026.278 -265.5533
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001072949 0 2026.29 -262.4176
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 3 3 0.001132032 0 2026.297 -262.493
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0009261882 0 2026.27 -266.764
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001032259 0 2026.282 -262.6961
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001138228 0 2026.293 -262.4368
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 3 3 0.001222884 0 2026.302 -262.8114
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0009552358 0 2026.273 -266.3751
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001032259 0 2026.285 -262.5424
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001138228 0 2026.293 -262.4368
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 3 3 0.001222884 0 2026.302 -262.8114
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.0009921389 0 2026.277 -265.8954
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001072949 0 2026.289 -262.422
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001132032 0 2026.297 -262.493
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.001032259 0 2026.285 -262.5416
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001138228 0 2026.293 -262.4368
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001222884 0 2026.302 -262.8114
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.001117727 0 2026.295 -262.4479
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001207773 0 2026.301 -262.6852
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.001300926 0 2026.309 -263.761
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001321338 0 2026.313 -265.1138
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 0 0 0.001427051 0 2026.32 -268.9745
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 1 1 0.001427051 0 2026.32 -268.9066
## ---
## iteration median_unadjusted_rate coef_of_var tmrca logLik
## 2 2 0.001427051 0 2026.32 -268.9066
## ---5. Outlier Removal
## [1] taxon q p loglik rates branch.length
## <0 rows> (or 0-length row.names)td1 <- td0
tr1 <- tr
while (length(ot) > 0) {
tr1 <- drop.tip(tr1, ot)
td1 <- dater(
tr1, sts, s = aln_len,
clock = "strict",
meanRateLimits = c(0.0005, 0.002),
numStartConditions = 20, quiet = TRUE, #FALSE
ncpu = 4
)
ot <- .otips(td1)
}
# Tips removed as outliers (empty = none removed)
print(setdiff(td0$tip.label, td1$tip.label))## character(0)6. Clock model selection
6.1 Relaxed clock test
rct <- relaxedClockTest(
tr1, sts, aln_len,
meanRateLimits = c(0.0005, 0.002),
quiet = FALSE, ncpu = 8, nreps=5000,
overrideTempConstraint = FALSE
) Best clock model: strict
Null distribution of rate coefficient of variation: 0.000571629763668337 0.000823308918992667
Returning best treedater fit6.3 Inspect best-fit model (strict)
##
## Phylogenetic tree with 126 tips and 125 internal nodes.
##
## Tip labels:
## PP_006XHKB.2, PP_00764EH.1, PP_00765UM.1, PP_0076617.1, PP_00765M1.1, PP_0076SYS.1, ...
## Node labels:
## , 0/72, 75/73, 0/21, 0/4, 0/7, ...
##
## Rooted; includes branch length(s).
##
## Time of common ancestor
## 2026.20852968182
##
## Time to common ancestor (before most recent sample)
## 0.265442920918304
##
## Weighted mean substitution rate (adjusted by branch lengths)
## 0.000633988705841585
##
## Unadjusted mean substitution rate
## 0.000633988705841585
##
## Clock model
## strict
##
## Coefficient of variation of rates
## 0## [1] "2026-03-18 02:43:12 UTC"
7. Bootstrap for confidence intervals (from iqtree ufboot distribution)
# Bootstrap
btd <- treedater::boot(
td1, btres,
ncpu = 8,
overrideTempConstraint = FALSE, quiet = TRUE
)## Running in quiet mode. To print progress, set quiet=FALSE.## pseudo ML 2.5 % 97.5 %
## Time of common ancestor 2.026209e+03 2.026089e+03 2.026291e+03
## Mean substitution rate 6.339887e-04 3.357899e-04 1.343159e-03
##
## For more detailed output, $trees provides a list of each fit to each simulation## 2.5% 97.5%
## "2026-02-02 13:02:11 UTC" "2026-04-17 01:34:46 UTC"8. BactDating
trbd1 <- tr1
trbd1$edge.length <- trbd1$edge.length * aln_len
bd1a <- bactdate(
trbd1, sts[tr1$tip.label],
initMu = aln_len * td1$mean.rate,
model = "arc", showProgress = FALSE # TRUE
)
bd1a## Phylogenetic tree dated using BactDating
##
## Phylogenetic tree with 126 tips and 125 internal nodes.
##
## Tip labels:
## PP_006XHKB.2, PP_00764EH.1, PP_00765UM.1, PP_0076617.1, PP_00765M1.1, PP_0076SYS.1, ...
## Node labels:
## , 0/72, 75/73, 0/21, 0/4, 0/7, ...
##
## Rooted; includes branch length(s).
## Probability of root branch=0.49
## likelihood=-2.62e+02 [-2.75e+02;-2.49e+02]
## prior=-2.28e+02 [-2.59e+02;-2.00e+02]
## mu=1.29e+01 [9.45e+00;1.67e+01]
## sigma=3.09e-01 [1.66e-02;6.20e-01]
## alpha=2.31e+00 [1.72e+00;3.10e+00]
## Root date=2025.86 [2025.63;2026.06]
## Root date for most likely root=2025.86 [2025.64;2026.05]Note: The estimated mean rate from BactDating is consistent with treedater, but the inferred root date is unusually early.
9. DiagnoDating
9.1 Root-to-tip (DiagnoDating)
NEEDS DiagnoDating >= 0.9.1, ie the fix-treedater-units branch (PR #1), which is where the additive clock, the seqlen argument and matching of named dates were fixed. Before that, treedater results always came back as model=‘poisson’ with relax=0:
devtools::install_github('xavierdidelot/DiagnoDating', ref = 'fix-treedater-units')
## $rate
## date
## 0.0002432144
##
## $ori
## (Intercept)
## 2025.296
##
## $pvalue
## [1] 0.2215DiagnoDating, like BactDating, wants branch lengths in substitutions rather than per site, so rescale the tree. seqlen then tells treedater how to convert back, and omega0 and meanRateLimits stay in treedater’s own units of substitutions per site per year.
9.2 Run & diagnose
ddtd <- runDating( trdd, sts[ tr1$tip.label ], algo = 'treedater', keepRoot = TRUE,
seqlen = aln_len, clock = 'strict', omega0 = td1$mean.rate, meanRateLimits = td1$mean.rate*c(1,1+1e-6) ) ### Result from treedater, model poisson, clock rate 11.67, relaxation parameter 0.00, root date 2026.23
##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: Normal distribution
## with parameters mean = 0.000, sd = 1.000
##
## data: n
## An = 2.4185, p-value = 0.0547310. Profile likelihood for the clock rate
ratetologlik <- function(omega){
capture.output( x <- sapply( btres, function(tr) dater(tr, sts, s=aln_len, clock="strict", omega0 = omega, meanRateLimits = omega*c(.98,1.02), quiet=TRUE, ncpu=1)$loglik ) )
x |> pomp::logmeanexp()
}
omegas <- seq( 0.0001, 0.002, length.out = 30)
llprof <- unlist( parallel::mclapply( omegas, ratetologlik, mc.cores=8 ) )
source( 'scripts/profplot-02b.R')
11. Effective population size (mlesky)
library(mlesky)
range_tr <- max(sts) - td1$timeOfMRCA
range_tr_weeks <- ceiling(range_tr * 52.1775)
range_tr_weeks## [1] 14mlsk <- mlskygrid(td1, sampleTimes=sts, res=NULL, tau=NULL, tau_lower=0.001, tau_upper=100,
ncross=10, ncpu=4, quiet=F, model=2)
plot(mlsk, logy=T)
mlsk_pboot <- mlesky::parboot(mlsk, nrep=1000, ncpu=8, dd=F)
# regular bootstrap NEEDS mlesky >= 0.1.9: # devtools::install_github('emvolz-phylodynamics/mlesky')
# bmlsk <- mlesky::boot( mlsk, btd$trees, ncpu = 8, sampleTimes = sts )
plot(mlsk_pboot, logy=T, ggplot=T)
12. Ancestral state estimation with tree figure
library(phangorn)
collapse_dated <- function(td, aln_len, min_subs = 0.5) {
phy <- as.phylo(td)
subs <- td$intree$edge.length * aln_len # substitutions per edge (aligned to phy$edge)
nd <- node.depth.edgelength(phy) # node heights (time) from root
tmp <- phy; tmp$edge.length <- subs
topo <- di2multi(tmp, tol = min_subs) # collapse by substitutions
Ntt <- Ntip(topo); dts <- numeric(Ntt + topo$Nnode) # re-time from the original treedater dates
dts[1:Ntt] <- (td$timeOfMRCA + nd[1:Ntip(phy)])[match(topo$tip.label, phy$tip.label)]
desc <- Descendants(topo, (Ntt+1):(Ntt+topo$Nnode), "tips")
for (i in seq_along(desc)) {
tips <- topo$tip.label[desc[[i]]]
mn <- if (length(tips)==Ntip(phy)) Ntip(phy)+1 else getMRCA(phy, tips)
dts[Ntt+i] <- td$timeOfMRCA + nd[mn]
}
topo$edge.length <- dts[topo$edge[,2]] - dts[topo$edge[,1]]
topo$root.time <- td$timeOfMRCA
topo
}
.norm <- function(x){ x <- toupper(trimws(x)); x[x=="" | x=="NOT PROVIDED"] <- "UNKNOWN"
dplyr::recode(x, "MONGWALU"="MONGBWALU","MUNGWALU"="MONGBWALU","NYAKUNDE"="NYANKUNDE") }
# shared locality set + palette (both trees have identical tips, so colours are comparable)
loc0 <- .norm(md_match$geoLocAdmin2[ match(as.phylo(td1)$tip.label, md_match$accessionVersion) ])
top5 <- names(sort(table(loc0), decreasing=TRUE))[1:5]
LEV <- c(top5, "Other")
pal <- setNames(c("#0072B2","#D55E00","#009E73","#E69F00","#CC79A7","grey72"), LEV)
mrsd <- as.Date(date_decimal(max(sts))); ci <- btd$timeOfMRCA_CI; break_gr <- 2026.35
mb <- decimal_date(seq(as.Date("2026-01-01"), as.Date("2026-08-01"), by="month"))
ml <- format(seq(as.Date("2026-01-01"), as.Date("2026-08-01"), by="month"), "%b %Y")
geo_tree_fig <- function(tree, root_date, title) {
phy <- as.phylo(tree); nt <- Ntip(phy)
loc <- .norm(md_match$geoLocAdmin2[ match(phy$tip.label, md_match$accessionVersion) ])
loc[!loc %in% top5] <- "Other"
dat <- phangorn::phyDat(as.matrix(data.frame(s=loc, row.names=phy$tip.label)), type="USER", levels=LEV)
mat <- do.call(rbind, lapply(phangorn::ancestral.pars(phy, dat, "MPR"), function(z) z[1,])); colnames(mat) <- LEV
map <- LEV[apply(mat,1,which.max)]; conf <- apply(mat,1,max)
ntrans <- sum(map[phy$edge[,1]] != map[phy$edge[,2]])
lab <- rep(NA_character_, nt+phy$Nnode) # UFBoot labels, only on collapsed trees
if (!is.binary(phy)) {
ufb <- as.numeric(sub(".*/", "", tr1$node.label))
for (n in (nt+1):(nt+phy$Nnode)) {
tips <- phy$tip.label[ unlist(phangorn::Descendants(phy, n, "tips")) ]
mn <- if (length(tips)==nt) Ntip(tr1)+1 else getMRCA(tr1, tips)
if (!is.null(mn) && !is.na(mn)) lab[n] <- as.character(round(ufb[mn - Ntip(tr1)]))
}
}
nodedf <- data.frame(node=1:(nt+phy$Nnode), state=factor(map,levels=LEV), conf=conf, lab=lab)
pg <- ggtree(tree, mrsd=mrsd, size=0.5) %<+% nodedf + aes(color=state)
# ggtree rounds mrsd to whole days, shifting its x-axis; off re-aligns absolute-date annotations
rootx <- pg$data$x[pg$data$node == nt+1]; rooty <- pg$data$y[pg$data$node == nt+1]
off <- rootx - root_date; cih <- 0.4 # half-height of the thin TMRCA CI band, in tip-row units
pg$layers <- c(annotate("rect", xmin=ci[1]+off, xmax=ci[2]+off, ymin=rooty-cih, ymax=rooty+cih, fill="grey30", alpha=0.45), pg$layers) # thin TMRCA 95% CI band at the root, behind the tree
pg +
geom_vline(xintercept=break_gr+off, linetype="dashed", color="grey55", linewidth=0.4) +
# annotate("label", x=break_gr+off, y=nt*0.965, label="growth-rate\nchange 2026.35", size=2.7,
# color="grey40", hjust=-0.04, vjust=1, lineheight=0.9, fill="white", alpha=0.7) +
annotate("point", x=rootx, y=rooty, shape=18, size=2.8, color="grey20") +
annotate("text", x=mean(ci)+off, y=rooty+cih+nt*0.02, label="TMRCA 95% CI", size=2.6, color="grey35", vjust=0) +
geom_nodepoint(aes(color=state, size=conf), alpha=0.85) +
geom_nodelab(aes(label=lab), geom="label", size=1.9, vjust=-0.6, color="grey20",
fill="white", label.size=0, label.padding=unit(0.04,"lines"), na.rm=TRUE) +
geom_tippoint(aes(fill=state), shape=21, size=2, stroke=0.25, color="white") +
scale_color_manual(values=pal, guide="none") +
scale_fill_manual(values=pal, name="geoLocAdmin2") +
scale_size_continuous(range=c(0.6,3), guide="none") +
scale_x_continuous(breaks=mb+off, labels=ml, limits=c(min(ci[1]+off, rootx)-0.02, max(sts)+off+0.02)) +
labs(title=title,
subtitle=glue("{nt} genomes | {phy$Nnode} internal nodes | MRCA {format(as.Date(date_decimal(root_date)),'%d %b %Y')} ")) +
theme_tree2() +
theme(plot.title=element_text(face="bold",size=12), plot.subtitle=element_text(size=8.5,color="grey40"),
axis.text.x=element_text(size=9,color="grey40"), legend.position=c(0.12,0.75),
legend.background=element_rect(fill="white",color="grey85"), legend.key.size=unit(0.9,"lines"))
}
td1_collapsed <- collapse_dated(td1, aln_len, min_subs = 0.5)
pg_coll <- geo_tree_fig(td1_collapsed, td1_collapsed$root.time, "")
print(pg_coll)