################################################################################ # # RevBayes Example: Bayesian inference of phylogeny using a Jukes-Cantor model # # This file: Runs the full MCMC on a single gene under the Jukes-Cantor # subsitution model using an unconstrained (unrooted) tree model. # # authors: Sebastian Hoehna, Tracy A. Heath, Michael Landis and Brian R. Moore # ################################################################################ ####################### # Reading in the Data # ####################### ###### This just defines a single model for all sites ####### ### Read in sequence data for both genes data <- readDiscreteCharacterData("primates_cytb.nex") # Get some useful variables from the data. We need these later on. n_species <- data.ntaxa() n_sites <- data.nchar() names <- data.names() n_branches <- 2 * n_species - 3 # set my move index mi = 0 ###################### # Substitution Model # ###################### #### specify the Jukes-Cantor substitution model applied uniformly to all sites ### Q := fnJC(4) ############## # Tree model # ############## #### Specify a uniform prior on the tree topology #### topology ~ dnUniformTopology(names) # moves on the tree moves[++mi] = mvNNI(topology) moves[++mi] = mvSPR(topology) #### Specify a prior and moves on the branch lengths #### # create a random variable for each branch length using a for loop for (i in 1:n_branches) { # We use here the exponential distribution with rate 1.0 as the branch length prior br_lens[i] ~ dnExponential(10.0) # Add a move for the branch length. We just take a simple scaling move since the value is a positive real number. moves[++mi] = mvScale(br_lens[i]) } TL := sum(br_lens) # Build the tree by combining the topology with the branch length. phylogeny := treeAssembly(topology, br_lens) ################### # PhyloCTMC Model # ################### # the sequence evolution model seq ~ dnPhyloCTMC(tree=phylogeny, Q=Q, type="DNA") # attach the data seq.clamp(data) ############# # THE Model # ############# # We define our model. # We can use any node of our model as a handle, here we chose to use the rate matrix. mymodel = model(Q) monitors[1] = mnModel(filename="output/primates_cytb_JC_posterior.log",printgen=10, separator = TAB) monitors[2] = mnFile(filename="output/primates_cytb_JC_posterior.trees",printgen=10, separator = TAB, phylogeny) monitors[3] = mnScreen(printgen=1000, TL) mymcmc = mcmc(mymodel, monitors, moves) mymcmc.burnin(generations=10000,tuningInterval=1000) mymcmc.run(generations=30000) # Now, we will analyze the tree output. # Let us start by reading in the tree trace treetrace = readTreeTrace("output/primates_cytb_JC_posterior.trees", treetype="non-clock") # and get the summary of the tree trace treetrace.summarize() mapTree(treetrace,"output/primates_cytb_JC.tree") # you may want to quit RevBayes now #q()