In stemangiola/rpharma2020_tidytranscriptomics: Introduction to Tidy Transcriptomics
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```{css, echo = FALSE}
Formatting for polls
.poll {
background-color: #fff4d4;
}
.poll code {
background-color: #fff4d4;
}
# Background
Background information and setup instructions for this workshop can be found [here](https://stemangiola.github.io/rpharma2020_tidytranscriptomics/articles/background.html).
# Part 1 Bulk RNA-seq Core
This workshop will present how to perform analysis of RNA sequencing data following the tidy data paradigm [@wickham2014tidy]. The tidy data paradigm provides a standard way to organise data values within a dataset, where each variable is a column, each observation is a row, and data is manipulated using an easy-to-understand vocabulary. Most importantly, the data structure remains consistent across manipulation and analysis functions.
This can be achieved for RNA sequencing data with the [tidybulk](https://stemangiola.github.io/tidybulk/), [tidyHeatmap](https://stemangiola.github.io/tidyHeatmap) [@mangiola2020tidyheatmap] and [tidyverse](https://www.tidyverse.org/) [@wickham2019welcome] packages. The tidybulk package provides a tidy data structure and a modular framework for bulk transcriptional analyses. tidyHeatmap provides a tidy implementation of ComplexHeatmap. These packages are part of the tidytranscriptomics suite that introduces a tidy approach to RNA sequencing data representation and analysis
### Acknowledgements
Some of the material in Part 1 was adapted from an R for RNA sequencing workshop first run [here](http://combine-australia.github.io/2016-05-11-RNAseq/). Use of the airway and pasilla datasets was inspired by the [DESeq2 vignette](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html).
```r
knitr::include_graphics("../inst/vignettes/tidybulk_logo.png")
```{poll_1 class.source="poll"}
Poll: Will you code along and/or try out exercises during this workshop?
```{poll_2 class.source="poll"}
Poll: Do you have experience with transcriptomic analyses?
```{poll_3 class.source="poll"}
Poll: Do you have experience with tidyverse?
## Introduction
Measuring gene expression on a genome-wide scale has become common practice over the last two decades or so, with microarrays predominantly used pre-2008. With the advent of next generation sequencing technology in 2008, an increasing number of scientists use this technology to measure and understand changes in gene expression in often complex systems. As sequencing costs have decreased, using RNA sequencing to simultaneously measure the expression of tens of thousands of genes for multiple samples has never been easier. The cost of these experiments has now moved from generating the data to storing and analysing it.
There are many steps involved in analysing an RNA sequencing dataset. Sequenced reads are aligned to a reference genome, then the number of reads mapped to each gene can be counted. This results in a table of counts, which is what we perform statistical analyses on in R. While mapping and counting are important and necessary tasks, today we will be starting from the count data and showing how differential expression analysis can be performed in a friendly way using the Bioconductor package, tidybulk.
First, let’s load all the packages we will need to analyse the data. Note that you should load the *tidybulk* library after the tidyverse core packages for best integration.
```r
# load libraries
# dataset
library(airway)
# tidyverse core packages
library(tibble)
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(ggplot2)
# tidyverse-friendly packages
library(plotly)
library(ggrepel)
library(tidyHeatmap)
library(tidybulk)
Plot settings. Set the colours and theme we will use for our plots.
# Use colourblind-friendly colours
friendly_cols <- dittoSeq::dittoColors()
# Set theme
custom_theme <-
list(
scale_fill_manual(values = friendly_cols),
scale_color_manual(values = friendly_cols),
theme_bw() +
theme(
panel.border = element_blank(),
axis.line = element_line(),
panel.grid.major = element_line(size = 0.2),
panel.grid.minor = element_line(size = 0.1),
text = element_text(size = 12),
legend.position = "bottom",
strip.background = element_blank(),
axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10)),
axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10)),
axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)
)
)
Airway RNA sequencing dataset
Here we will perform our analysis using the data from the airway package. The airway data comes from the paper by [@himes2014rna]; and it includes 8 samples from human airway smooth muscle cells, from 4 cell lines. For each cell line treated (with dexamethasone) and untreated (negative control) a sample has undergone RNA sequencing and gene counts have been generated.
Setting up the data
The airway data is stored as a Bioconductor RangedSummarizedExperiment object. We will convert this object into a tidybulk tibble. A tibble is the tidyverse table format.
In this workshop we will be using the tidyverse pipe %>%. This 'pipes' the output from the command on the left into the command on the right/below. Using the pipe is not essential but it reduces the amount of code we need to write when we have multiple steps (as we'll see later). It also can make the steps clearer and easier to see. For more details on the pipe see here.
# load airway RNA sequencing data
data(airway)
# convert to tidybulk tibble
counts_airway <-
airway %>%
tidybulk()
We can type the name of the object to view.
counts_airway
The counts_airway object contains information about genes and samples, the first column has the Ensembl gene identifier, the second column has the sample identifier and the third column has the gene transcription abundance. The abundance is the number of reads aligning to the gene in each experimental sample. The remaining columns include sample-wise information. The dex column tells us whether the samples are treated or untreated and the cell column tells us what cell line they are from.
We can shorten the sample names. We can remove the SRR1039 prefix that's present in all of them, as shorter names can fit better in some of the plots we will create. We can use mutate() together with str_replace() to remove the SRR1039 string from the sample column.
counts_format <- counts_airway %>%
mutate(sample=str_remove(sample, "SRR1039"))
Adding gene symbols
We can get the gene symbols for these Ensembl gene ids using the Bioconductor annotation package for human, org.Hs.eg.db and add them as a column using mutate again.
counts_tt <- counts_format %>%
mutate(symbol = AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys = as.character(feature), keytype = "ENSEMBL", column="SYMBOL", multiVals = "first"))
counts_tt
With tidyverse, all above operations can be linked with the %>%, as shown below. This has the benefits that
- no temporary variables need to be created
- less typing is required
- the steps can be seen more clearly.
counts_tt <-
airway %>%
tidybulk() %>%
mutate(sample=str_remove(sample, "SRR1039")) %>%
mutate(symbol = AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys = as.character(feature), keytype = "ENSEMBL", column="SYMBOL", multiVals = "first"))
From this tidybulk tibble, we can perform differential expression analysis with the tidybulk package.
Filtering lowly transcribed genes
Genes with very low counts across all libraries provide little evidence for differential expression and they can interfere with some of the statistical approximations that are used later in the pipeline. They also add to the multiple testing burden when estimating false discovery rates, reducing power to detect differentially expressed genes. These genes should be filtered out prior to further analysis.
We can perform the filtering using tidybulk keep_abundant or identify_abundant. These functions can use the edgeR filterByExpr function described in [@law2016rna] to automatically identify the genes with adequate abundance for differential expression testing. By default, this will keep genes with ~10 counts in a minimum number of samples, the number of the samples in the smallest group. In this dataset the smallest group size is four (as we have four dex-treated samples vs four untreated). Alternatively, we could use identify_abundant to identify which genes are abundant or not (TRUE/FALSE), rather than just keeping the abundant ones.
# Filtering counts
counts_filtered <- counts_tt %>% keep_abundant(factor_of_interest=dex)
# take a look
counts_filtered
After running keep_abundant we have a column called .abundant containing TRUE (identify_abundant would have TRUE/FALSE).
Scaling counts to normalise
Scaling of counts, normalisation, is performed to eliminate uninteresting differences between samples due to sequencing depth or composition. A more detailed explanation can be found here. In the tidybulk package the function scale_abundance generates scaled counts, with scaling factors calculated on abundant (filtered) transcripts and applied to all transcripts. We can choose from different normalisation methods. Here we will use the default, edgeR's trimmed mean of M values (TMM), [@robinson2010scaling]. TMM normalisation (and most scaling normalisation methods) scale relative to one sample.
# Scaling counts
counts_scaled <- counts_filtered %>% scale_abundance()
# take a look
counts_scaled
After we run scale_abundance we should see some columns have been added at the end. The counts_scaled column contains the scaled counts.
We can visualise the difference of abundance densities before and after scaling. As tidybulk output is compatible with tidyverse, we can simply pipe it into standard tidyverse functions such as filter, pivot_longer and ggplot. We can also take advantage of ggplot's facet_wrap to easily create multiple plots.
counts_scaled %>%
pivot_longer(cols = c("counts", "counts_scaled"), names_to = "source", values_to = "abundance") %>%
ggplot(aes(x=abundance + 1, color=sample)) +
geom_density() +
facet_wrap(~source) +
scale_x_log10() +
custom_theme
In this dataset the distributions of the counts are not very different to each other before scaling but scaling does make the distributions more similar. If we saw a sample with a very different distribution we may need to investigate it.
As tidybulk smoothly integrates with ggplot2 and other tidyverse packages it can save on typing and make plots easier to generate. Compare the code for creating density plots with tidybulk versus standard base R below (standard code adapted from [@law2016rna]).
tidybulk
# tidybulk
airway %>%
tidybulk() %>%
keep_abundant(factor_of_interest=dex) %>%
scale_abundance() %>%
pivot_longer(cols = c("counts", "counts_scaled"), names_to = "source", values_to = "abundance") %>%
ggplot(aes(x=abundance + 1, color=sample)) +
geom_density() +
facet_wrap(~source) +
scale_x_log10() +
custom_theme
base R using edgeR
# Example code, no need to run
# Prepare data set
dgList <- SE2DGEList(airway)
group <- factor(dgList$samples$dex)
keep.exprs <- filterByExpr(dgList, group=group)
dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE]
nsamples <- ncol(dgList)
logcounts <- log2(dgList$counts)
# Setup graphics
col <- RColorBrewer::brewer.pal(nsamples, "Paired")
par(mfrow=c(1,2))
# Plot raw counts
plot(density(logcounts[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title(main="Counts")
for (i in 2:nsamples){
den <- density(logcounts[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", legend=dgList$samples$Run, text.col=col, bty="n")
# Plot scaled counts
dgList_norm <- calcNormFactors(dgList)
lcpm_n <- cpm(dgList_norm, log=TRUE)
plot(density(lcpm_n[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")
title("Counts scaled")
for (i in 2:nsamples){
den <- density(lcpm_n[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", legend=dgList_norm$samples$Run, text.col=col, bty="n")
Exploratory analyses
Dimensionality reduction
By far, one of the most important plots we make when we analyse RNA sequencing data are principal-component analysis (PCA) or multi-dimensional scaling (MDS) plots. We reduce the dimensions of the data to identify the greatest sources of variation in the data. A principal components analysis is an example of an unsupervised analysis, where we don't need to specify the groups. If your experiment is well controlled and has worked well, what we hope to see is that the greatest sources of variation in the data are the treatments/groups we are interested in. It is also an incredibly useful tool for quality control and checking for outliers. We can use the reduce_dimensions function to calculate the dimensions.
# Get principal components
counts_scal_PCA <-
counts_scaled %>%
reduce_dimensions(method="PCA")
```{poll_4 class.source="poll"}
Poll: What fraction of variance is explained by PC3?
See ?reduce_dimensions for how to get additional dimensions.
This joins the result to the counts object.
```r
# Take a look
counts_scal_PCA
For plotting, we can select just the sample-wise information with pivot_sample.
# take a look
counts_scal_PCA %>% pivot_sample()
We can now plot the reduced dimensions.
# PCA plot
counts_scal_PCA %>%
pivot_sample() %>%
ggplot(aes(x=PC1, y=PC2, colour=dex, shape=cell)) +
geom_point() +
geom_text_repel(aes(label=sample), show.legend = FALSE) +
custom_theme
The samples separate by treatment on PC1 which is what we hope to see. PC2 separates the N080611 cell line from the other samples, indicating a greater difference between that cell line and the others.
Hierarchical clustering with heatmaps
An alternative to principal component analysis for examining relationships between samples is using hierarchical clustering. Heatmaps are a nice visualisation to examine hierarchical clustering of your samples. tidybulk has a simple function we can use, keep_variable, to extract the most variable genes which we can then plot with tidyHeatmap.
counts_scaled %>%
# extract 500 most variable genes
keep_variable(.abundance = counts_scaled, top = 500) %>%
# create heatmap
heatmap(
.column = sample,
.row = feature,
.value = counts_scaled,
transform = log1p
) %>%
add_tile(dex) %>%
add_tile(cell)
In the heatmap we can see the samples cluster into two groups, treated and untreated, for three of the cell lines, and the cell line (N080611) again is further away from the others.
Tidybulk enables a simplified way of generating a clustered heatmap of variable genes. Compare the code below for tidybulk versus a base R method.
base R using edgeR
# Example code, no need to run
dgList <- SE2DGEList(airway)
group <- factor(dgList$samples$dex)
keep.exprs <- filterByExpr(dgList, group=group)
dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE]
dgList <- calcNormFactors(dgList)
logcounts <- cpm(dgList, log=TRUE)
var_genes <- apply(logcounts, 1, var)
select_var <- names(sort(var_genes, decreasing=TRUE))[1:500]
highly_variable_lcpm <- logcounts[select_var,]
colours <- c("#440154FF", "#21908CFF", "#fefada" )
col.group <- c("red","grey")[group]
gplots::heatmap.2(highly_variable_lcpm, col=colours, trace="none", ColSideColors=col.group, scale="row")
Differential expression
Now that we are happy that the data looks good, we can continue to testing for differentially expressed (DE) genes. We will use the test_differential_abundance from tidybulk which currently uses edgeR [@robinson2010edger] to perform the differential expression analysis. We give test_differential_abundance our tidybulk counts object and a formula, specifying the column that contains our groups to be compared. If all our samples were from the same cell line, and there were no additional factors contributing variance such as batch differences, we could use the formula 0 + dex. However, each treated and untreated sample is from a different cell line so we add the cell line as an additional factor 0 + dex + cell. We also provide the names of the groups we want to compare to .contrasts (e.g. .contrasts = c("dextreat - dexuntreat")).
We only have one contrast here so we omit the suffix.
counts_de <- counts_filtered %>%
test_differential_abundance(
.formula = ~ 0 + dex + cell,
.contrasts = c("dextrt - dexuntrt"),
omit_contrast_in_colnames = TRUE
)
The results will be joined to our counts for every sample.
# take a look
counts_de
If we just want a table of differentially expressed genes we can select the transcript-wise information with pivot_transcript.
# take a look
counts_de %>% pivot_transcript()
Now we have columns with our log-fold change (logFC), false-discovery rate (FDR) and probability value (p-value). logFC is log2(treated/untreated).
Tidybulk enables a simplified way of performing an RNA sequencing differential expression analysis (with the benefit of smoothly integrating with ggplot2 and other tidyverse packages). Compare the code for a tidybulk edgeR analysis versus standard edgeR below.
standard edgeR
# Example code, no need to run
dgList <- SE2DGEList(airway)
group <- factor(dgList$samples$dex)
keep.exprs <- filterByExpr(dgList, group=group)
dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE]
dgList <- calcNormFactors(dgList)
cell <- factor(dgList$samples$cell)
design <- model.matrix(~ 0 + group + cell)
dgList <- estimateDisp(dgList, design)
fit <- glmQLFit(dgList, design)
TvsU <- makeContrasts(TvsU=grouptrt-groupuntrt, levels=design)
qlf <- glmQLFTest(fit, contrast=TvsU)
Table of differentially expressed genes
We can write out our differentially expressed genes to a file that can be loaded into e.g. Excel. write_tsv will create a tab-separated file.
# save results
counts_de %>%
pivot_transcript() %>%
write_tsv("de_results.tsv")
Counting differentially expressed genes
In order to decide which genes are differentially expressed, we usually take a cut-off of 0.05 on the FDR (or adjusted P value), NOT the raw p-value. This is because we are testing many genes (multiple testing), and the chances of finding differentially expressed genes is very high when you do that many tests. Hence we need to control the false discovery rate, which is the adjusted p-value column in the results table. What this means is that if 100 genes are significant at a 5% false discovery rate, we are willing to accept that 5 will be false positives.
We can count how many differentially expressed genes there are. We'll filter on FDR 0.05.
counts_de %>%
filter(FDR < 0.05) %>%
summarise(num_de = n_distinct(feature))
```{poll_5 class.source="poll"}
Poll: How many differentially expressed transcripts are there for FDR < 0.05 if we did not include cell line in the formula?
### Extracting top differentially expressed genes
We can see the top genes by smallest p-value. We'll take a look at the top 6.
```r
topgenes <-
counts_de %>%
pivot_transcript() %>%
arrange(PValue) %>%
head(6)
topgenes
We can extract the symbols for these top genes to use in some of the plots we will make.
topgenes_symbols <- topgenes %>% pull(symbol)
# take a look
topgenes_symbols
Plots after testing for differentially expressed
Volcano plots
Volcano plots are a useful genome-wide plot for checking that the analysis looks good. Volcano plots enable us to visualise the significance of change (p-value) versus the fold change (logFC). Highly significant genes are towards the top of the plot. We can also colour significant genes (e.g. genes with false-discovery rate < 0.05)
# volcano plot, minimal
counts_de %>%
ggplot(aes(x=logFC, y=PValue, colour=FDR < 0.05)) +
geom_point() +
scale_y_continuous(trans = "log10_reverse") +
custom_theme
A more informative plot, integrating some of the packages in tidyverse.
counts_de %>%
pivot_transcript() %>%
# Subset data
mutate(significant = FDR<0.05 & abs(logFC) >=2) %>%
mutate(symbol = ifelse(symbol %in% topgenes_symbols, as.character(symbol), "")) %>%
# Plot
ggplot(aes(x = logFC, y = PValue, label=symbol)) +
geom_point(aes(color = significant, size = significant, alpha=significant)) +
geom_text_repel() +
# Custom scales
custom_theme +
scale_y_continuous(trans = "log10_reverse") +
scale_color_manual(values=c("black", "#e11f28")) +
scale_size_discrete(range = c(0, 2))
Stripcharts
Before following up on the differentially expressed genes with further lab work, it is also recommended to have a look at the expression levels of the individual samples for the genes of interest. We can use stripcharts to do this. These will help show if expression is consistent amongst replicates in the groups.
With stripcharts we can see if replicates tend to group together and how the expression compares to the other groups. We'll also add a box plot to show the distribution.
strip_chart <-
counts_scaled %>%
# extract counts for top differentially expressed genes
filter(symbol %in% topgenes_symbols) %>%
# make faceted stripchart
ggplot(aes(x = dex, y = counts_scaled + 1, fill = dex, label = sample)) +
geom_boxplot() +
geom_jitter() +
facet_wrap(~symbol) +
scale_y_log10()+
custom_theme
strip_chart
We can also easily check the raw and scaled counts for these genes.
counts_scaled %>%
# extract counts for top differentially expressed genes
filter(symbol %in% topgenes_symbols) %>%
# reshape to create column ("source") containing the raw and scaled counts
pivot_longer(
c(counts, counts_scaled),
names_to = "source",
values_to = "count"
) %>%
# make faceted stripchart
ggplot(aes(x = source, y = count + 1, fill = dex)) +
geom_boxplot() +
facet_wrap(~symbol) +
scale_y_log10() +
custom_theme
Interactive Plots
A really nice feature of using tidyverse and ggplot2 is that we can make interactive plots quite easily using the plotly package. This can be very useful for exploring what genes or samples are in the plots. We can make interactive plots directly from our ggplot2 object (strip_chart). Having label in the aes is useful to visualise the identifier of the data point (here the sample id) or other variables when we hover over the plot.
We can also specify which parameters from the aes we want to show up when we hover over the plot with tooltip.
strip_chart %>% ggplotly(tooltip = c("label", "y"))
Automatic bibliography
Tidybulk provides a handy function called get_bibliography that you can use to obtain the references for the methods used. The references are in BibTeX format and can be imported into your reference manager.
get_bibliography(counts_de)
Key Points
- RNA sequencing data can be represented and analysed in a 'tidy' way using tidybulk and the tidyverse
- With the modularity offered by piping we don't need to create variables, unless an object is used more than one. This improves robustness of the code.
- The principles of tidy transcriptomics are to interface as much as possible with commonly known manipulation and visualisation tools, rather than creating custom functions.
- Some of the key steps in an RNA sequencing analysis are (i) filtering lowly abundant transcripts, (ii) adjusting for differences in sequencing depth and composition, (iii) testing for differential expression
- Dimensionality reduction (PCA or MDS) plots are very important for exploring the data
- Density plots, volcano plots, strip-charts and heatmaps are useful visualisation tools for evaluating the hypothesis testing.
Supplementary
Some things we don't have time to cover in Part 1 of this workshop can be found in the Supplementary material.
Exercises
Try to apply what you've learned to another dataset. This dataset was generated from the pasilla package, which obtained the data from the paper by [@brooks2011conservation]. Here we provide it as a SummarizedExperiment object. The dataset has 7 samples from Drosophila (fruitfly): 3 treated with siRNA knockdown of the pasilla gene and 4 untreated controls, noted in column "condition". Some of the samples have been sequenced with paired-end sequencing and some with single-end, noted in column "type".
Load the data and create the tidybulk object with:
data("pasilla", package = "rpharma2020tidytranscriptomics")
counts_tt <-
pasilla %>%
tidybulk() %>%
mutate(symbol = AnnotationDbi::mapIds(org.Dm.eg.db::org.Dm.eg.db, keys=as.character(feature), keytype = "FLYBASE", column="SYMBOL", multiVals = "first"))
```{poll_6 class.source="poll"}
Poll:
1) What is the Fraction of Variance for PC1?
2) How many differentially expressed genes are there for treated vs untreated (FDR < 0.05)?
3) What is the FBgn id of the 10th most differentially expressed gene (by smallest P value)?
Extra
1.4 What code can generate a heatmap of variable genes (starting from `count_scaled`)?
1.5 What code can you use to visualise expression of the pasilla gene (gene id: FBgn0261552)
1.6 What code can generate an interactive volcano plot that has gene ids showing on hover?
1.7 What code can generate a heatmap of the top 100 differentially expressed genes?
# Part 2 Bulk RNA-seq Extended
## Tidybulk ADD versus GET modes
In this Part 2 we will see `action="get"` being used, so we will explain here what it is doing.
Every tidybulk function takes a tidybulk tibble as input, and
* `action="add"` outputs the new information joined to the original input data frame (default)
* `action="get"` outputs the new information with only the sample or transcript information, depending on what the analysis is
For example, with `action="add"` (default), we can add the PCA dimensions to the original data set. So we still have a row for every transcript in every sample.
```r
counts_scaled %>%
reduce_dimensions(
method = "PCA",
action = "add")
Or with action= "get" we can add the PCA dimensions to the original data set selecting just the sample-wise columns. Note that we now have just one row for every sample.
counts_scaled %>%
reduce_dimensions(
method = "PCA",
action = "get")
Comparison of differential analysis methods
tidybulk integrates several popular methods for differential transcript abundance testing: the edgeR quasi-likelihood [@chen2016reads] (tidybulk default method), edgeR likelihood ratio [@mccarthy2012differential], limma-voom [@law2014voom] and DESeq2 [@love2014moderated]. A common question researchers have is which method to choose. Mike Love, DESeq2 author has this advice in his blog.
knitr::include_graphics("../inst/vignettes/which_method.png")
tidybulk can help you decide which method (or methods) to use, as it provides an easy way to run multiple and see how they compare.
We can perform differential analysis with several methods, and the results will be added to the original dataset.
# load additional libraries
library(forcats)
library(tidygate)
library(GGally)
As before, we first pre-process the data, creating a tibble and identifying abundant genes.
pasilla_de <-
rpharma2020tidytranscriptomics::pasilla %>%
# Convert SummarizedExperiment object to tibble
tidybulk() %>%
# Add gene symbols
mutate(symbol = AnnotationDbi::mapIds(org.Dm.eg.db::org.Dm.eg.db, keys=as.character(feature), keytype = "FLYBASE", column="SYMBOL", multiVals = "first")) %>%
# Filter counts
keep_abundant(factor_of_interest=condition)
This is an example for the default method for differential abundance testing as we saw previously with the airway dataset. It uses the edgeR quasi-likelihood method.
pasilla_de %>%
# Test differential composition
test_differential_abundance(
~ condition + type,
action="get"
) %>%
# Sort by P value
arrange(PValue)
Now let's try to perform multiple methods on the same dataset.
de_all <-
pasilla_de %>%
# edgeR QLT
test_differential_abundance(
~ condition + type,
method = "edger_quasi_likelihood",
prefix = "edgerQLT_"
) %>%
# edgeR LRT
test_differential_abundance(
~ condition + type,
method = "edger_likelihood_ratio",
prefix = "edgerLR_"
) %>%
# limma-voom
test_differential_abundance(
~ condition + type,
method = "limma_voom",
prefix = "voom_"
) %>%
# DESeq2
test_differential_abundance(
~ condition + type,
method = "deseq2",
prefix = "deseq2_"
)
# take a look
de_all
Note
You may notice that the methods produce columns with different names for similar outputs. If you wish to make these consistent you can do that with tidyverse rename. For example, to rename the p value adjusted columns you could run below.
de_all %>% rename(deseq2_FDR = deseq2_padj, voom_FDR = voom_adj.P.Val)
```{poll_7 class.source="poll"}
Which method detects the most differentially abundant transcripts, p value adjusted for multiple testing < 0.05 (FDR, adj.P.Val, padj)?
We can visually compare the log fold change (logFC) of transcript abundance for the comparison of interest (treated vs untreated) for all methods. We will notice that the consistency of the logFC is really high for the methods.
```r
de_all %>%
pivot_transcript() %>%
select(edgerQLT_logFC, edgerLR_logFC, voom_logFC, deseq2_log2FoldChange, feature) %>%
ggpairs(1:4)
Similarly, we can visually compare the significance for all methods. In this case the difference is larger.
de_all %>%
pivot_transcript() %>%
select(edgerQLT_PValue, edgerLR_PValue, voom_P.Value, deseq2_pvalue, feature) %>%
ggpairs(1:4)
We can select some of the transcripts for further analysis using the tidygate package.
With tidygate, we can interactively draw gates to select points we want using gate. We specify which columns we want to plot in the scatterplot, and how many gates we want to draw. We can also specify the opacity if we want to make it easier to see overlapping points.
de_gate <-
de_all %>%
gate(
feature,
edgerQLT_PValue,
deseq2_pvalue,
opacity=0.3,
how_many_gates = 2
)
We then click to draw gates around the points we want, for example as shown in the screenshot below.
knitr::include_graphics("../inst/vignettes/comparison_different_DE_methods_gates.png")
That will add a column called gate, specifying which gate the points (transcripts) are in.
# using pre-selected gates just to render the html version of this document
de_gate <-
de_all %>%
gate(
feature,
edgerQLT_PValue,
deseq2_pvalue, gate_list = rpharma2020tidytranscriptomics::de_gate_gates
)
de_gate
We can check how many transcripts we've got in each gate.
de_gate %>%
pivot_transcript() %>%
dplyr::count(gate)
We can now select the transcripts from our two gates i.e. more significant in edgeR (gate 1) and more significant in DESeq2 (gate 2).
de_gate %>%
# Generate scaled counts for plotting
scale_abundance() %>%
# Filter for transcripts within the gates
filter(gate > 0) %>%
# Rename for clarity
mutate(gate = case_when(
gate == 1 ~ "more in edgeR",
gate == 2 ~ "more in DESeq2",
TRUE ~ gate
)) %>%
# Order the plots for the transcripts
mutate(feature = fct_reorder(feature, edgerQLT_PValue, min)) %>%
# Plot
ggplot(aes(condition, counts_scaled, color=gate)) +
geom_point() +
facet_wrap(~feature, scale="free_y") +
custom_theme
This enables us to see, for example, that DESeq2 produces a more conservative logFC statistic for the transcript FBgn0052939 .
de_gate %>%
pivot_transcript %>%
filter(feature == "FBgn0052939")%>%
select(edgerQLT_logFC, deseq2_log2FoldChange)
Cell type composition analysis
If we are sequencing tissue samples, we may want to know what cell types are present and if there are differences in expression between them. tidybulk has a deconvolve_cellularity function that can help us do this.
For this example we will use a subset of the breast cancer dataset from The Cancer Genome Atlas (TCGA).
BRCA_tidy <-
rpharma2020tidytranscriptomics::BRCA %>%
tidybulk(patient, transcript, count)
BRCA_tidy
With tidybulk, we can easily infer the proportions of cell types within a tissue using one of several published methods (Cibersort [@newman2015robust], EPIC [@racle2017simultaneous] and llsr [@abbas2009deconvolution]). Here we will use Cibersort which provides a default signature called LM22 to define the cell types. LM22 contains 547 genes that identify 22 human immune cell types.
BRCA_cell_type <-
BRCA_tidy %>%
deconvolve_cellularity(action="get")
BRCA_cell_type
Cell type proportions are added to the tibble as new columns. The prefix makes it easy to reshape the data frame if needed, for visualisation or further analyses.
BRCA_cell_type_long <-
BRCA_cell_type %>%
# Reshape
pivot_longer(
contains("cibersort"),
names_prefix = "cibersort: ",
names_to = "cell_type",
values_to = "proportion"
)
BRCA_cell_type_long
We can plot the proportions of immune cell types for each patient.
BRCA_cell_type_long %>%
# Plot proportions
ggplot(aes(x=patient, y=proportion, fill=cell_type)) +
geom_bar(stat = "identity") +
custom_theme
We can visualise the similarity of the tissue composition for the patients by performing a dimensionality reduction on cell type and proportion (rather than on transcript and counts as we did previously).
BRCA_cell_type_long %>%
# Filter cell types with proportion zero in all patients
group_by(cell_type) %>%
filter(sum(proportion) > 0) %>%
ungroup() %>%
reduce_dimensions(
patient,
cell_type,
proportion,
method="PCA",
action="get"
) %>%
ggplot(aes(PC1, PC2, label=patient)) +
geom_point(color="red") +
ggrepel::geom_text_repel(size=3) +
custom_theme
```{poll_8 class.source="poll"}
Poll: What is the most abundant cell type overall in BRCA samples?
We can also perform differential tissue composition analyses, similar to how we performed differential transcript abundance analyses. We use tidybulk's `test_differential_cellularity` and can perform our analyses using a known factor of interest, such as tumour subtype, or using survival data. Here we use survival data available from TCGA [@liu2018integrated].
```r
library(survival)
BRCA_tidy_survival <-
BRCA_tidy %>%
test_differential_cellularity(Surv(time, event_occurred) ~ . ) %>%
arrange(p.value)
BRCA_tidy_survival %>%
dplyr::select(.cell_type, p.value, everything())
We can visualise the proportions for the cell types most associated with survival.
BRCA_tidy_survival %>%
dplyr::slice(1:2) %>%
unnest(cell_type_proportions) %>%
ggplot(aes(time, .proportion, color = factor(event_occurred))) +
geom_point() +
facet_wrap(~ .cell_type) +
scale_x_log10() +
scale_y_continuous(trans = "logit") +
custom_theme
Key Points
tidybulk allows streamlined multi-method analysestidygate allows the selection of arbitrary points in a two-dimensional plot, and add the gate information to the input tibbletidybulk allow easy analyses of cell type composition- Testing for differences in tissue composition between samples is analogous to the testing for differences in transcript abundance
Supplementary
Some things we don't have time to cover in Part 2 of this workshop can be found in the Supplementary material.
Q & A
Part 3 Single-cell RNA-seq
[Background information] (https://stemangiola.github.io/rpharma2020_tidytranscriptomics/articles/background.html#differences-between-bulk-and-single-cell-rna-sequencing-1)
In Part 2 we showed how we can study the cell-type composition of a biological sample using bulk RNA sequencing. Single cell sequencing enables a more direct estimation of cell-type composition and gives greater resolution. For bulk RNA sequencing we need to infer the cell types using the abundance of transcripts in the whole sample, with single-cell RNA sequencing we can directly measure the transcripts in each cell and then classify the cells into cell types.
Seurat is a very popular analysis toolkit for single cell RNA sequencing data [@butler2018integrating; @stuart2019comprehensive] .
tidyseurat provides a bridge between the Seurat single-cell package [@butler2018integrating; @stuart2019comprehensive] and the tidyverse [@wickham2019welcome]. It creates an invisible layer that enables viewing the
Seurat object as a tidyverse tibble, and provides Seurat-compatible dplyr, tidyr, ggplot and plotly functions.
# load additional libraries
library(rpharma2020tidytranscriptomics)
library(dplyr)
library(purrr)
library(stringr)
library(Seurat)
library(SingleR)
library(tidyseurat)
Create tidyseurat
This is a seurat object but it is evaluated as tibble. So it is fully compatible both with Seurat and tidyverse APIs.
pbmc_small_tidy <- rpharma2020tidytranscriptomics::pbmc_small %>% tidy()
It looks like a tibble
pbmc_small_tidy
But it is a Seurat object after all
Assays(pbmc_small_tidy)
Polish the data
We can interact with our object as we do with any tibble. In this case we want to polish an annotation column.
pbmc_small_tidy_clean <-
pbmc_small_tidy %>%
# Clean groups
mutate(groups = groups %>% str_remove("^g")) %>%
# Extract sample
extract(file, "sample", "../data/sample([a-z0-9]+)/outs.+")
pbmc_small_tidy_clean
Preprocess the dataset
pbmc_small_scaled <-
pbmc_small_tidy_clean %>%
SCTransform(verbose = FALSE) %>%
FindVariableFeatures(verbose = FALSE)
pbmc_small_scaled
Reduce dimensions
Beside PCA which is a linear dimensionality reduction, we can apply neighbour aware methods such as UMAP, to better define locally similar cells. We can calculate the first 3 UMAP dimensions using the Seurat framework.
pbmc_small_UMAP <-
pbmc_small_scaled %>%
RunPCA(verbose = FALSE) %>%
RunUMAP(reduction = "pca", dims = 1:15, n.components = 3L)
And we can plot them using 3D plot using plotly.
pbmc_small_UMAP %>%
plot_ly(
x = ~`UMAP_1`,
y = ~`UMAP_2`,
z = ~`UMAP_3`,
colors = friendly_cols[1:4]
)
knitr::include_graphics("../inst/vignettes/plotly.png")
Identify clusters
We proceed with cluster identification with Seurat.
pbmc_small_cluster <-
pbmc_small_UMAP %>%
FindNeighbors(verbose = FALSE) %>%
FindClusters(method = "igraph", verbose = FALSE)
pbmc_small_cluster
Now we can interrogate the object as if it was a regular tibble data frame.
pbmc_small_cluster %>%
count(sample, groups, seurat_clusters)
Manual cell type classification
We can identify cluster markers using Seurat.
# Identify top 10 markers per cluster
markers <-
pbmc_small_cluster %>%
FindAllMarkers(only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) %>%
group_by(cluster) %>%
top_n(10, avg_logFC)
markers
# Plot heatmap
pbmc_small_cluster %>%
DoHeatmap(
features = markers$gene,
group.colors = friendly_cols
)
Automatic cell type classification
We can infer cell type identities using SingleR [@aran2019reference] and manipulate the output using tidyverse. SingleR accepts any log-normalised transcript abundance matrix
# Get cell type reference data
hpca <- HumanPrimaryCellAtlasData()
# Infer cell identities
cell_type_df <-
# extracting counts from Seurat object
pbmc_small_cluster@assays[["SCT"]]@counts %>%
log1p() %>%
# SingleR
SingleR(
ref = hpca,
labels = hpca$label.main,
method = "cluster",
clusters = pbmc_small_cluster %>% pull(seurat_clusters)
) %>%
# Formatting results
as.data.frame() %>%
as_tibble(rownames = "seurat_clusters") %>%
select(seurat_clusters, first.labels)
# Join UMAP and cell type info
pbmc_small_cell_type <-
pbmc_small_cluster %>%
left_join(
rpharma2020tidytranscriptomics::cell_type_df,
by = "seurat_clusters"
)
# Reorder columns
pbmc_small_cell_type %>%
select(cell, first.labels, everything())
We can easily summarise the results. For example, we can see how cell type classification overlaps with cluster classification.
pbmc_small_cell_type %>%
count(seurat_clusters, first.labels)
Nested analyses
A powerful tool we can use with tidyseurat is nest. We can easily perform independent analyses on subsets of the dataset. First we classify cell types in lymphoid and myeloid; then, nest based on the new classification
pbmc_small_nested <-
pbmc_small_cell_type %>%
filter(first.labels != "Platelets") %>%
mutate(cell_class = if_else(`first.labels` %in% c("Macrophage", "Monocyte"), "myeloid", "lymphoid")) %>%
nest(data = -cell_class)
pbmc_small_nested
Now we can independently for the lymphoid and myeloid subsets (i) find variable features, (ii) reduce dimensions, and (iii) cluster using both tidyverse and Seurat seamlessly.
pbmc_small_nested_reanalysed <-
pbmc_small_nested %>%
mutate(data = map(
data, ~ .x %>%
FindVariableFeatures(verbose = FALSE) %>%
RunPCA(npcs = 10, verbose = FALSE) %>%
FindNeighbors(verbose = FALSE) %>%
FindClusters(method = "igraph", verbose = FALSE)
))
pbmc_small_nested_reanalysed
Key Points
- Some basic steps of a single-cell RNA sequencing analysis are dimensionality reduction, cluster identification and cell type classification
tidyseurat is an invisible layer that operates on a Seurat object and enables us to visualise and manipulate data as if it were a tidy data frame.tidyseurat object is a Seurat object so it can be used with any Seurat compatible method
Contributing
If you want to suggest improvements for this workshop or ask questions, you can do so as described here.
Reproducibility
Record package and version information with sessionInfo
sessionInfo()
References
stemangiola/rpharma2020_tidytranscriptomics documentation built on Oct. 9, 2020, 9:41 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
```{css, echo = FALSE}
Formatting for polls
.poll { background-color: #fff4d4; } .poll code { background-color: #fff4d4; }
# Background Background information and setup instructions for this workshop can be found [here](https://stemangiola.github.io/rpharma2020_tidytranscriptomics/articles/background.html). # Part 1 Bulk RNA-seq Core This workshop will present how to perform analysis of RNA sequencing data following the tidy data paradigm [@wickham2014tidy]. The tidy data paradigm provides a standard way to organise data values within a dataset, where each variable is a column, each observation is a row, and data is manipulated using an easy-to-understand vocabulary. Most importantly, the data structure remains consistent across manipulation and analysis functions. This can be achieved for RNA sequencing data with the [tidybulk](https://stemangiola.github.io/tidybulk/), [tidyHeatmap](https://stemangiola.github.io/tidyHeatmap) [@mangiola2020tidyheatmap] and [tidyverse](https://www.tidyverse.org/) [@wickham2019welcome] packages. The tidybulk package provides a tidy data structure and a modular framework for bulk transcriptional analyses. tidyHeatmap provides a tidy implementation of ComplexHeatmap. These packages are part of the tidytranscriptomics suite that introduces a tidy approach to RNA sequencing data representation and analysis ### Acknowledgements Some of the material in Part 1 was adapted from an R for RNA sequencing workshop first run [here](http://combine-australia.github.io/2016-05-11-RNAseq/). Use of the airway and pasilla datasets was inspired by the [DESeq2 vignette](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html). ```r knitr::include_graphics("../inst/vignettes/tidybulk_logo.png")
```{poll_1 class.source="poll"} Poll: Will you code along and/or try out exercises during this workshop?
```{poll_2 class.source="poll"}
Poll: Do you have experience with transcriptomic analyses?
```{poll_3 class.source="poll"} Poll: Do you have experience with tidyverse?
## Introduction Measuring gene expression on a genome-wide scale has become common practice over the last two decades or so, with microarrays predominantly used pre-2008. With the advent of next generation sequencing technology in 2008, an increasing number of scientists use this technology to measure and understand changes in gene expression in often complex systems. As sequencing costs have decreased, using RNA sequencing to simultaneously measure the expression of tens of thousands of genes for multiple samples has never been easier. The cost of these experiments has now moved from generating the data to storing and analysing it. There are many steps involved in analysing an RNA sequencing dataset. Sequenced reads are aligned to a reference genome, then the number of reads mapped to each gene can be counted. This results in a table of counts, which is what we perform statistical analyses on in R. While mapping and counting are important and necessary tasks, today we will be starting from the count data and showing how differential expression analysis can be performed in a friendly way using the Bioconductor package, tidybulk. First, let’s load all the packages we will need to analyse the data. Note that you should load the *tidybulk* library after the tidyverse core packages for best integration. ```r # load libraries # dataset library(airway) # tidyverse core packages library(tibble) library(dplyr) library(tidyr) library(readr) library(stringr) library(ggplot2) # tidyverse-friendly packages library(plotly) library(ggrepel) library(tidyHeatmap) library(tidybulk)
Plot settings. Set the colours and theme we will use for our plots.
# Use colourblind-friendly colours friendly_cols <- dittoSeq::dittoColors() # Set theme custom_theme <- list( scale_fill_manual(values = friendly_cols), scale_color_manual(values = friendly_cols), theme_bw() + theme( panel.border = element_blank(), axis.line = element_line(), panel.grid.major = element_line(size = 0.2), panel.grid.minor = element_line(size = 0.1), text = element_text(size = 12), legend.position = "bottom", strip.background = element_blank(), axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10)), axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10)), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1) ) )
Airway RNA sequencing dataset
Here we will perform our analysis using the data from the airway package. The airway data comes from the paper by [@himes2014rna]; and it includes 8 samples from human airway smooth muscle cells, from 4 cell lines. For each cell line treated (with dexamethasone) and untreated (negative control) a sample has undergone RNA sequencing and gene counts have been generated.
Setting up the data
The airway data is stored as a Bioconductor RangedSummarizedExperiment object. We will convert this object into a tidybulk tibble. A tibble is the tidyverse table format.
In this workshop we will be using the tidyverse pipe %>%. This 'pipes' the output from the command on the left into the command on the right/below. Using the pipe is not essential but it reduces the amount of code we need to write when we have multiple steps (as we'll see later). It also can make the steps clearer and easier to see. For more details on the pipe see here.
# load airway RNA sequencing data data(airway) # convert to tidybulk tibble counts_airway <- airway %>% tidybulk()
We can type the name of the object to view.
counts_airway
The counts_airway object contains information about genes and samples, the first column has the Ensembl gene identifier, the second column has the sample identifier and the third column has the gene transcription abundance. The abundance is the number of reads aligning to the gene in each experimental sample. The remaining columns include sample-wise information. The dex column tells us whether the samples are treated or untreated and the cell column tells us what cell line they are from.
We can shorten the sample names. We can remove the SRR1039 prefix that's present in all of them, as shorter names can fit better in some of the plots we will create. We can use mutate() together with str_replace() to remove the SRR1039 string from the sample column.
counts_format <- counts_airway %>% mutate(sample=str_remove(sample, "SRR1039"))
Adding gene symbols
We can get the gene symbols for these Ensembl gene ids using the Bioconductor annotation package for human, org.Hs.eg.db and add them as a column using mutate again.
counts_tt <- counts_format %>% mutate(symbol = AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys = as.character(feature), keytype = "ENSEMBL", column="SYMBOL", multiVals = "first")) counts_tt
With tidyverse, all above operations can be linked with the %>%, as shown below. This has the benefits that
- no temporary variables need to be created
- less typing is required
- the steps can be seen more clearly.
counts_tt <- airway %>% tidybulk() %>% mutate(sample=str_remove(sample, "SRR1039")) %>% mutate(symbol = AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db, keys = as.character(feature), keytype = "ENSEMBL", column="SYMBOL", multiVals = "first"))
From this tidybulk tibble, we can perform differential expression analysis with the tidybulk package.
Filtering lowly transcribed genes
Genes with very low counts across all libraries provide little evidence for differential expression and they can interfere with some of the statistical approximations that are used later in the pipeline. They also add to the multiple testing burden when estimating false discovery rates, reducing power to detect differentially expressed genes. These genes should be filtered out prior to further analysis.
We can perform the filtering using tidybulk keep_abundant or identify_abundant. These functions can use the edgeR filterByExpr function described in [@law2016rna] to automatically identify the genes with adequate abundance for differential expression testing. By default, this will keep genes with ~10 counts in a minimum number of samples, the number of the samples in the smallest group. In this dataset the smallest group size is four (as we have four dex-treated samples vs four untreated). Alternatively, we could use identify_abundant to identify which genes are abundant or not (TRUE/FALSE), rather than just keeping the abundant ones.
# Filtering counts counts_filtered <- counts_tt %>% keep_abundant(factor_of_interest=dex) # take a look counts_filtered
After running keep_abundant we have a column called .abundant containing TRUE (identify_abundant would have TRUE/FALSE).
Scaling counts to normalise
Scaling of counts, normalisation, is performed to eliminate uninteresting differences between samples due to sequencing depth or composition. A more detailed explanation can be found here. In the tidybulk package the function scale_abundance generates scaled counts, with scaling factors calculated on abundant (filtered) transcripts and applied to all transcripts. We can choose from different normalisation methods. Here we will use the default, edgeR's trimmed mean of M values (TMM), [@robinson2010scaling]. TMM normalisation (and most scaling normalisation methods) scale relative to one sample.
# Scaling counts counts_scaled <- counts_filtered %>% scale_abundance() # take a look counts_scaled
After we run scale_abundance we should see some columns have been added at the end. The counts_scaled column contains the scaled counts.
We can visualise the difference of abundance densities before and after scaling. As tidybulk output is compatible with tidyverse, we can simply pipe it into standard tidyverse functions such as filter, pivot_longer and ggplot. We can also take advantage of ggplot's facet_wrap to easily create multiple plots.
counts_scaled %>% pivot_longer(cols = c("counts", "counts_scaled"), names_to = "source", values_to = "abundance") %>% ggplot(aes(x=abundance + 1, color=sample)) + geom_density() + facet_wrap(~source) + scale_x_log10() + custom_theme
In this dataset the distributions of the counts are not very different to each other before scaling but scaling does make the distributions more similar. If we saw a sample with a very different distribution we may need to investigate it.
As tidybulk smoothly integrates with ggplot2 and other tidyverse packages it can save on typing and make plots easier to generate. Compare the code for creating density plots with tidybulk versus standard base R below (standard code adapted from [@law2016rna]).
tidybulk
# tidybulk airway %>% tidybulk() %>% keep_abundant(factor_of_interest=dex) %>% scale_abundance() %>% pivot_longer(cols = c("counts", "counts_scaled"), names_to = "source", values_to = "abundance") %>% ggplot(aes(x=abundance + 1, color=sample)) + geom_density() + facet_wrap(~source) + scale_x_log10() + custom_theme
base R using edgeR
# Example code, no need to run # Prepare data set dgList <- SE2DGEList(airway) group <- factor(dgList$samples$dex) keep.exprs <- filterByExpr(dgList, group=group) dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE] nsamples <- ncol(dgList) logcounts <- log2(dgList$counts) # Setup graphics col <- RColorBrewer::brewer.pal(nsamples, "Paired") par(mfrow=c(1,2)) # Plot raw counts plot(density(logcounts[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title(main="Counts") for (i in 2:nsamples){ den <- density(logcounts[,i]) lines(den$x, den$y, col=col[i], lwd=2) } legend("topright", legend=dgList$samples$Run, text.col=col, bty="n") # Plot scaled counts dgList_norm <- calcNormFactors(dgList) lcpm_n <- cpm(dgList_norm, log=TRUE) plot(density(lcpm_n[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="") title("Counts scaled") for (i in 2:nsamples){ den <- density(lcpm_n[,i]) lines(den$x, den$y, col=col[i], lwd=2) } legend("topright", legend=dgList_norm$samples$Run, text.col=col, bty="n")
Exploratory analyses
Dimensionality reduction
By far, one of the most important plots we make when we analyse RNA sequencing data are principal-component analysis (PCA) or multi-dimensional scaling (MDS) plots. We reduce the dimensions of the data to identify the greatest sources of variation in the data. A principal components analysis is an example of an unsupervised analysis, where we don't need to specify the groups. If your experiment is well controlled and has worked well, what we hope to see is that the greatest sources of variation in the data are the treatments/groups we are interested in. It is also an incredibly useful tool for quality control and checking for outliers. We can use the reduce_dimensions function to calculate the dimensions.
# Get principal components counts_scal_PCA <- counts_scaled %>% reduce_dimensions(method="PCA")
```{poll_4 class.source="poll"} Poll: What fraction of variance is explained by PC3? See ?reduce_dimensions for how to get additional dimensions.
This joins the result to the counts object. ```r # Take a look counts_scal_PCA
For plotting, we can select just the sample-wise information with pivot_sample.
# take a look counts_scal_PCA %>% pivot_sample()
We can now plot the reduced dimensions.
# PCA plot counts_scal_PCA %>% pivot_sample() %>% ggplot(aes(x=PC1, y=PC2, colour=dex, shape=cell)) + geom_point() + geom_text_repel(aes(label=sample), show.legend = FALSE) + custom_theme
The samples separate by treatment on PC1 which is what we hope to see. PC2 separates the N080611 cell line from the other samples, indicating a greater difference between that cell line and the others.
Hierarchical clustering with heatmaps
An alternative to principal component analysis for examining relationships between samples is using hierarchical clustering. Heatmaps are a nice visualisation to examine hierarchical clustering of your samples. tidybulk has a simple function we can use, keep_variable, to extract the most variable genes which we can then plot with tidyHeatmap.
counts_scaled %>% # extract 500 most variable genes keep_variable(.abundance = counts_scaled, top = 500) %>% # create heatmap heatmap( .column = sample, .row = feature, .value = counts_scaled, transform = log1p ) %>% add_tile(dex) %>% add_tile(cell)
In the heatmap we can see the samples cluster into two groups, treated and untreated, for three of the cell lines, and the cell line (N080611) again is further away from the others.
Tidybulk enables a simplified way of generating a clustered heatmap of variable genes. Compare the code below for tidybulk versus a base R method.
base R using edgeR
# Example code, no need to run dgList <- SE2DGEList(airway) group <- factor(dgList$samples$dex) keep.exprs <- filterByExpr(dgList, group=group) dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE] dgList <- calcNormFactors(dgList) logcounts <- cpm(dgList, log=TRUE) var_genes <- apply(logcounts, 1, var) select_var <- names(sort(var_genes, decreasing=TRUE))[1:500] highly_variable_lcpm <- logcounts[select_var,] colours <- c("#440154FF", "#21908CFF", "#fefada" ) col.group <- c("red","grey")[group] gplots::heatmap.2(highly_variable_lcpm, col=colours, trace="none", ColSideColors=col.group, scale="row")
Differential expression
Now that we are happy that the data looks good, we can continue to testing for differentially expressed (DE) genes. We will use the test_differential_abundance from tidybulk which currently uses edgeR [@robinson2010edger] to perform the differential expression analysis. We give test_differential_abundance our tidybulk counts object and a formula, specifying the column that contains our groups to be compared. If all our samples were from the same cell line, and there were no additional factors contributing variance such as batch differences, we could use the formula 0 + dex. However, each treated and untreated sample is from a different cell line so we add the cell line as an additional factor 0 + dex + cell. We also provide the names of the groups we want to compare to .contrasts (e.g. .contrasts = c("dextreat - dexuntreat")).
We only have one contrast here so we omit the suffix.
counts_de <- counts_filtered %>% test_differential_abundance( .formula = ~ 0 + dex + cell, .contrasts = c("dextrt - dexuntrt"), omit_contrast_in_colnames = TRUE )
The results will be joined to our counts for every sample.
# take a look
counts_de
If we just want a table of differentially expressed genes we can select the transcript-wise information with pivot_transcript.
# take a look counts_de %>% pivot_transcript()
Now we have columns with our log-fold change (logFC), false-discovery rate (FDR) and probability value (p-value). logFC is log2(treated/untreated).
Tidybulk enables a simplified way of performing an RNA sequencing differential expression analysis (with the benefit of smoothly integrating with ggplot2 and other tidyverse packages). Compare the code for a tidybulk edgeR analysis versus standard edgeR below.
standard edgeR
# Example code, no need to run dgList <- SE2DGEList(airway) group <- factor(dgList$samples$dex) keep.exprs <- filterByExpr(dgList, group=group) dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE] dgList <- calcNormFactors(dgList) cell <- factor(dgList$samples$cell) design <- model.matrix(~ 0 + group + cell) dgList <- estimateDisp(dgList, design) fit <- glmQLFit(dgList, design) TvsU <- makeContrasts(TvsU=grouptrt-groupuntrt, levels=design) qlf <- glmQLFTest(fit, contrast=TvsU)
Table of differentially expressed genes
We can write out our differentially expressed genes to a file that can be loaded into e.g. Excel. write_tsv will create a tab-separated file.
# save results counts_de %>% pivot_transcript() %>% write_tsv("de_results.tsv")
Counting differentially expressed genes
In order to decide which genes are differentially expressed, we usually take a cut-off of 0.05 on the FDR (or adjusted P value), NOT the raw p-value. This is because we are testing many genes (multiple testing), and the chances of finding differentially expressed genes is very high when you do that many tests. Hence we need to control the false discovery rate, which is the adjusted p-value column in the results table. What this means is that if 100 genes are significant at a 5% false discovery rate, we are willing to accept that 5 will be false positives.
We can count how many differentially expressed genes there are. We'll filter on FDR 0.05.
counts_de %>% filter(FDR < 0.05) %>% summarise(num_de = n_distinct(feature))
```{poll_5 class.source="poll"} Poll: How many differentially expressed transcripts are there for FDR < 0.05 if we did not include cell line in the formula?
### Extracting top differentially expressed genes We can see the top genes by smallest p-value. We'll take a look at the top 6. ```r topgenes <- counts_de %>% pivot_transcript() %>% arrange(PValue) %>% head(6) topgenes
We can extract the symbols for these top genes to use in some of the plots we will make.
topgenes_symbols <- topgenes %>% pull(symbol) # take a look topgenes_symbols
Plots after testing for differentially expressed
Volcano plots
Volcano plots are a useful genome-wide plot for checking that the analysis looks good. Volcano plots enable us to visualise the significance of change (p-value) versus the fold change (logFC). Highly significant genes are towards the top of the plot. We can also colour significant genes (e.g. genes with false-discovery rate < 0.05)
# volcano plot, minimal counts_de %>% ggplot(aes(x=logFC, y=PValue, colour=FDR < 0.05)) + geom_point() + scale_y_continuous(trans = "log10_reverse") + custom_theme
A more informative plot, integrating some of the packages in tidyverse.
counts_de %>% pivot_transcript() %>% # Subset data mutate(significant = FDR<0.05 & abs(logFC) >=2) %>% mutate(symbol = ifelse(symbol %in% topgenes_symbols, as.character(symbol), "")) %>% # Plot ggplot(aes(x = logFC, y = PValue, label=symbol)) + geom_point(aes(color = significant, size = significant, alpha=significant)) + geom_text_repel() + # Custom scales custom_theme + scale_y_continuous(trans = "log10_reverse") + scale_color_manual(values=c("black", "#e11f28")) + scale_size_discrete(range = c(0, 2))
Stripcharts
Before following up on the differentially expressed genes with further lab work, it is also recommended to have a look at the expression levels of the individual samples for the genes of interest. We can use stripcharts to do this. These will help show if expression is consistent amongst replicates in the groups.
With stripcharts we can see if replicates tend to group together and how the expression compares to the other groups. We'll also add a box plot to show the distribution.
strip_chart <- counts_scaled %>% # extract counts for top differentially expressed genes filter(symbol %in% topgenes_symbols) %>% # make faceted stripchart ggplot(aes(x = dex, y = counts_scaled + 1, fill = dex, label = sample)) + geom_boxplot() + geom_jitter() + facet_wrap(~symbol) + scale_y_log10()+ custom_theme strip_chart
We can also easily check the raw and scaled counts for these genes.
counts_scaled %>% # extract counts for top differentially expressed genes filter(symbol %in% topgenes_symbols) %>% # reshape to create column ("source") containing the raw and scaled counts pivot_longer( c(counts, counts_scaled), names_to = "source", values_to = "count" ) %>% # make faceted stripchart ggplot(aes(x = source, y = count + 1, fill = dex)) + geom_boxplot() + facet_wrap(~symbol) + scale_y_log10() + custom_theme
Interactive Plots
A really nice feature of using tidyverse and ggplot2 is that we can make interactive plots quite easily using the plotly package. This can be very useful for exploring what genes or samples are in the plots. We can make interactive plots directly from our ggplot2 object (strip_chart). Having label in the aes is useful to visualise the identifier of the data point (here the sample id) or other variables when we hover over the plot.
We can also specify which parameters from the aes we want to show up when we hover over the plot with tooltip.
strip_chart %>% ggplotly(tooltip = c("label", "y"))
Automatic bibliography
Tidybulk provides a handy function called get_bibliography that you can use to obtain the references for the methods used. The references are in BibTeX format and can be imported into your reference manager.
get_bibliography(counts_de)
Key Points
- RNA sequencing data can be represented and analysed in a 'tidy' way using tidybulk and the tidyverse
- With the modularity offered by piping we don't need to create variables, unless an object is used more than one. This improves robustness of the code.
- The principles of tidy transcriptomics are to interface as much as possible with commonly known manipulation and visualisation tools, rather than creating custom functions.
- Some of the key steps in an RNA sequencing analysis are (i) filtering lowly abundant transcripts, (ii) adjusting for differences in sequencing depth and composition, (iii) testing for differential expression
- Dimensionality reduction (PCA or MDS) plots are very important for exploring the data
- Density plots, volcano plots, strip-charts and heatmaps are useful visualisation tools for evaluating the hypothesis testing.
Supplementary
Some things we don't have time to cover in Part 1 of this workshop can be found in the Supplementary material.
Exercises
Try to apply what you've learned to another dataset. This dataset was generated from the pasilla package, which obtained the data from the paper by [@brooks2011conservation]. Here we provide it as a SummarizedExperiment object. The dataset has 7 samples from Drosophila (fruitfly): 3 treated with siRNA knockdown of the pasilla gene and 4 untreated controls, noted in column "condition". Some of the samples have been sequenced with paired-end sequencing and some with single-end, noted in column "type".
Load the data and create the tidybulk object with:
data("pasilla", package = "rpharma2020tidytranscriptomics") counts_tt <- pasilla %>% tidybulk() %>% mutate(symbol = AnnotationDbi::mapIds(org.Dm.eg.db::org.Dm.eg.db, keys=as.character(feature), keytype = "FLYBASE", column="SYMBOL", multiVals = "first"))
```{poll_6 class.source="poll"}
Poll:
1) What is the Fraction of Variance for PC1?
2) How many differentially expressed genes are there for treated vs untreated (FDR < 0.05)?
3) What is the FBgn id of the 10th most differentially expressed gene (by smallest P value)?
Extra 1.4 What code can generate a heatmap of variable genes (starting from `count_scaled`)? 1.5 What code can you use to visualise expression of the pasilla gene (gene id: FBgn0261552) 1.6 What code can generate an interactive volcano plot that has gene ids showing on hover? 1.7 What code can generate a heatmap of the top 100 differentially expressed genes? # Part 2 Bulk RNA-seq Extended ## Tidybulk ADD versus GET modes In this Part 2 we will see `action="get"` being used, so we will explain here what it is doing. Every tidybulk function takes a tidybulk tibble as input, and * `action="add"` outputs the new information joined to the original input data frame (default) * `action="get"` outputs the new information with only the sample or transcript information, depending on what the analysis is For example, with `action="add"` (default), we can add the PCA dimensions to the original data set. So we still have a row for every transcript in every sample. ```r counts_scaled %>% reduce_dimensions( method = "PCA", action = "add")
Or with action= "get" we can add the PCA dimensions to the original data set selecting just the sample-wise columns. Note that we now have just one row for every sample.
counts_scaled %>% reduce_dimensions( method = "PCA", action = "get")
Comparison of differential analysis methods
tidybulk integrates several popular methods for differential transcript abundance testing: the edgeR quasi-likelihood [@chen2016reads] (tidybulk default method), edgeR likelihood ratio [@mccarthy2012differential], limma-voom [@law2014voom] and DESeq2 [@love2014moderated]. A common question researchers have is which method to choose. Mike Love, DESeq2 author has this advice in his blog.
knitr::include_graphics("../inst/vignettes/which_method.png")
tidybulk can help you decide which method (or methods) to use, as it provides an easy way to run multiple and see how they compare.
We can perform differential analysis with several methods, and the results will be added to the original dataset.
# load additional libraries library(forcats) library(tidygate) library(GGally)
As before, we first pre-process the data, creating a tibble and identifying abundant genes.
pasilla_de <- rpharma2020tidytranscriptomics::pasilla %>% # Convert SummarizedExperiment object to tibble tidybulk() %>% # Add gene symbols mutate(symbol = AnnotationDbi::mapIds(org.Dm.eg.db::org.Dm.eg.db, keys=as.character(feature), keytype = "FLYBASE", column="SYMBOL", multiVals = "first")) %>% # Filter counts keep_abundant(factor_of_interest=condition)
This is an example for the default method for differential abundance testing as we saw previously with the airway dataset. It uses the edgeR quasi-likelihood method.
pasilla_de %>% # Test differential composition test_differential_abundance( ~ condition + type, action="get" ) %>% # Sort by P value arrange(PValue)
Now let's try to perform multiple methods on the same dataset.
de_all <- pasilla_de %>% # edgeR QLT test_differential_abundance( ~ condition + type, method = "edger_quasi_likelihood", prefix = "edgerQLT_" ) %>% # edgeR LRT test_differential_abundance( ~ condition + type, method = "edger_likelihood_ratio", prefix = "edgerLR_" ) %>% # limma-voom test_differential_abundance( ~ condition + type, method = "limma_voom", prefix = "voom_" ) %>% # DESeq2 test_differential_abundance( ~ condition + type, method = "deseq2", prefix = "deseq2_" ) # take a look de_all
Note
You may notice that the methods produce columns with different names for similar outputs. If you wish to make these consistent you can do that with tidyverse
rename. For example, to rename the p value adjusted columns you could run below.de_all %>% rename(deseq2_FDR = deseq2_padj, voom_FDR = voom_adj.P.Val)
```{poll_7 class.source="poll"} Which method detects the most differentially abundant transcripts, p value adjusted for multiple testing < 0.05 (FDR, adj.P.Val, padj)?
We can visually compare the log fold change (logFC) of transcript abundance for the comparison of interest (treated vs untreated) for all methods. We will notice that the consistency of the logFC is really high for the methods. ```r de_all %>% pivot_transcript() %>% select(edgerQLT_logFC, edgerLR_logFC, voom_logFC, deseq2_log2FoldChange, feature) %>% ggpairs(1:4)
Similarly, we can visually compare the significance for all methods. In this case the difference is larger.
de_all %>% pivot_transcript() %>% select(edgerQLT_PValue, edgerLR_PValue, voom_P.Value, deseq2_pvalue, feature) %>% ggpairs(1:4)
We can select some of the transcripts for further analysis using the tidygate package.
With tidygate, we can interactively draw gates to select points we want using gate. We specify which columns we want to plot in the scatterplot, and how many gates we want to draw. We can also specify the opacity if we want to make it easier to see overlapping points.
de_gate <- de_all %>% gate( feature, edgerQLT_PValue, deseq2_pvalue, opacity=0.3, how_many_gates = 2 )
We then click to draw gates around the points we want, for example as shown in the screenshot below.
knitr::include_graphics("../inst/vignettes/comparison_different_DE_methods_gates.png")
That will add a column called gate, specifying which gate the points (transcripts) are in.
# using pre-selected gates just to render the html version of this document de_gate <- de_all %>% gate( feature, edgerQLT_PValue, deseq2_pvalue, gate_list = rpharma2020tidytranscriptomics::de_gate_gates ) de_gate
We can check how many transcripts we've got in each gate.
de_gate %>% pivot_transcript() %>% dplyr::count(gate)
We can now select the transcripts from our two gates i.e. more significant in edgeR (gate 1) and more significant in DESeq2 (gate 2).
de_gate %>% # Generate scaled counts for plotting scale_abundance() %>% # Filter for transcripts within the gates filter(gate > 0) %>% # Rename for clarity mutate(gate = case_when( gate == 1 ~ "more in edgeR", gate == 2 ~ "more in DESeq2", TRUE ~ gate )) %>% # Order the plots for the transcripts mutate(feature = fct_reorder(feature, edgerQLT_PValue, min)) %>% # Plot ggplot(aes(condition, counts_scaled, color=gate)) + geom_point() + facet_wrap(~feature, scale="free_y") + custom_theme
This enables us to see, for example, that DESeq2 produces a more conservative logFC statistic for the transcript FBgn0052939 .
de_gate %>% pivot_transcript %>% filter(feature == "FBgn0052939")%>% select(edgerQLT_logFC, deseq2_log2FoldChange)
Cell type composition analysis
If we are sequencing tissue samples, we may want to know what cell types are present and if there are differences in expression between them. tidybulk has a deconvolve_cellularity function that can help us do this.
For this example we will use a subset of the breast cancer dataset from The Cancer Genome Atlas (TCGA).
BRCA_tidy <- rpharma2020tidytranscriptomics::BRCA %>% tidybulk(patient, transcript, count) BRCA_tidy
With tidybulk, we can easily infer the proportions of cell types within a tissue using one of several published methods (Cibersort [@newman2015robust], EPIC [@racle2017simultaneous] and llsr [@abbas2009deconvolution]). Here we will use Cibersort which provides a default signature called LM22 to define the cell types. LM22 contains 547 genes that identify 22 human immune cell types.
BRCA_cell_type <- BRCA_tidy %>% deconvolve_cellularity(action="get") BRCA_cell_type
Cell type proportions are added to the tibble as new columns. The prefix makes it easy to reshape the data frame if needed, for visualisation or further analyses.
BRCA_cell_type_long <- BRCA_cell_type %>% # Reshape pivot_longer( contains("cibersort"), names_prefix = "cibersort: ", names_to = "cell_type", values_to = "proportion" ) BRCA_cell_type_long
We can plot the proportions of immune cell types for each patient.
BRCA_cell_type_long %>% # Plot proportions ggplot(aes(x=patient, y=proportion, fill=cell_type)) + geom_bar(stat = "identity") + custom_theme
We can visualise the similarity of the tissue composition for the patients by performing a dimensionality reduction on cell type and proportion (rather than on transcript and counts as we did previously).
BRCA_cell_type_long %>% # Filter cell types with proportion zero in all patients group_by(cell_type) %>% filter(sum(proportion) > 0) %>% ungroup() %>% reduce_dimensions( patient, cell_type, proportion, method="PCA", action="get" ) %>% ggplot(aes(PC1, PC2, label=patient)) + geom_point(color="red") + ggrepel::geom_text_repel(size=3) + custom_theme
```{poll_8 class.source="poll"} Poll: What is the most abundant cell type overall in BRCA samples?
We can also perform differential tissue composition analyses, similar to how we performed differential transcript abundance analyses. We use tidybulk's `test_differential_cellularity` and can perform our analyses using a known factor of interest, such as tumour subtype, or using survival data. Here we use survival data available from TCGA [@liu2018integrated]. ```r library(survival) BRCA_tidy_survival <- BRCA_tidy %>% test_differential_cellularity(Surv(time, event_occurred) ~ . ) %>% arrange(p.value) BRCA_tidy_survival %>% dplyr::select(.cell_type, p.value, everything())
We can visualise the proportions for the cell types most associated with survival.
BRCA_tidy_survival %>% dplyr::slice(1:2) %>% unnest(cell_type_proportions) %>% ggplot(aes(time, .proportion, color = factor(event_occurred))) + geom_point() + facet_wrap(~ .cell_type) + scale_x_log10() + scale_y_continuous(trans = "logit") + custom_theme
Key Points
tidybulkallows streamlined multi-method analysestidygateallows the selection of arbitrary points in a two-dimensional plot, and add the gate information to the input tibbletidybulkallow easy analyses of cell type composition- Testing for differences in tissue composition between samples is analogous to the testing for differences in transcript abundance
Supplementary
Some things we don't have time to cover in Part 2 of this workshop can be found in the Supplementary material.
Q & A
Part 3 Single-cell RNA-seq
[Background information] (https://stemangiola.github.io/rpharma2020_tidytranscriptomics/articles/background.html#differences-between-bulk-and-single-cell-rna-sequencing-1)
In Part 2 we showed how we can study the cell-type composition of a biological sample using bulk RNA sequencing. Single cell sequencing enables a more direct estimation of cell-type composition and gives greater resolution. For bulk RNA sequencing we need to infer the cell types using the abundance of transcripts in the whole sample, with single-cell RNA sequencing we can directly measure the transcripts in each cell and then classify the cells into cell types.
Seurat is a very popular analysis toolkit for single cell RNA sequencing data [@butler2018integrating; @stuart2019comprehensive] .
tidyseurat provides a bridge between the Seurat single-cell package [@butler2018integrating; @stuart2019comprehensive] and the tidyverse [@wickham2019welcome]. It creates an invisible layer that enables viewing the Seurat object as a tidyverse tibble, and provides Seurat-compatible dplyr, tidyr, ggplot and plotly functions.
# load additional libraries library(rpharma2020tidytranscriptomics) library(dplyr) library(purrr) library(stringr) library(Seurat) library(SingleR) library(tidyseurat)
Create tidyseurat
This is a seurat object but it is evaluated as tibble. So it is fully compatible both with Seurat and tidyverse APIs.
pbmc_small_tidy <- rpharma2020tidytranscriptomics::pbmc_small %>% tidy()
It looks like a tibble
pbmc_small_tidy
But it is a Seurat object after all
Assays(pbmc_small_tidy)
Polish the data
We can interact with our object as we do with any tibble. In this case we want to polish an annotation column.
pbmc_small_tidy_clean <- pbmc_small_tidy %>% # Clean groups mutate(groups = groups %>% str_remove("^g")) %>% # Extract sample extract(file, "sample", "../data/sample([a-z0-9]+)/outs.+") pbmc_small_tidy_clean
Preprocess the dataset
pbmc_small_scaled <- pbmc_small_tidy_clean %>% SCTransform(verbose = FALSE) %>% FindVariableFeatures(verbose = FALSE) pbmc_small_scaled
Reduce dimensions
Beside PCA which is a linear dimensionality reduction, we can apply neighbour aware methods such as UMAP, to better define locally similar cells. We can calculate the first 3 UMAP dimensions using the Seurat framework.
pbmc_small_UMAP <- pbmc_small_scaled %>% RunPCA(verbose = FALSE) %>% RunUMAP(reduction = "pca", dims = 1:15, n.components = 3L)
And we can plot them using 3D plot using plotly.
pbmc_small_UMAP %>% plot_ly( x = ~`UMAP_1`, y = ~`UMAP_2`, z = ~`UMAP_3`, colors = friendly_cols[1:4] )
knitr::include_graphics("../inst/vignettes/plotly.png")
Identify clusters
We proceed with cluster identification with Seurat.
pbmc_small_cluster <- pbmc_small_UMAP %>% FindNeighbors(verbose = FALSE) %>% FindClusters(method = "igraph", verbose = FALSE) pbmc_small_cluster
Now we can interrogate the object as if it was a regular tibble data frame.
pbmc_small_cluster %>% count(sample, groups, seurat_clusters)
Manual cell type classification
We can identify cluster markers using Seurat.
# Identify top 10 markers per cluster markers <- pbmc_small_cluster %>% FindAllMarkers(only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) %>% group_by(cluster) %>% top_n(10, avg_logFC) markers # Plot heatmap pbmc_small_cluster %>% DoHeatmap( features = markers$gene, group.colors = friendly_cols )
Automatic cell type classification
We can infer cell type identities using SingleR [@aran2019reference] and manipulate the output using tidyverse. SingleR accepts any log-normalised transcript abundance matrix
# Get cell type reference data hpca <- HumanPrimaryCellAtlasData() # Infer cell identities cell_type_df <- # extracting counts from Seurat object pbmc_small_cluster@assays[["SCT"]]@counts %>% log1p() %>% # SingleR SingleR( ref = hpca, labels = hpca$label.main, method = "cluster", clusters = pbmc_small_cluster %>% pull(seurat_clusters) ) %>% # Formatting results as.data.frame() %>% as_tibble(rownames = "seurat_clusters") %>% select(seurat_clusters, first.labels)
# Join UMAP and cell type info pbmc_small_cell_type <- pbmc_small_cluster %>% left_join( rpharma2020tidytranscriptomics::cell_type_df, by = "seurat_clusters" ) # Reorder columns pbmc_small_cell_type %>% select(cell, first.labels, everything())
We can easily summarise the results. For example, we can see how cell type classification overlaps with cluster classification.
pbmc_small_cell_type %>% count(seurat_clusters, first.labels)
Nested analyses
A powerful tool we can use with tidyseurat is nest. We can easily perform independent analyses on subsets of the dataset. First we classify cell types in lymphoid and myeloid; then, nest based on the new classification
pbmc_small_nested <- pbmc_small_cell_type %>% filter(first.labels != "Platelets") %>% mutate(cell_class = if_else(`first.labels` %in% c("Macrophage", "Monocyte"), "myeloid", "lymphoid")) %>% nest(data = -cell_class) pbmc_small_nested
Now we can independently for the lymphoid and myeloid subsets (i) find variable features, (ii) reduce dimensions, and (iii) cluster using both tidyverse and Seurat seamlessly.
pbmc_small_nested_reanalysed <- pbmc_small_nested %>% mutate(data = map( data, ~ .x %>% FindVariableFeatures(verbose = FALSE) %>% RunPCA(npcs = 10, verbose = FALSE) %>% FindNeighbors(verbose = FALSE) %>% FindClusters(method = "igraph", verbose = FALSE) )) pbmc_small_nested_reanalysed
Key Points
- Some basic steps of a single-cell RNA sequencing analysis are dimensionality reduction, cluster identification and cell type classification
tidyseuratis an invisible layer that operates on aSeuratobject and enables us to visualise and manipulate data as if it were a tidy data frame.tidyseuratobject is aSeurat objectso it can be used with anySeuratcompatible method
Contributing
If you want to suggest improvements for this workshop or ask questions, you can do so as described here.
Reproducibility
Record package and version information with sessionInfo
sessionInfo()
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.