---
title: "Running `MultiRNAflow` on a RNA-seq raw counts with different time points and several biological conditions"
shorttitle: "MultiRNAflow"
author:
- name: "Rodolphe Loubaton"
affiliation: "Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France"
- name: "Nicolas Champagnat"
affiliation: "Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France"
- name: "Pierre Vallois"
affiliation: "Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France"
- name: "Laurent Vallat"
affiliation:
- "Université de Strasbourg, CNRS, UMR-7242 Biotechnology and cell signaling, F-67400 Illkirch, France"
- "Department of molecular genetic of cancers, Strasbourg University Hospital, F-67200 Strasbourg, France"
date: "`r format(Sys.Date(), '%d/%m/%Y')`"
output:
BiocStyle::html_document:
toc: true
number_sections: true
toc-title: "Contents"
bibliography: MultiRNAflowBiblio.bib
vignette: >
%\VignetteIndexEntry{Running_analysis_with_MultiRNAflow}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE, results="hide"}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
tidy = FALSE,
cache = FALSE,
#dev = "cairo_pdf",
message = FALSE, error = FALSE, warning = FALSE
)
```
# Introduction
## Quick description of the document
This document is a quick workflow describing how to use our R package
MultiRNAflow on one dataset (see [Dataset used as example]).\newline
For a more complete description of our package and complete outputs with graphs,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
## Dataset used as example
The **Mouse dataset 2** [@Weger2021TemporalMouseLiver] is accessible on the
Gene Expression Omnibus (GEO) database with the accession number GSE135898.
This dataset contains the temporal transcription profile of 16 mice with Bmal1
and Cry1/2 knocked-down under an ad libitum (AL) or night restricted feeding
(RF) regimen. Data were collected at 0, 4h, 8h, 16, 20h and 24h.
Therefore, there are six time points and eight biological conditions.
As there are only two mice per biological condition, we decided not to take
into account the effect of the regimen.
The dataset contains temporal expression data of 40327 genes. \newline
To illustrate the use of our package, we take 500 genes, over the global 40327
genes in the original dataset.
This sub dataset is saved in the file **RawCounts_Weger2021_MOUSEsub500**.
# Quick workflow
## Load package, example dataset and preamble
1. Load MultiRNAflow
```{r library, warning=FALSE, message=FALSE}
library(MultiRNAflow)
```
2. Load **Mouse dataset 2**
```{r DataWeger2021}
data("RawCounts_Weger2021_MOUSEsub500")
```
3. Structure of the dataset (Preamble)
The dataset must be a data.frame containing raw counts data.
If it is not the case, the functions `DATAnormalization()` and
`DEanalysisGlobal()` will stop and print an error.
Each line should correspond to a gene, each column to a sample,
except a particular column that may contain strings of characters
describing the names of the genes.
The first line of the data.frame should contain the names of the columns
(strings of characters) that must have the following structure.
```{r ColnamesDataWeger2021}
head(colnames(RawCounts_Weger2021_MOUSEsub500), n=5)
```
In this example, "Gene" indicates the column which contains the names of
the different genes.
The other column names contain all kind of information about the sample,
including the biological condition, the time of measurement and the name of
the individual (e.g patient, replicate, mouse, yeasts culture...).
Other kinds of information can be stored in the column names
(such as patient information), but they will not be used by the package.
The various information in the column names must be separated by underscores.
The order of these information is arbitrary but must be the same
for all columns. For instance, the sample "BmKo_t0_r1" corresponds to the first
replicate (r1) of the biological condition BmKo at time t0 (reference time).
The information located to the left of the first underscore will be considered
to be in position 1, the information located between the first underscore and
the second one will be considered to be in position 2, and so on.
In the previous example, the biological condition is in position 1,
the time is in position 2 and the replicate is in position 3.
In most of the functions of our package, the order of the previous information
in the column names will be indicated with the inputs `Group.position`,
`Time.position` and `Individual.position`.
Similarly the input `Column.gene` will indicate the number of
the column containing gene names.
## Preprocessing
4. Preprocessing with `DATAprepSE()`
```{r PreprocessingMouse2, eval=TRUE}
resDATAprepSE <- DATAprepSE(RawCounts=RawCounts_Weger2021_MOUSEsub500,
Column.gene=1,
Group.position=1,
Time.position=2,
Individual.position=3)
```
* **Output**: The function realizes the normalization step and returns
a SummarizeExperiment object containing
- all information about the raw counts data
- a DESeqDataSet object to be used by the function `DEanalysisGlobal`
for the statistical (supervised) analysis.
* **Input**:
- The dataset must be RNA-seq raw counts (`RawCounts`)
- The argument `Column.gene=1` means that the first column of the dataset
contain genes name, `Time.position=2` means that the time of measurements
is between the first and the second underscores in the columns names,
`Individual.position=3` means that the name of the individual is between
the second and the third underscores in the columns names and
`Group.position=NULL` means that there is only one biological condition
in the dataset. Similarly, `Time.position=NULL` would mean that
there is only one time of measurement for each individual and
`Column.gene=NULL` would mean that there is no column containing gene names.
* **Other**:
- Write `?DATAprepSE` in your console for more information about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled "MultiRNAflow_vignette-knitr.pdf".
## Exploratory data analysis
5. Normalization with `DATAnormalization()`
```{r NormalizationMouse2, eval=TRUE}
resNorm <- DATAnormalization(SEres=resDATAprepSE,
Normalization="vst",
Blind.rlog.vst=FALSE,
Plot.Boxplot=FALSE,
Colored.By.Factors=TRUE,
Color.Group=NULL,
path.result=NULL)
```
* **Output**: The function realizes the normalization step and
- returns the same SummarizedExperiment class object `resDATAprepSE` but with
the normalized data
- plots a boxplot (if `Plot.Boxplot=TRUE`) showing the distribution of
the normalized expression (`Normalization="vst"` means that the vst method
is used) of genes for each sample.
* **Input**:
- The results of the function `DATAprepSE()`.
- In order to display the output graph, set `Plot.Boxplot=TRUE`.
- If `Colored.By.Factors=TRUE`, the color of the boxplots would be different
for different biological conditions.
- In order to save the different results in a folder, select a folder path
in `path.result`.
* **Other**:
- Write `?DATAnormalization` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled "MultiRNAflow_vignette-knitr.pdf".
6. Principal component analysis (PCA) with `PCAanalysis()`
```{r PCAMus2, eval=TRUE}
resPCA <- PCAanalysis(SEresNorm=resNorm,
DATAnorm=TRUE,
gene.deletion=NULL,
sample.deletion=NULL,
Plot.PCA=FALSE,
Mean.Accross.Time=FALSE,
Color.Group=NULL,
Cex.label=0.6,
Cex.point=0.7, epsilon=0.2,
Phi=25,Theta=140,
motion3D=FALSE,
path.result=NULL)
```
* **Output**: When samples belong to different biological conditions and
different time points, the previous lines of code return
- the same SummarizedExperiment class object `resNorm` but with
the results of the function `PCA()` of the package `FactoMineR`.
- One 2D PCA graph, one 3D PCA graph and the same 3D PCA graph in a rgl window
(only if `motion3D=FALSE`) where samples are colored with different colors
for different biological conditions. Furthermore, lines are drawn
between each pair of consecutive points for each sample
(if `Mean.Accross.Time=FALSE`, otherwise it will be only between means).
- One 2D PCA graph, one 3D PCA graph and the same 3D PCA graph in a
rgl window (only if `motion3D=FALSE`) for each biological condition,
where samples are colored with different colors for different time points.
Furthermore, lines are drawn between each pair of consecutive points
for each sample (if `Mean.Accross.Time=FALSE`, otherwise it will be only
between means).
- The same graphs describe above but without lines.
* **Input**:
- The results of the function `DATAnormalization()`.
- We recommend the use of the normalized data (`DATAnorm=TRUE`) for the PCA
analysis.
- By default (if `Color.Group=NULL`), a color will be automatically applied
for each biological condition. If you want to change the colors,
read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
- If you want to delete, for instance, the genes ’ENSMUSG00000025921’ and
’ENSMUSG00000026113’ (respectively the second and sixth gene) and/or delete
the samples 'BmKo\_t2\_r1' and 'BmKo\_t5\_r2', set
- `gene.deletion=c("ENSMUSG00000025921", "ENSMUSG00000026113")` and/or
`sample.deletion=c("BmKo_t2_r1", "BmKo_t5_r2")`
- `gene.deletion=c(2,6)` and/or `sample.deletion=c(3,13)`.
The integers in `gene.deletion` and `sample.deletion` represent
respectively the row numbers and the column numbers of `ExprData`
where the selected genes and samples are located.
- In order to display the different output graph, set `Plot.PCA=TRUE`.
- In order to save the different results in a folder,
select a folder path in `path.result`.
* **Other**:
- Write `?PCAanalysis` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled "MultiRNAflow_vignette-knitr.pdf".
7. Hierarchical Clustering on Principle Components (HCPC) with `HCPCanalysis()`
```{r HCPCmus2500, warning=FALSE, message=FALSE, eval=TRUE}
resHCPC <- HCPCanalysis(SEresNorm=resNorm,
DATAnorm=TRUE,
gene.deletion=NULL,
sample.deletion=NULL,
Plot.HCPC=FALSE,
Phi=25,Theta=140,
Cex.point=0.6,
epsilon=0.2,
Cex.label=0.6,
motion3D=FALSE,
path.result=NULL)
```
* **Output**:
- the same SummarizedExperiment class object `resNorm` but with
the results of the function `HCPCanalysis()` of the package `FactoMineR`
- A dendrogram to illustrate how each cluster is composed
- A graph indicating by color for each sample, its cluster,
the biological condition and the time point associated to the sample.
- One 2D PCA graph, one 3D PCA graph and the same 3D PCA graph in a rgl window
(only if `motion3D=FALSE`). These PCA graphs are identical to the outputs
of `PCAanalysis()` with all samples but samples are colored with different
colors for different clusters.
* **Input**:
- The results of the function `DATAnormalization()`.
- We recommend the use of the normalized data (`DATAnorm=TRUE`)
for the HCPC analysis.
- In order to display the different output graph, set `Plot.HCPC=TRUE`.
- In order to save the different results in a folder, select a folder path in `path.result`.
* **Other**:
- The optimal number of clusters is automatically computed by the
R function `HCPC()`.
- Write `?HCPCanalysis` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled "MultiRNAflow_vignette-knitr.pdf".
8. Temporal clustering analysis with `MFUZZanalysis()`.
```{r MfuzzMus2, eval=TRUE, results='hide', fig.show='hide'}
resMFUZZ <- MFUZZanalysis(SEresNorm=resNorm,
DATAnorm=TRUE,
DataNumberCluster=NULL,
Method="hcpc",
Membership=0.5,
Min.std=0.1,
Plot.Mfuzz=FALSE,
path.result=NULL)
```
* **Output**:
- the same SummarizedExperiment class object `resNorm` but with
+ The summary of the results of the function `mfuzz()`
+ The number of clusters automatically computed
(if `DataNumberCluster=NULL`). If `Method="hcpc"`, the function plots
the scaled within-cluster inertia, but if `Method="kmeans"`,
the function plots the scaled within-cluster inertia.
As the number of genes can be very high, we recommend to select
`Method="hcpc"` which is by default.
- the output graph from the function `mfuzz()` which represents the most
common temporal behavior among all genes for the biological condition ‘BmKo’.
* **Input**:
- The results of the function `DATAnormalization()`.
- We recommend the use of the normalized data (`DATAnorm=TRUE`)
for the clustering analysis.
- For each cluster, genes with membership values below the threshold
`Membership` (numeric value between 0 and 1) will not be displayed.
The membership values correspond to the probability of gene to belong
to each cluster.
- All genes where their standard deviations are smaller than the threshold
`Min.std` (numeric positive value) will be excluded.
- In order to display the different output graph, set `Plot.Mfuzz=TRUE`.
- In order to save the different results in a folder, select a folder path in `path.result`.
* **Other**:
- Write `?MFUZZanalysis` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled "MultiRNAflow_vignette-knitr.pdf".
9. Plot temporal expression with with `DATAplotExpressionGenes()`
```{r DATAprofileGenesMusBmCrKoWt, eval=TRUE, results='hide', fig.show='hide'}
resEXPR <- DATAplotExpressionGenes(SEresNorm=resNorm,
DATAnorm=TRUE,
Vector.row.gene=c(17),
Color.Group=NULL,
Plot.Expression=FALSE,
path.result=NULL)
```
* **Output**: A graph plotting for each biological condition:
the evolution of the 17th gene expression of the three replicates across time
and the evolution of the mean and the standard deviation of the 17th gene
expression across time. The color of the different lines are different
for different biological conditions.
* **Input**:
- The results of the function `DATAnormalization()`.
- We recommend the use of the normalized data (`DATAnorm=TRUE`).
- If the user wants to select several genes, for instance the 97th, the 192th,
the 194th and the 494th, he needs to set `Vector.row.gene=c(97,192,194,494)`.
- By default (if `Color.Group=NULL`), a color will be automatically applied
for each biological condition. If you want to change the colors,
read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
- In order to display the different output graph, set `Plot.Expression=TRUE`.
- In order to save the different results in a folder,
select a folder path in `path.result`.
* **Other**:
- Write `?DATAplotExpressionGenes` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
## Supervised statistical analysis
10. Differential Expresion (DE) analysis with `DEanalysisGlobal()`
For the speed of the algorithm, we will take only three biological conditions
and three times
```{r SubMusBmCrKoWt}
Sub3bc3T <- RawCounts_Weger2021_MOUSEsub500
Sub3bc3T <- Sub3bc3T[, seq_len(73)]
SelectTime <- grep("_t0_", colnames(Sub3bc3T))
SelectTime <- c(SelectTime, grep("_t1_", colnames(Sub3bc3T)))
SelectTime <- c(SelectTime, grep("_t2_", colnames(Sub3bc3T)))
Sub3bc3T <- Sub3bc3T[, c(1, SelectTime)]
resSEsub <- DATAprepSE(RawCounts=Sub3bc3T,
Column.gene=1,
Group.position=1,
Time.position=2,
Individual.position=3)
```
```{r DEMusBmCrKoWt, eval=TRUE, echo=TRUE}
resDE <- DEanalysisGlobal(SEres=resSEsub,
pval.min=0.05,
log.FC.min=1,
Plot.DE.graph=FALSE,
path.result=NULL)
```
* **Output**:
- The results of the function `DESeq()`
- a data.frame (output `Results`) which contains
- gene names
- pvalues, log2 fold change and DE genes between each pairs of
biological conditions for each fixed time.
- pvalues, log2 fold change and DE genes between each time versus the
reference time t0 for each biological condition.
- Temporal pattern (succession of 0 and 1) for each biological condition.
The positions of 1 in one of these two columns, indicate the set of times ti
such that the gene is DE between ti and the reference time t0,
for the biological condition associated to the given column.
- Specific genes for each biological condition at each time.
A 1 in one of these columns means the gene is specific to the biological
condition at a fixed time associated to the given column. 0 otherwise.
A gene is called specific to a given biological condition BC1 at a time ti,
if the gene is DE between BC1 and any other biological conditions
at time ti, but not DE between any pair of other biological conditions
at time ti.
- Signature genes the signatures genes for each biological condition
at each time. A 1 in one of these columns means the gene is signature gene
to the biological condition at a fixed time associated to the given column.
0 otherwise. A gene is called signature of a biological condition BC1
at a given time ti, if the gene is specific to the biological condition BC1
at time ti and DE between ti versus the reference time t0
for the biological condition BC1.
- the following plots from the temporal statistical analysis
- a barplot which gives the number of DE genes between ti and the reference
time t0, for each time ti (except the reference time t0) and
biological condition.
- alluvial graphs of DE genes, one per biological condition.
- $N_{bc}$ graphs showing the number of DE genes as a function of time
for each temporal group, one per biological condition. By temporal group,
we mean the sets of genes which are first DE at the same time.
- $2\times N_{bc}$ UpSet plot showing the number of DE genes
belonging to each DE temporal pattern, for each biological condition.
By temporal pattern, we mean the set of times ti such that the gene is
DE between ti and the reference time t0 (see `DEplotVennBarplotTime()`).
- an alluvial graph for DE genes which are DE at least one time for
each group.
- the following plots from the statistical analysis by biological condition
- a barplot which gives the number of specific DE genes for each
biological condition and time (see `DEplotBarplotFacetGrid()`).
- $N_{bc}(N_{bc}-1)/2$ UpSet plot which give the number of genes
for each possible intersection (set of pairs of biological conditions),
one per time (see `DEplotVennBarplotGroup()`).
- an alluvial graph of genes which are specific at least one time
(see `DEplotAlluvial()`).
- the following plots from the combination of temporal and biological
statistical analysis
- a barplot which gives the number of signature genes for each biological
condition and time (see `DEplotBarplotFacetGrid()`).
- a barplot showing the number of genes which are DE at at least one time,
specific at at least one time and signature at at least one time,
for each biological condition.
- an alluvial graph of genes which are signature at least one time
(see `DEplotAlluvial()`).
* **Input**:
- A gene is considered as DE if the pvalue associated to the gene is inferior
to `pval.min` (numeric value between 0 and 1) and if the absolute
log fold change associated to the gene is superior to `log.FC.min`
(numeric positive value).
- In order to display the different output graph, set `Plot.DE.graph=TRUE`.
- In order to save the different results in a folder,
select a folder path in `path.result`.
- Set `RawCounts=RawCounts_Weger2021_MOUSEsub500` in order to use
the complete dataset.
* **Other**:
- Write `?DEanalysisGlobal` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
11. Volcano and ratio intensity (MA) plots with `DEplotVolcanoMA()`
```{r VolcanoMA_Mouse2, eval=TRUE, echo=TRUE}
resMAvolcano <- DEplotVolcanoMA(SEresDE=resDE,
NbGene.plotted=2,
SizeLabel=3,
Display.plots=FALSE,
Save.plots=FALSE)
```
* **Output**: The function returns Volcano plots and MA plots from the results
of our function `DEanalysisGlobal()`.
* **Input**:
- In order to display the different output graph, set `Display.plots=TRUE`.
- In order to save the different results in a folder, set `Save.plots=TRUE`
and and a folder path in the input `path.result`
in the function `DEanalysisGlobal()`.
* **Other**:
- Write `?DEplotVolcanoMA` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
12. Heatmaps with `DEplotHeatmaps()`
```{r Heatmaps_Mouse2, eval=TRUE, echo=TRUE}
resHEAT <- DEplotHeatmaps(SEresDE=resDE,
ColumnsCriteria=2,
Set.Operation="union",
NbGene.analysis=20,
SizeLabelRows=5,
SizeLabelCols=5,
Display.plots=FALSE,
Save.plots=TRUE)
```
* **Output**: The function returns two heatmaps: one heatmap of gene expressions
between samples and selected genes; and a correlation heatmap between samples.
* **Input**:
- If `Set.Operation="union"` then the rows extracted from the different
datasets included in `SEresDE` are those such that the sum of the selected
columns of `SummarizedExperiment::rowData(SEresDE)` by `ColumnsCriteria`
is >0. For example, the selected genes can be those DE at least at t1 or t2
(versus the reference time t0).
- In order to display the different output graph, set `Display.plots=TRUE`.
- In order to save the different results in a folder, set `Save.plots=TRUE`
and and a folder path in the input `path.result`
in the function `DEanalysisGlobal()`.
* **Other**:
- Write `?DEplotHeatmaps` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
13. GO enrichment analysis with `GSEAQuickAnalysis()` and `GSEApreprocessing()`
```{r GSEAquickAnalysis_Mouse2, eval=TRUE}
resRgprofiler2 <- GSEAQuickAnalysis(Internet.Connection=FALSE,
SEresDE=resDE,
ColumnsCriteria=2,
ColumnsLog2ordered=NULL,
Set.Operation="union",
Organism="mmusculus",
MaxNumberGO=20,
Background=FALSE,
Display.plots=FALSE,
Save.plots=TRUE)
```
* **Output**: The function realizes, from the outputs of `DEanalysisGlobal()`,
an enrichment analysis (GSEA) of a subset of genes with the R package
`gprofiler2`.
- A data.frame which contains the outputs of `gprofiler2::gost()`
- A Manhattan plot showing all GO names according to their pvalue.
- A lollipop graph showing the `MaxNumberGO` most important GO.
* **Input**:
- If `Set.Operation="union"` then the rows extracted from the different
datasets included in `SEresDE` are those such that the sum of the selected
columns of `SummarizedExperiment::rowData(SEresDE)` by `ColumnsCriteria`
is >0. For example, the selected genes can be those DE at least at t1 or t2
(versus the reference time t0).
- In order to display the different output graph, set `Display.plots=TRUE`.
- In order to save the different results in a folder, set `Save.plots=TRUE`
and a folder path in the input `path.result`
in the function `DEanalysisGlobal()`.
* **Other**:
- Write `?GSEAQuickAnalysis` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
```{r GSEAprepro_Mouse2, eval=TRUE}
resGSEApreprocess <- GSEApreprocessing(SEresDE=resDE,
ColumnsCriteria=2,
Set.Operation="union",
Rnk.files=FALSE,
Save.files=FALSE)
```
* **Output**:
- A vector of character containing gene names specified by `ColumnsCriteria`
and `Set.Operation`.
- A vector of character containing all gene names
- And, in case where `Save.files=TRUE` and the path.result of
`DEanalysisGlobal()` is not NULL, specific files designed to be used as input
for the following online tools and software : GSEA, DAVID, WebGestalt,
gProfiler, Panther, ShinyGO, Enrichr, GOrilla
* **Input**: If `Set.Operation="union"` then the rows extracted from the
different datasets included in `SEresDE` are those such that the sum of
the selected columns of `SummarizedExperiment::rowData(SEresDE)`
by `ColumnsCriteria` is >0.
For example, the selected genes can be those DE at least at t1 or t2
(versus the reference time t0).
* **Other**:
- Write `?GSEApreprocessing` in your console for more information
about the function.
- For a more complete description of the function and package,
the user must read our pdf file entitled 'MultiRNAflow_vignette-knitr.pdf".
# SessionInfo
Here is the output of `sessionInfo()` on the system on which this document
was compiled.
```{r SessionInfo}
sessionInfo()
```