RNA-Seq - differential expression using DESeq2

---

The Snf2 dataset

The RNA-Seq dataset we will use in this practical has been produced by Gierliński et al, 2015) and (Schurch et al, 2016)).

The dataset is composed of 48 samples of yeast wild-type (WT) strain, and 48 samples of Snf2 knock-out mutant cell line. The prepared RNA-Seq libraries (unstranded) were pooled and sequenced on seven lanes of a single flow-cell on an Illumina HiSeq 2000 resulting in a total of 1 billion 50-bp single-end reads across the 96 samples. RNA-Seq reads have been cleaned, mapped and counted to generated a count data matrix containing \(7126\) rows/genes.

In this tutorial, we will start from the count table generated by the authors of the study, and use different R tools in order to

• explore the distribution of the data;
• test alternative normalisation procedures;
• select differentially expressed genes (DEG);
• cluster the count table to regroup genes having similar transcription profiles;
• compare the lists of co-regulated genes with the annotations from the Gene Ontology, in order to test if they are statistically enriched for particular functional classes.

---

Loading the dataset

R enables to download data directly from the Web with the download.file function. After, we begin with some verification steps.

The expression matrix and phenotypic information will be loaded into R using the read.table function. Both table will be converted into a data.frame object when loaded into R.

• The countTable object should contain counts for each gene (row) and each sample (column).

• The phenoTable object should contain the designation of each genotype, that is 48 biological replicates of the wild type strain (WT) and the Snf2 mutant.

Note that the created tables include column names as factor object.

Reminder: a factor is a vector with levels (categories), which permits an efficient storage and indexing, but can in some cases lead to misleading effects. To circumvent this, we will sometimes have to convert the factor to a vector, with the R command as.vector().

            WT1   WT2   WT3   WT4   WT5
yal034w-a    78   116   103   113   113
yal035w    2815  4160  2969  4411  4836
yal036c     331   599   348   525   691
yal037c-a     0    37    23    15    51
yal037c-b    70   183   199   208   196
yal037w      22    69    32    58    51
yal038w   21080 18869 14649 28551 20208
yal039c    1541  1447  1318  2251  1501

[1] "strain" "color" 

# Load the count table
countTable <- read.table("http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data/counts.txt", header=TRUE, row.names=1)

# View part of the data/table
countTable[101:108, 1:5]

# Load the sample description file 
# (traditionally named "pheno" for "phenotype", athough in our case the conditions are genotypes)
phenoTable <- read.table("http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data/expDesign.txt", header=TRUE, row.names=1)


# Define a strain-specific color for each sample, and add it as a supplementary column
strainColor <- c("WT"="green","Snf2"="orange") # Choose one color per strain
print(strainColor)

## Add a column to the pheno type to specify the color of each sample according to its genotype
phenoTable$color[phenoTable$strain == "Snf"] <- strainColor["Snf2"]
phenoTable$color[phenoTable$strain == "WT"] <- strainColor["WT"]

# Check the dimensions of the pheno table
colnames(phenoTable)

Exercice

1. Find the commands to visualize the first and last line of a data.frame object (clue: similar to bash commands)

View solution|Hide solution

head(phenoTable)
tail(phenoTable)

2. Check the dimensions of the data (tables).

   1. How many genes do we have?
   2. How many samples?
   3. Check that the dimensions of the sample table and pheno table are consistent.

View solution|Hide solution

## Get the number of genes and sample from the dimensions of the count table
nb_genes <- nrow(countTable)
print(nb_genes)

nb.samples <- ncol(countTable)
print(nb.samples)

## Check that the dimension of the pheno table is consistent with the number of samples
dim(countTable)
dim(phenoTable)
## Each row of the pheno table corresponds to one sample.
## Each column of the count table corresponds to one sample.

3. Check the number of replicates in each class using the table function. To better understand what you could do with this function, play with it (varying the parameters) with the following test table.

## Quick test with a small table generated on the flight
testTable <- data.frame(
  ABC = sample(c("A", "B", "C", "D", "E"), replace=TRUE, 20),
  gender = sample(c("boy", "girl", "robot"), replace=TRUE, 20),
  yesno = sample(c("yes", "no"), replace=TRUE, 20))

View solution|Hide solution

## Check the content of the testTable
testTable

## Get stats about column content
table(testTable)
table(testTable[,2])
table(testTable[,2],testTable[,3])

## The same with column names instead of column numbers
table(testTable$gender)
table(testTable[, c("gender", "yesno")])

## Check the number of samples per strain in the yeast count table
table(phenoTable$strain)

## Check the consistency between strains and sample colors
table(phenoTable)

4. Draw a barplot showing the total number of reads in each sample (clue: use colSums function) with sample colored (according to the color column in phenoTable). What can you say from this diagram?

View solution|Hide solution

barplot(colSums(countTable)/1000000, 
        main="Total number of reads per sample (million)",
        col=phenoTable$color, 
        #        names.arg = "", 
        las=1,  horiz=TRUE,
        ylab="Samples", cex.names=0.5,
        xlab="Million counts")

Million reads per sample.

The barplot indicates the library sizes (total number of reads) for each sample. We can see important differences, ranging from 5.57 to 15.2 million reads per sample.

This has important consequences: any read count should be interpreted relative to the sequencing depth of the corresponding sample. We will thus need to normalise the counts. Consequently, it makes not much sense to compute summary statistics per gene (mean counts, standard deviation, …) before having normalized the data.

---

Descriptive statistics

Basic statistics

Before going further in the analysis, we will compute some descriptive statistics on the dataset. At this stage we only compute statistics per sample, since statistics per gene are meaningful only after library-wise normalization of the counts.

# Min, Max, median (...). 
# Here on the first 4 samples
summary(countTable[,1:4])

Among the reference functions in R, apply (and its derivates) allow to repeated execution of a function to a serie of data, at a noteworthy higher speed tha using for loops. Note that, if working on a matrix/data frame, the apply function requires that you indicate if you want to apply the function to the rows (MARGIN=1) or to the columns (MARGIN=2). Some built-in functions already exist as colMeans, colSums, rowMeans and rowSums. Let’s use apply to get the summary for every column of the count table. It produces a matrix as output that we transpose to obtain samples as rows and the statistics as columns, and move back to a data.frame object.

statsPerSample <- data.frame(t(apply(countTable, 2, summary)))
head(statsPerSample)

    Min. X1st.Qu. Median Mean X3rd.Qu.   Max.
WT1    0     49.0    224  837      561 189000
WT2    0     88.0    371 1110      853 197000
WT3    0     74.2    297  903      670 172000
WT4    0    104.0    430 1460     1020 329000
WT5    0     84.0    353 1120      819 225000
WT6    0    153.0    667 2050     1600 357000

# Redefine column nemes for readability
names(statsPerSample) <- c("min", "Q1", "median", "mean", "Q3", "max")
kable(head(statsPerSample))

	min	Q1	median	mean	Q3	max
WT1	0	49.0	224	837	561	189000
WT2	0	88.0	371	1110	853	197000
WT3	0	74.2	297	903	670	172000
WT4	0	104.0	430	1460	1020	329000
WT5	0	84.0	353	1120	819	225000
WT6	0	153.0	667	2050	1600	357000

We can now add some columns to the stats per sample. For example, here are the library sums and the 5% quantile.

statsPerSample$libsum <- colSums(countTable) # library sums
statsPerSample$perc05 <- apply(countTable, 2, quantile, probs=0.05)

Exercice

Integrate some more columns to statsPerSample: - the 10%, 90% and 95% quantile;

View solution|Hide solution

## Add columns with some indicative percentiles
statsPerSample[, "perc10"] <- apply(countTable, 2, quantile, 0.10)
statsPerSample[, "perc90"] <- apply(countTable, 2, quantile, 0.90)
statsPerSample[, "perc95"] <- apply(countTable, 2, quantile, 0.95)

## Alternatively, this can be done in a single command
statsPerSample[,c("perc10", "perc90", "perc95")] <- 
  apply(countTable, 2, quantile, c(0.10, 0.90, 0.95))

• the number of zeros for each column;

View solution|Hide solution

## Add a column with the number of zero counts per sample
statsPerSample$zeros <- apply(countTable==0, 2, sum)

• the percentage these zeros represent. View solution|Hide solution

## Add a column to the stats per sample, with the percentage of genes with zero counts
statsPerSample$percZeros <- 100*statsPerSample$zeros/nrow(countTable)

Finally, randomly choose (with the sample() function) 10 samples (rows) of the stats per sample table and look at their statistics.

statsPerSample[sample(1:ncol(countTable), size = 10),]

We show hereafter the statistics per sample (subset only).

	min	Q1	median	mean	Q3	max	libsum	perc05	perc10	perc90	perc95	zeros	percZeros
Snf26	0	124.0	506	1380	1140	210000	9812933	0	2392	2530	2526	467	6.55
Snf2	0	98.0	422	1220	960	200000	8694091	0	1977	2683	3066	552	7.75
Snf3	0	98.0	408	1230	954	210000	8773138	0	3316	4570	5280	533	7.48
WT8	0	80.0	360	1380	900	343000	9839118	0	1624	1738	1868	576	8.08
Snf43	0	119.0	476	1340	1060	219000	9568585	0	2	1	3	482	6.76
WT37	0	75.0	344	1550	901	468000	11009738	0	2	2	2	586	8.22
WT19	0	84.0	364	1310	879	311000	9326815	0	2	1	2	564	7.92
WT39	0	61.0	280	1040	689	247000	7408720	0	4056	5098	3496	627	8.80
WT10	0	80.2	360	1110	866	188000	7920887	0	2	2	2	547	7.68
Snf11	0	109.0	458	1390	1050	247000	9907483	0	1524	2802	2181	514	7.21

Interpretation of the stats per sample

• In all samples, several hundreds of genes are completely undetected (zeros).
• The huge difference between the central tendency (mean, median) and the maximum reflects the presence of outliers: a few genes have hundreds of thousands counts.
• The mean count per genes show important fluctuations, indicating the need for normalisation.
• The median is much lower than the mean, indicating right-skewed asymmetry.

Distributions

The summary only displays a few milestone values (mean, median, quartiles). In order to get a better intuition of the data, we can draw an histogram of all values.

par(mfrow=c(3,1))

hist(as.matrix(countTable), col="blue", border="white", breaks=100)

hist(as.matrix(countTable), col="blue", border="white",
     breaks=20000, xlim=c(0,2000), main="Counts per gene",
     xlab="Counts (truncated axis)", ylab="Number of genes", 
     las=1, cex.axis=0.7)

epsilon <- 1 # pseudo-count to avoid problems with log(0)
hist(as.matrix(log2(countTable + epsilon)), breaks=100, col="blue", border="white",
     main="Log2-transformed counts per gene", xlab="log2(counts+1)", ylab="Number of genes", 
     las=1, cex.axis=0.7)

Histogram of counts per genes. Top: raw counts. the scale is determined by the gene with the highest count, which is apparently an outlier. Middle: raw counts, with X axis truncated to 2000 in order to display a representative range despite outliers. Bottom: log2-transformed counts (bottom) per gene, with a pseudocount of 1 to avoid minus infinitevalues resulting from zero counts.

par(mfrow=c(1,1))

Interpretation

• The top histogram is not very informative so far, apparently due to the presence of a few very high count values, that impose a very large scale on the \(X\) axis.
• The middle histogram shows the representative range. Note the height of the first bin, which includes the zero counts.
• The logarithmic transformation (bottom histogram) improves the readability. Note that we added a pseudo-count of 1 to avoid problems with the log transformation of zero counts (which gives \(-\infty\)).

To get better insights into the distribution per sample, the boxplot offer a good perspective.

boxplot(log2(countTable + epsilon), col=phenoTable$color, pch=".", 
        horizontal=TRUE, cex.axis=0.5,
        las=1, ylab="Samples", xlab="log2(counts +1)")

Box plots of non-normalized log2(counts+1) per sample.

Another way to get an intuition of the value distributions is to use the plotDensity function, which draws one density curve for each sample.

# We will require one function from the affy package
if(!require("affy")){
  source("http://bioconductor.org/biocLite.R")
  biocLite("affy")  
}
library(affy)
plotDensity(log2(countTable + epsilon), lty=1, col=phenoTable$color, lwd=2)
grid()
legend("topright", legend=names(strainColor), col=strainColor, lwd=2)

Densities of log2(counts+1). Each curve corresponds to one sample.

Beware: the R function plotDensity does not display the actual distribution of your values, but a polynomial fit. The representation thus generally looks smoother than the actual data. It is important to realize that, in some particular cases, the fit can lead to extrapolated values which can be misleading.

Scatter plots

Let’s have a look at the scatter plots using the pairs function. We will only represent 10 randomly selected samples.

# Define a function to draw a scatter plot for a pair of variables (samples) with density colors
plotFun <- function(x,y){ 
  dns <- densCols(x,y); 
  points(x,y, col=dns, pch=".", panel.first=grid());  
#  abline(a=0, b=1, col="brown")
  }

# Plot the scatter plot for a few pairs of variables selected at random
set.seed(123) # forces the random number generator to produce fixed results
pairs(log2(countTable[,sample(ncol(countTable),5)] + epsilon), 
      panel=plotFun, lower.panel = NULL)

Scatter plot of log2-counts for a random selection of samples.

The command pairs draws a scatter plot for each pair of columns of the input dataset. The plot shows a fairly good reproducibility between samples of the same type (WT or KO, respectively): all points are aligned along te diagonal, with a relatively wider dispersion at the bottom, corresponding to small number fluctuations.

In contrast, on all the plots comparing a WT and a KO sample, we can see some points (genes) discarding from the diagonal.

---