Skip to main content

STEP 7: Analysis

Downstream analysis workflow

For the following steps, a fully worked example notebook is available on the next page with all outputs rendered. You can also download the original .Rmd from the GitHub repository to run it in RStudio on your own count matrix.

View worked example notebook | Download .Rmd from GitHub

It contains the following steps:

  • Quality control of the count matrix - PCA, to determine if samples cluster by condition.
  • Differential expression analysis - Using DESeq2 to identify differentially expressed genes.
  • Visualization of results - MAplot to visualize the results of the differential expression analysis.
  • Targeted analysis genes - Pulling out specific genes of interest to visualize their expression across conditions.
  • Exporting up and down regulated genes - Exporting the results of the differential expression analysis to a .bed file for further analysis.

However below is a brief/simplified overview of the steps involved in the analysis.

Differential Expression Analysis with DESeq2

DESeq2 is the standard tool for RNA-seq differential expression analysis (runs in R).

Basic DESeq2 Workflow

# Load libraries
library(DESeq2)

# Read count matrix
countData <- read.table("counts.txt", header=TRUE, row.names=1)

# Create sample metadata
colData <- data.frame(
  condition = c("control", "control", "treatment", "treatment"),
  row.names = colnames(countData))

# Create DESeq2 dataset
dds <- DESeqDataSetFromMatrix(
  countData = countData,
  colData = colData,
  design = ~ condition)

# Run differential expression analysis
dds <- DESeq(dds)

# Extract results
results <- results(dds, contrast=c("condition", "treatment", "control"))

# Get significant genes (adjusted p-value < 0.05)
sig_genes <- subset(results, padj < 0.05)

Quality Control Plots

# PCA plot
vsd <- vst(dds, blind=FALSE)
plotPCA(vsd, intgroup="condition")

# MA plot
plotMA(results, ylim=c(-2,2))

# Volcano plot
plot(results$log2FoldChange, -log10(results$padj),
     xlab="log2 Fold Change", ylab="-log10 padj")

Alternative Tools

  • edgeR - Another R package for differential expression
  • limma-voom - Versatile linear modeling approach
  • sleuth - For use with Kallisto/Salmon quantification

Typical Outputs

  • Count matrix - Gene expression levels per sample
  • Differential expression results - Log2 fold changes, p-values
  • QC plots - PCA, heatmaps, MA plots
  • Gene lists - Upregulated/downregulated genes

Expression Visualizations in R

Here are some other pointers for visualization steps that can be performed on the count matrix.

Sample Correlation Heatmap

library(DESeq2)
library(pheatmap)

# Calculate sample-to-sample distances
vsd <- vst(dds, blind=FALSE)
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)

# Plot heatmap
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists)

PCA Plot

Principal component analysis to visualize sample relationships:

library(DESeq2)
library(ggplot2)

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

# PCA plot
plotPCA(vsd, intgroup="condition") +
  theme_classic() +
  ggtitle("PCA of RNA-seq Samples")

MA Plot

Show relationship between mean expression and fold change:

# After running DESeq2
plotMA(results, ylim=c(-5,5), main="MA Plot")

Volcano Plot

Visualize differential expression results:

library(ggplot2)

# Create data frame from results
df <- as.data.frame(results)
df$significant <- df$padj < 0.05 & abs(df$log2FoldChange) > 1

# Plot
ggplot(df, aes(x=log2FoldChange, y=-log10(padj), color=significant)) +
  geom_point(alpha=0.5) +
  scale_color_manual(values=c("grey", "red")) +
  theme_classic() +
  labs(title="Volcano Plot",
       x="log2 Fold Change",
       y="-log10 Adjusted P-value") +
  geom_hline(yintercept=-log10(0.05), linetype="dashed") +
  geom_vline(xintercept=c(-1,1), linetype="dashed")

Gene-specific Visualizations

Expression Bar Plots

Show expression of specific genes across conditions:

# Get normalized counts
counts_norm <- counts(dds, normalized=TRUE)

# Plot specific gene
gene_of_interest <- "Actb"
barplot(counts_norm[gene_of_interest,],
        main=gene_of_interest,
        ylab="Normalized Counts",
        las=2)

Gene Expression Box Plots

library(ggplot2)

# Extract gene counts
gene_counts <- plotCounts(dds, gene="ENSMUSG00000000001", 
                          intgroup="condition", returnData=TRUE)

# Box plot
ggplot(gene_counts, aes(x=condition, y=count)) +
  geom_boxplot() +
  geom_jitter(width=0.1) +
  scale_y_log10() +
  theme_classic() +
  labs(title="Gene Expression")

Heatmap of Top Differentially Expressed Genes

library(pheatmap)

# Get top 50 DE genes
top_genes <- head(order(results$padj), 50)

# Get normalized counts
mat <- assay(vsd)[top_genes,]

# Scale by row (gene)
mat_scaled <- t(scale(t(mat)))

# Plot heatmap
pheatmap(mat_scaled,
         cluster_rows=TRUE,
         cluster_cols=TRUE,
         show_rownames=TRUE,
         annotation_col=colData(dds)[,"condition", drop=FALSE],
         main="Top 50 Differentially Expressed Genes")

Best Practices

  1. Use biological replicates - At least 3 per condition
  2. Pre-filter low counts - Remove genes with very low counts
  3. Check QC plots - PCA, sample distances before analysis
  4. Use appropriate thresholds - padj < 0.05, |log2FC| > 1
  5. Validate key findings - qPCR or other experimental methods
  6. Consider batch effects - Account for technical variation