In jtlovell/limmaDE2: A pipeline for differential expression analyses in LIMMA

knitr::opts_chunk$set(fig.width=8, fig.height=6,echo=TRUE, warning=FALSE, message=FALSE)

Overview

The R package limmaDE2 contains a set of functions designed to make running differential expression analyses in LIMMA more user friendly. There are several functions of this package:

• pipeLIMMA Run analyses of differnetial expression in LIMMA.
• makeBinarySig Find and count the number of significantly differentially expressed genes
• pqHists Plot the distribution of P and FDR-corrected P (Q) values
• voom2PCA Conduct principal component analyses on voom (or otherwise) normalized expression matrices
• volcanoPlot Generate log2 foldchange vs. P-value volcano plots
• volcanoPair Compare two sets of log2 fold changes, colored by p-values
• pipeDESeq2 Runs an analysis of differential expression similar to that of pipeLIMMA, except through the DESeq2 package.

1. Installation

The limmaDE2 packages is not on CRAN and can only be installed from github. Make sure the devtools package is installed on your system. The limmaDE2 package also requires several other R packages to be installed. These packages have functions upon which limmaDE2 depends:

Essential Packages:

• limma contains all of the basic statistical elements that are called here.
• qvalue fdr-corrected p-values
• edgeR Differential expression package used for some transformations

Data for this vignette comes from:

• SimSeq permits simulation of RNA-seq data

Other useful packages:

• plyr summarize and manipulate lists and dataframes
• vegan some distance-based methods
• Heatplus heatmap plotting
• topGO gene ontology (GO) analysis
• ggplot2 routines for complex plots
• hexbin methods to simplify data in plots

To install limmaDE2:

library(devtools)
install_github("jtlovell/limmaDE2")

Load all packages that will be needed for the tutorial:

library("limmaDE2")
library("ggplot2")
library("reshape2")
library("SimSeq")
library("DESeq2")

2. Importing data and checking

For limmaDE2 to work, two matrices are needed:

• counts: Raw transcript abundance counts
• info: Experimental design information These two matrices must match exactly, where each row in the info matrix corresponds to each column in counts. For best performance, the name of each gene should be stored in the rownames of the counts matrix.

data(kidney) # from simseq
counts<-kidney$counts
counts<-counts[sample(1:nrow(counts),1000),]
info<-with(kidney, 
           data.frame(id = paste(replic, treatment, sep = "_"), 
                      rep=replic, 
                      Treatment=ifelse(treatment == "Tumor","tumor","cntr"), 
                      stringsAsFactors=F))
colnames(counts)<-info$id

Add in another category that gives us a factorial design

group.a<-c("4619", "4712", "4863", "4865", "5452", "5453", "5454", "5455",
          "5456", "5457","5458", "5461", "5462", "5463", "5465", "5467",
          "5468", "5469", "5470", "5549","5552", "5580", "5641", "5672",
          "5689", "5701", "5703", "5706", "5989", "6088")
info$group<-ifelse(info$rep %in% group.a, "a","b")
with(info, table(group, Treatment))

info$trt.grp<-with(info, paste(Treatment, group, sep="_"))

The top of the example experimental design dataset

head(info)

What the example counts dataset looks like

counts[1:10,1:3]

In all statistical analyses, it is important to set the levels of the experimental factors. This is esspecially true in linear modelling, such as limma, where levels are tested against the base level. Here, we set the "Non-Tumor" treatment and "a" group as the base.

info$Treatment <- factor(info$Treatment,
                           levels = c("cntr", "tumor"))
info$group <- factor(info$group,
                           levels = c("a", "b"))

3. Basic analysis of differential expression

To test the plasticity (differential expression) of all genes between the tumor and control tissue:

stats <- pipeLIMMA(counts = counts, 
                   info = info, 
                   block = NULL, 
                   formula = "~ Treatment")
lmStats<-stats$stats
voom<-stats$voom$E

pipeLIMMA returns three elements: stats, voom and fstats. These are the statistical output of the limma functions: eBayes, voom and topTableF. Inspect the top few rows and columns of each.

stats: eBayes linear model statistics

knitr::kable(stats$stats[1:6,1:6])

voom: normalized expression matrix

knitr::kable(stats$voom[["E"]][1:6,1:6])

fstats: F-statistics for each factor in the model, in this case, just treatment

knitr::kable(stats$fstats[1:6,])

4. Differential expression using a blocking variable

Lets say that here, there is some sort of experimental blocking. In this case, we employ a routine in limma that calculates the duplicateCorrelation among blocks, then uses the blocking variable in the linear model fit.

info$block <- rep(1:2,each=nrow(info)/2)
stats.block <- pipeLIMMA(counts = counts, 
                   info = info, 
                   block = info$block,
                   formula = "~ Treatment")

5. Differential expression in a factorial experiment

stats.factorial <- pipeLIMMA(counts = counts, 
                   info = info, 
                   block = NULL, 
                   formula = "~ Treatment + group + Treatment*group")

6. Calculating differential expression using specific contrasts

Sometimes, it may make more sense to use specific contrasts to test for differential expression. For example, let's say we are interested in how the tumors and non-tumor tissues differ for each of the two groups.

First, we need to construct a design matrix.

design <- with(info, model.matrix(~ 0 + trt.grp))
colnames(design)<-gsub("trt.grp","",colnames(design))
head(design)

Then we construct a contrast matrix.

contrast.matrix <- makeContrasts(
  tumor_a - cntr_a , 
  tumor_b - cntr_b,
  levels = design)
head(contrast.matrix)

Finally, we fit the model with a contrast matrix as the design, overriding the formula argument.

stats <- pipeLIMMA(counts = counts, 
                   info = info, 
                   block = NULL, 
                   design = design, 
                   contrast.matrix = contrast.matrix)
stats.contrasts <- stats$stats

6. Count and plot the number of significantly differentially expressed genes

The function makeBinarySig looks for a provided string (e.g. "Q.Value") and outputs a matrix with whether or not those columns have values <= the provided alpha

sigs <- makeBinarySig(stats.contrasts, what = "Qvalue", alpha = 0.05)

Make a venn diagram of number of differentially expressed genes among the experimental factors

counts2Venn(x=sigs, cols=c(1,2), names=c("in.grpA","in.grpB"),
   colors=c("blue","red"),type="limma", legx=-3.3,legy=-3)

Make a euler diagram of number of differentially expressed genes among the experimental factors

counts2Venn(x=sigs, cols=c(1,2), names=c("in.grpA","in.grpB"),
   colors=c("blue","red"),type="Euler", legx=-3.3,legy=-3)

7. Make a volcano plot of the results

Volcano plots are a good way to look at the differences between two experimental levels. Here, we compare the extent of differential expression between the "high" treatment to the "low" treatment.

with(lmStats, volcanoPlot(pval = ebayesPvalue_Treatmenttumor,
                          lfc = Treatmenttumor_logFC,
                          sig = ebayesQvalue_Treatmenttumor,
                          main = "no tumor vs. tumor Volcano Plot", 
                          xlab = "tumor - no tumor Log2 Fold Change",
                          bty = "n", legpos = "top", leginset = c(0,-.1)))

8. Using two contrasts, make a pairwise volcano plot

Custom colors

sigs <- data.frame(makeBinarySig(stats.contrasts, what = "Qvalue", alpha = 0.05))
names(sigs)<-c("sig.a","sig.b")
sigs$sign.a<-sign(stats.contrasts$tumor_a...cntr_a_logFC)
sigs$sign.b<-sign(stats.contrasts$tumor_b...cntr_b_logFC)

cols<-with(sigs, ifelse(sig.a + sig.b == 0,  
                        "grey",
                        ifelse(sig.a + sig.b == 2 & sign.a*sign.b == 1, 
                               "pink",
                               ifelse(sig.a + sig.b == 2 & sign.a*sign.b == -1,
                                      "cornflowerblue",
                                      ifelse(sig.a == 1, "darkblue", "darkred")))))

Make the plot

with(stats.contrasts, volcanoPair(lfc1 = tumor_a...cntr_a_logFC,
                                  lfc2 = tumor_b...cntr_b_logFC,
                                  pt.col = cols, pt.pch = 19, pt.cex=.5,
                                  xlab = "Tumor - control LFC (group A)",
                                  ylab = "Tumor - control LFC (group B)"))

Principal component analysis of gene expression

It is often a good idea to get an idea of how the individuals (libraries) are structured - how similar are they to eachother. Principal component analyses allow for this kind of inference.

pca12 <- counts2PCA(counts=voom, info = info, ids = info$id, pcas2return = 6)
pca12.var <- pca12[[2]]
pca12 <- pca12[[1]]
gcols <- c("darkblue", "blue", "gold", "green", "pink", "red")
ggplot(pca12, aes(x = PC1, y = PC2, col = group, shape = Treatment)) +
  geom_vline(xintercept = 0, lty = 2)+
  geom_hline(yintercept = 0, lty = 2)+
  geom_point() +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank()) +
  labs(x = paste("PCA1 (", pca12.var[1],"%)", sep = ""),
       y = paste("PCA2 (", pca12.var[2],"%)", sep = ""),
       title = "PCA analysis of voom-normalized expression")

For more information, visit lovelleeb.weebly.com/analytics

jtlovell/limmaDE2 documentation built on May 20, 2019, 3:14 a.m.