Lab 5: Quantification with Salmon & Exploratory Analysis

Start an interactive R session

# Request an interactive session with enough RAM for R
srun --pty --time=01:00:00 --mem=16G --cpus-per-task=4 bash

# Load the R module
module load R

# Start R
R

Load required R packages

library(tximport)        # import Salmon/kallisto transcript-level estimates
library(DESeq2)          # differential expression and VST transformation
library(GenomicFeatures) # build transcript database from GTF
library(ggplot2)         # grammar-of-graphics plotting
library(pheatmap)        # annotated heatmaps

Build a transcript-to-gene mapping from the chr11 GTF

tximport needs a two-column data frame mapping each transcript ID (TXNAME) to its gene ID (GENEID) to summarise transcript-level Salmon estimates to gene level.

gtf_path <- "/ibex/user/$USER/workshop/references/chr11/GRCh38.chr11.gtf"

# Parse the GTF into a TxDb (transcript database) object
txdb    <- makeTxDbFromGFF(gtf_path, format = "gtf")

# Extract all transcript IDs as the key set
k       <- keys(txdb, keytype = "TXNAME")

# Build a two-column data frame: TXNAME → GENEID
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

head(tx2gene)
cat("Transcript-to-gene mappings:", nrow(tx2gene), "\n")

Import Salmon output with tximport

# Set base workshop directory using the current user's home path
WORKSHOP <- paste0("/ibex/user/", Sys.getenv("USER"), "/workshop")

# Define the six sample IDs matching the Salmon output directory names
sample_ids <- c(
    "KO_1_SRR10045016", "KO_2_SRR10045017", "KO_3_SRR10045018",
    "WT_1_SRR10045019", "WT_2_SRR10045020", "WT_3_SRR10045021"
)

# Build paths to each sample's quant.sf file
quant_files <- file.path(WORKSHOP, "counts/salmon", sample_ids, "quant.sf")
names(quant_files) <- sample_ids  # names are used as column names in the output matrix

# Stop early with a clear error if any file is missing
stopifnot(all(file.exists(quant_files)))

# Import and summarise from transcript to gene level using the tx2gene map
# ignoreTxVersion strips the version suffix (e.g. ENST00000123456.3 → ENST00000123456)
txi <- tximport(
    quant_files,
    type            = "salmon",
    tx2gene         = tx2gene,
    ignoreTxVersion = TRUE
)

cat("Genes imported:", nrow(txi$counts), "\n")
cat("Samples:", ncol(txi$counts), "\n")

Build sample metadata (coldata)

# Derive condition labels from sample names (KO_* → knockdown, WT_* → control)
condition <- ifelse(grepl("^KO", sample_ids), "knockdown", "control")

# Build a sample metadata table (coldata) required by DESeq2
# Levels set so "control" is the reference level in downstream contrasts
coldata <- data.frame(
    sample    = sample_ids,
    condition = factor(condition, levels = c("control", "knockdown")),
    row.names = sample_ids
)

print(coldata)

Create DESeqDataSet and apply variance-stabilising transformation

# Create a DESeqDataSet from tximport output, incorporating sample metadata
# design = ~ condition tells DESeq2 which variable to model
dds <- DESeqDataSetFromTximport(
    txi,
    colData = coldata,
    design  = ~ condition
)

# Remove genes with very low total counts — they carry no statistical power
# and slow down downstream steps; ≥10 counts across all samples is a common threshold
dds <- dds[rowSums(counts(dds)) >= 10, ]
cat("Genes after low-count filter:", nrow(dds), "\n")

# Apply variance-stabilising transformation (VST) for visualisation
# blind = TRUE means the transformation ignores the sample groups (fair QC)
vsd <- vst(dds, blind = TRUE)
cat("VST complete\n")

PCA plot

# Use DESeq2's built-in PCA function; returnData = TRUE gives us the coordinates
# instead of plotting directly so we can customise with ggplot2
pca_data <- plotPCA(vsd, intgroup = "condition", returnData = TRUE)

# Extract the percentage of variance explained by each PC for axis labels
pct_var  <- round(100 * attr(pca_data, "percentVar"), 1)

# Build a customised PCA scatter plot
p_pca <- ggplot(pca_data, aes(x = PC1, y = PC2, colour = condition, label = name)) +
    geom_point(size = 4) +                                         # one dot per sample
    geom_text(vjust = -0.8, hjust = 0.5, size = 3.5) +            # sample labels above dots
    scale_colour_manual(values = c(control = "#0d7377", knockdown = "#e07b4c")) +
    xlab(paste0("PC1: ", pct_var[1], "% variance")) +
    ylab(paste0("PC2: ", pct_var[2], "% variance")) +
    ggtitle("PCA — GSE136366  (VST, chr11)") +
    theme_bw(base_size = 13)

print(p_pca)

# Save to a PDF in the results directory
results_dir <- file.path(WORKSHOP, "results")
dir.create(results_dir, showWarnings = FALSE)  # create if it doesn't exist yet
ggsave(file.path(results_dir, "pca_plot_R.pdf"), p_pca, width = 7, height = 5)

Sample-distance heatmap

# Extract the VST-transformed count matrix (genes × samples)
vst_mat   <- assay(vsd)

# Transpose so samples are rows, then compute Euclidean distances between samples
dist_mat  <- dist(t(vst_mat))

# Convert to a symmetric data frame for pheatmap (requires a matrix/data frame input)
dist_df   <- as.data.frame(as.matrix(dist_mat))

# Build a one-column annotation data frame for the colour sidebar on rows and columns
ann_col <- data.frame(
    condition = coldata$condition,
    row.names = rownames(coldata)
)
# Define colours for the annotation sidebar to match the PCA plot palette
ann_colors <- list(condition = c(control = "#0d7377", knockdown = "#e07b4c"))

# Draw the heatmap and save directly to PDF via the filename argument
pheatmap(
    dist_df,
    annotation_col    = ann_col,           # colour sidebar on columns
    annotation_row    = ann_col,           # colour sidebar on rows
    annotation_colors = ann_colors,
    color             = colorRampPalette(c("#2166ac", "white"))(100),  # blue→white gradient
    display_numbers   = TRUE,              # show distance values in each cell
    number_format     = "%.1f",            # one decimal place
    main              = "Sample-to-sample distances (VST, chr11)",
    filename          = file.path(results_dir, "sample_distance_heatmap_R.pdf"),
    width = 8, height = 7
)