Rev Language Reference

mvSpeciesNarrow - Narrow-exchange joint move on species tree and gene trees for multispecies coalescent models.

Makes a narrow-exchange move both in the species tree and in the gene trees that contain nodes of the relevant populations.

Usage

mvSpeciesNarrow(TimeTree speciesTree, RealPos weight, Probability tuneTarget)

Arguments

speciesTree : TimeTree (<stochastic> pass by reference)
The species tree variable on which this move operates.
weight : RealPos (pass by value)
The weight how often on average this move will be used per iteration.
Default : 1
tuneTarget : Probability (pass by value)
The acceptance probability targeted by auto-tuning.
Default : 0.44

Details

The species tree must be ultrametric. All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. This move jointly performs narrow exchange moves (Nearest-Neighbor Interchanges without branch length alterations) on the species tree and on gene trees, all of which must be ultrametric.

Example

# We are going to save the trees we simulate in the folder simulatedTrees:
dataFolder = "simulatedTrees/"
# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:
n_species <- 10
n_genes <- 2
n_alleles <- 3
# we simulate an ultrametric species tree:
# Species names:
for (i in 1:n_species) {
        species[i] <- taxon(taxonName="Species_"+i, speciesName="Species_"+i)
}
spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=10, rho=1, samplingStrategy="uniform", condition="nTaxa", taxa=species)
print(spTree)
# let's pick a constant effective population size of 50:
popSize <- 50
# let's simulate gene trees now:
# taxa names:
for (g in 1:n_genes) {
  for (i in 1:n_species) {
    for (j in 1:n_alleles) {
        taxons[g][(i-1)*n_alleles+j] <- taxon(taxonName="Species_"+i+"_"+j, speciesName="Species_"+i)
    }
  }
  geneTrees[g] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxons[g])
  print(geneTrees[g])
}
# We can save the species tree and the gene trees:
write(spTree, filename=dataFolder+"speciesTree")
# Saving the gene trees
for (i in 1:(n_genes)) {
  write(geneTrees[i], filename=dataFolder+"geneTree_"+i+".tree")
}
# set my move index
mi = 0
move_species_narrow_exchange = mvSpeciesNarrow( speciesTree=spTree, weight=5 )
for (i in 1:n_genes) {
   move_species_narrow_exchange.addGeneTreeVariable( geneTrees[i] )
}
moves[++mi] = move_species_narrow_exchange
# We get a handle on our model.
# We can use any node of our model as a handle, here we choose to use the topology.
mymodel = model(spTree)
# Monitors to check the progression of the program
monitors[1] = mnScreen(printgen=10, spTree)
# Here we use a plain MCMC. You could also set nruns=2 for a replicated analysis
# or use mcmcmc with heated chains.
mymcmc = mcmc(mymodel, monitors, moves, nruns=4)
mymcmc.run(generations=1000)
mymcmc.operatorSummary()

Methods

  • addGeneTreeVariable(TimeTree geneTree)

See Also

Edit this page