| Title: | Computational inference of epimutation rates and spectra from high-throughput DNA methylation data in plants |
|---|---|
| Description: | AlphaBeta is a computational method for estimating epimutation rates and spectra from high-throughput DNA methylation data in plants. The method has been specifically designed to: 1. analyze 'germline' epimutations in the context of multi-generational mutation accumulation lines (MA-lines). 2. analyze 'somatic' epimutations in the context of plant development and aging. |
| Authors: | Yadollah Shahryary Dizaji [cre, aut], Frank Johannes [aut], Rashmi Hazarika [aut] |
| Maintainer: | Yadollah Shahryary Dizaji <[email protected]> |
| License: | GPL-3 |
| Version: | 1.21.0 |
| Built: | 2024-10-30 03:30:09 UTC |
| Source: | https://github.com/bioc/AlphaBeta |
This model assumes that heritable gains and losses in cytosine methylation are selectively neutral.
ABneutral(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)ABneutral(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
weight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABneutral RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABneutral" out <- ABneutral(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABneutral" out <- ABneutral(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)
This model assumes that somatically heritable gains and losses in cytosine methylation are selectively neutral.
ABneutralSOMA(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)ABneutralSOMA(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
weight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABneutralSoma RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABneutralSOMA_CG_estimates" out <- ABneutralSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABneutralSOMA_CG_estimates" out <- ABneutralSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)
Run model that considers no accumulation of epimutations (ABnull)
ABnull(pedigree.data, out.dir, out.name)ABnull(pedigree.data, out.dir, out.name)
pedigree.data |
Generation table name, you can find sample file in |
out.dir |
outputdirectory |
out.name |
name of file |
ABnull RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata out.name <- "CG_global_estimates_ABnull" out <- ABnull(pedigree.data = pedigree, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata out.name <- "CG_global_estimates_ABnull" out <- ABnull(pedigree.data = pedigree, out.dir=getwd(), out.name=out.name) summary(out)
Plotting Estimating epimutation
ABplot( pedigree.names, output.dir, out.name, alpha = 0.5, geom.point.size = 2, geom.line.size = 0.9, plot.height = 8, plot.width = 11, plot.type = "both", lsq.line = "theory", intract = FALSE )ABplot( pedigree.names, output.dir, out.name, alpha = 0.5, geom.point.size = 2, geom.line.size = 0.9, plot.height = 8, plot.width = 11, plot.type = "both", lsq.line = "theory", intract = FALSE )
pedigree.names |
Models output AB*.Rdata |
output.dir |
output directory |
out.name |
filename |
alpha |
ggplot parameters |
geom.point.size |
ggplot parameters |
geom.line.size |
ggplot parameters |
plot.height |
ggplot parameters |
plot.width |
ggplot parameters |
plot.type |
type of plot (data.only, fit.only, both) |
lsq.line |
Least Square Regression line (theory or pred) |
intract |
to see intarctive plot. (useing plotly) |
plot
# Get some toy data file <- system.file("extdata/dm/","Col_CG_global_estimates_ABneutral.Rdata", package="AlphaBeta") ABplot(pedigree.names=file, output.dir=getwd(), out.name="ABneutral")# Get some toy data file <- system.file("extdata/dm/","Col_CG_global_estimates_ABneutral.Rdata", package="AlphaBeta") ABplot(pedigree.names=file, output.dir=getwd(), out.name="ABneutral")
This model assumes that heritable losses of cytosine methylation are under negative selection.
ABselectMM(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)ABselectMM(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
nweight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABselectMM RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABselectMM" out <- ABselectMM(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABselectMM" out <- ABselectMM(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)
This model assumes that somatically heritable gains of cytosine methylation are under negative selection.
ABselectMMSOMA( pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name )ABselectMMSOMA( pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name )
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
weight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABneutralSoma RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABselectMMSOMA_CG_estimates" out <- ABselectMMSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABselectMMSOMA_CG_estimates" out <- ABselectMMSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)
This model assumes that heritable gains of cytosine methylation are under negative selection.
ABselectUU(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)ABselectUU(pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name)
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
weight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABselectMM RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABselectUU" out3 <- ABselectUU(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out3)#Get some toy data inFile <- readRDS(system.file("extdata/dm/","output.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 1 Nstarts <- 2 out.name <- "CG_global_estimates_ABselectUU" out3 <- ABselectUU(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out3)
This model assumes that somatically heritable gains of cytosine methylation are under negative selection.
ABselectUUSOMA( pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name )ABselectUUSOMA( pedigree.data, p0uu, eqp, eqp.weight, Nstarts, out.dir, out.name )
pedigree.data |
pedigree data. |
p0uu |
initial proportion of unmethylated cytosines. |
eqp |
equilibrium proportion of unmethylated cytosines. |
eqp.weight |
weight assigned to equilibrium function. |
Nstarts |
iterations for non linear LSQ optimization. |
out.dir |
output directory. |
out.name |
output file name. |
ABneutralSoma RData file.
#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABselectUUSOMA_CG_estimates" out <- ABselectUUSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)#Get some toy data inFile <- readRDS(system.file("extdata/soma/","outputSoma.rds", package="AlphaBeta")) pedigree <- inFile$Pdata p0uu_in <- inFile$tmpp0 eqp.weight <- 0.001 Nstarts <- 2 out.name <- "ABselectUUSOMA_CG_estimates" out <- ABselectUUSOMA(pedigree.data = pedigree, p0uu=p0uu_in, eqp=p0uu_in, eqp.weight=eqp.weight, Nstarts=Nstarts, out.dir=getwd(), out.name=out.name) summary(out)
Bootstrap analysis with the best model
BOOTmodel(pedigree.data, Nboot, out.dir, out.name)BOOTmodel(pedigree.data, Nboot, out.dir, out.name)
pedigree.data |
pedigree data. |
Nboot |
number of boot. |
out.dir |
output directory. |
out.name |
output file name. |
bootstrap result.
## Get some toy data inFile <- system.file("extdata/models/","ABneutral_CG_global_estimates.Rdata", package="AlphaBeta") Nboot <- 4 out.name <-"Boot_CG_global_estimates_ABneutral" Bout <- BOOTmodel(pedigree.data=inFile, Nboot=Nboot, out.dir=getwd(), out.name=out.name) summary(Bout)## Get some toy data inFile <- system.file("extdata/models/","ABneutral_CG_global_estimates.Rdata", package="AlphaBeta") Nboot <- 4 out.name <-"Boot_CG_global_estimates_ABneutral" Bout <- BOOTmodel(pedigree.data=inFile, Nboot=Nboot, out.dir=getwd(), out.name=out.name) summary(Bout)
calculate divergence times of the pedigree
buildPedigree(nodelist, edgelist, cytosine = "CG", posteriorMaxFilter = 0.99)buildPedigree(nodelist, edgelist, cytosine = "CG", posteriorMaxFilter = 0.99)
nodelist |
input file containing information on generation times and pedigree lineages "extdata" called "nodelist.fn" |
edgelist |
input file containing edges |
cytosine |
Type of cytosine (CHH/CHG/CG) |
posteriorMaxFilter |
Filter value, based on posteriorMax |
generating divergence matrices file.
# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") df<-read.csv(file) df$filename <- gsub("^", paste0(dirname(dirname(file)),"/"), df$filename ) write.csv(df, file = paste0(dirname(file),"/", "tmp_nodelist.fn"), row.names=FALSE, quote=FALSE) file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") file2 <- system.file("extdata/dm/","edgelist.fn", package="AlphaBeta") buildPedigree(nodelist = file, edgelist=file2, cytosine="CG", posteriorMaxFilter=0.99)# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") df<-read.csv(file) df$filename <- gsub("^", paste0(dirname(dirname(file)),"/"), df$filename ) write.csv(df, file = paste0(dirname(file),"/", "tmp_nodelist.fn"), row.names=FALSE, quote=FALSE) file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") file2 <- system.file("extdata/dm/","edgelist.fn", package="AlphaBeta") buildPedigree(nodelist = file, edgelist=file2, cytosine="CG", posteriorMaxFilter=0.99)
Estimating epimutation rates from high-throughput DNA methylation data
dMatrix(nodelist, cytosine, posteriorMaxFilter)dMatrix(nodelist, cytosine, posteriorMaxFilter)
nodelist |
list of samples, you can find sample file in "extdata" called "nodelist.fn" |
cytosine |
Type of cytosine (CHH/CHG/CG) |
posteriorMaxFilter |
Filter value, based on posteriorMax ex: >= 0.95 or 0.99 |
generating divergence matrices file.
# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") df<-read.csv(file) df$filename<-sub("^",paste0(dirname(file),"/"),df$filename ) write.csv(df, file = paste0(dirname(file),"tmp_nodelist.fn"),row.names=FALSE,quote=FALSE) file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") dMatrix(file, "CG", 0.99)# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") df<-read.csv(file) df$filename<-sub("^",paste0(dirname(file),"/"),df$filename ) write.csv(df, file = paste0(dirname(file),"tmp_nodelist.fn"),row.names=FALSE,quote=FALSE) file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") dMatrix(file, "CG", 0.99)
Comparison of different models and selection of best model
FtestRSS(pedigree.select, pedigree.null)FtestRSS(pedigree.select, pedigree.null)
pedigree.select |
pedigree model. |
pedigree.null |
ABnull pedigree. |
result of Ftest.
## Get some toy data file1 <- system.file("extdata/models/","ABneutral_CG_global_estimates.Rdata", package="AlphaBeta") file2 <- system.file("extdata/models/","ABnull_CG_global_estimates.Rdata", package="AlphaBeta") out <- FtestRSS(pedigree.select=file1, pedigree.null=file2)## Get some toy data file1 <- system.file("extdata/models/","ABneutral_CG_global_estimates.Rdata", package="AlphaBeta") file2 <- system.file("extdata/models/","ABnull_CG_global_estimates.Rdata", package="AlphaBeta") out <- FtestRSS(pedigree.select=file1, pedigree.null=file2)
Plotting Pedigree tree
plotPedigree( nodelist, edgelist, sampling.design, out.pdf = NULL, output.dir = NULL, plot.width = 11, plot.height = 8, vertex.label = NULL, vertex.size = 12, aspect.ratio = 2.5 )plotPedigree( nodelist, edgelist, sampling.design, out.pdf = NULL, output.dir = NULL, plot.width = 11, plot.height = 8, vertex.label = NULL, vertex.size = 12, aspect.ratio = 2.5 )
nodelist |
input file containing information on generation times and pedigree lineages "extdata" called "nodelist.fn" |
edgelist |
input file containing edges "edgelist.fn" |
sampling.design |
"progenitor.intermediate"; "sibling"; "progenitor.endpoint";"tree" |
out.pdf |
output file name |
output.dir |
output directory |
plot.width |
plotting with |
plot.height |
plotting height |
vertex.label |
label vertix |
vertex.size |
size of vertix |
aspect.ratio |
aspect.ration |
plot pedigree matrices file.
# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") file2 <- system.file("extdata/dm/","edgelist.fn", package="AlphaBeta") plotPedigree(nodelist = file, edgelist=file2, sampling.design="sibling",vertex.label=TRUE, out.pdf="Plot", output.dir=getwd() )# Get some toy data file <- system.file("extdata/dm/","nodelist.fn", package="AlphaBeta") file2 <- system.file("extdata/dm/","edgelist.fn", package="AlphaBeta") plotPedigree(nodelist = file, edgelist=file2, sampling.design="sibling",vertex.label=TRUE, out.pdf="Plot", output.dir=getwd() )
Estimating epimutation rates from high-throughput DNA methylation data
rc.meth.lvl(nodelist, cytosine, posteriorMaxFilter)rc.meth.lvl(nodelist, cytosine, posteriorMaxFilter)
nodelist |
List of samples, you can find sample file in "extdata" called "nodelist.fn" |
cytosine |
Type of cytosine (CHH/CHG/CG) |
posteriorMaxFilter |
Filter value, based on posteriorMax |
rc meth lvl.
## Get some toy data file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") rc.meth.lvl(file, "CG", 0.99)## Get some toy data file <- system.file("extdata/dm/","tmp_nodelist.fn", package="AlphaBeta") rc.meth.lvl(file, "CG", 0.99)