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 metabolomic data > datafilename <- system.file("extdata", "metabolo_sample.txt", package="DRomics") > (o <- continuousomicdata(datafilename, 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 0.69 1.1 1.79 2.92 4.78 7.76 10 6 2 2 2 6 2 Number of items: 109 Identifiers of the first 20 items: [1] "P_2" "P_4" "P_5" "P_6" "P_7" "P_10" "P_11" "P_12" "P_14" "P_16" [11] "P_19" "P_21" "P_22" "P_26" "P_32" "P_34" "P_35" "P_36" "P_37" "P_38" Warning message: In continuousomicdata(datafilename, check = TRUE) : We recommend you to check that your omic data were correctly pretreated before importation. In particular data (e.g. metabolomic signal) should have been log-transformed, without replacing 0 values by NA values (consider using the half minimum method instead for example). > if (visualize) + plot(o) > > # item selection using the quadratic method > # options to put in shiny : select.method (3 methods), FDR (numerical positive value < 1) > (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.001)) Removing intercept from test coefficients Number of selected items using a quadratic trend test with an FDR of 0.001: 15 Identifiers of the responsive items: [1] "P_120" "AP_M52" "P_135" "P_78" "AP_M78" "AP_M44" "AP_M54" "AP_M35" [9] "P_157" "P_177" "AP_M1" "P_39" "AP_M39" "P_62" "AP_M12" > if (visualize) + { + (s_lin <- itemselect(o, select.method = "linear", FDR = 0.001)) + (s_ANOVA <- itemselect(o, select.method = "ANOVA", FDR = 0.001)) + + } > > # no options in shiny > (f <- drcfit(s_quad, progressbar = TRUE)) The fitting may be long if the number of selected items is high. | | | 0% | |===== | 7% | |========= | 13% | |============== | 20% | |=================== | 27% | |======================= | 33% | |============================ | 40% | |================================= | 47% | |===================================== | 53% | |========================================== | 60% | |=============================================== | 67% | |=================================================== | 73% | |======================================================== | 80% | |============================================================= | 87% | |================================================================= | 93% | |======================================================================| 100% Results of the fitting using the AICc to select the best fit model Distribution of the chosen models among the 15 fitted dose-response curves: Hill linear exponential Gauss-probit 0 6 8 1 log-Gauss-probit 0 Distribution of the trends (curve shapes) among the 15 fitted dose-response curves: U dec inc 1 9 5 > f$fitres id irow adjpvalue model nbpar b c d 1 P_120 59 1.022669e-05 linear 2 0.14162429 NA 4.836753 2 AP_M52 100 1.022669e-05 exponential 3 0.44676314 NA 6.849541 3 P_135 69 1.182072e-05 exponential 3 0.04311356 NA 5.807941 4 P_78 44 1.746189e-05 exponential 3 0.04861713 NA 5.758516 5 AP_M78 109 2.452999e-05 linear 2 -0.04907406 NA 6.830002 6 AP_M44 96 3.850161e-05 linear 2 -0.05828339 NA 7.316757 7 AP_M54 101 2.019833e-04 exponential 3 0.40825132 NA 7.512451 8 AP_M35 93 2.168516e-04 exponential 3 0.57880272 NA 6.781761 9 P_157 74 2.756153e-04 linear 2 0.08013911 NA 4.332679 10 P_177 79 2.756153e-04 exponential 3 0.03223247 NA 5.283958 11 AP_M1 81 2.756153e-04 exponential 3 0.56848216 NA 5.880813 12 P_39 21 7.332645e-04 Gauss-probit 4 1.76061106 6.065256 6.065256 13 AP_M39 94 7.506425e-04 linear 2 -0.06388180 NA 6.691864 14 P_62 35 8.377951e-04 exponential 3 -0.01448606 NA 5.270084 15 AP_M12 85 9.166752e-04 linear 2 -0.04107552 NA 7.591279 e f SDres typology trend y0 yatdosemax 1 NA NA 0.26453114 L.inc inc 4.836753 5.935757 2 -0.8989528 NA 0.12907637 E.dec.convex dec 6.849541 6.402858 3 2.6724129 NA 0.10662804 E.inc.convex inc 5.807941 6.551317 4 3.3406581 NA 0.08320269 E.inc.convex inc 5.758516 6.206046 5 NA NA 0.07511639 L.dec dec 6.830002 6.449187 6 NA NA 0.10624269 L.dec dec 7.316757 6.864478 7 -0.4369577 NA 0.15227062 E.dec.convex dec 7.512451 7.104199 8 -0.4796435 NA 0.22895443 E.dec.convex dec 6.781761 6.202958 9 NA NA 0.21005208 L.inc inc 4.332679 4.954559 10 2.2690295 NA 0.17416352 E.inc.convex inc 5.283958 6.237017 11 -2.2721901 NA 0.20358047 E.dec.convex dec 5.880813 5.331017 12 0.5665977 -0.6516214 0.25511702 GP.U U 5.446520 6.065102 13 NA NA 0.15734348 L.dec dec 6.691864 6.196142 14 1.9335368 NA 0.18516350 E.dec.concave dec 5.270084 4.483012 15 NA NA 0.10465889 L.dec dec 7.591279 7.272533 yrange maxychange xextrem yextrem 1 1.0990045 1.0990045 NA NA 2 0.4466835 0.4466835 NA NA 3 0.7433767 0.7433767 NA NA 4 0.4475295 0.4475295 NA NA 5 0.3808147 0.3808147 NA NA 6 0.4522791 0.4522791 NA NA 7 0.4082513 0.4082513 NA NA 8 0.5788027 0.5788027 NA NA 9 0.6218795 0.6218795 NA NA 10 0.9530600 0.9530600 NA NA 11 0.5497964 0.5497964 NA NA 12 0.6514669 0.6185823 0.5665977 5.413635 13 0.4957227 0.4957227 NA NA 14 0.7870722 0.7870722 NA NA 15 0.3187460 0.3187460 NA NA > if (visualize) + plot(f) > > if (visualize) + { + # evaluate the impact of preventsfitsoutofrange, enablesfequal0inGP, enablesfequal0inlGP + data(Scenedesmus_metab) + (o1 <- continuousomicdata(Scenedesmus_metab, check = TRUE)) + + # datafilename <- system.file("extdata", "metabolo_sample.txt", package="DRomics") + # (o1 <- continuousomicdata(datafilename, check = TRUE)) + s_quad1 <- itemselect(o1, select.method = "quadratic", FDR = 0.1) + (f1 <- drcfit(s_quad1, + preventsfitsoutofrange = FALSE, + enablesfequal0inGP = FALSE, + enablesfequal0inLGP = FALSE, + progressbar = TRUE)) + (f1bis <- drcfit(s_quad1, + preventsfitsoutofrange = TRUE, + enablesfequal0inGP = FALSE, + enablesfequal0inLGP = FALSE, + progressbar = TRUE)) + (f1ter <- drcfit(s_quad1, + preventsfitsoutofrange = TRUE, + enablesfequal0inGP = TRUE, + enablesfequal0inLGP = TRUE, + progressbar = TRUE)) + + (idremovedinf1bis <- f1$fitres$id[!is.element(f1$fitres$id, f1bis$fitres$id)]) + # targetplot(items = idremovedinf1bis, f1) + + (idchanged <- f1bis$fitres$id[which(f1bis$fitres$model != f1ter$fitres$model | + f1bis$fitres$f != f1ter$fitres$f)]) + # targetplot(items = idchanged, f1bis, dose_log_transfo = TRUE) + # targetplot(items = idchanged, f1ter, dose_log_transfo = TRUE) + + f1bis$fitres[f1bis$fitres$id %in% idchanged, ] + f1ter$fitres[f1ter$fitres$id %in% idchanged, ] + # no impact on those data + + } > > > > if (visualize) + { + # various plot of fitted curves (without data) + 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)) 4 BMD-xfold values and 0 BMD-zSD values are not defined (coded NaN as the BMR stands outside the range of response values defined by the model). 5 BMD-xfold values and 0 BMD-zSD values could not be calculated (coded NA as the BMR stands within the range of response values defined by the model but outside the range of tested doses). > if (visualize) + (r.2 <- bmdcalc(f, z = 2, x = 50)) > > # plot of BMD > # options in shiny : BMDtype (2 possibilities), plottype (3 possibilities), by (3 possibilities) > # hist.bins (integer for hist only) > if (visualize) + { + plot(r, BMDtype = "zSD", plottype = "ecdf", by = "none") + plot(r, BMDtype = "xfold", plottype = "ecdf", by = "none") + + plot(r, plottype = "hist", by = "none", hist.bins = 10) + plot(r, plottype = "density", by = "none") + + plot(r, plottype = "ecdf", by = "trend", hist.bins = 10) + + } > > if (doboot) + { + niter <- 1000 + niter <- 10 + + # Calculation of confidence intervals on BMDs by Bootstrap + 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.87 0.89 10.70