Last edited: 2019-02-06
set.seed(42) # knitr::opts_chunk$set(comment=NA, # fig.align="center", # fig.width = 7, # fig.height = 7, # warning=FALSE)
r Biocpkg("ideal") is an R package distributed as part of the Bioconductor project.
To install the package, start R and enter:
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install("ideal")
The GitHub repository for
r Biocpkg("ideal") is https://github.com/federicomarini/ideal.
This is the place to file an issue, report a bug, or provide a pull request.
r Rpackage("ideal") is installed, it can be loaded by the following command.
r Biocpkg("ideal") is a Bioconductor package containing a Shiny application for analyzing RNA-Seq data in the context of differential expression.
This enables an interactive and at the same time analysis, keeping the functionality accessible, and yet providing a comprehensive selection of graphs and tables to mine the dataset at hand.
r Biocpkg("ideal") is an R package which fully leverages the infrastructure of the Bioconductor project in order to deliver an interactive yet reproducible analysis for the detection of differentially expressed genes in RNA-Seq datasets.
Graphs, tables, and interactive HTML reports can be readily exported and shared across collaborators.
The dynamic user interface displays a broad level of content and information, subdivided by thematic tasks.
All in all, it aims to enforce a proper analysis, by reaching out both life scientists and experienced bioinformaticians, and also fosters the communication between the two sides, offering robust statistical methods and high standard of accessible documentation.
It is structured in a similar way to the
r Biocpkg("pcaExplorer"), also designed as an interactive companion tool for RNA-seq analysis focused rather on the exploratory data analysis e.g. using principal components analysis as a main tool.
The interactive/reactive design of the app, with a dynamically generated user interface makes it easy and immediate to apply the gold standard methods (in the current implementation, based on
r Biocpkg("DESeq2")) in a way that is information-rich and accessible also to the bench biologist, while also providing additional insight also for the experienced data analyst.
Reproducibility is supported via state saving and automated report generation.
If you use
r Biocpkg("ideal") for your analysis, please cite it as here below:
There are different ways to use
r Biocpkg("ideal") for interactive differential expression analysis.
First load the library:
and then launch the app with the
This takes the following essential parameters as input:
DESeqDataSetobject. If not provided, then a
expdesignneed to be provided. If none of the above is provided, it is possible to upload the data during the execution of the Shiny App
DESeqResultsobject. If not provided, it can be computed during the execution of the application
data.frameobject, with row.names as gene identifiers (e.g. ENSEMBL ids) and a column,
gene_name, containing e.g. HGNC-based gene symbols. If not provided, it can be constructed during the execution via the
countmatrix- a count matrix, with genes as rows and samples as columns. If not provided, it is possible to upload the data during the execution of the Shiny App
data.framecontaining the info on the experimental covariates of each sample. If not provided, it is possible to upload the data during the execution of the Shiny App
Different modalities are supported to launch the application:
ideal(dds_obj = dds, res_obj = res, annotation_obj = anno), where the objects are precomputed in the current session and provided as parameters
ideal(dds_obj = dds), as in the command above, but where the result object is assembled at runtime
ideal(countmatrix = countmatrix, expdesign = expdesign), where instead of passing the defined
DESeqDataSetobject, its components are given, namely the count matrix (e.g. generated after a run of featureCounts or HTSeq-count) and a data frame with the experimental covariates. The design formula can be constructed interactively at runtime
ideal(), where the count matrix and experimental design can simply be uploaded at runtime, where all the derived objects can be extracted and computed live. These files have to be formatted as tabular text files, and a function in the package tries to guess the separator, based on heuristics of occurrencies per line of commonly used characters
r Biocpkg("ideal") without installing any additional software, you can access the public instance of the Shiny Server made available at the Institute of Medical Biostatistics, Epidemiology and Informatics (IMBEI) in Mainz.
This resource is accessible at this address:
A deployment-oriented version of the package is available at https://github.com/federicomarini/ideal_serveredition.
This repository contains also detailed instruction to setup the running instance of a Shiny Server, where
r Biocpkg("ideal") can be run without further installation for the end-users.
Please note that you still need
r Biocpkg("ideal") to be installed there once during the setup phase - for this operation, you might require root administrator permissions.
The user interface is dynamically displayed according to the provided and computed objects, with tabs that are actively usable only once the required input is effectively available.
Moreover, for some relevant UI widgets, the user can receive additional information by hovering over with the mouse, with the functionality powered by the
r CRANpkg("shinyBS") package.
For the user which is either new with the app UI/functionality, or not extensively familiar with the topic of differential expression, it is possible to obtain a small guided tour of the App by clicking on the respective help buttons, marked in the app like this - please note that this button is clickable but does not start any tour.
These trigger the start of a step-by-step guide and feature introduction, powered by the
r CRANpkg("rintrojs") package.
Some of the input controls which affect different tabs are located in the sidebar, while others are as well in the individual tabs of the app. By changing one or more of the input parameters, the user can get a fine control on what is computed and displayed.
Width and Height for the figures to export are input here in cm.
This displays a list of the underlying objects with which basically all of the analysis can be performed.
A green tick icon appears close to each when the respective component is either provided or calculated.
For obtaining the best analysis experience in
r Biocpkg("ideal"), it is recommended to provide all of them.
Clicking on this button activated the
intro.js based tour for getting to know the components and the structure of the app.
Dedicated step-by-step procedures are also available in each individual tab.
The task menu, accessible by clicking on the cog icon in the upper right part of the application, provides two functionalities:
Exit ideal & savewill close the application and store the content of the
valuesreactive objects in a list of two elements in the
ideal_envenvironment, respectively called
Save State as .RDatawill similarly store
r_datain a binary file named
idealState_YYYYMMDD_HHMMSS.Rdata, without closing the application
r Biocpkg("ideal") app is a one-paged dashboard, structured in different panels, where each of them is focused on a different aspect of the data exploration.
On top of the panels, three
valueBox objects serve as guiding elements for having an overview of the data at hand: how many genes and samples are in the data, how many entries are in the annotation object, and how many genes were found to be differentially expressed in the results.
Whenever each of the underlying objects is available, the background color turns from red to green.
For the main analysis, the available panels are described in the following subsections.
The landing page for the app is also where you might likely be reading this text (otherwise in the package vignette).
The Data Setup panel is where you can upload or inspect the required inputs for running the app.
This builds on the primary idea used by
r Biocpkg("pcaExplorer") and extends it with the following aspects:
data.framein advance, and is based on the widely adopted
r Biocpkg("ideal")are), the information flow can move from the data exploration to decisions taken at the moment of testing
A diagnostic mean-dispersion plot is also provided in a collapsible element at the bottom of the panel, shown when the
DESeqDataSet is generated and the
DESeq command from the
DESeq2 package has been applied.
r Biocpkg("pcaExplorer"), interactive tables for the raw, normalized, or variance stabilized transformed (preferred over regularized logarithm transformation, because of its speed during the dispersion estimation - especially with larger number of samples) counts are shown in this tab.
The user can also generate a sample-to-sample correlation scatter plot with the selected data.
r Biocpkg("ideal") has an option to include a filter step at the gene level by removing genes with low absolute or averages low values.
After this, it might be possible to have to re-run the analysis in step 3 from the Data Setup panel.
This tab is an interface for generating the summary tables after testing for DE.
It is usually based on the Wald test, as implemented in DESeq2, but when the factor of interest is assuming more than two levels, the user can also perform an ANOVA-like test across the groups with the likelihood ratio test.
Options for enabling/disabling automated independent filtering, adding the additional column of unshrunken log2 fold change values (instead of the moderated estimates used by default), as well as using the Independent Hypothesis Weighting (
r Biocpkg("IHW")) framework, are provided.
The False Discovery Rate (FDR) can be set from the sidebar panel, and a couple of diagnostic plots, such as the histogram of raw p-values and the distribution of log2fc, are shown below the interactive enhanced version of the table - with clickable elements to link to ENSEMBL database and NCBI website.
In this tab an interactive MA plot for the contrast selected in the Extract Results tab is displayed. Clicking on a single gene in the zoomed plot (enabled by brushing in the main plot), it is possible to obtain a boxplot for its expression values, flanked by an overview of information accessed live from the Entrez database. Alternatively, a volcano plot of -log10(p-value) versus log fold change can provide a slightly different perspective. The subset of selected genes are also here presented in static and interactive heatmaps, with the underlying data accessible from the collapsible box element.
The functionality in the Gene Finder builds upon the one provided by
r Biocpkg("pcaExplorer"), and allows to query up to four genes in the same view, which can here be selected from a dropdown input list which supports autocompletion.
A combined summary table (with both normalized counts and results statistics) is located below an MA plot where the selected genes are marked and annotated on the plot. To avoid repeating this manually, the user can also quickly upload a list of genes as text file (one gene identifier per line), such as members of gene families (e.g. all cytokines, all immunoglobulines, ...) or defined by common function (e.g. all housekeeping genes, or others based on any annotation).
The Functional Analysis tab takes the user from the simple lists of DE genes to insight on the affected biological pathways, with three approaches based on the Gene Ontology (GO) databases. This panel of ideal has a slim interface to the following methods for performing functional analysis:
limma::goanafor the quick yet standard implementation
topGO, particularly valuable for pruning terms which are topologically less meaningful than their specific nodes
goseq, which accounts for the specific length bias intrinsic in RNA-Seq assays (longer genes have higher chances of being called DE).
r Biocpkg("ideal") allows the user to work simultaneously with more gene lists, two of which can be uploaded in a custom way (e.g. list of gene families, or extracted from other existing publications).
The interaction among these lists can be visually represented in Venn diagrams, as well as with the appealing alternative from the UpSetR package, where all combination of sets are explicitly shown.
Each of the methods for GO enrichment delivers its own interactive
DT-based table, which can then be explored interactively with the display of a heatmap for all the (DE) genes annotated to a particular term, picking the normalized transformed values for comparing robustly the expression values.
This is simply triggered by clicking any of the rows for the results tables.
Another useful feature is provided by the clickable link to the AmiGO database on each of the GO term identifiers.
The Signatures Explorer tab allows the user to check the behavior of a number of provided gene signatures in the data at hand, displaying this as a heatmap.
This panel is composed by different well panels:
in the Setup Options, you can select and upload a gene signature file, in
gmt format (e.g. like the ones provided in the MSigDB database, or from WikiPathways), and quickly compute the variance stabilized transformed version of your data, which is more amenable for visualization than raw or normalized counts
in the Conversion options tab, you can create an annotation vector, used to bring the ids from your data and the ids the
gmt used for encoding the signature elements.
This works based on the
the lower well panels control the appearance of the heatmap, also with an option to display all genes annotated in that pathway, or only the ones detected as differentially expressed (for this you need to provide or compute the result object)
The Report Editor tab works in the same way of
r Biocpkg("pcaExplorer"), with the scope of providing an interface to full computational reproducibility of the analyses.
Markdown options and
Editor options are available, and the text editor, based on the
shinyAce package, contains a comprehensive template report, that can be edited to the best convenience of the user.
The code contained in the template report fetches the latest state of the reactive values in the ongoing session, and its output is a comprehensive HTML file that can be expanded, edited, previewed in the tab itself, downloaded, and shared with a few mouse clicks.
r Biocpkg("pcaExplorer"), the functionality to display the report preview is based on
knit2html, and some elements such as
DataTable objects might not render correctly.
To render them correctly, please install the PhantomJS executable before launching the app.
This can be done by using the
r CRANpkg("webshot") package and calling
webshot::install_phantomjs() - HTML widgets will be rendered automatically as screenshots.
Keep in mind that the fully rendered report (the one you can obtain with the "Generate & Save" button) is not affected by this, since it uses
The About tab contains the output of
sessionInfo, plus general information on
r Biocpkg("ideal"), including the link to the Github development version.
If requested, the modular structure of the app can be easily expanded, and many new operations on the same set of input data and derived results can be embedded in the same framework.
r Biocpkg("ideal")on an exemplary data set
We can run
r Biocpkg("ideal") for demonstration purpose on published datasets that are available as SummarizedExperiment in an experiment Bioconductor packages.
We will use the
r Biocpkg("airway") dataset, which can be installed with this command:
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager") BiocManager::install("airway")
This package provides a
RangedSummarizedExperiment object of read counts in genes for an RNA-Seq experiment on four human airway smooth muscle cell lines treated with dexamethasone. More details such as gene models and count quantifications can be found in the
r Biocpkg("airway") package vignette.
r Biocpkg("ideal") on this dataset, the following commands are required. First, prepare the objects to be passed as parameters of
library(airway) library(DESeq2) data(airway) dds_airway <- DESeqDataSet(airway,design= ~ cell + dex) dds_airway # run deseq on it dds_airway <- DESeq(dds_airway) # extract the results res_airway <- results(dds_airway, contrast = c("dex","trt","untrt"),alpha = 0.05)
Then launch the app itself.
ideal(dds_obj = dds_airway) # or also providing the results object ideal(dds_obj = dds_airway,res_obj = res_airway)
annotation for this dataset can be built manually by exploiting the
r Biocpkg("org.Hs.eg.db") package.
library(org.Hs.eg.db) genenames_airway <- mapIds(org.Hs.eg.db,keys = rownames(dds_airway),column = "SYMBOL",keytype="ENSEMBL") annotation_airway <- data.frame(gene_id = rownames(dds_airway), gene_name = genenames_airway, row.names = rownames(dds_airway), stringsAsFactors = FALSE) head(annotation_airway)
or alternatively, can be handily created at runtime in the optional step.
Then again, the app can be launched with:
ideal(dds_obj = dds_airway, annotation_obj = annotation_airway)
If desired, alternatives can be used. See the well written annotation workflow available at the Bioconductor site (https://bioconductor.org/help/workflows/annotation/annotation/).
Let's suppose you performed part of your analysis with
r Biocpkg("edgeR") or
voom and you want to continue inspecting your data interactively using
You can use the functionality provided by the
r Biocpkg("DEFormats") package to convert the object with:
library(DEFormats) library(edgeR) library(limma) dge_airway <- as.DGEList(dds_airway) # this is your initial object # your factors for the design: dex <- colData(dds_airway)$dex cell <- colData(dds_airway)$cell redo_dds_airway <- as.DESeqDataSet(dge_airway) # force the design to ~cell + dex design(redo_dds_airway) <- ~cell+group #TODO: this is due to the not 100% re-conversion via DEFormats ### with edgeR y <- calcNormFactors(dge_airway) design <- model.matrix(~ cell + dex) y <- estimateDisp(y,design) # If you performed quasi-likelihood F-tests fit <- glmQLFit(y,design) qlf <- glmQLFTest(fit) # contrast for dexamethasone treatment topTags(qlf) # If you performed likelihood ratio tests fit <- glmFit(y,design) lrt <- glmLRT(fit) topTags(lrt) # lrt to DESeqResults tbledger <- lrt$table colnames(tbledger)[colnames(tbledger) == 'PValue'] <- 'pvalue' colnames(tbledger)[colnames(tbledger) == 'logFC'] <- 'log2FoldChange' colnames(tbledger)[colnames(tbledger) == 'logCPM'] <- 'baseMean' # get from the logcpm to something more the baseMean for better tbledger$baseMean <- (2^tbledger$baseMean) * mean(dge_airway$samples$lib.size) / 1e6 # use the constructor for DESeqResults edger_resu <- DESeqResults(DataFrame(tbledger)) edger_resu <- DESeq2:::pvalueAdjustment(edger_resu, independentFiltering = TRUE, alpha = 0.05, pAdjustMethod = "BH") # cfr with res_airway # plot(res_airway$pvalue,edger_resu$pvalue) ### with limma-voom x <- calcNormFactors(dge_airway) design <- model.matrix(~0+dex+cell) contr.matrix <- makeContrasts(dextrt - dexuntrt,levels=colnames(design)) v <- voom(x, design) vfit <- lmFit(v, design) vfit <- contrasts.fit(vfit, contrasts=contr.matrix) efit <- eBayes(vfit) tbllimma <- topTable(efit,number= Inf, sort.by = "none") colnames(tbllimma)[colnames(tbllimma) == 'P.Value'] <- 'pvalue' colnames(tbllimma)[colnames(tbllimma) == 'logFC'] <- 'log2FoldChange' colnames(tbllimma)[colnames(tbllimma) == 'AveExpr'] <- 'baseMean' colnames(tbllimma)[colnames(tbllimma) == 'adj.P.Val'] <- 'padj' # get from the logcpm to something more the baseMean for better tbllimma$baseMean <- (2^tbllimma$baseMean) * mean(dge_airway$samples$lib.size) / 1e6 # use the constructor for DESeqResults limma_resu <- DESeqResults(DataFrame(tbllimma)) # cfr with res_airway # plot(res_airway$pvalue,limma_resu$pvalue)
After preparing the objects, you can launch the app with the following command:
ideal(redo_dds_airway,edger_resu) # or ... ideal(redo_dds_airway,limma_resu)
Basically, you need a
myresults_df) with 3 columns to create the
DESeqResults object to be passed to the app, with this names:
baseMeanas an average expression value
log2FoldChangeas an indicator for the effect size
pvalueas a measure of significance
Then, you just need to call
padj can be computed either by
p.adjust or using
DESeq2:::pvalueAdjustment as in the example above.
The functions exported by the
r Biocpkg("ideal") package can be also used in a standalone scenario, provided the required objects are in the working environment.
They are listed here for an overview, but please refer to the documentation for additional details.
Where possible, for each function a code snippet will be provided for its typical usage.
deseqresult2tbl generate a tidy table with the results of DESeq2, sorted by the values in the
summary(res_airway) res_airway head(deseqresult2tbl(res_airway))
deseqresult2DEgenes only includes genes detected as DE.
head(deseqresult2DEgenes(res_airway,FDR = 0.05)) # the output in the first lines is the same, but dim(res_airway) dim(deseqresult2DEgenes(res_airway))
This tables can be enhanced with clickable links to the ENSEMBL and NCBI gene databases by the following commands:
myde <- deseqresult2DEgenes(res_airway,FDR = 0.05) myde$symbol <- mapIds(org.Hs.eg.db,keys = as.character(myde$id),column = "SYMBOL",keytype="ENSEMBL") myde_enhanced <- myde myde_enhanced$id <- ideal:::createLinkENS(myde_enhanced$id,species = "Homo_sapiens") myde_enhanced$symbol <- ideal:::createLinkGeneSymbol(myde_enhanced$symbol) DT::datatable(myde_enhanced[1:100,], escape = FALSE)
ggplotCounts extends the functionality of the
plotCounts function of
r Biocpkg("DESeq2"), and plots the normalized counts of a single gene as a boxplot, with jittered points superimposed.
ggplotCounts(dds = dds_airway, gene = "ENSG00000103196", # CRISPLD2 in the original publication intgroup = "dex")
annotation_obj is provided, their gene name can also be included in the title.
ggplotCounts(dds = dds_airway, gene = "ENSG00000103196", # CRISPLD2 in the original publication intgroup = "dex", annotation_obj = annotation_airway)
When used in the context of the app, it is possible to seamless search for genes also by their gene name, making exploration even more immediate.
goseqTable is a wrapper to extract the functional GO terms enriched in in a list of (DE) genes, based on the algorithm and the implementation in the
r Biocpkg("goseq") package.
Its counterpart, based on the
r Biocpkg("topGO") package, can be found in the
r Biocpkg("pcaExplorer") package.
Please note that the following code chunks are not run only for reducing compilation time. The reader is invited to try out the functionality by simply executing these commands.
res_subset <- deseqresult2DEgenes(res_airway)[1:100,] # taking only a subset of the DE genes myde <- res_subset$id myassayed <- rownames(res_airway) mygo <- goseqTable(de.genes = myde, assayed.genes = myassayed, genome = "hg38", id = "ensGene", testCats = "GO:BP", addGeneToTerms = FALSE) head(mygo)
As for the results, this table can be enhanced by adding the links for each category to the AmiGO database.
mygo_enhanced <- mygo # using the unexported function to link the GO term to the entry in the AmiGO db mygo_enhanced$category <- ideal:::createLinkGO(mygo_enhanced$category) DT::datatable(mygo_enhanced,escape=FALSE)
The MA plot provided by
r Biocpkg("ideal") displays the gene-wise log2-fold changes (logFCs) versus the average expression value.
As a main input parameter, a
DESeqResults object is required.
Control on the appearance of the plot can be applied by selecting the False Discovery Rate (
FDR), the point transparency (
point_alpha), adding horizontal lines at particular logFC values (
The user can also decide to add rug plots in the margins (setting
To facilitate the inspection of a particular gene or gene set,
intgenes can assume the value of a vector of genes (either the IDs or the gene symbols if
symbol column is provided in
Labels can be added via
labels_intgenes, while classical labels/title can be also edited as preferred (see
plot_ma for all details).
plot_ma(res_airway, FDR = 0.05, hlines = 1, title ="Adding horizontal lines", add_rug = FALSE) plot_ma(res_airway, FDR = 0.1, intgenes = c("ENSG00000103196", # CRISPLD2 "ENSG00000120129", # DUSP1 "ENSG00000163884", # KLF15 "ENSG00000179094"), # PER1 title = "Providing a shortlist of genes", add_rug = FALSE )
res_airway2 <- res_airway res_airway2$symbol <- mapIds(org.Hs.eg.db,keys = rownames(res_airway2),column = "SYMBOL",keytype="ENSEMBL") plot_ma(res_airway2, FDR = 0.05, intgenes = c("CRISPLD2", # CRISPLD2 "DUSP1", # DUSP1 "KLF15", # KLF15 "PER1"), # PER1 annotation_obj = annotation_airway, hlines = 2, add_rug = FALSE, title = "Putting gene names..." )
The volcano plot gives a different flavor for the gene overview, displaying log2-fold changes and log p-values.
In a way similar to
plot_ma, genes can be annotated with
intgenes, and vertical lines can be added via
ylim_up controls the y axis upper limit to visualize better the bulk of genes - keep in mind that very small p-values due to robust differences/large effect sizes can be "cut out".
plot_volcano(res_airway2, FDR = 0.05, intgenes = c("CRISPLD2", # CRISPLD2 "DUSP1", # DUSP1 "KLF15", # KLF15 "PER1"), # PER1 title = "Selecting the handful of genes - and displaying the full range for the -log10(pvalue)", ylim_up = -log10(min(res_airway2$pvalue, na.rm =TRUE)))
sepguesser makes an educated guess on the separator character for the input text file (
The separator list can be provided as a vector in
sep_list (defaults to comma, tab, semicolon, and whitespace - which ideally could cover most of the cases).
The heuristics is based on the number of occurrencies of each separator in each line.
sepguesser(system.file("extdata/design_commas.txt",package = "ideal")) sepguesser(system.file("extdata/design_semicolons.txt",package = "ideal")) sepguesser(system.file("extdata/design_spaces.txt",package = "ideal")) anyfile <- system.file("extdata/design_tabs.txt",package = "ideal") # we know it is going to be TAB guessed_sep <- sepguesser(anyfile) guessed_sep # to be used for reading in the same file, without having to specify the sep read.delim(anyfile, header = TRUE, as.is = TRUE, sep = guessed_sep, quote = "", row.names = 1, check.names = FALSE)
While running the app, the user can
reactiveValuesin an environment, or in binary format
This functionality to retrieve and share the output is provided by action buttons that are placed close to each element of interest.
annotation_obj is a quick helper to make your data and results easier to read. For creating it, you can exploit the corresponding
org.XX.eg.db packages, available in Bioconductor.
Currently available are the following (in alphabetical order):
They can all be easily installed with
In a similar way, for using at best the
goseq package, the gene length information can also be computed if the corresponding
TxDb packages are installed. Currently, following
TxDb packages are available:
As for the
org.XX.eg.db above, these can all be quickly installed with
Additional functionality for the
r Biocpkg("ideal") will be added in the future, as it is tightly related to a topic under current development research.
Improvements, suggestions, bugs, issues and feedback of any type can be sent to firstname.lastname@example.org.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.