0.4 Solutions
0.4.1 Answer 1 : Load the Data
The data stored inside airway is a SummarizedExperiment object.
library("airway")
data(airway)
se <- airway
se## class: RangedSummarizedExperiment
## dim: 64102 8
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowRanges metadata column names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample0.4.2 Answer 2 : Create the DESeqDataSet
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.
library("DESeq2")
dds <- DESeqDataSet(se, design = ~ cell + dex)0.4.3 Answer 3 : Differential Expression Pipeline
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:
dds$dex <- relevel(dds$dex, "untrt")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:
dds <- DESeq(dds)## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testingThis 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.
0.4.4 Answer 4 : Build 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.
(res <- results(dds))## log2 fold change (MAP): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 64102 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003 708.60217 -0.37424998 0.09873107 -3.7906000 0.0001502838 0.001251416
## ENSG00000000005 0.00000 NA NA NA NA NA
## ENSG00000000419 520.29790 0.20215551 0.10929899 1.8495642 0.0643763851 0.192284345
## ENSG00000000457 237.16304 0.03624826 0.13684258 0.2648902 0.7910940556 0.910776144
## ENSG00000000460 57.93263 -0.08523371 0.24654402 -0.3457140 0.7295576905 0.878646940
## ... ... ... ... ... ... ...
## LRG_94 0 NA NA NA NA NA
## LRG_96 0 NA NA NA NA NA
## LRG_97 0 NA NA NA NA NA
## LRG_98 0 NA NA NA NA NA
## LRG_99 0 NA NA NA NA NAAs res is a DataFrame object, it carries metadata with information on the meaning of the columns:
mcols(res, use.names=TRUE)## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized counts for all samples
## log2FoldChange results log2 fold change (MAP): dex trt vs untrt
## lfcSE results standard error: dex trt vs untrt
## stat results Wald statistic: dex trt vs untrt
## pvalue results Wald test p-value: dex trt vs untrt
## padj results BH adjusted p-valuesThe 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} \approx 2.82\).
We can also summarize the results with the following line of code, which reports some additional information
summary(res)##
## out of 33469 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 2641, 7.9%
## LFC < 0 (down) : 2242, 6.7%
## outliers [1] : 0, 0%
## low counts [2] : 15441, 46%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results0.4.5 Answer 5 : Visualize Results
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.
topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene=topGene, intgroup=c("dex"))