#################### # Read in the data # #################### data_ITS = readDiscreteCharacterData("data/fagus_ITS.nex") data_matK = readDiscreteCharacterData("data/fagus_matK.nex") data_rbcL = readDiscreteCharacterData("data/fagus_rbcL.nex") # get some useful information about the data taxa = data_ITS.taxa() num_taxa = data_ITS.ntaxa() num_branches = 2 * num_taxa - 3 # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() # name the analysis name = "partition_2" ######################################### # Define the prior on the tree topology # ######################################### # We assume a uniform prior on topology. topology ~ dnUniformTopology(taxa) moves.append( mvNNI(topology, weight=10.0) ) moves.append( mvSPR(topology, weight=10.0) ) ########################################## # Define the prior on the branch lengths # ########################################## for(i in 1:num_branches){ br_lens[i] ~ dnExponential(10.0) moves.append( mvScale(br_lens[i], weight=1.0) ) } TL := sum(br_lens) ################################################ # Combine the tree topology and branch lengths # ################################################ phylogeny := treeAssembly(topology, br_lens) ############################################ # Define the substitution model parameters # ############################################ # ITS has GTR pi_ITS ~ dnDirichlet(v(1,1,1,1)) moves.append( mvBetaSimplex(pi_ITS, weight=1.0) ) er_ITS ~ dnDirichlet(v(1,1,1,1,1,1)) moves.append( mvBetaSimplex(er_ITS, weight=1.0) ) Q_ITS := fnGTR(er_ITS, pi_ITS) # matK and rbcL have GTR pi_matK_rbcL ~ dnDirichlet(v(1,1,1,1)) moves.append( mvBetaSimplex(pi_matK_rbcL, weight=1.0) ) er_matK_rbcL ~ dnDirichlet(v(1,1,1,1,1,1)) moves.append( mvBetaSimplex(er_matK_rbcL, weight=1.0) ) Q_matK_rbcL := fnGTR(er_matK_rbcL, pi_matK_rbcL) ################################################# # Define the model of among-site rate variation # ################################################# # ITS has Gamma alpha_ITS ~ dnUniform(0,1E8) moves.append( mvScale(alpha_ITS, weight=1.0) ) site_rates_ITS := fnDiscretizeGamma(alpha_ITS, alpha_ITS, 4) ######################################### # Define the among-locus rate variation # ######################################### num_sites[1] = data_ITS.nchar() num_sites[2] = data_matK.nchar() + data_rbcL.nchar() relative_rates ~ dnDirichlet(v(1,1)) moves.append( mvBetaSimplex(relative_rates, weight=1.0) ) subset_rates := relative_rates * sum(num_sites) / num_sites ################################# # Define the phyloCTMC model # # (AKA the likelihood function) # ################################# seq_ITS ~ dnPhyloCTMC(tree=phylogeny, branchRates=subset_rates[1], Q=Q_ITS, type="DNA", siteRates=site_rates_ITS) seq_ITS.clamp(data_ITS) # attach the observed data seq_matK ~ dnPhyloCTMC(tree=phylogeny, branchRates=subset_rates[2], Q=Q_matK_rbcL, type="DNA") seq_matK.clamp(data_matK) # attach the observed data seq_rbcL ~ dnPhyloCTMC(tree=phylogeny, branchRates=subset_rates[2], Q=Q_matK_rbcL, type="DNA") seq_rbcL.clamp(data_rbcL) # attach the observed data ######################### # Make the model object # ######################### my_model = model(phylogeny) ##################### # Make the monitors # ##################### monitors.append( mnModel(filename="output/" + name + "/posterior_samples.log",printgen=10, separator = TAB) ) monitors.append( mnFile(filename="output/" + name + "/tree_samples.trees",printgen=10, separator = TAB, phylogeny) ) monitors.append( mnScreen(printgen=100, TL) ) ################ # Run the MCMC # ################ analysis = mcmc(my_model, monitors, moves) analysis.burnin(generations=1000, tuningInterval=100) analysis.operatorSummary() analysis.run(generations=10000) ######################################### # Read in the posterior sample of trees # # and compute the MAP tree # ######################################### # start by reading in the tree trace treetrace = readTreeTrace("output/" + name + "/tree_samples.trees", treetype="non-clock") map_tree = mapTree(treetrace,"output/" + name + "/MAP_tree.tree") ################################### # Run the stepping-stone analysis # ################################### ss_analysis = powerPosterior(my_model, monitors, moves, "output/" + name + "/ss", cats=20, alpha=0.3) ss_analysis.burnin(generations=1000,tuningInterval=100) ss_analysis.run(generations=5000) ss = steppingStoneSampler("output/" + name + "/ss", "power", "likelihood", TAB) ss.marginal() # exit the program q()