data <- readDiscreteCharacterData("MOLECULAR_CONV_SOLAN.nex") n_taxa <- data.ntaxa() taxa <- data.taxa() n_branches <- n_taxa*2 - 2 topology <- readTrees("start_two.tre")[1] mvi = 0 mni = 0 ################### er_prior <- v(1,1,1,1,1,1) er ~ dnDirichlet(er_prior) moves[++mvi] = mvDirichletSimplex(er, weight=3.0) moves[++mvi] = mvBetaSimplex(er, weight=3.0) pi_prior <- v(1,1,1,1) pi ~ dnDirichlet(pi_prior) moves[++mvi] = mvDirichletSimplex(pi, weight=3.0) moves[++mvi] = mvBetaSimplex(pi, weight=3.0) Q := fnGTR(er,pi) alpha_prior_mean <- ln(5.0) alpha_prior_sd <- 0.587405 alpha ~ dnLognormal( alpha_prior_mean, alpha_prior_sd ) gamma_rates := fnDiscretizeGamma( alpha, alpha, 4 ) moves[++mvi] = mvScale(alpha, weight=3.0) pinvar ~ dnBeta(1,1) moves[++mvi] = mvSlide(pinvar, weight=5.0) ################### rho <- 0.000205 root_time ~ dnUniform(46, 86.9) moves[++mvi] = mvScale(root_time,lambda=0.01,weight=3.0) moves[++mvi] = mvScale(root_time,lambda=0.1,weight=3.0) moves[++mvi] = mvScale(root_time,lambda=1.0,weight=3.0) moves[++mvi] = mvSlide(root_time,weight=1.0) speciation ~ dnExponential(10) moves[++mvi] = mvScale(speciation,lambda=1,weight=3.0) moves[++mvi] = mvScale(speciation,lambda=0.1,weight=3.0) moves[++mvi] = mvScale(speciation,lambda=0.01,weight=3.0) relativeExtinction ~ dnBeta(1,1) moves[++mvi] = mvScale(relativeExtinction, lambda = 1, weight=3.0) moves[++mvi] = mvScale(relativeExtinction, lambda = 0.1, weight=3.0) moves[++mvi] = mvScale(relativeExtinction, lambda = 0.01, weight=3.0) extinction := speciation * relativeExtinction timetree ~ dnBDP(lambda=speciation, mu=extinction, rho=rho, rootAge=root_time, samplingStrategy="uniform", condition="nTaxa", taxa=taxa) timetree.setValue(topology) clade_solanoidea = clade("Nicotiana_glauca", "Solanum_jamesii") tmrca_solanoidea := tmrca(timetree, clade_solanoidea) min_age_solanoidea = 23 max_age_solanoidea = 86.9 width_age_prior_solanoidea = (max_age_solanoidea - min_age_solanoidea)/2 mean_age_prior_solanoidea = min_age_solanoidea + width_age_prior_solanoidea solanoidea_calibration ~ dnUniform(tmrca_solanoidea - width_age_prior_solanoidea, tmrca_solanoidea + width_age_prior_solanoidea) solanoidea_calibration.clamp(mean_age_prior_solanoidea) moves[++mvi] = mvSubtreeScale(timetree, weight=5.0) moves[++mvi] = mvNodeTimeSlideUniform(timetree, weight=15.0) moves[++mvi] = mvTreeNodeAgeSlide(timetree, weight=40.0) moves[++mvi] = mvTreeScale(tree=timetree, delta=1.0, weight=3.0) moves[++mvi] = mvTreeScale(tree=timetree, delta=0.1, weight=3.0) moves[++mvi] = mvTreeScale(tree=timetree, delta=0.01, weight=3.0) ################### clock_mean ~ dnUnif(1E-6, 10) moves[++mvi] = mvScale(clock_mean, lambda = 1, weight = 3) moves[++mvi] = mvScale(clock_mean, lambda = 0.1, weight = 3) moves[++mvi] = mvScale(clock_mean, lambda = 0.01, weight = 3) for(i in 1:n_branches){ branch_rates[i] ~ dnLnorm(ln(clock_mean), 0.2937025) moves[++mvi] = mvScale(branch_rates[i], lambda = 1, weight=2.0) } 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, siteRates=gamma_rates, pInv=pinvar, branchRates=branch_rates, type="DNA") seq.clamp(data) ################### mymodel = model(timetree) monitors[++mni] = mnModel(filename="tight_relaxed/output.log", printgen=10, separator = TAB) monitors[++mni] = mnFile(filename="tight_relaxed/output.trees", printgen=10, timetree) monitors[++mni] = mnScreen(printgen=10, root_time, tmrca_solanoidea) mymcmc = mcmc(mymodel, monitors, moves, nruns=2) mymcmc.burnin(generations=0, tuningInterval=250) mymcmc.run(generations=75000)