GTR_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 GTR+G substitution model applied uniformly to all sites ###
er_prior <- v(1,1,1,1,1,1)
er ~ dnDirichlet(er_prior)
moves[++mi] = mvSimplexElementScale(er,weight=3)


pi_prior <- v(1,1,1,1) 
pi ~ dnDirichlet(pi_prior)
moves[++mi] = mvSimplexElementScale(pi,weight=2)


#### create a deterministic variable for the rate matrix ####
Q := fnGTR(er,pi) 


#############################
# Among Site Rate Variation #
#############################

alpha_prior <- 0.05
alpha ~ dnExponential( alpha_prior )
gamma_rates := fnDiscretizeGamma( alpha, alpha, 4, false )

# add moves for the stationary frequencies, exchangeability rates and the shape parameter
moves[++mi] = mvScale(alpha,weight=2)



# the probability of a site being invariable
pinvar ~ dnBeta(1,1)
# add moves for the probability if a site being invariable
moves[++mi] = mvScale(pinvar)
moves[++mi] = mvSlide(pinvar)


##############
# 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, siteRates=gamma_rates, pInv=pinvar, 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)