In stemangiola/zhejiang2020: Introduction to Tidy Transcriptomics

options(digits=3)
options(width=90)

{width=100px}

Abstract

This is the tidy version of the material in the session RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR of this workshop

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.

Set-up

library(zhejiang2020)

# tidyverse core packages
library(tibble)
library(dplyr)
library(tidyr)
library(readr)
library(magrittr)
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(),
        text = element_text(size = 12),
        legend.position = "bottom",
        strip.background = element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1)
      )
  )

Data packaging

Organising information within a tibble (user-friendly table)

We can create a tibble with our gene-transcript abundance and sample annotation information

counts = 

  # Create the counts table
  left_join(

    # Transcript abundance
    zhejiang2020::dge_list$counts %>%
      as_tibble(rownames="entrez") %>% 
      pivot_longer(-entrez, names_to="sample", values_to="count"),

    # Sample annotation
     zhejiang2020::dge_list$samples %>%
      as_tibble(rownames="sample") %>%
      select(-lib.size, -norm.factors, -files)

  ) %>%

  # Add gene symbols
  mutate(symbol = AnnotationDbi::mapIds(
      org.Mm.eg.db::org.Mm.eg.db, 
      keys = entrez, 
      keytype = "ENTREZID", 
      column="SYMBOL", 
      multiVals = "first"
    )) %>%

  # Filter for empty symbols
  filter(!is.na(symbol)) %>%

  # Memorise key column names for streamlined analyses using `tidybulk`
  tidybulk(sample, symbol, count)


counts

Aggregate duplicate gene symbols

Our gene annotation contains 28 genes that map to multiple chromosomes, in this case we will combine all chromosome information from the multi-mapped genes

counts_aggregated = 
  counts %>%
  aggregate_duplicates(aggregation_function = median)

counts_aggregated

Data pre-processing

Identifying genes that are lowly expressed

As before, identifying lowly transcribed genes is necessary for several downstream analyses. We can specify the factor of interest for a more informed filtering. This function uses the edgeR utility filterByExpr.

counts_abundant = 
  counts_aggregated %>%
  identify_abundant(factor_of_interest = group)

counts_abundant

Scaling data for visualisation

We can compensate for technical differences in sequencing depth, scaling the data (also called normalisation). By default the TMM [@robinson2010scaling] method is used. The scaling will be calculated on the highly-transcribed genes and applied on all genes.

counts_scaled = 
  counts_abundant %>%
  scale_abundance()

counts_scaled %>% select(sample, symbol, contains("count"), everything())

We can reproduce the log-transcript-abundance density of unfiltered and filtered data (seen in the previous session of the workshop) using tidyverse tools

bind_rows(
  counts_scaled %>% mutate(label = "1.Unfiltered"),
  counts_scaled %>% filter(.abundant) %>% mutate(label = "2.Filtered")
) %>%
  ggplot(aes(count_scaled +1, color = sample)) + 
  geom_density() +
  facet_wrap(~label) +
  scale_x_log10() +
  custom_theme

We can reproduce the log-transcript-abundance density of unscaled and scaled data (seen in the previous session of the workshop) using tidyverse tools

counts_scaled %>%
  filter(.abundant) %>%

  # We reshape the data in order to build a faceted plot
  pivot_longer(contains("count"), names_to="processing", values_to="value") %>%

  # Build the plot
  ggplot(aes(sample, value + 1, fill=sample)) +
  geom_boxplot() +
  facet_wrap(~processing) +
  scale_y_log10() +
  custom_theme

Dimensionality reduction

As previously shown, we can perform dimensionality reduction to further explore our data. The reduce_dimensions function, will perform calculations only on highly-transcribed genes (i.e. .abundant == TRUE)

counts_scaled %>%
  reduce_dimensions(method = "MDS", action="get") %>%

  # We reshape the data in order to build a faceted plot
  pivot_longer(c(group, lane), names_to="annotation", values_to="value") %>%

  # Build the plot
  ggplot(aes(Dim1, Dim2, color=value, label=value)) +
  geom_text() +
  facet_wrap(~annotation) +
  custom_theme

Differential expression analysis

We can replicate the differential expression analyses using tidybulk

model.matrix(
  ~0+group+lane,
  data = pivot_sample(counts_scaled)
)

counts_test = 
  counts_scaled %>%
  test_differential_abundance(
    .formula = ~0+group+lane,
    .contrasts = c("groupBasal-groupLP", "groupBasal - groupML", "groupLP - groupML"),
    method = "limma_voom",
    action="get"
  )

counts_test

We can reproduce the fitted means (x-axis) and variances (y-axis) relationship of each gene, using the raw results from limma-voom.

counts_test %>% 
  attr("internals") %$% 
  voom %>% 
  limma::eBayes() %>% 
  limma::plotSA(main="Final model: Mean-variance trend")

Useful graphical representations of differential expression results

With ggplot2 we We can reproduce and customise the plot for the association between fold-change and average log-abundance

counts_test %>%

  filter(.abundant) %>%

  # Label significant
  mutate(significant = `adj.P.Val___groupBasal-groupLP`<0.05) %>%

  # Subset labels
  mutate(symbol = ifelse(abs(`logFC___groupBasal-groupLP`) >=8, as.character(symbol), "")) %>%

  ggplot(aes(
    x=`AveExpr___groupBasal-groupLP`, 
    y=`logFC___groupBasal-groupLP`, 
    label=symbol
  )) +
  geom_point(aes(color = significant, size = significant, alpha=significant)) +

  # Customisation
  geom_text_repel() +
  scale_color_manual(values=c("black", "#e11f28")) +
  scale_size_discrete(range = c(0, 1)) +
  theme_bw()

sessionInfo()

References

stemangiola/zhejiang2020 documentation built on Dec. 31, 2020, 7:33 a.m.