%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{1. Introduction to Bioconductor}

# Bioconductor for Sequence Analysis -- Introduction

Epigenomics 2014<br>
Author: Martin Morgan (<a href="mailto:mtmorgan@fhcrc.org">mtmorgan@fhcrc.org</a>)<br>
Date: 24 August, 2014

```{r setup, echo=FALSE}
knitr::opts_chunk$set(cache=TRUE)
```

## Overall work flow

1. Experimental design
    - Keep it simple, e.g., 'control' and 'treatment' groups
    - Replicate within treatments!
2. Wet-lab sequence preparation
    - Record covariates, including processing day -- likely 'batch effects'
3. (Illumina) Sequencing (Bentley et al., 2008, doi:10.1038/nature07517)

    <img src="figures/Solexa-bridge-pcr.jpg" height="400px" width="300px" />

4. Alignment
    - Choose to match task, e.g., [Rsubread][], Bowtie2 good for ChIPseq,
      some forms of RNAseq; BWA, GMAP better for variant calling
5. Analysis
     a. Reduction, e.g., to 'count table'
     b. Differential expression, peak identification, ...
6. Comprehension

## 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

![Alt Sequencing Ecosystem](figures/SequencingEcosystem.png)

- 824 packages
- Discover and navigate via [biocViews][]
- Package 'landing page'
  - 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
- 'Release' (every six months) and 'devel' branches

Objects

- Represent complicated data types
- Foster interoperability
- S4 object system
  - Introspection: `getClass()`, `showMethods(..., where=search())`,
    `selectMethod()`
  - 'accessors' and other documented functions / methods for
    manipulation, rather than direct access to the object structure
- Interactive help
  - `method?"substr,<tab>"` to select help on methods, `class?D<tab>`
    for help on classes

Example

```{r Biostrings, message=FALSE}
require(Biostrings)                     # Biological sequences
data(phiX174Phage)                      # sample data, see ?phiX174Phage
phiX174Phage
m <- consensusMatrix(phiX174Phage)[1:4,] # nucl. x position counts
polymorphic <- which(colSums(m != 0) > 1)
m[, polymorphic]
```
```{r showMethods, eval=FALSE}
showMethods(class=class(phiX174Phage), where=search())
```

### Case study: Working with DNA sequence data 

1. Load the Biostrings package and phiX174Phage data set. What class
   is phiX174Phage? Find the help page for the class, and identify
   interesting functions that apply to it.
2. Discover vignettes in the Biostrings package with
   `vignette(package="Biostrings")`. Add another argument to the
   `vignette` function to view the 'BiostringsQuickOverview' vignette.
3. Navigate to the Biostrings landing page on
   http://bioconductor.org. Do this by visiting the biocViews
   page. Can you find the BiostringsQuickOverview vignette on the web
   site?
4. The following code loads some sample data, 6 versions of the
   phiX174Phage genome as a DNAStringSet object.
```{r phiX}
   library(Biostrings)
   data(phiX174Phage)
```
   Explain what the following code does, and how it works
```{r consensusMatrix}
   m <- consensusMatrix(phiX174Phage)[1:4,]
   polymorphic <- which(colSums(m != 0) > 1)
   mapply(substr, polymorphic, polymorphic, MoreArgs=list(x=phiX174Phage))
```

<a name="classes-and-methods"></a>

## 'S4' Classes, methods, and packages

This section focuses on classes, methods, and packages, with the goal
being to learn to navigate the help system and interactive discovery
facilities.

### Motivation

Sequence analysis is specialized
- Large data needs to be processed in a memory- and time-efficient manner
- Specific algorithms have been developed for the unique
  characteristics of sequence data

Additional considerations
- Re-use of existing, tested code is easier to do and less error-prone
  than re-inventing the wheel.
- Interoperability between packages is easier when the packages share
  similar data structures.

Solution: use well-defined _classes_ to represent complex data;
_methods_ operate on the classes to perform useful functions.  Classes
and methods are placed together and distributed as _packages_ so that
we can all benefit from the hard work and tested code of others.

<a name="i-granges"></a>

### Case study: IRanges and GRanges

The [IRanges][] package defines an important class for specifying
integer ranges, e.g.,
```{r iranges}
library(IRanges)
ir <- IRanges(start=c(10, 20, 30), width=5)
ir
```

There are many interesting operations to be performed on ranges, e.g,
`flank()` identifies adjacent ranges
```{r iranges-flank}
flank(ir, 3)
```

Consult the help page for flank, `?flank`, and explore other
range-based operations.

The `IRanges` class is part of a class hierarchy. To see this, ask R for
the class of `ir`, and for the class definition of the `IRanges` class
```{r iranges-class}
class(ir)
getClassDef(class(ir))
```

Notice that `IRanges` extends the `Ranges` class. Now try entering
`?"flank,<tab>`, where `<tab>` means to press the tab key to ask for
tab completion (may not be necessary in Rstudio). You can see that
there are help pages for several different classes. Tab-complete to

```{r iranges-flank-method, eval=FALSE}
?"flank,Ranges-method" 
```

and verify that you're at the page that describes the method relevant
to an `IRanges` instance.

The [GenomicRanges][] package extends the notion of ranges to include
features relevant to application of ranges in sequence analysis,
particularly the ability to associate a range with a sequence name
(e.g., chromosome) and a strand. Create a `GRanges` instance based on
our `IRanges` instance, as follows
```{r granges}
library(GenomicRanges)
gr <- GRanges(c("chr1", "chr1", "chr2"), ir, strand=c("+", "-", "+"))
gr
```

The notion of flanking sequence has a more nuanced meaning in
biology. In particular we might expect that flanking sequence on the
`+` strand would precede the range, but on the minus strand would
follow it. Verify that `flank` applied to a `GRanges` object has this
behavior.
```{r granges-flank}
flank(gr, 3)
```

Discover what classes `GRanges` extends, find the help page
documenting the behavior of `flank` when applied to a GRanges object,
and verify that the help page documents the behavior we just observed.
```{r granges-class}
class(gr)
getClassDef(class(gr))
```
```{r granges-flank-method, eval=FALSE}
?"flank,GenomicRanges-method"
```
Notice that the available `flank()` methods have been augmented by the
methods defined in the _GenomicRanges_ package.

It seems like there might be a number of helpful methods available for
working with genomic ranges; we can discover some of these from the
command line, indicating that the methods should be on the current
`search()` path

```{r granges-methods, eval=FALSE}
showMethods(class="GRanges", where=search())
```

Use `help()` to list the help pages in the `GenomicRanges` package,
and `vignettes()` to view and access available vignettes; these are
also available in the Rstudio 'Help' tab.
```{r granges-man-and-vignettes, eval=FALSE}
help(package="GenomicRanges")
vignette(package="GenomicRanges")
vignette(package="GenomicRanges", "GenomicRangesHOWTOs")
```

<a name="package-tour"></a>

## 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
```{r require}
  library(GenomicRanges)
```
  and the help system queried interactively, as outlined above:
```{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.
- Popular ChIP-seq packages include [DiffBind][] for comparison of
  peaks across samples, [ChIPQC][] for quality assessment, and
  [ChIPpeakAnno][] 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")`
- [rtracklayer][]'s `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. Explore the
  [ShortRead vignette](http://bioconductor.org/packages/release/bioc/vignettes/ShortRead/inst/doc/Overview.pdf)
  and Scalable Genomics labs to see approaches to effectively
  processing the large files.

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.hg19.knownGene][])
    contain gene models (exon coordinates, exon / transcript
    relationships, etc) derived from common sources such as the hg19
    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.
  - **BSgenome** packages (e.g., [BSgenome.Hsapiens.UCSC.hg19][])
    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 (we'll do more of this in Thursday's lab) at the
  [annotation work flow](http://bioconductor.org/help/workflows/annotation/annotation/)
  on the Bioconductor web site.

## Summary

Bioconductor is a large collection of R packages for the analysis and
comprehension of high-throughput genomic data. Bioconductor relies on
formal classes to represent genomic data, so it is important to
develop a rudimentary comfort with classes, including seeking help for
classes and methods. Bioconductor uses vignettes to augment
traditional help pages; these can be very valuable in illustrating
overall package use.


[AnnotationData]: http://bioconductor.org/packages/release/BiocViews.html#___AnnotationData
[AnnotationDbi]: http://bioconductor.org/packages/release/bioc/html/AnnotationDbi.html
[AnnotationHub]: http://bioconductor.org/packages/release/bioc/html/AnnotationHub.html
[BSgenome.Hsapiens.UCSC.hg19]: http://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg19.html
[BSgenome]: http://bioconductor.org/packages/release/bioc/html/BSgenome.html
[BiocParallel]: http://bioconductor.org/packages/release/bioc/html/BiocParallel.html
[Biostrings]: http://bioconductor.org/packages/release/bioc/html/Biostrings.html
[Bsgenome.Hsapiens.UCSC.hg19]: http://bioconductor.org/packages/release/data/annotation/html/Bsgenome.Hsapiens.UCSC.hg19.html
[CNTools]: http://bioconductor.org/packages/release/bioc/html/CNTools.html
[ChIPQC]: http://bioconductor.org/packages/release/bioc/html/ChIPQC.html
[ChIPpeakAnno]: http://bioconductor.org/packages/release/bioc/html/ChIPpeakAnno.html
[DESeq2]: http://bioconductor.org/packages/release/bioc/html/DESeq2.html
[DiffBind]: http://bioconductor.org/packages/release/bioc/html/DiffBind.html
[GenomicAlignments]: http://bioconductor.org/packages/release/bioc/html/GenomicAlignments.html
[GenomicRanges]: http://bioconductor.org/packages/release/bioc/html/GenomicRanges.html
[Homo.sapiens]: http://bioconductor.org/packages/release/data/annotation/html/Homo.sapiens.html
[IRanges]: http://bioconductor.org/packages/release/bioc/html/IRanges.html
[KEGGREST]: http://bioconductor.org/packages/release/bioc/html/KEGGREST.html
[PSICQUIC]: http://bioconductor.org/packages/release/bioc/html/PSICQUIC.html
[Rsamtools]: http://bioconductor.org/packages/release/bioc/html/Rsamtools.html
[ShortRead]: http://bioconductor.org/packages/release/bioc/html/ShortRead.html
[TxDb.Hsapiens.UCSC.hg19.knownGene]: http://bioconductor.org/packages/release/data/annotation/html/TxDb.Hsapiens.UCSC.hg19.knownGene.html
[VariantAnnotation]: http://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html
[VariantFiltering]: http://bioconductor.org/packages/release/bioc/html/VariantFiltering.html
[VariantTools]: http://bioconductor.org/packages/release/bioc/html/VariantTools.html
[biocViews]: http://bioconductor.org/packages/release/BiocViews.html#___Software
[biomaRt]: http://bioconductor.org/packages/release/bioc/html/biomaRt.html
[cn.mops]: http://bioconductor.org/packages/release/bioc/html/cn.mops.html
[edgeR]: http://bioconductor.org/packages/release/bioc/html/edgeR.html
[ensemblVEP]: http://bioconductor.org/packages/release/bioc/html/ensemblVEP.html 
[h5vc]: http://bioconductor.org/packages/release/bioc/html/h5vc.html
[limma]: http://bioconductor.org/packages/release/bioc/html/limma.html
[metagenomeSeq]: http://bioconductor.org/packages/release/bioc/html/metagenomeSeq.html
[org.Hs.eg.db]: http://bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html
[org.Sc.sgd.db]: http://bioconductor.org/packages/release/data/annotation/html/org.Sc.sgd.db.html
[phyloseq]: http://bioconductor.org/packages/release/bioc/html/phyloseq.html
[rtracklayer]: http://bioconductor.org/packages/release/bioc/html/rtracklayer.html
[snpStats]: http://bioconductor.org/packages/release/bioc/html/snpStats.html