---
title: "Differential expression, manipulation, and visualization of RNA-seq reads"
author: "Michael Love, Simon Anders, Wolfgang Huber"
date: "21 July 2015"
vignette: >
%\VignetteIndexEntry{Differential expression, manipulation, and visualization of RNA-seq reads}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r opts, echo=FALSE, results="hide"}
library("knitr")
opts_chunk$set(error=FALSE)
```
```{r style, echo=FALSE, results="asis"}
library("BiocStyle")
BiocStyle::markdown()
```
# Differential expression, manipulation, and visualization of RNA-seq reads
Michael Love [1], Simon Anders [2], Wolfgang Huber [3]
[1] Dept. of Biostatistics, Dana-Farber Cancer Institute &
Harvard T.H.Chan School of Public Health, Boston, US;
[2] Institute for Molecular Medicine Finland (FIMM), Helsinki, Finland;
[3] European Molecular Biology Laboratory (EMBL), Heidelberg, Germany.
### Link to more material online
This workflow is an abridged/modified version of the longer
[RNA-seq workflow](http://bioconductor.org/help/workflows/rnaseqGene/)
available under [Workflows](http://bioconductor.org/help/workflows/)
on the Bioconductor website.
## Contents
Here we will cover basic steps in RNA-seq analysis
using a variety of Bioconductor packages:
* [Loading gene annotations](#gene) (GenomicFeatures, AnnotationHub)
* [Creating an RNA-seq count table](#count) (GenomicAlignments, Rsubread)
* [Exploratory analysis and visualization](#eda) (DESeq2, PoiClaClu)
* [Differential expression testing](#de) (DESeq2)
* [Annotation of results tables](#annotate) (AnnotationDbi)
* [Generation of HTML reports](#report) (ReportingTools)
* [Examining RNA-seq alignments](#align) (GenomicAlignments, Gviz)
## Introduction
This lab will walk you through an RNA-Seq differential expression workflow,
using `r Biocpkg("DESeq2")` along with other Bioconductor packages.
We note that a number of other Bioconductor packages can also be used for
statistical inference of differential expression at the gene level including
`r Biocpkg("edgeR")`, `r Biocpkg("BaySeq")`, `r Biocpkg("DSS")` and
`r Biocpkg("limma")`.
## Finding help
Remember, if you have a question about a particular function,
for example, `DESeq()`, you can read about the arguments and details with:
```{r, eval=FALSE}
?DESeq
```
All Bioconductor packages also come with descriptive vignettes.
The DESeq2 vignette can be accessed with:
```{r, eval=FALSE}
browseVignettes("DESeq2")
```
If you have questions about this workflow, please post them to the
[Bioconductor support site](https://support.bioconductor.org/).
If the questions concern a specific package, you can tag the post with the name of the package,
or for general questions about the workflow, tag the post with `deseq2`. Note the
[posting guide](http://www.bioconductor.org/help/support/posting-guide/)
for crafting an optimal question for the support site.
## Experimental data
The data used in this workflow consists of RNA-seq samples of
human **airway smooth muscle cells** treated with dexamethasone,
a synthetic **glucocorticoid steroid** with
anti-inflammatory effects. Glucocorticoids
are used, for example, in asthma patients to prevent or reduce
inflammation of the airways. In this experiment, we have two
samples each from **four cell lines**:
one sample was untreated, while another sample was treated with 1 micromolar
dexamethasone for 18 hours. The publication for the experiment is:
Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM,
Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr,
Tantisira KG, Weiss ST, Lu Q. "RNA-Seq Transcriptome Profiling
Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates
Cytokine Function in Airway Smooth Muscle Cells." PLoS One. 2014 Jun
13;9(6):e99625.
PMID: [24926665](http://www.ncbi.nlm.nih.gov/pubmed/24926665).
GEO: [GSE52778](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778).
## The biological question
Our aim in this workflow is to identify **differences in
gene expression due to treatment**, compared to control samples,
while controlling for the differences across cell line.
## Count matrix
As input, the *DESeq2* package expects *count data* as obtained,
e.g., from RNA-seq or another high-throughput sequencing experiment,
in the form of a matrix of integer values. The value in the *i*-th row
and the *j*-th column of the matrix tells how many reads
(or fragments in a paired-end experiment) can be
uniquely assigned to gene *i* in sample *j*.
The count values must be raw counts of sequencing reads. This is
important for *DESeq2*'s statistical model to hold, as only the
actual counts allow assessing the measurement precision
correctly.
For more details on why we work with counts, or the statistical
model underlying the methods, see the
[DESeq2 paper](http://www.genomebiology.com/2014/15/12/550).
## Aligning reads
We will assume we have already aligned the reads to the genome
and have BAM files for each of the samples. In order to see
the code which was used for this experiment to align the paired-end
reads to the genome, see the
[full workflow](http://bioconductor.org/help/workflows/rnaseqGene/)
or the [vignette](http://bioconductor.org/packages/release/data/experiment/vignettes/airway/inst/doc/airway.html)
for the airway data package.
## Loading gene annotations
First we will need to define the genes in which we will
count the RNA-seq fragments. Our initial goal is to
construct a *TxDb* object, which contains all the information
about the exons, transcripts and genes of an organism.
From this we can pull out the exons associated with each
gene, which will be used for counting.
We have different options for how we read in the gene model information.
One option is to read the gene model from a
[GTF file](http://www.ensembl.org/info/website/upload/gff.html), using
*makeTxDbFromGFF* from the `r Biocpkg("GenomicFeatures")`
package. GTF files can be downloaded from
Ensembl's FTP site or other gene model repositories.
We first need to load our data package, which contains the
GTF file (and small example BAM files and count matrix).
```{r message=FALSE}
library("airway")
```
The R function *system.file* can be used to find out where on your
computer the files from a package have been installed. Here we ask for
the full path to the `extdata` directory, which is part of the
`r Biocexptpkg("airway")` package.
Note that this GTF file only contains a subset of all the genes
from the 75th Ensembl release.
```{r}
dir <- system.file("extdata", package="airway", mustWork=TRUE)
gtffile <- file.path(dir, "Homo_sapiens.GRCh37.75_subset.gtf")
```
Now we can create the *TxDb*.
```{r message=FALSE}
library("GenomicFeatures")
txdb <- makeTxDbFromGFF(gtffile, format="gtf")
```
Now we can produce a *GRangesList* of all the exons grouped by gene.
Each element of the list defines all the exons associated with a given
gene.
```{r}
ebg <- exonsBy(txdb, by="gene")
ebg[[1]]
```
An even faster option for creating a *TxDb*,
compared to reading in GTF files, is to
pull down the pre-processed GTF file using `r Biocpkg("AnnotationHub")`.
Below we show an example of how you could construct a
*TxDb* for Saccharomyces cerevisiae from the Ensembl resource.
Browse the *AnnotationHub* vignettes for more details on usage.
```{r annohub, message=FALSE, warning=FALSE}
library("AnnotationHub")
ah <- AnnotationHub()
```
```{r}
qu <- query(ah, c("Ensembl","gtf","Saccharomyces cerevisiae","release-80"))
stopifnot(length(qu) == 1) # check there was only one match
gtf <- qu[[1]] # this actually downloads the gene model
txdbFromAHub <- makeTxDbFromGRanges(gtf)
```
There are still other options for obtaining a *TxDb*.
For the *known genes* track from the UCSC Genome Browser,
one can use the pre-built Bioconductor package:
`r Biocannopkg("TxDb.Hsapiens.UCSC.hg19.knownGene")`.
The *makeTxDbFromBiomart* and *makeTxDbFromUCSC* functions
can be used to automatically pull
gene models from Biomart and UCSC, respectively.
See the `r Biocpkg("GenomicFeatures")` package for more details.
## Creating an RNA-seq count table
We will now show how to count reads/fragments which can be uniquely assigned
to the exons of a gene, across all the genes and all the samples.
Besides the main count matrix, which we will use later, the
`r Biocexptpkg("airway")` package also contains a small subset
of reads from eight *BAM files* (the reads which overlap the
genes in our subset GTF file).
We will use these BAM files to demonstrate how a count matrix can be constructed.
Afterwards, we will load the full count matrix
corresponding to all samples and all data, which is already provided
in the same package, and will continue the analysis with that full
table.
In our directory specified earlier, we find the eight BAM files (and some other files):
```{r}
list.files(dir)
```
Typically, we have a **sample table** with experimental metadata for our samples
(which samples were in treatment vs control, which were in batch 1 vs batch 2, etc.).
For your own project, you might create such a comma-separated
value (CSV) file using a text editor or spreadsheet software such as Excel.=
We can load such a file with the function *read.csv*.
(The parentheses around the last line are used
to print the result in addition to storing it to the
`sampleTable` object.)
```{r}
csvfile <- file.path(dir,"sample_table.csv")
( sampleTable <- read.csv(csvfile,row.names=1) )
```
Using the `Run` column in the sample table, we construct the full
paths to the files we want to perform the counting operation on:
```{r}
filenames <- file.path(dir, paste0(sampleTable$Run, "_subset.bam"))
```
We indicate in Bioconductor that these files are BAM files using the
*BamFileList* function. Here we can optionally also
specify details about how the BAM files should
be treated, e.g., only process 2 million reads at a time.
```{r message=FALSE}
library("Rsamtools")
bamfiles <- BamFileList(filenames, yieldSize=2000000)
```
**Note:** make sure that the chromosome names of the genomic features
in the annotation you use are consistent with the chromosome names of
the reference used for read alignment. Otherwise, the scripts might
fail to count any reads to features due to the mismatching names.
We can check the chromosome names in the alignment files like so:
```{r}
seqinfo(bamfiles[1])
```
After these preparations, the actual counting is easy. The function
*summarizeOverlaps* from the `r Biocpkg("GenomicAlignments")`
package will do this. This produces a *SummarizedExperiment*
object, which contains a variety of information about
an experiment, and will be described in more detail below.
**Note:** If it is desired to perform counting using multiple cores, one can use
the *register* and *MulticoreParam* functions from the
`r Biocpkg("BiocParallel")` package before the counting call below.
```{r message=FALSE}
library("GenomicAlignments")
```
The result of the counting call will be a *SummarizedExperiment*
object, explained in a section below.
```{r}
se <- summarizeOverlaps(features=ebg, reads=bamfiles,
mode="Union",
singleEnd=FALSE,
ignore.strand=TRUE,
fragments=TRUE )
```
We specify a number of arguments besides the `features` and the
`reads`. The `mode` argument describes what kind of read overlaps will
be counted as a hit. These modes are shown in Figure 1 of the "Counting reads with
summarizeOverlaps" vignette for the `r Biocpkg("GenomicAlignments")`
package. Setting `singleEnd` to `FALSE` indicates that the experiment
produced paired-end reads, and we want to count a pair of reads only once
toward the read count for a gene.
In order to produce correct counts, it is important to know if the
RNA-Seq experiment was strand-specific or not. This experiment was not
strand-specific so we set `ignore.strand` to `TRUE`.
The `fragments` argument can be used when `singleEnd=FALSE`
to specify if unpaired reads should be counted (yes if `fragments=TRUE`).
## SummarizedExperiment
Here we show the component parts of a *SummarizedExperiment* object,
and also its subclasses, such as the *DESeqDataSet* which is
explained in the next section. The `assay(s)` (pink block) contains
the matrix (or matrices) of summarized values, the `rowRanges` (blue
block) contains information about the genomic ranges, and the
`colData` (green block) contains information about the samples or
experiments. The highlighted line in each block represents the first
row (note that the first row of `colData` lines up with the first
column of the `assay`.
```{r sumexp, echo=FALSE}
par(mar=c(0,0,0,0))
plot(1,1,xlim=c(0,100),ylim=c(0,100),bty="n",
type="n",xlab="",ylab="",xaxt="n",yaxt="n")
polygon(c(45,80,80,45),c(10,10,70,70),col="pink",border=NA)
polygon(c(45,80,80,45),c(68,68,70,70),col="pink3",border=NA)
text(62.5,40,"assay(s)")
text(62.5,30,"e.g. 'counts'")
polygon(c(20,40,40,20),c(10,10,70,70),col="skyblue",border=NA)
polygon(c(20,40,40,20),c(68,68,70,70),col="skyblue3",border=NA)
text(30,40,"rowRanges")
polygon(c(45,80,80,45),c(75,75,90,90),col="palegreen",border=NA)
polygon(c(45,47,47,45),c(75,75,90,90),col="palegreen3",border=NA)
text(62.5,82.5,"colData")
```
The example code above actually only counts a small subset of reads
from the original experiment. Nevertheless, we
can still investigate the resulting *SummarizedExperiment* by looking at
the counts in the `assay` slot, the phenotypic data about the samples in
`colData` slot (in this case an empty *DataFrame*), and the data about the
genes in the `rowRanges` slot.
```{r}
se
head(assay(se))
colSums(assay(se))
colData(se)
rowRanges(se)
```
Note that the `rowRanges` slot is a *GRangesList*, which
contains all the information about the exons for each gene, i.e., for each row
of the count matrix. It also contains metadata about the construction
of the gene model in the `metadata` slot.
```{r}
str(metadata(rowRanges(se)))
```
The `colData` slot, so far empty, should contain all the metadata.
We hence assign our sample table to it:
```{r}
colData(se) <- DataFrame(sampleTable)
```
## featureCounts from Rsubread
We could have also counted the fragments which align to each gene
using the *featureCounts* function from the `r Biocpkg("Rsubread")`
package. **Note:** this package only exists for Linux and Mac
machines, not for Windows machines.
Here we show an example of counting, using only the
first file from the vector of 8 files. This function
produces a *list* with the count matrix called *counts*.
```{r message=FALSE}
library("Rsubread")
fc <- featureCounts(filenames[1], annot.ext=gtffile,
isGTFAnnotationFile=TRUE, isPairedEnd=TRUE)
names(fc)
head(fc$counts)
```
## Branching point for statistical analysis
At this point, we have counted the reads which overlap the genes in
the gene model we specified. This is a branching point where we could
use a variety of Bioconductor packages for exploration and
differential expression of the counts, including
`r Biocpkg("edgeR")`, `r Biocpkg("BaySeq")`,
`r Biocpkg("DSS")` and `r Biocpkg("limma")`.
We will continue, using `r Biocpkg("DESeq2")`.
The *SummarizedExperiment* object (or count matrix and sample table)
is all we need to start our analysis. In the following section we will
show how to use it to create the data object used by `r Biocpkg("DESeq2")`.
## The *DESeqDataSet*, sample information, and the design formula
Bioconductor software packages often define and use a custom class for
their data object, which makes sure that all the needed data slots are
consistently provided and fulfill the requirements. In addition,
Bioconductor has general data classes (such as the
*SummarizedExperiment*) that can be used to move data between
packages. In `r Biocpkg("DESeq2")`, the custom class is called
*DESeqDataSet*. It is built on top of the *SummarizedExperiment* class
(in technical terms, *DESeqDataSet* is a subclass), and it is easy
to convert *SummarizedExperiment* instances into *DESeqDataSet* and vice versa.
One of the main differences is that the `assay` slot is instead
accessed using the *count* accessor, and the class enforces that the
values in this matrix are non-negative integers.
A second difference is that the *DESeqDataSet* has an associated
**design formula**. The experimental design is specified at the
beginning of the analysis, as it will inform many of the *DESeq2*
functions how to treat the samples in the analysis (one exception is
the size factor estimation, i.e., the adjustment for differing library
sizes, which does not depend on the design formula). The design
formula tells which variables in the column metadata (`colData`)
specify the experimental design and how these factors should be used
in the analysis. The column metadata is where the information about the
samples (columns of the count matrix) are stored.
The simplest design formula for differential expression would be
`~ condition`, where `condition` is a column in `colData(dds)` which
specifies which of two (or more groups) the samples belong to. For
the airway experiment, we will specify `~ cell + dex`, which
means that we want to test for the effect of dexamethasone (the last
factor), controlling for the effect of different cell line (the first
factor).
You can use R's formula notation to express any experimental design
that can be described within an ANOVA-like framework. Note that
*DESeq2* uses the same formula notation as, for instance, the
linear model (*lm*)
function of base R. If the question of interest is whether a fold
change due to treatment is different across groups, interaction terms
can be included using models such as
`~ group + treatment + group:treatment`. See the manual page for
`?results` for examples of extracting contrasts from more complex
designs such as these.
In the following sections, we will demonstrate the construction of the
*DESeqDataSet* from two starting points:
* from a *SummarizedExperiment* object created by, e.g.,
*summarizeOverlaps* in the above example
* more generally, from a count matrix and a sample information table
which have been loaded into R
For a full example of using the *HTSeq* Python package for read
counting, please see the `r Biocexptpkg("pasilla")` vignette. For an
example of generating the *DESeqDataSet* from files produced by
*htseq-count*, please see the `r Biocpkg("DESeq2")` vignette.
## Starting from *SummarizedExperiment*
We now use R's *data* command to load a prepared
*SummarizedExperiment* that was generated from the publicly available
sequencing data files associated with the Himes et al. paper,
described above. The steps we used to produce this object were
equivalent to those you worked through in the previous sections,
except that we used all the reads and all the genes. For more details
on the exact steps used to create this object type
`browseVignettes("airway")` into your R session.
```{r}
data("airway")
se <- airway
```
We can quickly check the millions of fragments which uniquely aligned
to the genes (the second argument of *round* tells how many decimal
points to keep).
```{r}
round( colSums(assay(se)) / 1e6, 1 )
```
Supposing we have constructed a *SummarizedExperiment* using
one of the methods described in the previous section, we now need to
make sure that the object contains all the necessary information about
the samples, i.e., a table with information on the count matrix's columns
stored in the `colData` slot:
```{r}
colData(se)
```
Here we see that this object already contains an informative
`colData` slot -- because we have already prepared it for you, as
described in the `r Biocexptpkg("airway")` vignette.
However, when you work with your own data, you will have to add the
pertinent sample / phenotypic information for the experiment at this stage.
We highly recommend keeping this information in a comma-separated
value (CSV) or tab-separated value (TSV) file, which can be exported
from an Excel spreadsheet, and the assign this to the `colData` slot,
making sure that the rows correspond to the columns of the
*SummarizedExperiment*. We made sure of this correspondence by
specifying the BAM files using a column of the sample table.
Once we have our fully annotated *SummarizedExperiment* object,
we can construct a *DESeqDataSet* object from it, which will then form
the starting point of the actual *DESeq2* package, described in the
following sections. We add an appropriate design for the analysis.
```{r message=FALSE}
library("DESeq2")
```
```{r}
dds <- DESeqDataSet(se, design = ~ cell + dex)
```
If we only wanted to perform transformations and exploratory data analysis
we could use a `~ 1` for the design, but be careful, because a true
experimental design, e.g. `~ condition` would need to be added
later before differential expression (or else we would only be
testing the intercept).
Note that there are two alternative functions,
*DESeqDataSetFromMatrix* and *DESeqDataSetFromHTSeq*, which allow you
to get started in case you have your data not in the form of a
*SummarizedExperiment* object, but either as a simple matrix of count
values or as output files from the *htseq-count* script from the
*HTSeq* Python package.
Below we demonstrate using *DESeqDataSetFromMatrix*.
## Starting from count matrices
In this section, we will show how to build an *DESeqDataSet* supposing
we only have a count matrix and a table of sample information.
**Note:** if you have prepared a *SummarizedExperiment* you should skip this
section. While the previous section would be used to contruct a
*DESeqDataSet* from a *SummarizedExperiment*, here we first extract
the individual object (count matrix and sample info) from the
*SummarizedExperiment* in order to build it back up into a new object
-- only for demonstration purposes.
In practice, the count matrix would either be read in from a file or
perhaps generated by an R function like *featureCounts* from the
`r Biocpkg("Rsubread")` package.
The information in a *SummarizedExperiment* object can be
accessed with accessor functions. For example, to see the actual data,
i.e., here, the read counts, we use the *assay* function. (The *head*
function restricts the output to the first few lines.)
```{r}
countdata <- assay(se)
head(countdata)
```
In this count matrix, each row represents an Ensembl gene, each column
a sequenced RNA library, and the values give the raw numbers of
sequencing reads that were mapped to the respective gene in each
library. We also have information on each of the samples (the columns of the
count matrix). If you've counted reads with some other software, you need to check at
this step that the columns of the count matrix correspond to the rows
of the column metadata.
```{r}
coldata <- colData(se)
```
We now have all the ingredients to prepare our data object in a form
that is suitable for analysis, namely:
* `countdata`: a table with the read counts
* `coldata`: a table with information about the samples
To now construct the data object from the matrix of counts and the
sample information table, we use:
```{r}
(ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ cell + dex))
```
We will continue with the object generated from the
*SummarizedExperiment* section.
# Exploratory analysis and visualization
There are two separate paths in this workflow, the one we will
see first involves **transformations of the counts**
in order to visualizally explore sample relationships.
In the second part, we will go back to the orignal raw counts
for **statistical testing**. This is critical because
the statistical testing methods rely on original data
(not scaled or transformed) for calculating the precision of measurements.
## The rlog transformation
Many common statistical methods for exploratory analysis of
multidimensional data, especially methods for clustering and
ordination (e.g., principal-component analysis and the like), work
best for (at least approximately) homoskedastic data; this means that
the variance of an observed quantity (here, the expression
strength of a gene) does not depend on the mean. In RNA-Seq data,
however, variance grows with the mean. For example, if one performs
PCA (principal components analysis) directly on a matrix of normalized
read counts, the result typically depends only on the few most
strongly expressed genes because they show the largest absolute
differences between samples. A simple and often used strategy to avoid
this is to take the logarithm of the normalized count values plus a
small pseudocount; however, now the genes with low counts tend to
dominate the results because, due to the strong Poisson noise inherent
to small count values, they show the strongest relative differences
between samples.
As a solution, *DESeq2* offers the *regularized-logarithm transformation*,
or *rlog* for short. For genes with high counts, the rlog
transformation differs not much from an ordinary log2
transformation. For genes with lower counts, however, the values are
shrunken towards the genes' averages across all samples. Using
an empirical Bayesian prior on inter-sample differences in the form of
a *ridge penalty*, this is done such that the rlog-transformed data
are approximately homoskedastic. See the help for `?rlog` for more
information and options. Another transformation, the *variance
stabilizing transformation*, is discussed alongside the *rlog* in the
*DESeq2* vignette.
**Note:** the rlog transformation is provided for applications *other*
than differential testing. For differential testing we recommend the
*DESeq* function applied to raw counts, as described later
in this workflow, which also takes into account the dependence of the
variance of counts on the mean value during the dispersion estimation
step.
The function *rlog* returns a *SummarizedExperiment*
object which contains the rlog-transformed values in its *assay* slot:
```{r}
rld <- rlog(dds)
head(assay(rld))
```
To show the effect of the transformation, we plot the first sample
against the second, first simply using the *log2* function (after adding
1, to avoid taking the log of zero), and then using the rlog-transformed
values. For the *log2* method, we need estimate size factors to
account for sequencing depth (this is done automatically for the
*rlog* method).
```{r rldplot, fig.width=8, fig.height=4}
par( mfrow = c( 1, 2 ) )
dds <- estimateSizeFactors(dds)
plot(log2( 1 + counts(dds, normalized=TRUE)[ , 1:2] ),
pch=16, cex=0.3)
plot(assay(rld)[ , 1:2],
pch=16, cex=0.3)
```
Note that, in order to make it easier to see where several points are
plotted on top of each other, we set the plotting color to a
semi-transparent black and changed the points to solid circles
(`pch=16`) with reduced size (`cex=0.3`).
We can see how genes with low counts seem to be excessively variable
on the ordinary logarithmic scale, while the rlog transform compresses
differences for genes for which the data cannot provide good information anyway.
## Sample distances
A useful first step in an RNA-Seq analysis is often to assess overall
similarity between samples: Which samples are similar to each other,
which are different? Does this fit to the expectation from the
experiment's design?
We use the R function *dist* to calculate the Euclidean distance
between samples. To avoid that the distance measure is dominated by a
few highly variable genes, and have a roughly equal contribution from
all genes, we use it on the rlog-transformed data:
```{r}
sampleDists <- dist( t( assay(rld) ) )
sampleDists
```
Note the use of the function *t* to transpose the data matrix. We need
this because *dist* calculates distances between data *rows* and
our samples constitute the columns.
We visualize the distances in a heatmap, using the function
*pheatmap* from the `r CRANpkg("pheatmap")` package.
```{r message=FALSE}
library("pheatmap")
library("RColorBrewer")
```
In order to plot the sample distance matrix with the rows/columns
arranged by those distances in the matrix,
we manually provide the `sampleDists` to the `clustering_distance`
argument of the *pheatmap* function.
Otherwise the *pheatmap* function would assume that the matrix contains
the data values themselves, and would calculate distances between the
rows/columns of the distance matrix, which is not desired.
```{r distheatmap, fig.width=8}
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- paste( rld$dex, rld$cell, sep="-" )
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors)
```
Note that we have changed the row names of the distance matrix to
contain treatment type and patient number instead of sample ID, so
that we have all this information in view when looking at the heatmap.
Another option for calculating sample distances is to use the Poisson
Distance, implemented in the CRAN package
`r CRANpkg("PoiClaClu")`. Similar to the transformations offered in
*DESeq2*, this measure of dissimilarity also takes the variance
structure of counts into consideration when calculating the distances
between samples. The *PoissonDistance* function takes the original
count matrix (not normalized) with samples as rows instead of columns,
so we need to transpose the counts in `dds`.
```{r message=FALSE}
library("PoiClaClu")
poisd <- PoissonDistance(t(counts(dds)))
```
We can plot the heatmap as before:
```{r poisdistheatmap, fig.width=8}
samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- paste( rld$dex, rld$cell, sep="-" )
colnames(samplePoisDistMatrix) <- NULL
pheatmap(samplePoisDistMatrix,
clustering_distance_rows=poisd$dd,
clustering_distance_cols=poisd$dd,
col=colors)
```
## PCA plot
Another way to visualize sample-to-sample distances is a
principal-components analysis (PCA). In this ordination method, the
data points (i.e., here, the samples) are projected onto the 2D plane
such that they spread out in the two directions which explain most of
the differences in the data. The x-axis is the direction (or principal
component) which separates the data points the most. The amount of the
total variance which is contained in the direction is printed in the
axis label.
```{r plotpca, fig.width=6, fig.height=4.5}
plotPCA(rld, intgroup = c("dex", "cell"))
```
Here, we have used the function *plotPCA* which comes with *DESeq2*.
The two terms specified by `intgroup` are the interesting groups for
labeling the samples; they tell the function to use them to choose
colors. We can also build the PCA plot from scratch using
`r CRANpkg("ggplot2")`. This is done by asking the *plotPCA* function
to return the data used for plotting rather than building the plot.
See the *ggplot2* [documentation](http://docs.ggplot2.org/current/)
for more details on using *ggplot*.
```{r}
(data <- plotPCA(rld, intgroup = c( "dex", "cell"), returnData=TRUE))
percentVar <- round(100 * attr(data, "percentVar"))
```
We can then use this data to build up the plot, specifying that the
color of the points should reflect dexamethasone treatment and the
shape should reflect the cell line.
```{r message=FALSE}
library("ggplot2")
```
```{r ggplotpca, fig.width=6, fig.height=4.5}
ggplot(data, aes(PC1, PC2, color=dex, shape=cell)) + geom_point(size=3) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance"))
```
From both visualizations, we see that the differences between cells are
considerable, though not stronger than the differences due to
treatment with dexamethasone. This shows why it will be important to
account for this in differential testing by using a paired design
("paired", because each dex treated sample is paired with one
untreated sample from the *same* cell line). We are already set up for
this by using the design formula `~ cell + dex` when setting up the
data object in the beginning.
## MDS plot
Another plot, very similar to the PCA plot, can be made using the
*multidimensional scaling* (MDS) function in base R. This is useful when we
don't have the original data, but only a matrix of distances. Here we
have the MDS plot for the distances calculated from the *rlog*
transformed counts:
```{r mdsrlog, fig.width=6, fig.height=4.5}
mds <- data.frame(cmdscale(sampleDistMatrix))
mds <- cbind(mds, as.data.frame(colData(rld)))
ggplot(mds, aes(X1,X2,color=dex,shape=cell)) + geom_point(size=3)
```
And here from the *PoissonDistance*:
```{r mdspois, fig.width=6, fig.height=4.5}
mds <- data.frame(cmdscale(samplePoisDistMatrix))
mds <- cbind(mds, as.data.frame(colData(dds)))
ggplot(mds, aes(X1,X2,color=dex,shape=cell)) + geom_point(size=3)
```
# Differential gene expression testing
Now, suppose we want to do statistical testing on this dataset
to investigate which, if any, genes show evidence of
differential expression due to treatment, and controlling
for cell line.
It will be convenient to make sure that `untrt` is the first level in
the `dex` factor, so that the default log2 fold changes are calculated
as treated over untreated (by default R will chose the first
alphabetical level, remember: computers don't know what to do unless
you tell them). The function *relevel* achieves this:
```{r}
dds$dex <- relevel(dds$dex, "untrt")
```
## Running the pipeline
Finally, we are ready to run the differential expression pipeline.
With the data object prepared, the *DESeq2* analysis can now be run
with a single call to the function *DESeq*:
```{r}
dds <- DESeq(dds)
```
This function will print out a message for the various steps it
performs. These are described in more detail in the manual page for
*DESeq*, which can be accessed by typing `?DESeq`. Briefly these are:
the estimation of size factors (which control for differences in the
library size of the sequencing experiments), the estimation of
dispersion for each gene, and fitting a generalized linear model.
A *DESeqDataSet* is returned which contains all the fitted
information within it, and the following section describes how to
extract out results tables of interest from this object.
## Building the results table
Calling *results* without any arguments will extract the estimated
log2 fold changes and *p* values for the last variable in the design
formula. If there are more than 2 levels for this variable, *results*
will extract the results table for a comparison of the last level over
the first level.
```{r}
( res <- results(dds) )
```
As `res` is a *DataFrame* object, it carries metadata
with information on the meaning of the columns:
```{r}
mcols(res, use.names=TRUE)
```
The first column, `baseMean`, is a just the average of the normalized
count values, dividing by size factors, taken over all samples. The
remaining four columns refer to a specific contrast, namely the
comparison of the `trt` level over the `untrt` level for the factor
variable `dex`. See the help page for *results* (by typing `?results`)
for information on how to obtain other contrasts.
The column `log2FoldChange` is the effect size estimate. It tells us
how much the gene's expression seems to have changed due to treatment
with dexamethasone in comparison to untreated samples. This value is
reported on a logarithmic scale to base 2: for example, a log2 fold
change of 1.5 means that the gene's expression is increased by a
multiplicative factor of 2^1.5 or approximately 2.82.
Of course, this estimate has an uncertainty associated with it, which
is available in the column `lfcSE`, the standard error estimate for
the log2 fold change estimate. We can also express the uncertainty of
a particular effect size estimate as the result of a statistical
test. The purpose of a test for differential expression is to test
whether the data provides sufficient evidence to conclude that this
value is really different from zero. *DESeq2* performs for each gene a
*hypothesis test* to see whether evidence is sufficient to decide
against the *null hypothesis* that there is no effect of the treatment
on the gene and that the observed difference between treatment and
control was merely caused by experimental variability (i.e., the type
of variability that you can just as well expect between different
samples in the same treatment group). As usual in statistics, the
result of this test is reported as a *p* value, and it is found in the
column `pvalue`. (Remember that a *p* value indicates the probability
that a fold change as strong as the observed one, or even stronger,
would be seen under the situation described by the null hypothesis.)
We can also summarize the results with the following line of code,
which reports some additional information, which will be covered in
later sections.
```{r}
summary(res)
```
Note that there are many genes with differential expression due to
dexamethasone treatment at the FDR level of 10%. This makes sense, as
the smooth muscle cells of the airway are known to react to
glucocorticoid steroids. However, there are two ways to be more strict
about which set of genes are considered significant:
* lower the false discovery rate threshold (the threshold on `padj` in
the results table)
* raise the log2 fold change threshold from 0 using the `lfcThreshold`
argument of *results*. See the *DESeq2* vignette for a demonstration
of the use of this argument.
If we lower the false discovery rate threshold, we should also
tell this value to `results()`, so that the function will use an
alternative threshold for the optimal independent filtering step:
```{r}
res.05 <- results(dds, alpha=.05)
table(res.05$padj < .05)
```
Sometimes a subset of the *p* values in `res` will be `NA` ("not
available"). This is *DESeq*'s way of reporting that all counts for
this gene were zero, and hence not test was applied. In addition, *p*
values can be assigned `NA` if the gene was excluded from analysis
because it contained an extreme count outlier. For more information,
see the outlier detection section of the vignette.
## Other comparisons
In general, the results for a comparison of any two levels of a
variable can be extracted using the `contrast` argument to
*results*. The user should specify three values: the name of the
variable, the name of the level in the numerator, and the name of the
level in the denominator. Here we extract results for the log2 of the
fold change of one cell line over another:
```{r}
( res.cell <- results(dds, contrast=c("cell", "N061011", "N61311")) )
```
If results for an interaction term are desired, the `name`
argument of *results* should be used. Please see the
help for the *results* function for more details.
## Plotting counts for individual genes
A quick way to visualize the counts for a particular gene is to use
the *plotCounts* function, which takes as arguments the
*DESeqDataSet*, a gene name, and the group over which to plot the
counts.
```{r plotcounts, fig.width=5, fig.height=5}
topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene=topGene, intgroup=c("dex"))
```
We can also make more customizable plots using the *ggplot* function from the
`r CRANpkg("ggplot2")` package:
```{r ggplotcountsjitter, fig.height=5}
data <- plotCounts(dds, gene=topGene, intgroup=c("dex","cell"), returnData=TRUE)
ggplot(data, aes(x=dex, y=count, color=cell)) +
scale_y_log10() +
geom_point(position=position_jitter(width=.1,height=0), size=3)
```
Here we use a more structural arrangement instead of random jitter,
and color by the treatment.
```{r ggplotcountsdot, fig.height=5}
ggplot(data, aes(x=dex, y=count, fill=dex)) +
scale_y_log10() +
geom_dotplot(binaxis="y", stackdir="center")
```
Note that the *DESeq* test actually takes into account the cell line
effect, so a more detailed plot would also show the cell lines.
```{r ggplotcountsgroup, fig.height=5}
ggplot(data, aes(x=dex, y=count, color=cell, group=cell)) +
scale_y_log10() + geom_point(size=3) + geom_line()
```
## MA-plot: fold changes for all genes
An "MA-plot" provides a useful overview for an experiment with a
two-group comparison. The log2 fold change for a particular
comparison is plotted on the y-axis and the average of the counts
normalized by size factor is shown on the x-axis ("M" for minus,
because a log ratio is equal to log minus log, and "A" for average).
```{r plotma}
plotMA(res, ylim=c(-5,5))
```
Each gene is represented with a dot. Genes with an adjusted $p$ value
below a threshold (here 0.1, the default) are shown in red. The
*DESeq2* package incorporates a prior on log2 fold changes, resulting
in moderated log2 fold changes from genes with low counts and highly
variable counts, as can be seen by the narrowing of spread of points
on the left side of the plot. This plot demonstrates that only genes
with a large average normalized count contain sufficient information
to yield a significant call.
We can label individual points on the MA plot as well. Here we use the
*with* R function to plot a circle and text for a selected row of the
results object. Within the *with* function, only the `baseMean` and
`log2FoldChange` values for the selected rows of `res` are used.
```{r plotma2}
plotMA(res, ylim=c(-5,5))
with(res[topGene, ], {
points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
```
## Gene clustering
In the sample distance heatmap made previously, the dendrogram at the
side shows us a hierarchical clustering of the samples. Such a
clustering can also be performed for the genes. Since the clustering
is only relevant for genes that actually carry signal, one usually
carries it out only for a subset of most highly variable genes. Here,
for demonstration, let us select the 35 genes with the highest
variance across samples. We will work with the *rlog* transformed
counts:
```{r message=FALSE}
library("genefilter")
topVarGenes <- head(order(-rowVars(assay(rld))),20)
```
The heatmap becomes more interesting if we do not look at absolute
expression strength but rather at the amount by which each gene
deviates in a specific sample from the gene's average across all
samples. Hence, we center each genes' values across samples,
and plot a heatmap. We provide the column side colors to help identify
the treated samples (in blue) from the untreated samples (in grey).
```{r genescluster}
mat <- assay(rld)[ topVarGenes, ]
mat <- mat - rowMeans(mat)
df <- as.data.frame(colData(rld)[,c("cell","dex")])
pheatmap(mat, annotation_col=df)
```
We can now see blocks of genes which covary across patients. Note that
a set of genes at the top of the heatmap are separating the N061011
cell line from the others. At the bottom of the heatmap, we see a set
of genes for which the treated samples have higher gene expression.
## Annotation of results tables
Our result table only uses Ensembl gene IDs, but gene names may be
more informative. Bioconductor's annotation packages help with mapping
various ID schemes to each other.
We load the `r Biocpkg("AnnotationDbi")` package and the annotation package
`r Biocannopkg("org.Hs.eg.db")`:
```{r message=FALSE}
library("AnnotationDbi")
library("org.Hs.eg.db")
```
This is the organism annotation package ("org") for
*Homo sapiens* ("Hs"), organized as an *AnnotationDbi*
database package ("db"), using Entrez Gene IDs ("eg") as primary key.
To get a list of all available key types, use:
```{r}
columns(org.Hs.eg.db)
```
We can use the *mapIds* function to add invidual columns to our results
table. To add the gene symbol and Entrez ID, we call *mapIds* twice:
```{r}
res$symbol <- mapIds(org.Hs.eg.db,
keys=row.names(res),
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
keys=row.names(res),
column="ENTREZID",
keytype="ENSEMBL",
multiVals="first")
```
Now the results have the desired external gene ids:
```{r}
resOrdered <- res[order(res$padj),]
head(resOrdered)
```
## Exporting plain text results
You can easily save the results table in a CSV file, which you can
then load with a spreadsheet program such as Excel. The call to
*as.data.frame* is necessary to convert the *DataFrame* object
(`r Biocpkg("IRanges")` package) to a *data.frame* object which can be
processed by *write.csv*.
```{r eval=FALSE}
write.csv(as.data.frame(resOrdered), file="results.csv")
```
## Generation of HTML reports
Another more sophisticated package for exporting results from various
Bioconductor analysis packages is the `r Biocpkg("ReportingTools")` package.
*ReportingTools* will automatically generate dynamic HTML documents,
including links to external databases using gene identifiers
and boxplots summarizing the normalized counts across groups.
```{r message=FALSE}
library("ReportingTools")
```
The most basic report is just a nicer HTML interface to the results table in R. We convert our results object into a dataframe and pass this to the ReportingTools functions. Here we will take just the top 1000 genes:
```{r}
resTop <- head(resOrdered,1000)
resTop <- as(resTop, "data.frame")
resTop <- resTop[,c("padj","baseMean","log2FoldChange","pvalue",
"symbol","entrez")]
resTop$ensembl <- rownames(resTop)
```
There are three steps: (1) initializing a report, (2) publishing (pushing) the results to the report, and (3) finishing the report. The report will be written in a directory `reports` created in the current working directory.
```{r eval=FALSE}
getwd()
htmlRep <- HTMLReport(shortName = "airway_DE_report",
title = "Airway differential expression analysis 06/2015",
reportDirectory = "./reports")
publish(resTop, htmlRep)
url <- finish(htmlRep)
browseURL(url)
```
We can now open in a web browser the file `./reports/airway_DE_report.html`.
This gives a readable report with pagination, as well as sortable,
filterable numeric columns, and a search function.
We can also use some of the built-in functionality of *ReportingTools*
for *DESeq2* objects (there are similar functions for *edgeR* and *limma*
objects).
By including the *DESeqDataSet* along with a `factor` indicating the
condition of the samples, the output will include boxplots of normalized
counts. Remember, though that our statistical test was not just of
the differences across dexamethasone treatment, but the differences
of the pairs (taking account of the cell line information).
Additionally, if we change the primary ID of our dataset to be
Entrez IDs, then ReportingTools will also include the gene
symbol, name and a link to an external database. In order to see
what this looks like, let's create a subset of the results table
and *DESeqDataSet* in the same order: the top 50 genes by adjusted p-value.
We then add the Entrez ID as the rownames of both objects.
```{r}
resReport <- head(resOrdered, 50)
ddsReport <- dds[ rownames(resReport), ]
rownames(resReport) <- resReport$entrez
rownames(ddsReport) <- rownames(resReport)
```
Now we can pass these objects to ReportingTools:
```{r eval=FALSE}
htmlRep <- HTMLReport(shortName = "airway_DE_report",
title = "Airway differential expression analysis 06/2015",
reportDirectory = "./reports")
publish(resReport, htmlRep, DataSet=ddsReport,
annotation.db="org.Hs.eg.db",
n=50, factor=dds$dex)
url <- finish(htmlRep)
browseURL(url)
```
## Examining RNA-seq alignments
Here we will show some basic techniques for examining RNA-seq
alignments in R, using the `r Biocpkg("GenomicAlignments")`
and the `r Biocpkg("Gviz")` packages. First we will define a
*GRange* which covers the subset of genes we loaded in the
beginning of this workflow. In preparation for loading the
alignments into R, we also need to create an index file for the
BAM file we will be reading alignments from.
```{r}
total.range <- GRanges("1", IRanges(11e6, 11.4e6))
indexBam(filenames[1])
```
Now we can read in the alignments from the BAM file, and
we can specify that we only want to read in the alignments
from a specific region using `ScamBamParam`.
```{r}
ga <- readGAlignmentPairs(filenames[1], param=ScanBamParam(which=total.range))
ga
```
Note that without specifying a region, or without setting `yieldSize`
while defining a `BamFile`, this would start to load all the
alignments into your R session, which would likely take up too much
memory.
There are now a couple of methods defined for *GAlignmentPairs* objects,
which can be found on the help page `?GAlignmentPairs`. One of these is
to compute the genomic coverage. We can then look at the coverage over
a small area, convert this to a numeric vector and plot it using base R
graphics.
```{r simplecov}
cov <- coverage(ga)
zoom <- GRanges("1", IRanges(11.07e6, 11.09e6))
cov[zoom]
cov.numeric <- as.numeric(cov[zoom][[1]])
plot(cov.numeric, type="h")
```
Now we will use the `r Biocpkg("Gviz")` package to make a more sophisticated
version of this plot. See the vignette and help files for *Gviz* for more
details on each function and its arguments.
```{r gviz, message=FALSE}
library("Gviz")
options(ucscChromosomeNames=FALSE)
grtrack <- GeneRegionTrack(txdb)
gtrack <- GenomeAxisTrack()
altrack <- AlignmentsTrack(filenames[1])
plotTracks(list(gtrack, grtrack, altrack),
from=11.07e6, to=11.09e6, chromosome="1")
```
The above plot includes diagrams for each alignment. We can simplify
this visualization by just showing the coverage and the exon junction
coverage, similar to the [sashimi plot](https://www.broadinstitute.org/igv/Sashimi)
which can be viewed in IGV.
```{r gvizSashimi}
altrack <- AlignmentsTrack(filenames[1], type=c("coverage","sashimi"))
plotTracks(list(gtrack, grtrack, altrack),
from=11.07e6, to=11.09e6, chromosome="1")
```
Suppose now we want to know, of the transcripts which can be
viewed in this window, for which are most of the RNA-seq fragments
compatible? Here, *compatible* means that the junctions in
the aligned reads pairs must be compatible with the exon-exon junctions
of the transcript.
```{r}
ebt <- exonsBy(txdb, by="tx")
ebt.sub <- ebt[ebt %over% zoom]
fco <- findCompatibleOverlaps(ga, ebt.sub)
fco
```
We can investigate this *Hits* object, and tally the
number of fragments comptabible with each transcript.
We then save the top transcript, by fragment count.
```{r}
(tab <- table(subjectHits(fco)))
ebt.top <- ebt.sub[names(tab)[which.max(tab)]]
```
Now we can use *GeneRegionTrack* again, this time
providing a *GRangesList* with the single top
transcript. The *Gviz* package has many options
for visualization and different data input.
```{r gvizebttop}
grtrack2 <- GeneRegionTrack(ebt.top)
plotTracks(list(gtrack, grtrack2, altrack),
from=11.07e6, to=11.09e6, chromosome="1")
```
## Session information
As last part of this document, we call the function *sessionInfo*,
which reports the version numbers of R and all the packages used in
this session. It is good practice to always keep such a record as it
will help to trace down what has happened in case that an R script
ceases to work because the functions have been changed in a newer
version of a package. The session information should also **always**
be included in any emails to the
[Bioconductor support site](https://support.bioconductor.org) along
with all code used in the analysis.
```{r}
sessionInfo()
```