Rev Language Reference

dnMultiSpeciesCoalescentUniformPrior - Multispecies coalescent distribution with uniform prior on effective population sizes

Multispecies coalescent distribution describing how gene trees can be generated from within a species tree given effective population sizes. Requires an ultrametric species tree, an upper bound for the uniform prior on effective population sizes (a single real positive), and taxa with species and individual names.

Usage

dnMultiSpeciesCoalescentUniformPrior(TimeTree speciesTree, RealPos max, Taxon[] taxa)

Arguments

speciesTree : TimeTree (pass by const reference)
The species tree in which the gene trees evolve.
max : RealPos (pass by const reference)
The maximum effective population size.
taxa : Taxon[] (pass by value)
The vector of taxa which have species and individual names.

Domain Type

TimeTree

Details

The species tree must be ultrametric. This distribution uses a uniform prior on effective population sizes. As a consequence, effective population sizes are analytically integrated out and treated as nuisance parameters (Hey & Nielsen 2007).

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 maximum effective population size of 50:
popSize <- 50

# Let's simulate gene trees now.
# Taxon names:
for (g in 1:n_genes) {
    for (i in 1:n_species) {
        for (j in 1:n_alleles) {
            taxa[g][(i-1)*n_alleles+j] <- taxon(taxonName="Species_"+i+"_"+j, speciesName="Species_"+i)
        }
    }
    geneTrees[g] ~ dnMultiSpeciesCoalescentUniformPrior(speciesTree=spTree, max=popSize, taxa=taxa[g])
    print(geneTrees[g])
}

# We can save the species tree:
write(spTree, filename=dataFolder+"speciesTree")

# Saving the gene trees:
for (i in 1:(n_genes)) {
    write(geneTrees[i], filename=dataFolder+"geneTree_"+i+".tree")
}

See Also

Edit this page