Task Overview

We are going to re-analyse a dataset containing transcriptomic samples from healthy patients and patients infected with mycobacterium tuberculosis, which causes tuberculosis disease. This data has previously been analysed by Singhania et al, 2018. Tuberculosis is a heterogeneous disease with high global burden. Patients infected with mycobacterium tuberculosis can be classified as active or latent infections, depending on whether clinical symptoms are present or not. Diagnosing tuberculosis infection can be difficult, especially in asymptomatic patients with latent tuberculosis infection (LTBI).

Therefore, it would be beneficial to identify some biomarkers which can act as sensitive and specific indicators of the disease state. The data is open-access and available through the National Center for Biotechnology Information’s Gene Expression Omnibus project. Data were collected from three cohorts (the London, South Africa, and Leicester cohorts).

The goal of this practical will be to perform an analysis of these data, using the principles learned in the lecture. Specifically, we will practice data pre-processing, normalisation, batch effect correction, data exploration with principal component analysis and hierarchical clustering, and differential gene expression analysis.

Data Management

0a) Set up your project

Why do we need a project?

When you analyse data in R, you need to tell R where files are. For example, you might load data into your environment by defining a path "C:/Users/me/Documents/data/file.xlsx". However, this is problematic, since it is specific to your own computer setup. The code will then break as soon as you move the files, share them with a colleague, or re-run the analysis on a different machine.

To avoid this, we will introduce some basic organisational skills for an R project. Firstly, we will contain the work in a single R project folder (the root), and refer to every file relative to that folder. This makes specifying paths for loading data and saving outputs more robust.

Using a project folder also helps to organise files in a clean way.

• raw data : the original data files, which you should never modify;
• analysis : the code (such as this .Rmd file) that reads the raw data and produces results;
• output : everything the code creates (figures, tables, processed data), which should be able to be perfectly regenerated by re-running the code.

Create the project folder

1. In RStudio, click File > New Project... > New Directory > New Project.
2. Give the project a name, for example omics_practicals, and choose where to save it on your computer. Click Create Project.
3. RStudio restarts inside your new project. You can check this in the top-right corner of RStudio, which now shows the project name, and in the Files pane, which shows a single file called e.g. omics_practicals.Rproj.

Next, create three sub-folders inside the project. In the Files pane of RStudio, click New Folder and create data-raw, analysis, and output. Alternatively, run the following code once in the R console.

dir.create("data-raw")
dir.create("analysis")
dir.create("output")

Your project should now look like this:

omics_practicals/
├── omics_practicals.Rproj    <- marks the "root" of the project
├── data-raw/                 <- the 3 raw data files
├── analysis/                 <- this practical (.Rmd file)
└── output/                   <- figures and results that you save

Add the files

1. Copy the three Excel files I sent you into the data-raw folder: GSE107991_Raw_counts_Berry_London.xlsx, GSE107992_Raw_counts_Berry_SouthAfrica.xlsx, and GSE107993_Raw_counts_Leicester_non_progressor_longitudnal_only.xlsx.
2. Copy this practical (the .Rmd file) into the analysis folder.
3. In the Files pane of RStudio, open the analysis folder and click on the .Rmd file to open it.

Important: always start your work by opening the project (double-click omics_practicals.Rproj, or use File > Open Project), and open the .Rmd file from within RStudio’s Files pane. If you open the .Rmd file directly from your Downloads folder, it will not be part of the project and the paths below will not work.

How do we refer to files?

By default, R Markdown runs code from the folder where the .Rmd file is saved (here, analysis/), which would make paths confusing. Instead, we will use the here package. The function here() always returns the path to the root of your project (the folder containing the .Rproj file), wherever the code is run from. We then build paths from this root. For example, here("data-raw", "my_file.xlsx") points to my_file.xlsx inside the data-raw folder, and here("output", "my_plot.pdf") points to a file in the output folder. The same code works on any computer and operating system, and you never have to type a full path.

0b) Install and Load Required Packages

Using install.packages(), install the following packages from the CRAN : BiocManager, readxl, dplyr, flashClust, UpSetR, tidyr, ggplot2, stringr, dendextend, ggdendro, ggrepel, VennDiagram, grid, here.

Using BiocManager::install(), install the following packages from BioConductor : GEOquery, edgeR, limma, pvca, dearseq, sva.

The following chunk will check if you already have any of these packages and then install the missing ones.

cran_packages <- c(
  "BiocManager",
  "readxl",
  "dplyr",
  "flashClust",
  "UpSetR",
  "tidyr",
  "ggplot2",
  "stringr",
  "dendextend",
  "ggdendro",
  "ggrepel",
  "VennDiagram",
  "grid",
  "here"
)

# Install missing CRAN packages
install.packages(setdiff(cran_packages, installed.packages()[, "Package"]))

# ---- BioConductor Packages ----
bioc_packages <-
  c("GEOquery", "edgeR", "limma", "pvca", "dearseq", "sva")

# Load BiocManager and install missing Bioconductor packages
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install(setdiff(bioc_packages, installed.packages()[, "Package"]))

0c) Check your project setup

Run the following chunk to check that everything is in the right place. The first line should print the path to your project folder (the one containing omics_practicals.Rproj), and all three files should be found.

# Path to the root of the project
here()

## [1] "/home/ah3/Desktop/Work/Teaching/PHDS_omics_data_course"

# The three raw data files that should be in the data-raw folder
raw_files <- c(
  "GSE107991_Raw_counts_Berry_London.xlsx",
  "GSE107992_Raw_counts_Berry_SouthAfrica.xlsx",
  "GSE107993_Raw_counts_Leicester_non_progressor_longitudnal_only.xlsx"
)

# Check that each file exists
data.frame(file = raw_files,
           found = file.exists(here("data-raw", raw_files)))

##                                                                  file found
## 1                              GSE107991_Raw_counts_Berry_London.xlsx  TRUE
## 2                         GSE107992_Raw_counts_Berry_SouthAfrica.xlsx  TRUE
## 3 GSE107993_Raw_counts_Leicester_non_progressor_longitudnal_only.xlsx  TRUE

If any file is not found (FALSE), or if the path printed above is not your project folder, go back to the previous steps: check the file names, that the files are in data-raw, and that you opened the project and then this .Rmd file from the Files pane.

The data come from the NCBI’s Gene Expression Omnibus. Here is the link for the series A modular transcriptional signature identifies phenotypic heterogeneity of human tuberculosis infection, which consists of the SubSeries GSE107991, GSE107992, and GSE107993. These SubSeries correspond to three different patient cohorts, called the London, South Africa, and Leicester cohorts, respectively. You do not need to download anything yourself: the raw count data are the files in your data-raw folder, and we will download the metadata directly in R below.

In this first step, I will give you the code to load the data and merge it together to make it easier to work with.

0d) Load data

Load the raw counts data with the read_excel() function, using here() to build the path to each file.

london_raw <-
  read_excel(here("data-raw", "GSE107991_Raw_counts_Berry_London.xlsx"))

southafrica_raw <-
  read_excel(here("data-raw", "GSE107992_Raw_counts_Berry_SouthAfrica.xlsx"))

leicester_raw <-
  read_excel(here(
    "data-raw",
    "GSE107993_Raw_counts_Leicester_non_progressor_longitudnal_only.xlsx"
  ))

0e) Download metadata

Use the getGEO() function from the GEOquery package to download the metadata corresponding to these files.

geo_london <- getGEO("GSE107991", GSEMatrix = FALSE)
geo_southafrica <- getGEO("GSE107992", GSEMatrix = FALSE)
geo_leicester <- getGEO("GSE107993", GSEMatrix = FALSE)

0f) Process the metadata

Extract the relevant information (cohort, location, set, disease status and sample name) from the metadata of each cohort.

# Extract the information of sample i from the GEO metadata, as a vector
extract_sample_info <- function(i, geo) {
  sample_label <- geo@gsms[[i]]@header$source_name_ch1

  # Remove underscores within values, so that splitting on "_" separates the fields
  sample_label <- gsub("Active_TB", "ActiveTB", sample_label)
  sample_label <- gsub("Test_set", "TestSet", sample_label)
  sample_label <- gsub("Validation_set", "ValidationSet", sample_label)
  sample_label <- gsub("Non_progressor", "NonProgressor", sample_label)

  c(unlist(strsplit(sample_label, split = "_")), geo@gsms[[i]]@header$title)
}

# Build a data frame of the sample information (one row per sample) of a cohort
build_cohort_metadata <- function(geo) {
  sample_info <- sapply(seq_along(geo@gsms), FUN = extract_sample_info, geo = geo)
  cohort_metadata <- cbind.data.frame(GSMID = names(geo@gsms), t(sample_info))
  rownames(cohort_metadata) <- names(geo@gsms)
  cohort_metadata
}

metadata_london <- build_cohort_metadata(geo_london)
colnames(metadata_london)[-1] <-
  c("Cohort", "Location", "Set", "Status", "SampleName")

metadata_southafrica <- build_cohort_metadata(geo_southafrica)
colnames(metadata_southafrica)[-1] <-
  c("Cohort", "Location", "Set", "Status", "SampleName")

# Note : in the Leicester metadata, the order of Status and Set is reversed
metadata_leicester <- build_cohort_metadata(geo_leicester)
colnames(metadata_leicester)[-1] <-
  c("Cohort", "Location", "Status", "Set", "SampleName")

0g) Merge data

Merge the three cohorts. We also change the format of the counts to wide, so that each row is a sample and each column is a gene.

# Stack the metadata of the three cohorts (columns are matched by name)
metadata_cohorts <-
  rbind.data.frame(metadata_london, metadata_southafrica, metadata_leicester)

# Merge the counts of the three cohorts. The first 3 columns identify the genes
# (Genes, Gene_name, Gene_biotype) and are used to match them.
counts_by_gene <-
  merge(
    x = merge(x = london_raw, y = southafrica_raw, by = 1:3),
    y = leicester_raw,
    by = 1:3
  ) %>%
  dplyr::select(-Gene_name, -Gene_biotype)

# Transpose, so that rows are samples and columns are genes
# (the first row, which contains the gene names, is dropped)
counts_wide <- t(counts_by_gene)[-1, ] %>%
  as.data.frame()
colnames(counts_wide) <- counts_by_gene$Genes
rownames(counts_wide) <- seq_len(nrow(counts_wide))

# Merge the counts with the metadata by sample name
counts_metadata <- merge(
  x = metadata_cohorts,
  y = bind_cols(SampleName = colnames(counts_by_gene)[-1], counts_wide),
  by = "SampleName"
) %>%
  distinct()

# Separate the counts (stored as numbers) from the metadata
metadata_cols <- c("SampleName", "GSMID", "Cohort", "Location", "Set", "Status")

counts_raw <- counts_metadata %>%
  dplyr::select(-all_of(metadata_cols)) %>%
  sapply(as.numeric) %>%
  as.data.frame()

metadata_all <- counts_metadata %>%
  dplyr::select(all_of(metadata_cols))

# Visualise the start of the merged data
counts_metadata[, 1:8] %>%
  head(5)

##              SampleName      GSMID Cohort Location     Set   Status
## 1  Berry_London_Sample1 GSM2886035  Berry   London TestSet ActiveTB
## 2 Berry_London_Sample10 GSM2886044  Berry   London TestSet ActiveTB
## 3 Berry_London_Sample11 GSM2886045  Berry   London TestSet     LTBI
## 4 Berry_London_Sample12 GSM2886046  Berry   London TestSet     LTBI
## 5 Berry_London_Sample13 GSM2886047  Berry   London TestSet     LTBI
##   ENSG00000000003 ENSG00000000005
## 1               5               0
## 2               7               0
## 3               2               0
## 4              14               0
## 5              26               1

1) Data pre-processing

1a) Remove unexpressed genes

A typical pre-processing step is to remove genes which are unexpressed (i.e. which have 0 counts in every sample), since these provide no information. Remove these genes from the data.

How many genes do you remove with this step?

# Total counts of each gene across all samples
gene_total_counts <- colSums(counts_raw)

# Unexpressed genes have no counts in any sample
unexpressed_genes <- colnames(counts_raw)[gene_total_counts == 0]

counts_expressed <- counts_raw %>%
  dplyr::select(-all_of(unexpressed_genes))

message(length(unexpressed_genes), " unexpressed genes removed from the data.")

## 6443 unexpressed genes removed from the data.

Answer: 6,443 genes were removed from the data.

1b) Remove non-protein coding genes

Many investigators only include protein-coding genes in their analysis. Information on whether a gene is protein coding or not may be found in the raw counts under the Gene_biotype column. For example, the first biotype values for the London data are:

london_raw$Gene_biotype %>%
  head()

## [1] "protein_coding" "protein_coding" "protein_coding" "protein_coding"
## [5] "protein_coding" "protein_coding"

Using this information, filter the counts to only include genes that are protein-coding. How many genes were removed in this stage? Store the filtered counts in an object called counts_protein_coding, which is used in the following sections.

# Names of the non-protein-coding genes (we use the London data to find them)
non_protein_coding_genes <- london_raw %>%
  filter(Gene_biotype != "protein_coding") %>%
  pull(Genes)

counts_protein_coding <- counts_expressed %>%
  dplyr::select(!any_of(non_protein_coding_genes))

message(
  ncol(counts_expressed) - ncol(counts_protein_coding),
  " non-protein-coding genes removed from the data."
)

## 32221 non-protein-coding genes removed from the data.

Answer: 32,221 genes were removed from the data.

Now that we have finished loading and merging the data, we can remove the unnecessary objects from the global environment with the rm() function. This is good practice with high-dimensional data, as memory can quickly be filled by large, unnecessary objects.

# List the objects to keep
objects_to_keep <- c("metadata_all", "counts_protein_coding")

# Remove everything else from the global environment
rm(list = setdiff(ls(), objects_to_keep))

# Garbage collection to free unused memory
invisible(gc())

2) Normalisation

Normalisation is a crucial step in RNA-seq analysis to ensure that samples are comparable. An important point here is that our data comes from 3 different studies (London, SA, and Leicester). Let’s investigate the impact of this study variable on the data.

2a) Investigate study effects

The library size of a sample is the total number of counts across all genes. Compute the library sizes and visualise them on a bar chart using ggplot2. Colour the samples by the study from which they came. What do you notice?

Save this figure as a PDF in your output folder using ggsave(), with here() to build the path (e.g. here("output", "library_sizes.pdf")). Do the same for every figure you produce in this practical, ensuring that the saved figure has reasonable text sizes, axis names, etc.

# Library size of each sample, with its name and study
library_size_df <- data.frame(
  library_size = rowSums(counts_protein_coding),
  SampleName = metadata_all$SampleName,
  Location = metadata_all$Location
)

library_size_plot <-
  ggplot(library_size_df,
         aes(x = SampleName, y = library_size, fill = Location)) +
  geom_col() +
  theme(axis.text.x = element_blank()) +
  labs(x = "Sample", y = "Library size", title = "Library sizes by sample and location")

# Save the figure as a PDF in the output folder
ggsave(here("output", "library_sizes.pdf"), library_size_plot, width = 12, height = 8)

# Display the figure
library_size_plot

Answer: There are differences in library size both between individuals and between studies. Therefore, these samples are not comparable without normalisation.

2b) Log2-CPM normalisation

In the lecture, we introduced a simple way to normalise data based on library size: log2-counts-per-million (log2-CPM). This is calculated by first dividing each count by its sample’s library size, then multiplying by a million, and finally taking a log2 transformation.

Using the cpm() function from the edgeR package with the log = TRUE argument, compute the log2-CPM of these data.

Tip: If you are unsure how this function works, run ?cpm to bring up its documentation in the Help panel of RStudio.

# cpm() expects genes in rows and samples in columns, so we transpose the counts
# (and transpose the result back, so that samples are rows again)
counts_log2_cpm <- cpm(t(counts_protein_coding), log = TRUE) %>%
  t() %>%
  as.data.frame()

As a sanity check, I recommend that you inspect the data and compare it to the non-normalised data to convince yourself that everything looks OK.

# First 5 samples and genes, before and after normalisation
counts_protein_coding[1:5, 1:5]

##   ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457
## 1               5               0             681             699
## 2               7               0             643             820
## 3               2               0             514             483
## 4              14               0             419             596
## 5              26               1             472             548
##   ENSG00000000460
## 1             141
## 2             116
## 3             118
## 4             190
## 5             141

counts_log2_cpm[1:5, 1:5]

##   ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457
## 1      -1.9161564       -3.592338        4.636814        4.674328
## 2      -1.8660662       -3.592338        4.142860        4.492203
## 3      -2.6531811       -3.592338        4.295027        4.205673
## 4      -0.7631771       -3.592338        3.929357        4.435385
## 5       0.2600683       -2.991041        4.344663        4.559236
##   ENSG00000000460
## 1        2.383151
## 2        1.702604
## 3        2.192350
## 4        2.797838
## 5        2.615326

In the lecture, we introduced principal component analysis (PCA) as a useful tool for representing high-dimensional data in a lower-dimensional space, while preserving the signals which contribute to between-sample variability.

I will first demonstrate how to use the prcomp() function to perform a PCA on the non-normalised data. Since PCA is sensitive to the scale of the data, it is important to center and scale the data (such that every variable has mean 0 and standard deviation 1) using the center and scale. arguments of prcomp().

pca_raw <- prcomp(counts_protein_coding, center = TRUE, scale. = TRUE)

# The coordinates of the samples in the principal component (PC) space.
# Each column is a PC, e.g. the coordinates of the first samples on PC1 are
pc_scores_raw <- pca_raw$x
head(pc_scores_raw[, 1])

## [1] -19.66209 -53.24977 -26.31920 -33.88317 -12.67034 -53.31632

# The proportion of variability explained by each PC
pc_variance_raw <- summary(pca_raw)$importance["Proportion of Variance", ]

# e.g. the proportion explained by PC1 is
pc_variance_raw[1]

##     PC1 
## 0.38552

I will now demonstrate how to visualise the data in a reduced dimension using the results of the PCA. In particular, we might be interested in how the study location (London, Leicester, or South Africa) contributes to variability in the data. To do this, we can plot the samples on the first 2 PCs and colour the points by location.

# Coordinates of the samples on PC1 and PC2, with their location
pca_raw_df <- data.frame(
  PC1 = pc_scores_raw[, 1],
  PC2 = pc_scores_raw[, 2],
  Location = metadata_all$Location
)

# The axis labels show the percentage of variability explained by each PC
pca_raw_plot <- ggplot(pca_raw_df, aes(x = PC1, y = PC2, color = Location)) +
  geom_point(size = 3, alpha = 0.5) +
  labs(
    x = paste0("PC1 (", round(pc_variance_raw[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(pc_variance_raw[2] * 100, 1), "%)"),
    title = "PCA of raw counts"
  )

ggsave(here("output", "pca_raw.pdf"), pca_raw_plot, width = 12, height = 8)

pca_raw_plot

What does this visualisation tell us?

Answer: A large amount of variability in the data can be explained just by one variable. This is likely due to the strong influence of the Location variable, showing that a batch effect is present.

Now, use the example above to perform a PCA on the log2-CPM normalised data and plot the samples on the first two PCs, coloured by location. Has this normalisation step removed the study effect?

pca_log2_cpm <- prcomp(counts_log2_cpm, center = TRUE, scale. = TRUE)
pc_variance_log2_cpm <- summary(pca_log2_cpm)$importance["Proportion of Variance", ]

pca_log2_cpm_df <- data.frame(
  PC1 = pca_log2_cpm$x[, 1],
  PC2 = pca_log2_cpm$x[, 2],
  Location = metadata_all$Location
)

pca_log2_cpm_plot <- ggplot(pca_log2_cpm_df, aes(x = PC1, y = PC2, color = Location)) +
  geom_point(size = 3, alpha = 0.5) +
  labs(
    x = paste0("PC1 (", round(pc_variance_log2_cpm[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(pc_variance_log2_cpm[2] * 100, 1), "%)"),
    title = "PCA of log2-CPM data"
  )

ggsave(here("output", "pca_log2_cpm.pdf"), pca_log2_cpm_plot, width = 12, height = 8)

pca_log2_cpm_plot

Answer: The batch effect is still present. In particular, the samples from the Leicester cohort are very separated from the others on PC1.

2c) Log2-CPM-TMM normalisation

Another typical normalisation approach is to use the trimmed mean of M-values (TMM) method. This calculates normalisation factors which exclude extreme values. The factors are computed on a DGEList object (the edgeR format for count data, with genes in rows) using the calcNormFactors() function, and cpm() then takes them into account. This can be done with the following code.

# Create a DGEList object and compute the TMM normalisation factors
dge_tmm <- DGEList(t(counts_protein_coding)) %>%
  calcNormFactors(method = "TMM")

counts_log2_cpm_tmm <- cpm(dge_tmm, log = TRUE) %>%
  t() %>%
  as.data.frame()

Using your previous code, plot the samples on the first two PCs, coloured by location. Has this normalisation step solved the problem?

pca_log2_cpm_tmm <- prcomp(counts_log2_cpm_tmm, center = TRUE, scale. = TRUE)
pc_variance_log2_cpm_tmm <- summary(pca_log2_cpm_tmm)$importance["Proportion of Variance", ]

pca_log2_cpm_tmm_df <- data.frame(
  PC1 = pca_log2_cpm_tmm$x[, 1],
  PC2 = pca_log2_cpm_tmm$x[, 2],
  Location = metadata_all$Location
)

pca_log2_cpm_tmm_plot <- ggplot(pca_log2_cpm_tmm_df, aes(x = PC1, y = PC2, color = Location)) +
  geom_point(size = 3, alpha = 0.5) +
  labs(
    x = paste0("PC1 (", round(pc_variance_log2_cpm_tmm[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(pc_variance_log2_cpm_tmm[2] * 100, 1), "%)"),
    title = "PCA of log2-CPM-TMM data"
  )

ggsave(here("output", "pca_log2_cpm_tmm.pdf"), pca_log2_cpm_tmm_plot, width = 12, height = 8)

pca_log2_cpm_tmm_plot

Answer: Even after this normalisation, there is still a study effect, especially related to the Leicester location.

2d) Batch effect correction

We can view the study as a batch variable, since we do not expect it to contribute much to the heterogeneity between samples. To correct for batch effects in RNA-seq data, we can use the ComBat-Seq method, with the ComBat_seq() function from the sva package.

In this method, we must specify two things. Firstly, the batch variable (Location), and secondly, a variable whose effects we want to remain intact after correction (in this case, the Status variable, which provides information on the disease states of the patients). We apply ComBat-Seq to the raw, filtered counts to correct for batch = Location, while preserving the effects of group = Status, using the following code (this may take a couple of minutes).

# ComBat_seq() expects genes in rows, so we transpose the counts
counts_combat <- ComBat_seq(
  counts = t(counts_protein_coding),
  batch = metadata_all$Location,
  group = metadata_all$Status
)

## Found 3 batches
## Using full model in ComBat-seq.
## Adjusting for 2 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

Now that we have attempted to correct for the batch effect in the raw counts, use the previous code on this corrected matrix to normalise the data using log2-CPM-TMM. Then, visualise the data on the first two PCs, as before, and conclude whether a study effect remains in the data. Store the normalised counts in an object called counts_combat_log2_cpm_tmm, which is used in the following sections.

# The corrected counts already have genes in rows
dge_combat <- DGEList(counts_combat) %>%
  calcNormFactors(method = "TMM")

counts_combat_log2_cpm_tmm <- cpm(dge_combat, log = TRUE) %>%
  t() %>%
  as.data.frame()

pca_combat <- prcomp(counts_combat_log2_cpm_tmm, center = TRUE, scale. = TRUE)
pc_variance_combat <- summary(pca_combat)$importance["Proportion of Variance", ]

pca_combat_df <- data.frame(
  PC1 = pca_combat$x[, 1],
  PC2 = pca_combat$x[, 2],
  Location = metadata_all$Location
)

pca_combat_plot <- ggplot(pca_combat_df, aes(x = PC1, y = PC2, color = Location)) +
  geom_point(size = 3, alpha = 0.5) +
  labs(
    x = paste0("PC1 (", round(pc_variance_combat[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(pc_variance_combat[2] * 100, 1), "%)"),
    title = "PCA of ComBat-Seq corrected log2-CPM-TMM data"
  )

ggsave(here("output", "pca_combat_location.pdf"), pca_combat_plot, width = 12, height = 8)

pca_combat_plot

Answer: The batch effect appears to have been mitigated.

Since we have finished with the pre-processing, we can remove the unnecessary objects from the global environment again.

objects_to_keep <- c("metadata_all", "counts_combat_log2_cpm_tmm", "counts_combat")

rm(list = setdiff(ls(), objects_to_keep))
invisible(gc())

3) Data Exploration

Now that we have pre-processed the data and (hopefully) removed any batch effects, we can start to explore it. The goal is to use tools like PCA and hierarchical clustering to explore any obvious signals in the data related to the patient status.

3a) Principal Component Analysis

Perform a PCA to explore whether the variability in the gene expression data can be explained by the disease status of patients. Like before, perform a PCA on the batch-corrected, log2-CPM-TMM normalised data and plot the samples on the first two PCs. However, this time, colour the points by their disease status. What does this suggest?

pca_combat <- prcomp(counts_combat_log2_cpm_tmm, center = TRUE, scale. = TRUE)
pc_variance_combat <- summary(pca_combat)$importance["Proportion of Variance", ]

pca_combat_status_df <- data.frame(
  PC1 = pca_combat$x[, 1],
  PC2 = pca_combat$x[, 2],
  Status = metadata_all$Status
)

pca_combat_status_plot <- ggplot(pca_combat_status_df, aes(x = PC1, y = PC2, color = Status)) +
  geom_point(size = 3, alpha = 0.5) +
  labs(
    x = paste0("PC1 (", round(pc_variance_combat[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(pc_variance_combat[2] * 100, 1), "%)"),
    title = "PCA of ComBat-Seq corrected log2-CPM-TMM data, by disease status"
  )

ggsave(here("output", "pca_combat_status.pdf"), pca_combat_status_plot, width = 12, height = 8)

pca_combat_status_plot

Answer: The samples with active tuberculosis seem to systematically vary from the others, suggesting that this disease state is contributing to some of the variability observed in the data.

3b) Hierarchical clustering

In addition, we can explore signal in the data using hierarchical clustering algorithms. Here, the goal is to determine which samples are similar and which are different. If these similarities and/or differences can be explained by the disease state, it suggests a signal might be present.

To perform hierarchical clustering, we need three ingredients:

• the dataset;
• a distance metric, defining the similarity between data points;
• a linkage function, defining the similarity between clusters.

In this first example, I will demonstrate how to perform agglomerative hierarchical clustering on the gene expression data using Euclidean distance and complete linkage, and plot the results as a dendrogram. First, we extract the sample names and disease statuses, which we will need later for the plot.

sample_names <- metadata_all$SampleName
sample_status <- metadata_all$Status

We can now compute the distances between the samples using the dist() function.

sample_distances <- dist(counts_combat_log2_cpm_tmm, method = "euclidean")

Then, we can use the hclust() function from the flashClust package to perform agglomerative hierarchical clustering based on these distances.

hclust_complete <- hclust(sample_distances, method = "complete")

A dendrogram can be created from this result with the as.dendrogram() function. I will also provide some code to plot it with the sample labels coloured by the Status variable.

Do you notice any patterns? Is there any obvious way to determine the number of clusters in this case?

dendrogram_complete <- as.dendrogram(hclust_complete)

# Order of the samples along the dendrogram, used to label the leaves in the right order
leaf_order <- order.dendrogram(dendrogram_complete)
labels(dendrogram_complete) <- sample_names[leaf_order]

# Convert the dendrogram into a format suitable for ggplot2
dendrogram_complete_data <- dendro_data(dendrogram_complete)

leaf_labels_complete <- data.frame(
  label = sample_names[leaf_order],
  Status = sample_status[leaf_order],
  x = dendrogram_complete_data$labels$x
)

dendrogram_complete_plot <- ggplot() +
  geom_segment(
    data = dendrogram_complete_data$segments,
    aes(x = x, y = y, xend = xend, yend = yend)
  ) +
  geom_text(
    data = leaf_labels_complete,
    aes(x = x, y = 0, label = label, color = Status),
    hjust = 1,
    vjust = 0.5,
    angle = 90,
    size = 3
  ) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Hierarchical clustering of samples (Euclidean distance, complete linkage)",
    x = "Samples",
    y = "Height"
  ) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

ggsave(here("output", "dendrogram_complete.pdf"), dendrogram_complete_plot, width = 15, height = 9)

dendrogram_complete_plot

Answer: Some of the red samples (active TB) tend to cluster together, although these clusters are not very distinct. This makes sense given what we previously observed with the PCA, which revealed that the status contributes to the variability in the data, but not hugely so.

Now change the linkage function by setting method = "ward" in the hclust() call. Use the previous code to plot the resulting dendrogram. A criterion to select the number of clusters which often accompanies Ward’s method is to choose the number of clusters such that the next jump in the tree is the greatest. How many clusters would you choose with this criterion? What is your conclusion now?

hclust_ward <- hclust(sample_distances, method = "ward")

dendrogram_ward <- as.dendrogram(hclust_ward)
leaf_order <- order.dendrogram(dendrogram_ward)
labels(dendrogram_ward) <- sample_names[leaf_order]
dendrogram_ward_data <- dendro_data(dendrogram_ward)

leaf_labels_ward <- data.frame(
  label = sample_names[leaf_order],
  Status = sample_status[leaf_order],
  x = dendrogram_ward_data$labels$x
)

dendrogram_ward_plot <- ggplot() +
  geom_segment(
    data = dendrogram_ward_data$segments,
    aes(x = x, y = y, xend = xend, yend = yend)
  ) +
  geom_text(
    data = leaf_labels_ward,
    aes(x = x, y = 0, label = label, color = Status),
    hjust = 1,
    vjust = 0.5,
    angle = 90,
    size = 3
  ) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Hierarchical clustering of samples (Euclidean distance, Ward linkage)",
    x = "Samples",
    y = "Height"
  ) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

ggsave(here("output", "dendrogram_ward.pdf"), dendrogram_ward_plot, width = 15, height = 9)

dendrogram_ward_plot

Answer: Now the samples appear to cluster together more strongly by status. A large jump in the tree can be observed between around 600 and 800 distance. Therefore, using this method, we would cut the tree at the last jump, resulting in two clusters.

Now clear the environment of unnecessary objects again.

objects_to_keep <- c("metadata_all", "counts_combat_log2_cpm_tmm", "counts_combat")

rm(list = setdiff(ls(), objects_to_keep))
invisible(gc())

4) Differential Gene Expression Analysis

We are interested in identifying genes whose expression differs between disease states. For the purposes of this exercise, we will group the control and latent tuberculosis (LTBI) patients together, and compare them to the active infection group. We will use multiple methods to do this, and compare their results at the end. In each case, we define the differentially expressed genes (DEGs) as those with adjusted p-value < 0.05 and |log2-fold-change| > 1.

4a) EdgeR

We will be following the edgeR pipeline, found here.

edgeR starts with raw integer counts. We can use the counts from ComBat-Seq (prior to normalisation). Starting from a DGEList object, follow the edgeR pipeline. This starts with estimating TMM normalisation factors with calcNormFactors(). Then, we estimate a common dispersion using estimateCommonDisp(). The negative binomial model can be fit with glmQLFit(). P-values can be extracted with glmQLFTest(), and the top genes can be extracted using the topTags() function.

First, I will give you the code to define the groups we want to compare. The variable testGroup groups the control and LTBI patients into one control group, and the design matrix design describes this comparison (its coefficient testGroupActiveTB compares the ActiveTB group to the control group). We will reuse design in the next sections.

# New variable grouping the control and LTBI patients into one control group
metadata_all <- metadata_all %>%
  mutate(testGroup = recode_factor(Status, "LTBI" = "Control")) %>%
  relocate(testGroup, .after = Status)

# Design matrix : the coefficient "testGroupActiveTB" compares ActiveTB to control
design <- model.matrix(~ testGroup, data = metadata_all)

Extract a vector of the names of the DEGs. How many genes are differentially expressed according to this criterion?

# DGEList object with the ComBat-Seq corrected counts, and TMM normalisation factors
dge_edgeR <- DGEList(counts = counts_combat, group = metadata_all$testGroup) %>%
  calcNormFactors(method = "TMM")

# Estimate a common dispersion for all genes
dge_common_disp <- estimateCommonDisp(dge_edgeR, design)

# Fit the quasi-likelihood GLM and test the coefficient of interest (coef = 2, ActiveTB vs control)
fit_common_disp <- glmQLFit(dge_common_disp, design)
qlf_common_disp <- glmQLFTest(fit_common_disp, coef = 2)

# Results for all genes, and names of the DEGs
results_edgeR_common <- topTags(qlf_common_disp, n = Inf) %>%
  as.data.frame()

deg_edgeR_common <- results_edgeR_common %>%
  filter(FDR < 0.05, abs(logFC) > 1) %>%
  rownames()

message(length(deg_edgeR_common), " differentially expressed genes found.")

## 2561 differentially expressed genes found.

Answer: 2,561 differentially expressed genes were found.

Repeat the previous experiment, but instead of estimating a common dispersion, estimate one for each gene using the estimateDisp() function. Do the results change significantly?

# Estimate a dispersion for each gene
dge_genewise_disp <- estimateDisp(dge_edgeR, design)

fit_genewise_disp <- glmQLFit(dge_genewise_disp, design)
qlf_genewise_disp <- glmQLFTest(fit_genewise_disp, coef = 2)

results_edgeR_genewise <- topTags(qlf_genewise_disp, n = Inf) %>%
  as.data.frame()

deg_edgeR_genewise <- results_edgeR_genewise %>%
  filter(FDR < 0.05, abs(logFC) > 1) %>%
  rownames()

message(length(deg_edgeR_genewise), " differentially expressed genes found.")

## 2537 differentially expressed genes found.

# DEGs found with both common and gene-wise dispersions
message(
  length(intersect(deg_edgeR_common, deg_edgeR_genewise)),
  " differentially expressed genes found with both methods."
)

## 2534 differentially expressed genes found with both methods.

Answer: 2,534 differentially expressed genes were found, which is not too different.

4b) Voom-Limma

Now, we will perform the same task using voom-limma. A tutorial can be found here.

Following the guide, apply the voom-limma pipeline with the functions voom(), lmFit(), eBayes() and decideTests(), reusing the design matrix defined above. The output should be a vector giving the names of the DEGs. How many genes are differentially expressed?

# DGEList object with the ComBat-Seq corrected counts, and TMM normalisation factors
dge_voom <- DGEList(counts = counts_combat, group = metadata_all$testGroup) %>%
  calcNormFactors(method = "TMM")

# Estimate the mean-variance relationship with voom(), and save the diagnostic plot
pdf(here("output", "voom_mean_variance.pdf"), width = 8, height = 6)
voom_data <- voom(dge_voom, design, plot = TRUE)
invisible(dev.off())

# Fit the linear models, then compute the empirical Bayes moderated tests
fit_voom <- lmFit(voom_data, design)
fit_voom <- eBayes(fit_voom, robust = TRUE)

# Number of down- and up-regulated genes (adjusted p-value < 0.05 and |log2FC| > 1)
voom_decisions <- decideTests(fit_voom, adjust.method = "BH", p.value = 0.05, lfc = 1)
summary(voom_decisions)

##        (Intercept) testGroupActiveTB
## Down          6902               357
## NotSig        1435             18649
## Up           11050               381

# Results for all genes, and names of the DEGs
results_voom <- topTable(fit_voom, coef = 2, n = Inf)

deg_voom <- results_voom %>%
  filter(adj.P.Val < 0.05, abs(logFC) > 1) %>%
  rownames()

message(length(deg_voom), " differentially expressed genes found.")

## 738 differentially expressed genes found.

Answer: 738 differentially expressed genes were found.

4c) dearseq

Finally we will apply the dearseq pipeline. A guide can be found here.

Use the dear_seq() function to perform differential expression analysis. Since log2-fold-change values are not directly available from the dearseq package, you may use the ones calculated by edgeR.

dearseq_test <- dear_seq(
  # Log2-CPM-TMM data, with genes in rows
  exprmat = t(counts_combat_log2_cpm_tmm),
  # The variable to test (ActiveTB vs control)
  variables2test = design[, "testGroupActiveTB", drop = FALSE],
  # The variables to adjust for (here, only the intercept)
  covariates = design[, "(Intercept)", drop = FALSE],
  which_test = "asymptotic",
  # The data are already normalised and log2-transformed
  preprocessed = TRUE
)

# Adjusted p-values from dearseq, with the log2FC from edgeR (the genes are in the same order)
results_dearseq <- data.frame(
  adj.P.Val = dearseq_test$pvals$adjPval,
  logFC = qlf_common_disp$table$logFC,
  row.names = rownames(counts_combat)
)

deg_dearseq <- results_dearseq %>%
  filter(adj.P.Val < 0.05, abs(logFC) > 1) %>%
  rownames()

message(length(deg_dearseq), " differentially expressed genes found.")

## 618 differentially expressed genes found.

Answer: 618 differentially expressed genes were found.

4d) Volcano plot

A common representation of the results of differential gene expression analysis is the so-called volcano plot. The x-axis of a volcano plot is the log2FC value and the y-axis is the adjusted p-value on the -log10 scale, such that smaller p-values are higher on the y-axis. Typically, the thresholds which define differential expression (e.g. adjusted p-value < 0.05 and |log2FC| > 1 in our case) are represented by lines on the plot. Significantly up-regulated genes are coloured red and significantly down-regulated genes blue. Sometimes, the few most significant genes are labelled with their name.

Choose one of your analysis pipelines, and make a volcano plot to represent the results. Annotate the top 10 most significantly differentially expressed genes (by adjusted p-value) on the plot. Save the figure.

# Volcano plot of a results table with columns logFC and adj.P.Val (genes in rownames)
plot_volcano <- function(results, title) {
  results <- results %>%
    mutate(
      gene = rownames(results),
      expression = case_when(
        adj.P.Val < 0.05 & logFC > 1 ~ "Up-regulated",
        adj.P.Val < 0.05 & logFC < -1 ~ "Down-regulated",
        TRUE ~ "Not significant"
      )
    )

  # The 10 most significant DEGs, to be labelled
  top_genes <- results %>%
    filter(expression != "Not significant") %>%
    arrange(adj.P.Val) %>%
    head(10)

  ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = expression)) +
    geom_point(alpha = 0.7, size = 1.5) +
    scale_color_manual(values = c(
      "Up-regulated" = "red",
      "Down-regulated" = "blue",
      "Not significant" = "grey"
    )) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    geom_text_repel(
      data = top_genes,
      aes(label = gene),
      size = 3,
      box.padding = 0.5,
      point.padding = 0.5,
      force = 10,
      segment.color = "grey50"
    ) +
    labs(
      title = title,
      x = "log2 fold change",
      y = "-log10 adjusted p-value",
      color = "Expression"
    ) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"))
}

The function above builds the volcano plot of any results table, so we can use it for the three pipelines.

# edgeR (common dispersion). edgeR calls the adjusted p-value "FDR"
volcano_edgeR_plot <- results_edgeR_common %>%
  rename(adj.P.Val = FDR) %>%
  plot_volcano(title = "Volcano plot : edgeR")
ggsave(here("output", "volcano_edgeR.pdf"), volcano_edgeR_plot, width = 10, height = 8)
volcano_edgeR_plot

# voom-limma
volcano_voom_plot <- plot_volcano(results_voom, title = "Volcano plot : voom-limma")
ggsave(here("output", "volcano_voom.pdf"), volcano_voom_plot, width = 10, height = 8)
volcano_voom_plot

# dearseq
volcano_dearseq_plot <- plot_volcano(results_dearseq, title = "Volcano plot : dearseq")
ggsave(here("output", "volcano_dearseq.pdf"), volcano_dearseq_plot, width = 10, height = 8)
volcano_dearseq_plot

4e) Results comparison

Finally, let’s make some visualisations to compare the results of the 3 pipelines (voom-limma, edgeR and dearseq).

First, use the VennDiagram package to visualise the agreement between the methods.

# The DEGs of each pipeline
deg_lists <- list(
  VoomLimma = deg_voom,
  EdgeR = deg_edgeR_common,
  dearseq = deg_dearseq
)

# Stop VennDiagram from writing a log file in the working directory
futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")

## NULL

venn_plot <- venn.diagram(
  x = deg_lists,
  filename = NULL,
  fill = c("red", "blue", "green"),  # Colour of each set
  alpha = 0.5,                       # Transparency of the fills
  lwd = 2,                           # Line width of the circles
  col = "black",                     # Line colour of the circles
  cex = 2,                           # Text size
  fontface = "bold",
  fontfamily = "sans"
)

# Save the Venn diagram as a PDF in the output folder
pdf(here("output", "venn_diagram.pdf"), width = 8, height = 8)
grid.newpage()
grid.draw(venn_plot)
invisible(dev.off())

# Display the Venn diagram
grid.newpage()
grid.draw(venn_plot)

How many genes are identified by all three methods? Do the methods broadly agree?

deg_all_methods <- Reduce(intersect, deg_lists)

message(length(deg_all_methods), " genes were identified by all three methods.")

## 455 genes were identified by all three methods.

Answer: 455 genes were identified by all three methods. In particular, edgeR identifies a large number of distinct genes.

An alternative way to visualise agreement is with a so-called UpSet plot. This representation works well when the number of sets to compare is greater than 3, such that a Venn diagram becomes hard to read. Using the UpSetR package, make an UpSet plot. (Hint : the fromList() function converts a list of gene names into the input format of upset()). Do our 3 pipelines have concordant results? Make a conclusion about the sensitivity of the results to the choice of method.

upset_plot <- upset(
  fromList(deg_lists),
  sets = names(deg_lists),
  order.by = "freq",              # Order the intersections by size
  empty.intersections = "on",     # Include empty intersections
  main.bar.color = "black",
  sets.bar.color = c("red", "blue", "green")
)

# Save the UpSet plot as a PDF in the output folder
pdf(here("output", "upset_plot.pdf"), width = 10, height = 6)
print(upset_plot)
invisible(dev.off())

# Display the UpSet plot
upset_plot

Answer: Although some of the results contain many shared genes, there are also a significant number which are unique to a given method, suggesting that the results are sensitive to the choice of method. It may be beneficial to apply multiple methods to the same data and compare the results.