Using Limma to find differentially expressed genes

This is a lama, and I am learning limma. (Photo Credit: Brett Sayles)

Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies.Nucleic Acids Research 43(7), e47.

limma is an R package hosted on Bioconductor which finds differentially expressed genes for RNA-seq or microarray. Recently I’ve been working on a PCR-based low-density array and noticed that I forgot how to use limma for the one hundredth time, so I decided to make a note.

Material

• Log-transformed expression data in a matrix:
  Each column represents an experiement, and each row represents a detected gene/probe.
• Design matrix:
  Each column represents a status of your data (e.g., wild-type, mutant, rescued…), and each row corresponds to an experiment in the expression matrix. This could be generated from an expressionDataset object, or we could create this manually if you only have a raw expression dataset like me.

# Create a design matrix
design <- model.matrix(~ 0+factor(c(1,1,2,2,3,3)))
# assign the column names
colnames(design) <- c("Wild_type", "Mutant", "Rescue")

The code above should generate something like this:

• Contrast matrix: Then, we must tell limma whom we are going to compare with whom.

cont_matrix <- makeContrasts(MutvsWT = Mutant-Wild_type, MutvsRes = Mutant-Rescue, levels=design)

The code above would give:

Analysis

Now, I am pretty much ready (except for the need to understand the math).

# Fit the expression matrix to a linear model
fit <- lmFit(exp_matrix, design)
# Compute contrast
fit_contrast <- contrasts.fit(fit, cont_matrix)
# Bayes statistics of differential expression
# *There are several options to tweak!*
fit_contrast <- eBayes(fit_contrast)
# Generate a vocalno plot to visualize differential expression
volcanoplot(fit_contrast)
# Generate a list of top 100 differentially expressed genes
top_genes <- topTable(fit_contrast, number = 100, adjust = "BH")
# Summary of results (number of differentially expressed genes)
result <- decideTests(fit_contrast)
summary(result)

Have fun with the result, and enjoy some experiment time.