---
title: Overview of the DRomics package
author: Marie Laure Delignette Muller, Elise Billoir, Floriane Larras
date: '`r Sys.Date()`'
output:
  html_vignette:
  # pdf_document:
    toc: yes
    number_sections: yes

vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Overview of the DRomics package} 
  %!\VignetteEncoding{UTF-8}
  \usepackage[utf8]{inputenc}
---

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
# 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.

```{r}
# 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:

+ **RNAseqdata() should be used to import RNAseq as counts of reads**,
+ **microarraydata() should be used to import single-channel microarray data in log2 scale**,
+ **continuousomicdata() should be used to import other continuous omics data** such as metabolomics, proteomics,..., **in a scale that enables the use of a normal error model** 
in Steps 2 and 3. metabolomicdata() is the former name, but still available, of this function.

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 

```{r}
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.


```{r}
(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).


```{r}
microarrayfilename <- system.file("extdata", "transcripto_sample.txt", package = "DRomics")
```


```{r}
(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

```{r}
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.

```{r}
(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 

```{r}
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. 

```{r}
(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.

```{r}
(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). 

```{r}
(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) 

```{r}
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).

```{r}
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

```{r}
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).

```{r, fig.width = 7, fig.height = 5}
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}
(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).

```{r}
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.

```{r}
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.


```{r}
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).

```{r}
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).

```{r}
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
```{r}
(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).

```{r}
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.


```{r}
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).

```{r}
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.)

```{r}
# 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)
```

2. **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). 

```{r}
# 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)
```

3. **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.  

```{r}
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).

```{r}
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.)

```{r, eval = FALSE}
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).

```{r}
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. 

```{r}
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).

```{r}
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. 

```{r}
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).

```{r}
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).

```{r}
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](#enrich), below is an example of R code to import the DRomics results for microarray data, and to merge them with information on functional annotation.

```{r}
# 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](#enrich) metabolomics dataframe (extended results with functional annotation) is renamed for the sake of homogeneity.

```{r}
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).

```{r}
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.

```{r}
(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.

```{r}
(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.

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

```{r}
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).

```{r}
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.

```{r}
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).

```{r}
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).


```{r}
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.

```{r}
# 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).

```{r}
# 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).

```{r}
# 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 

+ Burnham, KP, Anderson DR (2004). Multimodel inference: understanding AIC and BIC in model selection. Sociological methods & research, 33(2), 261-304.

+ EFSA Scientific Committee, Hardy A, Benford D, Halldorsson T, Jeger MJ, Knutsen KH, ...
& Schlatter JR  (2017). Update: use of the benchmark dose approach in risk assessment.
EFSA Journal, 15(1), e04658.

+ Hurvich, CM, Tsai, CL (1989). Regression and time series model selection in small samples. Biometrika, 76(2), 297-307.

+ Larras F, Billoir E, Baillard V, Siberchicot A, Scholz S, Wubet T, Tarkka M, Schmitt-Jansen M and Delignette-Muller ML (2018). DRomics : a turnkey tool to support the use of the dose-response framework for omics data in ecological risk assessment. Environmental Science & Technology. [https://pubs.acs.org/doi/10.1021/acs.est.8b04752](https://pubs.acs.org/doi/10.1021/acs.est.8b04752).
You can also find this article at : [https://hal.archives-ouvertes.fr/hal-02309919](https://hal.archives-ouvertes.fr/hal-02309919)

+ Larras F, Billoir E, Scholz S, Tarkka M, Wubet T, Delignette-Muller ML, Schmitt-Jansen M (2020). A multi-omics concentration-response framework uncovers novel understanding of triclosan effects in the chlorophyte Scenedesmus vacuolatus. Journal of Hazardous Materials. [https://doi.org/10.1016/j.jhazmat.2020.122727](https://doi.org/10.1016/j.jhazmat.2020.122727).

