2.1 CPM/TPM

For analysis other than DGE, the data set must be corrected before use. The most basic correction required is sequencing depth. This is achieved using rescaling the counts to counts per million.

We read in our filtered count data and metadata. You can use the files you have created or read from the provided URL.

Filtered data is here https://raw.githubusercontent.com/NBISweden/course_rnaseq_slu/master/data/counts_filtered.txt and the metadata is here https://raw.githubusercontent.com/NBISweden/course_rnaseq_slu/master/data/metadata_raw.csv.

cf <- read.table("./data/counts_filtered.txt",header=TRUE)
mr <- read.csv2("./data/metadata_raw.csv",header=TRUE,stringsAsFactors=F)
rownames(mr) <- mr$Sample_ID

We can use the function cpm() from R package edgeR to convert to counts-per-million.

R

cc <- edgeR::cpm(cf)
boxplot(log10(as.matrix(cc)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="CPM")

But, CPM data has some drawbacks. It is not suitable for within-sample comparisons. The total number of reads per sample varies from sample to sample. This also makes it harder to compare one experiment to another. In addition, gene length is not controlled for in this correction. RPKM/FPKM normalisations correct for gene length, but they are not recommended because they are not comparable between samples.

A better correction method that resolves these issues is TPM (transcripts-per-million). The code for computing TPM is simple.

#' @title Compute TPM from a read count matrix
#' @param counts A numeric data.frame of read counts with samples (columns) and genes (rows).
#' @param len A vector of gene cds length equal to number of rows of dfr.
#'
#' https://support.bioconductor.org/p/91218/
#'
tpm <- function(counts,len) {
  x <- counts/len
  return(t(t(x)*1e6/colSums(x)))
}

We read in the annotation data from URL: https://raw.githubusercontent.com/NBISweden/course_rnaseq_slu/master/data/human_genes.txt. Then, we remove duplicated ensembl IDs and compute gene lengths.

g <- read.delim("./data/human_genes.txt",header=T,stringsAsFactors=F)
g <- g[!duplicated(g$ensembl_gene_id),]
g$len <- g$end_position-g$start_position
rownames(g) <- g$ensembl_gene_id

Next, we find shared genes between count data and annotation data and match their order.

igenes <- intersect(rownames(cf),g$ensembl_gene_id)
g1 <- g[igenes,]
cf1 <- cf[igenes,]
all.equal(rownames(cf1),g1$ensembl_gene_id)

## [1] TRUE

R

ct <- tpm(cf1,g1$len)
boxplot(log10(as.matrix(ct)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="TPM")

2.2 DESeq2

DESeq2 internally corrects counts for sequencing depth and RNA compositional bias using Median of ratios method. The details of this method are described further below under DGE size factors. To run this method, we create a DESeq2 object using the count data and metadata.

R

library(DESeq2)
mr$Time <- factor(mr$Time)
d <- DESeqDataSetFromMatrix(countData=cf,colData=mr,design=~Time)
d <- DESeq2::estimateSizeFactors(d,type="ratio")
cd <- counts(d,normalized=TRUE)
boxplot(log10(as.matrix(cd)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="DESeq2")

2.3 VST

For the purpose of exploratory analysis such as MDS, PCA, clustering etc, VST (variance-stabilizing-transformation) is recommended. VST is also run using DESeq2. As in the previous step, a DESeq2 object is created.

R

library(DESeq2)
mr$Time <- factor(mr$Time)
d <- DESeqDataSetFromMatrix(countData=cf,colData=mr,design=~Time)
d <- DESeq2::estimateSizeFactors(d,type="ratio")
d <- DESeq2::estimateDispersions(d)
cv <- as.data.frame(assay(varianceStabilizingTransformation(d,blind=T)),check.names=F)

# write vst counts to working directory
if(!file.exists("./data/counts_vst.txt")) write.table(cv,"./data/counts_vst.txt",sep="\t",dec=".",quote=FALSE)

# plot
boxplot(log10(as.matrix(cv)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="VST")

The effect of VST transformation can be clearly seen in a mean vs variance plot.

R

rowVar <- function(x) apply(x,1,var)
par(mfrow=c(2,2))
plot(log10(rowMeans(cf)+1),log10(rowVar(cf)+1),xlab=expression('Log'[10]~'Mean count'),ylab=expression('Log'[10]~'Variance'),main="Filtered")
plot(log10(rowMeans(ct)+1),log10(rowVar(ct)+1),xlab=expression('Log'[10]~'Mean count'),ylab=expression('Log'[10]~'Variance'),main="TPM")
plot(log10(rowMeans(cd)+1),log10(rowVar(cd)+1),xlab=expression('Log'[10]~'Mean count'),ylab=expression('Log'[10]~'Variance'),main="DESeq2")
plot(rowMeans(cv),rowVar(cv),xlab='Mean count',ylab='Variance',main="VST")
par(mfrow=c(1,1))

For RNA-seq data, as the mean count value increases, the variance increases. There is a strong almost linear relationship as seen in the figures. The statistical methods such as PCA expects similar variance across the range of mean values. If not, the higher variance genes will contribute more than the lower variance genes. Such data is said to be heteroscedastic and needs to be corrected. Correction using log transformation (with pseudocount) now gives inflates the contribution of the low variance genes. To obtain similar variance across the whole range of mean values, DESeq2 offers two methods VST (variance stabilising transformation) and RLOG (regularised log transformation).

As the name suggests, VST transformation stabilizes variance across the whole range of count values. VST is recommended for clustering or visualisation. It is not intended for differential gene expression. If the size factors vary dramatically between samples, then RLOG transformation is recommended.

Finally, we can compare all of the various transformations in a single plot.

R

par(mfrow=c(1,4))
boxplot(log10(as.matrix(cf)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="Filtered")
boxplot(log10(as.matrix(ct)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="TPM")
boxplot(log10(as.matrix(cd)+1),ylab=expression('Log'[10]~'Read counts'),las=2,main="DESeq2")
boxplot(as.matrix(cv),ylab='Read counts',las=2,main="VST")
par(mfrow=c(1,1))

Would it be possible to have one perfect normalisation method for all types of analyses? Is there any drawback to using gene length corrected counts in differential gene expression analyses?

3.1 Correlation

It is a good idea to start by checking correlation between samples. RNA-Seq samples generally have very high correlation (R^2 > 0.9). R^2 values below 0.8 may be an indication of an outlier sample. Depending on other QC metric, these low correlation samples may be discarded from analyses.

We use the function cor() for computing sample-to-sample Spearman correlation. Note that the input matrix has genes as rows and samples as columns. This generates a sample-to-sample pairwise correlation matrix. This matrix can be plotted as a heatmap using the pheatmap() function from the pheatmap R package.

R

dmat <- as.matrix(cor(cv,method="spearman"))

library(pheatmap)
pheatmap(dmat,border_color=NA,annotation_col=mr[,"Time",drop=F],
         annotation_row=mr[,"Time",drop=F],annotation_legend=T)

In the matrix, red colour denotes higher correlation (higher similarity) and blue denotes lower correlation (lower similarity). pheatmap() also hierarchically clusters rows and columns based on correlation values. The dendrograms show how samples cluster. Annotation colours denote Time groups. Notice that samples group by Time.

3.2 PCA

To run PCA, we use the R function prcomp(). It takes in a matrix where samples are rows and variables are columns. Therefore, we transpose our count matrix using the function t(). If we do not transpose, then PCA is run on the genes rather than the samples. The next line of code plots the variance explained by the top PCs.

R

pcaobj <- prcomp(x=t(cv))
{barplot(round(pcaobj$sdev^2/sum(pcaobj$sdev^2)*100,2),las=2,
        names.arg=colnames(pcaobj$x),ylab="% Variance explained",
        xlab="PCA principal components")
abline(h=2, lty=2)}

The first two principal components in total explain 85% (75%+10%) of the variance in the data set. This means a scatterplot of PC1 vs PC2 can help to visualise the most important trend in the data. Then we merge the rotated data with the metadata and plot a scatterplot coloured by our variable of interest (Time).

R

pcamat1 <- as.data.frame(pcaobj$x)
pcamat2 <- merge(pcamat1,mr,by=0)

ggplot(pcamat2,aes(PC1,PC2,colour=Time))+
  geom_point()+
  theme_bw()

Samples cluster by the Time variable as expected.

Sometimes the first two PCs may not be the ones that will best separate the sample groups, so it is a good idea to look at more PCs.

R

p1 <- ggplot(pcamat2,aes(PC1,PC3,colour=Time))+
  geom_point()+
  theme_bw()
p2 <- ggplot(pcamat2,aes(PC2,PC3,colour=Time))+
  geom_point()+
  theme_bw()
ggpubr::ggarrange(p1,p2,nrow=1,ncol=2)

An alternative way to create a PCA plot is directly from the DESeq2 object using the DESeq2 function plotPCA().

R

plotPCA(varianceStabilizingTransformation(d),intgroup="Time")

3.3 Clustering

For clustering, we create a sample-to-sample pairwise distance matrix (here we use euclidean distance). The rows and columns of this matrix is then hierarchically clustered. We use the function dist() to compute the distance. Note that for a sample-to-sample matrix, the rows need to be samples and columns should be genes. Therefore, we use the function t() to transpose our VST normalised count matrix.

R

dmat <- as.matrix(dist(t(cv)))

library(pheatmap)
pheatmap(dmat,border_color=NA,annotation_col=mr[,"Time",drop=F],
         annotation_row=mr[,"Time",drop=F],annotation_legend=T)

Hierarchically clustered sample-to-sample euclidean distance matrix. Larger distances mean lower similarity and vice-versa. In the matrix, red colour denotes larger distance (lower similarity) and blue denotes small distance (higher similarity). Annotation colours denote Time groups. The dendrogram helps to visualise sample clustering. Notice that samples group by Time.

Try to run the PCA using one of the other normalisation methods, say logCPM and/or DESeq2 normalised counts.
How different are the results? How do samples group? Can these differences be explained?

4.1 Size factors

The first step is estimating size factors. The data is normalised for sequencing depth and compositional bias as done for the VST step. DESeq2 uses a method called median-of-ratios for this step.

d <- DESeq2::estimateSizeFactors(d,type="ratio")

Optional

For those interested in the details of the median-of-ratios method, click below.

This is a step-by-step guide to computing normalisation factors (size factors) using the median-of-ratios method.

1. The geometric mean is computed across all samples for each gene.

R

gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

gmean <- apply(cf,1,gm_mean)
head(gmean)

## ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 
##       479.11181       820.01611        81.68888       319.29571 
## ENSG00000001036 ENSG00000001084 
##       920.13961       520.55995

2. Read counts for each gene is divided by the geometric mean of that gene to create a ratio.

R

ratio <- cf/gmean
head(ratio)[,1:5]

##                  Sample_1  Sample_2  Sample_3 Sample_4  Sample_5
## ENSG00000000003 0.6699898 0.6324202 0.4257879 1.026900 0.9496739
## ENSG00000000419 0.8487638 0.8048622 0.5755984 1.159733 1.1743672
## ENSG00000000457 0.7222525 0.6610447 0.5386290 1.334331 0.8936344
## ENSG00000000460 1.2496253 1.2684166 0.7391267 1.393692 1.4218794
## ENSG00000001036 0.9172521 0.8781276 0.5379618 1.157433 1.1270029
## ENSG00000001084 1.0277395 0.8798218 0.5570924 1.386968 1.3312588

3. The median ratio for each sample (across all genes) is taken as the size factor for that sample.

R

sf <- apply(ratio,2,median)
sf

##  Sample_1  Sample_2  Sample_3  Sample_4  Sample_5  Sample_6  Sample_7 
## 0.8975253 0.8407580 0.5086763 1.1263621 1.0931051 0.8123073 0.7558324 
##  Sample_8  Sample_9 Sample_10 Sample_11 Sample_12 
## 1.1745210 1.0193067 1.3720949 1.2387549 1.8679810

We can verify that these values are correct by comparing with size factors generated by DESeq2.

R

# deseq2 size factors
sizeFactors(d)

##  Sample_1  Sample_2  Sample_3  Sample_4  Sample_5  Sample_6  Sample_7 
## 0.9003753 0.8437393 0.5106445 1.1276451 1.0941383 0.8133849 0.7553903 
##  Sample_8  Sample_9 Sample_10 Sample_11 Sample_12 
## 1.1744008 1.0189325 1.3642797 1.2325485 1.8555904

If we plot the size factors for each sample against the total counts for each sample, we get the plot below. We can see that they correlate very well. Size factors are mostly correcting for total counts, ie; sequencing depth.

R

plot(sizeFactors(d),colSums(cf),xlab="Size factors",ylab="Total counts")

The raw counts can then be divided by the size factor to yield normalised read counts.

R

# custom
head(t(t(cf)/sf))[,1:5]
# deseq2
head(counts(d,normalized=TRUE))[,1:5]

##                  Sample_1  Sample_2  Sample_3  Sample_4  Sample_5
## ENSG00000000003 357.65009 360.38907 401.04088 436.80445 416.24541
## ENSG00000000419 775.46562 785.00590 927.89851 844.31105 880.97655
## ENSG00000000457  65.73631  64.22776  86.49901  96.77172  66.78223
## ENSG00000000460 444.55572 481.70816 463.94925 395.07720 415.33059
## ENSG00000001036 940.36349 961.03752 973.11390 945.52183 948.67361
## ENSG00000001084 596.08349 544.74652 570.10713 641.00166 633.97378
##                  Sample_1  Sample_2  Sample_3  Sample_4  Sample_5
## ENSG00000000003 356.51801 359.11565 399.49511 436.30747 415.85236
## ENSG00000000419 773.01102 782.23211 924.32203 843.35041 880.14467
## ENSG00000000457  65.52823  64.00081  86.16561  96.66161  66.71917
## ENSG00000000460 443.14856 480.00607 462.16101 394.62769 414.93840
## ENSG00000001036 937.38692 957.64173 969.36314 944.44605 947.77780
## ENSG00000001084 594.19669 542.82168 567.90972 640.27234 633.37514

The function estimateSizeFactors() has options to set a custom locfunc other than the median. Why is this useful? What happens if you change it to “shorth”. Check out the help page for estimateSizeFactors().

4.2 Gene dispersion

When it comes to comparing values between groups, some measure of variation is needed to estimate the variability in gene counts within groups. Dispersion is a measure of variation in a data set. Variance and standard deviation are not a good measure to estimate variability because it correlates with the mean.

Optional

For some more discussion on dispersion, click below.

We can create a mean counts vs variance plot for all genes in our data set.

R

dm <- apply(cd,1,mean)
dv <- apply(cd,1,var)

ggplot(data.frame(mean=log10(dm+1),var=log10(dv+1)),
       aes(mean,var))+
  geom_point(alpha=0.2)+
  geom_smooth(method="lm")+
  labs(x=expression('Log'[10]~'Mean counts'),y=expression('Log'[10]~'Variance'))+
  theme_bw()

We see a mostly linear relationship on the log scale. The blue line denotes a linear fit. Genes that have larger read counts show higher variance. It’s hard to say which genes are more variable based on this alone. Therefore, variance is not a good measure to identify variation in read counts. A measure that controls for this mean-variance relationship is what we need.

One option is the coefficient of variation (CV).

R

cva <- function(x) sd(x)/mean(x)
dc <- apply(cd,1,cva)

ggplot(data.frame(mean=log10(dm+1),var=dc),
       aes(mean,var))+
  geom_point(alpha=0.2)+
  geom_smooth()+
  labs(x=expression('Log'[10]~'Mean counts'),y="Coefficient of variation")+
  theme_bw()

Now, we see that genes with lower counts have higher variability and genes with larger counts have lower variability. A measure like CV is taking the ratio of ‘variation’ to mean. cv=sd(x)/mean(x).

This becomes even more apparent if we compute the CV and mean over replicates within sample groups (Time).

R

dx1 <- data.frame(t0=apply(cd[,1:3],1,cva),t2=apply(cd[,4:6],1,cva),
            t6=apply(cd[,7:9],1,cva),t24=apply(cd[,10:12],1,cva))
dx1$gene <- rownames(dx1)
dx1 <- tidyr::gather(dx1,key=sample,value=cv,-gene)
rownames(dx1) <- paste0(dx1$gene,"-",dx1$sample)

dx2 <- data.frame(t0=apply(cd[,1:3],1,mean),t2=apply(cd[,4:6],1,mean),
            t6=apply(cd[,7:9],1,mean),t24=apply(cd[,10:12],1,mean))
dx2$gene <- rownames(dx2)
dx2 <- tidyr::gather(dx2,key=sample,value=mean,-gene)
rownames(dx2) <- paste0(dx2$gene,"-",dx2$sample)

dx3 <- merge(dx1,dx2,by=0)

ggplot(dx3,aes(x=log10(mean+1),y=cv))+
  geom_point(alpha=0.2)+
  geom_smooth()+
  facet_wrap(~sample.x)+
  labs(x=expression('Log'[10]~'Mean counts'),y="Coefficient of variation")+
  theme_bw()

We find that CV strongly declines with increasing counts. Genes with low counts show higher dispersion. For the sake of completeness, we can also plot the relationship between CV and variance for the same sample groups.

R

dx1 <- data.frame(t0=apply(cd[,1:3],1,cva),t2=apply(cd[,4:6],1,cva),
            t6=apply(cd[,7:9],1,cva),t24=apply(cd[,10:12],1,cva))
dx1$gene <- rownames(dx1)
dx1 <- tidyr::gather(dx1,key=sample,value=cv,-gene)
rownames(dx1) <- paste0(dx1$gene,"-",dx1$sample)

dx2 <- data.frame(t0=apply(cd[,1:3],1,var),t2=apply(cd[,4:6],1,var),
            t6=apply(cd[,7:9],1,var),t24=apply(cd[,10:12],1,var))
dx2$gene <- rownames(dx2)
dx2 <- tidyr::gather(dx2,key=sample,value=var,-gene)
rownames(dx2) <- paste0(dx2$gene,"-",dx2$sample)

dx3 <- merge(dx1,dx2,by=0)

ggplot(dx3,aes(x=log10(var+1),y=cv))+
  geom_point(alpha=0.2)+
  geom_smooth()+
  facet_wrap(~sample.x)+
  labs(x=expression('Log'[10]~'Variance'),y="Coefficient of variation")+
  theme_bw()

DESeq2 computes it’s own version of dispersion in a more robust manner taking into account low count values. The DESeq2 dispersion estimates are inversely related to the mean and directly related to variance. The dispersion estimate is a good measure of the variation in gene expression for a certain mean value.

Now, the variance or dispersion estimate for genes with low counts is unreliable when there are too few replicates. To overcome this, DESeq2 borrows information from other genes. DESeq2 assumes that genes with similar expression levels have similar dispersion values. Dispersion estimates are computed for each gene separately using maximum likelihood estimate. A curve is fitted to these gene-wise dispersion estimates. The gene-wise estimates are then ‘shrunk’ to the fitted curve.

Gene-wide dispersions, fitted curve and shrinkage can be visualised using the plotDispEsts() function.

R

d <- DESeq2::estimateDispersions(d)
plotDispEsts(d)

The black points denote the maximum likelihood dispersion estimate for each gene. The red curve denote the fitted curve. The blue points denote the new gene dispersion estimates after shrunk towards the curve. The circled blue points denote estimates that are not shrunk as they are too far away from the curve. Thus, shrinkage method is important to reduce false positives in DGE analysis involving too few replicates.

It is a good idea to visually check the dispersion shrinkage plot to verify that the modelling works for your data set.

4.3 Testing

Overdispersion is the reason why RNA-Seq data is better modelled as negative-binomial distribution rather than poisson distribution. Poisson distribution has a mean = variance relationship, while negative-binomial distribution has a variance > mean relationship. The last step in the DESeq2 workflow is fitting the Negative Binomial model for each gene and performing differential expression testing. This is based on the log fold change values computed on the corrected count estimates between groups.

logFC = log2 (corrected counts group A / corrected counts group B)

The most commonly used testing for comparing two groups in DESeq2 is the Walds’s test. The null hypothesis is that the groups are not different and logFC=0. The list of contrasts can be seen using resultsNames(). Then we can pick our comparisons of interest.

dg <- nbinomWaldTest(d)
resultsNames(dg)

## [1] "Intercept"      "Time_t2_vs_t0"  "Time_t24_vs_t0" "Time_t6_vs_t0"

And we can get the result tables for the three different comparisons. The summary of the result object shows the number of genes that are differentially expressed with positive or negative fold-change and outliers.

res1 <- results(dg,name="Time_t2_vs_t0",alpha=0.05)
summary(res1)
res2 <- results(dg,name="Time_t6_vs_t0",alpha=0.05)
summary(res2)
res3 <- results(dg,name="Time_t24_vs_t0",alpha=0.05)
summary(res3)

## 
## out of 14578 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 413, 2.8%
## LFC < 0 (down)     : 696, 4.8%
## outliers [1]       : 0, 0%
## low counts [2]     : 2261, 16%
## (mean count < 26)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## 
## out of 14578 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2484, 17%
## LFC < 0 (down)     : 2856, 20%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## 
## out of 14578 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 5066, 35%
## LFC < 0 (down)     : 5093, 35%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

You can also build up the comparison using contrasts. Contrats need the condition, level to compare and the reference level (base level). For example, results(dg,contrast=c("Time","t2","t0"),alpha=0.05).

Note that in both cases above, t0 is the reference level and other levels are compared to this. Therefore, a fold-change of 2 would mean that, a gene is 2 fold higher expressed in the test level compared to t0.

The results() function has many useful arguments. One can set a threshold on the logFC values using lfcThreshold. By default, no filtering is performed based on logFC values. Outliers are detected and p-values are set to NA automatically using cooksCutoff. independentFiltering remove genes with low counts.

head(res1)

## log2 fold change (MLE): Time t2 vs t0 
## Wald test p-value: Time t2 vs t0 
## DataFrame with 6 rows and 6 columns
##                         baseMean     log2FoldChange              lfcSE
##                        <numeric>          <numeric>          <numeric>
## ENSG00000000003 490.017213130775  0.220619809383656  0.112761083788383
## ENSG00000000419 817.780657499214 0.0592719659601406  0.101481284662652
## ENSG00000000457  82.078767893562  0.207748623172909  0.220404896514001
## ENSG00000000460  356.07160020304 -0.129186381676072  0.115139182335103
## ENSG00000001036 919.606750211436 0.0288826917736355 0.0851500691928604
## ENSG00000001084 529.593965301543  0.211964772147193 0.0929810563117334
##                              stat             pvalue              padj
##                         <numeric>          <numeric>         <numeric>
## ENSG00000000003  1.95652437854969 0.0504034135093116 0.263505451695327
## ENSG00000000419 0.584067950629267  0.559174596319379 0.830262316260423
## ENSG00000000457 0.942577168015466   0.34589722283592 0.689945926100409
## ENSG00000000460 -1.12200190288033    0.2618616308297 0.612624613846088
## ENSG00000001036  0.33919751384133  0.734460942036292 0.909638554355053
## ENSG00000001084  2.27965545408033 0.0226281312376226 0.159263252830742

The results table contains mean expression value (baseMean), log2 fold change (log2FoldChange), log2 fold change standard error (lfcSE), wald test statistic (stat), wald test p-value (pvalue) and BH adjusted p-value (padj) for each gene.

It is a good idea to look at the distribution on unadjusted p-values.

R

hist(res1$pvalue[res1$baseMean>1],main="res1 Pval distribution",xlab="P-values")

This is the kind of distribution to be expected when the p-values are “well-behaved”. For more explanation on p-value distributions, see here. If this distribution is very different or strange, then it might indicate an underlying problem.

Note that the results object is a DESeqResults class object and not a data.frame. It can be converted to a data.frame using as.data.frame() for exporting to a file.

We can filter the results table as needed.

# convert res to data.frame
df_res1 <- na.omit(as.data.frame(res1))
df_res2 <- na.omit(as.data.frame(res2))
df_res3 <- na.omit(as.data.frame(res3))

# all genes
nrow(df_res1)
# only genes with padj <0.05
nrow(df_res1[df_res1$padj<0.05,])
# only genes with padj <0.05 and an absolute fold change >2
nrow(df_res1[df_res1$padj<0.05 & abs(df_res1$log2FoldChange)>2,])

## [1] 12317
## [1] 1109
## [1] 11

Note that manually filtering by log2FoldChange on the results table is not exactly the same as setting the lfcThreshold argument in the results() function.

Finally, we can add additional information to our results to make the interpretation and downstream analyses easier. We read in the gene info downloaded through biomaRt.

# add row names as a new column
df_res1$ensembl_gene_id <- rownames(df_res1)
df_res2$ensembl_gene_id <- rownames(df_res2)
df_res3$ensembl_gene_id <- rownames(df_res3)

# read genes info
hg <- read.delim("./data/human_genes.txt",header=T,sep="\t",stringsAsFactors=F)
hg <- hg[!duplicated(hg$ensembl_gene_id),]

# merge dataframes
df_res1 <- merge(df_res1,hg,by="ensembl_gene_id")
df_res2 <- merge(df_res2,hg,by="ensembl_gene_id")
df_res3 <- merge(df_res3,hg,by="ensembl_gene_id")

head(df_res1)

##   ensembl_gene_id  baseMean log2FoldChange      lfcSE       stat
## 1 ENSG00000000003 490.01721     0.22061981 0.11276108  1.9565244
## 2 ENSG00000000419 817.78066     0.05927197 0.10148128  0.5840680
## 3 ENSG00000000457  82.07877     0.20774862 0.22040490  0.9425772
## 4 ENSG00000000460 356.07160    -0.12918638 0.11513918 -1.1220019
## 5 ENSG00000001036 919.60675     0.02888269 0.08515007  0.3391975
## 6 ENSG00000001084 529.59397     0.21196477 0.09298106  2.2796555
##       pvalue      padj entrezgene external_gene_name chromosome_name
## 1 0.05040341 0.2635055       7105             TSPAN6               X
## 2 0.55917460 0.8302623       8813               DPM1              20
## 3 0.34589722 0.6899459      57147              SCYL3               1
## 4 0.26186163 0.6126246      55732           C1orf112               1
## 5 0.73446094 0.9096386       2519              FUCA2               6
## 6 0.02262813 0.1592633       2729               GCLC               6
##   start_position end_position strand   gene_biotype
## 1      100627109    100639991     -1 protein_coding
## 2       50934867     50958555     -1 protein_coding
## 3      169849631    169894267     -1 protein_coding
## 4      169662007    169854080      1 protein_coding
## 5      143494811    143511690     -1 protein_coding
## 6       53497341     53616970     -1 protein_coding
##                                                                                      description
## 1                                              tetraspanin 6 [Source:HGNC Symbol;Acc:HGNC:11858]
## 2 dolichyl-phosphate mannosyltransferase subunit 1, catalytic [Source:HGNC Symbol;Acc:HGNC:3005]
## 3                                   SCY1 like pseudokinase 3 [Source:HGNC Symbol;Acc:HGNC:19285]
## 4                        chromosome 1 open reading frame 112 [Source:HGNC Symbol;Acc:HGNC:25565]
## 5                                        alpha-L-fucosidase 2 [Source:HGNC Symbol;Acc:HGNC:4008]
## 6                 glutamate-cysteine ligase catalytic subunit [Source:HGNC Symbol;Acc:HGNC:4311]

Now we have annotated our DEG list with information such as gene names, genomic coordinates, biotype etc.

Save the data to a text file for use in functional anlayses.

write.table(df_res1,"./data/deg-time_t2_vs_t0.txt",sep="\t",col.names=TRUE,row.names=FALSE,quote=F)
write.table(df_res2,"./data/deg-time_t6_vs_t0.txt",sep="\t",col.names=TRUE,row.names=FALSE,quote=F)
write.table(df_res3,"./data/deg-time_t24_vs_t0.txt",sep="\t",col.names=TRUE,row.names=FALSE,quote=F)

How many DE genes do you find for each test? Is there a trend? How much do they change of you set a threshold of 1 on fold-change?

nrow(df_res1[df_res1$padj<0.05,])
nrow(df_res2[df_res2$padj<0.05,])
nrow(df_res3[df_res3$padj<0.05,])

nrow(df_res1[df_res1$padj<0.05 & abs(df_res1$log2FoldChange)>1,])
nrow(df_res2[df_res2$padj<0.05 & abs(df_res2$log2FoldChange)>1,])
nrow(df_res3[df_res3$padj<0.05 & abs(df_res3$log2FoldChange)>1,])

## [1] 1086
## [1] 5222
## [1] 9909
## [1] 104
## [1] 788
## [1] 4082

4.3.1 lfcShrink

This is an extra step to generate more accurate log2 fold changes. This step corrects the log2 fold changes for genes with high dispersion.

R

lres1 <- lfcShrink(dg,coef="Time_t2_vs_t0",res=res1)
summary(lres1)

## 
## out of 14578 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 413, 2.8%
## LFC < 0 (down)     : 696, 4.8%
## outliers [1]       : 0, 0%
## low counts [2]     : 2261, 16%
## (mean count < 26)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

R

head(lres1)

## log2 fold change (MAP): Time t2 vs t0 
## Wald test p-value: Time t2 vs t0 
## DataFrame with 6 rows and 6 columns
##                         baseMean      log2FoldChange              lfcSE
##                        <numeric>           <numeric>          <numeric>
## ENSG00000000003 490.017213130775   0.166374521328742 0.0903392402749302
## ENSG00000000419 817.780657499214  0.0494916503724324 0.0846342204609171
## ENSG00000000457  82.078767893562  0.0925123104527917  0.114348354878349
## ENSG00000000460  356.07160020304 -0.0895216223688513  0.091747791502412
## ENSG00000001036 919.606750211436  0.0257185248646338 0.0746622716252995
## ENSG00000001084 529.593965301543   0.189029922610361 0.0795917180410762
##                              stat             pvalue              padj
##                         <numeric>          <numeric>         <numeric>
## ENSG00000000003  1.95652437854969 0.0504034135093116 0.263505451695327
## ENSG00000000419 0.584067950629267  0.559174596319379 0.830262316260423
## ENSG00000000457 0.942577168015466   0.34589722283592 0.689945926100409
## ENSG00000000460 -1.12200190288033    0.2618616308297 0.612624613846088
## ENSG00000001036  0.33919751384133  0.734460942036292 0.909638554355053
## ENSG00000001084  2.27965545408033 0.0226281312376226 0.159263252830742

We see that the number of genes up and down has not changed. But, we can see that the logFC distribution has changed in the density plot below.

R

par(mfrow=c(1,2))
plot(density(na.omit(res1$log2FoldChange)),main="Default log2FoldChange")
plot(density(na.omit(lres1$log2FoldChange)),main="lfcShrink log2FoldChange")
par(mfrow=c(1,1))

This step does not change the total number of DE genes. This may be useful in downstream steps where genes need to be filtered down based on fold change or if fold change values are used in functional analyses such as GSEA.

5.1 MA plot

The MA plot shows mean expression vs log fold change for all genes. The plotMA() function from DESeq2 takes a results object as input. Differentially expressed genes are marked in red.

R

DESeq2::plotMA(res1)

How does this plot change if you set log fold change filtering to minimum value of 1. How does the plot change when you use lfcShrink()?

5.2 Volcano plot

A volcano plot is similar to the MA plot. It plots log fold change vs adjusted p-values.

R

df_res1_1 <- df_res1[df_res1$padj<0.05,]

ggplot()+
  geom_point(data=df_res1,aes(x=log2FoldChange,y=-log10(padj)),col="grey80",alpha=0.5)+
  geom_point(data=df_res1_1,aes(x=log2FoldChange,y=-log10(padj)),col="red",alpha=0.7)+
  geom_hline(aes(yintercept=-log10(0.05)),alpha=0.5)+
  theme_bw()

X axis denotes log fold change and the y axis denotes -log10 adjusted p-value. The adjusted p-value is transformed so that the smallest p-values appears at the top. The horizontal grey line denotes the significance threshold line. All genes above this line (coloured red as well) are considered significant.

Why is the y-axis (p-value) on a -log scale?

5.3 Counts plot

It can be a good idea to manually verify some of these genes by plotting out it’s actual read count values. We can use the function plotCounts() to visualise the data points for a gene of interest. Below, we see the counts before and after normalisation.

R

plotCounts(d,gene=rownames(res1)[1],intgroup="Time",normalized=F)
plotCounts(d,gene=rownames(res1)[1],intgroup="Time",normalized=T)

By looking at the count plots, do you agree that the top DE genes differ significantly between the groups compared?

5.4 Heatmap

A heatmap is a good way to visualise the clustered relative expression levels of all DEGs.

We first select DE genes from the DGE results table using some criteria. Here we select genes with adjusted-p-value < 0.05 and absolute log fold change >1. A log fold change of 1 translates to a fold change of 2 (2^1=2). Then we can use this list of genes to filter the VST table for only these genes.

df_res3 <- na.omit(df_res3)
deg <- df_res3[df_res3$padj<0.05 & abs(df_res3$log2FoldChange)>1,]
cv1 <- cv[deg$ensembl_gene_id,]

We can use a heatmap function to explore the visual consequences of clustering. Here, we can make use of the pheatmap() function, which by default will do the clustering of the rows and columns.

R

cl <- pheatmap(cv1, scale="row",color=colorRampPalette(c("navy","white","firebrick"))(90),border_color=NA,cluster_cols=F,cutree_rows=2,show_rownames=FALSE)

We can see that up and down genes cluster distinctly. By setting scale="row", the data is are converted to z-scores which means that they are scaled row-wise. Therefore, the colours show fold-change up or down relative to the mean expression value for the gene. The rows are then clustered based on these z-scores.

Why do samples 10, 11 and 12 show the strongest values for up and down?