Workshop

Transcripttranscriptomics 2: Statistical Analysis

Author

Emma Rand

Published

24 August, 2025

Introduction

Session overview

In the workshop, you will learn how to perform differential expression analysis on raw counts using DESeq2 (Love, Huber, and Anders 2014) or on logged normalised expression values using scran (Lun, McCarthy, and Marioni 2016) or both.

Set up

🎄 Arabidopisis

🎬 Open the arab-88H RStudio Project and the cont-low-root.R script.

💉 Leishmania

🎬 Open the leish-88H RStudio Project and the pro-meta.R script.

🐭 Stem cells

🎬 Open the mice-88H RStudio Project and the hspc-prog.R script.

Everyone

🎬 Make a new folder results in the project directory.

This is where we will save our results.

🎬 Load tidyverse (Wickham et al. 2019) You most likely have this code at the top of your script already.

library(tidyverse)

── Attaching core tidyverse packages ─────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package to force all conflicts to become errors

Have you ever stopped to think about this message? It is telling us that there are functions in the dplyr package that have the same name as functions in the stats package and that R will use the dplyr version. As this is what you want, this has always been fine. It still is fine in this case. However, as you start to load more packages, you will want to know if you are using a function from a package that has the same name as a function in another loaded package. This is where the conflicted (Wickham 2023) package comes in. Conflicted will warn you when you are using a function that has the same name as a function in another package. You can then choose which function to use.

🎬 Load the conflicted package:

library(conflicted)

Instead of getting a warning every time you are using a function that has a function with the same name in another package, we can declare a preference for one function over another. This is useful for the functions you use a lot or ones where you are certain you always want to use a particular function.

For example, to always use the dplyr version of filter() by default you can add this to the top of your script:

conflicts_prefer(dplyr::filter)

We will also want to ensure that we are using the setdiff() function from the GenomicRanges package which is used by DESeq2.

conflicts_prefer(GenomicRanges::setdiff)

Import

🎄 Arabidopisis

We need to import the root tissue data that were filtered to remove genes with 3 or 4 zeros and those where the total counts was less than 20.

🎬 Import the data from the data-processed folder.

Now go to Differential Expression Analysis.

💉 Leishmania

We need to import the procyclic- and metacyclic-promastigote data that were filtered to remove genes with 4, 5 or 6 zeros and those where the total counts was less than 20.

🎬 Import the data from the data-processed folder.

Now go to Differential Expression Analysis.

🐭 Stem cells

We need to import the progenitor and HSPC cell data that were combined.

🎬 Import the data from the data-processed folder..

Now go to Differential Expression Analysis.

Differential Expression Analysis

🎄 Arabidopisis

These are the steps we will take

Find the genes that are expressed in only one treatment group.
Create a DESeqDataSet object. This is a special object that is used by the DESeq2 package
Prepare the normalised counts from the DESeqDataSet object.
Do differential expression analysis on the genes. This needs to be done on the raw counts.

All but the first step are done with the DESeq2 package

1. Genes expressed in one treatment

The genes expressed in only one treatment group are those with zeros in all replicates in one group and non-zero values in all replicates in the other group. For example, those shown here:

gene_id	gene_name	LWR1	LWR2	LWR3	LWR4	LWR5	LWR6
AT1G63820	AT1G63820	13	4	4	7	10	9
AT1G27090	AT1G27090	16	6	6	9	13	14
AT1G14950	AT1G14950	14	7	4	9	13	17
AT1G18350	MKK7	6	8	4	4	6	8
AT1G65910	NAC028	15	11	5	1	3	15

We will use filter() to find these genes.

🎬 Find the genes that are expressed only in the control group:

root_lowni_only <- root_filtered |>
  filter(if_all(starts_with("CTR"), \(x) x == 0),
         if_all(starts_with("LWR"), \(x) x > 0) )

This code creates a new data frame called root_lowni_only that contains only the rows from root_filtered where:

All CTR columns have a value of 0
All LWR columns have a value greater than 0

filter() keeps rows that match the conditions you give it. Here we give two conditions:

All columns whose names start with “CTR” must have a value of 0.
All columns whose names start with “LWR” must have a value greater than 0.

The if_all() function applies a test to each selected column (those starting with “CTR” or “LWR”) for each row. It returns TRUE only if all the selected columns pass the test.

The bit \(x) x == 0 is an anonymous function — a function without a name. In base R syntax, it’s equivalent to: function(x) x == 0 Here, x will be the values from one column at a time, and the function checks whether they are equal to 0. We use an anonymous function because the condition needs to be applied separately to each column, but it’s simple enough that there’s no need to define a full named function elsewhere in the script.

❓ How many genes are expressed only in the low nickel group?

Note: There may have been more genes with counts in the low nickel group that were removed when we filtered on the total number of counts being less that 50.

🎬 Now you find any genes that are expressed only in the control group.

❓ How many genes are expressed only in the control group?

❓ Do the results make sense to you in light of what you know about the biology?

🎬 Write all the genes that are expressed one group only to file (saved in results)

2. Create DESeqDataSet object

🎬 Load the DESeq2 package:

A DEseqDataSet object is a custom data type that is used by DESeq2. Custom data types are common in the Bioconductor¹ packages. They are used to store data in a way that is useful for the analysis. These data types typically have data, transformed data, metadata and experimental designs within them.

To create a DESeqDataSet object, we need to provide three things:

The raw counts - these are in root_filtered
The meta data which gives information about the samples and which treatment groups they belong to
A design matrix which captures the design of the statistical model.

The counts must in a matrix rather than a dataframe. Unlike a dataframe, a matrix has columns of all the same type. That is, it will contain only the counts. The gene ids are given as row names rather than a column. The matrix() function will create a matrix from a dataframe of columns of the same type and the select() function can be used to remove the gene ids column.

🎬 Create a matrix of the counts:

root_count_mat <- root_filtered |>
  select(-gene_id, -gene_name) |>
  as.matrix()

🎬 Add the gene ids as row names to the matrix:

# add the row names to the matrix
rownames(root_count_mat) <- root_filtered$gene_id

You might want to view the matrix (click on it in your environment pane).

The metadata are in a file, arab_meta_data.txt. This is a tab-delimited file. The first column is the sample name and the other columns give the “treatments”. In this case, the treatments are tissue (with two levels) and nickel (with two levels).

🎬 Make a folder called meta and save the file to it.

🎬 Read the metadata into a dataframe:

meta <- read_table("meta/arab_meta_data.txt")

🎬 Examine the resulting dataframe.

We need to add the sample names as row names to the metadata dataframe. This is because the DESeqDataSet object will use the row names to match the samples in the metadata to the samples in the counts matrix.

🎬 Add the sample names as row names to the metadata dataframe:

meta <- meta |>
  column_to_rownames("sample_id")

We are dealing only with the root data so we need to remove the samples that are not in the root data.

🎬 Filter the metadata to keep only the root information:

meta_root <- meta |>
  filter(tissue == "root")

We can now create the DESeqDataSet object. The design formula describes the statistical model. You should recognise the form from previous work. The ~ can be read as “explain by” and on its right hand side are the explanatory variables. That is, the model is counts explained by nickel status.

Note that:

The names of the columns in the count matrix have to exactly match the names of the rows in the metadata dataframe. They also need to be in the same order.
The names of the explanatory variables in the design formula have to match the names of columns in the metadata.

🎬 Create the DESeqDataSet object:

dds <- DESeqDataSetFromMatrix(root_count_mat,
                              colData = meta_root,
                              design = ~ nickel)

The warning “Warning: some variables in design formula are characters, converting to factors” just means that the variable type of nickel in the metadata dataframe is “char” and it has been converted into a factor type.

To help you understand what the DESeqDataSet object we have called dds contains, we can look its contents

The counts are in dds@assays@data@listData[["counts"]] and the metadata are in dds@colData but the easiest way to see them is to use the counts() and colData() functions from the DESeq2 package.

🎬 View the counts:

counts(dds) |> View()

You should be able to see that this is the same as in root_count_mat.

🎬 View the column information:

colData(dds)

DataFrame with 12 rows and 2 columns
          tissue   nickel
     <character> <factor>
CTR1        root  control
CTR2        root  control
CTR3        root  control
CTR4        root  control
CTR5        root  control
...          ...      ...
LWR2        root   low_ni
LWR3        root   low_ni
LWR4        root   low_ni
LWR5        root   low_ni
LWR6        root   low_ni

You should be able to see this is the same as in meta_root.

3. Prepare the normalised counts

The normalised counts are the counts that have been transformed to account for the library size (i.e., the total number of reads in a sample) and the gene length. We have to first estimate the normalisation factors and store them in the DESeqDataSet object and then we can get the normalised counts.

🎬 Estimate the factors for normalisation and store them in the DESeqDataSet object:

dds <- estimateSizeFactors(dds)

🎬 Look at the factors (just for information):

sizeFactors(dds)

     CTR1      CTR2      CTR3      CTR4      CTR5      CTR6      LWR1      LWR2 
1.0241774 1.3312269 1.0264881 1.1672503 1.0274692 1.0104721 0.9390444 0.8343701 
     LWR3      LWR4      LWR5      LWR6 
0.8624386 0.9117912 1.0260560 0.9740408

The normalised counts will be useful to use later. To get the normalised counts we again used the counts() function but this time we use the normalized=TRUE argument.

🎬 Save the normalised to a matrix:

normalised_counts <- counts(dds, normalized = TRUE)

🎬 Make a dataframe of the normalised counts, adding a column for the gene ids at the same time:

root_normalised_counts <- data.frame(normalised_counts,
                                    gene_id = row.names(normalised_counts))

4. Differential expression analysis

We use the DESeq() function to do the differential expression analysis. This function fits the statistical model to the data and then uses the model to calculate the significance of the difference between the treatments. It again stores the results in the DESseqDataSet object. Note that the differential expression needs the raw (unnormalised counts) as it does its own normalisation as part of the process.

🎬 Run the differential expression analysis and store the results in the same object:

dds <- DESeq(dds)

The function will take only a few moments to run on this data but can take longer for bigger datasets.

We need to define the contrasts we want to test. We want to test the difference between the treatments so we will define the contrast as control and low_ni.

🎬 Define the contrast:

contrast_suf <- c("nickel", "control", "low_ni")

Note that nickel is the name of the column in the metadata dataframe and control and low_ni are the names of the levels in the nickel column. By putting them in the order control , low_ni we are saying the fold change will be control / low_ni. This means:

positive log fold changes indicate control > low_ni and
negative log fold changes indicates low_ni > control.

If we had put them in the order low_ni, control we would have the reverse.

🎬 Extract the results from the DESseqDataSet object:

results_suf <- results(dds,
                       contrast = contrast_suf)

This will give us the log₂ fold change, the p-value and the adjusted p-value for the comparison between the control and low_ni for each gene.

🎬 Put the results in a dataframe and add the gene ids as a column:

root_results <- data.frame(results_suf,
                          gene_id = row.names(results_suf))

It is useful to have the normalised counts and the statistical results in one dataframe.

🎬 Merge the two dataframes:

# merge the results with the normalised counts
root_results <- root_normalised_counts |>
  left_join(root_results, by = "gene_id")

Now go to Add gene information.

💉 Leishmania

These are the steps we will take

Find the genes that are expressed in only one treatment group.
Create a DESeqDataSet object. This is a special object that is used by the DESeq2 package
Prepare the normalised counts from the DESeqDataSet object.
Do differential expression analysis on the genes. This needs to be done on the raw counts.

All but the first step are done with the DESeq2 package

1. Genes expressed in one treatment

The genes expressed in only one treatment group are those with zeros in all replicates in one group and non-zero values in all replicates in the other group.

We will use filter() to find these genes.

🎬 Find the genes that are expressed only at the procyclic-promastigote stage:

pro_meta_pro_only <- pro_meta_filtered  |>
  filter(lm_pro_1 > 0,
         lm_pro_2 > 0,
         lm_pro_3 > 0,
         lm_meta_1 == 0,
         lm_meta_2 == 0,
         lm_meta_2 == 0)

❓ How many genes are expressed only in the procyclic-promastigote stage group?

🎬 Now you find any genes that are expressed only at the metacyclic stage

❓ How many genes are expressed only at the metacyclic stage?

❓ Do the results make sense to you in light of what you know about the biology?

🎬 Write all the genes that are expressed one group only to file (saved in results)

2. Create DESeqDataSet object

🎬 Load the DESeq2 package:

A DEseqDataSet object is a custom data type that is used by DESeq2. Custom data types are common in the Bioconductor² packages. They are used to store data in a way that is useful for the analysis. These data types typically have data, transformed data, metadata and experimental designs within them.

To create a DESeqDataSet object, we need to provide three things:

The raw counts - these are in pro_meta_filtered
The meta data which gives information about the samples and which treatment groups they belong to
A design matrix which captures the design of the statistical model.

🎬 Create a matrix of the counts:

pro_meta_count_mat <- pro_meta_filtered  |>
  select(-gene_id) |>
  as.matrix()

🎬 Add the gene ids as row names to the matrix:

# add the row names to the matrix
rownames(pro_meta_count_mat) <- pro_meta_filtered$gene_id

You might want to view the matrix (click on it in your environment pane).

The metadata are in a file, leish_meta_data.txt. This is a tab-delimited file. The first column is the sample name and the other columns give the “treatments”. In this case, the treatment is stage (with three levels).

🎬 Make a folder called meta and save the file to it.

🎬 Read the metadata into a dataframe:

meta <- read_table("meta/leish_meta_data.txt")

🎬 Examine the resulting dataframe.

🎬 Add the sample names as row names to the metadata dataframe:

meta <- meta |>
  column_to_rownames("sample_id")

We are dealing only with the wild data so we need to remove the samples that are not in the wild data.

🎬 Filter the metadata to keep only the procyclic and metacyclic information:

meta_pro_meta <- meta |>
  filter(stage != "amastigotes")

Note that:

The names of the columns in the count matrix have to exactly match the names of the rows in the metadata dataframe. They also need to be in the same order.
The names of the explanatory variables in the design formula have to match the names of columns in the metadata.

🎬 Create the DESeqDataSet object:

dds <- DESeqDataSetFromMatrix(pro_meta_count_mat,
                              colData = meta_pro_meta,
                              design = ~ stage)

The warning “Warning: some variables in design formula are characters, converting to factors” just means that the variable type of stage in the metadata dataframe is “char” and it has been converted into a factor type.

To help you understand what the DESeqDataSet object we have called dds contains, we can look its contents

🎬 View the counts:

counts(dds) |> View()

You should be able to see that this is the same as in pro_meta_count_mat.

🎬 View the column information:

colData(dds)

DataFrame with 6 rows and 2 columns
               stage replicate
            <factor> <numeric>
lm_pro_1  procyclic          1
lm_pro_2  procyclic          2
lm_pro_3  procyclic          3
lm_meta_1 metacyclic         1
lm_meta_2 metacyclic         2
lm_meta_3 metacyclic         3

You should be able to see this is the same as in meta_pro_meta.

3. Prepare the normalised counts

🎬 Estimate the factors for normalisation and store them in the DESeqDataSet object:

dds <- estimateSizeFactors(dds)

🎬 Look at the factors (just for information):

sizeFactors(dds)

 lm_pro_1  lm_pro_2  lm_pro_3 lm_meta_1 lm_meta_2 lm_meta_3 
1.3029351 0.9158157 0.9943186 0.7849299 0.8443586 1.3250409

The normalised counts will be useful to use later. To get the normalised counts we again used the counts() function but this time we use the normalized=TRUE argument.

🎬 Save the normalised to a matrix:

normalised_counts <- counts(dds, normalized = TRUE)

🎬 Make a dataframe of the normalised counts, adding a column for the gene ids at the same time:

pro_meta_normalised_counts <- data.frame(normalised_counts,
                                    gene_id = row.names(normalised_counts))

4. Differential expression analysis

🎬 Run the differential expression analysis and store the results in the same object:

dds <- DESeq(dds)

The function will take only a few moments to run on this data but can take longer for bigger datasets.

We need to define the contrasts we want to test. We want to test the difference between the treatments so we will define the contrast as procyclic and metacyclic.

🎬 Define the contrast:

contrast_pro_meta <- c("stage", "procyclic", "metacyclic")

Note that stage is the name of the column in the metadata dataframe and procyclic and metacyclic are the names of the levels in the stage column. By putting them in the order procyclic , metacyclic we are saying the fold change will be procyclic / metacyclic. This means:

positive log fold changes indicate procyclic > metacyclic and
negative log fold changes indicates metacyclic > procyclic.

If we had put them in the order metacyclic, procyclic we would have the reverse.

🎬 Extract the results from the DESseqDataSet object:

results_pro_meta <- results(dds,
                       contrast = contrast_pro_meta)

This will give us the log₂ fold change, the p-value and the adjusted p-value for the comparison between procyclic and metacyclic stage for each gene

🎬 Put the results in a dataframe and add the gene ids as a column:

pro_meta_results <- data.frame(results_pro_meta,
                          gene_id = row.names(results_pro_meta))

It is useful to have the normalised counts and the statistical results in one dataframe.

🎬 Merge the two dataframes:

# merge the results with the normalised counts
pro_meta_results <- pro_meta_normalised_counts |>
  left_join(pro_meta_results, by = "gene_id")

Now go to Add gene information.

🐭 Stem cells

These are the steps we will take

Find the genes that are expressed in only one cell type (the prog or the hspc).
Prepare the data for differential expression analysis with the scran package.
Do differential expression analysis on the genes using the scran package. This needs to be done on the logged normalised counts.

1. Genes expressed in one cell type

The genes expressed in only cell type are those with zeros in all the cells of the other type. We can find these by summing the expression values for each gene across the cells of one type and filtering for those that are zero.

To do row wise aggregates such as the sum across rows we can use the rowwise() function. c_across() allows us to use the colon notation to select columns. This is very useful when you have a lot of columns because it would be annoying to have to list all of them. HSPC_001:HSPC_852 means all the columns from HSPC_001 to HSPC_852.

🎬 Find the genes expressed only in progenitor cells (i.e., those that are 0 in every HSPC cell):

hspc_prog |> 
  rowwise() |> 
  filter(sum(c_across(HSPC_001:HSPC_852)) == 0)

# A tibble: 0 × 1,500
# Rowwise: 
# ℹ 1,500 variables: ensembl_gene_id <chr>, HSPC_001 <dbl>, HSPC_002 <dbl>,
#   HSPC_003 <dbl>, HSPC_004 <dbl>, HSPC_006 <dbl>, HSPC_008 <dbl>,
#   HSPC_009 <dbl>, HSPC_011 <dbl>, HSPC_012 <dbl>, HSPC_014 <dbl>,
#   HSPC_015 <dbl>, HSPC_016 <dbl>, HSPC_017 <dbl>, HSPC_018 <dbl>,
#   HSPC_020 <dbl>, HSPC_021 <dbl>, HSPC_022 <dbl>, HSPC_023 <dbl>,
#   HSPC_024 <dbl>, HSPC_025 <dbl>, HSPC_026 <dbl>, HSPC_027 <dbl>,
#   HSPC_028 <dbl>, HSPC_030 <dbl>, HSPC_031 <dbl>, HSPC_033 <dbl>, …

We already know from last week’s work that there are no genes that are zero across all the cells (both types). If we did not know that, we would need to add |> filter(sum(c_across(Prog_001:Prog_852)) != 0)

meaning zero in all the HSPC but not zero in all the Prog

❓ How many genes are expressed only in the progenitor cells only

🎬 Now you find any genes that are expressed only in the HSPC cells.

❓ How many genes are expressed only in the HSPC cells?

❓ Do the results make sense to you in light of what you know about the biology?

🎬 Write all the genes that are expressed one cell type only to file (saved in results)

2. Prepare the data for analysis with `scran`

scran can use a matrix or a dataframe of counts but these must be log normalised counts. If using a dataframe, the columns must only contain the expression values (not the gene ids). The rows can be named to retain the gene ids.

hspc_prog is a dataframe so we will use the ensembl gene ids to name the rows and remove the gene ids from the dataframe.

🎬 Add the gene ids as the row names:

hspc_prog <- hspc_prog |>
  column_to_rownames("ensembl_gene_id")

Like DESeq2, scran needs metadata to define which columns were in which group. Instead of having this is a file, we will create a vector that indicates which column belongs to which cell type.

🎬 Create a vector that indicates which column belongs to which cell type:

n_hspc <- 701
n_prog <- 798

cell_type <- rep(c("hspc","prog"), 
                 times = c(n_hspc, n_prog))

The number of times each cell type is repeated is the number of cells of that type. Do check that the length of the cell_type vector is the same as the number of columns in the hspc_prog dataframe.

3. Differential expression analysis

🎬 Load the scran package:

Differential expression is carried out with the findMarkers() function. It takes two arguments. The first argument is the dataframe containing the data and the second argument is the vector indicating which columns are in which cell type. Make sure that the order of the cell types in the vector is appropriate for the order of the columns in the dataframe.

🎬 Run the differential expression analysis:

results_hspc_prog <- findMarkers(hspc_prog, 
                             cell_type)

The output is a list object which, rather unnecessarily, includes two dataframes. This is not really necessary as the results are the same except for the fold change having a different sign.

The dataframe results_hspc_prog$prog is log prog - log hspc (i.e.,Prog/HSPC). This means:
- Positive fold change: prog is higher than hspc
- Negative fold change: hspc is higher than prog

The results_hspc_prog$prog dataframe
	Top	summary.logFC	logFC.hspc	ensembl_gene_id
ENSMUSG00000024399	1	-4.3900364	-4.3900364	ENSMUSG00000024399
ENSMUSG00000079563	2	-4.0005511	-4.0005511	ENSMUSG00000079563
ENSMUSG00000042817	3	-3.6990453	-3.6990453	ENSMUSG00000042817
ENSMUSG00000063856	4	0.9115206	0.9115206	ENSMUSG00000063856
ENSMUSG00000024053	5	3.0351647	3.0351647	ENSMUSG00000024053
ENSMUSG00000069792	6	-3.4659544	-3.4659544	ENSMUSG00000069792

The dataframe results_hspc_prog$hspc is log hspc - log prog (i.e., HSPC/Prog). This means:
- Positive fold change: hspc is higher than prog
- Negative fold change: prog is higher than hspc

The results_hspc_prog$hspc dataframe. Notice the sign of the fold change is the other way
	Top	summary.logFC	logFC.prog	ensembl_gene_id
ENSMUSG00000024399	1	4.3900364	4.3900364	ENSMUSG00000024399
ENSMUSG00000079563	2	4.0005511	4.0005511	ENSMUSG00000079563
ENSMUSG00000042817	3	3.6990453	3.6990453	ENSMUSG00000042817
ENSMUSG00000063856	4	-0.9115206	-0.9115206	ENSMUSG00000063856
ENSMUSG00000024053	5	-3.0351647	-3.0351647	ENSMUSG00000024053
ENSMUSG00000069792	6	3.4659544	3.4659544	ENSMUSG00000069792

We only need one of these results dataframes and we will use the prog one. It does not matter which one you use but you do need keep the direction of the foldchange in mind when interpreting the results.

Notice that the dataframes are ordered by significance rather than the original gene order.

It is useful to have the normalised counts and the statistical results in one dataframe to which we will add the gene information from Ensembl. Having all the information together will make it easier to interpret the results and select genes of interest. We will need to extract the results from the list object, add the gene ids and then join with the normalised counts.

🎬 Extract the results dataframe from the list object and add the gene ids as a column:

hspc_prog_results <- data.frame(results_hspc_prog$prog, 
                                ensembl_gene_id = row.names(results_hspc_prog$prog))

🎬 Return the ensembl gene ids as a column to the normalised counts:

hspc_prog <- hspc_prog |>
  rownames_to_column(var = "ensembl_gene_id")

🎬 Merge the results dataframe with the normalised counts:

# merge the results with the normalised counts
hspc_prog_results <- hspc_prog_results |>
  left_join(hspc_prog, by = "ensembl_gene_id")

Now go to Add gene information.

Add gene information

🎄 Arabidopisis

Ensembl (Martin et al. 2023; Birney et al. 2004)is a bioinformatics project to organise all the biological information around the sequences of large genomes. The are a large number of databases and BioMart (Smedley et al. 2009) provides a consistent interface to the material. There are web-based tools to use these but the R package biomaRt (Durinck et al. 2009, 2005) gives you programmatic access making it easier to integrate information into R dataframes.

🎬 Load the biomaRt (Durinck et al. 2009, 2005) package:

library(biomaRt)

The biomaRt package includes a function to list all the available datasets

🎬 List the Ensembl “marts” available:

listEnsemblGenomes()

              biomart                        version
1       protists_mart      Ensembl Protists Genes 61
2 protists_variations Ensembl Protists Variations 61
3          fungi_mart         Ensembl Fungi Genes 61
4    fungi_variations    Ensembl Fungi Variations 61
5        metazoa_mart       Ensembl Metazoa Genes 61
6  metazoa_variations  Ensembl Metazoa Variations 61
7         plants_mart        Ensembl Plants Genes 61
8   plants_variations   Ensembl Plants Variations 61

plants_mart looks like the one we want. We can see what genomes are available with names like “Arabidopsis” in this mart using the searchDatasets() function.

🎬

searchDatasets(useEnsemblGenomes(biomart = "plants_mart"), 
               pattern = "Arabidopsis")

             dataset                         description version
4   ahalleri_eg_gene Arabidopsis halleri genes (Ahal2.2) Ahal2.2
6    alyrata_eg_gene    Arabidopsis lyrata genes (v.1.0)   v.1.0
11 athaliana_eg_gene Arabidopsis thaliana genes (TAIR10)  TAIR10

athaliana_eg_gene is the Arabidopsis thaliana genes (TAIR10) dataset we want.

🎬 Connect to the athaliana_eg_gene database in plants_mart:

ensembl <- useEnsemblGenomes(biomart = "plants_mart",
                             dataset = "athaliana_eg_gene")

🎬 See the the types of information we can retrieve:

listAttributes(mart = ensembl) |> View()

There are many (1,796!) possible bits of information (attributes) that can be obtained.

We use the getBM() function to retrieve information from the database. The filters argument is used to specified what kind of identifier we are supplying in values to retrieve information. The attributes argument is used to select the information we want to retrieve. The values argument is used to specify the identifiers. The mart argument is used to specify the connection we created.

🎬 Get the the gene name and a description. We also retreive the gene id so we can later join the information with the results:

gene_info <- getBM(filters = "ensembl_gene_id",
                   attributes = c("ensembl_gene_id",
                                  "external_gene_name",
                                  "description"),
                   values = root_results$gene_id,
                   mart = ensembl)

You should view the resulting dataframe to see what information is available. You can use glimpse() or View().

🎬 Merge the gene information with the results:

# join the gene info with the results
root_results <- root_results |>
  left_join(gene_info,
            by = join_by(gene_id == ensembl_gene_id))

🎬 Save the results to a file:

write_csv(root_results, file = "results/root_results.csv")

Now go to Look after future you.

💉 Leishmania

I got the information from TriTrypDB
which is a functional genomic resource for the Trypanosomatidae and Plasmodidae
https://tritrypdb.org/tritrypdb/app/downloads section
I downloaded the L. mexicana MHOM/GT/2001/U1103 Full GFF and extracted the gene information and saved it as leishmania_mex.xlsx

We will import this file and join it to the results dataframe.

🎬 Load the readxl (Wickham and Bryan 2023) package:

library(readxl)

🎬 Import the Xenbase gene information file:

gene_info <- read_excel("meta/leishmania_mex.xlsx")

Error in read_excel("meta/leishmania_mex.xlsx"): could not find function "read_excel"

You should view the resulting dataframe to see what information is available. You can use glimpse() or View().

🎬 Merge the gene information with the results:

# join the gene info with the results
pro_meta_results <- pro_meta_results |>
  left_join(gene_info, by = "gene_id")

Error in `left_join()`:
! Join columns in `y` must be present in the data.
✖ Problem with `gene_id`.

🎬 Save the results to a file:

write_csv(pro_meta_results, file = "results/pro_meta_results.csv")

Now go to Look after future you.

🐭 Stem cells

Ensembl (Martin et al. 2023; Birney et al. 2004)is a bioinformatics project to organise all the biological information around the sequences of large genomes. The are a large number of databases but BioMart (Smedley et al. 2009) provides a consistent interface to the material. There are web-based tools to use these but the R package biomaRt (Durinck et al. 2009, 2005) gives you programmatic access making it easier to integrate information into R dataframes

🎬 Load the biomaRt (Durinck et al. 2009, 2005) package:

library(biomaRt)

🎬 Connect to the mouse database

# Connect to the mouse database

ensembl <- useMart(biomart = "ensembl", 
                   dataset = "mmusculus_gene_ensembl")

🎬 See information we can retrieve:

# See what information we can retrieve
listAttributes(mart = ensembl) |> View()

There are many (~3000!) possible bits of information (attributes) that can be obtained.

🎬 Get the the gene name and a description. We also retrieve the gene id so we can later join the information with the results:

gene_info <- getBM(filters = "ensembl_gene_id",
                   attributes = c("ensembl_gene_id",
                                  "external_gene_name",
                                  "description"),
                   values = hspc_prog_results$ensembl_gene_id,
                   mart = ensembl)

You should view the resulting dataframe to see what information is available. You can use glimpse() or View(). Notice the dataframe returned only has 422 rows - one of the ids does not have information.

🎬 Merge the gene information with the results:

# join the gene info with the results
hspc_prog_results <- hspc_prog_results |> 
  left_join(gene_info, by = "ensembl_gene_id")

🎬 Save the results to a file:

write_csv(hspc_prog_results, file = "results/hspc_prog_results.csv")

Now go to Look after future you.

🤗 Look after future you!

🎬 Go through you script and tidy up. You might need to :

collect together library statements at the top of the script
remove code that you needed to start today but which wouldn’t be needed if you were running the script from the top (e.g., importing)
edit your comments for clarity
rename variables for consistency or clarity
remove house keeping or exploratory code
restyle code, add code section headers etc

🥳 Finished

Well Done!

Independent study following the workshop

Consolidate

The Code file

These contain all the code needed in the workshop even where it is not visible on the webpage.

The workshop.qmd file is the file I use to compile the practical. Qmd stands for Quarto markdown. It allows code and ordinary text to be interleaved to produce well-formatted reports including webpages. View the source code for this workshop using the </> Code button at the top of the page. Coding and thinking answers are marked with #---CODING ANSWER--- and #---THINKING ANSWER---

Pages made with R (R Core Team 2024), Quarto (Allaire et al. 2024), knitr (Xie 2024, 2015, 2014), kableExtra (Zhu 2021)

References

Allaire, J. J., Charles Teague, Carlos Scheidegger, Yihui Xie, and Christophe Dervieux. 2024. “Quarto.” https://doi.org/10.5281/zenodo.5960048.

Birney, Ewan, T. Daniel Andrews, Paul Bevan, Mario Caccamo, Yuan Chen, Laura Clarke, Guy Coates, et al. 2004. “An Overview of Ensembl.” Genome Research 14 (5): 925–28. https://doi.org/10.1101/gr.1860604.

Durinck, Steffen, Yves Moreau, Arek Kasprzyk, Sean Davis, Bart De Moor, Alvis Brazma, and Wolfgang Huber. 2005. “BioMart and Bioconductor: A Powerful Link Between Biological Databases and Microarray Data Analysis.” Bioinformatics 21: 3439–40.

Durinck, Steffen, Paul T. Spellman, Ewan Birney, and Wolfgang Huber. 2009. “Mapping Identifiers for the Integration of Genomic Datasets with the r/Bioconductor Package biomaRt.” Nature Protocols 4: 1184–91.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15: 550. https://doi.org/10.1186/s13059-014-0550-8.

Lun, Aaron T. L., Davis J. McCarthy, and John C. Marioni. 2016. “A Step-by-Step Workflow for Low-Level Analysis of Single-Cell RNA-Seq Data with Bioconductor.” F1000Res. 5: 2122. https://doi.org/10.12688/f1000research.9501.2.

Martin, Fergal J, M Ridwan Amode, Alisha Aneja, Olanrewaju Austine-Orimoloye, Andrey G Azov, If Barnes, Arne Becker, et al. 2023. “Ensembl 2023.” Nucleic Acids Research 51 (D1): D933–41. https://doi.org/10.1093/nar/gkac958.

R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Smedley, Damian, Syed Haider, Benoit Ballester, Richard Holland, Darin London, Gudmundur Thorisson, and Arek Kasprzyk. 2009. “BioMart Biological Queries Made Easy.” BMC Genomics 10 (1): 22. https://doi.org/10.1186/1471-2164-10-22.

Wickham, Hadley. 2023. Conflicted: An Alternative Conflict Resolution Strategy. https://conflicted.r-lib.org/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, and Jennifer Bryan. 2023. Readxl: Read Excel Files. https://readxl.tidyverse.org.

Xie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2024. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Zhu, Hao. 2021. “kableExtra: Construct Complex Table with ’Kable’ and Pipe Syntax.” https://CRAN.R-project.org/package=kableExtra.

Footnotes

Bioconductor is a project that develops and supports R packages for bioinformatics.↩︎
Bioconductor is a project that develops and supports R packages for bioinformatics.↩︎

--- title: "Workshop" subtitle: "Transcripttranscriptomics 2: Statistical Analysis" author: "Emma Rand" toc: true toc-depth: 4 toc-location: right execute: echo: true include: true error: true bibliography: ../../references.bib editor: markdown: wrap: 72 params: mouse: secretome # or: surfaceome --- # Introduction ## Session overview In the workshop, you will learn how to perform differential expression analysis on raw counts using **`DESeq2`** [@DESeq2] or on logged normalised expression values using **`scran`** [@scran] or both. # Set up ## 🎄 *Arabidopisis* 🎬 Open the `arab-88H` RStudio Project and the `cont-low-root.R` script. ## 💉 *Leishmania* 🎬 Open the `leish-88H` RStudio Project and the `pro-meta.R` script. ## 🐭 Stem cells 🎬 Open the `mice-88H` RStudio Project and the `hspc-prog.R` script. ## Everyone 🎬 Make a new folder `results` in the project directory. This is where we will save our results. 🎬 Load **`tidyverse`** [@tidyverse] You most likely have this code at the top of your script already. ```{r} library(tidyverse) ``` ``` ── Attaching core tidyverse packages ─────────────────────────────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.1.3 ✔ readr 2.1.4 ✔ forcats 1.0.0 ✔ stringr 1.5.0 ✔ ggplot2 3.4.3 ✔ tibble 3.2.1 ✔ lubridate 1.9.3 ✔ tidyr 1.3.0 ✔ purrr 1.0.2 ── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package to force all conflicts to become errors ``` Have you ever stopped to think about this message? It is telling us that there are functions in the `dplyr` package that have the same name as functions in the `stats` package and that R will use the `dplyr` version. As this is what you want, this has always been fine. It still is fine in this case. However, as you start to load more packages, you will want to know if you are using a function from a package that has the same name as a function in another loaded package. This is where the **`conflicted`** [@conflicted] package comes in. **`Conflicted`** will warn you when you are using a function that has the same name as a function in another package. You can then choose which function to use. 🎬 Load the **`conflicted`** package: ```{r} library(conflicted) ``` Instead of getting a warning every time you are using a function that has a function with the same name in another package, we can declare a preference for one function over another. This is useful for the functions you use a lot or ones where you are certain you always want to use a particular function. For example, to always use the **`dplyr`** version of `filter()` by default you can add this to the top of your script: ```{r} conflicts_prefer(dplyr::filter) ``` We will also want to ensure that we are using the `setdiff()` function from the **`GenomicRanges`** package which is used by **`DESeq2`**. ```{r} conflicts_prefer(GenomicRanges::setdiff) ``` # Import ## 🎄 *Arabidopisis* We need to import the root tissue data that were filtered to remove genes with 3 or 4 zeros and those where the total counts was less than 20. 🎬 Import the data from the `data-processed` folder. ```{r} #| echo: false #---CODING ANSWER--- root_filtered <- read_csv("data-processed/root_filtered.csv") ``` Now go to [Differential Expression Analysis](#arabidopisis-differential). ## 💉 *Leishmania* We need to import the procyclic- and metacyclic-promastigote data that were filtered to remove genes with 4, 5 or 6 zeros and those where the total counts was less than 20. 🎬 Import the data from the `data-processed` folder. ```{r} #| echo: false #---CODING ANSWER--- pro_meta_filtered <- read_csv("data-processed/pro_meta_filtered.csv") ``` Now go to [Differential Expression Analysis](#leishmania-differential). ## 🐭 Stem cells We need to import the progenitor and HSPC cell data that were combined. 🎬 Import the data from the `data-processed` folder.. ```{r} #| echo: false #---CODING ANSWER--- hspc_prog <- read_csv("data-processed/hspc_prog.csv") ``` Now go to [Differential Expression Analysis](#stem-cells-differential). # Differential Expression Analysis ## 🎄 *Arabidopisis* {#arabidopisis-differential} These are the steps we will take 1. Find the genes that are expressed in only one treatment group. 2. Create a DESeqDataSet object. This is a special object that is used by the **`DESeq2`** package 3. Prepare the normalised counts from the DESeqDataSet object. 4. Do differential expression analysis on the genes. This needs to be done on the raw counts. All but the first step are done with the **`DESeq2`** package ### 1. Genes expressed in one treatment The genes expressed in only one treatment group are those with zeros in all replicates in one group and non-zero values in all replicates in the other group. For example, those shown here: ```{r} #| echo: false # sort the data by the values in the control column root_filtered |> filter(CTR1 == 0) |> filter(CTR2 == 0) |> filter(CTR3 == 0) |> filter(CTR4 == 0) |> filter(CTR5 == 0) |> filter(CTR6 == 0) |> head(5) |> knitr::kable() ``` We will use `filter()` to find these genes. 🎬 Find the genes that are expressed only in the control group: ```{r} root_lowni_only <- root_filtered |> filter(if_all(starts_with("CTR"), \(x) x == 0), if_all(starts_with("LWR"), \(x) x > 0) ) ``` This code creates a new data frame called `root_lowni_only` that contains only the rows from root_filtered where: - All CTR columns have a value of 0 - All LWR columns have a value greater than 0 `filter()` keeps rows that match the conditions you give it. Here we give two conditions: - All columns whose names start with "CTR" must have a value of 0. - All columns whose names start with "LWR" must have a value greater than 0. The `if_all()` function applies a test to each selected column (those starting with "CTR" or "LWR") for each row. It returns `TRUE` only if all the selected columns pass the test. The bit `\(x) x == 0` is an anonymous function — a function without a name. In base R syntax, it’s equivalent to: `function(x) x == 0` Here, `x` will be the values from one column at a time, and the function checks whether they are equal to 0. We use an anonymous function because the condition needs to be applied separately to each column, but it’s simple enough that there’s no need to define a full named function elsewhere in the script. ❓ How many genes are expressed only in the low nickel group?   Note: There may have been more genes with counts in the low nickel group that were removed when we filtered on the total number of counts being less that 50. 🎬 Now you find any genes that are expressed only in the control group. ```{r} #| echo: false #---CODING ANSWER--- root_cont_only <- root_filtered |> filter(if_all(starts_with("CTR"), \(x) x > 0), if_all(starts_with("LWR"), \(x) x == 0) ) ``` ❓ How many genes are expressed only in the control group?   ❓ Do the results make sense to you in light of what you know about the biology?       🎬 Write all the genes that are expressed one group only to file (saved in `results`) ```{r} #| echo: false #---CODING ANSWER--- write_csv(root_lowni_only, "results/root_lowni_only.csv") write_csv(root_cont_only, "results/root_cont_only.csv") ``` ### 2. Create DESeqDataSet object 🎬 Load the DESeq2 package: ```{r} #| echo: false #---CODING ANSWER--- library(DESeq2) ``` A DEseqDataSet object is a custom data type that is used by **`DESeq2`**. Custom data types are common in the Bioconductor[^1] packages. They are used to store data in a way that is useful for the analysis. These data types typically have data, transformed data, metadata and experimental designs within them. [^1]: [Bioconductor](https://www.bioconductor.org/) is a project that develops and supports R packages for bioinformatics. To create a DESeqDataSet object, we need to provide three things: - The raw counts - these are in `root_filtered` - The meta data which gives information about the samples and which treatment groups they belong to - A design matrix which captures the design of the statistical model. The counts must in a *matrix* rather than a dataframe. Unlike a dataframe, a matrix has columns of all the same type. That is, it will contain only the counts. The gene ids are given as row names rather than a column. The `matrix()` function will create a matrix from a dataframe of columns of the same type and the `select()` function can be used to remove the gene ids column. 🎬 Create a matrix of the counts: ```{r} root_count_mat <- root_filtered |> select(-gene_id, -gene_name) |> as.matrix() ``` 🎬 Add the gene ids as row names to the matrix: ```{r} # add the row names to the matrix rownames(root_count_mat) <- root_filtered$gene_id ``` You might want to view the matrix (click on it in your environment pane). The metadata are in a file, [arab_meta_data.txt](meta/arab_meta_data.txt). This is a tab-delimited file. The first column is the sample name and the other columns give the "treatments". In this case, the treatments are `tissue` (with two levels) and `nickel` (with two levels). 🎬 Make a folder called `meta` and save the file to it. 🎬 Read the metadata into a dataframe: ```{r} meta <- read_table("meta/arab_meta_data.txt") ``` 🎬 Examine the resulting dataframe. We need to add the sample names as row names to the metadata dataframe. This is because the DESeqDataSet object will use the row names to match the samples in the metadata to the samples in the counts matrix. 🎬 Add the sample names as row names to the metadata dataframe: ```{r} meta <- meta |> column_to_rownames("sample_id") ``` We are dealing only with the root data so we need to remove the samples that are not in the root data. 🎬 Filter the metadata to keep only the root information: ```{r} meta_root <- meta |> filter(tissue == "root") ``` We can now create the DESeqDataSet object. The design formula describes the statistical model. You should recognise the form from previous work. The `~` can be read as "explain by" and on its right hand side are the explanatory variables. That is, the model is counts explained by nickel status. Note that: - The names of the columns in the count matrix have to exactly match the names of the rows in the metadata dataframe. They also need to be in the same order. - The names of the explanatory variables in the design formula have to match the names of columns in the metadata. 🎬 Create the DESeqDataSet object: ```{r} dds <- DESeqDataSetFromMatrix(root_count_mat, colData = meta_root, design = ~ nickel) ``` The warning "Warning: some variables in design formula are characters, converting to factors" just means that the variable type of nickel in the metadata dataframe is "char" and it has been converted into a factor type. To help you understand what the `DESeqDataSet` object we have called `dds` contains, we can look its contents The counts are in `dds@assays@data@listData[["counts"]]` and the metadata are in `dds@colData` but the easiest way to see them is to use the `counts()` and `colData()` functions from the **`DESeq2`** package. 🎬 View the counts: ```{r} #| eval: false counts(dds) |> View() ``` You should be able to see that this is the same as in `root_count_mat`. 🎬 View the column information: ```{r} colData(dds) ``` You should be able to see this is the same as in `meta_root`. ### 3. Prepare the normalised counts The normalised counts are the counts that have been transformed to account for the library size (i.e., the total number of reads in a sample) and the gene length. We have to first estimate the normalisation factors and store them in the DESeqDataSet object and then we can get the normalised counts. 🎬 Estimate the factors for normalisation and store them in the DESeqDataSet object: ```{r} dds <- estimateSizeFactors(dds) ``` 🎬 Look at the factors (just for information): ```{r} sizeFactors(dds) ``` The normalised counts will be useful to use later. To get the normalised counts we again used the `counts()` function but this time we use the `normalized=TRUE` argument. 🎬 Save the normalised to a matrix: ```{r} normalised_counts <- counts(dds, normalized = TRUE) ``` 🎬 Make a dataframe of the normalised counts, adding a column for the gene ids at the same time: ```{r} root_normalised_counts <- data.frame(normalised_counts, gene_id = row.names(normalised_counts)) ``` ### 4. Differential expression analysis We use the `DESeq()` function to do the differential expression analysis. This function fits the statistical model to the data and then uses the model to calculate the significance of the difference between the treatments. It again stores the results in the DESseqDataSet object. Note that the differential expression needs the raw (unnormalised counts) as it does its own normalisation as part of the process. 🎬 Run the differential expression analysis and store the results in the same object: ```{r} dds <- DESeq(dds) ``` The function will take only a few moments to run on this data but can take longer for bigger datasets. We need to define the contrasts we want to test. We want to test the difference between the treatments so we will define the contrast as `control` and `low_ni`. 🎬 Define the contrast: ```{r} contrast_suf <- c("nickel", "control", "low_ni") ``` Note that `nickel` is the name of the column in the metadata dataframe and `control` and `low_ni` are the names of the levels in the `nickel` column. By putting them in the order `control` , `low_ni` we are saying the fold change will be control / low_ni. This means: - positive log fold changes indicate control \> low_ni and - negative log fold changes indicates low_ni \> control. If we had put them in the order `low_ni`, `control` we would have the reverse. 🎬 Extract the results from the DESseqDataSet object: ```{r} results_suf <- results(dds, contrast = contrast_suf) ``` This will give us the log~2~ fold change, the *p*-value and the adjusted *p*-value for the comparison between the control and low_ni for each gene. 🎬 Put the results in a dataframe and add the gene ids as a column: ```{r} root_results <- data.frame(results_suf, gene_id = row.names(results_suf)) ``` It is useful to have the normalised counts and the statistical results in one dataframe. 🎬 Merge the two dataframes: ```{r} # merge the results with the normalised counts root_results <- root_normalised_counts |> left_join(root_results, by = "gene_id") ``` Now go to [Add gene information](#arabidopisis-add-info). ## 💉 *Leishmania* {#leishmania-differential} These are the steps we will take 1. Find the genes that are expressed in only one treatment group. 2. Create a DESeqDataSet object. This is a special object that is used by the **`DESeq2`** package 3. Prepare the normalised counts from the DESeqDataSet object. 4. Do differential expression analysis on the genes. This needs to be done on the raw counts. All but the first step are done with the **`DESeq2`** package ### 1. Genes expressed in one treatment The genes expressed in only one treatment group are those with zeros in all replicates in one group and non-zero values in all replicates in the other group. We will use `filter()` to find these genes. 🎬 Find the genes that are expressed only at the procyclic-promastigote stage: ```{r} pro_meta_pro_only <- pro_meta_filtered |> filter(lm_pro_1 > 0, lm_pro_2 > 0, lm_pro_3 > 0, lm_meta_1 == 0, lm_meta_2 == 0, lm_meta_2 == 0) ``` ❓ How many genes are expressed only in the procyclic-promastigote stage group?   🎬 Now you find any genes that are expressed only at the metacyclic stage ```{r} #| echo: false #---CODING ANSWER--- pro_meta_meta_only <- pro_meta_filtered |> filter(lm_pro_1 == 0, lm_pro_2 == 0, lm_pro_3 == 0, lm_meta_1 > 0, lm_meta_2 > 0, lm_meta_2 > 0) ``` ❓ How many genes are expressed only at the metacyclic stage?   ❓ Do the results make sense to you in light of what you know about the biology?     🎬 Write all the genes that are expressed one group only to file (saved in `results`) ```{r} #| echo: false #| eval: false #---CODING ANSWER--- # If we had any genes expressed in only one group we would write # them to a file like this. write_csv(pro_meta_pro_only, "results/pro_meta_pro_only.csv") write_csv(pro_meta_meta_only, "results/pro_meta_meta_only.csv") ``` ### 2. Create DESeqDataSet object 🎬 Load the DESeq2 package: ```{r} #| echo: false #| eval: false #---CODING ANSWER--- library(DESeq2) ``` A DEseqDataSet object is a custom data type that is used by **`DESeq2`**. Custom data types are common in the Bioconductor[^2] packages. They are used to store data in a way that is useful for the analysis. These data types typically have data, transformed data, metadata and experimental designs within them. [^2]: [Bioconductor](https://www.bioconductor.org/) is a project that develops and supports R packages for bioinformatics. To create a DESeqDataSet object, we need to provide three things: - The raw counts - these are in `pro_meta_filtered` - The meta data which gives information about the samples and which treatment groups they belong to - A design matrix which captures the design of the statistical model. The counts must in a *matrix* rather than a dataframe. Unlike a dataframe, a matrix has columns of all the same type. That is, it will contain only the counts. The gene ids are given as row names rather than a column. The `matrix()` function will create a matrix from a dataframe of columns of the same type and the `select()` function can be used to remove the gene ids column. 🎬 Create a matrix of the counts: ```{r} pro_meta_count_mat <- pro_meta_filtered |> select(-gene_id) |> as.matrix() ``` 🎬 Add the gene ids as row names to the matrix: ```{r} # add the row names to the matrix rownames(pro_meta_count_mat) <- pro_meta_filtered$gene_id ``` You might want to view the matrix (click on it in your environment pane). The metadata are in a file, [leish_meta_data.txt](meta/leish_meta_data.txt). This is a tab-delimited file. The first column is the sample name and the other columns give the "treatments". In this case, the treatment is stage (with three levels). 🎬 Make a folder called `meta` and save the file to it. 🎬 Read the metadata into a dataframe: ```{r} meta <- read_table("meta/leish_meta_data.txt") ``` 🎬 Examine the resulting dataframe. We need to add the sample names as row names to the metadata dataframe. This is because the DESeqDataSet object will use the row names to match the samples in the metadata to the samples in the counts matrix. 🎬 Add the sample names as row names to the metadata dataframe: ```{r} meta <- meta |> column_to_rownames("sample_id") ``` We are dealing only with the wild data so we need to remove the samples that are not in the wild data. 🎬 Filter the metadata to keep only the procyclic and metacyclic information: ```{r} meta_pro_meta <- meta |> filter(stage != "amastigotes") ``` We can now create the DESeqDataSet object. The design formula describes the statistical model. You should recognise the form from previous work. The `~` can be read as "explain by" and on its right hand side are the explanatory variables. That is, the model is counts explained by stage status. Note that: - The names of the columns in the count matrix have to exactly match the names of the rows in the metadata dataframe. They also need to be in the same order. - The names of the explanatory variables in the design formula have to match the names of columns in the metadata. 🎬 Create the DESeqDataSet object: ```{r} dds <- DESeqDataSetFromMatrix(pro_meta_count_mat, colData = meta_pro_meta, design = ~ stage) ``` The warning "Warning: some variables in design formula are characters, converting to factors" just means that the variable type of stage in the metadata dataframe is "char" and it has been converted into a factor type. To help you understand what the `DESeqDataSet` object we have called `dds` contains, we can look its contents The counts are in `dds@assays@data@listData[["counts"]]` and the metadata are in `dds@colData` but the easiest way to see them is to use the `counts()` and `colData()` functions from the **`DESeq2`** package. 🎬 View the counts: ```{r} #| eval: false counts(dds) |> View() ``` You should be able to see that this is the same as in `pro_meta_count_mat`. 🎬 View the column information: ```{r} colData(dds) ``` You should be able to see this is the same as in `meta_pro_meta`. ### 3. Prepare the normalised counts The normalised counts are the counts that have been transformed to account for the library size (i.e., the total number of reads in a sample) and the gene length. We have to first estimate the normalisation factors and store them in the DESeqDataSet object and then we can get the normalised counts. 🎬 Estimate the factors for normalisation and store them in the DESeqDataSet object: ```{r} dds <- estimateSizeFactors(dds) ``` 🎬 Look at the factors (just for information): ```{r} sizeFactors(dds) ``` The normalised counts will be useful to use later. To get the normalised counts we again used the `counts()` function but this time we use the `normalized=TRUE` argument. 🎬 Save the normalised to a matrix: ```{r} normalised_counts <- counts(dds, normalized = TRUE) ``` 🎬 Make a dataframe of the normalised counts, adding a column for the gene ids at the same time: ```{r} pro_meta_normalised_counts <- data.frame(normalised_counts, gene_id = row.names(normalised_counts)) ``` ### 4. Differential expression analysis We use the `DESeq()` function to do the differential expression analysis. This function fits the statistical model to the data and then uses the model to calculate the significance of the difference between the treatments. It again stores the results in the DESseqDataSet object. Note that the differential expression needs the raw (unnormalised counts) as it does its own normalisation as part of the process. 🎬 Run the differential expression analysis and store the results in the same object: ```{r} dds <- DESeq(dds) ``` The function will take only a few moments to run on this data but can take longer for bigger datasets. We need to define the contrasts we want to test. We want to test the difference between the treatments so we will define the contrast as `procyclic` and `metacyclic`. 🎬 Define the contrast: ```{r} contrast_pro_meta <- c("stage", "procyclic", "metacyclic") ``` Note that `stage` is the name of the column in the metadata dataframe and `procyclic` and `metacyclic` are the names of the levels in the `stage` column. By putting them in the order `procyclic` , `metacyclic` we are saying the fold change will be procyclic / metacyclic. This means: - positive log fold changes indicate procyclic \> metacyclic and - negative log fold changes indicates metacyclic \> procyclic. If we had put them in the order `metacyclic`, `procyclic` we would have the reverse. 🎬 Extract the results from the DESseqDataSet object: ```{r} results_pro_meta <- results(dds, contrast = contrast_pro_meta) ``` This will give us the log~2~ fold change, the *p*-value and the adjusted *p*-value for the comparison between procyclic and metacyclic stage for each gene 🎬 Put the results in a dataframe and add the gene ids as a column: ```{r} pro_meta_results <- data.frame(results_pro_meta, gene_id = row.names(results_pro_meta)) ``` It is useful to have the normalised counts and the statistical results in one dataframe. 🎬 Merge the two dataframes: ```{r} # merge the results with the normalised counts pro_meta_results <- pro_meta_normalised_counts |> left_join(pro_meta_results, by = "gene_id") ``` Now go to [Add gene information](#leishmania-add-info). ## 🐭 Stem cells {#stem-cells-differential} These are the steps we will take 1. Find the genes that are expressed in only one cell type (the prog or the hspc). 2. Prepare the data for differential expression analysis with the **`scran`** package. 3. Do differential expression analysis on the genes using the **`scran`** package. This needs to be done on the logged normalised counts. ### 1. Genes expressed in one cell type The genes expressed in only cell type are those with zeros in all the cells of the other type. We can find these by summing the expression values for each gene across the cells of one type and filtering for those that are zero. To do row wise aggregates such as the sum across rows we can use the `rowwise()` function. `c_across()` allows us to use the colon notation to select columns. This is very useful when you have a lot of columns because it would be annoying to have to list all of them. `HSPC_001:HSPC_852` means all the columns from `HSPC_001` to `HSPC_852`. 🎬 Find the genes expressed only in progenitor cells (*i.e.*, those that are 0 in every HSPC cell): ```{r} hspc_prog |> rowwise() |> filter(sum(c_across(HSPC_001:HSPC_852)) == 0) ``` We already know from last week's work that there are no genes that are zero across all the cells (both types). If we did not know that, we would need to add `|> filter(sum(c_across(Prog_001:Prog_852)) != 0)` meaning zero in all the HSPC but not zero in all the Prog ❓ How many genes are expressed only in the progenitor cells only   🎬 Now you find any genes that are expressed only in the HSPC cells. ```{r} #| include: false #---CODING ANSWER--- hspc_prog |> rowwise() |> filter(sum(c_across(Prog_001:Prog_852)) == 0) ``` ❓ How many genes are expressed only in the HSPC cells?   ❓ Do the results make sense to you in light of what you know about the biology?        🎬 Write all the genes that are expressed one cell type only to file (saved in `results`) ```{r} #| echo: false #| eval: false #---CODING ANSWER--- # If we had any genes expressed in only one cell type we would write # them to a file like this. write_csv(hspc_prog_hspc_only, "results/hspc_prog_hspc_only.csv") write_csv(hspc_prog_prog_only, "results/hspc_prog_prog_only.csv") ``` ### 2. Prepare the data for analysis with **`scran`** **`scran`** can use a matrix or a dataframe of counts but these must be log normalised counts. If using a dataframe, the columns must only contain the expression values (not the gene ids). The rows can be named to retain the gene ids. `hspc_prog` is a dataframe so we will use the ensembl gene ids to name the rows and remove the gene ids from the dataframe. 🎬 Add the gene ids as the row names: ```{r} hspc_prog <- hspc_prog |> column_to_rownames("ensembl_gene_id") ``` Like **`DESeq2`**, **`scran`** needs metadata to define which columns were in which group. Instead of having this is a file, we will create a vector that indicates which column belongs to which cell type. 🎬 Create a vector that indicates which column belongs to which cell type: ```{r} n_hspc <- 701 n_prog <- 798 cell_type <- rep(c("hspc","prog"), times = c(n_hspc, n_prog)) ``` The number of times each cell type is repeated is the number of cells of that type. Do check that the length of the `cell_type` vector is the same as the number of columns in the `hspc_prog` dataframe. ### 3. Differential expression analysis 🎬 Load the **`scran`** package: ```{r} #| echo: false library(scran) ``` Differential expression is carried out with the `findMarkers()` function. It takes two arguments. The first argument is the dataframe containing the data and the second argument is the vector indicating which columns are in which cell type. Make sure that the order of the cell types in the vector is appropriate for the order of the columns in the dataframe. 🎬 Run the differential expression analysis: ```{r} results_hspc_prog <- findMarkers(hspc_prog, cell_type) ``` The output is a list object which, rather unnecessarily, includes two dataframes. This is not really necessary as the results are the same except for the fold change having a different sign. - The dataframe `results_hspc_prog$prog` is log prog - log hspc (*i.e.*,Prog/HSPC). This means: - Positive fold change: prog is higher than hspc - Negative fold change: hspc is higher than prog ```{r} #| echo: false data.frame(results_hspc_prog$prog, ensembl_gene_id = row.names(results_hspc_prog$prog)) |> head() |> knitr::kable(cap = "The results_hspc_prog$prog dataframe") ``` - The dataframe `results_hspc_prog$hspc` is log hspc - log prog (*i.e.*, HSPC/Prog). This means: - Positive fold change: hspc is higher than prog - Negative fold change: prog is higher than hspc ```{r} #| echo: false data.frame(results_hspc_prog$hspc, ensembl_gene_id = row.names(results_hspc_prog$hspc)) |> head() |> knitr::kable(cap = "The results_hspc_prog$hspc dataframe. Notice the sign of the fold change is the other way") ``` We only need one of these results dataframes and we will use the `prog` one. It does not matter which one you use but you do need keep the direction of the foldchange in mind when interpreting the results. Notice that the dataframes are ordered by significance rather than the original gene order. It is useful to have the normalised counts and the statistical results in one dataframe to which we will add the gene information from Ensembl. Having all the information together will make it easier to interpret the results and select genes of interest. We will need to extract the results from the list object, add the gene ids and then join with the normalised counts. 🎬 Extract the results dataframe from the list object and add the gene ids as a column: ```{r} hspc_prog_results <- data.frame(results_hspc_prog$prog, ensembl_gene_id = row.names(results_hspc_prog$prog)) ``` 🎬 Return the ensembl gene ids as a column to the normalised counts: ```{r} hspc_prog <- hspc_prog |> rownames_to_column(var = "ensembl_gene_id") ``` 🎬 Merge the results dataframe with the normalised counts: ```{r} # merge the results with the normalised counts hspc_prog_results <- hspc_prog_results |> left_join(hspc_prog, by = "ensembl_gene_id") ``` Now go to [Add gene information](#stem-add-info). # Add gene information ## 🎄 *Arabidopisis* {#arabidopisis-add-info} [Ensembl](https://www.ensembl.org/index.html) [@martin2023; @birney2004]is a bioinformatics project to organise all the biological information around the sequences of large genomes. The are a large number of databases and [BioMart](https://www.ensembl.org/info/data/biomart/index.html) [@smedley2009] provides a consistent interface to the material. There are web-based tools to use these but the R package **`biomaRt`** [@biomaRt1; @biomaRt2] gives you programmatic access making it easier to integrate information into R dataframes. 🎬 Load the **`biomaRt`** [@biomaRt1; @biomaRt2] package: ```{r} library(biomaRt) ``` The **`biomaRt`** package includes a function to list all the available datasets 🎬 List the Ensembl "marts" available: ```{r} listEnsemblGenomes() ``` `plants_mart` looks like the one we want. We can see what genomes are available with names like "Arabidopsis" in this mart using the `searchDatasets()` function. 🎬 ```{r} searchDatasets(useEnsemblGenomes(biomart = "plants_mart"), pattern = "Arabidopsis") ``` `athaliana_eg_gene` is the Arabidopsis thaliana genes (TAIR10) dataset we want. 🎬 Connect to the `athaliana_eg_gene` database in `plants_mart`: ```{r} ensembl <- useEnsemblGenomes(biomart = "plants_mart", dataset = "athaliana_eg_gene") ``` 🎬 See the the types of information we can retrieve: ```{r} #| eval: false listAttributes(mart = ensembl) |> View() ``` There are many (1,796!) possible bits of information (attributes) that can be obtained. We use the `getBM()` function to retrieve information from the database. The `filters` argument is used to specified what kind of identifier we are supplying in `values` to retrieve information. The `attributes` argument is used to select the information we want to retrieve. The `values` argument is used to specify the identifiers. The mart argument is used to specify the connection we created. 🎬 Get the the gene name and a description. We also retreive the gene id so we can later join the information with the results: ```{r} gene_info <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "external_gene_name", "description"), values = root_results$gene_id, mart = ensembl) ``` You should view the resulting dataframe to see what information is available. You can use `glimpse()` or `View()`. 🎬 Merge the gene information with the results: ```{r} # join the gene info with the results root_results <- root_results |> left_join(gene_info, by = join_by(gene_id == ensembl_gene_id)) ``` 🎬 Save the results to a file: ```{r} write_csv(root_results, file = "results/root_results.csv") ``` ```{r} #| echo: false #---UNDER THE HOOD CODE--- write_csv(root_results, file = "../week-5/results/root_results.csv") ``` Now go to [Look after future you](#look-after-future-you). ## 💉 *Leishmania* {#leishmania-add-info} - I got the information from [TriTrypDB](https://tritrypdb.org/tritrypdb/app/)\ which is a functional genomic resource for the Trypanosomatidae and Plasmodidae - <https://tritrypdb.org/tritrypdb/app/downloads> section - I downloaded the [*L. mexicana* MHOM/GT/2001/U1103 Full GFF](https://tritrypdb.org/common/downloads/release-68/LmexicanaMHOMGT2001U1103/gff/data/TriTrypDB-68_LmexicanaMHOMGT2001U1103.gff) and extracted the gene information and saved it as [leishmania_mex.xlsx](meta/leishmania_mex.xlsx) We will import this file and join it to the results dataframe. 🎬 Load the **`readxl`** [@readxl] package: ```{r} #| eval: false library(readxl) ``` 🎬 Import the Xenbase gene information file: ```{r} gene_info <- read_excel("meta/leishmania_mex.xlsx") ``` You should view the resulting dataframe to see what information is available. You can use `glimpse()` or `View()`. 🎬 Merge the gene information with the results: ```{r} # join the gene info with the results pro_meta_results <- pro_meta_results |> left_join(gene_info, by = "gene_id") ``` 🎬 Save the results to a file: ```{r} write_csv(pro_meta_results, file = "results/pro_meta_results.csv") ``` ```{r} #| echo: false #---UNDER THE HOOD CODE--- write_csv(pro_meta_results, file = "../week-5/results/pro_meta_results.csv") ``` Now go to [Look after future you](#look-after-future-you). ## 🐭 Stem cells {#stem-cells-add-info} [Ensembl](https://www.ensembl.org/index.html) [@martin2023; @birney2004]is a bioinformatics project to organise all the biological information around the sequences of large genomes. The are a large number of databases but [BioMart](https://www.ensembl.org/info/data/biomart/index.html) [@smedley2009] provides a consistent interface to the material. There are web-based tools to use these but the R package **`biomaRt`** [@biomaRt1; @biomaRt2] gives you programmatic access making it easier to integrate information into R dataframes 🎬 Load the **`biomaRt`** [@biomaRt1; @biomaRt2] package: ```{r} #| eval: false library(biomaRt) ``` 🎬 Connect to the mouse database ```{r} # Connect to the mouse database ensembl <- useMart(biomart = "ensembl", dataset = "mmusculus_gene_ensembl") ``` 🎬 See information we can retrieve: ```{r} #| eval: false # See what information we can retrieve listAttributes(mart = ensembl) |> View() ``` There are many (\~3000!) possible bits of information (attributes) that can be obtained. We use the `getBM()` function to retrieve information from the database. The `filters` argument is used to specified what kind of identifier we are supplying in `values` to retrieve information. The `attributes` argument is used to select the information we want to retrieve. The `values` argument is used to specify the identifiers. The mart argument is used to specify the connection we created. 🎬 Get the the gene name and a description. We also retrieve the gene id so we can later join the information with the results: ```{r} gene_info <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "external_gene_name", "description"), values = hspc_prog_results$ensembl_gene_id, mart = ensembl) ``` You should view the resulting dataframe to see what information is available. You can use `glimpse()` or `View()`. Notice the dataframe returned only has `r nrow(gene_info)` rows - one of the ids does not have information. 🎬 Merge the gene information with the results: ```{r} # join the gene info with the results hspc_prog_results <- hspc_prog_results |> left_join(gene_info, by = "ensembl_gene_id") ``` 🎬 Save the results to a file: ```{r} write_csv(hspc_prog_results, file = "results/hspc_prog_results.csv") ``` ```{r} #| echo: false #---UNDER THE HOOD CODE--- write_csv(hspc_prog_results, file = "../week-5/results/hspc_prog_results.csv") ``` Now go to [Look after future you](#look-after-future-you). # 🤗 Look after future you! {#look-after-future-you} 🎬 Go through you script and tidy up. You might need to : - collect together library statements at the top of the script - remove code that you needed to start today but which wouldn't be needed if you were running the script from the top (e.g., importing) - edit your comments for clarity - rename variables for consistency or clarity - remove house keeping or exploratory code - restyle code, add code section headers etc # 🥳 Finished Well Done! # Independent study following the workshop [Consolidate](study_after_workshop.qmd) # The Code file These contain all the code needed in the workshop even where it is not visible on the webpage. The [workshop.qmd](workshop.qmd) file is the file I use to compile the practical. Qmd stands for Quarto markdown. It allows code and ordinary text to be interleaved to produce well-formatted reports including webpages. View the source code for this workshop using the `</> Code` button at the top of the page. Coding and thinking answers are marked with `#---CODING ANSWER---` and `#---THINKING ANSWER---` Pages made with R [@R-core], Quarto [@Allaire_Quarto_2024], `knitr` [@knitr1; @knitr2; @knitr3], `kableExtra` [@kableExtra] # References

Introduction

Session overview

Set up

🎄 Arabidopisis

💉 Leishmania

🐭 Stem cells

Everyone

Import

🎄 Arabidopisis

💉 Leishmania

🐭 Stem cells

Differential Expression Analysis

🎄 Arabidopisis

1. Genes expressed in one treatment

2. Create DESeqDataSet object

3. Prepare the normalised counts

4. Differential expression analysis

💉 Leishmania

1. Genes expressed in one treatment

2. Create DESeqDataSet object

3. Prepare the normalised counts

4. Differential expression analysis

🐭 Stem cells

1. Genes expressed in one cell type

2. Prepare the data for analysis with scran

3. Differential expression analysis

Add gene information

🎄 Arabidopisis

💉 Leishmania

🐭 Stem cells

🤗 Look after future you!

🥳 Finished

Independent study following the workshop

The Code file

References

Footnotes

2. Prepare the data for analysis with `scran`