R Under development (unstable) (2024-01-29 r85841 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 RNAseq data and normalization > # with respect to library size and transformation > # options to put in shiny : transfo.method (2 methods, rlog or vst) > datafilename <- system.file("extdata", "RNAseq_sample.txt", package="DRomics") > # small data set 'less than 1000 items (999) > (o.vst <- RNAseqdata(datafilename, check = TRUE, transfo.method = "vst")) converting counts to integer mode 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.22 0.67 2 6 2 3 3 3 3 Number of items: 999 Identifiers of the first 20 items: [1] "NM_144958" "NR_102758" "NM_172405" "NM_029777" "NM_001130188" [6] "NM_207141" "NM_001162368" "NM_008117" "NM_001168290" "NM_010910" [11] "NM_001004147" "NM_001146318" "NM_145597" "NM_001161797" "NM_021483" [16] "NR_002862" "NR_033520" "NM_134027" "NM_010381" "NM_019388" Data were normalized with respect to library size and tranformed using the following method: vst Warning message: In RNAseqdata(datafilename, check = TRUE, transfo.method = "vst") : To optimize the dose-response modelling, it is recommended to use a dose-response design with at least six different tested doses. > plot(o.vst) > > if (visualize) # too long computation ! + { + plot(o.vst, range = 1.5) # boxplot visualizing outliers + (o.vst.notblind <- RNAseqdata(datafilename, check = TRUE, transfo.method = "vst", + transfo.blind = FALSE)) + plot(o.vst.notblind) + + (o.rlog <- RNAseqdata(datafilename, check = TRUE, transfo.method = "rlog")) + plot(o.rlog) + + (o.rlog.notblind <- RNAseqdata(datafilename, check = TRUE, transfo.method = "rlog", + transfo.blind = FALSE)) + plot(o.rlog.notblind) + + } > > if(visualize) # too long computation ! + { + data(Zhou_kidney_pce) + + # variance stabilizing tranformation + (o1 <- RNAseqdata(Zhou_kidney_pce, check = TRUE, transfo.method = "vst")) + plot(o1) + + # regularized logarithm + (o2 <- RNAseqdata(Zhou_kidney_pce, check = TRUE, transfo.method = "rlog")) + plot(o2) + + # variance stabilizing tranformation (not blind to the experimental design) + (o3 <- RNAseqdata(Zhou_kidney_pce, check = TRUE, transfo.method = "vst", + transfo.blind = FALSE)) + plot(o3) + + # regularized logarithm (not blind to the experimental design) + (o4 <- RNAseqdata(Zhou_kidney_pce, check = TRUE, transfo.method = "rlog", + transfo.blind = FALSE)) + plot(o4) + + } > > > if (visualize) + { + # item selection using the quadratic method + # options to put in shiny : select.method (3 methods), FDR (numerical positive value < 1) + o <- o.rlog + (s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01)) + (s_lin <- itemselect(o, select.method = "linear", FDR = 0.01)) + (s_ANOVA <- itemselect(o, select.method = "ANOVA", FDR = 0.01)) + + (f <- drcfit(s_quad, progressbar = TRUE)) + f$fitres + plot(f) + + # 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 == "bell", ], xmax = max(f$omicdata$dose), + facetby = "id") + + curvesplot(f$fitres[f$fitres$trend == "U", ], xmax = max(f$omicdata$dose), + facetby = "id") + } > > if (visualize) + { + # evaluate the impact of preventsfitsoutofrange, enablesfequal0inGP, enablesfequal0inlGP + data(Zhou_kidney_pce) + (o1 <- RNAseqdata(Zhou_kidney_pce, check = TRUE, transfo.method = "rlog")) + 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, dose_log_transfo = FALSE) + + (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, ] + + + } > > > # calculation of benchmark doses > # options in shiny : z (numerical positive value), x (numerical positive value : percentage) > if (visualize) + { + (r <- bmdcalc(f, z = 1, x = 10)) + (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) + 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 = "hist", by = "trend", hist.bins = 10) + } > > if (doboot) + { + # Calculation of confidence intervals on BMDs by Bootstrap + 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) + + } > > if(visualize) # too long computation ! + { + data(Zhou_kidney_pce) + + # exploration of an extreme case (BMD at 0) + d <- Zhou_kidney_pce + (o <- RNAseqdata(d)) + plot(o) + (s <- itemselect(o, select.method = "quadratic", FDR = 0.01)) + (f <- drcfit(s, progressbar = TRUE)) + head(f$fitres) + + r <- bmdcalc(f, z = 1) + plot(r) + + bmdplotwithgradient(r$res, BMDtype = "zSD") + bmdplotwithgradient(r$res, BMDtype = "zSD", BMD_log_transfo = FALSE) + + # no more 0 BMD values using argument minBMD + # res0 <- r$res[r$res$BMD.zSD == 0, ] + # curvesplot(res0, xmin =0.0000000001, xmax = max(f$omicdata$dose), + # colorby = "model", dose_log_transfo = TRUE) + # plot(f, items = r$res[r$res$BMD.zSD == 0, ]$id) + # plot(f, items = r$res[r$res$BMD.zSD == 0, ]$id, dose_log_trans = TRUE) + } > > proc.time() user system elapsed 11.23 1.07 12.31