---
title: "Lab 1.3: Starting an analysis"
output:
BiocStyle::html_document:
toc: true
vignette: >
% \VignetteIndexEntry{Lab 1.3: Starting an analysis}
% \VignetteEngine{knitr::rmarkdown}
---
```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```
```{r setup, echo=FALSE, messages=FALSE, warnings=FALSE}
knitr::opts_chunk$set(
eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE"))
)
suppressPackageStartupMessages({
library(Biostrings)
library(GenomicRanges)
})
```
Original Authors: Martin Morgan<br />
Presenting Authors: [Martin Morgan][], [Lori Shepherd][]</br >
Date: 22 July, 2019</br >
Back: [Monday labs](lab-1-intro-to-r-bioc.html)
[Martin Morgan]: mailto: Martin.Morgan@RoswellPark.org
[Lori Shepherd]: mailto: Lori.Shepherd@RoswellPark.org
**Objective**: An overview of software available in _Bioconductor_.
**Lessons learned**:
- How to discover _CRAN_ and _Bioconductor_ packages
- Finding vignettes, workflows, and tutorials
- Exploring an existing analysis
- Approaching your own data
- How to get help
# _Bioconductor_
Analysis and comprehension of high-throughput genomic data
- Statistical analysis: large data, technological artifacts, designed
experiments; rigorous
- Comprehension: biological context, visualization, reproducibility
- High-throughput
- Sequencing: RNASeq, ChIPSeq, variants, copy number, ...
- Microarrays: expression, SNP, ...
- Flow cytometry, proteomics, images, ...
## Packages, vignettes, work flows
- 1750 software packages
- Discover and navigate via [biocViews][]
- Package 'landing page', e.g., `r Biocpkg("Gviz")`
- Title, author / maintainer, short description, citation,
installation instructions, ..., download statistics
- All user-visible functions have help pages, most with runnable
examples
- 'Vignettes' an important feature in _Bioconductor_ -- narrative
documents illustrating how to use the package, with integrated code
- Example: `AnnotationHub` [landing page][AH] references
[HOW-TO vignette][AH-howto] illustrating some fun use cases.
- Some are extensive; check out [Gviz][], [limma][], [edgeR][],
[DESeq2][]!
- 'Release' (every six months) and 'devel' branches
[AH]: https://bioconductor.org/packages/AnnotationHub
[AH-howto]: http://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/AnnotationHub-HOWTO.html
# Packages
## Finding packages
'Domain-specific' analysis
**Exercise**
- Visit the listing of [All Packages][].
- Use the 'Search biocViews' box in the upper left to identify
packages that have been tagged for RNASeq analysis. Explore other
analysis like ChIPSeq, Epigenetics, VariantAnnotation, Proteomics,
SingleCell, etc. Explore the graph of software packages by
expanding and contracting individual terms.
- Return to RNASeq. Two very popular packages are DESeq2 and
edgeR. Visit the 'landing page' of one of these packages. The
landing page has a title, authors, instructions for citing use the
package, etc.
- Breifly explore the vignette and reference manual links. When
would you consult the vignette? When would the reference manual be
helpful?
[All Packages]: https://bioconductor.org/packages/release/BiocViews.html#___Software
_Bioconductor_ provides 'infrastructure' for working with genomic
data. We'll explore some of these in more detail in a later part of
this lab. For now...
**Exercise**
- Visit the landing pages of each of the [Biostrings][],
[GenomicRanges][], [VariantAnnotation][], and [GenomicAlignments][]
packages. Create a short summary describing what each package does,
and when it might be useful.
- Visit the landing page for the [SummarizedExperiment][]
package. This package is meant to provide a data representation that
helps manage, in a coordinated fashion, 'assays' (e.g., gene x
sample matrix of RNASeq counts) and the row (e.g., genomic
coordinates, P-values) and column (e.g., sample sheets). Briefly
review the user-oriented vignette ('SummarizedExperiment for
Coordinating Experimental Assays, Samples, and Regions of Interest')
to get a sense for how this package might be used.
- Visit the landing page for the [rtracklayer][] package. From the
reference manual, when would the various `import()` and `export()`
functions be useful?
Annotation packages are data-, rather than software-, centric,
providing information about the relationship between different
identifiers, gene models, reference genomes, etc.
**Exercise**
- On the page listing [All Packages][], click on the AnnotationData
top-level term.
- Search, using the box on the right-hand side, for annotation
packages that start with the following letters to get a sense of the
packages and organisms available.
- `org.*`: symbol mapping
- `TxDb.*` and `EnsDb.*`: gene models
- `BSgenome.*`: reference genomes
We'll see in a subsequent lab that a wealth of additional annotation
resources, including updated EnsDb and reference genomes, are
available through [AnnotationHub][].
Workflow packages are meant to provide a comprehensive introduction to
work flows that require several different packages. These can be quite
extensive documents, providing a very rich source of information.
**Exercise**
- Briefly explore the 'Simple Single Cell' workflow (or other workflow
relevant to your domain of interest) to get a sense of the material
the workflow covers.
## Installing packages
Likely the packages needed for this course are already
installed. Nonetheless it is useful to know how to install other
packages.
_Bioconductor_ has a particular approach to making packages
available. We have a 'devel' branch where new packages and features
are introduced, and a 'release' branch where users have access to
stable packages. Each six months, in spring and fall, the current
'devel' version of packages is branched to become the next
'release'. Packages within a release are tested with one another, so
it is important to install packages from the same release. The
[BiocManager][] package tries to make it easy to do this.
The first step to package installation is to make sure that the
[BiocManager][] package has been installed using standard _R_
procedures.
```{r, eval = FALSE}
if (!require(BiocManager))
install.packages("BiocManager", repos = "https://cran.r-project.org")
```
Then, install the package(s) you would like to use
```{r, eval = FALSE}
BiocManager::install(c("Biostrings", "GenomicRanges"))
```
[BiocManager][] knows how to install CRAN and github packages, too.
There are several common problems encountered with package
installation. Often, packages have been installed using methods
different from the one recommended here, and the packages are from
different _Bioconductor_ releases. This leads to problems when
packages from different releases are incompatible with one another.
**Exercise** Verify that your packages are current and installed from
the same _Bioconductor_ release with
```{r, eval = FALSE}
BiocManager::valid()
```
Two common problems are that some packages are too old (a newer
version of the package exists) or too new (some packages have been
installed using a method other than [BiocManager][]). If there are
packages that are too old or too new, it is almost always a good idea
to follow the instructions from `BiocManager::valid()` to correct the
situation.
## Loading and using packages
Packages need to be installed only once for each version of _R_ you
use, but need to be _loaded_ into each new _R_ session that you
start. Packages are loaded using
```{r, messages = FALSE}
library(Biostrings)
```
Whan a package is loaded, it can sometimes generate messages that are
informational only, if you are confident this is the case for the
packages you're loading, use `suppressPackageStartupMessages()` for a
quieter experience:
```{r}
suppressPackageStartupMessages({
library(GenomicRanges)
library(GenomicAlignments)
})
```
**Exercise** It is usually very helpful to explore package vignettes.
- Visit the vignette of the [DESeq2][] package, and walk through a few
steps to understand what the vignette provides in terms of
instructions for starting with the package, functionality the
package provides, mathematical and statistical details of the
implementation, and how the analysis provided by the package might
be extended by other packages in the _Bioconductor_ ecosystem. One
can visit vignettes through RStudio, or by running commands such as
```{r, eval = FALSE}
vignette(package = "DESeq2")
browseVignettes("DESeq2")
```
- Most vignettes are written in such a way that the _R_ code of the
vignette _must_ be correct for the vignette to be produced. The code
itself is available in the package. Find the code for the DESeq2
vignette
```{r}
dir(system.file(package="DESeq2", "doc"))
vign <- system.file(package="DESeq2", "doc", "DESeq2.R")
```
open it in RStudio (e.g., using File -> Open File... menu), step
through the first few lines of _R_ code and compare your output to
the output in the vignette. Alternatively, run the entire analysis
in the vignette with the command
```{r, eval = FALSE}
source(vign, echo = TRUE, max.lines = Inf)
```
**Exercise** Help pages provide more focused instructions for use of
particular functions. It is often con
- Load the [Biostrings][] package
```{r, message = FALSE}
library(Biostrings)
```
- Look for help on the function `letterFrequency()` using the command
```
?letterFrequency
```
note that there is tab completion after the `?` and first few
letters of the command.
- The help page is quite complicated, documenting several different
functions. In the 'Description' section, find a description of what
`letterFrequency()` does. In the 'Usage' section, find the arguments
that can be used with `letterFrequency()`, and try to understand,
from the `Arguments` section what each argument might be or how it
influences the computation. The `Value` section attempts to describe
the return value of the `letterFrequency()` function.
- Sometimes an example is worth a thousand words. Can you run the
first two sections of the example at the end of the help page (for
`alphabetFrequency()` and `letterFrequency()` to arrive at a better
understanding of how the `letterFrequency()` function works?
# Getting help
Where to get help?
- Most questions from users: https://support.bioconductor.org
What can you get help on?
- Problems on how to use software
- Statistical interpretation of results
- But not: extensive statistical consultation on, e.g., advanced
experimental design. Best source: local experts and collaborators.
How to ask a good question
- Simplify to just a few lines of _R_ code.
- Must be able to be run by someone else
- Use built-in data sets or simple simulations
- Don't reference files on your system!
- include the output of `sessionInfo()`, which often shows problems
with out-of-date packages.
**Exercise** Visit the support site and review the five most recent
questions. Which do you think are 'good', from the guidelines offered
above? Which have received helpful answers? Can you figure out who the
person answering the question is, i.e., why do they think they have an
answer?
# A sequence analysis package tour
This very open-ended topic points to some of the most prominent
_Bioconductor_ packages for sequence analysis. Use the opportunity in
this lab to explore the package vignettes and help pages highlighted
below; many of the material will be covered in greater detail in
subsequent labs and lectures.
Basics
- _Bioconductor_ packages are listed on the [biocViews][] page. Each
package has 'biocViews' (tags from a controlled vocabulary)
associated with it; these can be searched to identify appropriately
tagged packages, as can the package title and author.
- Each package has a 'landing page', e.g., for
[GenomicRanges][]. Visit this landing page, and note the
description, authors, and installation instructions. Packages are
often written up in the scientific literature, and if available the
corresponding citation is present on the landing page. Also on the
landing page are links to the vignettes and reference manual and, at
the bottom, an indication of cross-platform availability and
download statistics.
- A package needs to be installed once, using the instructions on the
landing page. Once installed, the package can be loaded into an R
session and the help system queried interactively, as outlined
above:
```{r require}
library(GenomicRanges)
```
```{r help, eval=FALSE}
help(package="GenomicRanges")
vignette(package="GenomicRanges")
vignette(package="GenomicRanges", "GenomicRangesHOWTOs")
?GRanges
```
Domain-specific analysis -- explore the landing pages, vignettes, and
reference manuals of two or three of the following packages.
- Important packages for analysis of differential expression include
[edgeR][] and [DESeq2][]; both have excellent vignettes for
exploration. Additional research methods embodied in _Bioconductor_
packages can be discovered by visiting the [biocViews][] web page,
searching for the 'DifferentialExpression' view term, and narrowing
the selection by searching for 'RNA seq' and similar.
- Single-cell 'omics are increasingly important. From the
[biocViews][] page, enter 'single cell' in the 'search table'
field. Check out the [single-cell workflow][sc-workflow] for an
overview of key packages.
- Popular ChIP-seq packages include [DiffBind][] and [csaw][] for
comparison of peaks across samples, [ChIPQC][] for quality
assessment, and [ChIPpeakAnno][] and [ChIPseeker][] for annotating
results (e.g., discovering nearby genes). What other ChIP-seq
packages are listed on the [biocViews][] page?
- Working with called variants (VCF files) is facilitated by packages
such as [VariantAnnotation][], [VariantFiltering][], and
[ensemblVEP][]; packages for calling variants include, e.g.,
[h5vc][] and [VariantTools][].
- Several packages identify copy number variants from sequence data,
including [cn.mops][]; from the [biocViews][] page, what other copy
number packages are available? The [CNTools][] package provides some
useful facilities for comparison of segments across samples.
- Microbiome and metagenomic analysis is facilitated by packages such
as [phyloseq][] and [metagenomeSeq][].
- Metabolomics, chemoinformatics, image analysis, and many other
high-throughput analysis domains are also represented in
_Bioconductor_; explore these via biocViews and title searches.
Working with sequences, alignments, common web file formats, and raw
data; these packages rely very heavily on the [IRanges][] /
[GenomicRanges][] infrastructure that we will encounter later in the
course.
- The [Biostrings][] package is used to represent DNA and other
sequences, with many convenient sequence-related functions. Check
out the functions documented on the help page `?consensusMatrix`,
for instance. Also check out the [BSgenome][] package for working
with whole genome sequences, e.g., `?"getSeq,BSgenome-method"`
- The [GenomicAlignments][] package is used to input reads aligned to
a reference genome. See for instance the `?readGAlignments` help
page and `vigentte(package="GenomicAlignments",
"summarizeOverlaps")`
- The [rtracklayer][] `import` and `export` functions can read in many
common file types, e.g., BED, WIG, GTF, ..., in addition to querying
and navigating the UCSC genome browser. Check out the `?import` page
for basic usage.
- The [ShortRead][] and [Rsamtools][] packages can be used for
lower-level access to FASTQ and BAM files, respectively.
- Many genomics data files are very large. We'll explore strategies of
_restriction_ (only input some of the data in the file) and
_iteration_ (read the file in chunks, rather than its entirety) for
processing large data in other labs.
Annotation: _Bioconductor_ provides extensive access to 'annotation'
resources (see the [AnnotationData][] biocViews hierarchy); these are
covered in greater detail in Thursday's lab, but some interesting
examples to explore during this lab include:
- [biomaRt][], [PSICQUIC][], [KEGGREST][] and other packages for
querying on-line resources; each of these have informative vignettes.
- [AnnotationDbi][] is a cornerstone of the
[Annotation Data][AnnotationData] packages provided by _Bioconductor_.
- **org** packages (e.g., [org.Hs.eg.db][]) contain maps between
different gene identifiers, e.g., ENTREZ and SYMBOL. The basic
interface to these packages is described on the help page `?select`
- **TxDb** packages (e.g., [TxDb.Hsapiens.UCSC.hg38.knownGene][])
contain gene models (exon coordinates, exon / transcript
relationships, etc) derived from common sources such as the hg38
knownGene track of the UCSC genome browser. These packages can be
queried, e.g., as described on the `?exonsBy` page to retrieve all
exons grouped by gene or transcript.
- **EnsDb** packages and databases (e.g. [EnsDb.Hsapiens.v86][]) provide,
similar to TxDb packages, gene models, but also protein annotations (protein
sequences and protein domains within these) and additional annotation
columns such as `"gene_biotype"` or `"tx_biotype"` defining the *biotype* of
the features (e.g. lincRNA, protein_coding, miRNA etc). EnsDb databases are
designed for [Ensembl][] annotations and contain annotations for all genes
(protein coding and non-coding) for a specific Ensembl release.
- **BSgenome** packages (e.g., [BSgenome.Hsapiens.UCSC.hg38][])
contain whole genomes of model organisms.
- [VariantAnnotation][] and [ensemblVEP][] provide access to sequence
annotation facilities, e.g., to identify coding variants; see the
[Introduction to VariantAnnotation](http://bioconductor.org/packages/release/bioc/vignettes/ShortRead/inst/doc/Overview.pdf)
vignette for a brief introduction; we'll re-visit this during the
Thursday lab.
- Take a quick look (there are more activites in other labs) at the
[annotation work flow](http://bioconductor.org/help/workflows/annotation/annotation/)
on the _Bioconductor_ web site.
A number of _Bioconductor_ packages help with visualization and
reporting, in addition to functions provided by indiidual packages.
- [Gviz][] provides a track-like visualization of genomic regions;
it's got an amazing vignette.
- [ComplexHeatmap][] does an amazing job of all sorts of heatmaps,
including OncoPrint-style summaries.
- [ReportingTools][] provides a flexible way to generate static and
dynamic HTML-based reports.
# End matter
## Session Info
```{r}
sessionInfo()
```
## Acknowledgements
Research reported in this tutorial was supported by the National Human
Genome Research Institute and the National Cancer Institute of the
National Institutes of Health under award numbers U41HG004059 and
U24CA180996.
This project has received funding from the European Research Council
(ERC) under the European Union's Horizon 2020 research and innovation
programme (grant agreement number 633974)
[biocViews]: http://bioconductor.org/packages/release/BiocViews.html#___Software
[AnnotationData]: http://bioconductor.org/packages/release/BiocViews.html#___AnnotationData
[aprof]: http://cran.r-project.org/web/packages/aprof/index.html
[hexbin]: http://cran.r-project.org/web/packages/hexbin/index.html
[lineprof]: https://github.com/hadley/lineprof
[microbenchmark]: http://cran.r-project.org/web/packages/microbenchmark/index.html
[AnnotationDbi]: http://bioconductor.org/packages/AnnotationDbi
[AnnotationHub]: https://bioconductor.org/packages/AnnotationHub
[BSgenome]: http://bioconductor.org/packages/BSgenome
[BiocManager]: https://cran.r-project.org/package=BiocManager
[Biostrings]: http://bioconductor.org/packages/Biostrings
[CNTools]: http://bioconductor.org/packages/CNTools
[ChIPQC]: http://bioconductor.org/packages/ChIPQC
[ChIPpeakAnno]: http://bioconductor.org/packages/ChIPpeakAnno
[ChIPseeker]: http://bioconductor.org/packages/ChIPseeker
[ComplexHeatmap]: http://bioconductor.org/packages/ComplexHeatmap
[csaw]: http://bioconductor.org/packages/csaw
[DESeq2]: http://bioconductor.org/packages/DESeq2
[DiffBind]: http://bioconductor.org/packages/DiffBind
[GenomicAlignments]: http://bioconductor.org/packages/GenomicAlignments
[GenomicRanges]: http://bioconductor.org/packages/GenomicRanges
[Gviz]: http://bioconductor.org/packages/Gviz
[IRanges]: http://bioconductor.org/packages/IRanges
[KEGGREST]: http://bioconductor.org/packages/KEGGREST
[PSICQUIC]: http://bioconductor.org/packages/PSICQUIC
[rtracklayer]: http://bioconductor.org/packages/rtracklayer
[Rsamtools]: http://bioconductor.org/packages/Rsamtools
[ReportingTools]: http://bioconductor.org/packages/ReportingTools
[ShortRead]: http://bioconductor.org/packages/ShortRead
[SummarizedExperiment]: http://bioconductor.org/packages/SummarizedExperiment
[VariantAnnotation]: http://bioconductor.org/packages/VariantAnnotation
[VariantFiltering]: http://bioconductor.org/packages/VariantFiltering
[VariantTools]: http://bioconductor.org/packages/VariantTools
[biomaRt]: http://bioconductor.org/packages/biomaRt
[cn.mops]: http://bioconductor.org/packages/cn.mops
[h5vc]: http://bioconductor.org/packages/h5vc
[edgeR]: http://bioconductor.org/packages/edgeR
[ensemblVEP]: http://bioconductor.org/packages/ensemblVEP
[limma]: http://bioconductor.org/packages/limma
[metagenomeSeq]: http://bioconductor.org/packages/metagenomeSeq
[phyloseq]: http://bioconductor.org/packages/phyloseq
[snpStats]: http://bioconductor.org/packages/snpStats
[org.Hs.eg.db]: http://bioconductor.org/packages/org.Hs.eg.db
[TxDb.Hsapiens.UCSC.hg38.knownGene]: http://bioconductor.org/packages/TxDb.Hsapiens.UCSC.hg38.knownGene
[BSgenome.Hsapiens.UCSC.hg38]: http://bioconductor.org/packages/BSgenome.Hsapiens.UCSC.hg38
[EnsDb.Hsapiens.v86]: http://bioconductor.org/packages/EnsDb.Hsapiens.v86
[Ensembl]: http://www.ensembl.org
[sc-workflow]: http://bioconductor.org/packages/release/workflows/html/simpleSingleCell.html