From 30b8015debc6de793cf7d8c38f039d3c3834bcc0 Mon Sep 17 00:00:00 2001
From: sjspielman <4701111+sjspielman@users.noreply.github.com>
Date: Thu, 8 Aug 2024 16:26:54 +0000
Subject: [PATCH 1/2] Live and rendered notebooks
---
RNA-seq/02-gastric_cancer_tximeta-live.Rmd | 16 +-
RNA-seq/02-gastric_cancer_tximeta.nb.html | 76 +-
.../03-gastric_cancer_exploratory-live.Rmd | 20 +-
RNA-seq/03-gastric_cancer_exploratory.nb.html | 86 +-
RNA-seq/04-nb_cell_line_tximeta.nb.html | 4 +-
RNA-seq/05-nb_cell_line_DESeq2-live.Rmd | 54 +-
RNA-seq/05-nb_cell_line_DESeq2.nb.html | 99 +-
RNA-seq/06-openpbta_heatmap-live.Rmd | 6 +-
RNA-seq/06-openpbta_heatmap.nb.html | 339 +-
.../01-intro_to_base_R-live.Rmd | 2 +-
.../01-intro_to_base_R.nb.html | 14 +-
.../02-intro_to_ggplot2-live.Rmd | 2 +-
.../02-intro_to_ggplot2.nb.html | 14 +-
.../03-intro_to_tidyverse-live.Rmd | 4 +-
.../03-intro_to_tidyverse.nb.html | 16 +-
.../01-overrepresentation_analysis-live.Rmd | 62 +-
.../01-overrepresentation_analysis.nb.html | 3459 +++++++++++++--
.../02-gene_set_enrichment_analysis-live.Rmd | 37 +-
.../02-gene_set_enrichment_analysis.nb.html | 199 +-
.../03-gene_set_variation_analysis-live.Rmd | 104 +-
.../03-gene_set_variation_analysis.nb.html | 3760 +++++++++++++----
.../01-read_filter_normalize_scRNA.nb.html | 14 +-
.../02-dataset_integration-live.Rmd | 2 +-
.../02-dataset_integration.nb.html | 22 +-
.../03-differential_expression.nb.html | 14 +-
.../04-overrepresentation_analysis.nb.html | 14 +-
.../05-gene_set_enrichment_analysis.nb.html | 18 +-
scRNA-seq/01-scRNA_quant_qc.nb.html | 6 +-
scRNA-seq/02-filtering_scRNA.nb.html | 31 +-
scRNA-seq/03-normalizing_scRNA.nb.html | 12 +-
.../04-dimension_reduction_scRNA.nb.html | 12 +-
.../05-clustering_markers_scRNA-live.Rmd | 2 +-
scRNA-seq/05-clustering_markers_scRNA.nb.html | 16 +-
scRNA-seq/06-celltype_annotation.nb.html | 10 +-
34 files changed, 6835 insertions(+), 1711 deletions(-)
diff --git a/RNA-seq/02-gastric_cancer_tximeta-live.Rmd b/RNA-seq/02-gastric_cancer_tximeta-live.Rmd
index 645c251d..e9a8bec1 100644
--- a/RNA-seq/02-gastric_cancer_tximeta-live.Rmd
+++ b/RNA-seq/02-gastric_cancer_tximeta-live.Rmd
@@ -53,8 +53,12 @@ We'll need the `quant.sf` files for all the samples in an experiment which we ha
```{r input-names}
# the quant files themselves
-sf_files <- list.files(quant_dir, recursive = TRUE, full.names = TRUE,
- pattern = "quant.sf")
+sf_files <- list.files(
+ quant_dir,
+ recursive = TRUE,
+ full.names = TRUE,
+ pattern = "quant.sf"
+)
```
```{r metadata-file}
@@ -65,7 +69,7 @@ meta_file <- file.path(data_dir, "gastric-cancer_metadata.tsv")
**Output**
```{r output-names, live = TRUE}
-# Name the output gastric-cancer_tximeta.RDS and use the directory created
+# Name the output gastric-cancer_tximeta.rds and use the directory created
# above as the rest of the path
```
@@ -98,8 +102,10 @@ sample_names
- a `names` column with the sample names
```{r names_sf_files}
-coldata <- data.frame(files = sf_files,
- names = sample_names)
+coldata <- data.frame(
+ files = sf_files,
+ names = sample_names
+)
```
We have more information about these samples stored in the metadata file that we will also want stored in `coldata`.
diff --git a/RNA-seq/02-gastric_cancer_tximeta.nb.html b/RNA-seq/02-gastric_cancer_tximeta.nb.html
index 381101a7..5ad82e07 100644
--- a/RNA-seq/02-gastric_cancer_tximeta.nb.html
+++ b/RNA-seq/02-gastric_cancer_tximeta.nb.html
@@ -3040,7 +3040,7 @@
Libraries and functions
Directories and files
-
+
# directory where the data are located
data_dir <- file.path("data", "gastric-cancer")
@@ -3050,9 +3050,7 @@
Directories and files
# create a directory to hold the tximeta results if it doesn't exist yet
txi_dir <- file.path(data_dir, "txi")
-if (!dir.exists(txi_dir)) {
- dir.create(txi_dir, recursive = TRUE)
-}
# Name the output gastric-cancer_tximeta.RDS and use the directory created
+
+
# Name the output gastric-cancer_tximeta.rds and use the directory created
# above as the rest of the path
-txi_out_file <- file.path(txi_dir, "gastric-cancer_tximeta.RDS")
diff --git a/RNA-seq/03-gastric_cancer_exploratory-live.Rmd b/RNA-seq/03-gastric_cancer_exploratory-live.Rmd
index 0277ed47..534795cf 100644
--- a/RNA-seq/03-gastric_cancer_exploratory-live.Rmd
+++ b/RNA-seq/03-gastric_cancer_exploratory-live.Rmd
@@ -43,7 +43,7 @@ data_dir <- file.path("data", "gastric-cancer")
# directory with the tximeta processed data
txi_dir <- file.path(data_dir, "txi")
-txi_file <- file.path(txi_dir, "gastric-cancer_tximeta.RDS")
+txi_file <- file.path(txi_dir, "gastric-cancer_tximeta.rds")
```
We'll create a directory to hold our plots.
@@ -77,8 +77,7 @@ First, let's read in the data we processed with `tximeta`.
We use the tissue of origin in the design formula because that will allow us to model this variable of interest.
```{r ddset}
-ddset <- DESeqDataSet(gene_summarized,
- design = ~ tissue)
+ddset <- DESeqDataSet(gene_summarized, design = ~tissue)
```
### Variance stabilizing transformation
@@ -91,7 +90,7 @@ Since different samples are usually sequenced to different depths, we want to tr
We also want to deal with the fact that genes with low counts are also likely to have higher variance (on the log2 scale), as that could bias our clustering.
To handle both of these considerations, we can calculate a Variance Stabilizing Transformation of the count data, and work with that transformed data for our analysis.
-See [this section of the `DESeq2` vignette](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization) for more on this topic.
+See [this section of the `DESeq2` vignette](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization) for more on this topic.
```{r vst}
vst_data <- vst(ddset)
@@ -112,15 +111,16 @@ Visualizing PC1 and PC2 can give us insight into how different variables (e.g.,
```{r plotPCA, live = TRUE}
-# DESeq2 built in function is called plotPCA and we want to color points by
-# tissue
+# DESeq2 built in function is called plotPCA and we want to
+# color points by tissue
```
Save the most recent plot to file with `ggsave` from `ggplot2`
```{r save-pdf}
-# Save the PDF file
+# Save plot as a PDF file
+# Note that `last_plot()` is the default, but we are making it explicit here
ggplot2::ggsave(pca_plot_file, plot = ggplot2::last_plot())
```
@@ -140,8 +140,8 @@ Now we can add a new column with toy batch information and re-store the `colData
```{r add-batch}
# Add batch information
-sample_info$batch <- c("batch1", "batch1", "batch1", "batch1", "batch2",
- "batch1", "batch2", "batch1")
+sample_info$batch <- c("batch1", "batch1", "batch1", "batch1",
+ "batch2", "batch1", "batch2", "batch1")
```
If this batch information were real we would have included it with the sample metadata when we made the original `SummarizedExperiment` object with `tximeta`.
@@ -149,7 +149,7 @@ We would then include it in the model stored in our DESeq2 object using the `des
Here we will take a bit of a shortcut and add it directly to the `colData()` for our `vst()`-transformed data.
```{r coldata-vst, live = TRUE}
-# Add coldata() with batch info to vst_data
+# Add updated colData() with batch info to vst_data
```
diff --git a/RNA-seq/03-gastric_cancer_exploratory.nb.html b/RNA-seq/03-gastric_cancer_exploratory.nb.html
index b4d26555..a953eeb0 100644
--- a/RNA-seq/03-gastric_cancer_exploratory.nb.html
+++ b/RNA-seq/03-gastric_cancer_exploratory.nb.html
@@ -3031,25 +3031,23 @@
Libraries and functions
Directories and files
-
+
# Main data directory
data_dir <- file.path("data", "gastric-cancer")
# directory with the tximeta processed data
txi_dir <- file.path(data_dir, "txi")
-txi_file <- file.path(txi_dir, "gastric-cancer_tximeta.RDS")
# Create a plots directory if it does not exist yet
plots_dir <- file.path("plots", "gastric-cancer")
-if (!dir.exists(plots_dir)) {
- dir.create(plots_dir, recursive = TRUE)
-}
and help us spot any technical effects (more on that below).
-
-
# DESeq2 built in function is called plotPCA and we want to color points by
-# tissue
+
+
# DESeq2 built in function is called plotPCA and we want to
+# color points by tissue
plotPCA(vst_data, intgroup = "tissue")
using ntop=500 top features by variance
-
+
@@ -3168,8 +3165,9 @@
Principal component analysis
ggplot2
-
-
# Save the PDF file
+
+
# Save plot as a PDF file
+# Note that `last_plot()` is the default, but we are making it explicit here
ggplot2::ggsave(pca_plot_file, plot = ggplot2::last_plot())
# Add updated colData() with batch info to vst_data
colData(vst_data) <- sample_info
-
+
# PCA plot - tissue *and* batch
# We want plotPCA to return the data so we can have more control about the plot
-pca_data <- plotPCA(vst_data,
- intgroup = c("tissue", "batch"),
- returnData = TRUE)
data. The quant.sf files for each sample can be found in
data/NB-cell/salmon_quant/<SAMPLE>.
Save the tximeta output as
-data/NB-cell/txi/NB-cell_tximeta.RDS. Note that
+data/NB-cell/txi/NB-cell_tximeta.rds. Note that
data/NB-cell/txi/ is a new directory.
diff --git a/RNA-seq/05-nb_cell_line_DESeq2-live.Rmd b/RNA-seq/05-nb_cell_line_DESeq2-live.Rmd
index 5e9afa62..73b470b6 100644
--- a/RNA-seq/05-nb_cell_line_DESeq2-live.Rmd
+++ b/RNA-seq/05-nb_cell_line_DESeq2-live.Rmd
@@ -50,7 +50,7 @@ library(EnhancedVolcano)
```{r input-files}
# directory with the tximeta processed data
txi_dir <- file.path("data", "NB-cell", "txi")
-txi_file <- file.path(txi_dir, "NB-cell_tximeta.RDS")
+txi_file <- file.path(txi_dir, "NB-cell_tximeta.rds")
```
@@ -61,9 +61,7 @@ We'll create a results directory to hold our results.
```{r results-dir}
# Create a results directory if it doesn't already exist
results_dir <- file.path("results", "NB-cell")
-if (!dir.exists(results_dir)) {
- dir.create(results_dir, recursive = TRUE)
-}
+fs::dir_create(results_dir)
```
We will also need a directory to store our plots.
@@ -77,7 +75,7 @@ We will also need a directory to store our plots.
```{r output-files}
# RDS for the output of DESeq analysis
deseq_file <- file.path(results_dir,
- "NB-cell_DESeq_amplified_v_nonamplified.RDS")
+ "NB-cell_DESeq_amplified_v_nonamplified.rds")
# DESeq2 results table
deseq_df_file <- file.path(results_dir,
@@ -151,8 +149,11 @@ Genes that have very low counts are not likely to yield reliable differential ex
We will keep only genes with total counts of at least 10 across all samples.
```{r filter_ddset}
+# create a vector of TRUE and FALSE values where
+# TRUE corresponds to genes with counts of at least 10
genes_to_keep <- rowSums(counts(ddset)) >= 10
-ddset <- ddset[genes_to_keep, ]
+# use which() to prevent any NAs sneaking through
+ddset <- ddset[which(genes_to_keep), ]
```
@@ -210,11 +211,9 @@ So for this data set we will be using `ashr` ([Stephens 2017](https://doi.org/10
```{r lfc_shrink}
# calculate shrunken log2 fold change estimates
deseq_shrunken <- lfcShrink(deseq_object,
- # the coefficient we want to reestimate
- coef = 2,
- # We will use `ashr` for estimation
- type = "ashr"
- )
+ coef = 2, # the coefficient we want to reestimate
+ type = "ashr" # We will use `ashr` for estimation
+)
```
Let's compare the log2 fold change estimates from the two results tables by creating a plot.
@@ -226,20 +225,23 @@ comparison_df <- data.frame(
lfc_original = deseq_results$log2FoldChange,
lfc_shrunken = deseq_shrunken$log2FoldChange,
logmean = log10(deseq_results$baseMean)
- )
+)
```
Now we can plot the original and shrunken log2 fold change values to see what happened after shrinkage.
```{r plot_comparison}
ggplot(comparison_df,
- aes(x = lfc_original,
- y = lfc_shrunken,
- color = logmean)) +
+ aes(
+ x = lfc_original,
+ y = lfc_shrunken,
+ color = logmean
+ )
+) +
geom_point(alpha = 0.1) +
theme_bw() +
scale_color_viridis_c() +
- coord_cartesian(xlim = c(-10,10), ylim = c(-10,10)) # zoom in on the middle
+ coord_cartesian(xlim = c(-10, 10), ylim = c(-10, 10)) # zoom in on the middle
```
We will now do a bit of manipulation to store the results in a data frame and add the gene symbols.
@@ -279,16 +281,16 @@ Even better, it outputs a `ggplot2` object, so if we want to customize it furthe
```{r volcano}
EnhancedVolcano(deseq_df,
- x = 'log2FoldChange', # fold change statistic to plot
- y = 'pvalue', # significance values
- lab = deseq_df$gene_symbol, # labels for points
- pCutoff = 1e-05, # The p value cutoff we will use (default)
- FCcutoff = 1, # The fold change cutoff (default)
- title = NULL, # no title
- subtitle = NULL, # or subtitle
- caption = NULL, # or caption
- labSize = 3 # smaller labels
- ) +
+ x = "log2FoldChange", # fold change statistic to plot
+ y = "pvalue", # significance values
+ lab = deseq_df$gene_symbol, # labels for points
+ pCutoff = 1e-05, # The p value cutoff we will use (default)
+ FCcutoff = 1, # The fold change cutoff (default)
+ title = NULL, # no title
+ subtitle = NULL, # or subtitle
+ caption = NULL, # or caption
+ labSize = 3 # smaller labels
+) +
# change the overall theme
theme_classic() +
# move the legend to the bottom
diff --git a/RNA-seq/05-nb_cell_line_DESeq2.nb.html b/RNA-seq/05-nb_cell_line_DESeq2.nb.html
index 75271f50..69b713f5 100644
--- a/RNA-seq/05-nb_cell_line_DESeq2.nb.html
+++ b/RNA-seq/05-nb_cell_line_DESeq2.nb.html
@@ -3056,10 +3056,10 @@
Directories and files
Input
-
+
# directory with the tximeta processed data
txi_dir <- file.path("data", "NB-cell", "txi")
-txi_file <- file.path(txi_dir, "NB-cell_tximeta.RDS")
We’ll create a results directory to hold our results.
-
+
# Create a results directory if it doesn't already exist
results_dir <- file.path("results", "NB-cell")
-if (!dir.exists(results_dir)) {
- dir.create(results_dir, recursive = TRUE)
-}
+fs::dir_create(results_dir)
We will also need a directory to store our plots.
-
+
# Create a plots directory if it doesn't already exist
plots_dir <- file.path("plots", "NB-cell")
-if (!dir.exists(plots_dir)) {
- dir.create(plots_dir, recursive = TRUE)
-}
+fs::dir_create(plots_dir)
-
+
# RDS for the output of DESeq analysis
deseq_file <- file.path(results_dir,
- "NB-cell_DESeq_amplified_v_nonamplified.RDS")
+ "NB-cell_DESeq_amplified_v_nonamplified.rds")
# DESeq2 results table
deseq_df_file <- file.path(results_dir,
@@ -3243,11 +3239,11 @@
Preparation
DESeq Dataset creation
-
+
# Create a DESeq2 dataset from `gene_summarized`
# remember that `status` is the variable of interest here
ddset <- DESeqDataSet(gene_summarized,
- design = ~ status)
+ design = ~status)
using counts and average transcript lengths from tximeta
# create a vector of TRUE and FALSE values where
+# TRUE corresponds to genes with counts of at least 10
+genes_to_keep <- rowSums(counts(ddset)) >= 10
+# use which() to prevent any NAs sneaking through
+ddset <- ddset[which(genes_to_keep), ]
@@ -3429,14 +3428,12 @@
Shrinking log2 fold change estimates
2017)
-
+
# calculate shrunken log2 fold change estimates
deseq_shrunken <- lfcShrink(deseq_object,
- # the coefficient we want to reestimate
- coef = 2,
- # We will use `ashr` for estimation
- type = "ashr"
- )
+ coef = 2, # the coefficient we want to reestimate
+ type = "ashr" # We will use `ashr` for estimation
+)
using 'ashr' for LFC shrinkage. If used in published research, please cite:
@@ -3450,12 +3447,12 @@
Shrinking log2 fold change estimates
First we will combine the results into a new data frame.
ggplot(comparison_df,
- aes(x = lfc_original,
- y = lfc_shrunken,
- color = logmean)) +
+ aes(
+ x = lfc_original,
+ y = lfc_shrunken,
+ color = logmean
+ )
+) +
geom_point(alpha = 0.1) +
theme_bw() +
scale_color_viridis_c() +
- coord_cartesian(xlim = c(-10,10), ylim = c(-10,10)) # zoom in on the middle
+ coord_cartesian(xlim = c(-10, 10), ylim = c(-10, 10)) # zoom in on the middle
@@ -3535,18 +3535,18 @@
Making a Volcano Plot
have used before.
-
+
EnhancedVolcano(deseq_df,
- x = 'log2FoldChange', # fold change statistic to plot
- y = 'pvalue', # significance values
- lab = deseq_df$gene_symbol, # labels for points
- pCutoff = 1e-05, # The p value cutoff we will use (default)
- FCcutoff = 1, # The fold change cutoff (default)
- title = NULL, # no title
- subtitle = NULL, # or subtitle
- caption = NULL, # or caption
- labSize = 3 # smaller labels
- ) +
+ x = "log2FoldChange", # fold change statistic to plot
+ y = "pvalue", # significance values
+ lab = deseq_df$gene_symbol, # labels for points
+ pCutoff = 1e-05, # The p value cutoff we will use (default)
+ FCcutoff = 1, # The fold change cutoff (default)
+ title = NULL, # no title
+ subtitle = NULL, # or subtitle
+ caption = NULL, # or caption
+ labSize = 3 # smaller labels
+) +
# change the overall theme
theme_classic() +
# move the legend to the bottom
@@ -3578,8 +3578,8 @@
Session Info
sessionInfo()
-
-
R version 4.4.0 (2024-04-24)
+
+
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS
@@ -3615,16 +3615,16 @@
diff --git a/RNA-seq/06-openpbta_heatmap-live.Rmd b/RNA-seq/06-openpbta_heatmap-live.Rmd
index 7de8399a..8da2be99 100644
--- a/RNA-seq/06-openpbta_heatmap-live.Rmd
+++ b/RNA-seq/06-openpbta_heatmap-live.Rmd
@@ -20,7 +20,7 @@ This notebook will demonstrate how to:
---
In this notebook, we cluster RNA-seq data from the Open Pediatric Brain Tumor Atlas (OpenPBTA) project and create a heatmap.
-OpenPBTA is a collaborative project organized by the CCDL and the Center for Data-Driven Discovery in Biomedicine (D3b) at the Children's Hospital of Philadelphia conducted openly on GitHub.
+OpenPBTA is a collaborative project organized by the Data Lab and the Center for Data-Driven Discovery in Biomedicine (D3b) at the Children's Hospital of Philadelphia conducted openly on GitHub.
You can read more about the project [here](https://github.com/alexslemonade/openpbta-analysis/#openpbta-analysis).
@@ -64,7 +64,7 @@ We have stored the data we'll use in this notebook in `data/open-pbta`.
histologies_file <- file.path(data_dir, "pbta-histologies-subset.tsv")
# The RNA-seq counts table
-rnaseq_file = file.path(data_dir, "pbta-rsem-expected_count-subset.rds")
+rnaseq_file <- file.path(data_dir, "pbta-rsem-expected_count-subset.rds")
```
#### Output files
@@ -302,7 +302,7 @@ Heatmap(zscores_mat,
name = "z-score")
```
-`ComplexHeatmap` gives a few warning here, but nothing to be concerned about.
+`ComplexHeatmap` gives a few warnings here, but nothing to be concerned about.
One is complaining that the color scale may be truncated relative to our data.
The other warns that our very large heatmap itself is being drawn at a lower resolution to speed things up.
diff --git a/RNA-seq/06-openpbta_heatmap.nb.html b/RNA-seq/06-openpbta_heatmap.nb.html
index c27172c4..d6e6cfa8 100644
--- a/RNA-seq/06-openpbta_heatmap.nb.html
+++ b/RNA-seq/06-openpbta_heatmap.nb.html
@@ -2920,26 +2920,135 @@
Set up
Libraries
-
-
# We will manipulate RNASeq data with DESeq2 at the start
+
+
# We will manipulate RNASeq data with DESeq2 at the start
+library(DESeq2)
-
-
Warning message:
-replacing previous import ‘S4Arrays::makeNindexFromArrayViewport’ by ‘DelayedArray::makeNindexFromArrayViewport’ when loading ‘SummarizedExperiment’
-
-
-
library(DESeq2)
-
-# Then we'll be doing a bit of data wrangling with the Tidyverse
+
+
Loading required package: S4Vectors
+
+
+
Loading required package: stats4
+
+
+
Loading required package: BiocGenerics
+
+
+
+Attaching package: 'BiocGenerics'
+
+
+
The following objects are masked from 'package:stats':
+
+ IQR, mad, sd, var, xtabs
Welcome to Bioconductor
+
+ Vignettes contain introductory material; view with
+ 'browseVignettes()'. To cite Bioconductor, see
+ 'citation("Biobase")', and for packages 'citation("pkgname")'.
+
+
+
+Attaching package: 'Biobase'
+
+
+
The following object is masked from 'package:MatrixGenerics':
+
+ rowMedians
+
+
+
The following objects are masked from 'package:matrixStats':
+
+ anyMissing, rowMedians
+
+
+
Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
+'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'
+
+
+
# Then we'll be doing a bit of data wrangling with the Tidyverse
library(tidyverse)
✖ lubridate::second() masks S4Vectors::second()
✖ lubridate::second<-() masks S4Vectors::second<-()
✖ dplyr::slice() masks IRanges::slice()
-ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
-
-
+ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
+
+
# ComplexHeatmap is the package we'll use for making a heatmap
# It will do the hierarchical clustering for us as well
library(ComplexHeatmap)
This message can be suppressed by:
suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
-
+
@@ -2992,7 +3101,7 @@
Directories and files
data/open-pbta.
-
+
data_dir <- file.path("data", "open-pbta")
# We'll store the heatmap in plots/open-pbta - create directory if it doesn't exist yet
@@ -3005,7 +3114,7 @@
Directories and files
Input files
-
+
# The metadata describing the samples
histologies_file <- file.path(data_dir, "pbta-histologies-subset.tsv")
@@ -3019,7 +3128,7 @@
Rows: 607 Columns: 33
+── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
-chr (31): Kids_First_Biospecimen_ID, CNS_region, sample_id, aliquot_id, Kids_First_Participant_ID, experimental_strategy, sample_type, composition, tumor_descriptor, primary_site, reported_gender, race...
+chr (31): Kids_First_Biospecimen_ID, CNS_region, sample_id, aliquot_id, Kids...
dbl (2): OS_days, age_last_update_days
+
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
-
+
Use the chunk below to explore the metadata data frame.
@@ -3062,7 +3173,7 @@
Metadata
the Tidyverse.
-
+
histology_count_df <- histologies_df |>
# Count how many samples are in each short_histology and name the column
# with that number n
@@ -3073,13 +3184,11 @@
Metadata
histology_count_df
-
-
@@ -3088,7 +3197,7 @@
RNA-seq data
Read in the expression count matrix (stored as a data frame).
-
+
# Read in and examine the RNA-seq data
rnaseq_exp <- read_rds(rnaseq_file)
@@ -3104,7 +3213,7 @@
Convert and round
convert from a data frame to a matrix.
-
+
rnaseq_mat <- rnaseq_exp |>
# move gene_id to the rownames
tibble::column_to_rownames("gene_id") |>
@@ -3128,7 +3237,7 @@
# The Okabe Ito palette is recommended for those with color vision deficiencies
histology_colors <- palette.colors(palette = "Okabe-Ito")[2:5]
# `palette.colors()` returns a named vector, which can cause trouble
@@ -3316,7 +3427,7 @@
Heatmap(zscores_mat,
# The distance metric used for clustering the rows
# This is different from the default (Euclidean)
@@ -3372,17 +3483,23 @@
Heatmap itself!
# of the cells of the heatmap itself
name = "z-score")
-
-
The automatically generated colors map from the minus and plus 99^th of the absolute values in the matrix. There are outliers in the matrix whose patterns might be hidden by this
-color mapping. You can manually set the color to `col` argument.
+
+
The automatically generated colors map from the minus and plus 99^th of
+the absolute values in the matrix. There are outliers in the matrix
+whose patterns might be hidden by this color mapping. You can manually
+set the color to `col` argument.
-Use `suppressMessages()` to turn off this message.
-`use_raster` is automatically set to TRUE for a matrix with more than 2000 rows. You can control `use_raster` argument by explicitly setting TRUE/FALSE to it.
+Use `suppressMessages()` to turn off this message.
+
+
+
`use_raster` is automatically set to TRUE for a matrix with more than
+2000 rows. You can control `use_raster` argument by explicitly setting
+TRUE/FALSE to it.
Set `ht_opt$message = FALSE` to turn off this message.
name = "z-score",
use_raster = FALSE) # higher resolution for output (be careful with PDF output!)
-
-
The automatically generated colors map from the minus and plus 99^th of the absolute values in the matrix. There are outliers in the matrix whose patterns might be hidden by this
-color mapping. You can manually set the color to `col` argument.
+
+
The automatically generated colors map from the minus and plus 99^th of
+the absolute values in the matrix. There are outliers in the matrix
+whose patterns might be hidden by this color mapping. You can manually
+set the color to `col` argument.
Use `suppressMessages()` to turn off this message.
-
-
+
+
# Shut down current graphics device
dev.off()
-
-
null device
- 1
+
+
png
+ 2
@@ -3473,17 +3592,17 @@
PCA as an alternative to clustering
plotPCA() function from DESeq2.
-
+
# Use plotPCA, but return the data for custom plotting
pca_df <- plotPCA(vst_data,
ntop = 5000, # use the top 5000 genes by variance
intgroup = "short_histology",
returnData = TRUE)
diff --git a/intro-to-R-tidyverse/01-intro_to_base_R-live.Rmd b/intro-to-R-tidyverse/01-intro_to_base_R-live.Rmd
index 4d94a5fd..be1e6fa7 100644
--- a/intro-to-R-tidyverse/01-intro_to_base_R-live.Rmd
+++ b/intro-to-R-tidyverse/01-intro_to_base_R-live.Rmd
@@ -24,7 +24,7 @@ This notebook will demonstrate how to:
#### *More resources for learning R*
- [Swirl, an interactive tutorial](https://swirlstats.com/)
-- [_R for Data Science_ book](https://r4ds.had.co.nz/)
+- [_R for Data Science_ book](https://r4ds.hadley.nz/)
- [Tutorial on R, RStudio and R Markdown](https://ismayc.github.io/rbasics-book/)
- [Handy R cheatsheets](https://www.posit.co/resources/cheatsheets/)
- [R Markdown website](https://rmarkdown.rstudio.com)
diff --git a/intro-to-R-tidyverse/01-intro_to_base_R.nb.html b/intro-to-R-tidyverse/01-intro_to_base_R.nb.html
index efd95791..ce772647 100644
--- a/intro-to-R-tidyverse/01-intro_to_base_R.nb.html
+++ b/intro-to-R-tidyverse/01-intro_to_base_R.nb.html
@@ -2910,7 +2910,7 @@
diff --git a/intro-to-R-tidyverse/03-intro_to_tidyverse-live.Rmd b/intro-to-R-tidyverse/03-intro_to_tidyverse-live.Rmd
index eeb867ea..a0414c65 100644
--- a/intro-to-R-tidyverse/03-intro_to_tidyverse-live.Rmd
+++ b/intro-to-R-tidyverse/03-intro_to_tidyverse-live.Rmd
@@ -27,7 +27,7 @@ It is a pre-processed [astrocytoma microarray dataset](https://www.refine.bio/ex
**More tidyverse resources:**
-- [R for Data Science](https://r4ds.had.co.nz/)
+- [R for Data Science](https://r4ds.hadley.nz/)
- [tidyverse documentation](https://tidyverse.org/)
- [`dplyr` documentation](https://dplyr.tidyverse.org/)
- [`readr` documentation](https://readr.tidyverse.org/)
@@ -175,7 +175,7 @@ What information is contained in `gene_df`?
One nifty feature that was added to `R` in version 4.1 is the pipe: `|>`.
Pipes are very handy things that allow you to funnel the result of one expression to the next, making your code more streamlined and fluently expressing the flow of data through a series of operations.
-_Note:_ If you are using a version of `R` prior to 4.1 (or looking at older code), pipe functionality was available through the `magrittr` package , which used a pipe that looked like this: `%>%`.
+_Note:_ If you are using a version of `R` prior to 4.1 (or looking at older code), pipe functionality was available through the `magrittr` package, which used a pipe that looked like this: `%>%`.
That pipe was the inspiration for the native R pipe we are using here.
While there are some minor differences, you can mostly treat them interchangeably as long as you load the `magrittr` package or `dplyr`, which also loads that version of the pipe.
diff --git a/intro-to-R-tidyverse/03-intro_to_tidyverse.nb.html b/intro-to-R-tidyverse/03-intro_to_tidyverse.nb.html
index b2459ee8..e04143a4 100644
--- a/intro-to-R-tidyverse/03-intro_to_tidyverse.nb.html
+++ b/intro-to-R-tidyverse/03-intro_to_tidyverse.nb.html
@@ -2903,7 +2903,7 @@
# Filter the mouse data frame to the KEGG pathways that are included in the
# curated gene sets
-mm_kegg_df <- mm_msigdb_df %>%
- dplyr::filter(gs_cat == "C2", # curated gene sets
- gs_subcat == "CP:KEGG") # KEGG pathways
Note: We could specified that we wanted the KEGG gene sets using the category and subcategory arguments of msigdbr(), but we’re going for general steps!
+
Note: We could specified that we wanted the KEGG gene sets using
+the category and subcategory arguments of
+msigdbr(), but we’re going for general steps!
The clusterProfiler function we will use requires a data frame with two columns, where one column contains the term identifier or name and one column contains gene identifiers that match our gene lists we want to check for enrichment. Our data frame with KEGG terms contains Entrez IDs and gene symbols.
+
The clusterProfiler function we will use requires a data
+frame with two columns, where one column contains the term identifier or
+name and one column contains gene identifiers that match our gene lists
+we want to check for enrichment. Our data frame with KEGG terms contains
+Entrez IDs and gene symbols.
+Delimiter: "\t"
+chr (1): Gene
+dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
+
+ℹ Use `spec()` to retrieve the full column specification for this data.
+ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
@@ -657,92 +3249,85 @@
Read in DGE results and prep
Gene identifier conversion
-
Our data frame of DGE results contains Ensembl gene identifiers. So we will need to convert from these identifiers into either the gene symbols or Entrez IDs that are present in the data we extracted with msigdbr().
-
We’re going to convert our identifiers to gene symbols because they are a bit more human readable, but you can, with the change of a single argument, use the same code to convert to many other types of identifiers!
-
The annotation package org.Mm.eg.db contains information for different identifiers. org.Mm.eg.db is specific to Mus musculus – this is what the Mm in the package name is referencing. To perform gene identifier conversion in human (Homo sapiens) we could use org.Hs.eg.db; we would use org.Dr.eg.db for zebrafish (Danio rerio).
-
We can see what types of IDs are available to us in an annotation package with keytypes().
+
Our data frame of DGE results contains Ensembl gene identifiers. So
+we will need to convert from these identifiers into either the gene
+symbols or Entrez IDs that are present in the data we extracted with
+msigdbr().
+
We’re going to convert our identifiers to gene symbols because they
+are a bit more human readable, but you can, with the change of a single
+argument, use the same code to convert to many other types of
+identifiers!
+
The annotation package org.Mm.eg.db contains information
+for different identifiers. org.Mm.eg.db is specific to
+Mus musculus – this is what the Mm in the package
+name is referencing. To perform gene identifier conversion in human
+(Homo sapiens) we could use org.Hs.eg.db; we would
+use org.Dr.eg.db for zebrafish (Danio rerio).
+
We can see what types of IDs are available to us in an annotation
+package with keytypes().
Even though we’ll use this package to convert from Ensembl gene IDs (ENSEMBL) to gene symbols (SYMBOL), we could just as easily use it to convert from an Ensembl transcript ID (ENSEMBLTRANS) to Entrez IDs (ENTREZID).
-
The function we will use to map from Ensembl gene IDs to gene symbols is called mapIds().
+
Even though we’ll use this package to convert from Ensembl gene IDs
+(ENSEMBL) to gene symbols (SYMBOL), we could
+just as easily use it to convert from an Ensembl transcript ID
+(ENSEMBLTRANS) to Entrez IDs (ENTREZID).
+
The function we will use to map from Ensembl gene IDs to gene symbols
+is called mapIds().
-
+
# This returns a named vector which we can convert to a data frame, where
# the keys (Ensembl IDs) are the names
symbols_vector <- mapIds(org.Mm.eg.db, # Specify the annotation package
- # The vector of gene identifiers we want to
+ # The vector of gene identifiers we want to
# map
- keys = vs_low_df$Gene,
+ keys = vs_low_df$Gene,
# What type of gene identifiers we're starting
# with
- keytype = "ENSEMBL",
+ keytype = "ENSEMBL",
# The type of gene identifier we want returned
- column = "SYMBOL",
+ column = "SYMBOL",
# In the case of 1:many mappings, return the
# first one. This is default behavior!
- multiVals = "first")
+ multiVals = "first")
'select()' returned 1:many mapping between keys and columns
@@ -756,13 +3341,19 @@
Gene identifier conversion
-
This message is letting us know that sometimes Ensembl gene identifiers will map to multiple gene symbols. In this case, it’s also possible that a gene symbol will map to multiple Ensembl IDs.
-
Now we are ready to add the gene symbols to our data frame with the DGE results. We can use a join function from the dplyr package to do this, which will use the Ensembl gene IDs in both data frames to determine how to join the rows.
-
Let’s do this first for the comparison to the low stem cell capacity population.
+
This message is letting us know that sometimes Ensembl gene
+identifiers will map to multiple gene symbols. In this case, it’s also
+possible that a gene symbol will map to multiple Ensembl IDs.
+
Now we are ready to add the gene symbols to our data frame with the
+DGE results. We can use a join function from the
+dplyr package to do this, which will use the Ensembl gene
+IDs in both data frames to determine how to join the rows.
+
Let’s do this first for the comparison to the low stem cell capacity
+population.
-
-
vs_low_df <- symbols_df %>%
+
+
vs_low_df <- symbols_df |>
# An *inner* join will only return rows that are in both data frames
dplyr::inner_join(vs_low_df,
# The name of the column that contains the Ensembl gene IDs
@@ -774,37 +3365,43 @@
Gene identifier conversion
Drop NA values
-
Some of these rows have NA values in padj, which can happen for a number of reasons including when all samples have zero counts or a gene has low mean expression.
-
Let’s filter to rows that do not have any NA using a function tidyr::drop_na(). This will also drop genes that have an Ensembl gene identifier but no gene symbol!
+
Some of these rows have NA values in padj,
+which can
+happen for a number of reasons including when all samples have zero
+counts or a gene has low mean expression.
+
Let’s filter to rows that do not have any NA
+using a function tidyr::drop_na(). This will also drop
+genes that have an Ensembl gene identifier but no gene symbol!
-
-
# Remove rows that are not complete (e.g., contain NAs) by filtering to only
+
+
# Remove rows that are not complete (e.g., contain NAs) by filtering to only
# complete rows
-vs_low_df <- vs_low_df %>%
+vs_low_df <- vs_low_df |>
tidyr::drop_na()
-
Now we’ll read in our data frame of DGE results from another comparison. To save us some time during instruction, we’ve already done the gene identifier conversion and filtering to remove NA values in this notebook. We took a different series of steps to achieve the same thing, which is often possible in R!
+
Now we’ll read in our data frame of DGE results from another
+comparison. To save us some time during instruction, we’ve
+already done the gene identifier conversion and filtering to remove
+NA values in this
+notebook. We took a different series of steps to achieve the same
+thing, which is often possible in R!
+Delimiter: "\t"
+chr (2): Gene, gene_symbol
+dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
+
+ℹ Use `spec()` to retrieve the full column specification for this data.
+ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
@@ -812,14 +3409,27 @@
Drop NA values
Over-representation Analysis (ORA)
-
To test for over-representation, we can calculate a p-value with a hypergeometric test (ref).
Where N is the number of genes in the background distribution, M is the number of genes in a pathway, n is the number of genes we are interested in (our marker genes), and k is the number of genes that overlap between the pathway and our marker genes.
Where N is the number of genes in the background
+distribution, M is the number of genes in a pathway,
+n is the number of genes we are interested in (our marker
+genes), and k is the number of genes that overlap between
+the pathway and our marker genes.
Example: Suppose we have 17,980 genes detected in a Microarray study and 57 genes were differentially expressed. Among the differential expressed genes, 28 are annotated to a gene set.
+
Example: Suppose we have 17,980 genes detected in a
+Microarray study and 57 genes were differentially expressed. Among the
+differential expressed genes, 28 are annotated to a gene set.
-
We’ll call genes that are differentially expressed gene_in_interest and genes that are in the gene set in_gene_set.
+
We’ll call genes that are differentially expressed
+gene_in_interest and genes that are in the gene set
+in_gene_set.
@@ -838,7 +3448,10 @@
Over-representation Analysis (ORA)
-
We can assess if the 28 overlapping genes mean that the differentially expressed genes are over-represented in the gene set with the hypergeometric distribution. This corresponds to a one-sided Fisher’s exact test.
+
We can assess if the 28 overlapping genes mean that the
+differentially expressed genes are over-represented in the gene set with
+the hypergeometric distribution. This corresponds to a one-sided
+Fisher’s exact test.
@@ -859,56 +3472,88 @@
Over-representation Analysis (ORA)
-
When we test multiple pathways or gene sets, the p-values then need to be adjusted for multiple hypothesis testing.
+
When we test multiple pathways or gene sets, the
+p-values then need to be adjusted for multiple
+hypothesis testing.
High stem cell capacity ORA
-
Our DGE results are from data published as part of Sachs et al. (2014). The authors sorted populations of primary leukemia cells and examined the stem cell capacity of these cell populations. (This study may sound familiar if you’ve worked on one of our bulk RNA-seq exercise notebooks in the past!)
-
We compared the population that the authors identified as having high stem cell capacity to a low stem cell capacity population. We also compared the high stem cell capacity cells to a mix of populations (e.g., unsorted cells). You can see the code in here.
-
We’re interested in what pathways are over-represented in genes that specifically distinguish the high capacity population from the low capacity population.
-
Let’s generate a list of genes that have higher expression in the high stem cell capacity population compared to the low stem cell capacity population, but we’ll also want to exclude genes that show up in our other comparison to unsorted cells.
-
We’ll start with the high stem cell capacity vs. low stem cell capacity population comparison. Genes with positive log2 fold-changes (LFC) will be more highly expressed in the high stem cell capacity cells based on how we set up the analysis.
+
Our DGE results are from data published as part of Sachs et
+al. (2014). The authors sorted populations of primary leukemia
+cells and examined the stem cell capacity of these cell populations.
+(This study may sound familiar if you’ve worked on one of our bulk
+RNA-seq exercise notebooks in the past!)
+
We compared the population that the authors identified as having high
+stem cell capacity to a low stem cell capacity population. We also
+compared the high stem cell capacity cells to a mix of populations
+(e.g., unsorted cells). You can see the code in here.
+
We’re interested in what pathways are over-represented in genes that
+specifically distinguish the high capacity population from the low
+capacity population.
+
Let’s generate a list of genes that have higher expression in the
+high stem cell capacity population compared to the low stem cell
+capacity population, but we’ll also want to exclude genes that show up
+in our other comparison to unsorted cells.
+
We’ll start with the high stem cell capacity vs. low stem cell
+capacity population comparison. Genes with positive log2 fold-changes
+(LFC) will be more highly expressed in the high stem cell capacity cells
+based on how we set up the analysis.
-
-
vs_low_genes <- vs_low_df %>%
+
+
vs_low_genes <- vs_low_df |>
# Filter to the positive LFC and filter based on significance too (padj)
dplyr::filter(log2FoldChange > 0,
- padj < 0.05) %>%
+ padj < 0.05) |>
# Return a vector of gene symbols
dplyr::pull(gene_symbol)
-
Although we’re picking a commonly used cutoff (FDR < 0.05), it’s still arbitrary and we could just as easily pick a different threshold for our LFC values. When we generate lists of genes of interest for ORA, we typically pick an arbitrary cutoff. This is one of the approach’s weaknesses – we’ve removed all other context.
+
Although we’re picking a commonly used cutoff (FDR <
+0.05), it’s still arbitrary and we could just as easily pick a different
+threshold for our LFC values. When we generate lists of genes of
+interest for ORA, we typically pick an arbitrary cutoff. This is one of
+the approach’s weaknesses – we’ve removed all other context.
Now, we’ll take the same steps for our other results.
We want genes that are in the first comparison but not in the second! We can use setdiff(), a base R function for set operations, to get the list that we want.
+
We want genes that are in the first comparison but not in the second!
+We can use setdiff(), a base R function for set operations,
+to get the list that we want.
-
+
# What genes are in the first set but *not* in the second set
-genes_for_ora <- setdiff(vs_low_genes, vs_unsorted_genes)
As we saw above, calculating the p-value relies on the number of genes in the background distribution. Sometimes folks consider genes from the entire genome to comprise the background, but in the example borrowed from the clusterProfiler authors, they state:
+
As we saw above, calculating the p-value relies on the number of
+genes in the background distribution. Sometimes folks consider genes
+from the entire genome to comprise the background, but in the example
+borrowed from the clusterProfiler authors, they state:
17,980 genes detected in a Microarray study
Where the key phrase is genes detected.
-
If we were unable to include a gene in one of our differential expression comparisons because, for example, it had low mean expression in our experiment and therefore was filtered out in our tidyr::drop_na() step, we shouldn’t include in our background set.
-
We can use another function for set operations, intersect(), to get our background set of genes that were included in both comparisons.
+
If we were unable to include a gene in one of our differential
+expression comparisons because, for example, it had low mean expression
+in our experiment and therefore was filtered out in our
+tidyr::drop_na() step, we shouldn’t include in our
+background set.
+
We can use another function for set operations,
+intersect(), to get our background set of genes that were
+included in both comparisons.
@@ -927,25 +3572,29 @@
Background set
Run enricher()
-
Now that we have our background set, our genes of interest, and our pathway information, we’re ready to run ORA using the enricher() function.
+
Now that we have our background set, our genes of interest, and our
+pathway information, we’re ready to run ORA using the
+enricher() function.
-
+
kegg_ora_results <- enricher(
gene = genes_for_ora, # Genes of interest
- pvalueCutoff = 0.05,
+ pvalueCutoff = 0.05,
pAdjustMethod = "BH", # FDR
universe = background_set, # Background set
- # The pathway information should be a data frame with a term name or
+ # The pathway information should be a data frame with a term name or
# identifier and the gene identifiers
- TERM2GENE = dplyr::select(mm_kegg_df,
+ TERM2GENE = dplyr::select(mm_kegg_df,
gs_name,
gene_symbol)
)
-
Note: using enrichKEGG() is a shortcut for doing ORA using KEGG, but the approach we covered here can be used with any gene sets you’d like!
+
Note: using enrichKEGG() is a shortcut for doing ORA
+using KEGG, but the approach we covered here can be used with any gene
+sets you’d like!
What is returned by enricher()?
@@ -954,7 +3603,9 @@
Run enricher()
-
The information we’re most likely interested in is in the results slot. Let’s convert this into a data frame that we can write to file.
+
The information we’re most likely interested in is in the
+results slot. Let’s convert this into a data frame that we
+can write to file.
@@ -965,43 +3616,47 @@
Run enricher()
Visualizing results
-
We can use a dot plot to visualize our significant enrichment results.
+
We can use a dot plot to visualize our significant enrichment
+results.
enrichplot::dotplot(kegg_ora_results)
-
-
wrong orderBy parameter; set to default `orderBy = "x"`
-
-
+
-
We can use an UpSet plot to visualize the overlap between the gene sets that were returned as significant.
+
We can use an UpSet
+plot to visualize the overlap between the gene sets
+that were returned as significant.
enrichplot::upsetplot(kegg_ora_results)
-
+
-
We can see that some of the DNA repair pathways share genes. Gene sets or pathways aren’t independent, either! Sometimes multiple pathways that show up in our results as significant are indicative of only a handful of genes in our gene list.
-
We can look at the geneID column of our results to see what genes overlap; it’s a good idea to take a look.
+
We can see that some of the DNA repair pathways share genes. Gene
+sets or pathways aren’t independent, either! Sometimes multiple pathways
+that show up in our results as significant are indicative of only a
+handful of genes in our gene list.
+
We can look at the geneID column of our results to see
+what genes overlap; it’s a good idea to take a look.
-
-
kegg_result_df %>%
+
+
kegg_result_df |>
# Use dplyr::select() - the name of the pathway is in the ID column
dplyr::select(ID, geneID)
@@ -1026,67 +3681,76 @@
Session Info
sessionInfo()
-
-
R version 4.0.3 (2020-10-10)
-Platform: x86_64-pc-linux-gnu (64-bit)
-Running under: Ubuntu 20.04 LTS
+
+
We can see that the associated values vary for each row.
@@ -3185,7 +3189,7 @@
Removing duplicates
value.
-
+
filtered_dge_df <- dge_results_df |>
# Sort so that the highest absolute values of the log2 fold change are at the
# top
@@ -3198,20 +3202,18 @@
Removing duplicates
Let’s see what happened to our duplicate identifiers.
-
+
# Subset to & arrange by gene symbols that were duplicated in the original
# data frame of results
filtered_dge_df |>
dplyr::filter(gene_symbol %in% duplicated_gene_symbols) |>
dplyr::arrange(gene_symbol)
-
-
Now we’re ready to prep our pre-ranked list for GSEA.
@@ -3223,7 +3225,7 @@
Pre-ranked list
identifiers. This is step 1 – gene-level statistics.
-
+
lfc_vector <- filtered_dge_df |>
# Extract a vector of `log2FoldChange` named by `gene_symbol`
dplyr::pull(log2FoldChange, name = gene_symbol)
@@ -3234,7 +3236,7 @@
Pre-ranked list
Let’s look at the top ranked values.
-
+
# Look at first entries of the log2 fold change vector
head(lfc_vector)
@@ -3247,7 +3249,7 @@
Pre-ranked list
And the bottom of the list.
-
+
# Look at the last entries of the log2 fold change vector
tail(lfc_vector)
@@ -3267,7 +3269,7 @@
Run GSEA
specific, commonly used gene sets (e.g., gseKEGG()).
-
+
gsea_results <- GSEA(geneList = lfc_vector, # ordered ranked gene list
minGSSize = 25, # minimum gene set size
maxGSSize = 500, # maximum gene set set
@@ -3277,15 +3279,25 @@
Run GSEA
gs_name,
gene_symbol))
-
-
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
-
-preparing geneSet collections...
-GSEA analysis...
-Warning: There are ties in the preranked stats (0.1% of the list).
-The order of those tied genes will be arbitrary, which may produce unexpected results.leading edge analysis...
-done...
-
+
+
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
+
+
+
preparing geneSet collections...
+
+
+
GSEA analysis...
+
+
+
Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (0.1% of the list).
+The order of those tied genes will be arbitrary, which may produce unexpected results.
+
+
+
leading edge analysis...
+
+
+
done...
+
The warning about ties means that there are multiple genes that have
@@ -3309,7 +3321,7 @@
diff --git a/pathway-analysis/03-gene_set_variation_analysis-live.Rmd b/pathway-analysis/03-gene_set_variation_analysis-live.Rmd
index 13369210..84aa07d3 100644
--- a/pathway-analysis/03-gene_set_variation_analysis-live.Rmd
+++ b/pathway-analysis/03-gene_set_variation_analysis-live.Rmd
@@ -1,11 +1,11 @@
---
title: "Pathway analysis: Gene Set Variation Analysis (GSVA)"
-output:
+output:
html_notebook:
toc: true
toc_float: true
author: CCDL for ALSF
-date: 2020
+date: 2024
---
## Objectives
@@ -42,8 +42,6 @@ Note that these scores will depend on the samples included in the dataset when y
### Libraries
```{r libraries}
-# Pipes
-library(magrittr)
# Gene Set Variation Analysis
library(GSVA)
```
@@ -53,12 +51,12 @@ library(GSVA)
#### Directories
```{r directories}
-# We have some medulloblastoma data from the OpenPBTA project that we've
+# We have some medulloblastoma data from the OpenPBTA project that we've
# prepared ahead of time
input_dir <- file.path("data", "open-pbta")
# Create a directory specifically for the results using this dataset
-output_dir <- file.path("results", "open-pbta")
+output_dir <- file.path("results", "open-pbta")
if (!dir.exists(output_dir)) {
dir.create(output_dir, recursive = TRUE)
}
@@ -139,7 +137,7 @@ For GSVA, we need a matrix.
You may notice that GSVA has some commonalities with GSEA.
Rather than ranking genes based on some statistic _we_ selected ahead of time, GSVA fits a model and ranks genes based on their expression level relative to the sample distribution.
This is a way of asking if a gene _i_ is highly or lowly expressed in a sample _j_ in the context of this experiment and ranking accordingly ([Hänzelmann _et al._ 2013](https://doi.org/10.1186/1471-2105-14-7)).
-The pathway-level score calculated is a way of asking how genes _within_ a gene set vary as compared to genes that are _outside_ of that gene set ([Malhotra. 2018](https://towardsdatascience.com/decoding-gene-set-variation-analysis-8193a0cfda3)).
+The pathway-level score calculated is a way of asking how genes _within_ a gene set vary as compared to genes that are _outside_ of that gene set ([Malhotra. 2018](https://towardsdatascience.com/decoding-gene-set-variation-analysis-8193a0cfda3)).
(This is sometimes called a competitive test in gene set enrichment literature.)
The intuition here is that we will get pathway-level scores for each sample that indicate if genes in a pathway vary concordantly in one direction (overexpressed or underexpressed relative to the overall population) ([Hänzelmann _et al._ 2013](https://doi.org/10.1186/1471-2105-14-7)).
@@ -147,21 +145,30 @@ The output is a gene set by sample matrix of GSVA scores.
### Perform GSVA
+The main work for GSVA is done by the `gsva()` function, which can actually do a variety of different enrichment score calculations.
+To specify that we want the algorithm used by [Hänzelmann _et al._ (2013)](https://doi.org/10.1186/1471-2105-14-7), we will pass as the first argument a call to the `gsvaParam()` function, which is where we will put our data, gene sets, and other parameters specific to that algorithm.
+
```{r run_gsva}
-gsva_results <- gsva(rnaseq_mat,
- hallmarks_list,
- method = "gsva",
- # Appropriate for our transformed data on log2-like scale
- kcdf = "Gaussian",
- # Minimum gene set size
- min.sz = 15,
- # Maximum gene set size
- max.sz = 500,
- # Compute Gaussian-distributed scores
- mx.diff = TRUE)
+gsva_results <- gsva(
+ gsvaParam(
+ exprData = rnaseq_mat,
+ geneSets = hallmarks_list,
+ # Minimum gene set size
+ minSize = 15,
+ # Maximum gene set size
+ maxSize = 500,
+ # Kernel for estimation
+ # Gaussian for our transformed data on log2-like scale
+ kcdf = "Gaussian",
+ # enrichment score is the difference between largest positive and negative
+ maxDiff = TRUE
+ ),
+ # Use 4 cores for multiprocessing
+ BPPARAM = BiocParallel::MulticoreParam(4)
+)
```
-**Note: the `gsva()` documentation says we can use `kcdf = "Gaussian"` if we had RNA-seq log-CPMs, log-RPKMs or log-TPMs, but we would use `kcdf = "Poisson"` on integer counts.**
+**Note: the `gsvaParam()` documentation says we can use `kcdf = "Gaussian"` if we had RNA-seq log-CPMs, log-RPKMs or log-TPMs, but we would use `kcdf = "Poisson"` on integer counts.**
```{r gsva_peek}
# Let's explore what the output of gsva() looks like
@@ -169,7 +176,7 @@ gsva_results[1:5, 1:5]
```
### A note on gene set size
-
+
Often the scores of gene set enrichment methods are not comparable between gene sets of different sizes.
Let's do an experiment using randomly generated gene sets to explore this idea a bit more.
@@ -184,7 +191,7 @@ all_genes <- rownames(rnaseq_mat)
set.seed(2020)
```
-Our minimum gene set size earlier was 15 genes and our maximum gene set size was 500 genes.
+Our minimum gene set size earlier was 15 genes and our maximum gene set size was 500 genes.
We'll use the same minimum and maximum values for our random gene sets and some values in between.
```{r gene_set_sizes}
@@ -198,7 +205,7 @@ For each gene set size, we will generate 100 random gene sets.
# Set number of replicates
nreps <- 100
# Generate 100 random gene sets of each size
-random_gene_sets <- rep(gene_set_size, nreps) %>% # Repeat gene sizes so we run `nreps` times
+random_gene_sets <- rep(gene_set_size, nreps) |> # Repeat gene sizes so we run `nreps` times
purrr::map(
# Sample the vector of all genes, choosing the number of items specified
# in the element of gene set size
@@ -210,35 +217,28 @@ random_gene_sets <- rep(gene_set_size, nreps) %>% # Repeat gene sizes so we run
The Hallmarks list we used earlier stored the gene set names as the name of the list, so let's add names to our random gene sets that indicate what size they are and so `gsva()` doesn't get upset.
```{r name_random_gene_sets}
-# We will include the size of the gene set in the gene set name
+# We will include the size of the gene set in the gene set name
# Start by taking the length of each pathway and appending "pathway_" to that
# number
-lengths_vector <- random_gene_sets %>%
+lengths_vector <- random_gene_sets |>
# Get the length of each gene set (number of genes)
- purrr::map(~ length(.x)) %>%
+ purrr::map(~ length(.x)) |>
# Make it "pathway_"
- purrr::map(~ paste0("pathway_", .x)) %>%
+ purrr::map(~ paste0("pathway_", .x)) |>
# Return a vector
purrr::flatten_chr()
# Add the names in lengths_vector to the list - "pathway_"
-random_gene_sets <- random_gene_sets %>%
+random_gene_sets <- random_gene_sets |>
# make.names() appends a "version" if something is not unique
purrr::set_names(nm = make.names(lengths_vector, unique = TRUE))
```
-
+
Run GSVA on our dataset with the same parameters as before, but now with random gene sets.
```{r random_gsva, live = TRUE}
- # Appropriate for our transformed data on
- # log2-like scale
-
- # Minimum gene set size
-
- # Maximum gene set size
-
- # Compute Gaussian-distributed scores
+ # Use 4 cores for multiprocessing
```
@@ -247,20 +247,22 @@ First we need to get this data in an appropriate format for `ggplot2`.
```{r longer_random_gsva}
# The random results are a matrix
-random_long_df <- random_gsva_results %>%
- data.frame() %>%
+random_long_df <- random_gsva_results |>
+ data.frame() |>
# Gene set names are rownames
- tibble::rownames_to_column("gene_set") %>%
+ tibble::rownames_to_column("gene_set") |>
# Get into long format
- tidyr::pivot_longer(cols = -gene_set,
- names_to = "Kids_First_Biospecimen_ID",
- values_to = "gsva_score") %>%
+ tidyr::pivot_longer(
+ cols = -gene_set,
+ names_to = "Kids_First_Biospecimen_ID",
+ values_to = "gsva_score"
+ ) |>
# Remove the .version added by make.names()
- dplyr::mutate(gene_set = stringr::str_remove(gene_set, "\\..*")) %>%
+ dplyr::mutate(gene_set = stringr::str_remove(gene_set, "\\..*")) |>
# Add a column that keeps track of the gene set size
- dplyr::mutate(gene_set_size = stringr::word(gene_set, 2, sep = "_")) %>%
+ dplyr::mutate(gene_set_size = stringr::word(gene_set, 2, sep = "_")) |>
# We want to plot smallest no. genes -> largest no. genes
- dplyr::mutate(gene_set_size = factor(gene_set_size,
+ dplyr::mutate(gene_set_size = factor(gene_set_size,
levels = c(15, 25, 50, 100, 250, 500)))
```
@@ -269,8 +271,8 @@ Let's make a violin plot so we can look at the distribution of scores by gene se
```{r random_violin}
# Violin plot comparing GSVA scores of different random gene set sizes
-random_long_df %>%
- ggplot2::ggplot(ggplot2::aes(x = gene_set_size,
+random_long_df |>
+ ggplot2::ggplot(ggplot2::aes(x = gene_set_size,
y = gsva_score)) +
# Make a violin plot that's a pretty blue!
ggplot2::geom_violin(fill = "#99CCFF", alpha = 0.5) +
@@ -297,7 +299,7 @@ How might you use this information to inform your interpretation of GSVA scores?
### How can you use these scores?
If you did have groups information for your samples, you could test for pathway scores differences between groups.
-Here's [an example](https://htmlpreview.github.io/?https://github.com/AlexsLemonade/OpenPBTA-analysis/blob/9b44bf1c186b3126b16dbe5b87756b3eae3feec2/analyses/gene-set-enrichment-analysis/02-model-gsea.nb.html) of that from the OpenPBTA project itself!
+Here's [an example](https://htmlpreview.github.io/?https://github.com/AlexsLemonade/OpenPBTA-analysis/blob/9b44bf1c186b3126b16dbe5b87756b3eae3feec2/analyses/gene-set-enrichment-analysis/02-model-gsea.nb.html) of that from the OpenPBTA project itself!
You can also visualize this matrix in a heatmap.
Here's a figure from the OpenPBTA project, where the middle panel is a heatmap of GSVA scores that were significantly different between histologies.
@@ -309,9 +311,9 @@ Here's a figure from the OpenPBTA project, where the middle panel is a heatmap o
### Write results to file
```{r write_gsva_results}
-gsva_results %>%
- as.data.frame() %>%
- tibble::rownames_to_column("pathway") %>%
+gsva_results |>
+ as.data.frame() |>
+ tibble::rownames_to_column("pathway") |>
readr::write_tsv(file = gsva_results_file)
```
diff --git a/pathway-analysis/03-gene_set_variation_analysis.nb.html b/pathway-analysis/03-gene_set_variation_analysis.nb.html
index c3733022..eefefc0f 100644
--- a/pathway-analysis/03-gene_set_variation_analysis.nb.html
+++ b/pathway-analysis/03-gene_set_variation_analysis.nb.html
@@ -11,43 +11,2661 @@
-
+
Pathway analysis: Gene Set Variation Analysis (GSVA)
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+code{white-space: pre-wrap;}
+span.smallcaps{font-variant: small-caps;}
+span.underline{text-decoration: underline;}
+div.column{display: inline-block; vertical-align: top; width: 50%;}
+div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
+ul.task-list{list-style: none;}
+
-
