R Under development (unstable) (2024-01-23 r85822 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(DRomics) Loading required package: limma Loading required package: DESeq2 Loading required package: S4Vectors Loading required package: stats4 Loading required package: BiocGenerics Attaching package: 'BiocGenerics' The following object is masked from 'package:limma': plotMA The following objects are masked from 'package:stats': IQR, mad, sd, var, xtabs The following objects are masked from 'package:base': Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted, lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply, union, unique, unsplit, which.max, which.min Attaching package: 'S4Vectors' The following object is masked from 'package:utils': findMatches The following objects are masked from 'package:base': I, expand.grid, unname Loading required package: IRanges Attaching package: 'IRanges' The following object is masked from 'package:grDevices': windows Loading required package: GenomicRanges Loading required package: GenomeInfoDb Loading required package: SummarizedExperiment Loading required package: MatrixGenerics Loading required package: matrixStats Attaching package: 'MatrixGenerics' The following objects are masked from 'package:matrixStats': colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars Loading required package: Biobase Welcome to Bioconductor Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'. Attaching package: 'Biobase' The following object is masked from 'package:MatrixGenerics': rowMedians The following objects are masked from 'package:matrixStats': anyMissing, rowMedians DRomics has been loaded. !!!! IMPORTANT CHANGES IN DEFAULT PLOT ARGUMENTS !!!! Now all the plot functions use by default a log10 scale for dose and BMD values, except curvesplot(), for which the use of a dose log scale requires the specification by the user of a non null minimal value (xmin). We also put the default value of the argument scaling at TRUE in curvesplot() and bmdplotwithgradient(), to focus on shapes of dose-responses rather than on their amplitude. > visualize <- FALSE ## put to TRUE for a manual check of plots > doboot <- FALSE > > # importation and check of apical anchoring data > # datafilename <- system.file("extdata", "apical_anchoring.txt", package="DRomics") > # (o <- continuousanchoringdata(datafilename, backgrounddose = 0.1, check = TRUE)) > data("Scenedesmus_apical") > (o <- continuousanchoringdata(Scenedesmus_apical, backgrounddose = 0.1, check = TRUE)) Elements of the experimental design in order to check the coding of the data: Tested doses and number of replicates for each dose: 0 2.4 3.8 6.2 10.1 16.5 26.8 43.5 70.7 12 6 2 2 2 6 2 2 2 Number of endpoints: 2 Names of the endpoints: [1] "growth" "photosynthesis" Warning message: In continuousanchoringdata(Scenedesmus_apical, backgrounddose = 0.1, : We recommend you to check that your anchoring data are continuous and defined in a scale that enable the use of a normal error model (needed at each step of the workflow including the selection step). > > > # Use of only one endpoint > #(o <- continuousanchoringdata(Scenedesmus_apical[c(1,2), ], > # backgrounddose = 0.1,check = TRUE)) # growth > # (o <- continuousanchoringdata(Scenedesmus_apical[c(1,3), ], > # backgrounddose = 0.1, check = TRUE)) # photosynthesis > if (visualize) + plot(o) > > # item selection using the three methods > (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.05)) Number of selected items using a quadratic trend test with an FDR of 0.05: 2 Identifiers of the responsive items: [1] "growth" "photosynthesis" > if (visualize) + { + (s_lin <- itemselect(o, select.method = "linear", FDR = 0.05)) + (s_ANOVA <- itemselect(o, select.method = "ANOVA", FDR = 0.05)) + } > > # fit > (f <- drcfit(s_quad, progressbar = TRUE)) The fitting may be long if the number of selected items is high. | | | 0% | |=================================== | 50% | |======================================================================| 100% Results of the fitting using the AICc to select the best fit model Distribution of the chosen models among the 2 fitted dose-response curves: Hill linear exponential Gauss-probit 0 0 0 2 log-Gauss-probit 0 Distribution of the trends (curve shapes) among the 2 fitted dose-response curves: U 2 > if (visualize) + { + f$fitres + plot(f) + plot(f, dose_log_transfo = FALSE) + plot(f, plot.type = "dose_residuals") + + } > > # various plot of fitted curves (without data) > if (visualize) + { + curvesplot(f$fitres, xmax = max(f$omicdata$dose), + facetby = "model", colorby = "model") + curvesplot(f$fitres, xmax = max(f$omicdata$dose), + facetby = "typology") + + # plot of selection of curves + curvesplot(f$fitres[f$fitres$trend == "U", ], xmax = max(f$omicdata$dose), + facetby = "id") + + } > > > # calculation of benchmark doses > # options in shiny : z (numerical positive value), x (numerical positive value : percentage) > (r <- bmdcalc(f, z = 1, x = 10)) BMD-xfold and BMD-SD values could be calculated on all the curves (the BMR always stands within the range of response values defined by the model and within the range of tested doses). > if (visualize) + { + r$res + + # plot of BMD with gradient + bmdplotwithgradient(r$res, xmax = max(f$omicdata$dose)) + } > > # Calculation of confidence intervals on BMDs by Bootstrap > if (doboot) + { + niter <- 1000 + niter <- 10 + (b <- bmdboot(r, niter = niter)) # niter should be fixed at least at 1000 to get a reasonable precision + if (visualize) + plot(b) + + } > > proc.time() user system elapsed 9.07 0.71 9.78