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")) ############################################################################################################# mean_for_all_loci ~ dnLnorm(-3.00667, 0.587405/4) moves[++mvi] = mvScale(mean_for_all_loci, lambda = 1, weight = 3) moves[++mvi] = mvScale(mean_for_all_loci, lambda = 0.1, weight = 3) moves[++mvi] = mvScale(mean_for_all_loci, lambda = 0.01, weight = 3) loci_partitions_prior <- rep(1, marker_number[z]) loci_partitions_proportions ~ dnDirichlet(loci_partitions_prior) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 0.1, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 0.5, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 1, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 5, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 10, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 50, weight=3.0) moves[++mvi] = mvDirichletSimplex(loci_partitions_proportions, alpha = 100, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 0.1, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 0.5, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 1, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 5, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 10, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 50, weight=3.0) moves[++mvi] = mvBetaSimplex(loci_partitions_proportions, alpha = 100, weight=3.0) for (i in 1:marker_number[z]){ locus_rate[i] := ln(loci_partitions_proportions[i] * (mean_for_all_loci*marker_number[z])) } for (i in 1:marker_number[z]){ branch_rates[i] <- v(0, 0, 0, 0) } for (j in 1:marker_number[z]){ for (i in 1:n_branches){ branch_rates[j][i] ~ dnLnorm(locus_rate[j], 0.587405/4) moves[mvi++] = mvScale(branch_rates[j][i], lambda = 1, weight=2) } } moves[mvi++] = mvVectorScale(branch_rates[j], lambda = 1.0, weight = 2.0) moves[mvi++] = mvVectorSingleElementScale(branch_rates[j], lambda=30.0, weight = 1.0) ############################################################################################################# for (i in 1:marker_number[z]){ seq[i] ~ dnPhyloCTMC(timetree, Q=Q, branchRates=branch_rates[i], type="DNA") seq[i].clamp(data[i]) } ############################################################################################################# mymodel = model(timetree) monitors[++mni] = mnModel(filename = results[z] + concat[z] + "dirichlet_sampling.log", printgen=(5000 + (dirich_scale[z]*4000))/500, separator = TAB) monitors[++mni] = mnFile(filename= results[z] + concat[z] + "dirichlet_trees.trees", printgen=(5000 + (dirich_scale[z]*4000))/500, timetree) monitors[++mni] = mnScreen(printgen=500) mymcmc = mcmc(mymodel, monitors, moves, nruns=1) mymcmc.run(generations=5000 + (dirich_scale[z]*4000)) treetrace = readTreeTrace(results[z] + concat[z] + "dirichlet_trees.trees", treetype="clock") mapTree(treetrace, results[z] + concat[z] + "dirichlet_tree.tre") clear(monitors, moves, Q, root_time, rho, timetree, A_B, unknown_time, branch_rates, seq, mean_for_all_loci, loci_partitions_prior, loci_partitions_proportions, locus_rate)