Learn R Programming

HyPhy (version 1.0)

rgenetree: Samples gene family trees that evolved on a species tree with various conditions.

Description

Samples gene family trees that evolved on the branches of a species tree under the birth-death process. Trees may be conditioned on various priors for the number of reconstructed gene lineages at the root of the species tree, on the total number of genes in the gene tree, on the number of genes in each tip of the species tree, on at least one gene in each tip of the species tree or some combination of those conditions.

Usage

rgenetree(n, spec.phy, lams, mus, root = NULL, genetips = NULL, alltips = FALSE)

Arguments

n
An integer. The number of trees to sample.
spec.phy
The species tree on which to evolve the gene trees of class “phylo” with branch lengths.
lams
The gene duplication rate in units of duplications per spec.phy branch lengths. Either a single numeric for the whole tree or a vector of length dim(spec.phy$edge)[1] for each branch of the species tree.
mus
The gene loss rate in units of losses per spec.phy branch lengths. Either a single numeric for the whole tree or a vector of length dim(spec.phy$edge)[1] for each branch of the species tree.
root
The prior distribution of reconstructed gene lieages at the root of the species tree. Either NULL, a single positive integer, or a positive numeric vector summing to 1. See details.
genetips
The number of genes at the tips of the species tree. Either NULL, a single positive integer or a vector of positive integers of length length(spec.phy$tip.label). See details.
alltips
A logical. Should the sample only include gene trees with at least one gene in each species. Ignored if length(genetips==length(spec.phy$tip.label)

Value

If n==1, a tree of class “phylo”. Otherwise, multiple trees in a list of class “multiPhylo”.

Details

root sets the prior distribution of the number of reconstructed gene lineages at the root node of the species tree. If is.null(root), then the prior on the number of lineages is flat, but if is.null(genetips) then root must have some value. If root is a single positive integer, then all gene trees will start with root reconstructed gene lineages at the root. To establish a more complex prior root should be a numerical vector for which the ith element represents the prior probability that there are i gene lineages at the root, so that sum(root)==1. In that case root can be of any length and the prior probability of any number of lineages greater than length(root) will be 0. It should be noted that this is the prior probability assuming is.null(genetips) & !alltips, changes to those values will affect the probabilities (Hallinan 2012?).

genetips sets the conditions on the number of genes in the tips of the spec.phy. If is.null(genetips) then there are no conditions. If genetips is a single integer, then that value will be the total number of genes in each sampled gene tree, although the number of genes in each species may vary from sample to sample. One can set the number of genes in each tip of the spec.phy with a genetips of length length(spec.phy$tip.label), where the ith member of genetips will be the number of genes in spec.phy$tip.label[i].

References

N. Hallinan. Null models for gene family trees, Math. Biosci. (In review).

See Also

recon.score,plot.phylo,plot.recon(coming soon)

Examples

Run this code
##First we need a simple species tree
spec<-read.tree(text="((A:0.5,B:0.5):0.5,C:1);")
##Now we sample ten gene trees starting with 3 reconstructed gene lineages
phy.all<-rgenetree(10,spec,0.5,0.5,3)
plot(phy.all)

##Now let's make sure that every tip has at least one gene and set an exponential prior on the root
phy.full<-rgenetree(10,spec,0.5,0.5,exp(-(1:20))/sum(exp(-(1:20))),NULL,TRUE)
plot(phy.full)

##Now lets force the whole gene tree to end in 10 genes and set a flat prior for the root
phy.10<-rgenetree(10,spec,0.5,0.5,NULL,10)
plot(phy.10)

##Now lets start with 3 genes, set the number of genes at each tip of spec and vary mu between the branches of spec 
phy.253<-rgenetree(10,spec,0.5,c(0,1,0.2,0.5),3,c(2,5,3))
plot(phy.253)

Run the code above in your browser using DataLab