## ----setup, include = FALSE------------------------------------------------------------------------------------------- options(width=120) knitr::opts_chunk$set( collapse = FALSE, eval = FALSE, comment = " ", tidy.opts =list(width.cutoff=80), tidy = TRUE, size="small" ) ## ---- include=FALSE--------------------------------------------------------------------------------------------------- # library(AlphaBeta) # library(data.table) ## --------------------------------------------------------------------------------------------------------------------- # # Sample "generations.fn" file # generation.fn <- system.file("extdata","generations.fn", package="AlphaBeta") # file <- fread(generation.fn) # head(file) # ## ----include=FALSE---------------------------------------------------------------------------------------------------- # df<-read.csv(generation.fn) # df$filename<-sub("^",paste0(dirname(generation.fn),"/"),df$filename ) # write.csv(df, file = paste0(dirname(generation.fn),"/tmp_generation.fn"),row.names=FALSE,quote=FALSE) # generation.fn<- system.file("extdata","tmp_generation.fn", package="AlphaBeta") # ## ---- eval=FALSE, r,paged.print=FALSE--------------------------------------------------------------------------------- # dMatrix(genTable =generation.fn, # cytosine = "CG", # posteriorMaxFilter = 0.99 ) ## ----paged.print=FALSE------------------------------------------------------------------------------------------------ # #Sample output from dMatrix function # head(fread(system.file("extdata/dm","AB-dMatrix-CG-0.99.csv", package="AlphaBeta"))) ## ----eval=FALSE, paged.print=FALSE------------------------------------------------------------------------------------ # rc.meth.lvl(genTable = generation.fn , # cytosine = "CG", # posteriorMaxFilter = 0.99 ) ## ---- paged.print=FALSE----------------------------------------------------------------------------------------------- # #Sample output from proportions function # head(fread(system.file("extdata/dm","AB-methprop-CG-0.99.csv", package="AlphaBeta"))) # ## ---- paged.print=FALSE---------------------------------------------------------------------------------------------- # #Sample file # head(fread(system.file("extdata/dm","sampleInfo.csv", package="AlphaBeta"))) # ## ---- paged.print=FALSE----------------------------------------------------------------------------------------------- # #Sample file # head(fread(system.file("extdata/dm","branchPoints.csv", package="AlphaBeta"))) ## --------------------------------------------------------------------------------------------------------------------- # props.name <- read.table(system.file("extdata/dm","AB-methprop-CG-0.99.csv", package="AlphaBeta"), sep="\t", header=TRUE) # sample.info <-read.table(system.file("extdata/dm","sampleInfo.csv", package="AlphaBeta"), sep="\t", header=TRUE) # branch.points <-read.table(system.file("extdata/dm","branchPoints.csv", package="AlphaBeta"), sep="\t", header=TRUE) # dmatrix <-read.table(system.file("extdata/dm","AB-dMatrix-CG-0.99.csv", package="AlphaBeta"), sep="\t", header=TRUE) # context <- "CG" ## --------------------------------------------------------------------------------------------------------------------- # pedigree <- convertDMATRIX(sample.info=sample.info, # branch.points=branch.points, # dmatrix=dmatrix, # design="sibling") # head(pedigree) ## ----echo=FALSE------------------------------------------------------------------------------------------------------- # dt <- pedigree[,2] + pedigree[,3] - 2 * pedigree[,1] # plot(dt, pedigree[,"D.value"], ylab="Divergence value", xlab=expression(paste(Delta, " t"))) ## --------------------------------------------------------------------------------------------------------------------- # outliers <- "none" # dmatrix <- dmatrix[which(dmatrix[,1] != outliers), ] # dmatrix <- dmatrix[which(dmatrix[,2] != outliers), ] # pedigree <- pedigree[c(as.numeric(rownames(dmatrix))),] # # props <- props.name[which(as.character(props.name[,2]) == context),] # props <- props.name[which(!is.element(props.name[,1], outliers) == TRUE),] ## --------------------------------------------------------------------------------------------------------------------- # p0uu_in <- 1-mean(as.numeric(as.character(props[,3]))) # p0uu_in ## ----output.lines=5--------------------------------------------------------------------------------------------------- # # output directory # output.data.dir<-paste0(getwd(),"/") # # output <- ABneutral( # pedigree.data = pedigree, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 1, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "CG_global_estimates_ABneutral" # ) ## ----output.lines=5--------------------------------------------------------------------------------------------------- # summary(output) ## ----output.lines=5--------------------------------------------------------------------------------------------------- # head(output$pedigree) ## --------------------------------------------------------------------------------------------------------------------- # output <- ABselectMM( # pedigree.data = pedigree, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 1, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "CG_global_estimates_ABselectMM" # ) # # summary(output) ## --------------------------------------------------------------------------------------------------------------------- # output <- ABselectUU( # pedigree.data = pedigree, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 1, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "CG_global_estimates_ABselectUU" # ) # # summary(output) ## --------------------------------------------------------------------------------------------------------------------- # output <- ABnull(pedigree.data = pedigree, # out.dir = output.data.dir, # out.name = "CG_global_estimates_ABnull") # # summary(output) ## --------------------------------------------------------------------------------------------------------------------- # file1 <- # system.file("extdata/models/", # "CG_global_estimates_ABneutral.Rdata", # package = "AlphaBeta") # file2 <- # system.file("extdata/models/", # "CG_global_estimates_ABnull.Rdata", # package = "AlphaBeta") # # out <- FtestRSS(pedigree.select = file1, # pedigree.null = file2) # # out$Ftest ## --------------------------------------------------------------------------------------------------------------------- # file1 <- # system.file("extdata/models/", # "CG_global_estimates_ABselectMM.Rdata", # package = "AlphaBeta") # file2 <- # system.file("extdata/models/", # "CG_global_estimates_ABnull.Rdata", # package = "AlphaBeta") # # out <- FtestRSS(pedigree.select = file1, # pedigree.null = file2) # # out$Ftest ## --------------------------------------------------------------------------------------------------------------------- # file1 <- # system.file("extdata/models/", # "CG_global_estimates_ABselectUU.Rdata", # package = "AlphaBeta") # file2 <- # system.file("extdata/models/", # "CG_global_estimates_ABnull.Rdata", # package = "AlphaBeta") # # out <- FtestRSS(pedigree.select = file1, # pedigree.null = file2) # # out$Ftest ## --------------------------------------------------------------------------------------------------------------------- # # inputModel <- system.file("extdata/models/", # "CG_global_estimates_ABneutral.Rdata", # package = "AlphaBeta") # # # Bootstrapping models CG # output.data.dir <-paste0(getwd(),"/") # # Boutput <- BOOTmodel( # pedigree.data = inputModel, # Nboot = 2, # out.dir = output.data.dir, # out.name = "Boot_CG_global_estimates_ABneutral" # ) # # summary(Boutput) ## --------------------------------------------------------------------------------------------------------------------- # Boutput$standard.errors ## --------------------------------------------------------------------------------------------------------------------- # props.name <- read.table(system.file("extdata/soma/","AB-methprop-CG-0.99.csv", package="AlphaBeta"), sep="\t", header=TRUE,stringsAsFactors = FALSE) # sample.info <-read.table(system.file("extdata/soma/","sampleInfo.csv", package="AlphaBeta"), sep="\t", header=TRUE,stringsAsFactors = FALSE) # dmatrix <-read.table(system.file("extdata/soma/","AB-dMatrix-CG-0.99.csv", package="AlphaBeta"), sep="\t", header=TRUE,stringsAsFactors = FALSE) ## --------------------------------------------------------------------------------------------------------------------- # head(props.name) # head(sample.info) # head(dmatrix) ## --------------------------------------------------------------------------------------------------------------------- # # 'tall' is total age of tree # pedigree.out <- makePHYLO(tall=330, pedigree = dmatrix, sample.info = sample.info) # pedigree.out <- pedigree.out[[1]] # head(pedigree.out) ## --------------------------------------------------------------------------------------------------------------------- # p0uu_in <- mean(props[,3]) # p0uu_in ## --------------------------------------------------------------------------------------------------------------------- # outneutral <- ABneutralSOMA( # pedigree.data = pedigree.out, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 0.001, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "ABneutralSOMA_CG_estimates" # ) # summary(outneutral) ## --------------------------------------------------------------------------------------------------------------------- # head(outneutral$pedigree) ## --------------------------------------------------------------------------------------------------------------------- # outselectMM <- ABselectMMSOMA( # pedigree.data = pedigree.out, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 0.001, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "ABselectMMSOMA_CG_estimates" # ) # # summary(outselectMM) ## --------------------------------------------------------------------------------------------------------------------- # outselectUU <- ABselectUUSOMA( # pedigree.data = pedigree.out, # p0uu = p0uu_in, # eqp = p0uu_in, # eqp.weight = 0.001, # Nstarts = 2, # out.dir = output.data.dir, # out.name = "ABselectUUSOMA_CG_estimates" # ) # summary(outselectUU) ## --------------------------------------------------------------------------------------------------------------------- # sessionInfo()