JC_Model.Rev

################################################################################
#
# 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 #
#######################

# 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)