Differential expression analysis for studies with more than 2 groups

Last updated: 2018-01-12

Code version: 20cdf3d

Analyze RNA-seq data from Hoban et al., 2016, which measured gene expression in the prefrontal cortex of 3 groups of mice: conventional (con), germ-free (gf), and germ-free colonized with normal microbiota (exgf).

library("Biobase")
library("edgeR")

Warning: package 'limma' was built under R version 3.3.3

library("ggplot2")
library("limma")
dge <- readRDS("../data/ch03.rds")

Studies with more than 2 groups (Video)

Describe the scientific question, the experimental design, and the data collected for the 3-group study. Review RNA-seq technique and the need to standardize by library size.

Calculate the library size

Use colSums to calculate the library sizes and run summary to appreciate the variability.

Calculate counts per million

Use edgeR::calcNormFactors and edgeR::cpm to calculate TMM-normalized counts per million.

Create density plots of counts versus counts per million

Use limma::plotDensities to visualize the distribution of counts and the distribution of counts per million.

Filtering genes and samples (Video)

Describe the process of filtering lowly expressed genes (show density plots from previous exercise to demonstrate bimodal distribution) and identifying outlier samples with PCA. Discuss the nuances of removing outlier samples.

Remove lowly expressed genes

Remove any gene that does not have at least a log2 cpm > 0 in at least 4 samples (n = 4 per group). Use limma::plotDensities to confirm the distribution of counts per million is unimodal post-filtering.

Use PCA to identify outlier sample

Use prcomp to calculate PCA and then plot PC1 vs. PC2. Remove the exgf sample that clusters with the gf samples.

Differential expression for RNA-seq with voom (Video)

Review how to use model.matrix and limma::makeContrasts. Describe how limma::voom corrects for the mean-variance relationship of the count data.

Create design matrix for study with 3 groups

Use model.matrix to create a linear model with three binary variables (group-means parametrization).

Create contrasts matrix for study with 3 groups

Use limma::makeContrasts to test all pairwise comparisions.

Account for mean-variance relationship with voom

Use limma::voom to calculate weights for the linear model to account for the mean-variance relationship of the counts. Set plot = TRUE to visualize the mean-variance trend. Explore the EList object returned to confirm that the expression values in E are simply the counts per million.

Fit and test model for study with 3 groups

Use limma::lmFit, limma::contrasts.fit, limma::eBayes, and limma::decideTests to fit and test the model. Use limma::vennDiagram to visualize overlap of differentially expressed genes.

Session information

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X Yosemite 10.10.5

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] ggplot2_2.2.1       edgeR_3.16.5        limma_3.30.13      
[4] Biobase_2.34.0      BiocGenerics_0.20.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.11     knitr_1.18       magrittr_1.5     munsell_0.4.3   
 [5] colorspace_1.3-2 lattice_0.20-34  rlang_0.1.6      stringr_1.2.0   
 [9] plyr_1.8.4       tools_3.3.2      grid_3.3.2       gtable_0.2.0    
[13] git2r_0.20.0     htmltools_0.3.6  lazyeval_0.2.1   yaml_2.1.16     
[17] rprojroot_1.2    digest_0.6.13    tibble_1.3.3     evaluate_0.10.1 
[21] rmarkdown_1.8    stringi_1.1.5    scales_0.4.1     backports_1.1.0 
[25] locfit_1.5-9.1

This R Markdown site was created with workflowr