## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", crop = NULL ## Related to https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016656.html ) ## ----vignetteSetup, echo=FALSE, warning = FALSE------------------------------- ## Track time spent on making the vignette startTime <- Sys.time() ## Bib setup suppressPackageStartupMessages(library("RefManageR")) ## Write bibliography information bib <- c( R = citation(), BiocStyle = citation("BiocStyle")[1], knitr = citation("knitr")[1], RefManageR = citation("RefManageR")[1], rmarkdown = citation("rmarkdown")[1], sessioninfo = citation("sessioninfo")[1], testthat = citation("testthat")[1], microSTASIS = citation("microSTASIS")[1] ) ## ----"install", eval = FALSE-------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # BiocManager::install("microSTASIS") # ## Check that you have a valid Bioconductor installation # BiocManager::valid() ## ----"citation"--------------------------------------------------------------- ## Citation info citation("microSTASIS") ## ----"start"------------------------------------------------------------------ ## Load the package library("microSTASIS") ## ----"input"------------------------------------------------------------------ ## Show the input data data(clr) clr[1:8, 1:6] ## ----"pairedTimes"------------------------------------------------------------ ## Subseting the initial data to a list of multiple paired times times <- pairedTimes(data = clr, sequential = TRUE, common = "_0_") ## ----"TSE"-------------------------------------------------------------------- ## Loading the TSE package suppressPackageStartupMessages(library(TreeSummarizedExperiment)) ## ----"makeObject"------------------------------------------------------------- ## Creating a TreeSummarizedExperiment object ID_common_time <- rownames(clr) samples <- 1:dim(clr)[1] metadata <- data.frame(samples, stringr::str_split(ID_common_time, "_", simplify = TRUE)) colnames(metadata)[2:4] <- c("ID", "common", "time_point") both_tse <- TreeSummarizedExperiment(assays = list(counts = t(clr)), colData = metadata) altExp(both_tse) <- SummarizedExperiment(assays = SimpleList(data = t(clr)), colData = metadata) assayNames(altExp(both_tse)) <- "data" ## ----"pairedTimesTSE"--------------------------------------------------------- ## Subseting the data as two different TSE objects to two list of multiple paired times TSE_times <- pairedTimes(data = both_tse, sequential = TRUE, assay = "counts", alternativeExp = NULL, ID = "ID", timePoints = "time_point") TSE_altExp_times <- pairedTimes(data = both_tse, sequential = TRUE, assay = "data", alternativeExp = "unnamed1", ID = "ID", timePoints = "time_point") ## ----"iterativeClustering"---------------------------------------------------- ## Main algorithm applied to the matrix-derived object mS <- iterativeClustering(pairedTimes = times, common = "_") ## ----"iterativeClusteringTSE"------------------------------------------------- ## Main algorithm applied to the TSE-derived object mS_TSE <- iterativeClustering(pairedTimes = TSE_times, common = "_") ## ----"visualization"---------------------------------------------------------- ## Prepare the result for the later visualization functions results <- mSpreviz(results = mS, times = times) ## Scatter plot of the stability scores per patient and time plotmSscatter(results = results, order = "median", times = c("t1_t25", "t25_t26"), gridLines = TRUE, sideScale = 0.2) ## Heatmap of the stability scores per patient and time plotmSheatmap(results = results, order = "mean", times = c("t1_t25", "t25_t26"), label = TRUE) ## ----"CV"--------------------------------------------------------------------- ## Cross validation in a parallelized way with different k values cv_klist_k2 <- BiocParallel::bpmapply(iterativeClusteringCV, name = names(times), k = rep(2L, 3), MoreArgs = list(pairedTimes = times, results = mS, common = "_0_")) ## ----"error"------------------------------------------------------------------ ## Calculate the mean absolute error after computing the cross validation MAE_t1_t25 <- mSerrorCV(pairedTime = times$t1_t25, CVklist = cv_klist_k2[[1]], k = 2L) ## ----"error viz"-------------------------------------------------------------- ## Prepare the result for the later visualization functions MAE <- mSpreviz(results = list(MAE_t1_t25), times = list(t1_t25 = times$t1_t25)) ## Heatmap of the mean absolute error values per patient and time plotmSheatmap(results = MAE, times = c("t1_t25", "t25_t26"), label = TRUE, high = 'red2', low = 'forestgreen', midpoint = 5) ## ----"lines"------------------------------------------------------------------ ## Line plot of the stability scores per patient at a given paired times ## Both the proper score after iterative clustering and the cross validation ones plotmSlinesCV(pairedTime = times$t1_t25, CVklist = cv_klist_k2[[1]], k = 2L) ## ----"metadata"--------------------------------------------------------------- ## Create a metadata data frame metadata <- data.frame(Sample = rownames(clr), age = c(rep("youth", 65), rep("old", 131-65))) ## Modify the metadata to the proper format that match with results group <- mSmetadataGroups(metadata = metadata, samples = metadata$Sample, common = "_0_", individuals = results$individual, variable = "age") ## Boxplot with individual stability scores split by groups plotmSdynamics(results, groups = group, points = TRUE, linetype = 0) ## ----createVignette, eval=FALSE----------------------------------------------- # ## Create the vignette # library("rmarkdown") # system.time(render("microSTASIS.Rmd", "BiocStyle::html_document")) # ## Extract the R code # library("knitr") # knit("microSTASIS.Rmd", tangle = TRUE) ## ----reproduce1, echo=FALSE--------------------------------------------------- ## Date the vignette was generated Sys.time() ## ----reproduce2, echo=FALSE--------------------------------------------------- ## Processing time in seconds totalTime <- diff(c(startTime, Sys.time())) round(totalTime, digits = 3) ## ----reproduce3, echo=FALSE------------------------------------------------------------------------------------------- ## Session info library("sessioninfo") options(width = 120) session_info() ## ----vignetteBiblio, results = "asis", echo = FALSE, warning = FALSE, message = FALSE--------------------------------- ## Print bibliography PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))