In ndukler/epiAllele: Tests for signals of variable rates of evolution between gene classes and lineages for probabalistic alleles

First create some data and a tree. Pretend that sites are split into two classes, "Enhancer" and "Promoter" and also split by ontology classes "A" and "B".

library(epiAllele)
## Set seed for consistency
set.seed(123)
## Simulate random data
species=c("A","B","C","D","E")
aData=lapply(as.list(species),function(x) matrix(runif(120),ncol=3))
names(aData)=species

## Create test trees, both one that will fail, and one that will pass
tree=ape::rtree(n = length(species),tip.label = species)

## Create a data frame with labels
siteLabels=data.frame(cre.class=rep(c("Enhancer","Promoter"),each=nrow(aData[[1]])/2),
                      go.class=rep(c("A","B"),times=nrow(aData[[1]])/2))

Now create the allele data object:

ad=alleleData(data=aData,tree=tree,siteInfo = siteLabels)

Now create a rate model specifying that genomic regions should be seperated by both cre.class and go.class.

rateMod=rateModel(ad,siteLabelCriteria = c("cre.class","go.class"))

Note that the model assumes that all different genomic regions have different rate parameters. This can be seen using the plotTree command.

plotTree(obj = rateMod,colorByRate="index")

Now we can fit the model and plot the tree with values

fittedModel=fit(rateMod)
plotTree(obj = fittedModel$model,colorByRate="value")

Say we want to allow for lineage specific rates. Then we can extract and edge table from the allele data object and add a column called edgeGroup to

egt=getEdgeTable(ad)
egt[edgeID %in% c("8-4","8-3"),edgeGroup:="e1"]
egt[!edgeID %in% c("8-4","8-3"),edgeGroup:="e2"]
rateMod=rateModel(ad,siteLabelCriteria = c("cre.class","go.class"),lineageTable = egt)
plotTree(rateMod,colorByRate="index")

Then say we want to tie rates between two sets of sites for a specific group of edges.

rateModTied=tieRates(obj=rateMod,siteLabelA = "Enhancer_A",edgeGroupA = "e1",siteLabelB = "Enhancer_B",edgeGroupB = "e1")
plotTree(obj = rateModTied,colorByRate="index")

Then say we want to tie rates between two sets edge groups within a siteLabel

rateModelTied2=tieRates(rateMod,siteLabelA = "Enhancer_B",edgeGroupA = "e1",siteLabelB = "Enhancer_B",edgeGroupB = "e2")
plotTree(obj = rateModelTied2,colorByRate="index")

ndukler/epiAllele documentation built on May 5, 2019, 4:50 p.m.