Suggested solution - Introduction to edgeR GLMs

Here we demonstrate the use of edgeR to perform a differential expression analysis using data from Tuch et al. (PLOS) as detailed in the edgeR manual.

The data is a set of RNA-seq samples of oral squamous cell carcinomas and matched normal tissue from three patients that were previously quantified into raw counts.

We will use edgeR to do a differential expression analysis of Tumor vs Non-Tumor samples. We will start with a simple pairwise comparison of the Tumor and Non-Tumor samples, and then repeat the analysis adding the patient pairing information to the model design.


Load the count data

We start by importing the counts table into R using the read.delim function. Other functions to import tables include read.table and read.csv. We also specify that the values in the tables are separated by a TAB. You can type ?read.delim in the R console to display the documentation of the function.

rawdata <- read.delim("edgeR_example1_Tuch.tab", sep = "\t")

To check that the data was loaded properly we can use functions such as head (to displays the first lines of the table), dim (to display the dimensions of the table) and summary (to display summary statistics for each column). In RStudio you can also type View(rawdata) to view the full table on a separate window.

head(rawdata)
##       idRefSeq   N8 T8   N33 T33   N51  T51
## 1    NM_182502 2592  3  7805 321  3372    9
## 2    NM_003280 1684  0  1787   7  4894  559
## 3    NM_152381 9915 15 10396  48 23309 7181
## 4    NM_022438 2496  2  3585 239  1596    7
## 5 NM_001100112 4389  7  7944  16  9262 1818
## 6    NM_017534 4402  7  7943  16  9244 1815
dim(rawdata)
## [1] 15668     7
summary(rawdata)
##       idRefSeq           N8                 T8                N33        
##  NM_000014:    1   Min.   :     2.0   Min.   :     0.0   Min.   :     1  
##  NM_000016:    1   1st Qu.:   119.0   1st Qu.:    88.0   1st Qu.:   143  
##  NM_000017:    1   Median :   256.0   Median :   219.0   Median :   291  
##  NM_000018:    1   Mean   :   771.6   Mean   :   646.2   Mean   :  1270  
##  NM_000019:    1   3rd Qu.:   562.0   3rd Qu.:   503.0   3rd Qu.:   622  
##  NM_000020:    1   Max.   :393801.0   Max.   :330105.0   Max.   :581364  
##  (Other)  :15662                                                         
##       T33              N51               T51        
##  Min.   :     0   Min.   :      5   Min.   :     0  
##  1st Qu.:   171   1st Qu.:    317   1st Qu.:   223  
##  Median :   379   Median :    679   Median :   494  
##  Mean   :  1186   Mean   :   2162   Mean   :  1394  
##  3rd Qu.:   828   3rd Qu.:   1479   3rd Qu.:  1100  
##  Max.   :365430   Max.   :1675945   Max.   :633871  
##

For convenience, we separate the table in two: one containing the counts for all samples (columns 2 to 7), and another containing only the list of gene names (column 1).

rawcounts <- rawdata[, 2:7]
genes <- rawdata[, 1]


Simple pairwise differential expression analysis with edgeR GLMs

We need to import edgeR into the R environment.

library(edgeR)
## Loading required package: limma

We start by telling edgeR where our raw counts are, and calculate normalization factors.

y <- DGEList(counts=rawcounts, genes=genes)
y <- calcNormFactors(y)

y$samples
##     group lib.size norm.factors
## N8      1 12090121    1.0746801
## T8      1 10123913    1.1050377
## N33     1 19890767    0.7446805
## T33     1 18590376    1.0271063
## N51     1 33878462    0.9385600
## T51     1 21832978    1.1729954

After normalization, we can now produce a Multidimensional Scaling Plot (MDS) using the function plotMDS. This visualization, a type of dimensional reduction technique, places the samples on a plane such that the distance between samples approximates the typical log2 fold-changes between them.

plotMDS(y)

We now define the design of our comparison. We want to compare Tumor to Non-Tumor samples. So we first create a variable indicating which samples are from normal (N) or tumor (T) tissue. Then we define the design for the genewise linear models. Here, the ~ Tissue design is equivalent to a simple pairwise test of Tumor vs Non-tumor (i.e. the model only takes into account the originating tissue).

Tissue <- factor(c("N","T","N","T","N","T"))

design <- model.matrix(~ Tissue)
rownames(design) <- colnames(y)

design
##     (Intercept) TissueT
## N8            1       0
## T8            1       1
## N33           1       0
## T33           1       1
## N51           1       0
## T51           1       1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$Tissue
## [1] "contr.treatment"

Next we use this design to conduct the test of differential expression. In edgeR, this is done in 3 steps: estimation of the negative binomial dispersions (estimateDisp), fitting of the negative binomial model to the count data (glmFit) and hypothesis testing (glmLRT).

y <- estimateDisp(y, design, robust=TRUE)
fit <- glmFit(y, design)
lrt <- glmLRT(fit)

We now check how many genes were differentially expressed.

summary(decideTestsDGE(lrt))
##        TissueT
## Down       956
## NotSig   14335
## Up         377

In edgeR we make an MA-plot with the plotMD function. Up-regulated genes are indicated in red, and down-regulated genes are indicated in blue. The horizontal lines indicate 2x fold-changes.

plotMD(lrt)
abline(h=c(-1, 1), col="blue")

We can retrieve a table with all the results of differential expression using the topTags function. We also save it to a file so we can latter open it in Excel.

result <- as.data.frame(topTags(lrt, n = nrow(rawcounts)))

head(result)
##           genes     logFC   logCPM       LR       PValue          FDR
## 15660 NM_198964  3.969929 5.653315 70.20133 5.355024e-17 8.390251e-13
## 102   NM_004320 -4.472638 5.464585 66.53939 3.429713e-16 2.176337e-12
## 103   NM_173201 -4.462351 5.454729 66.15553 4.167099e-16 2.176337e-12
## 200   NM_031469 -4.022319 5.038015 64.88296 7.948160e-16 3.033213e-12
## 53    NM_005609 -5.241483 5.485800 64.49462 9.679643e-16 3.033213e-12
## 15654 NM_198966  3.872656 5.256142 63.37970 1.704659e-15 4.451432e-12
write.table(result, file = "edgeR_Tuch_Tumor_vs_NonTumor.csv", sep="\t", row.names = FALSE)


A more complex design: adding patient pairing information

Recall that tumor and non-samples were collected from 3 patients. Until now we have ignored this information in our design. Here we repeat the analysis by adding the sample pairing information to our model design, that will allow us to adjust for differences between patients.

For this we only have to change the design definition. We create a new Patient variable, and then include it as a blocking factor in the GLM design.

Patient <- factor(c(8, 8, 33, 33, 51, 51))
Tissue <- factor(c("N","T","N","T","N","T"))

design <- model.matrix(~ Patient + Tissue)
rownames(design) <- colnames(y)

design
##     (Intercept) Patient33 Patient51 TissueT
## N8            1         0         0       0
## T8            1         0         0       1
## N33           1         1         0       0
## T33           1         1         0       1
## N51           1         0         1       0
## T51           1         0         1       1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Patient
## [1] "contr.treatment"
##
## attr(,"contrasts")$Tissue
## [1] "contr.treatment"
y <- DGEList(counts=rawcounts, genes=genes)
y <- calcNormFactors(y)
y <- estimateDisp(y, design, robust=TRUE)
fit <- glmFit(y, design)
lrt <- glmLRT(fit)

de <- decideTestsDGE(lrt)
summary(de)
##        TissueT
## Down      1430
## NotSig   13740
## Up         498
plotMD(lrt)
abline(h=c(-1, 1), col="blue")

result_paired <- as.data.frame(topTags(lrt, n = nrow(rawcounts)))
write.table(result_paired, file = "edgeR_Tuch_Tumor_vs_NonTumor_paired.csv", sep="\t", row.names = FALSE)


Back

Back to previous page.