Analyzing RNA-Seq with DESeq2

Code chunk options

The knitr package provides a lot of options for customizing nearly all components of code chunk, including the text outputs, source codes, and plots. Under the text outputs, you can find echo. When echo = TRUE this means to display the source code in the output documents. When warning = FALSE this means that all warnings will be suppressed.

Install the required R libraries

# Install CRAN packages
#install.packages(c("knitr", "RColorBrewer", "pheatmap"))

library(knitr)
library(DESeq2) 
library(RColorBrewer)
library(pheatmap)
library(ggplot2)
library(tidyverse)
library(patchwork)
library(ggplotify)
library(ggpubr)

Libraries loaded

Loading the packages required for subsequent analyses. In the future, we will be using pacman::p_load. Please look to have these packages installed for use by next week.

# loading required packages 
#pacman::p_load(DESeq2, dplyr, tidyr, ggplot2, genekitr, pheatmap, 
               #RColorBrewer, EnhancedVolcano, clusterProfiler, DOSE, 
               #ggnewscale, enrichplot, knitr, cowplot, msigdbr, pathview, 
               #patchwork, ggplotify)

Overview

To use DESeq2 it requires the following:

1. counts matrix

2. table of sample information

3. design indicates how to model the samples

Why un-normalized counts?

As input, the DESeq2 package expects count data to be in the form of a matrix of integer values. The values in the matrix should be un-normalized counts. The DESeq2 model internally corrects for library size, so transformed or normalized values should not be used as input. To normalize for sequencing depth and RNA composition, DESeq2 uses the median of ratios method. On the user-end there is only one step, but on the back-end there are multiple steps involved, including:

Step 1: creating a pseudo-reference sample (row-wise geometric mean)

Step 2: calculates ratio of each sample to the reference

• For every gene in a sample, the ratios (sample/ref) are calculated (as shown below). This is performed for each sample in the dataset. Since the majority of genes are not differentially expressed, the majority of genes in each sample should have similar ratios within the sample.

Step 3: calculate the normalization factor for each sample (size factor)

• The median value (column-wise for the above table) of all ratios for a given sample is taken as the normalization factor (size factor) for that sample, as calculated below.

Step 4: calculate the normalized count values using the normalization factor

The DESeqDataSet

The object used by the DESeq2 package to store the (1) read counts, (2) table of sample information , and (3) design is the DESeqDataSet, this is usually represented as dds.

dds = 1 + 2 + 3

dds = readcounts, table of sample, design

There are multiple ways of constructing a DESeqDataSet, depending on what pipeline was used upstream of DESeq2 to generated counts or estimated counts. For example, if you used Salmon to generate read-counts, you would use DESeqDataSetFromTximport()

1. From transcript abundance files and tximport - DESeqDataSetFromTximport()
2. From a count matrix - DESeqDataSetFromMatrix
3. From htseq-count files - DESeqDataSetFromHTSeqCount

We will be discussing how to input files output from htseq-count program.

htseq-count input

Let’s take a closer look at the contents of the ENSG_counts folder:

• img folder contains images
• counts folder contains outputs from HTSeq2 in the form of a .count file
• The metadata table called GSE164713_tcf1_anno.csv

To input our metadata file:

sampleTable <- read.csv("GSE164713_tcf1_anno.csv", header = T)
head(sampleTable)

##         samplename               filename   condition replicate
## 1      WT_CD8_rep1 SRR13423162.ENSG.count      WT_CD8         1
## 2      WT_CD8_rep2 SRR13423163.ENSG.count      WT_CD8         2
## 3      WT_CD8_rep3 SRR13423164.ENSG.count      WT_CD8         3
## 4 TCF1_KO_CD8_rep1 SRR13423165.ENSG.count TCF1_KO_CD8         1
## 5 TCF1_KO_CD8_rep2 SRR13423166.ENSG.count TCF1_KO_CD8         2
## 6 TCF1_KO_CD8_rep3 SRR13423167.ENSG.count TCF1_KO_CD8         3

The order here actually matters. Let’s check why by turning to the DESeqDataSetFromHTSeqCount manual:

#?DESeqDataSetFromHTSeqCount

#Convert to data frame 
sampleTable <- data.frame(sampleTable)
head(sampleTable)

##         samplename               filename   condition replicate
## 1      WT_CD8_rep1 SRR13423162.ENSG.count      WT_CD8         1
## 2      WT_CD8_rep2 SRR13423163.ENSG.count      WT_CD8         2
## 3      WT_CD8_rep3 SRR13423164.ENSG.count      WT_CD8         3
## 4 TCF1_KO_CD8_rep1 SRR13423165.ENSG.count TCF1_KO_CD8         1
## 5 TCF1_KO_CD8_rep2 SRR13423166.ENSG.count TCF1_KO_CD8         2
## 6 TCF1_KO_CD8_rep3 SRR13423167.ENSG.count TCF1_KO_CD8         3

str(sampleTable)

## 'data.frame':    6 obs. of  4 variables:
##  $ samplename: chr  "WT_CD8_rep1" "WT_CD8_rep2" "WT_CD8_rep3" "TCF1_KO_CD8_rep1" ...
##  $ filename  : chr  "SRR13423162.ENSG.count" "SRR13423163.ENSG.count" "SRR13423164.ENSG.count" "SRR13423165.ENSG.count" ...
##  $ condition : chr  "WT_CD8" "WT_CD8" "WT_CD8" "TCF1_KO_CD8" ...
##  $ replicate : int  1 2 3 1 2 3

Next, we want to specify the condition column as a factor. By doing so, this will allow us to categorize the data and store it as a level.

#specify as factor 
condition <- factor(sampleTable$condition)

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.

Below we perform a minimal pre-filtering to keep only rows that have at least 10 reads total.

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

Note on factor levels

By default, R will choose a reference level for factors based on alphabetical order. Therefore, 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.

use relevel to specify the reference level:

dds$condition <- relevel(dds$condition, ref = "WT_CD8")

Library size differences

With DESeq2, size factors are calculated using the estimateSizeFactors() function. The size factors estimated by this function combines an adjustment for differences in library sizes with an adjustment for differences in the RNA composition of the samples. The latter is important due to the compositional nature of RNA-seq data. There is a fixed number of reads to distribute between the genes, and if a single (or a few) very highly expressed gene consume a large part of the reads, all other genes will consequently receive very low counts.

dds <- estimateSizeFactors(dds)

# Plot the size factors against library size
# and look for any patterns by group:

ggplot(data.frame(libSize = colSums(assay(dds)),
                  sizeFactor = sizeFactors(dds),
                  Group = dds$condition),
       aes(x = libSize, y = sizeFactor, col = condition)) + 
  geom_point(size = 5) + theme_bw() + 
  labs(x = "Library size", y = "Size factor")

Principal components analysis (PCA)

In order to test for differential analysis, we used raw counts. These matrices are so large that we need convenient ways to extract the important information from these large data matrices therefore it is useful to transform the data.

In using principal component analysis (PCA):

• The data is reduced to smaller matrices so they can more easily be examined, plotted and interpreted.
• The most important factors are extracted (principal components).

Below we are using an variance stabilizing transformation (vst) (Tibshirani 1988; Huber et al. 2003; Anders and Huber 2010) or regularized logarithm (rlog) (Love, Huber, and Anders 2014). Both transformations will produce transformed data on the log2 scale which has been normalized with respect to library size or other normalization factors.

Using the transformed data we will generate a few outputs.

• PCA plot: The PCA plot will show the samples in a 2D plane spanned by the first two principal components. This type of plot is useful for visualizing the overall effect of experimental covariates (ex. age) and batch effects (ex. day sample was prepared).

vstcounts <- vst(dds, blind = FALSE)

It is also possible to customize the PCA plot using the ggplot function. Below we are adding shapes to represent the different days the samples were created.

pcaData <- plotPCA(vstcounts, intgroup=c("condition"), returnData=TRUE)

percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color=condition)) +
  geom_point(size=2) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() 

Now, we are making sure the x and y axis are similar using ylim argument

#pdf(file = "output/PCA_plot_TCF1_RNAseq.pdf", height = 6, width = 6)

pcaData <- plotPCA(vstcounts, intgroup=c("condition"), returnData=TRUE)

percentVar <- round(100 * attr(pcaData, "percentVar"))

pca <- ggplot(pcaData, aes(PC1, PC2, color=condition, shape = condition)) +
  geom_point(size=4, alpha=0.8) + # adding transparency
  ylim(-10,10) + 
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() + 
  ggtitle("PCA of TCF1 RNASeq data") + 
  scale_color_brewer(palette = "Set1") + 
  theme_classic(base_size = 10) + 
  theme(legend.position = "right", legend.title = element_blank(),
        plot.title = element_text(hjust=0.5, face = "bold"))

#dev.off()

#saveRDS(pca, file = "pca.RDS")

Now that we got it looking the way we want it, we can take some time to interpret the plot. Each data point represents one sample. The x-axis typically PC1, explains the largest source of variation in the dataset. We can see that PC1 is primarily explained by knocking out the gene Tcf7 (encodes Tcf1) in CD8 T cells. The next largest source of variance within the dataset, or PC2, is explained by one sample created on day1 in the WT group. Overall, the groups are clustering as expected.

Finally, we can generate a correlation heatmap. This is a type of plot used to visualize the strength of relationships between samples.

#pdf(file = "output/correlation_plot_TCF1_RNAseq.pdf", height = 3, width = 4)

vst_mat <- assay(vstcounts)
vst_cor <- cor(vst_mat)
annotation_heatmap <- data.frame(row.names=colnames(vst_mat), condition)

cor_heatmap <- as.ggplot(pheatmap(vst_cor,
                                  annotation_col = annotation_heatmap,
                                  main = "sample correlation", 
                                  show_rownames = F, 
                                  show_colnames = F,
                                  fontsize = 6,  # Decrease overall font size
                                  fontsize_col = 6,  
                                  fontsize_row = 6,
                                  width = 3,
                                  height = 4))  # Decrease legend size

#dev.off()

saveRDS(cor_heatmap, "cor_heatmap.RDS")

library(ggplot2)
library(pheatmap)
library(grid)
library(cowplot)
library(magick)  


# Save heatmap to a PDF with controlled size
pdf("temp_heatmap.pdf", width = 4, height = 2)

pheatmap(vst_cor,
         annotation_col = annotation_heatmap,
         show_rownames = FALSE, 
         show_colnames = FALSE,
         fontsize = 8, 
         fontsize_col = 6, 
         fontsize_row = 6,
         legend_cex = 0.6)
dev.off()

# Convert PDF to image
heatmap_img <- image_read_pdf("temp_heatmap.pdf")

# Convert image to grob
heatmap_grob <- rasterGrob(as.raster(heatmap_img), interpolate = TRUE)

# Combine PCA plot and heatmap, ensuring heatmap is half the height
FigAB <- plot_grid(pca, 
                   heatmap_grob, 
                   labels = c("(A)", "(B)"), 
                   ncol = 2,
                   rel_widths = c(1, 1),    
                   rel_heights = c(1, 0.5), # Heatmap is half the height
                   align = "hv")  

# Save final figure
ggsave("QC_figure.pdf", FigAB, height = 6, width = 8, dpi = 300, units = "in", scale = 1)

Session info

sessionInfo()

## R version 4.4.3 (2025-02-28)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggpubr_0.6.0                ggplotify_0.1.2            
##  [3] patchwork_1.3.0             lubridate_1.9.4            
##  [5] forcats_1.0.0               stringr_1.5.1              
##  [7] dplyr_1.1.4                 purrr_1.0.4                
##  [9] readr_2.1.5                 tidyr_1.3.1                
## [11] tibble_3.2.1                tidyverse_2.0.0            
## [13] ggplot2_3.5.1               pheatmap_1.0.12            
## [15] RColorBrewer_1.1-3          DESeq2_1.46.0              
## [17] SummarizedExperiment_1.36.0 Biobase_2.66.0             
## [19] MatrixGenerics_1.18.1       matrixStats_1.5.0          
## [21] GenomicRanges_1.58.0        GenomeInfoDb_1.42.3        
## [23] IRanges_2.40.1              S4Vectors_0.44.0           
## [25] BiocGenerics_0.52.0         knitr_1.50                 
## 
## loaded via a namespace (and not attached):
##  [1] rlang_1.1.5             magrittr_2.0.3          compiler_4.4.3         
##  [4] vctrs_0.6.5             pkgconfig_2.0.3         crayon_1.5.3           
##  [7] fastmap_1.2.0           backports_1.5.0         XVector_0.46.0         
## [10] rmdformats_1.0.4        labeling_0.4.3          rmarkdown_2.29         
## [13] tzdb_0.5.0              UCSC.utils_1.2.0        xfun_0.51              
## [16] zlibbioc_1.52.0         cachem_1.1.0            jsonlite_2.0.0         
## [19] DelayedArray_0.32.0     BiocParallel_1.40.0     broom_1.0.8            
## [22] parallel_4.4.3          R6_2.6.1                bslib_0.9.0            
## [25] stringi_1.8.7           car_3.1-3               jquerylib_0.1.4        
## [28] numDeriv_2016.8-1.1     Rcpp_1.0.14             bookdown_0.42          
## [31] Matrix_1.7-3            timechange_0.3.0        tidyselect_1.2.1       
## [34] rstudioapi_0.17.1       abind_1.4-8             yaml_2.3.10            
## [37] codetools_0.2-20        lattice_0.22-6          plyr_1.8.9             
## [40] withr_3.0.2             coda_0.19-4.1           evaluate_1.0.3         
## [43] gridGraphics_0.5-1      pillar_1.10.1           carData_3.0-5          
## [46] generics_0.1.3          emdbook_1.3.13          hms_1.1.3              
## [49] munsell_0.5.1           scales_1.3.0            glue_1.8.0             
## [52] tools_4.4.3             apeglm_1.28.0           locfit_1.5-9.12        
## [55] ggsignif_0.6.4          fs_1.6.5                mvtnorm_1.3-3          
## [58] grid_4.4.3              bbmle_1.0.25.1          bdsmatrix_1.3-7        
## [61] colorspace_2.1-1        GenomeInfoDbData_1.2.13 Formula_1.2-5          
## [64] cli_3.6.4               S4Arrays_1.6.0          gtable_0.3.6           
## [67] rstatix_0.7.2           yulab.utils_0.2.0       sass_0.4.9             
## [70] digest_0.6.37           SparseArray_1.6.2       farver_2.1.2           
## [73] htmltools_0.5.8.1       lifecycle_1.0.4         httr_1.4.7             
## [76] MASS_7.3-65

Citation

title: “Analyzing RNA-Seq with DESeq2” author: “Michael I. Love, Simon Anders, and Wolfgang Huber”

Note: if you use DESeq2 in published research, please cite:

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8