Onecut3 in mouse hypothalamus development
Evgenii Tretiakov
2021-01-28
Last updated: 2022-10-04
Checks: 7 0
Knit directory: Zupancic_2022/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20220105) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 7f4395e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rproj.user/
Ignored: data/oc3_fin.h5seurat
Ignored: data/rar2020.srt.cont.oc2or3.raw.h5seurat
Unstaged changes:
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Gad1_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Gad2_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Nav1_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Nav2_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Onecut2_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Slc17a6_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Slc32a1_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Th_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Trh_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Onecut3-Trio_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Slc32a1-Onecut3_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Slc32a1-Th_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Slc32a1-Trh_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Th-Onecut3_.pdf
Modified: output/figures/stat-corr-plt_oc3-rna-data-Th-Trh_.pdf
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/dynamic_oc3-pop.Rmd) and
HTML (docs/dynamic_oc3-pop.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 7f4395e | Evgenii O. Tretiakov | 2022-10-04 | workflowr::wflow_publish("analysis/dynamic_oc3-pop.Rmd", update = T, |
| html | e928499 | EugOT | 2022-01-29 | Build site. |
| html | da268b6 | EugOT | 2022-01-29 | rebuild with fixed output device to png and pdf |
| Rmd | 3049d38 | EugOT | 2022-01-29 | fix output device to png and pdf |
| Rmd | b21b61a | EugOT | 2022-01-29 | save final subset of cells |
| html | b21b61a | EugOT | 2022-01-29 | save final subset of cells |
| html | fc91ebf | EugOT | 2022-01-29 | Build site. |
| Rmd | 1df6d48 | EugOT | 2022-01-29 | wflow_publish("*") |
| html | d106e47 | EugOT | 2022-01-29 | Build site. |
| Rmd | 6a77e95 | EugOT | 2022-01-29 | wflow_publish("*") |
| html | 6a77e95 | EugOT | 2022-01-29 | wflow_publish("*") |
# Load tidyverse infrastructure packages
library(here)
# here() starts at /home/etretiakov/src/Zupancic_2022
library(tidyverse)
# ── Attaching packages
# ───────────────────────────────────────
# tidyverse 1.3.2 ──
# ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
# ✔ tibble 3.1.8 ✔ dplyr 1.0.10
# ✔ tidyr 1.2.1 ✔ stringr 1.4.1
# ✔ readr 2.1.2 ✔ forcats 0.5.2
# ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
# ✖ dplyr::filter() masks stats::filter()
# ✖ dplyr::lag() masks stats::lag()
library(magrittr)
#
# Attaching package: 'magrittr'
#
# The following object is masked from 'package:purrr':
#
# set_names
#
# The following object is masked from 'package:tidyr':
#
# extract
library(zeallot)
library(future)
#
# Attaching package: 'future'
#
# The following objects are masked from 'package:zeallot':
#
# %->%, %<-%
# Load packages for scRNA-seq analysis and visualisation
library(sctransform)
library(Seurat)
# Attaching SeuratObject
# Attaching sp
library(SeuratWrappers)
library(SeuratDisk)
# Registered S3 method overwritten by 'SeuratDisk':
# method from
# as.sparse.H5Group Seurat
library(UpSetR)
library(patchwork)
library(Nebulosa)src_dir <- here("code")
data_dir <- here("data")
output_dir <- here("output")
plots_dir <- here(output_dir, "figures")
tables_dir <- here(output_dir, "tables")source(here(src_dir, "functions.R"))
source(here(src_dir, "genes.R"))reseed <- 42
set.seed(seed = reseed)# available cores
n_cores <- available_cores(prop2use = .5)
# Parameters for parallel execution
plan("multicore", workers = n_cores)
# Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
# processing ('multicore') is not supported when running R from RStudio
# because it is considered unstable. For more details, how to control forked
# processing or not, and how to silence this warning in future R sessions, see ?
# parallelly::supportsMulticore
options(future.globals.maxSize = Inf,
future.rng.onMisuse = "ignore")
plan()
# multicore:
# - args: function (..., workers = 16, envir = parent.frame())
# - tweaked: TRUE
# - call: plan("multicore", workers = n_cores)Read data
rar2020_ages_all <- c("E15", "E17", "P00", "P02", "P10", "P23")
rar2020_ages_postnat <- c("P02", "P10", "P23")
samples_df <- read_tsv(here("data/samples.tsv"))
# Rows: 8 Columns: 14
# ── Column specification ────────────────────────────────────────────────────────
# Delimiter: "\t"
# chr (8): sample, age, condition, fullname, name, sex, date, date10x
# dbl (6): ncells, libbatch, seqbatch, perfussed, nt, sn
#
# ℹ 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.
colours_wtree <- setNames(read_lines(here(data_dir, "colours_wtree.tsv")),
1:45)
onecut <- LoadH5Seurat(here(data_dir, "rar2020.srt.cont.oc2or3.raw.h5seurat"))
# Validating h5Seurat file
# Initializing RNA with data
# Adding counts for RNA
# Adding scale.data for RNA
# Adding variable feature information for RNA
# Adding reduction pca
# Adding cell embeddings for pca
# Adding feature loadings for pca
# Adding miscellaneous information for pca
# Adding reduction tsne
# Adding cell embeddings for tsne
# Adding miscellaneous information for tsne
# Adding reduction umap
# Adding cell embeddings for umap
# Adding miscellaneous information for umap
# Adding command information
# Adding cell-level metadata
# Adding miscellaneous information
# Adding tool-specific results
onecut3 <- subset(onecut, subset = Onecut3 > 0)Derive and filter matrix of Onecut3
mtx_oc3 <-
onecut3 %>%
GetAssayData("data", "RNA") %>%
as.data.frame() %>%
t()
rownames(mtx_oc3) <- colnames(onecut3)
# Filter features
filt_low_genes <-
colSums(mtx_oc3) %>%
.[. > quantile(., 0.4)] %>%
names()
mtx_oc3 %<>% .[, filt_low_genes]
min_filt_vector <-
mtx_oc3 %>%
as_tibble() %>%
select(all_of(filt_low_genes)) %>%
summarise(across(.fns = ~ quantile(.x, .1))) %>%
as.list %>%
map(as.double) %>%
simplify %>%
.[colnames(mtx_oc3)]
# Prepare table of intersection sets analysis
content_mtx_oc3 <-
(mtx_oc3 > min_filt_vector) %>%
as_tibble() %>%
mutate_all(as.numeric)Correlation analysis visualisation between different genes
p_corrs <- list(
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Onecut3, y = Trh, xfill = "#ffc400", yfill = "#e22ee2"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Slc32a1, y = Onecut3, xfill = "#0000da", yfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Th, y = Onecut3, xfill = "#006eff", yfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Slc32a1, y = Trh, xfill = "#0000da", yfill = "#e22ee2"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Th, y = Trh, xfill = "#006eff", yfill = "#e22ee2"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
x = Slc32a1, y = Th, xfill = "#0000da", yfill = "#006eff"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Slc32a1, x = Onecut3, yfill = "#0000da", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Th, x = Onecut3, yfill = "#006eff", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Slc17a6, x = Onecut3, yfill = "#ff0000", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Gad1, x = Onecut3, yfill = "#a50202", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Gad2, x = Onecut3, yfill = "#4002a5", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Onecut2, x = Onecut3, yfill = "#6402a5", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Nav1, x = Onecut3, yfill = "#2502a5", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Nav2, x = Onecut3, yfill = "#4002a5", xfill = "#ffc400"),
ggstatsplot::ggscatterstats(as.data.frame(mtx_oc3),
y = Trio, x = Onecut3, yfill = "#2502a5", xfill = "#ffc400")
)
# Registered S3 method overwritten by 'ggside':
# method from
# +.gg ggplot2
n_corrs <- list(
"oc3-rna-data-Onecut3-Trh",
"oc3-rna-data-Slc32a1-Onecut3",
"oc3-rna-data-Th-Onecut3",
"oc3-rna-data-Slc32a1-Trh",
"oc3-rna-data-Th-Trh",
"oc3-rna-data-Slc32a1-Th",
"oc3-rna-data-Onecut3-Slc32a1",
"oc3-rna-data-Onecut3-Th",
"oc3-rna-data-Onecut3-Slc17a6",
"oc3-rna-data-Onecut3-Gad1",
"oc3-rna-data-Onecut3-Gad2",
"oc3-rna-data-Onecut3-Onecut2",
"oc3-rna-data-Onecut3-Nav1",
"oc3-rna-data-Onecut3-Nav2",
"oc3-rna-data-Onecut3-Trio"
)
walk2(n_corrs, p_corrs, save_my_plot, type = "stat-corr-plt")
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.Visualise intersections sets that we are going to use (highlighted)
upset(as.data.frame(content_mtx_oc3),
order.by = "freq",
sets.x.label = "Number of cells",
text.scale = c(2, 1.6, 2, 1.3, 2, 3),
nsets = 15,
sets = c("Gad1", "Gad2", "Slc32a1", "Slc17a6"),
queries = list(
list(
query = intersects,
params = list("Gad1", "Gad2", "Slc32a1"),
active = T
),
list(
query = intersects,
params = list("Slc17a6"),
active = T
)
),
nintersects = 60,
empty.intersections = "on"
)
upset(as.data.frame(content_mtx_oc3),
order.by = "freq",
sets.x.label = "Number of cells",
text.scale = c(2, 1.6, 2, 1.3, 2, 3),
nsets = 15,
sets = c("Th", "Trh", "Slc32a1", "Slc17a6"),
queries = list(
list(
query = intersects,
params = list("Th", "Slc32a1"),
active = T
),
list(
query = intersects,
params = list("Trh", "Slc17a6"),
active = T
)
),
nintersects = 60,
empty.intersections = "on"
)
Regroup factor by stages for more balanced groups
onecut3$age %>% forcats::fct_count()
# # A tibble: 6 × 2
# f n
# <fct> <int>
# 1 E15 156
# 2 E17 104
# 3 P00 77
# 4 P02 74
# 5 P10 174
# 6 P23 7
onecut3$stage <-
onecut3$age %>%
forcats::fct_collapse(`Embrionic day 15` = "E15",
`Embrionic day 17` = "E17",
Neonatal = c("P00", "P02"),
Postnatal = c("P10", "P23"))
onecut3$stage %>% forcats::fct_count()
# # A tibble: 4 × 2
# f n
# <fct> <int>
# 1 Embrionic day 15 156
# 2 Embrionic day 17 104
# 3 Neonatal 151
# 4 Postnatal 181Make subset of stable neurons
onecut3$gaba_status <-
content_mtx_oc3 %>%
select(Gad1, Gad2, Slc32a1) %>%
mutate(gaba = if_all(.fns = ~ .x > 0)) %>%
.$gaba
onecut3$gaba_occurs <-
content_mtx_oc3 %>%
select(Gad1, Gad2, Slc32a1) %>%
mutate(gaba = if_any(.fns = ~ .x > 0)) %>%
.$gaba
onecut3$glut_status <-
content_mtx_oc3 %>%
select(Slc17a6) %>%
mutate(glut = Slc17a6 > 0) %>%
.$glut
oc3_fin <-
subset(onecut3,
cells = union(
WhichCells(onecut3,
expression = gaba_status == TRUE & glut_status == FALSE),
WhichCells(onecut3,
expression = glut_status == TRUE & gaba_occurs == FALSE)))Check contingency tables for neurotransmitter signature
oc3_fin@meta.data %>%
janitor::tabyl(glut_status, gaba_status)
# glut_status FALSE TRUE
# FALSE 0 102
# TRUE 92 0By age
oc3_fin@meta.data %>%
janitor::tabyl(age, gaba_status)
# age FALSE TRUE
# E15 29 9
# E17 12 4
# P00 12 6
# P02 16 2
# P10 21 81
# P23 2 0By stage
oc3_fin@meta.data %>%
janitor::tabyl(stage, gaba_status)
# stage FALSE TRUE
# Embrionic day 15 29 9
# Embrionic day 17 12 4
# Neonatal 28 8
# Postnatal 23 81Make splits of neurons by neurotransmitter signature
oc3_fin$status <- oc3_fin$gaba_status %>%
if_else(true = "GABAergic",
false = "glutamatergic")
Idents(oc3_fin) <- "status"
SaveH5Seurat(
object = oc3_fin,
filename = here(data_dir, "oc3_fin"),
overwrite = TRUE,
verbose = TRUE
)
# Warning: Overwriting previous file /home/etretiakov/src/Zupancic_2022/data/
# oc3_fin.h5seurat
# Creating h5Seurat file for version 3.1.5.9900
# Adding counts for RNA
# Adding data for RNA
# Adding scale.data for RNA
# Adding variable features for RNA
# No feature-level metadata found for RNA
# Adding cell embeddings for pca
# Adding loadings for pca
# No projected loadings for pca
# Adding standard deviations for pca
# No JackStraw data for pca
# Adding cell embeddings for tsne
# No loadings for tsne
# No projected loadings for tsne
# No standard deviations for tsne
# No JackStraw data for tsne
# Adding cell embeddings for umap
# No loadings for umap
# No projected loadings for umap
# No standard deviations for umap
# No JackStraw data for umap
## Split on basis of neurotrans and test for difference
oc3_fin_neurotrans <- SplitObject(oc3_fin, split.by = "status")
## Split on basis of age and test for difference
oc3_fin_ages <- SplitObject(oc3_fin, split.by = "age")DotPlots grouped by age
Expression of GABA receptors in GABAergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$GABAergic,
features = gabar,
group.by = "age",
cols = c("#adffff", "#0084ff"),
col.min = -1, col.max = 1
) + RotatedAxis()
Expression of GABA receptors in glutamatergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$glutamatergic,
features = gabar,
group.by = "age",
cols = c("#ffc2c2", "#ff3c00"),
col.min = -1, col.max = 1
) + RotatedAxis()
Expression of glutamate receptors in GABAergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$GABAergic,
features = glutr,
group.by = "age",
cols = c("#adffff", "#0084ff"),
col.min = -1, col.max = 1
) + RotatedAxis()
Expression of glutamate receptors in glutamatergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$glutamatergic,
features = glutr,
group.by = "age",
cols = c("#ffc2c2", "#ff3c00"),
col.min = -1, col.max = 1
) + RotatedAxis()
DotPlots grouped by stage
Expression of GABA receptors in GABAergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$GABAergic,
features = gabar,
group.by = "stage",
cols = c("#adffff", "#0084ff"),
col.min = -1, col.max = 1
) + RotatedAxis()
# Warning: Scaling data with a low number of groups may produce misleading results
Expression of GABA receptors in glutamatergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$glutamatergic,
features = gabar,
group.by = "stage",
cols = c("#ffc2c2", "#ff3c00"),
col.min = -1, col.max = 1
) + RotatedAxis()
# Warning: Scaling data with a low number of groups may produce misleading results
Expression of glutamate receptors in GABAergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$GABAergic,
features = glutr,
group.by = "stage",
cols = c("#adffff", "#0084ff"),
col.min = -1, col.max = 1
) + RotatedAxis()
# Warning: Scaling data with a low number of groups may produce misleading results
Expression of glutamate receptors in glutamatergic Onecut3 positive cells
DotPlot(object = oc3_fin_neurotrans$glutamatergic,
features = glutr,
group.by = "stage",
cols = c("#ffc2c2", "#ff3c00"),
col.min = -1, col.max = 1
) + RotatedAxis()
# Warning: Scaling data with a low number of groups may produce misleading results
sbs_mtx_oc <-
onecut %>%
GetAssayData("data", "RNA") %>%
as.data.frame() %>%
t()
rownames(sbs_mtx_oc) <- colnames(onecut)
# Filter features
filt_low_genes2 <-
colSums(sbs_mtx_oc) %>%
.[. > quantile(., 0.4)] %>%
names()
sbs_mtx_oc %<>% .[, filt_low_genes2]
min_filt_vector2 <-
sbs_mtx_oc %>%
as_tibble() %>%
select(all_of(filt_low_genes2)) %>%
summarise(across(.fns = ~ quantile(.x, .005))) %>%
as.list %>%
map(as.double) %>%
simplify %>%
.[filt_low_genes2]
# Prepare table of intersection sets analysis
content_sbs_mtx_oc <-
(sbs_mtx_oc > min_filt_vector2) %>%
as_tibble() %>%
mutate_all(as.numeric)
onecut$gaba_status <-
content_sbs_mtx_oc %>%
select(Gad1, Gad2, Slc32a1) %>%
mutate(gaba = if_all(.fns = ~ .x > 0)) %>%
.$gaba
onecut$gaba_occurs <-
content_sbs_mtx_oc %>%
select(Gad1, Gad2, Slc32a1) %>%
mutate(gaba = if_any(.fns = ~ .x > 0)) %>%
.$gaba
onecut$glut_status <-
content_sbs_mtx_oc %>%
select(Slc17a6) %>%
mutate(glut = Slc17a6 > 0) %>%
.$glut
oc_fin <-
subset(onecut,
cells = union(
WhichCells(onecut,
expression = gaba_status == TRUE & glut_status == FALSE),
WhichCells(onecut,
expression = glut_status == TRUE & gaba_occurs == FALSE)))
oc_fin$status <- oc_fin$gaba_status %>%
if_else(true = "GABAergic",
false = "glutamatergic")
Idents(oc_fin) <- "status"
sbs_mtx_oc <-
oc_fin %>%
GetAssayData("data", "RNA") %>%
as.data.frame() %>%
t()
rownames(sbs_mtx_oc) <- colnames(oc_fin)
# Filter features
filt_low_genes2 <-
colSums(sbs_mtx_oc) %>%
.[. > quantile(., 0.4)] %>%
names()
sbs_mtx_oc %<>% .[, filt_low_genes2]
min_filt_vector2 <-
sbs_mtx_oc %>%
as_tibble() %>%
select(all_of(filt_low_genes2)) %>%
summarise(across(.fns = ~ quantile(.x, .005))) %>%
as.list %>%
map(as.double) %>%
simplify %>%
.[filt_low_genes2]
# Prepare table of intersection sets analysis
content_sbs_mtx_oc <-
(sbs_mtx_oc > min_filt_vector2) %>%
as_tibble() %>%
mutate_all(as.numeric)upset(as.data.frame(content_sbs_mtx_oc),
order.by = "freq",
sets.x.label = "Number of cells",
text.scale = c(2, 1.6, 2, 1.3, 2, 3),
nsets = 15,
sets = c("Gad1", "Gad2", "Slc32a1", "Slc17a6", "Onecut3"),
queries = list(
list(
query = intersects,
params = list("Gad1", "Gad2", "Slc32a1", "Onecut3"),
active = T
),
list(
query = intersects,
params = list("Slc17a6", "Onecut3"),
active = T
)
),
nintersects = 20,
empty.intersections = NULL
)
upset(as.data.frame(content_sbs_mtx_oc),
order.by = "freq",
sets.x.label = "Number of cells",
text.scale = c(2, 1.6, 2, 1.3, 2, 3),
nsets = 15,
sets = c("Th", "Trh", "Slc32a1", "Slc17a6", "Onecut3"),
queries = list(
list(
query = intersects,
params = list("Th", "Slc32a1", "Onecut3"),
active = T
),
list(
query = intersects,
params = list("Trh", "Slc17a6", "Onecut3"),
active = T
)
),
nintersects = 20,
empty.intersections = NULL
)
sbs_mtx_oc_full <- content_sbs_mtx_oc |>
select(any_of(c(neurotrans, glutr, gabar, "Trh", "Th", "Onecut3"))) |>
dplyr::bind_cols(oc_fin@meta.data)
sbs_mtx_oc_full |> glimpse()
# Rows: 1,068
# Columns: 74
# $ Slc17a6 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Slc17a8 <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
# $ Slc1a1 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,…
# $ Slc1a2 <dbl> 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0,…
# $ Slc1a6 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
# $ Gad1 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
# $ Slc32a1 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
# $ Slc6a1 <dbl> 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,…
# $ Gria1 <dbl> 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0,…
# $ Gria2 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,…
# $ Gria3 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0,…
# $ Gria4 <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
# $ Grid1 <dbl> 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0,…
# $ Grid2 <dbl> 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,…
# $ Grik1 <dbl> 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1,…
# $ Grik2 <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
# $ Grik3 <dbl> 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1,…
# $ Grik4 <dbl> 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0,…
# $ Grik5 <dbl> 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0,…
# $ Grin1 <dbl> 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0,…
# $ Grin2a <dbl> 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
# $ Grin2b <dbl> 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0,…
# $ Grin2d <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Grin3a <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0,…
# $ Grm1 <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0,…
# $ Grm5 <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1,…
# $ Grm2 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Grm3 <dbl> 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0,…
# $ Grm4 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Grm7 <dbl> 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0,…
# $ Grm8 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0,…
# $ Gabra1 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,…
# $ Gabra2 <dbl> 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0,…
# $ Gabra3 <dbl> 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0,…
# $ Gabra4 <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Gabra5 <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0,…
# $ Gabrb1 <dbl> 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0,…
# $ Gabrb2 <dbl> 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0,…
# $ Gabrb3 <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,…
# $ Gabrg1 <dbl> 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0,…
# $ Gabrg2 <dbl> 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1,…
# $ Gabrg3 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,…
# $ Gabre <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Gabrq <dbl> 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
# $ Gabbr1 <dbl> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1,…
# $ Gabbr2 <dbl> 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0,…
# $ Trh <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
# $ Th <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
# $ Onecut3 <dbl> 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0,…
# $ nGene <int> 3456, 2277, 1660, 2308, 1478, 2216, 955, 2135, 3856, …
# $ nUMI <dbl> 8525, 4649, 2794, 5118, 3123, 4215, 1355, 3898, 9751,…
# $ orig.ident <chr> "FC2P23", "FC2P23", "FC2P23", "FC2P10", "FC2P10", "FC…
# $ res.0.2 <chr> "5", "5", "5", "5", "5", "5", "0", "5", "3", "0", "3"…
# $ res.0.4 <chr> "5", "5", "5", "5", "5", "5", "0", "5", "21", "4", "7…
# $ res.0.8 <chr> "19", "19", "19", "26", "19", "19", "1", "19", "17", …
# $ res.1.2 <chr> "19", "19", "19", "27", "19", "19", "6", "19", "28", …
# $ res.2 <chr> "20", "20", "20", "27", "20", "20", "13", "20", "34",…
# $ tree.ident <int> 9, 9, 9, 9, 9, 9, 9, 9, 11, 11, 11, 11, 11, 11, 12, 1…
# $ pro_Inter <chr> "13", "13", "13", "13", "13", "13", "13", "13", "18",…
# $ pro_Enter <chr> "13", "13", "13", "13", "13", "13", "13", "13", "18",…
# $ tree_final <fct> 11, 11, 11, 11, 11, 11, 11, 11, 8, 8, 8, 8, 8, 8, 6, …
# $ subtree <fct> 10, 10, 10, 10, 10, 10, 10, 10, 12, 12, 12, 12, 12, 1…
# $ prim_walktrap <fct> 1, 1, 1, 1, 1, 1, 3, 1, 3, 3, 3, 15, 3, 3, 3, 12, 11,…
# $ umi_per_gene <dbl> 2.466725, 2.041722, 1.683133, 2.217504, 2.112991, 1.9…
# $ log_umi_per_gene <dbl> 0.3921207, 0.3099965, 0.2261183, 0.3458645, 0.3248976…
# $ nCount_RNA <dbl> 8525, 4649, 2794, 5118, 3123, 4215, 1355, 3898, 9751,…
# $ nFeature_RNA <int> 3456, 2277, 1660, 2308, 1478, 2216, 955, 2135, 3856, …
# $ wtree <fct> 10, 10, 10, 10, 10, 10, 11, 10, 11, 11, 11, 21, 11, 1…
# $ age <fct> P23, P23, P23, P10, P10, P10, E17, E15, P23, P23, P10…
# $ postnat <int> 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,…
# $ gaba_status <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,…
# $ gaba_occurs <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,…
# $ glut_status <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALS…
# $ status <chr> "GABAergic", "GABAergic", "GABAergic", "GABAergic", "…
sbs_mtx_oc_full$modulator <-
sbs_mtx_oc_full %>%
select(Trh, Th) %>%
mutate(modulator = if_any(.fns = ~ .x > 0)) %>%
.$modulator
sbs_mtx_oc_full$oc3 <-
(sbs_mtx_oc_full$Onecut3 > 0)
library(ggstatsplot)
# You can cite this package as:
# Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
# Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
#
# Attaching package: 'ggstatsplot'
# The following object is masked from 'package:future':
#
# %<-%
# for reproducibility
set.seed(123)
## plot
# ggpiestats(
# data = sbs_mtx_oc_full,
# x = modulator,
# y = status,
# perc.k = 1,
# package = "ggsci",
# palette = "category10_d3",
# ) + # further modification with `{ggplot2}` commands
# ggplot2::theme(
# plot.title = ggplot2::element_text(
# color = "black",
# size = 14,
# hjust = 0
# )
# )
# plot
grouped_ggpiestats(
# arguments relevant for `ggpiestats()`
data = sbs_mtx_oc_full,
x = modulator,
y = oc3,
grouping.var = status,
perc.k = 1,
package = "ggsci",
palette = "category10_d3",
# arguments relevant for `combine_plots()`
title.text = "Neuromodulator specification of onecut-driven hypothalamic neuronal lineages by Onecut3 and main neurotransmitter expression",
caption.text = "Asterisks denote results from proportion tests; \n***: p < 0.001, ns: non-significant",
plotgrid.args = list(nrow = 2)
)
sessionInfo()
# R version 4.2.1 (2022-06-23)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 22.04.1 LTS
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
# LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
# [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
# [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
# [9] LC_ADDRESS=C LC_TELEPHONE=C
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] ggstatsplot_0.9.4 Nebulosa_1.6.0 patchwork_1.1.2
# [4] UpSetR_1.4.0 SeuratDisk_0.0.0.9020 SeuratWrappers_0.3.0
# [7] sp_1.5-0 SeuratObject_4.1.2 Seurat_4.2.0
# [10] sctransform_0.3.5 future_1.28.0 zeallot_0.1.0
# [13] magrittr_2.0.3 forcats_0.5.2 stringr_1.4.1
# [16] dplyr_1.0.10 purrr_0.3.4 readr_2.1.2
# [19] tidyr_1.2.1 tibble_3.1.8 ggplot2_3.3.6
# [22] tidyverse_1.3.2 here_1.0.1 workflowr_1.7.0
#
# loaded via a namespace (and not attached):
# [1] estimability_1.4.1 scattermore_0.8
# [3] R.methodsS3_1.8.2 coda_0.19-4
# [5] ragg_1.2.3 bit64_4.0.5
# [7] knitr_1.40 irlba_2.3.5
# [9] multcomp_1.4-20 DelayedArray_0.22.0
# [11] R.utils_2.12.0 data.table_1.14.2
# [13] rpart_4.1.16 RCurl_1.98-1.8
# [15] generics_0.1.3 BiocGenerics_0.42.0
# [17] callr_3.7.2 cowplot_1.1.1
# [19] TH.data_1.1-1 RANN_2.6.1
# [21] correlation_0.8.2 bit_4.0.4
# [23] tzdb_0.3.0 spatstat.data_2.2-0
# [25] xml2_1.3.3 lubridate_1.8.0
# [27] httpuv_1.6.6 SummarizedExperiment_1.26.1
# [29] assertthat_0.2.1 gargle_1.2.1
# [31] xfun_0.33 hms_1.1.2
# [33] jquerylib_0.1.4 evaluate_0.16
# [35] promises_1.2.0.1 fansi_1.0.3
# [37] dbplyr_2.2.1 readxl_1.4.1
# [39] igraph_1.3.5 DBI_1.1.3
# [41] htmlwidgets_1.5.4 spatstat.geom_2.4-0
# [43] googledrive_2.0.0 stats4_4.2.1
# [45] paletteer_1.4.1 ellipsis_0.3.2
# [47] ks_1.13.5 backports_1.4.1
# [49] insight_0.18.4 prismatic_1.1.1
# [51] deldir_1.0-6 MatrixGenerics_1.8.1
# [53] vctrs_0.4.2 SingleCellExperiment_1.18.1
# [55] Biobase_2.56.0 remotes_2.4.2
# [57] ROCR_1.0-11 abind_1.4-5
# [59] cachem_1.0.6 withr_2.5.0
# [61] progressr_0.11.0 rgdal_1.5-32
# [63] vroom_1.5.7 emmeans_1.8.1-1
# [65] mclust_5.4.10 goftest_1.2-3
# [67] cluster_2.1.4 lazyeval_0.2.2
# [69] crayon_1.5.2 hdf5r_1.3.6
# [71] labeling_0.4.2 pkgconfig_2.0.3
# [73] GenomeInfoDb_1.32.4 nlme_3.1-159
# [75] statsExpressions_1.3.3 rlang_1.0.6
# [77] globals_0.16.1 lifecycle_1.0.2
# [79] miniUI_0.1.1.1 MatrixModels_0.5-1
# [81] sandwich_3.0-2 modelr_0.1.9
# [83] rsvd_1.0.5 cellranger_1.1.0
# [85] rprojroot_2.0.3 polyclip_1.10-0
# [87] matrixStats_0.62.0 lmtest_0.9-40
# [89] datawizard_0.6.1 Matrix_1.5-1
# [91] zoo_1.8-11 reprex_2.0.2
# [93] whisker_0.4 ggridges_0.5.4
# [95] processx_3.7.0 googlesheets4_1.0.1
# [97] png_0.1-7 viridisLite_0.4.1
# [99] parameters_0.18.2 bitops_1.0-7
# [101] getPass_0.2-2 R.oo_1.25.0
# [103] KernSmooth_2.23-20 parallelly_1.32.1
# [105] spatstat.random_2.2-0 S4Vectors_0.34.0
# [107] scales_1.2.1 plyr_1.8.7
# [109] ica_1.0-3 zlibbioc_1.42.0
# [111] compiler_4.2.1 RColorBrewer_1.1-3
# [113] fitdistrplus_1.1-8 snakecase_0.11.0
# [115] cli_3.4.1 XVector_0.36.0
# [117] listenv_0.8.0 pbapply_1.5-0
# [119] ps_1.7.1 ggside_0.2.1
# [121] MASS_7.3-58.1 mgcv_1.8-40
# [123] tidyselect_1.1.2 stringi_1.7.8
# [125] textshaping_0.3.6 highr_0.9
# [127] yaml_2.3.5 ggrepel_0.9.1.9999
# [129] grid_4.2.1 sass_0.4.2
# [131] tools_4.2.1 future.apply_1.9.1
# [133] parallel_4.2.1 rstudioapi_0.14
# [135] git2r_0.30.1 janitor_2.1.0.9000
# [137] gridExtra_2.3 farver_2.1.1
# [139] Rtsne_0.16 digest_0.6.29
# [141] BiocManager_1.30.18 rgeos_0.5-9
# [143] shiny_1.7.2 pracma_2.4.2
# [145] Rcpp_1.0.9 GenomicRanges_1.48.0
# [147] broom_1.0.1 BayesFactor_0.9.12-4.4
# [149] performance_0.9.2 later_1.3.0
# [151] RcppAnnoy_0.0.19 httr_1.4.4
# [153] effectsize_0.7.0.5 colorspace_2.0-3
# [155] rvest_1.0.3 fs_1.5.2
# [157] tensor_1.5 reticulate_1.26
# [159] IRanges_2.30.1 splines_4.2.1
# [161] uwot_0.1.14 rematch2_2.1.2
# [163] spatstat.utils_2.3-1 systemfonts_1.0.4
# [165] plotly_4.10.0 xtable_1.8-4
# [167] jsonlite_1.8.1 R6_2.5.1
# [169] pillar_1.8.1 htmltools_0.5.3
# [171] mime_0.12 glue_1.6.2
# [173] fastmap_1.1.0 codetools_0.2-18
# [175] mvtnorm_1.1-3 utf8_1.2.2
# [177] lattice_0.20-45 bslib_0.4.0
# [179] spatstat.sparse_2.1-1 leiden_0.4.3
# [181] survival_3.4-0 rmarkdown_2.16
# [183] munsell_0.5.0 GenomeInfoDbData_1.2.8
# [185] haven_2.5.1 reshape2_1.4.4
# [187] gtable_0.3.1 bayestestR_0.13.0
# [189] spatstat.core_2.4-4