start_tree <- readTrees("three_taxon_tree.tre")[1] n_taxa <- 3 taxa <- start_tree.taxa() n_branches <- n_taxa*2 - 2 mvi = 0 mni = 0 ############################################################################################################# Q := fnJC(4) ############################################################################################################# root_time <- 1 rho <- 1 lambda ~ dnExponential(10) moves[++mvi] = mvScale(lambda,lambda=1,weight=3) timetree ~ dnBDP(lambda, mu=0, rho=rho, rootAge=root_time, samplingStrategy="uniform", condition="nTaxa", taxa=taxa) timetree.setValue(start_tree) moves[++mvi] = mvNodeTimeSlideUniform(timetree, weight=30.0) moves[++mvi] = mvTreeScale(timetree, delta = 1, weight=3.0) A_B <- clade("A", "B") unknown_time := tmrca(timetree, clade("A", "B")) ############################################################################################################# ucln_log_mean <- -3.00667 for (i in 1:n_branches){ branch_rates[i] ~ dnLnorm(ucln_log_mean, 0.587405/4) moves[++mvi] = mvScale(branch_rates[i], lambda = 1, weight=5) } moves[++mvi] = mvVectorScale(branch_rates, lambda = 1.0, weight = 2.0) moves[++mvi] = mvVectorSingleElementScale(branch_rates, lambda=30.0, weight = 1.0) ############################################################################################################# seq ~ dnPhyloCTMC(timetree, Q=Q, branchRates=branch_rates, type="DNA") seq.clamp(data) ############################################################################################################# mymodel = model(timetree) monitors[++mni] = mnModel(filename = results[z] + concat[z] + "standard_relaxed_sampling.log", printgen=(5000 + (scale_number[z]*4000))/1000, separator = TAB) monitors[++mni] = mnFile(filename= results[z] + concat[z] + "standard_relaxed_trees.trees", printgen=(5000 + (scale_number[z]*4000))/1000, timetree) monitors[++mni] = mnScreen(printgen=500) mymcmc = mcmc(mymodel, monitors, moves, nruns=1) mymcmc.run(generations=5000 + (scale_number[z]*4000)) treetrace = readTreeTrace(results[z] + concat[z] + "standard_relaxed_trees.trees", treetype="clock") mapTree(treetrace, results[z] + concat[z] + "standard_relaxed_tree.tre") clear(monitors, moves, Q, root_time, rho, timetree, A_B, unknown_time, lambda, ucln_log_mean, branch_rates, seq)