```{r, include = FALSE} knitr::opts_chunk$set(error = TRUE) ``` # 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](https://www.nature.com/articles/s41467-018-04579-w). 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. ```{r, eval = F} 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. ```{r, echo = F, results = 'hide', message = F} # This sets the options for the outputs when knitting the rmd file knitr::opts_chunk$set(warning = FALSE, message = TRUE, fig.width = 12, fig.height = 8) ``` ```{r, eval = F, results = 'hide', message = FALSE, warning = FALSE} 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"])) ``` ```{r, echo = F, results = 'hide', message = F} library(BiocManager) library(GEOquery) library(edgeR) library(limma) library(pvca) library(tidyr) library(readxl) library(dplyr) library(flashClust) library(UpSetR) library(ggplot2) library(stringr) library(dendextend) library(ggdendro) library(VennDiagram) library(grid) library(ggrepel) library(dearseq) library(sva) library(here) # Default look of all the ggplot2 figures theme_set(theme_minimal(base_size = 14)) ``` ## 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. ```{r} # Path to the root of the project here() # 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))) ``` > 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](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE107995) 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. ```{r, message = F} 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. ```{r, results = 'hide', message = F} 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. ```{r} # 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. ```{r} # 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) ``` # 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? ```{r} # Your code here ``` ## 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: ```{r} london_raw$Gene_biotype %>% head() ``` > 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. ```{r} # Your code here ``` > 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. ```{r} # 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. ```{r} # Your code here ``` ## 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. ```{r} # Your code here ``` > 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. ```{r} # Your code here ``` > 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()`. ```{r} 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]) # 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] ``` > 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. ```{r} # 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? > 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? ```{r} # Your code here ``` ## 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. ```{r} # 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? ```{r} # Your code here ``` ## 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). ```{r} # 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 ) ``` > 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. ```{r} # Your code here ``` > Since we have finished with the pre-processing, we can remove the unnecessary objects from the global environment again. ```{r} 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? ```{r} # Your code here ``` ## 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. ```{r} sample_names <- metadata_all$SampleName sample_status <- metadata_all$Status ``` > We can now compute the distances between the samples using the `dist()` function. ```{r} 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. ```{r} 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? ```{r} 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 ``` > 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? ```{r} # Your code here ``` > Now clear the environment of unnecessary objects again. ```{r} 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](https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf). > 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. ```{r} # 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? ```{r} # Your code here ``` > 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? ```{r} # Your code here ``` ## 4b) Voom-Limma > Now, we will perform the same task using voom-limma. A tutorial can be [found here](https://ucdavis-bioinformatics-training.github.io/2018-June-RNA-Seq-Workshop/thursday/DE.html). > 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? ```{r} # Your code here ``` ## 4c) dearseq > Finally we will apply the dearseq pipeline. A guide can be [found here](https://www.bioconductor.org/packages/release/bioc/vignettes/dearseq/inst/doc/dearseqUserguide.html). > 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`. ```{r} # Your code here ``` ## 4d) Volcano plot > A common representation of the results of differential gene expression analysis is the so-called [*volcano plot*](https://olvtools.com/en/rnaseq/help/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. ```{r} # Your code here ``` ## 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. ```{r} # Your code here ``` > How many genes are identified by all three methods? Do the methods broadly agree? ```{r} # Your code here ``` > 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. ```{r} # Your code here ```