---
title: "doseR"
author: "AJ Vaestermark, JR Walters."
date: "`r Sys.Date()`"
output: BiocStyle::html_document
vignette: >
%\VignetteIndexEntry{"doseR"}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
# Table of Contents
1. Introduction
2. Overview
3. Inputs
4. Normalization
5. Filtering
6. Comparing absolute expression
7. Comparing relative expression
8. Example workflow
9. Installation
10. Normalization & RPKM
11. Filtering "unexpressed" loci
12. Comparing absolute expression
13. Comparing relative expression
14. Workflow conclusion
# Introduction
The doseR package aims to provide a consistent, efficient, and generalizable
analytical framework for assessing sex chromosome dosage compensation using
RNA-seq based read-count data. The study of sex chromosome dosage compensation
has garnered increasing interest in the last 15 years or so. This growth has
been driven in part by the rise of RNA-seq as a method for assaying genome-wide
patterns of gene expression, along with the discovery of a surprising diversity
of dosage compensation patterns across species. Yet despite the growing
interest in the topic of dosage compensation, with many similarities among
the analyses performed in various studies, each subsequent study tends to
employ its own bespoke scripting to analyze and visualize the data. There
has yet to emerge a common computational framework specifically developed
to streamline and increase reproducibility among studies of dosage
compensation. Our hope is for doseR to help fill that void.
Functions provided in doseR offer support for several analytical steps
commonly done in studies of dosage compensation, including:
* Normalizing data relative to library size
* Removing unexpressed or low-abundance transcripts
* Visualizing absolute and relative expression levels across gene sets
* Statistical comparisons between distributions of gene expressions
Many of these elements may be of interest for studies of gene expression
not specifically concerning dosage compensation, and the doseR package is
not meant to be limited only to studies contrasting sex chromosomes versus
autosomes or males versus females. doseR supports any number of comparisons
between biological treatments and groupings of loci.
# Overview
The standard goal of most sex chromosome dosage compensation studies is
to assess whether the sex chromosome (i.e. the X or Z chromosome, depending
on male or female heterogamety) shows a globally distinct pattern of gene
expression between the sexes. On average, autosomal expression is expected to
be comparable between the sexes because autosomes are diploid in both sexes. In
contrast, differences in gene dose between the sexes due to heterogamety may
influence average expression on the X or Z, depending on what compensating
mechanisms have evolved. Thus most contemporary studies of sex chromosome
dosage compensation proceed by using sex-specific RNA-seq data to compare
distributions of normalized gene expression across chromosomes.
## Inputs
We assume users are starting with read count data from an RNA-seq experiment,
as is typically generated after employing one of many commonly used alignment
pipelines (e.g. RSEM, Salmon, etc). The doseR workflow starts with read count
data, and not RPKM (or, comparably, FPKM), because RPKM values are calculated
during the analysis after selecting the desired library-size normalization
factor.
RNA-seq data will typically come from distinct treatments (e.g. male, female)
and may have one or more replicates per treatment. Users also provide lengths
(in basepairs) for each transcript assayed. Finally, annotations indicating
the groupings of loci (e.g. chromosome) are also needed.
## Normalization
Accounting for differences in library size (i.e. depth of sequencing) between
samples is an important step in most statistical analyses involving RNA-seq
data.
doseR provides a variety of methods for defining the library size for each
sample.
After defining library sizes, read counts are scaled both by library size
and transcript length to give RPKM values that can be meaningfully compared
between loci and between samples.
## Filtering
One notable concern and potential artifact in dosage compensation analysis
is the filtering of “unexpressed” loci. In theory, it is ideal to remove
loci prior to analysis since we should be primarily concerned only with
actively expressed genes. In practice, however, it is far from clear how
best to do this or what criteria should be employed. When “unexpressed”
loci are disproportionately represented among chromosomes (or other gene
groups), failing to effectively remove them can cause spurious results.
At the other extreme, overly aggressive removal of low-abundance transcripts
can artificially compress distributions of gene expression, masking true
differences in the data. A variety of approaches have been implement in doseR
for filtering out loci considered to be “unexpressed”.
## Comparing absolute expression
One informative way to assess dosage compensation is to consider the
distribution of expression levels within sex between groups of genes,
typically for genes on the X or Z versus the autosomes. This is done
separately in each sex (or other treatments), providing a direct
inference of absolute differences in gene expression among chromosomes
(or other gene groupings). However, this approach is limited in that
it does not easily support a direct comparison between sexes (or
other treatments). doseR includes functions for visualization and
statistical tests of differences in absolute gene expression.
## Comparing relative expression
Examining the distribution of expression ratios (e.g. Male:Female)
between treatments across chromosomes (or other gene groupings) is
an alternative approach to assess whether global expression biases
exist for particular sets of genes, such as might be expected for
uncompensated X or Z chromosomes. While this approach can provide
nuanced assessment of dosage effects, it does not indicate what
absolute changes in gene expression may be occurring that underlie
shifts in expression ratio. doseR provides functions for
visualizing and statistically testing variation in expression
ratios across chromosomes (or other gene groupings), complementing
functions for examining absolute differences in expression among
groupings within treatment.
----------------------------------------
# Example workflow
To demonstrate a typical workflow using doseR, we will use a subset
of data from a study in Heliconius butterflies. H. melpomene
has 20 autosomes with chromosome 21 being the "Z" sex chromosome.
Like all Lepidoptera, this is a female-heterogametic species, so males
are diploid for the Z while females have a single Z and W chromosomes.
The relevant data are included in the doseR package, which must be installed.
## Installation
doseR can be installed from Bioconductor directly using the `BiocManager`.
```{r, eval = F}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("doseR")
```
Once installed, you must load the package.
```{r}
library(doseR)
library(SummarizedExperiment)
```
## Example data: _Heliconius_ RNA-seq
The primary data used by doseR are read counts from RNA-seq data,
in this case generated from legs of H. melpomene, with three
biological replicates each for male and female. Read count data
should be formatted as a matrix, with replicates in columns and
transcripts in rows. We also need a vector of corresponding transcript
lengths. Finally, data indicating gene groupings are needed. Here we
have a vector giving chromosomal linkage for each transcript (note
that these data are limited to only one representative transcript per gene).
```{r}
data(hmel.data.doser)
```
These three data elements are used to populate a new `SummarizedExperiment`,
which is the core data structure for the doseR package. This object
contains and organizes the read counts, RPKM values, segment lengths,
and annotations of gene groupings. However, there are some further data
manipulations required before we can create a new object.
First, we need to create a vector indicating to the replicate structure
in the data, in this case three male and three female replicates, as
organized in the counts matrix.
```{r}
reps <- c("Male", "Male", "Male", "Female", "Female", "Female")
```
Second, we must create a dataframe to hold the annotation of gene
groupings. In doing this, we will also convert the numeric
annotation vector into a factor, which facilitates
proper ordering of gene groups in later plotting.
Additionally, we will recode the chromosomal linkage information as
being Z-linked or autosomal.
```{r}
annotxn <- data.frame("Chromosome" = factor(hmel.dat$chromosome,
levels = 1:21))
annotxn$ZA <- factor(ifelse(hmel.dat$chromosome == 21, "Z", "A"),
levels = c("A", "Z"))
```
Now we are ready to populate a new `SummarizedExperiment` object.
```{r}
counts <- hmel.dat$readcounts
colData <- S4Vectors::DataFrame(Treatment=as.factor(reps),
row.names=colnames(hmel.dat$readcounts))
rowData <- S4Vectors::DataFrame(annotation = annotxn,
seglens = hmel.dat$trxLength,
row.names=rownames(hmel.dat$readcounts) )
se2 <- SummarizedExperiment(assays=list(counts=counts),
colData=colData, rowData=rowData)
```
## Normalization & RPKM
Next we select a method for calculating the library size (sometimes
called scaling factors) so that data are normalized for depth of
sequencing. Various methods for calculating a library size are provided,
but we will use the most straightforward method of summing total read counts.
```{r}
SummarizedExperiment::colData(se2)$Libsizes <- getLibsizes3(se2,
estimationType = "total")
```
With library sizes now stored in the `SummExperiment` object, we can calculate
RPKM values for each locus in each sample and store them in the RPKM slot.
```{r}
SummarizedExperiment::assays(se2)$rpkm <- make_RPKM(se2)
# se2 now equals se...
data(hmel.se) # loads hmel.dat list
# MD5 checksum equivalence:
#library(digest)
digest::digest(se) == digest::digest(se2)
```
It is always wise when first working with comparisons between genome-wide
expression data sets like this to check that there are not clear biases or
other artifacts showing up in the data, and that the normalization seems
reasonable. Making [MA plots](https://en.wikipedia.org/wiki/MA_plot) is one
important visualization to guard agains such problems.
```{r, fig.align="center", fig.width=6}
plotMA.se(se, samplesA ="Male", samplesB = "Female", cex = .2 ,
pch = 19, col = rgb(0,0,0, .2), xlab = "Log2(Average RPKM)",
ylab = "Log2(Male:Female)")
```
## Filtering "unexpressed" loci
It is generally good practice to remove loci considered to be "unexpressed".
In this case we will apply a filter so that included loci must have a mean RPKM
value across samples greather than 0.01. Many other approaches to filtering
are provided in doseR, and it is generally wise to explore how filtering
choices affect assessments of dosage compensation.
```{r}
f_se <- simpleFilter(se, mean_cutoff = 0.01, counts = FALSE)
```
## Comparing absolute expression
Having generated filtered RPKM values, we can proceed to comparing
distributions of expression among different groupings of genes, as
indicated in the annotation slot of the `SummExperiment` object. Biological
replicates are averaged to give a single representative RPKM value per gene.
In this case it is clear that the gene expression is reduced on the Z
chromosome relative to autosomes in both sexes.
```{r, fig.align="center", fig.width=5}
plotExpr(f_se, groupings = "annotation.ZA", clusterby_grouping = FALSE,
col=c("grey80","red","grey80","red"), notch=TRUE, ylab = "Log2(RPKM)")
```
Beyond visualizing with boxplots, we will also want statistical
summaries and tests applied to these data, which can be done with
the `generateStats()` function. Here we will want to
subset the `SummarizedExperiment` object so that tests are applied
using only data from a single treatment.
The function returns a list of three items:
1) A distributional summary for each grouping of genes analyzed
2) The results from applying a Kruskal-Wallis test, a non-parametric
test of differences among groups
3) The values used for generating 1 & 2, i.e. the average values
across replicates for each gene in each group.
Here we demonstrate this only for male samples, but a full analysis
would also require the same be done for the female samples.
```{r}
se.male <- f_se[, colData(f_se)$Treatment == "Male"]
male_ZvA <- generateStats(se.male , groupings = "annotation.ZA", LOG2 = FALSE)
male_ZvA$summary # distributional summary statistics
male_ZvA$kruskal # htest class output from kruskal.test()
lapply(male_ZvA$data, head) # a record of values used for statistics.
```
It may be of interest to examine other groupings of genes. For instance,
how variable are the distributions of expression across chromosomes?
Again, we'll limit this only to males for the sake of demonstration.
```{r, fig.align="center"}
plotExpr(f_se, groupings = "annotation.Chromosome", col=c(rep("grey80", 20),
"red"), notch=TRUE, ylab = "Log2(RPKM)", las = 2, treatment = "Male",
clusterby_grouping = TRUE )
```
## Comparing relative expression
Similarly to the plots and statistics given above for absolute expression,
doseR provides functions for examing relative expression between pairs of
treatments. For instance, we may want to assess whether there is a dosage
effect on gene expression for the Z chromosome. This would manifest as a
positive shift in the distribution of Male:Female expression levels on the
Z relative to autosomes.
We can examine this using boxplots, similar to absolute expression, above.
```{r, fig.align="center", fig.width=4}
plotRatioBoxes(f_se, groupings = "annotation.ZA", treatment1 = "Male",
treatment2 = "Female", outline = FALSE, col = c("grey80", "red"),
ylab = "Log(Male:Female)" )
```
Alternatively, there is a widely used convention in the dosage compensation
literature to visualize these distributions using a density plot, so this
functionality is also provided.
```{r, fig.align="center", fig.width = 5}
plotRatioDensity(f_se, groupings = "annotation.ZA", treatment1 = "Male",
treatment2 = "Female", type = "l", xlab = "Log(Male:Female)", ylab = "Density")
```
As above, we can apply the same plotting to the chromosomal
annotation if desired.
```{r, fig.align="center", fig.width=10}
par(mfrow = c(1,2))
plotRatioBoxes(f_se, groupings = "annotation.Chromosome", treatment1 =
"Male", treatment2 = "Female", outline = FALSE, col=c(rep("grey80", 20),
"red"), ylab = "Log(Male:Female)", xlab = "Chromosome" )
plotRatioDensity(f_se, groupings = "annotation.Chromosome", treatment1 =
"Male", treatment2 = "Female", type = "l", xlab = "Log(Male:Female)",
ylab = "Density", col=c(rep("grey80", 20), "red"), lty = 1)
```
Beyond plotting, we also want to know whether there are statistical
differences between the groupings indicated. The `test_diffs()` function
summarizes the distributions of expression ratios among groups and applies
a non-parametric Kruskal-Wallis test.
```{r}
za.ratios.test <- test_diffs(f_se, groupings = "annotation.ZA",
treatment1 = "Male", treatment2 = "Female", LOG2 = FALSE )
za.ratios.test$summary # summary statistics for each grouping
za.ratios.test$kruskal # htest class output from kruskal.test()
lapply(za.ratios.test$data, head) # values used for summaries and tests
```
## Workflow conclusion
As can be seen from these plots and statistics, these data reveal
a distinct reduction in absolute expression for Z linked loci, relative
to autosomes, for both sexes. Furthermore, there is a modest but
statistically significant dosage effect causing a male-bias in the
male:female ratio expression on the Z chromosome. The functions
support various comparisons between treatments and among any set
of gene groupings, as specified in the annotations.
In some cases, when analyzing more than two gene groupings (e.g.
across individual chromosomes), users may want to extend statistical
analyses to specific pairwise comparisons. For instance, is the
average Z expression specifically different from chromosome 14 in
particular? Because the Kruskal-Wallis test is a "global" test for
inequality among groups (analogous to an ANOVA), you may want to apply
a follow-up non-parametric post-hoc test for specific pairwise comparisons
after a significant Kruskal-Wallis test. Several such tests, such as the
[Nemenyi](https://en.wikipedia.org/wiki/Nemenyi_test), are implemented
in the [PMCMRplus]
(https://cran.r-project.org/web/packages/PMCMRplus/index.html) package.
