RNA-seq differential expression analysis in R

Using edgeR and limma

Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law Based on the course RNAseq analysis in R delivered on May 11/12th 2016

Resources and data files

This material has been created using the following resources:
http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf [@Lun2016]
http://monashbioinformaticsplatform.github.io/RNAseq-DE-analysis-with-R/99-RNAseq_DE_analysis_with_R.html

Data files downloaded from:
ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60450/suppl/GSE60450_Lactation-GenewiseCounts.txt.gz http://bioinf.wehi.edu.au/software/MSigDB/mouse_c2_v5.rdata http://bioinf.wehi.edu.au/software/MSigDB/mouse_H_v5.rdata

Data files:
sampleinfo.txt
GSE60450_Lactation-GenewiseCounts.txt
mouse_c2_v5.rdata
mouse_H_v5.rdata

Data files available from: https://figshare.com/s/1d788fd384d33e913a2a You should download these files and place them in your /data directory.

Differential expression with edgeR

Now that we are happy that we have normalised the data and that the quality looks good, we can continue to testing for differentially expressed genes. There are a number of packages to analyse RNA-Seq data. Most people use DESEQ2 or edgeR. We will use edgeR for the rest of this practical.

**First make sure we have all the objects and libraries loaded*

library(edgeR)
library(limma)
library(Glimma)
library(gplots)
library(org.Mm.eg.db)
load("Robjects/preprocessing.Rdata")

Recap of pre-processing

The previous section walked-through the pre-processing and transformation of the count data. Here, for completeness, we list the minimal steps required to process the data prior to differential expression analysis.

## Read the counts from the downloaded data
seqdata <- read.delim("data/GSE60450_Lactation-GenewiseCounts.txt", stringsAsFactors = FALSE)
#
# Remove first two columns from seqdata

countdata <- seqdata[,-(1:2)]

# Store EntrezGeneID as rownames
rownames(countdata) <- seqdata[,1]
countdata
colnames(countdata) <- substr(colnames(countdata), 1, 7)
countdata
## Calculate the Counts Per Million measure
myCPM <- cpm(countdata)
## Identify genes with at least 0.5 cpm in at least 2 samples
thresh <- myCPM > 0.5
keep <- rowSums(thresh) >= 2
# Subset the rows of countdata to keep the more highly expressed genes
counts.keep <- countdata[keep,]
## Convert to an edgeR object
dgeObj <- DGEList(counts.keep)
## Perform TMM normalisation
dgeObj <- calcNormFactors(dgeObj)
## Obtain corrected sample information
sampleinfo <- read.delim("data/SampleInfo_Corrected.txt")
group <- paste(sampleinfo$CellType,sampleinfo$Status,sep=".")
group

Create the design matrix

First we need to create a design matrix for the groups, as we have seen in the linear models lecture. We have two variables, status and cell type. We will fit two models under two assumptions; no interaction and interaction of these two factors.

Let’s start with the model with only main effects, that is no interaction. The main assumption here is that the effect of the status is the same in both type of cells.

# Create the two variables
group <- as.character(group)
type <- sapply(strsplit(group, ".", fixed=T), function(x) x[1])
status <- sapply(strsplit(group, ".", fixed=T), function(x) x[2])
# Specify a design matrix with an intercept term
design <- model.matrix(~ type + status)
design

Data exploration

An MDS plot shows distances, in terms of biological coefficient of variation (BCV), between samples. What do you think of the quality of the data? Can you anticipate if the interaction term will be important?

plotMDS(dgeObj, labels=group, cex=0.75, xlim=c(-4, 5))

Estimating the dispersion

The common dispersion estimates the overall BCV of the dataset, averaged over all genes:

dgeObj <- estimateCommonDisp(dgeObj)

Then we estimate gene-wise dispersion estimates, allowing a possible trend with averge count size:

dgeObj <- estimateGLMTrendedDisp(dgeObj)
dgeObj <- estimateTagwiseDisp(dgeObj)

Plot the estimated dispersions:

plotBCV(dgeObj)

Testing for differential expression

First, we fit genewise glms:

# Fit the linear model
fit <- glmFit(dgeObj, design)
names(fit)
head(coef(fit))

Conduct likelihood ratio tests for luminal vs basal and show the top genes:

lrt.BvsL <- glmLRT(fit, coef=2)
topTags(lrt.BvsL)

Challenge

Conduct likelihood ratio tests for virgin vs lactate and show the top genes.

Contrasts

Suppose we want to find differentially expressed genes between pregnant and virgin. We don’t have a parameter that explicitly will allow us to test that hypothesis. We need to build a contrast:

PvsV <- makeContrasts(statuspregnant-statusvirgin, levels=design)
lrt.pVsV <- glmLRT(fit, contrast=PvsV)
topTags(lrt.pVsV)

Challenge

1.Fit a model with interaction: What is the rationale to include the interaction (What assumption are you relaxing?) 2. Is the number of replicates good enough to include the interaction? 3. Is the interaction needed in the model?

Solution

Testing relative to a threshold

Statistical tests can sometimes find small log-fold changes to be highly significant. In many biological applications, small fold-change differences may not be interesting, in which case we would like to find “significant” DE genes that also have a minimum fold-change. The glmTreat function to test for differential expression relative to a fold-change threshold that we apply.

In this example, we will say that we are not interesting in any gene that has a fold-change less than 1.5, so we set the threshold to log2(1.5).

trt.BvsL <- glmTreat(fit, coef = 2, lfc = log2(1.5))
topTags(trt.BvsL)

Here, the top DE genes have large fold-changes anyway, but observe the difference in p-values between this TREAT test and the standard LRT above.

Challenge

1. Experiment with different lfc thresholds in glmTreat. How do these changes affect the results and any conclusions you might draw?

Solution

save(lrt.BvsL, dgeObj, group,file="Robjects/DE.Rdata")

Alternative: limma analysis

If the sequencing depth is reasonably consistent across the RNA samples, then the simplest and most robust approach to differential exis to use limma-trend. This approach will usually work well if the ratio of the largest library size to the smallest is not more than about 3-fold.

limma-trend

In the limma-trend approach, the counts are converted to logCPM values using edgeR’s cpm function:

logCPM <- cpm(dgeObj, log = TRUE, prior.count = 3)

The prior count is used here to damp down the variances of logarithms of low counts. The logCPM values can then be used in any standard limma pipeline, using the trend=TRUE argument when running eBayes. For example:

lmfit <- lmFit(logCPM, design)
lmfit <- eBayes(lmfit, trend=TRUE)
topTable(lmfit, coef = 2)

How does this gene ranking compare to that obtained with edgeR?

voom

When the library sizes are quite variable between samples, then the voom approach is theoretically more powerful than limma-trend. In this approach, the voom transformation is applied to the normalized and filtered DGEList object:

v <- voom(dgeObj, design, plot = TRUE)

After this, the usual limma pipelines for differential expression can be applied, for example:

vfit <- lmFit(v, design)
vfit <- eBayes(vfit)
topTable(vfit, coef = 2)

Or, to give more weight to fold-changes in the ranking, one could use a treat approach as introduced earlier:

tvfit <- treat(vfit, lfc = log2(1.5))
topTreat(tvfit, coef = 2)

Challenge

1. Compare the results obtained with edgeR, treat, limma-trend and limma-voom. How do they differ? Are these differences important?

Solution