Skip to content

Instantly share code, notes, and snippets.

@tomsing1
Created November 19, 2023 17:32
Show Gist options
  • Save tomsing1/06966801bb6f2c60a10c9891f65b0807 to your computer and use it in GitHub Desktop.
Save tomsing1/06966801bb6f2c60a10c9891f65b0807 to your computer and use it in GitHub Desktop.
Quarto markdown file that uses the quarto-webr extension to run R code in the browser
---
title: "Embedding R into Quarto documents with quarto-webr"
subtitle: "Example: intersecting differential expression results"
author: "Thomas Sandmann"
date: '2023/11/18'
format: html
engine: knitr
webr:
show-startup-message: true
packages: ['ggvenn', 'huxtable']
channel-type: "post-message"
filters:
- webr
editor_options:
chunk_output_type: console
---
```{webr-r}
#| context: setup
kUrlRoot <- paste0(
"https://tomsing1.github.io/blog/posts/quarto-webr"
)
kResults <- c(
"Day1" = "day1_results",
"Day2" = "day2_results",
"Day3" = "day3_results"
)
results <- lapply(kResults, \(result) {
url <- paste(kUrlRoot, paste0(result, ".csv"), sep = "/")
download.file(url, basename(url))
read.csv(basename(url))[, c("gene_name", "adj.P.Val", "logFC")]
})
names(results) <- names(kResults)
# keep only genes present in all results
keep <- Reduce(intersect, lapply(results, \(x) x$gene_name))
results <- lapply(results, \(result) result[result$gene_name %in% keep, ])
# Helper function to apply an FDR cutoff to a list of dataframes
apply_fdr <- function(x, fdr_cutoff = 0.01){
lapply(x, function(df){
df[df$adj.P.Val < fdr_cutoff, ]$gene_name
})
}
```
## Introduction
This interactive website allows users to intersect differential expression
results of three different comparisons from an experiment published by
[Xia et al](https://www.biorxiv.org/content/10.1101/2021.01.19.426731v1).
The experiment was performed in three different batches (e.g. samples were
collected on three different days). For this report, the differences
between samples from homozygous APP-SAA and wiltdype animals were estimated
_separately_ for each batch/day: `Day 1`, `Day 2` and `Day 3`.[^1]
[^1]: Analyzing each batch separately is _not_ the best way to estimate genotype
effects across the full experiment. For this exercise, I just want to obtain
three different contrasts to compare in a Venn diagram. You can find more
details about a more realistic analysis of this experiment in
[this blog post](https://tomsing1.github.io/blog/posts/nextflow-core-quantseq-3-xia/).
By manipulating the code cells include in this website, you can select
a set of genes that passes a user-defined FDR threshold on all three days.[^2]
[^2]: The FDR does not reveal the direction of differential expression. So it
is possible that we will pick up genes that are statistically significantly
differentially expressed in _opposite directions_.
## Step 1: Choose a false-discovery (FDR) threshold
1. Please choose an FDR threshold by changing the `FDR_threshold` the code
box below.
- It will be applied to results from all three comparisons.
- For example, setting `FDR_threshold < 0.1` will retain all genes with a
false-discovery rate < 10% in _all three experiments_.
2. Press the "Run Code" button to create a Venn diagram.
- You can save the plot by right-clicking and selecting `Save image as`.
3. Repeat until you like the results.
```{webr-r}
#| context: interactive
FDR_threshold <- 0.1
filtered <- apply_fdr(results, FDR_threshold)
ggvenn(filtered)
selected <- sort(Reduce(intersect, filtered))
```
## Step 2: List the genes that pass the FDR threshold in all three comparisons
Next, press the `Run Code` button to see the list of genes that pass your
selected FDR threshold.
- If you triple-click on the list of gene symbols, you can copy them into
another document.
```{webr-r}
#| context: interactive
paste(selected, collapse = ", ")
```
### Step 3: See the differential expression results for the selected genes
For more context, the following sections show you tables with the results
from your three experiments for the genes you selected in Step 1.
## Day 1
Hit `Run Code` button to see the expression of the selected genes
or download [results for all genes](https://tomsing1.github.io/blog/posts/quarto-webr/day1_results.csv.gz)
```{webr-r}
#| context: interactive
results[["Day1"]][match(selected, results[["Day1"]]$gene_name), ] |>
huxtable::hux()
```
### Day 2
Hit `Run Code` button to see the expression of the selected genes
or download [results for all genes](https://tomsing1.github.io/blog/posts/quarto-webr/day2_results.csv.gz)
```{webr-r}
#| context: interactive
results[["Day2"]][match(selected, results[["Day2"]]$gene_name), ] |>
huxtable::hux()
```
### Day 3
Hit `Run Code` button to see the expression of the selected genes
or download [results for all genes](https://tomsing1.github.io/blog/posts/quarto-webr/day3_results.csv.gz)
```{webr-r}
#| context: interactive
results[["Day3"]][match(selected, results[["Day3"]]$gene_name), ] |>
huxtable::hux()
```
```{webr-r}
#| context: interactive
sessionInfo()
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment