taxa <- readTaxonData("osmund_taxa.tsv") data_molecular <- readDiscreteCharacterData("MOLECULAR_OSMUND.nex") data_morphology <- readDiscreteCharacterData("MORPHOLOGY.nex") start <- readTrees("start_3_osmund.tre")[1] n_taxa_extant <- data_molecular.ntaxa() n_taxa <- taxa.size() n_branches <- n_taxa*2 - 2 data_molecular.addMissingTaxa(taxa) data_morphology.addMissingTaxa(taxa) mvi = 0 mni = 0 ########## er_prior <- v(1,1,1,1,1,1) er ~ dnDirichlet(er_prior) moves[++mvi] = mvDirichletSimplex(er, alpha = 0.01, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(er, alpha = 0.01, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(er, alpha = 0.1, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(er, alpha = 0.1, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(er, alpha = 1, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(er, alpha = 1, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(er, alpha = 10, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(er, alpha = 10, weight=3.0, tune=true) pi_prior <- v(1,1,1,1) pi ~ dnDirichlet(pi_prior) moves[++mvi] = mvDirichletSimplex(pi, alpha = 0.01, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(pi, alpha = 0.01, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(pi, alpha = 0.1, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(pi, alpha = 0.1, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(pi, alpha = 1, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(pi, alpha = 1, weight=3.0, tune=true) moves[++mvi] = mvDirichletSimplex(pi, alpha = 10, weight=3.0, tune=true) moves[++mvi] = mvBetaSimplex(pi, alpha = 10, weight=3.0, tune=true) 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, tune=true) pinvar ~ dnBeta(1,1) moves[++mvi] = mvSlide(pinvar, weight=5.0, tune=true) ############ Q_morpho := fnJC(3) ############ rho <- 0.62 origin_time ~ dnUnif(227, 472) moves[++mvi] = mvSlide(origin_time, delta=0.01, weight=5.0, tune=true) moves[++mvi] = mvSlide(origin_time, delta=0.1, weight=5.0, tune=true) moves[++mvi] = mvSlide(origin_time, delta=1, weight=5.0, tune=true) speciation ~ dnExponential(10) moves[++mvi] = mvScale(speciation,lambda=1,weight=3.0, tune=true) moves[++mvi] = mvScale(speciation,lambda=0.1,weight=3.0, tune=true) moves[++mvi] = mvScale(speciation,lambda=0.01,weight=3.0, tune=true) relativeExtinction ~ dnBeta(1,1) moves[++mvi] = mvScale(relativeExtinction, lambda = 1, weight=3.0, tune=true) moves[++mvi] = mvScale(relativeExtinction, lambda = 0.1, weight=3.0, tune=true) moves[++mvi] = mvScale(relativeExtinction, lambda = 0.01, weight=3.0, tune=true) extinction := speciation * relativeExtinction psi ~ dnExponential(50) moves[++mvi] = mvScale(psi, lambda=0.01, weight=1, tune=true) moves[++mvi] = mvScale(psi, lambda=0.1, weight=1, tune=true) moves[++mvi] = mvScale(psi, lambda=1, weight=1, tune=true) fbd_dist = dnFBDP(origin=origin_time, lambda=speciation, mu=extinction, psi=psi, rho=rho, taxa=taxa) ############CLADES all_with_fossils = clade("Todea_barbara", "Todea_papuana", "Leptopteris_fraseri", "Leptopteris_hymenophylloides", "Leptopteris_wilkesiana", "Osmunda_cinnamomea", "Osmunda_claytoniana", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica", "Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis", "Todea_tidwellii", "Osmunda_cinnamomea_Neogene_USA", "Osmunda_cinnamomea_Neogene_Japan", "Osmunda_cinnamomea_Cretaceous_Canada", "Osmunda_precinnamomea", "Osmunda_chengii", "Osmunda_bromeliaefolia", "Osmunda_arnoldii", "Osmunda_dowkeri", "Osmunda_dowkeri_no_morpho", "Osmunda_ilianensis", "Osmunda_nathorstii", "Osmunda_oregonensis", "Osmunda_pluma", "Osmunda_shimokawaensis", "Osmunda_wehrii", "Osmunda_pulchella", "Osmunda_liaoningensis", "Osmunda_wangii", "Osmunda_plumites", "Osmunda_claytoniites", "Osmunda_delawarensis", "Osmunda_diamensis", "Osmunda_japonica_fossilis", "Osmunda_lignita_one", "Osmunda_lignita_two", "Osmunda_macrophylla_one", "Osmunda_macrophylla_two", "Osmunda_parschlugiana", "Osmunda_strozzii", "Osmunda_palaeobanksiifolia", "Osmunda_macrophylla", "Osmunda_vancouverensis", "Osmundopsis_sturii", "Todea_amissa", "Osmundopsis_plectrophora") all_without_fossils = clade("Todea_barbara", "Todea_papuana", "Leptopteris_fraseri", "Leptopteris_hymenophylloides", "Leptopteris_wilkesiana", "Osmunda_cinnamomea", "Osmunda_claytoniana", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica", "Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis") todea_with_fossils = clade("Todea_barbara", "Todea_papuana", "Todea_amissa", "Leptopteris_fraseri", "Leptopteris_hymenophylloides", "Leptopteris_wilkesiana") todea_without_fossils = clade("Todea_barbara", "Todea_papuana", "Leptopteris_fraseri", "Leptopteris_hymenophylloides", "Leptopteris_wilkesiana") osmunda_overall_with_fossils = clade("Osmunda_cinnamomea", "Osmunda_claytoniana", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica", "Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis", "Osmunda_cinnamomea_Neogene_USA", "Osmunda_cinnamomea_Neogene_Japan", "Osmunda_cinnamomea_Cretaceous_Canada", "Osmunda_precinnamomea", "Osmunda_bromeliaefolia", "Osmunda_arnoldii", "Osmunda_dowkeri", "Osmunda_dowkeri_no_morpho", "Osmunda_ilianensis", "Osmunda_nathorstii", "Osmunda_oregonensis", "Osmunda_pluma", "Osmunda_shimokawaensis", "Osmunda_wehrii", "Osmunda_liaoningensis", "Osmunda_wangii", "Osmunda_delawarensis", "Osmunda_japonica_fossilis", "Osmunda_lignita_one", "Osmunda_lignita_two", "Osmunda_macrophylla_one", "Osmunda_macrophylla_two", "Osmunda_parschlugiana", "Osmunda_strozzii", "Osmunda_palaeobanksiifolia", "Osmunda_macrophylla", "Osmunda_vancouverensis") osmunda_overall_without_fossils = clade("Osmunda_cinnamomea", "Osmunda_claytoniana", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica", "Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis") osmunda_JAP_VAC_with_fossils = clade("Osmunda_bromeliaefolia", "Osmunda_delawarensis", "Osmunda_japonica_fossilis", "Osmunda_lignita_one", "Osmunda_lignita_two", "Osmunda_macrophylla_one", "Osmunda_macrophylla_two", "Osmunda_parschlugiana", "Osmunda_strozzii", "Osmunda_palaeobanksiifolia", "Osmunda_macrophylla", "Osmunda_ilianensis", "Osmunda_arnoldii", "Osmunda_dowkeri", "Osmunda_dowkeri_no_morpho", "Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica", "Osmunda_shimokawaensis") osmunda_JAP_VAC_without_fossils = clade("Osmunda_japonica", "Osmunda_lancea", "Osmunda_regalis", "Osmunda_banksiaefolia", "Osmunda_vachellii", "Osmunda_javanica") osmunda_JAP_REG_with_fossils = clade("Osmunda_shimokawaensis", "Osmunda_regalis", "Osmunda_japonica", "Osmunda_lancea") osmunda_JAP_REG_without_fossils = clade("Osmunda_regalis", "Osmunda_japonica", "Osmunda_lancea") ############# constraints = v(todea_with_fossils, osmunda_overall_with_fossils, osmunda_JAP_VAC_with_fossils, osmunda_JAP_REG_with_fossils) fbd_tree ~ dnConstrainedTopology(fbd_dist, constraints=constraints) moves[++mvi] = mvFNPR(fbd_tree, weight=150) moves[++mvi] = mvCollapseExpandFossilBranch(fbd_tree, origin_time, weight=6.0) moves[++mvi] = mvNodeTimeSlideUniform(fbd_tree, weight=40.0) moves[++mvi] = mvRootTimeSlideUniform(fbd_tree, origin_time, weight=5.0) all_with_fossils_age := tmrca(fbd_tree, all_with_fossils) all_without_fossils_age := tmrca(fbd_tree, all_without_fossils) all_difference := all_with_fossils_age - all_without_fossils_age all_calibration ~ dnUniform(all_difference-0.0001, all_difference+0.0001) all_calibration.clamp(0) todea_with_fossils_age := tmrca(fbd_tree, todea_with_fossils) todea_without_fossils_age := tmrca(fbd_tree, todea_without_fossils) todea_difference := todea_with_fossils_age - todea_without_fossils_age todea_calibration ~ dnUniform(todea_difference-0.0001, todea_difference+0.0001) todea_calibration.clamp(0) osmunda_overall_with_fossils_age := tmrca(fbd_tree, osmunda_overall_with_fossils) osmunda_overall_without_fossils_age := tmrca(fbd_tree, osmunda_overall_without_fossils) osmunda_overall_difference := osmunda_overall_with_fossils_age - osmunda_overall_without_fossils_age osmunda_overall_calibration ~ dnUniform(osmunda_overall_difference-0.0001, osmunda_overall_difference+0.0001) osmunda_overall_calibration.clamp(0) osmunda_JAP_VAC_with_fossils_age := tmrca(fbd_tree, osmunda_JAP_VAC_with_fossils) osmunda_JAP_VAC_without_fossils_age := tmrca(fbd_tree, osmunda_JAP_VAC_without_fossils) osmunda_JAP_VAC_difference := osmunda_JAP_VAC_with_fossils_age - osmunda_JAP_VAC_without_fossils_age osmunda_JAP_VAC_calibration ~ dnUniform(osmunda_JAP_VAC_difference-0.0001, osmunda_JAP_VAC_difference+0.0001) osmunda_JAP_VAC_calibration.clamp(0) osmunda_JAP_REG_with_fossils_age := tmrca(fbd_tree, osmunda_JAP_REG_with_fossils) osmunda_JAP_REG_without_fossils_age := tmrca(fbd_tree, osmunda_JAP_REG_without_fossils) osmunda_JAP_REG_difference := osmunda_JAP_REG_with_fossils_age - osmunda_JAP_REG_without_fossils_age osmunda_JAP_REG_calibration ~ dnUniform(osmunda_JAP_REG_difference-0.0001, osmunda_JAP_REG_difference+0.0001) osmunda_JAP_REG_calibration.clamp(0) taxa <- readTaxonData("osmund_taxa.tsv") ############# fossils = fbd_tree.getFossils() for(i in 1:fossils.size()) { t[i] := tmrca(fbd_tree, clade(fossils[i])) a_i = fossils[i].getMinAge() b_i = fossils[i].getMaxAge() F[i] ~ dnUniform(t[i] - b_i, t[i] - a_i) F[i].clamp( 0 ) } moves[++mvi] = mvFossilTimeSlideUniform(fbd_tree, origin_time, weight=5.0) moves[++mvi] = mvTipTimeSlideUniform(fbd_tree, origin_time, weight=5.0) fbd_tree.setValue(start) ################### clock_mean ~ dnUnif(10e-6, 10) moves[++mvi] = mvScale(clock_mean, lambda = 1, weight = 3, tune=true) moves[++mvi] = mvScale(clock_mean, lambda = 0.1, weight = 3, tune=true) moves[++mvi] = mvScale(clock_mean, lambda = 0.01, weight = 3, tune=true) ########### morpho_clock_mean ~ dnUnif(10e-6, 10) moves[++mvi] = mvScale(morpho_clock_mean, lambda = 1, weight = 3, tune=true) moves[++mvi] = mvScale(morpho_clock_mean, lambda = 0.1, weight = 3, tune=true) moves[++mvi] = mvScale(morpho_clock_mean, lambda = 0.01, weight = 3, tune=true) ########### phyMol ~ dnPhyloCTMC(tree=fbd_tree, Q=Q, siteRates=gamma_rates, pInv=pinvar, branchRates=clock_mean, type="DNA") phyMol.clamp(data_molecular) phyMorph ~ dnPhyloCTMC(tree=fbd_tree, Q=Q_morpho, branchRates=morpho_clock_mean, type="Standard") phyMorph.clamp(data_morphology) ########### mymodel = model(Q) monitors[++mni] = mnModel(filename="fbd_strict_high_rate/output.log", printgen=10, separator = TAB) monitors[++mni] = mnFile(filename="fbd_strict_high_rate/output.trees", printgen=10, fbd_tree) monitors[++mni] = mnScreen(printgen=10, origin_time) mymcmc = mcmc(mymodel, monitors, moves, nruns=2) mymcmc.burnin(generations=10000, tuningInterval=100) mymcmc.run(generations=5000000)