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 data and normalization if needed > # options to put in shiny : norm.method (4 methods) > ## sample of the transcripto data set > datafilename <- system.file("extdata", "transcripto_very_small_sample.txt", package="DRomics") > (o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess")) Just wait, the normalization using cyclicloess may take a few minutes. 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.223 2.148 3.774 6.631 5 5 5 5 5 5 Number of items: 100 Identifiers of the first 20 items: [1] "1" "2" "3" "4" "5.1" "6.1" "7.1" "8.1" "9.1" "10.1" [11] "11.1" "12.1" "13.1" "14.1" "15" "16.1" "17.1" "18.1" "19.1" "20.1" Data were normalized between arrays using the following method: cyclicloess > plot(o) > if (visualize) + { + (o.2 <- microarraydata(datafilename, check = TRUE, norm.method = "none")) + (o.3 <- microarraydata(datafilename, check = TRUE, norm.method = "quantile")) + (o.4 <- microarraydata(datafilename, check = TRUE, norm.method = "scale")) + + plot(o.2) + plot(o.3) + plot(o.4) + } > > > # 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: 10 Identifiers of the responsive items: [1] "15" "12.1" "27.1" "25.1" "4" "70" "7.1" "88.1" "92" "81" > > if (visualize) + { + (s_lin <- itemselect(o, select.method = "linear", FDR = 0.001)) + (s_ANOVA <- itemselect(o, select.method = "ANOVA", FDR = 0.001)) + } > > # fit of dose response models and choice of the best fit for each item > # no options in shiny > (f <- drcfit(s_quad, progressbar = TRUE)) The fitting may be long if the number of selected items is high. | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100% Results of the fitting using the AICc to select the best fit model Distribution of the chosen models among the 10 fitted dose-response curves: Hill linear exponential Gauss-probit 0 1 4 5 log-Gauss-probit 0 Distribution of the trends (curve shapes) among the 10 fitted dose-response curves: U bell dec inc 2 3 4 1 > f$fitres id irow adjpvalue model nbpar b c d 1 15 15 1.546048e-05 exponential 3 0.071422338 NA 7.740153 2 12.1 12 2.869315e-05 Gauss-probit 5 0.414151082 8.903415 7.564374 3 27.1 27 3.087292e-05 linear 2 -0.108446801 NA 15.608419 4 25.1 25 1.597308e-04 exponential 3 -0.128807373 NA 15.142111 5 4 4 2.302448e-04 Gauss-probit 4 3.079187979 9.851210 9.851210 6 70 70 2.323292e-04 exponential 3 -0.007515144 NA 6.682254 7 7.1 7 2.712029e-04 Gauss-probit 4 2.384260966 9.122630 9.122630 8 88.1 88 4.566344e-04 Gauss-probit 4 2.103260393 11.157946 11.157946 9 92 92 4.566344e-04 Gauss-probit 4 8.383709408 -27.819656 -27.819656 10 81 81 6.448977e-04 exponential 3 -0.025717121 NA 6.713592 e f SDres typology trend y0 yatdosemax 1 2.276376 NA 0.3183292 E.inc.convex inc 7.740153 8.983706 2 1.131878 0.7204105 0.3802228 GP.bell bell 7.585780 8.903415 3 NA NA 0.1648041 L.dec dec 15.608419 14.889308 4 3.404153 NA 0.2142472 E.dec.concave dec 15.142111 14.367457 5 1.959459 -1.6121602 0.2994936 GP.U U 8.534547 9.341173 6 1.074078 NA 1.1263294 E.dec.concave dec 6.682254 3.082935 7 1.735801 -0.9436037 0.2594717 GP.U U 8.398699 9.007963 8 1.755034 0.7172904 0.2131311 GP.bell bell 11.664352 11.206771 9 2.557103 37.5472852 1.1492863 GP.bell bell 8.021110 5.546332 10 1.357555 NA 1.1925115 E.dec.concave dec 6.713592 3.338828 yrange maxychange xextrem yextrem 1 1.2435537 1.2435537 NA NA 2 1.5579862 1.5579862 1.438980 9.143766 3 0.7191107 0.7191107 NA NA 4 0.7746535 0.7746535 NA NA 5 1.1021236 0.8066266 1.959459 8.239050 6 3.5993189 3.5993189 NA NA 7 0.8289366 0.6092641 1.735801 8.179027 8 0.6684651 0.4575814 1.755034 11.875236 9 4.1812974 2.4747781 2.557103 9.727629 10 3.3747644 3.3747644 NA NA > f$unfitres [1] id irow adjpvalue cause <0 rows> (or 0-length row.names) > if (visualize) + { + plot(f) + # Alternative plots + # with a chosen number of first items + plot(f, items = 12) + # with chosen items in a specified order + plot(f, items = c("301.2", "363.1", "383.1")) + # residual plots + plot(f, items = 12, plot.type = "fitted_residuals") + plot(f, items = 12, plot.type = "dose_residuals") + # plot with dose in log + plot(f, items = 12, plot.type = "dose_fitted", dose_log_transfo = FALSE) + plot(f, items = 12, plot.type = "dose_residuals", dose_log_transfo = FALSE) + plot(f, items = 12, plot.type = "fitted_residuals", dose_log_transfo = FALSE) + } > > if (visualize) + { + # evaluate the impact of preventsfitsoutofrange + datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics") + (o1 <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess")) + (s_quad1 <- itemselect(o1, select.method = "quadratic", FDR = 0.001)) + + (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)]) + + (idchanged <- f1bis$fitres$id[which(f1bis$fitres$model != f1ter$fitres$model | + f1bis$fitres$f != f1ter$fitres$f)]) + # no impact on this dataset + } > > > > # calculation of benchmark doses > # options in shiny : z (numerical positive value), x (numerical positive value : percentage) > (r <- bmdcalc(f, z = 1, x = 10)) 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") + plot(r, plottype = "hist", by = "none", hist.bins = 10) + plot(r, plottype = "density", by = "none") + + plot(r, plottype = "hist", by = "trend", hist.bins = 10) + plot(r, plottype = "hist", by = "model", hist.bins = 10) + plot(r, plottype = "hist", by = "typology", hist.bins = 10) + } > > # 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.84 0.92 10.78