DESeq2 exampleDESeq2 is one of the most commonly used packages for performing differential expression analysis. It has been for a while on the top 5% of downloaded Bioconductor packages. The vignette explains how to run a DESeq2 analysis and describes in further detail the example data that we are going to use here.
This example will show you how to use regionReport to make interactive HTML reports from DESeq2 results. It will also cover how to create a PDF report.
First we need to load the example data from the pasilla package. Once we have loaded the data, we need to create a DESeqDataSet object and then perform the differential expression analysis using the DESeq() function.
## Load example data from pasilla package
library('pasilla')
library('DESeq2')## Loading required package: S4Vectors## Loading required package: methods## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colnames,
## do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, lengths, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff,
## sort, table, tapply, union, unique, unsplit##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:base':
##
## colMeans, colSums, expand.grid, rowMeans, rowSums## Loading required package: IRanges## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Get the necessary data from the pasilla package
data('pasillaGenes')
countData <- counts(pasillaGenes)
colData <- pData(pasillaGenes)[, c('condition', 'type')]
## Create DESeqDataSet object from the pasilla package
dds <- DESeqDataSetFromMatrix(countData = countData,
colData = colData,
design = ~ condition)
## Perform the differential expression analysis
dds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingOnce we have the input data we can use the DESeq2Report() function from regionReport. First we create a directory where we’ll have the data.
## The output will be saved in the 'DESeq2-example' directory
dir.create('DESeq2-example', showWarnings = FALSE, recursive = TRUE)Next, we can create the HTML report. In this case we’ll change the default theme by using the theme_bw() function from the ggplot2 package.
## Use ggplot2::theme_bw()
library('ggplot2')
library('regionReport')
## Create the HTML report
report <- DESeq2Report(dds, project = 'DESeq2 HTML report',
intgroup = c('condition', 'type'), outdir = 'DESeq2-example',
output = 'index', theme = theme_bw())## Warning in c.bibentry(knitcitations = citation("knitcitations"),
## regionReport = citation("regionReport")[1], : method is only applicable to
## 'bibentry' objects## Writing 9 Bibtex entries ...## OK
## Results written to file 'DESeq2-example/index.bib'##
##
## processing file: index.Rmd##
|
| | 0%
|
|. | 2%
## inline R code fragments
##
##
|
|... | 4%
## label: docSetup (with options)
## List of 3
## $ bootstrap.show.code : logi FALSE
## $ dev : symbol device
## $ bootstrap.show.message: logi FALSE
##
##
|
|.... | 6%
## ordinary text without R code
##
##
|
|...... | 9%
## label: setup (with options)
## List of 1
## $ bootstrap.show.message: logi FALSE
##
##
|
|....... | 11%
## ordinary text without R code
##
##
|
|........ | 13%
## label: PCA##
|
|.......... | 15%
## inline R code fragments
##
##
|
|........... | 17%
## label: sampleDist##
|
|............ | 19%
## ordinary text without R code
##
##
|
|.............. | 21%
## label: MAplotalpha##
|
|............... | 23%
## inline R code fragments
##
##
|
|................. | 26%
## label: MAplotalphaHalf##
|
|.................. | 28%
## inline R code fragments
##
##
|
|................... | 30%
## label: MAplotalpha-nBest##
|
|..................... | 32%
## inline R code fragments
##
##
|
|...................... | 34%
## label: pvalueHistogram##
|
|........................ | 36%
## ordinary text without R code
##
##
|
|......................... | 38%
## label: pvalueSumm
##
|
|.......................... | 40%
## ordinary text without R code
##
##
|
|............................ | 43%
## label: pvalueTable (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|............................. | 45%
## ordinary text without R code
##
##
|
|.............................. | 47%
## label: padjHistogram##
|
|................................ | 49%
## inline R code fragments
##
##
|
|................................. | 51%
## label: padjSumm
##
|
|................................... | 53%
## inline R code fragments
##
##
|
|.................................... | 55%
## label: padjTable (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|..................................... | 57%
## inline R code fragments
##
##
|
|....................................... | 60%
## label: topFeatures (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|........................................ | 62%
## inline R code fragments
##
##
|
|......................................... | 64%
## label: plotCounts##
|
|........................................... | 66%
## inline R code fragments
##
##
|
|............................................ | 68%
## label: edgeR-BCV (with options)
## List of 2
## $ eval: symbol isEdgeR
## $ echo: language isEdgeR & outputIsHTML
##
##
|
|.............................................. | 70%
## inline R code fragments
##
##
|
|............................................... | 72%
## label: edgeR-MDS (with options)
## List of 2
## $ eval: symbol isEdgeR
## $ echo: language isEdgeR & outputIsHTML
##
##
|
|................................................ | 74%
## inline R code fragments
##
##
|
|.................................................. | 77%
## label: unnamed-chunk-1 (with options)
## List of 2
## $ child: symbol customCode
## $ eval : symbol hasCustomCode
##
##
|
|................................................... | 79%
## inline R code fragments
##
##
|
|..................................................... | 81%
## label: thecall (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...................................................... | 83%
## ordinary text without R code
##
##
|
|....................................................... | 85%
## label: reproducibility1 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|......................................................... | 87%
## ordinary text without R code
##
##
|
|.......................................................... | 89%
## label: reproducibility2 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|........................................................... | 91%
## ordinary text without R code
##
##
|
|............................................................. | 94%
## label: reproducibility3 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............................................................. | 96%
## inline R code fragments
##
##
|
|................................................................ | 98%
## label: bibliography (with options)
## List of 3
## $ results: chr "asis"
## $ echo : logi FALSE
## $ warning: logi FALSE
##
##
|
|.................................................................| 100%
## ordinary text without R code## output file: index.knit.md## /usr/local/bin/pandoc +RTS -K512m -RTS index.utf8.md --to html --from markdown+autolink_bare_uris+ascii_identifiers+tex_math_single_backslash-implicit_figures --output index.html --smart --email-obfuscation none --self-contained --standalone --section-divs --table-of-contents --toc-depth 3 --variable toc_float=1 --variable toc_selectors=h1,h2,h3 --variable toc_collapsed=1 --variable toc_smooth_scroll=1 --template /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/default.html --variable 'theme:spacelab' --include-in-header /var/folders/cx/n9s558kx6fb7jf5z_pgszgb80000gn/T//RtmpuTKlSb/rmarkdown-str399f2971fd33.html --mathjax --variable 'mathjax-url:https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --no-highlight --variable highlightjs=/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/highlight --variable navigationjs=/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/navigation-1.0 --variable code_folding=hide##
## Output created: index.htmlYou can view the final HTML report at DESeq2-example.
The HTML report has an interactive table with the top features (in this case genes). It allows you to re-order these top features by different criteria or even search for your feature of interest. However, sometimes you might prefer to create a PDF report. The following code will allow you to create such a report.
## Create PDF version of the same report
report <- DESeq2Report(dds, project = 'DESeq2 PDF report',
intgroup = c('condition', 'type'), outdir = 'DESeq2-example',
output = 'DESeq2Report', theme = theme_bw(), output_format = 'pdf_document',
device = 'pdf')## '/Library/Frameworks/R.framework/Resources/bin/R' --no-site-file \
## --no-environ --no-save --no-restore --quiet \
## --file='DESeq2-report-isolated-PDF.R' --slave## You can view the final PDF report at DESeq2-example PDF.
## Date generated:
Sys.time()## [1] "2016-04-12 07:32:16 EDT"## Time spent making this page:
proc.time()## user system elapsed
## 87.654 4.009 93.208## R and packages info:
options(width = 120)
devtools::session_info()## Session info -----------------------------------------------------------------------------------------------------------## setting value
## version R version 3.3.0 alpha (2016-03-23 r70368)
## system x86_64, darwin13.4.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## tz America/New_York
## date 2016-04-12## Packages ---------------------------------------------------------------------------------------------------------------## package * version date source
## acepack 1.3-3.3 2014-11-24 CRAN (R 3.3.0)
## annotate 1.49.1 2016-02-06 Bioconductor
## AnnotationDbi 1.33.8 2016-04-10 Bioconductor
## backports 1.0.2 2016-03-18 CRAN (R 3.3.0)
## bibtex 0.4.0 2014-12-31 CRAN (R 3.3.0)
## Biobase * 2.31.3 2016-01-14 Bioconductor
## BiocGenerics * 0.17.4 2016-04-07 Bioconductor
## BiocParallel 1.5.21 2016-03-23 Bioconductor
## biomaRt 2.27.2 2016-01-14 Bioconductor
## Biostrings 2.39.12 2016-02-21 Bioconductor
## bitops 1.0-6 2013-08-17 CRAN (R 3.3.0)
## BSgenome 1.39.4 2016-02-21 Bioconductor
## bumphunter 1.11.5 2016-03-29 Bioconductor
## checkmate 1.7.4 2016-04-08 CRAN (R 3.3.0)
## cluster 2.0.3 2015-07-21 CRAN (R 3.3.0)
## codetools 0.2-14 2015-07-15 CRAN (R 3.3.0)
## colorout * 1.1-2 2016-03-24 Github (jalvesaq/colorout@f96c00c)
## colorspace 1.2-6 2015-03-11 CRAN (R 3.3.0)
## DBI 0.3.1 2014-09-24 CRAN (R 3.3.0)
## DEFormats 0.99.8 2016-03-31 Bioconductor
## derfinder 1.5.30 2016-03-25 Bioconductor
## derfinderHelper 1.5.3 2016-03-23 Bioconductor
## DESeq 1.23.1 2016-01-14 Bioconductor
## DESeq2 * 1.11.42 2016-04-10 Bioconductor
## devtools * 1.10.0 2016-01-23 CRAN (R 3.3.0)
## digest 0.6.9 2016-01-08 CRAN (R 3.3.0)
## doRNG 1.6 2014-03-07 CRAN (R 3.3.0)
## DT * 0.1 2015-06-09 CRAN (R 3.3.0)
## edgeR 3.13.8 2016-04-08 Bioconductor
## evaluate 0.8.3 2016-03-05 CRAN (R 3.3.0)
## foreach 1.4.3 2015-10-13 CRAN (R 3.3.0)
## foreign 0.8-66 2015-08-19 CRAN (R 3.3.0)
## formatR 1.3 2016-03-05 CRAN (R 3.3.0)
## Formula 1.2-1 2015-04-07 CRAN (R 3.3.0)
## genefilter 1.53.3 2016-03-23 Bioconductor
## geneplotter 1.49.0 2016-01-14 Bioconductor
## GenomeInfoDb * 1.7.6 2016-01-29 Bioconductor
## GenomicAlignments 1.7.20 2016-02-25 Bioconductor
## GenomicFeatures 1.23.29 2016-04-05 Bioconductor
## GenomicFiles 1.7.9 2016-02-22 Bioconductor
## GenomicRanges * 1.23.25 2016-03-31 Bioconductor
## ggplot2 * 2.1.0 2016-03-01 CRAN (R 3.3.0)
## gridExtra 2.2.1 2016-02-29 CRAN (R 3.3.0)
## gtable 0.2.0 2016-02-26 CRAN (R 3.3.0)
## highr 0.5.1 2015-09-18 CRAN (R 3.3.0)
## Hmisc 3.17-3 2016-04-03 CRAN (R 3.3.0)
## htmltools 0.3.5 2016-03-21 CRAN (R 3.3.0)
## htmlwidgets 0.6 2016-02-25 CRAN (R 3.3.0)
## httr 1.1.0 2016-01-28 CRAN (R 3.3.0)
## IRanges * 2.5.43 2016-04-10 Bioconductor
## iterators 1.0.8 2015-10-13 CRAN (R 3.3.0)
## jsonlite 0.9.19 2015-11-28 CRAN (R 3.3.0)
## knitcitations 1.0.7 2015-10-28 CRAN (R 3.3.0)
## knitr * 1.12.3 2016-01-22 CRAN (R 3.3.0)
## knitrBootstrap 1.0.0 2016-03-24 Github (jimhester/knitrBootstrap@cdaa4a9)
## labeling 0.3 2014-08-23 CRAN (R 3.3.0)
## lattice 0.20-33 2015-07-14 CRAN (R 3.3.0)
## latticeExtra 0.6-28 2016-02-09 CRAN (R 3.3.0)
## limma 3.27.14 2016-03-23 Bioconductor
## locfit 1.5-9.1 2013-04-20 CRAN (R 3.3.0)
## lubridate 1.5.6 2016-04-06 CRAN (R 3.3.0)
## magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
## markdown 0.7.7 2015-04-22 CRAN (R 3.3.0)
## Matrix 1.2-4 2016-03-02 CRAN (R 3.3.0)
## matrixStats 0.50.1 2015-12-15 CRAN (R 3.3.0)
## memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
## munsell 0.4.3 2016-02-13 CRAN (R 3.3.0)
## nnet 7.3-12 2016-02-02 CRAN (R 3.3.0)
## pasilla * 0.11.0 2016-03-31 Bioconductor
## pheatmap * 1.0.8 2015-12-11 CRAN (R 3.3.0)
## pkgmaker 0.22 2014-05-14 CRAN (R 3.3.0)
## plyr 1.8.3 2015-06-12 CRAN (R 3.3.0)
## qvalue 2.3.2 2016-01-14 Bioconductor
## R6 2.1.2 2016-01-26 CRAN (R 3.3.0)
## RColorBrewer * 1.1-2 2014-12-07 CRAN (R 3.3.0)
## Rcpp 0.12.4 2016-03-26 CRAN (R 3.3.0)
## RCurl 1.95-4.8 2016-03-01 CRAN (R 3.3.0)
## RefManageR 0.10.13 2016-04-04 CRAN (R 3.3.0)
## regionReport * 1.5.47 2016-04-12 Bioconductor
## registry 0.3 2015-07-08 CRAN (R 3.3.0)
## reshape2 1.4.1 2014-12-06 CRAN (R 3.3.0)
## RJSONIO 1.3-0 2014-07-28 CRAN (R 3.3.0)
## rmarkdown 0.9.5 2016-02-22 CRAN (R 3.3.0)
## rngtools 1.2.4 2014-03-06 CRAN (R 3.3.0)
## rpart 4.1-10 2015-06-29 CRAN (R 3.3.0)
## Rsamtools 1.23.8 2016-04-10 Bioconductor
## RSQLite 1.0.0 2014-10-25 CRAN (R 3.3.0)
## rstudioapi 0.5 2016-01-24 CRAN (R 3.3.0)
## rtracklayer 1.31.10 2016-04-07 Bioconductor
## S4Vectors * 0.9.46 2016-04-07 Bioconductor
## scales 0.4.0 2016-02-26 CRAN (R 3.3.0)
## stringi 1.0-1 2015-10-22 CRAN (R 3.3.0)
## stringr 1.0.0 2015-04-30 CRAN (R 3.3.0)
## SummarizedExperiment * 1.1.23 2016-04-06 Bioconductor
## survival 2.38-3 2015-07-02 CRAN (R 3.3.0)
## VariantAnnotation 1.17.23 2016-04-07 Bioconductor
## withr 1.0.1 2016-02-04 CRAN (R 3.3.0)
## XML 3.98-1.4 2016-03-01 CRAN (R 3.3.0)
## xtable 1.8-2 2016-02-05 CRAN (R 3.3.0)
## XVector 0.11.8 2016-04-06 Bioconductor
## yaml 2.1.13 2014-06-12 CRAN (R 3.3.0)
## zlibbioc 1.17.1 2016-03-19 Bioconductor