From Gene Counts to Differential Expression - DESeq2 Tutorial
note: this tutorial is loosely based off the DESeq2 and bioconductor documentation found here. Some explanation text is taken directly from the documentation and referenced.
In this tutorial, we will walk through going from gene counts to differential expression results.
Introduction
The true end goal in RNA Sequencing analysis is determining the genes that are differentially expressed between a control and treatment group. In this case, if you followed the previous tutorial, you know we are investigating the transcriptomic response of TCPMOH on the developing zebrafish. In this dataset subset we have two control biological replicates and two TCPMOH biological replicates. To perform differential expression analysis, we will be using using DESeq2 (Love, Huber, and Anders 2014). DESeq2 is a great tool for dealing with RNA-seq data and running Differential Gene Expression (DGE) analysis. It takes read count files from different samples, combines them into a big table (with genes in the rows and samples in the columns) and applies normalization for sequencing depth and library composition. We do not need to account for gene length normalization does because we are comparing the counts between sample groups for the same gene.
Load necessary packages
The first step is to make sure we have the proper packages loaded. We will be using DESeq2, tidyverse, and ggplot2. If you haven’t already installed these before, you can run the following commands:
# install DESeq2
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
# install tidyverse (ggplot2 is in tidyverse)
install.packages("tidyverse")
Now we can load the packages into our workspace by running the following commands:
# package setup
library( "DESeq2" )
library(tidyverse)
library(ggplot2)
Getting Started - Data
This tutorial assumes you have a file, per sample, of gene counts. This would be a typical output from a tool such as featureCounts which we used in the previous tutorial. To give an idea of what this looks like, lets peak into one of the files:
# define the path to the file
# note this path is relative to the current working directory and the files are stored in a folder called data
dmso1_data_file_path <- 'data/DMSO_1_counts.tsv'
# read in the file
dmso_1_count_data <- read.csv(dmso1_data_file_path, header = TRUE, sep = "\t")
# show the first few lines
head(dmso_1_count_data)
## Geneid DMSO_1.RNA.STAR..mapped.bam
## 1 CYTB 136445
## 2 ND6 1443
## 3 ND5 77429
## 4 ND4 71269
## 5 ND4L 12122
## 6 ND3 10141
As you can see, the file has two columns. One identifying the gene name (Geneid) and the second identifying the number of counts per gene. It is useful to go ahead and rename the second column as it is currently cumbersome.
# rename the second column
colnames(dmso_1_count_data)[2] ="DMSO_1"
Lets go ahead and do this for the rest of the samples:
# dmso2
dmso_2_count_data <- read.csv('data/DMSO_2_counts.tsv', header = TRUE, sep = "\t")
colnames(dmso_2_count_data)[2] ="DMSO_2"
# tcpmoh1
tcpmoh_1_count_data <- read.csv('data/TCPMOH_1_counts.tsv', header = TRUE, sep = "\t")
colnames(tcpmoh_1_count_data)[2] ="TCPMOH_1"
# tcpmoh2
tcpmoh_2_count_data <- read.csv('data/TCPMOH_2_counts.tsv', header = TRUE, sep = "\t")
colnames(tcpmoh_2_count_data)[2] ="TCPMOH_2"
Organize Data
To actually run DESeq2, we want all the gene count data in one data frame. Therefore, we need to merge the data sets based on the Geneid.
# put all data frames into list
df_list <- list(dmso_1_count_data, dmso_2_count_data, tcpmoh_1_count_data, tcpmoh_2_count_data)
# merge all data frames in list
count_data <- df_list %>% reduce(full_join, by='Geneid')
# make the Geneid the row id
count_data <- count_data %>% remove_rownames %>% column_to_rownames(var="Geneid")
head(count_data)
## DMSO_1 DMSO_2 TCPMOH_1 TCPMOH_2
## CYTB 136445 176209 137586 193845
## ND6 1443 2208 1317 1695
## ND5 77429 93110 74762 93454
## ND4 71269 94885 79189 100194
## ND4L 12122 14628 11453 14168
## ND3 10141 14687 10184 13766
We also need a meta data table. This can be done in R or by filling out a csv/tsv file and reading it in. I will just generate it here.
# get the sample names from the count_data matrix
SampleName <- c(colnames(count_data))
# specify the conditions for each sample
# In my sample names, DMSO_1 and DMSO_2 are control replicated and TCPMOH_1 and TCPMOH_2 are treated replicates
condition <- c("control", "control", "treated", "treated")
# generate the metadata data frame
meta_data <- data.frame(SampleName, condition)
# make the sample name the row id
meta_data <- meta_data %>% remove_rownames %>% column_to_rownames(var="SampleName")
meta_data
## condition
## DMSO_1 control
## DMSO_2 control
## TCPMOH_1 treated
## TCPMOH_2 treated
We want to double check that the name of the columns in the counts matrix is the same as the names in the meta data file. We can do this using the following two lines of code. Ideally, you should see TRUE and TRUE returned. If not, the rows in meta_data do not match the columns in counts_data and you should trouble shoot before moving forward. This is usually just a misorganization somewhere or a typo.
all(colnames(count_data) %in% rownames(meta_data))
## [1] TRUE
all(colnames(count_data) == rownames(meta_data))
## [1] TRUE
Differential Expression Analysis
Set Up DESeq Objects
Now that we have the data organized, we can run the differential expression analysis. The first step is to create a DESeq data object to complete the downstream analysis.
# create deseq data set object
dds <- DESeqDataSetFromMatrix(countData = count_data,
colData = meta_data,
design = ~ condition)
We can investigate the DESeq object:
# run the below to see the object elements
dds
## class: DESeqDataSet
## dim: 30421 4
## metadata(1): version
## assays(1): counts
## rownames(30421): CYTB ND6 ... cep97 rpl24
## rowData names(0):
## colnames(4): DMSO_1 DMSO_2 TCPMOH_1 TCPMOH_2
## colData names(1): condition
We see here that the dimension is 36351x4 meaning we have 36351 genes and 4 samples. The row names are the gene ids. Everything looks as it should! An optional step is to pre fliter some of the genes with low counts.
Pre-Filtering
While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.
# filter any counts less than 10
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds
## class: DESeqDataSet
## dim: 24356 4
## metadata(1): version
## assays(1): counts
## rownames(24356): CYTB ND6 ... cep97 rpl24
## rowData names(0):
## colnames(4): DMSO_1 DMSO_2 TCPMOH_1 TCPMOH_2
## colData names(1): condition
You can see that when we pre-filter, we end up with 26859 genes left. The only downside to this step is these genes will no longer be included in the final output file. This isn’t an issue for most cases but if you are like me and often aggregate results from multiple studies, this can make your life difficult as the number of genes will now be different. Anyway, pre-filtering is great if you are just running a ‘typical’ workflow.
Set Reference Factor
By default, R will choose a reference level for factors based on alphabetical order. Then, if you never tell the DESeq2 functions which level you want to compare against (e.g. which level represents the control group), the comparisons will be based on the alphabetical order of the levels. Here, we just want to be sure the level we are comparing against is truly the control using the relevel function.
# set the reference to be the control
dds$condition <- relevel(dds$condition, ref = 'control')
Getting Normalized Counts
Sometimes, we want the normalized counts as an output for downstream analysis. To obtain them, you can run the following commands.
# get normalized counts
dds <- estimateSizeFactors(dds)
normalized_counts <- counts(dds, normalized = TRUE)
# If you would like to write the normalized counts to a file, you can run the following command. Note that you need to specify the file path.
# write.table(normalized_counts, file=normalized_counts_data_file_path, sep = '\t', quote=F, col.names = NA)
Run Differential Expression
The standard differential expression analysis steps are wrapped into a single function, DESeq. Results tables are generated using the function results, which extracts a results table with log2 fold changes, p values and adjusted p values. With no additional arguments to results, the log2 fold change and Wald test p value will be for the last variable in the design formula, and if this is a factor, the comparison will be the last level of this variable over the reference level (see previous note on factor levels). However, the order of the variables of the design do not matter so long as the user specifies the comparison to build a results table for, using the name or contrast arguments of results.
dds <- DESeq(dds)
res <- results(dds)
res
## log2 fold change (MLE): condition treated vs control
## Wald test p-value: condition treated vs control
## DataFrame with 24356 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## CYTB 158703.52 0.0480368 0.105504 0.455308 0.6488875 0.962519
## ND6 1643.65 -0.2958785 0.155837 -1.898644 0.0576113 0.487636
## ND5 83990.25 -0.0468041 0.102839 -0.455123 0.6490211 0.962519
## ND4 85437.84 0.0876639 0.110877 0.790639 0.4291546 0.917258
## ND4L 12990.62 -0.0880488 0.107552 -0.818662 0.4129792 0.908968
## ... ... ... ... ... ... ...
## hikeshi 1064.229 0.0618543 0.146767 0.421445 0.673430 0.966905
## eed 875.412 -0.0960908 0.166265 -0.577937 0.563307 0.948176
## nfkbiz 948.290 0.0612779 0.172558 0.355116 0.722503 0.977389
## cep97 661.991 -0.0895011 0.212066 -0.422043 0.672993 0.966905
## rpl24 35114.027 0.1329105 0.123686 1.074584 0.282561 0.844292
We can also get a summary:
summary(res)
##
## out of 24356 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 201, 0.83%
## LFC < 0 (down) : 270, 1.1%
## outliers [1] : 0, 0%
## low counts [2] : 9916, 41%
## (mean count < 204)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
If you would like to know the exact number of adjusted p-values that are <0.01, you can run the following command:
sum(res$padj < 0.1, na.rm=TRUE)
## [1] 471
As you can see, the results function defaults to a p-value <0.1. For stronger statistical significance, we often want to look at those <0.05. This can be specified in the results function:
res05 <- results(dds, alpha=0.05)
summary(res05)
##
## out of 24356 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 148, 0.61%
## LFC < 0 (down) : 197, 0.81%
## outliers [1] : 0, 0%
## low counts [2] : 6611, 27%
## (mean count < 75)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
And then the number of genes with an adjusted p-value <0.05 are:
sum(res$padj < 0.05, na.rm=TRUE)
## [1] 348
As you can see, this number is less as an adjusted p-value <0.05 is a stronger statistical significance. At this point, it is useful to download the results into one location. We can do that by writing the results to a csv file.
# write results to file, notice you need to specify the file path.
# write.csv(as.data.frame(res), file=deseq_results_file_path)
Visualizing Results
At this point, we can visualize and make conclusions for our data. The first question is what is considered a differentially expressed gene? There are two things to consider: (1) the lof2FC value and (2) the adjusted p-value.
DESeq2 uses the so-called Benjamini-Hochberg (BH) adjustment; in brief, this method calculates for each gene an adjusted p value which answers the following question: if one called significant all genes with a p value less than or equal to this gene’s p value threshold, what would be the fraction of false positives (the false discovery rate, FDR) among them (in the sense of the calculation outlined above)? These values, called the BH-adjusted p values, are given in the column padj of the results object. Hence, if we consider a fraction of 10% false positives acceptable, we can consider all genes with an adjusted p value below 10%=0.1 as significant.
Mathematically speaking, any log2FC value different than zero with low variability. Typically speaking, you will see people identify a differentially expressed gene for log2FC cutoff value between 1-1.5 (positive or negative). Additionally, the adjusted p-value is typically chosen to be <0.05 but is sometimes chosen to be <0.1 (DESeq default) or even <0.01 for stronger statistical confidence. For this tutorial, I will choose genes with an adjusted p-value <0.05 and a |Log2FC|<1.
# convert results data to basic dataframe
data <- data.frame(res)
head(data)
## baseMean log2FoldChange lfcSE stat pvalue padj
## CYTB 158703.518 0.04803679 0.1055039 0.4553082 0.64888750 0.9625189
## ND6 1643.651 -0.29587855 0.1558368 -1.8986437 0.05761133 0.4876364
## ND5 83990.255 -0.04680413 0.1028385 -0.4551225 0.64902109 0.9625189
## ND4 85437.841 0.08766391 0.1108773 0.7906391 0.42915461 0.9172576
## ND4L 12990.616 -0.08804881 0.1075521 -0.8186621 0.41297922 0.9089675
## ND3 12017.929 -0.07758008 0.1247725 -0.6217724 0.53409151 0.9438768
Notice the column names! The columns we are typically interested in are log2FoldChange and padj.
PCA Plot
Portions of this section were taken from this tutorial. DESeq2 has a built-in function for plotting PCA plots, that uses ggplot2 under the hood. This is great because it saves us having to type out lines of code and having to fiddle with the different ggplot2 layers. In addition, it takes the rlog object as an input directly, hence saving us the trouble of extracting the relevant information from it.
The function plotPCA() requires two arguments as input: an rlog object and the intgroup (the column in our metadata that we are interested in).
rld <- rlog(dds)
plotPCA(rld)
By default the function uses the top 500 most variable genes. Ideally, we would want to see that our samples cluster together based on treatment group.
Volcano Plots
A volcano plot is a type of scatterplot that shows statistical significance (P value) versus magnitude of change (fold change). It enables quick visual identification of genes with large fold changes that are also statistically significant. These may be the most biologically significant genes. We can get as simple or as fancy as we would like when it comes to volcano plots!
Identify Up and Down Regulated Genes
# add an additional column that identifies a gene as unregulated, downregulated, or unchanged
# note the choice of pvalue and log2FoldChange cutoff.
data <- data %>%
mutate(
Expression = case_when(log2FoldChange >= log(1) & padj <= 0.05 ~ "Up-regulated",
log2FoldChange <= -log(1) & padj <= 0.05 ~ "Down-regulated",
TRUE ~ "Unchanged")
)
head(data)
## baseMean log2FoldChange lfcSE stat pvalue padj
## CYTB 158703.518 0.04803679 0.1055039 0.4553082 0.64888750 0.9625189
## ND6 1643.651 -0.29587855 0.1558368 -1.8986437 0.05761133 0.4876364
## ND5 83990.255 -0.04680413 0.1028385 -0.4551225 0.64902109 0.9625189
## ND4 85437.841 0.08766391 0.1108773 0.7906391 0.42915461 0.9172576
## ND4L 12990.616 -0.08804881 0.1075521 -0.8186621 0.41297922 0.9089675
## ND3 12017.929 -0.07758008 0.1247725 -0.6217724 0.53409151 0.9438768
## Expression
## CYTB Unchanged
## ND6 Unchanged
## ND5 Unchanged
## ND4 Unchanged
## ND4L Unchanged
## ND3 Unchanged
Get the Top 10 Most Variable Genes
top <- 10
# we are getting the top 10 up and down regulated genes by filtering the column Up-regulated and Down-regulated and sorting by the adjusted p-value.
top_genes <- bind_rows(
data %>%
filter(Expression == 'Up-regulated') %>%
arrange(padj, desc(abs(log2FoldChange))) %>%
head(top),
data %>%
filter(Expression == 'Down-regulated') %>%
arrange(padj, desc(abs(log2FoldChange))) %>%
head(top)
)
# create a datframe just holding the top 10 genes
Top_Hits = head(arrange(data,pvalue),10)
Top_Hits
## baseMean log2FoldChange lfcSE stat pvalue
## pdk2b 9748.7672 1.329178 0.1346152 9.873903 5.402079e-23
## si:ch211-132b12.7 4790.7498 -1.239026 0.1256370 -9.861952 6.085478e-23
## trim63a 622.7489 1.887092 0.1942775 9.713383 2.644124e-22
## eevs 2600.6421 -1.623744 0.1680947 -9.659699 4.471781e-22
## zgc:113054 3347.4297 -1.115578 0.1315768 -8.478527 2.280618e-17
## si:ch211-161h7.5 2632.1364 -1.424045 0.1735748 -8.204213 2.321063e-16
## hebp2 2600.6605 -1.098667 0.1359765 -8.079829 6.485740e-16
## ndrg1b 716.7809 1.302238 0.1659977 7.844918 4.332367e-15
## psme4a 1946.1987 1.279067 0.1718202 7.444217 9.752090e-14
## hsp70.3 817.8511 -1.287905 0.1735650 -7.420303 1.168528e-13
## padj Expression
## pdk2b 4.393715e-19 Up-regulated
## si:ch211-132b12.7 4.393715e-19 Down-regulated
## trim63a 1.272705e-18 Up-regulated
## eevs 1.614313e-18 Down-regulated
## zgc:113054 6.586426e-14 Down-regulated
## si:ch211-161h7.5 5.586024e-13 Down-regulated
## hebp2 1.337916e-12 Down-regulated
## ndrg1b 7.819922e-12 Up-regulated
## psme4a 1.564669e-10 Up-regulated
## hsp70.3 1.687355e-10 Down-regulated
Basic Volcano Plot
data$label = if_else(rownames(data) %in% rownames(Top_Hits), rownames(data), "")
# basic plot
p1 <- ggplot(data, aes(log2FoldChange, -log(pvalue,10))) + # -log10 conversion
geom_point( size = 2/5) +
xlab(expression("log"[2]*"FC")) +
ylab(expression("-log"[10]*"P-Value")) +
xlim(-4.5, 4.5)
p1
Basic Volcano Plot - Red Significant Genes
# basic plot with line + red for p < 0.05
p2 <- ggplot(data, aes(log2FoldChange, -log(pvalue,10))) + # -log10 conversion
geom_point(aes(color = Expression), size = 2/5) +
#geom_hline(yintercept= -log10(0.05), linetype="dashed", linewidth = .4) +
xlab(expression("log"[2]*"FC")) +
ylab(expression("-log"[10]*"P-Value")) +
scale_color_manual(values = c("firebrick3", "black", "firebrick3")) +
xlim(-4.5, 4.5) +
theme(legend.position = "none")
p2
# with labels for top 10 sig overall
library(ggrepel)
p3 <- ggplot(data, aes(log2FoldChange, -log(pvalue,10))) + # -log10 conversion
geom_point(aes(color = Expression), size = 2/5) +
# geom_hline(yintercept=-log10(0.05), linetype="dashed", linewidth = .4) +
xlab(expression("log"[2]*"FC")) +
ylab(expression("-log"[10]*"P-Value")) +
scale_color_manual(values = c("firebrick3", "black", "firebrick3")) +
xlim(-4.5, 4.5) +
theme(legend.position = "none") +
geom_text_repel(aes(label = label), size = 2.5)
p3
Up/Down Regulated Genes Identified
# plot with up/down
p4 <- ggplot(data, aes(log2FoldChange, -log(pvalue,10))) + # -log10 conversion
geom_point(aes(color = Expression), size = 2/5) +
#geom_hline(yintercept=-log10(0.05), linetype="dashed", linewidth = .4) +
xlab(expression("log"[2]*"FC")) +
ylab(expression("-log"[10]*"P-Value")) +
scale_color_manual(values = c("dodgerblue3", "black", "firebrick3")) +
xlim(-4.5, 4.5)
p4
Up/Down Regulated Genes Identified and Labels
p5 <- ggplot(data, aes(log2FoldChange, -log(pvalue,10))) + # -log10 conversion
geom_point(aes(color = Expression), size = 2/5) +
# geom_hline(yintercept=-log10(0.05), linetype="dashed", linewidth = .4) +
xlab(expression("log"[2]*"FC")) +
ylab(expression("-log"[10]*"P-Value")) +
scale_color_manual(values = c("dodgerblue3", "black", "firebrick3")) +
xlim(-4.5, 4.5) +
geom_text_repel(aes(label = label), size = 2.5)
p5
Conclusion
Hope this tutorial was helpful! There are many other things that you can accomplish while analyzing differentially expressed genes. In fact, the next logical step would be to perform a gene enrichment analysis. Tutorial coming soon…