Bulk RNA-seq Analysis: From Raw Read Counts to Publication-Style Plots

Overview

This is Differential Gene Expression Analysis guideline for LS6014 Project based on:

For further detail, please check these resources.

By the end you will have:

Loaded a gene × sample count matrix
Run a DESeq2 differential expression (DE) test between two conditions
Normalized counts so samples are comparable
Explored sample structure with PCA and clustering
Produced the standard set of RNA-seq figures: PCA plot, volcano plot, MA plot, DE heatmap, and single-gene boxplots
Exported a tidy results table and figures for your report

As you will use DESeq2, 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

Let’s get started!!

Let’s begin by opening up RStudio and setting up a new project for this analysis.

Go to the File menu and select New Project.
In the New Project window, choose New Directory. Then, choose Empty Project. Name your new directory DEanalysis and then “Create the project as subdirectory of:” the Desktop (or location of your choice).
The new project should automatically open in RStudio.
Check whether or not you are in the correct working directory, using getwd().
Go to the File menu and select New File, then select R Script. Save the file as de_script.R.
Create 2 new folders: “data” and “results”.
Save the full counts matrix file in the data directory.

Now let’s install the packages:

Firstly (only in the first time):

install.packages(c("BiocManager", "tidyverse","RColorBrewer","pheatmap","ggrepel","cowplot","clusterProfiler"))

Then, install the below packages from Bioconductor - (only in the first time):

library(BiocManager)
install("DESeq2")
install("clusterProfiler")
install("DOSE")
install("org.Hs.eg.db")
install("org.Mm.eg.db")
install("pathview")
install("DEGreport")
install("tximport")
install("AnnotationHub")
install("ensembldb")
install("apeglm")

Load the packages one at a time. If you want to, you can run sessionInfo() afterwards.

# Load packages
library(DESeq2)
library(tidyverse)
library(RColorBrewer)
library(pheatmap)
library(ggrepel)
library(cowplot)
library(clusterProfiler)
library(DEGreport)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(DOSE)
library(pathview)
library(tximport)
library(AnnotationHub)
library(ensembldb)
library(apeglm)
library(dplyr)

Now load your data:

# Load in data
# Read in the raw read counts and transform to matrix
rawCounts_1 = read.delim('data/full_counts.tsv')

rawCounts = rawCounts_1[,-1]
rownames(rawCounts) = rawCounts_1[,1]

count_matrix = as.matrix(rawCounts)

DESeq2 needs sample information (metadata) for performing DGE analysis. Let’s create the sample information

# Create sample information/metadata
coldata = data.frame(
  sample = c( "Eoma_1", "Eoma_3", "Eoma_4", "Eoma_5", "CCOC_1", "CCOC_2", "CCOC_4","CCOC_5" ),
  condition = c( "benign_endometrioma", "benign_endometrioma",  "benign_endometrioma", "benign_endometrioma", "concurrent_endometriosis", "concurrent_endometriosis","concurrent_endometriosis","concurrent_endometriosis"), 
  row.names = "sample")

coldata$condition = as.factor(coldata$condition)

It is essential to have the name of the columns in the count matrix in the same order as that in name of the samples (rownames in coldata).

# Check all is TRUE
all(rownames(coldata) %in% colnames(count_matrix))

all(rownames(coldata) == colnames(count_matrix))

Construct DESeqDataSet for DGE analysis

# Construct DESeqDataSet for DGE analysis
dds = DESeqDataSetFromMatrix(countData = count_matrix, 
                             colData = coldata, 
                             design = ~ condition)

If you have additional feature data, it can be added to the DESeqDataSet by adding to the metadata columns of a newly constructed object. In our situation, it would be helpful to add a new column with the gene names (we only have now the gene identifiers)

In bulk RNA-seq data, the way you convert gene IDs to gene names depends on which gene identifiers you currently have:

Ensembl IDs (e.g., ENSG00000141510)
Entrez IDs (e.g., 7157)
RefSeq IDs (e.g., NM_000546)
Other organism-specific IDs

The most common approach in R is to use Bioconductor’s AnnotationDbi and organism annotation packages

# Check your ID type first
head(rownames(count_matrix))

# Convert Ensembl IDs to gene symbols (Human)
genes = rawCounts_1$Geneid


# Remove Ensembl version numbers
rownames(count_matrix) = sub("\\..*", "", rownames(count_matrix))

genes = rownames(count_matrix)

gene_annotation = mapIds(
  org.Hs.eg.db,
  keys = genes,
  column = "SYMBOL",
  keytype = "ENSEMBL",
  multiVals = "first")

ID = data.frame(genes,gene_annotation)

ID = ID %>%
  mutate(gene_ID = coalesce(gene_annotation,genes)) # replace values NA with the original FB name

mcols(dds) = DataFrame(mcols(dds),ID)

rownames(dds) = rowData(dds)$gene_ID

Pre-filter the genes which have low counts. Low count genes may not have sufficient evidence for differential gene expression. Furthermore, removing low count genes reduce the load of multiple hypothesis testing corrections.

We will remove the genes which have < 10 reads (this can vary based on research goal) in total across all the samples. Pre-filtering helps to remove genes that have very few mapped reads, reduces memory, and increases the speed of the DESeq2 analysis.

# Pre-filtering
dds = dds[rowSums(counts(dds)) >= 10,]

Select the reference level for condition comparisons.

The reference level can set using ref parameter. The comparisons of other conditions will be compared against this reference i.e, the log2 fold changes will be calculated based on ref value (infected/control). If this parameter is not set, comparisons will be based on alphabetical order of the levels.

# Set control condition as reference
dds$condition = relevel(dds$condition, ref = "benign_endometrioma")

Differential expression analysis

Performing differential gene expression analysis using DESeq2 is the process of comparing gene expression levels between different sample groups (such as control vs. treatment) using the DESeq2 R package.

DESeq2 statistically identifies genes that are significantly upregulated or downregulated while accounting for biological variability and differences in sequencing depth, helping researchers understand the molecular changes associated with a condition or treatment.

# Perform differential gene expression analysis
dds = DESeq(dds)

# See all comparisons
resultsNames(dds)

# **Get gene expression table**
res = results(dds)
res

Order gene expression table by adjusted p value (Benjamini-Hochberg FDR method):

# Order gene expression table based on adjusted p-value
res[order(res$padj),]

Export differential gene expression analysis table to CSV file:

# Export gene expression analysis table
write.csv(as.data.frame(res[order(res$padj),] ), file="condition_concurrent_endometriosis_vs_control_benign_endometrioma.csv")

Get summary of differential gene expression with adjusted p value cut-off at 0.05:

# Get data summary when p value cut-off at 0.05
summary(results(dds, alpha=0.05))

Raw counts are not directly comparable across samples: a sample sequenced twice as deeply will have roughly twice the counts for every gene even with identical biology, and a handful of very highly expressed genes can distort the apparent expression of everything else (a “composition” effect). Get normalized counts.

# Normalize counts
normalized_counts = counts(dds, normalized=TRUE)
head(normalized_counts)

Exploratory Data Analysis

Before trusting any DE result, check whether samples group the way biology predicts:

PCA on the most variable genes — samples should broadly separate along PC1/PC2 by condition. If they instead separate by something else (batch, sequencing lane, extraction date), that variable is a confounder to model.

### Transform counts for data visualization
rld = rlog(dds, blind=TRUE)

# The blind=TRUE argument results in a transformation unbiased to sample condition information
# The rlog function returns a DESeqTransform object, another type of DESeq-specific object.
library("ggplot2")

p = plotPCA(rld, intgroup="condition") 
p

p + #geom_text(aes(label=name)) +
  ggtitle("Principal Component Analysis") +
  labs(colour  = "Conditions") +
  geom_point(size = 5) +
  theme(aspect.ratio=1) 

# Save PCA plot
ggsave("PCA_plot.jpeg")

# Optional #
# Identify genes that contribute most to the PCs
rld_mat = assay(rld)
pca = prcomp(t(rld_mat))

PC1 = data.frame(sort(abs(pca$rotation[,"PC1"]), decreasing=TRUE)[1:50])
genes_PC1 = rownames(PC1)

PC2 = data.frame(sort(abs(pca$rotation[,"PC2"]), decreasing=TRUE)[1:50])
genes_PC2 = rownames(PC2)

Visualizing the results

Four plots cover almost everything you need to report an RNA-seq DE result:

Volcano plot — effect size (log2 fold change) vs. statistical confidence (-log10 adjusted p-value) for every gene.
MA plot — fold change vs. average expression level; a healthy MA plot is centered on zero fold-change across the full range of expression.
Heatmap of top DE genes — how the top hits behave across every individual sample, not just on average.
Single-gene boxplot — the most defensible way to show one specific gene of interest, with every replicate visible.

Volcano Plots

Volcano plots represent the results of a differential expression test. While DESeq2 has an integrated volcano plot, the packages EnhancedVolcano draws nicer and more customisable plots. It takes the results of DESeq2 as input.

# Convert to dataframe
volcano_data = as.data.frame(res) %>%
na.omit() %>%
mutate(
  significance = case_when(
    padj < 0.05 & log2FoldChange > 1 ~ "Upregulated",
    padj < 0.05 & log2FoldChange < -1 ~ "Downregulated",
    TRUE ~ "Not Significant"
  )
)

# Volcano plot
ggplot(volcano_data,
       aes(x = log2FoldChange,
           y = -log10(padj),
           color = significance)) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c(
  "Upregulated" = "red",
  "Downregulated" = "blue",
  "Not Significant" = "grey"
)) +
geom_vline(xintercept = c(-1, 1),
           linetype = "dashed",
           color = "black") +
geom_hline(yintercept = -log10(0.05),
           linetype = "dashed",
           color = "black") +
labs(
  title = "Volcano Plot",
  x = "Log2 Fold Change",
  y = "-Log10 Adjusted P-value"
) +
theme_minimal()


# Label the top significant genes
library(ggrepel)

top_genes = volcano_data %>%
  arrange(padj) %>%
  head(10)

ggplot(volcano_data,
       aes(log2FoldChange, -log10(padj), color = significance)) +
  geom_point(alpha = 0.7) +
  geom_text_repel(
    data = top_genes,
    aes(label = rownames(top_genes)),
    size = 3
  ) +
  theme_minimal()

MA plots

An MA plot displays the relationship between average gene expression and log₂ fold change between experimental conditions. Genes with significant differential expression appear above or below the zero line. Fold-change shrinkage was applied to reduce the influence of low-count genes with highly variable expression estimates, resulting in more accurate and interpretable effect sizes.

# Basic MA plot
plotMA(res,
       ylim = c(-5, 5),
       main = "DESeq2 MA Plot")


# MA plot with shrunk log2 fold changes (recommended)
# Shrink log2 fold changes
resLFC = lfcShrink(dds,
                    coef = 2,
                    type = "apeglm")

# MA plot
plotMA(resLFC,
       ylim = c(-5, 5),
       main = "MA Plot (Shrunken Log2FC)")

# Greater control over appearance
ma_data = as.data.frame(res) %>%
  na.omit() %>%
  mutate(
    significant = padj < 0.05
  )

ggplot(ma_data,
       aes(x = baseMean,
           y = log2FoldChange,
           color = significant)) +
  geom_point(alpha = 0.6, size = 1.5) +
  scale_x_log10() +
  scale_color_manual(values = c("grey", "red")) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(
    title = "MA Plot",
    x = "Mean Expression (baseMean)",
    y = "Log2 Fold Change"
  ) +
  theme_minimal()

Heatmaps

A heatmap was generated using variance-stabilized expression values obtained with DESeq2’s variance stabilizing transformation (VST). The top differentially expressed genes were selected based on adjusted p-values and visualized after row-wise scaling (z-score transformation). Hierarchical clustering was applied to both genes and samples to identify patterns of expression and sample similarity. The color scale represents relative expression levels, with red indicating higher expression and blue indicating lower expression relative to each gene’s mean expression across samples.

# Heatmap of the 50 top differentially expressed genes

# Variance stabilizing transformation
vsd = vst(dds, blind = FALSE)

# Select top 50 DE genes
top_genes = rownames(res[order(res$padj), ])[1:50]

# Extract transformed counts
mat = assay(vsd)[top_genes, ]

# Center each gene (row)
mat = t(scale(t(mat)))

# Sample information
annotation_col = as.data.frame(colData(dds)[, "condition", drop = FALSE])

# Heatmap
pheatmap(
mat,
annotation_col = annotation_col,
show_rownames = TRUE,
cluster_rows = TRUE,
cluster_cols = TRUE,
fontsize_row = 6,
main = "Top 50 Differentially Expressed Genes"
)

Single-gene expression plots

Visualizes the normalized expression levels of one gene across the experimental groups, allowing comparison of its expression between conditions and assessment of sample-to-sample variation. DESeq2’s plotCounts() function is commonly used for this purpose.

# Using DESeq2's built-in plotCounts()
d = plotCounts(
  dds,
  gene = "FGF7",
  intgroup = "condition",
  returnData = TRUE
)

ggplot(d, aes(condition, count, color = condition)) +
  geom_point(position = position_jitter(width = 0.1), size = 3) +
  stat_summary(fun = mean, geom = "crossbar",
               width = 0.5, color = "black") +
  scale_y_log10() +
  theme_bw() +
  labs(
    title = "TP53 expression",
    y = "Normalized counts"
  )

Functional enrichment

A ranked gene list is usually the input to a functional enrichment analysis, asking whether your DE genes concentrate in particular biological pathways rather than being a random subset of the genome. Gene Set Enrichment Analysis (GSEA) uses the entire ranked gene list, so it doesn’t require a significance cutoff first.

Dot plot of enrichment analysis results. Each dot represents an enriched term. Dot size corresponds to the number of genes associated with the term, while color indicates the adjusted p-value. Larger and darker-colored dots represent more significantly enriched terms supported by a greater number of genes.

Gene Set Enrichment Analysis (GSEA)

library(clusterProfiler)
library(org.Hs.eg.db)
library(dplyr)

# Create ranked gene list from DESeq2 results
res_df = as.data.frame(res) %>%
  na.omit()

gene_list = res_df$log2FoldChange
names(gene_list) = rownames(res_df)

# Convert gene symbols to Entrez IDs
gene_df = bitr(
  names(gene_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db
)

# Keep only mapped genes
gene_list = gene_list[gene_df$SYMBOL]
names(gene_list) = gene_df$ENTREZID

# Sort ranked list
gene_list = sort(gene_list, decreasing = TRUE)

# Run GSEA (GO Biological Process)
gsea_go = gseGO(
  geneList = gene_list,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  minGSSize = 10,
  maxGSSize = 500,
  pvalueCutoff = 0.05,
  verbose = FALSE
)

# View results
head(as.data.frame(gsea_go))

# Plot
dotplot(gsea_go, showCategory = 10)