################################################################################ # # RevBayes Example: Inferring a Phylogeny of Fossil Bears Using the Mk Model # # This file: Runs the full MCMC ... # # authors: April M. Wright, Michael J. Landis, Sebastian Höhna # ################################################################################ ####################### # Reading in the Data # ####################### # Import the morphological character matrix # # this file contains only the taxa for which morphological characters are available # morpho <- readDiscreteCharacterData("data/bears.nex") ## helpers taxa <- morpho.names() num_taxa <- taxa.size() num_branches <- 2 * num_taxa - 2 # Create some vector for the moves and monitors of this analysis moves = VectorMoves() monitors = VectorMonitors() ############## # Tree model # ############## # Set up branch length hyperprior with a move br_len_lambda ~ dnExp(0.2) moves.append( mvScale(br_len_lambda, weight=2) ) # Define the tree parameter. # First, we generate the topology. # We assume a uniform prior on topology. phylogeny ~ dnUniformTopologyBranchLength(taxa, branchLengthDistribution=dnExponential(br_len_lambda)) # compute the tree length from the phylogeny tree_length := phylogeny.treeLength() moves.append( mvNNI(phylogeny, weight=num_branches/2.0) ) moves.append( mvSPR(phylogeny, weight=num_branches/10.0) ) moves.append( mvBranchLengthScale(phylogeny, weight=num_branches) ) ########################################### # Binary morphological substitution model # ########################################### # Create the Q matrix. These data are binary, so we initialize the Jukes-Cantor matrix with # two states Q_morpho <- fnJC(2) # Set up Gamma-distributed rate variation. alpha_morpho ~ dnUniform( 0.0, 1E6 ) rates_morpho := fnDiscretizeGamma( alpha_morpho, alpha_morpho, 4 ) # Moves on the parameters to the Gamma distribution. moves.append( mvScale(alpha_morpho,lambda=1, weight=2.0) ) #Combine all of our elements into a CTMC. Because we have not observed any invariant sites, # we specify the coding is 'variable'. phyMorpho ~ dnPhyloCTMC(tree=phylogeny, siteRates=rates_morpho, Q=Q_morpho, type="Standard") phyMorpho.clamp(morpho) ######## # MCMC # ######## # initialize the model object # mymodel = model(phylogeny) # Create a vector of monitors # # 1. for the full model # monitors.append( mnModel(filename="output/mk.log", printgen=10) ) # 2. the tree # monitors.append( mnFile(filename="output/mk.trees", printgen=10, phylogeny) ) # 3. and a few select parameters to be printed to the screen # monitors.append( mnScreen(printgen=100) ) # Initialize the MCMC object # mymcmc = mcmc(mymodel, monitors, moves, nruns=2, combine="mixed") # Run the MCMC # mymcmc.run(generations=20000, tuningInterval=200) # check the performance of the MCMC/moves mymcmc.operatorSummary() # Read in the tree trace and construct the consensus tree tree # trace = readTreeTrace("output/mk.trees", treetype="non-clock") trace.setBurnin(0.25) # Summarize tree trace and the consensus tree to file mapTree(trace, file="output/mk.map.tre") consensusTree(trace, file="output/mk.majrule.tre") # Quit RevBayes # q()