Overview of the DRomics package

# output:
#  html_vignette
require(DRomics)
require(ggplot2)
set.seed(1234)
options(digits = 3)
knitr::opts_chunk$set(echo = TRUE,
                      eval = TRUE, 
                      message=FALSE, 
                      warning=FALSE,
                      cache=FALSE,
                      fig.width = 7, 
                      fig.height = 5)

Introduction

This vignette is intended to help users to start using the DRomics package. It is complementary to the reference manual where you can find more details on each function of the package. The first part of this vignette (Main workflow, steps 1 to 4) could also help users of the Shiny application. If you do not want to use R functions, you can skip the pieces of code and focus on the explanations and on the outputs that are also given in the Shiny application. And if one day you want go further using the R functions, we recommend you to start from the whole R code corresponding to your analysis that is provided on the last page of the Shiny application.

Main workflow

Step 1: importation, check and normalization / transformation of data if needed

General format of imported data

Whatever the type of data imported in DRomics, data can be imported from a .txt file (e.g. "mydata.txt") containing one row per item, after the first row coding for doses or concentrations for each sample (each column except the first one), with the first column corresponding to the identifier of each item (identifier of the probe, transcript, metabolite, ..., or name of the endpoint for anchoring data), and the other columns giving the responses of the item for each sample. In the first row, after a name for the identifier column, we must have the tested doses or concentrations in a numeric format for the corresponding sample (for example, if there are triplicates for each treatment, the first line could be "item", 0, 0, 0, 0.1, 0.1, 0.1, etc.). This file is imported within DRomics using the function read.table() with its default field separator (sep argument).

Alternatively an R object of class data.frame can be directly given in input, corresponding to the output of read.table(file, header = FALSE) on a file described as above.

You can see below an example of a RNAseq data set that is available in DRomics both as an R object (named Zhou_kidney_pce) and as a text file named "RNAseq_sample.txt" containing just a sample of the previous one.

# Load and look at the first line of the R object
data(Zhou_kidney_pce)
nrow(Zhou_kidney_pce)
head(Zhou_kidney_pce)

# Import the text file just to see what will be automatically imported
datafilename <- system.file("extdata", "RNAseq_sample.txt", package = "DRomics")
# for your local file datafilename choose the name: "yourchosenname.txt"
d <- read.table(file = datafilename, header = FALSE)
nrow(d)
head(d)

Types of data that may be imported in DRomics

DRomics offers the possibility to work on different types of omics data (see next paragraph for their description) but also on continuous anchoring data. When working on omics data, all the lines of the dataframe (except the first one coding for the doses or concentrations) correspond the same type of data (e.g. raw counts for RNAseq data). When working on anchoring data, the different lines (except the first one coding for the doses or concentrations) correspond to different endpoints that may correspond to different types of data (e.g. biomass, length,..), but all are assumed continuous data compatible with a normal error distribution (e.g. after logarithmic transformation for metabolomic data) for the selection and modelling steps.

Three types of omics data may be imported in DRomics using the following functions:

In Steps 1 and 2 count data are internally analysed using functions of the Bioconductor package DESeq2 while continuous data (microarray data and other continuous omics data) are internally analysed using functions of the Bioconductor package limma.

An example with RNAseq data

RNAseqfilename <- system.file("extdata", "RNAseq_sample.txt", package = "DRomics")

For RNAseq data, imperatively imported in raw counts, you have to choose the transformation method used to stabilize the variance ("rlog" or "vst"). In the example below "vst" was used only to make this vignette quick to compile, but "rlog" is strongly recommended and chosen by default even if more computer intensive than "vst" (see ?RNAseqdata for details). Whatever the chosen method, data are automatically normalized with respect to library size and transformed in a log2 scale.

(o.RNAseq <- RNAseqdata(RNAseqfilename, transfo.method = "vst"))
plot(o.RNAseq, cex.main = 0.8, col = "green", range = 1e6)

In the previous example the argument range (internally passed to boxplot) is put to a high value just to plot true minimum and maximum values and to prevent the automatic plot of many outliers as individual points in the plot of raw counts.

An example with microarray data

For single-channel microarray data, imperatively imported in log scale (classical and recommended log2 scale), you can choose the between array normalization method ("cyclicloess", "quantile", "scale" or "none"). In the example below "quantile" was used only to make this vignette quick to compile, but "cyclicloess" is strongly recommended and chosen by default even if more computer intensive than the others (see ?microarraydata for details).

microarrayfilename <- system.file("extdata", "transcripto_sample.txt", package = "DRomics")
(o.microarray <- microarraydata(microarrayfilename, norm.method = "quantile"))
plot(o.microarray, cex.main = 0.8, col = "green", range = 1e6)

In the previous example the argument range is intended to be internally passed by to the boxplot() function in order to enable long whiskers to be plotted, without individualizing extreme values, just to make the obtained figure lighter.

An example with metabolomic data

metabolofilename <- system.file("extdata", "metabolo_sample.txt", package = "DRomics")

No normalization nor transformation is provided in function continuousomicdata(). The pre-treatment of metabolomic data must be done before importation of data, and data must be imported in log scale, so that they can be directly modelled using a normal error model. This strong hypothesis is required both for selection of items and for dose-reponse modelling. In the context of a multi-omics approach we recommend you the use of a log2 transformation, instead of the classical log10 for such data, so as to facilitate the comparison of results obtained with transcriptomics data generally handled in a log2 scale.

As an example, a basic procedure for this pre-treatment of metabolomic data could follow the three steps described thereafter: i) removing of metabolites for which the proportion of missing data (non detections) across all the samples is too high (more than 20 to 50 percents according to your tolerance level); ii) retrieving of missing values data using half minimum method (i.e. half of the minimum value found for a metabolite across all samples); iii) log-transformation of values. If a scaling to the total intensity (normalization by sum of signals in each sample) or another normalization is necessary and pertinent, we recommend to do it before those three previously decribed steps.

(o.metabolo <- continuousomicdata(metabolofilename))
plot(o.metabolo, col = "green", range = 1e6)

We renamed metabolomicdata() to continuousomicdata() (while keeping the first name available) to offer its use to other continuous omic data such as proteomics data or RT-QPCR data. As for metabolomic data, the pretreatment of other continuous omic data data must be done before importation of data, and data must be imported in a scale that enables the use of a normal error model. This strong hypothesis is required both for selection of items and for dose-reponse modelling.

An example with continuous anchoring apical data

anchoringfilename <- system.file("extdata", "apical_anchoring.txt", package = "DRomics")

No transformation is provided in function continuousanchoringdata(). If needed the pretreatment of data must be done before importation of data, so that they can be directly modelled using a normal error model. This strong hypothesis is required both for selection of responsive endpoints and for dose-reponse modelling.

(o.anchoring <- continuousanchoringdata(anchoringfilename, backgrounddose = 0.1))
# Here the argument backgrounddose is used to specify that doses below  or equal to 0.1
# are considered as 0 in the DRomics workflow.
# Specifying this argument is necessary when there is no doses at 0 in the data 
plot(o.anchoring)

For such data the plot() function provides a dose-response plot for each endpoint.

Step 2: selection of significantly responding items

For the second step of the workflow, function itemselect() must be used with the output of the function used in Step 1 as first argument (output of RNAseqdata(), microarraydata(), continuousomicdata() or continuousanchoringdata()). Below is an example with microarray data.

The false discovery rate (FDR) corresponds to the expected proportion of items that will be falsely detected as responsive. With a very large data set it is important to define a selection step based on an FDR not only to reduce the number of items to be further processed, but also to remove too noisy dose-response signals that may impair the quality of the results. We recommend to set a value between 0.001 and 0.1 depending of the initial number of items. When this number is very high (more than several tens of thousands), we recommend a FDR less than 0.05 (0.001 to 0.01) to increase the robustness of the results (Larras et al. 2018).

Concerning the method used for selection, we recommend the default choice ("quadratic") for a typical omics dose-response design (many doses/concentrations with few replicates per condition). It enables the selection of both monotonic and biphasic dose-responses. If you want to focus on monotonic dose-responses, the "linear" method could be chosen. For a design with a small number of doses/concentrations and many replicates (not an optimal for dose-response modelling), the "ANOVA" method could be preferable.

See ?itemselect and Larras et al. 2018 for details.

(s_quad <- itemselect(o.microarray, select.method = "quadratic", FDR = 0.01))

Step 3: fit of dose-response models, choice of the best fit for each curve

Fit

For Step 3 the function drcfit() must be simply used with the output of itemselect() as first argument. Description of the fitted models and of the procedure to select the best fit are described in Larras et al. 2018 and in ?drcfit. The former use of the AIC (Akaike criterion- default information criterion used for the selection of the best fit model in DRomics versions < 2.2-0) was replaced by the use of the AICc (second-order Akaike criterion) in order to prevent the overfitting that may occur with dose-response designs with a small number of data points, as recommended and now classically done in regression (Hurvich and Tsai, 1989; Burnham and Anderson DR, 2004).

As the call to this function may take time, by default a progressbar is provided. Some arguments of this function can be used to specify parallel computing to accelerate the computation (see ?drcfit for details).

(f <- drcfit(s_quad, progressbar = FALSE))

In the following you can see the first ten lines of the output dataframe on our example (see ?drcfit for a complete description of the columns of the output dataframe.) This output dataframe provides information such as best-fit model, parameter value, coordinates of particular points, and the trend of the curve (among increasing, decreasing, U-shaped, bell-shaped)

head(f$fitres, 10)

Plot of fitted curves

By default the plot() function used on the output of the drcfit() function provides the first 20 fitted curves (or the ones you specify using the argument items) with observed points. Fitted curves are represented in red, replicates are represented in open circles and means of replicates at each dose/concentration are represented by solid circles. All the fitted curves may be saved in a pdf file using the plotfit2pdf() function (see ?drcfit).

plot(f) 

The fitted curves may be represented using a log scale for the dose/concentration using argument dose_log_transfo (see ?drcfit for details and examples).

Another specific plot function named targetplot() can be used to plot targeted items, whether they were or not selected in step 2 and fitted in step 3. See an example below and details in ?targetplot

targetitems <- c("88.1", "1", "3", "15")
targetplot(targetitems, f = f)

Plot of residuals

To check the assumption of the normal error model, two types of residual plots can be used ("dose_residuals" or "fitted_residuals"). The residual plots for all items may also be saved in a pdf file using the plotfit2pdf() function (see ?drcfit).

plot(f, plot.type = "dose_residuals")

Step 4: calculation of x-fold and z-SD benchmark doses

Calculation of BMD

The two types of benchmark doses (BMD-zSD and BMD-xfold) proposed by the EFSA (2017) are systematically calculated for each fitted dose-response curve using the function bmdcalc() with the output of the drcfit() function as a first argument (see Larras et al. 2018 or ?drcfit for details).

The argument z, by default at 1, is used to define the BMD-zSD as the dose at which the response is reaching y0 +/- z * SD, with y0 the level at the control given by the dose-response fitted model and SD the residual standard deviation of the dose-response fitted model.

The argument x, by default at 10 (for 10%), is used to define the BMD-xfold as the dose at which the response is reaching y0 +/- (x/100) * y0.

(r <- bmdcalc(f, z = 1, x = 10))

In the following you can see the first ten lines of the output dataframe of the function bmdcalc() on our example (see ?bmdcalc for a complete description of the columns of the output dataframe).

head(r$res, 10)

Various plots of the BMD distribution

The default plot of the output of the bmdcalc() function provides the distribution of benchmark doses as an ECDF (Empirical Cumulative Density Function) plot for the chosen BMD ("zSD"" or "xfold"). See an example below.

plot(r, BMDtype = "zSD", plottype = "ecdf") 

Different alternative plots are proposed (see ?bmdcalc for details) that can be obtained using the argument plottype to choose the type of plot ("ecdf", "hist" or "density") and the argument by to split the plot by "trend", "model" or "typology". Below is an example of a density plot of BMD-zSD split by trend of dose-response curves.

plot(r, BMDtype = "zSD", plottype = "density", by = "trend") 

Plot of BMD distribution with a color gradient for signal intensity

On a BMD ECDF plot one can add of a color gradient for each item coding for the intensity of the signal (after shift of the control signal at 0) as a function of the dose (see ?bmdplotwithgradient for details and an example below).

bmdplotwithgradient(r$res, BMDtype = "zSD",
                    facetby = "trend", 
                    shapeby = "model",
                    line.size = 1.2) + labs(shape = "model")

As in the previous example, you can use the argument line.size to manually adjust the width of lines in that plot if the default value does not give a visual result that suits you.

Plot of fitted curves with BMD values

It is possible to add the output of bmdcalc() in the argument BMDoutput of the plot() function of drcfit() to add BMD values (when defined) as a vertical line on each fitted curve. Horizontal dotted lines corresponding to the two BMR potential values will be also added. All the fitted curves may also be saved in the same way in a pdf file using the plotfit2pdf() function (see ?drcfit).

plot(f, BMDoutput = r) 

Step 5: calculation of confidence intervals on the BMDs by bootstrap

Confidence intervals on BMD values can be calculated by bootstrap. As the call to this function may take much time, by default a progressbar is provided and some arguments can be used to specify parallel computing to accelerate the computation (see ?bmdboot for details).

In the example below a small number of iterations was used just to make this vignette quick to compile, but the default value of the argument niter (1000) should be considered as a minimal value to obtain stable results.

Bootstrap calculation

(b <- bmdboot(r, niter = 50, progressbar = FALSE))

This function gives an output corresponding to the output of the bmdcalc() function completed with bounds of BMD confidence intervals (by default 95% confidence intervals).

head(b$res, 10)

Add of confidence intervals on BMD ECDF plots

The plot() function applied on the output the bmdboot() function gives an ECDF plot of the chosen BMD with the confidence interval of each BMD (see an example below). By default BMDs with an infinite confidence interval bound are not plotted.

plot(b, BMDtype = "zSD", by = "trend") 

Plot of fitted curves with BMD values and confidence intervals

It is possible to add the output of bmdboot() in the argument BMDoutput of the plot() function of drcfit() to add BMD values (when defined) as a vertical line on each fitted curve and bounds of the confidence intervals (when successfully calculated) as two dashed lines. Horizontal dotted lines corresponding to the two BMR potential values will be also added. All the fitted curves may also be saved in the same way in a pdf file using the plotfit2pdf() function (see ?drcfit).

plot(f, BMDoutput = b) 

Help for biological interpretation of DRomics outputs

This section illustrates functions of DRomics that are meant to help the biological interpretation of outputs. The idea is to augment the output dataframe with new column(s) bringing biological information such as provided by functional annotation of the items (e.g. KEGG pathway classes or GO terms) then to use this information to organize the visualisation of the DRomics output.

Below is used an example from a metabolomic data set previously analysed using DRomics.

Enrichment of the dataframe of DRomics results with functional annotation {#enrich}

This enrichment is not done using DRomics functions, but with relevant R functions such as merge().

An example of how to proceed:

  1. Import the dataframe with DRomics results to be used: the output $res of bmdcalc() or bmdboot() functions from step 4 or 5 of the main DRomics workflow.

(This step will not be necessary if previous steps are done directly in R using the DRomics package as described previously in this vignette. We did it to take a real example that took a long time to run but from which results are stored in the package.)

# code to import the file for this example in our package
resfilename <- system.file("extdata", "triclosanSVmetabres.txt", package = "DRomics")
res <- read.table(resfilename, header = TRUE, stringsAsFactors = TRUE)

# to see the structure of this file
str(res)
  1. Import the dataframe with functional annotation (or any other descriptor/category you want to use, here KEGG pathway classes) of each item present in the 'res' file.

Examples are embedded in the DRomics package, but be cautious, generally this file must be produced by the user. Each item may have more than one annotation (i.e. more than one line).

# code to import the file for this example in our package
annotfilename <- system.file("extdata", "triclosanSVmetabannot.txt", package = "DRomics")
annot <- read.table(annotfilename, header = TRUE, stringsAsFactors = TRUE)

# to see the structure of this file
str(annot)
  1. Merging of both previous dataframes in order to obtain a so-called 'extendedres' dataframe gathering, for each item, metrics derived from the DRomics workflow and functional annotation.

Arguments by.x and by.y of the merge() function indicate the column name in res and annot dataframes, respectively, that must be used for the merging.

annotres <- merge(x = res, y = annot, by.x = "id", by.y = "metab.code")
head(annotres)

BMD ECDF plot by functional group

BMD ECDF plot split by group defined from functional annotation

Using the function bmdplot() and its argument facetby, the BMD ECDF plot can be split by group (here by KEGG pathway class). Confidence intervals can be added on this plot and color coding for trend in this example (See ?bmdplot for more options).

bmdplot(annotres, BMDtype = "zSD", add.CI = TRUE, 
                    facetby = "path_class", 
                    colorby = "trend") + labs(col = "trend")

Function ecdfplotwithCI() can also used as an alternative as below to provide the same plot differing just by the coloring of intervals only. (See ?ecdfplotwithCI for more options.)

ecdfplotwithCI(variable = annotres$BMD.zSD, 
               CI.lower = annotres$BMD.zSD.lower, 
               CI.upper = annotres$BMD.zSD.upper, 
               by = annotres$path_class,
               CI.col = annotres$trend) + labs(col = "trend")

BMD ECDF plot with color gradient split by group defined from functional annotation

Using the function bmdplotwithgradient() and its argument facetby, the BMD plot with color gradient can be split here by KEGG pathway class. (See ?bmdplotwithgradient for more options).

bmdplotwithgradient(annotres, BMDtype = "zSD",
                    facetby = "path_class", 
                    shapeby = "trend") + labs(shape = "trend")

the same representation with labels of items (so without shapeby)

The argument add.label set at TRUE will display item identifiers instead of points.

bmdplotwithgradient(annotres, BMDtype = "zSD",
                    facetby = "path_class", 
                    add.label = TRUE) +
  theme(strip.text.x = element_text(size = 6))

To increase the visibility of this plot, you can plot it one by one for each group of interest as below for the group "Lipid metabolism". In that case in can be useful to control the limits of the color gradient and the limits on the x-axis in order to use the same x-scale and signal-scale, as in the following example (see ?bmdplotwithgradient for details).

annotres_lipid <- annotres[annotres$path_class == "Lipid metabolism",] 

bmdplotwithgradient(annotres_lipid, BMDtype = "zSD",
                    facetby = "path_class", 
                    add.label = TRUE,
                    limits4colgradient = c(-0.8, 0.8),
                    xmin = 0, xmax = 6.5,
                    label.size = 3) 

Sensitivity plot of functional groups

It is also possible to show a summary of BMD values in each pathway/category as a given summary (argument BMDsummary which can be "first.quartile", "median" or "median.and.IQR" for medians with the interquartile range as an interval) using the function sensitivityplot(). Moreover, this function will provide information on the number of items involved in each pathway/category. (See ?sensitivityplot for more options).

As an example, below is an ECDF plot of 25th quantiles of BMD-zSD calculated here by pathway class.

sensitivityplot(annotres, BMDtype = "zSD",
                    group = "path_class",
                 BMDsummary = "first.quartile") 

Trend plot per functional group

It is possible to represent the repartition of trends in each functional group using function trendplot() (see ?trendplot for details).

trendplot(annotres, group = "path_class") 

Plot of dose-response curves for a specific functional group

The function curvesplot() can show the dose-response curves for a specific pathway/category with arguments left to the choice of the user. In this example, only results related to the "lipid metabolism" pathway class are kept. The plot is split by id (argument facetby) and colored by trend (argument colorby). (See ?curvesplot for more options).

LMres <- annotres[annotres$path_class == "Lipid metabolism", ]
curvesplot(LMres, facetby = "id", npoints = 100, line.size = 1,
           colorby = "trend",
           xmin = 0, xmax = 8) + labs(col = "trend")

Help for multi-omics approaches

This section illustrates functions of DRomics that are meant to help the interpretation of outputs by linking several omics levels. The idea is to augment the output dataframe with new column(s) bringing information on the molecular level, then to use this information to organize the visualisation of the DRomics output.

Below is used an example linking a transcriptomic (microarray) and a metabolomic data set issued from the same experiment.

An example with metabolomics and transcriptomics data for Scenedesmus and triclosan (cf. Larras et al. 2020)

Enrichment of the dataframes of DRomics results with functional annotation

Following the same steps as described before for metabolomics, below is an example of R code to import the DRomics results for microarray data, and to merge them with information on functional annotation.

# 1. Import the dataframe with DRomics results to be used
contigresfilename <- system.file("extdata", "triclosanSVcontigres.txt", package = "DRomics")
contigres <- read.table(contigresfilename, header = TRUE, stringsAsFactors = TRUE)
str(contigres)

# 2. Import the dataframe with functional annotation (or any other descriptor/category 
# you want to use, here KEGG pathway classes) 
contigannotfilename <- system.file("extdata", "triclosanSVcontigannot.txt", package = "DRomics")
contigannot <- read.table(contigannotfilename, header = TRUE, stringsAsFactors = TRUE)
str(contigannot)

# 3. Merging of both previous dataframes   
contigextendedres <- merge(x = contigres, y = contigannot, by.x = "id", by.y = "contig")
# to see the structure of this dataframe
str(contigextendedres)

The previouly created metabolomics dataframe (extended results with functional annotation) is renamed for the sake of homogeneity.

metabextendedres <- annotres

Binding of transcriptomics and metabolomics dataframes

The next step is the binding of dataframes at both levels adding a variable (named level) coding for the level (here a factor with two levels, metabolites and contigs).

extendedres <- rbind(metabextendedres, contigextendedres)
extendedres$level <- factor(c(rep("metabolites", nrow(metabextendedres)),
                              rep("contigs", nrow(contigextendedres))))
str(extendedres)

Comparison of results obtained at different molecular levels using basic R functions

Below are examples of illustrations that can be used to compare the results obtained at several levels of biological organization using basic R functions.

Frequencies of pathways per molecular level

Here basic R functions are used to compute and plot frequencies of pathways by molecular levels.

(t.pathways <- table(extendedres$path_class, extendedres$level)) 
original.par <- par()
par(las = 2, mar = c(4,13,1,1))
barplot(t(t.pathways), beside = TRUE, horiz = TRUE, 
        cex.names = 0.7, legend.text = TRUE, 
        main = "Frequencies of pathways")
par(original.par)

Proportions of pathways per molecular level

Here basic R functions are used to compute and plot proportions of pathways by molecular levels.

(t.prop.pathways <- prop.table(t.pathways, margin = 2)) 
original.par <- par()
par(las = 2, mar = c(4,13,1,1))
barplot(t(t.prop.pathways), beside = TRUE, horiz = TRUE, 
        cex.names = 0.7, legend.text = TRUE, 
        main = "Proportion of pathways")
par(original.par)

ECDF plot of BMD_zSD by pathway using different colors or facets for molecular levels

Here the ggplot2 grammar is used to plot the ECDF of BMD_zSD using different colors or facets for the different molecular levels.

if (require(ggplot2))
{
   ggplot(extendedres, aes(x = BMD.zSD, color = level)) +
      stat_ecdf(geom = "step") + ylab("ECDF")

}
if (require(ggplot2))
{
   ggplot(extendedres, aes(x = BMD.zSD)) + facet_wrap(~ level) +
      stat_ecdf(geom = "step") + ylab("ECDF")

}

Comparison of results obtained at different molecular levels using DRomics functions

With DRomics functions bmdplot(), bmdplotwithgradient(), sensitivityplot(), trendplot() and curvesplot(), it is also possible to easily integrate different molecular levels (or another grouping level).

ECDF plot of BMD_zSD per pathway and molecular level

Using the function bmdplot() and its argument facetby, the BMD ECDF plot can be split by group (here by KEGG pathway class) and colored or split by molecular level. (See ?bmdplot for more options).

bmdplot(extendedres, BMDtype = "zSD", 
                    facetby = "path_class", 
                    colorby = "level") + labs(col = "molecular level")

Below is a focus on few pathways with split on pathway and molecular level.

chosen_path_class <- c("Membrane transport", "Lipid metabolism", "Energy metabolism")
selectedres <- extendedres[extendedres$path_class %in% chosen_path_class, ]
bmdplot(selectedres, BMDtype = "zSD", add.CI = TRUE, 
                    facetby = "path_class", 
                    facetby2 = "level") + labs(col = "trend")

BMD ECDF plot with color gradient split by pathway and molecular level

Using the function bmdplotwithgradient() and its arguments facetby and facetby2, the BMD plot with color gradient can be split here by KEGG pathway class and molecular level, here on a selection of pathway classes present at both molecular levels.(See ?bmdplotwithgradient for more options). For a better visualization, it should have been preferable to enter metabolic data in DRomics in log2 instead of log10, to be in same signal scale as for transcriptomic data (indeed the amplitude of the metabolic signal is lower due to the use the log10).

chosen_path_class <- c("Membrane transport", "Lipid metabolism", "Energy metabolism")
selectedres <- extendedres[extendedres$path_class %in% chosen_path_class, ]
bmdplotwithgradient(selectedres, BMDtype = "zSD",
               facetby = "path_class", facetby2 = "level")

Sensitivity plot per pathway and molecular level

Using the function sensitivityplot() and its arguments group and colorby, it is possible to show a summary of BMD values with size of points coding for the number of items in each group as in the example, where the 25th quartiles of BMD values are represented per KEGG pathway class for each molecular level. (See ?sensitivityplot for more options).

sensitivityplot(extendedres, BMDtype = "zSD",
                group = "path_class", colorby = "level",
                BMDsummary = "first.quartile") 

Plot of the trend repartition per pathway and molecular level

Using the function trendplot() and its arguments facetby it is possible to show the repartition of trends of responses in each pathway class for both molecular levels.

# Preliminary optional alphabetic ordering of path_class groups
extendedres$path_class <- factor(extendedres$path_class, 
                levels = sort(levels(extendedres$path_class), decreasing = TRUE))
trendplot(extendedres, group = "path_class", facetby = "level") 

Plot of the dose-response curves for a specific metabolic pathway

Using the function curvesplot(), specific dose-response curves can be shown. In this example, first only results related to the "lipid metabolism" pathway class are kept. Then, the plot is split by molecular level (argument facetby) and colored by trend (argument colorby). (See ?curvesplot for more options).

# Plot of the dose-response curves for a specific metabolic pathway
# in this example the "lipid metabolism" pathclass
LMres <- extendedres[extendedres$path_class == "Lipid metabolism", ]
curvesplot(LMres, facetby = "level", free.y.scales = TRUE, npoints = 100, line.size = 1,
           colorby = "trend",
           xmin = 0, xmax = 8) + labs(col = "DR_trend")

Plot of dose-response curves for a selection of pathways

The same plot can be drawn for several chosen pathways with the use of argument facetby2. (See ?curvesplot for more options).

# Plot of the dose-response curves for a specific metabolic pathway
# in this example the "lipid metabolism" pathclass
curvesplot(selectedres, facetby = "path_class", facetby2 = "level",
           free.y.scales = TRUE, npoints = 100, line.size = 1,
           colorby = "trend",
           xmin = 0, xmax = 8) + labs(col = "DR_trend")

References



Try the DRomics package in your browser

Any scripts or data that you put into this service are public.

DRomics documentation built on Jan. 6, 2022, 5:07 p.m.