## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----installation, echo=TRUE, eval=FALSE-------------------------------------- # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # # BiocManager::install("SOMNiBUS") ## ---- eval=FALSE, echo=FALSE-------------------------------------------------- # ROOT_PACKAGE_PATH <- paste(getwd(), "/", sep = "") # devtools::document(ROOT_PACKAGE_PATH) # devtools::load_all(ROOT_PACKAGE_PATH) ## ----load package, echo=TRUE, message=FALSE, warning=FALSE-------------------- library(SOMNiBUS) ## ----------------------------------------------------------------------------- head(RAdat) ## ----------------------------------------------------------------------------- RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ]) ## ----------------------------------------------------------------------------- out <- binomRegMethModel(data = RAdat.f, n.k = rep(5, 3), p0 = 0.003, p1 = 0.9, Quasi = FALSE, RanEff = FALSE) ## ---- eval=FALSE-------------------------------------------------------------- # out.ctype <- binomRegMethModel(data = RAdat.f, n.k = rep(5, 3), p0 = 0.003, p1 = 0.9, covs = "T_cell") ## ----------------------------------------------------------------------------- as.integer(length(unique(RAdat.f$Position)) / 20) ## ----------------------------------------------------------------------------- out$reg.out ## ---- fig.cap="The estimates (solid red lines) and 95% pointwise confidence intervals (dashed red lines) of the intercept, the smooth effect of cell type and RA on methylation levels. ", fig.height= 4, fig.width=9---- binomRegMethModelPlot(out) ## ---- fig.cap="The estimates (solid red lines) and 95% pointwise confidence intervals (dashed red lines) of the intercept, the smooth effect of cell type and RA on methylation levels. (Same ranges of Y axis.)", fig.height= 4, fig.width=9---- binomRegMethModelPlot(out, same.range = TRUE) ## ----------------------------------------------------------------------------- pos <- out$uni.pos my.p <- length(pos) newdata <- expand.grid(pos, c(0, 1), c(0, 1)) colnames(newdata) <- c("Position", "T_cell", "RA") ## ----------------------------------------------------------------------------- my.pred <- binomRegMethModelPred(out, newdata, type = "link.scale") ## ---- fig.height= 6, fig.width=6, fig.cap="The predicted methylation levels in the logit scale for the 4 groups of samples with different disease and cell type status."---- plot(pos[order(pos)], (my.pred[(newdata$RA == 0 & newdata$T_cell == 0)])[order(pos)], type = "l", xlab = "Position", ylab = "Predicted methylation levels (in logit scale)", col = "blue", main = "Logit scale", ylim = c(min(my.pred), max(my.pred)), lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 0 & newdata$T_cell == 1)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "green", lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 1 & newdata$T_cell == 0)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "red", lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 1 & newdata$T_cell == 1)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "black", lwd = 2 ) legend("top", c("RA MONO", "RA TCELL", "CTRL MONO", "CTRL TCELL"), fill = c("red", "black", "blue", "green"), title = "Disease and Cell Type", bty = "n", cex = 0.8 ) ## ---- fig.height= 6, fig.width=6 , fig.cap="The predicted methylation levels proportion scale (right) for the 4 groups of samples with different disease and cell type status."---- my.pred <- binomRegMethModelPred(out, newdata, type = "proportion") plot(pos[order(pos)], (my.pred[(newdata$RA == 0 & newdata$T_cell == 0)])[order(pos)], type = "l", xlab = "Position", ylab = "Predicted methylation levels (in logit scale)", col = "blue", main = "Proportion scale", ylim = c(min(my.pred), max(my.pred)), lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 0 & newdata$T_cell == 1)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "green", lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 1 & newdata$T_cell == 0)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "red", lwd = 2 ) lines(pos[order(pos)], (my.pred[(newdata$RA == 1 & newdata$T_cell == 1)])[order(pos)], type = "l", xlab = "Position", ylab = "predicted", col = "black", lwd = 2 ) legend("top", c("RA MONO", "RA TCELL", "CTRL MONO", "CTRL TCELL"), fill = c("red", "black", "blue", "green"), title = "Disease and Cell Type", bty = "n", cex = 0.8 ) ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()