RNASeq in parathyroid adenomas

This work flow derives from a more comprehensive example developed in
the parathyroidSE package and the
CSAMA 2014
workshop. It uses DESeq2 to estimate gene-level differential
expression. Analysis picks up after counting reads aligning to
Ensembl gene regions; counting can be accomplished using
summarizeOverlaps(), as described in the parathyroidSE vignette.

The data is from
Haglund et al.,
Evidence of a Functional Estrogen Receptor in Parathyroid Adenomas, J
Clin Endocrin Metab, Sep 2012. The experiment investigated the role
of the estrogen receptor in parathyroid tumors. The investigators
derived primary cultures of parathyroid adenoma cells from 4
patients. These primary cultures were treated with diarylpropionitrile
(DPN), an estrogen receptor β agonist, or with 4-hydroxytamoxifen
(OHT). RNA was extracted at 24 hours and 48 hours from cultures under
treatment and control.

Start by loading the count data, which has been carefully curatd as a
SummarizedExperiment object. Explore the object with accessors such
as rowData(se), colData(se), exptData(se)$MIAME, and
head(assay(se))

The first step is to add an experimental design and perform additional
preparatory steps. The experiment contains several technical
replicates, and while these could be incorporated in the model, here
we simply pool them (pooling technical replicates is often appropriate
for RNASeq data).

The entire work flow is executed with a call to the DESeq()
function. This includes library size factor estimation, dispersion
estimates, model fitting, independent filtering, and formulating test
statistics based on the specified design. Here we generate the results
for a comparison between DPN and control.

The remaining steps provide some basic checks and visualizations of
the data. We identify genes for which the adjusted (for multiple
comparison) p value is less than 0.1, and order these by their log2
fold change.

An 'MA' plot represents each gene as a point, with average expression
over all samples on the x-axis, and log2 fold change between treatment
groups on the y-axis; highlighted values are genes with adjusted p
values less than 0.1.