## ----echo=FALSE, include=FALSE------------------------------------------------
library(knitr)
knitr::opts_chunk$set(
cache=TRUE, warning=FALSE,
message=FALSE, cache.lazy=FALSE)
## ----eval=FALSE---------------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly=TRUE))
# install.packages("BiocManager")
#
# BiocManager::install("tidySingleCellExperiment")
## ----message=FALSE------------------------------------------------------------
# Bioconductor single-cell packages
library(scran)
library(scater)
library(igraph)
library(celldex)
library(SingleR)
library(SingleCellSignalR)
# Tidyverse-compatible packages
library(purrr)
library(GGally)
library(tidyHeatmap)
# Both
library(tidySingleCellExperiment)
# Other
library(Matrix)
library(dittoSeq)
## -----------------------------------------------------------------------------
data(pbmc_small, package="tidySingleCellExperiment")
pbmc_small_tidy <- pbmc_small
## -----------------------------------------------------------------------------
pbmc_small_tidy
## -----------------------------------------------------------------------------
counts(pbmc_small_tidy)[1:5, 1:4]
## -----------------------------------------------------------------------------
# Turn off the tibble visualisation
options("restore_SingleCellExperiment_show" = TRUE)
pbmc_small_tidy
## -----------------------------------------------------------------------------
# Turn on the tibble visualisation
options("restore_SingleCellExperiment_show" = FALSE)
## -----------------------------------------------------------------------------
pbmc_small_tidy$file[1:5]
## -----------------------------------------------------------------------------
# Create sample column
pbmc_small_polished <-
pbmc_small_tidy %>%
extract(file, "sample", "../data/([a-z0-9]+)/outs.+", remove=FALSE)
# Reorder to have sample column up front
pbmc_small_polished %>%
select(sample, everything())
## -----------------------------------------------------------------------------
# Use colourblind-friendly colours
friendly_cols <- dittoSeq::dittoColors()
# Set theme
custom_theme <- list(
scale_fill_manual(values=friendly_cols),
scale_color_manual(values=friendly_cols),
theme_bw() + theme(
aspect.ratio=1,
legend.position="bottom",
axis.line=element_line(),
text=element_text(size=12),
panel.border=element_blank(),
strip.background=element_blank(),
panel.grid.major=element_line(linewidth=0.2),
panel.grid.minor=element_line(linewidth=0.1),
axis.title.x=element_text(margin=margin(t=10, r=10, b=10, l=10)),
axis.title.y=element_text(margin=margin(t=10, r=10, b=10, l=10))))
## ----plot1--------------------------------------------------------------------
pbmc_small_polished %>%
ggplot(aes(nFeature_RNA, fill=groups)) +
geom_histogram() +
custom_theme
## ----plot2--------------------------------------------------------------------
pbmc_small_polished %>%
ggplot(aes(groups, nCount_RNA, fill=groups)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter(width=0.1) +
custom_theme
## -----------------------------------------------------------------------------
pbmc_small_polished %>%
join_features(features=c("HLA-DRA", "LYZ")) %>%
ggplot(aes(groups, .abundance_counts + 1, fill=groups)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter(aes(size=nCount_RNA), alpha=0.5, width=0.2) +
scale_y_log10() +
custom_theme
## ----preprocess---------------------------------------------------------------
# Identify variable genes with scran
variable_genes <-
pbmc_small_polished %>%
modelGeneVar() %>%
getTopHVGs(prop=0.1)
# Perform PCA with scater
pbmc_small_pca <-
pbmc_small_polished %>%
runPCA(subset_row=variable_genes)
pbmc_small_pca
## ----pc_plot------------------------------------------------------------------
# Create pairs plot with 'GGally'
pbmc_small_pca %>%
as_tibble() %>%
select(contains("PC"), everything()) %>%
GGally::ggpairs(columns=1:5, aes(colour=groups)) +
custom_theme
## ----cluster------------------------------------------------------------------
pbmc_small_cluster <- pbmc_small_pca
# Assign clusters to the 'colLabels'
# of the 'SingleCellExperiment' object
colLabels(pbmc_small_cluster) <-
pbmc_small_pca %>%
buildSNNGraph(use.dimred="PCA") %>%
igraph::cluster_walktrap() %$%
membership %>%
as.factor()
# Reorder columns
pbmc_small_cluster %>%
select(label, everything())
## ----cluster count------------------------------------------------------------
# Count number of cells for each cluster per group
pbmc_small_cluster %>%
count(groups, label)
## -----------------------------------------------------------------------------
# Identify top 10 markers per cluster
marker_genes <-
pbmc_small_cluster %>%
findMarkers(groups=pbmc_small_cluster$label) %>%
as.list() %>%
map(~ .x %>%
head(10) %>%
rownames()) %>%
unlist()
# Plot heatmap
pbmc_small_cluster %>%
join_features(features=marker_genes) %>%
group_by(label) %>%
heatmap(
.row=.feature, .column=.cell,
.value=.abundance_counts, scale="column")
## ----umap---------------------------------------------------------------------
pbmc_small_UMAP <-
pbmc_small_cluster %>%
runUMAP(ncomponents=3)
## ----umap plot, eval=FALSE----------------------------------------------------
# pbmc_small_UMAP %>%
# plot_ly(
# x=~`UMAP1`,
# y=~`UMAP2`,
# z=~`UMAP3`,
# color=~label,
# colors=friendly_cols[1:4])
## ----eval=FALSE---------------------------------------------------------------
# # Get cell type reference data
# blueprint <- celldex::BlueprintEncodeData()
#
# # Infer cell identities
# cell_type_df <-
# logcounts(pbmc_small_UMAP) %>%
# Matrix::Matrix(sparse = TRUE) %>%
# SingleR::SingleR(
# ref=blueprint,
# labels=blueprint$label.main,
# method="single") %>%
# as.data.frame() %>%
# as_tibble(rownames="cell") %>%
# select(cell, first.labels)
## -----------------------------------------------------------------------------
# Join UMAP and cell type info
data(cell_type_df)
pbmc_small_cell_type <-
pbmc_small_UMAP %>%
left_join(cell_type_df, by="cell")
# Reorder columns
pbmc_small_cell_type %>%
select(cell, first.labels, everything())
## -----------------------------------------------------------------------------
# Count number of cells for each cell type per cluster
pbmc_small_cell_type %>%
count(label, first.labels)
## -----------------------------------------------------------------------------
pbmc_small_cell_type %>%
# Reshape and add classifier column
pivot_longer(
cols=c(label, first.labels),
names_to="classifier", values_to="label") %>%
# UMAP plots for cell type and cluster
ggplot(aes(UMAP1, UMAP2, color=label)) +
facet_wrap(~classifier) +
geom_point() +
custom_theme
## -----------------------------------------------------------------------------
pbmc_small_cell_type %>%
# Add some mitochondrial abundance values
mutate(mitochondrial=rnorm(dplyr::n())) %>%
# Plot correlation
join_features(features=c("CST3", "LYZ"), shape="wide") %>%
ggplot(aes(CST3+1, LYZ+1, color=groups, size=mitochondrial)) +
facet_wrap(~first.labels, scales="free") +
geom_point() +
scale_x_log10() +
scale_y_log10() +
custom_theme
## -----------------------------------------------------------------------------
pbmc_small_nested <-
pbmc_small_cell_type %>%
filter(first.labels != "Erythrocytes") %>%
mutate(cell_class=if_else(
first.labels %in% c("Macrophages", "Monocytes"),
true="myeloid", false="lymphoid")) %>%
nest(data=-cell_class)
pbmc_small_nested
## ----warning=FALSE------------------------------------------------------------
pbmc_small_nested_reanalysed <-
pbmc_small_nested %>%
mutate(data=map(data, ~ {
# feature selection
variable_genes <- .x %>%
modelGeneVar() %>%
getTopHVGs(prop=0.3)
# dimension reduction
.x <- .x %>%
runPCA(subset_row=variable_genes) %>%
runUMAP(ncomponents=3)
# clustering
colLabels(.x) <- .x %>%
buildSNNGraph(use.dimred="PCA") %>%
cluster_walktrap() %$%
membership %>%
as.factor()
return(.x)
}))
pbmc_small_nested_reanalysed
## -----------------------------------------------------------------------------
pbmc_small_nested_reanalysed %>%
# Convert to 'tibble', else 'SingleCellExperiment'
# drops reduced dimensions when unifying data sets.
mutate(data=map(data, ~as_tibble(.x))) %>%
unnest(data) %>%
# Define unique clusters
unite("cluster", c(cell_class, label), remove=FALSE) %>%
# Plotting
ggplot(aes(UMAP1, UMAP2, color=cluster)) +
facet_wrap(~cell_class) +
geom_point() +
custom_theme
## ----eval=FALSE---------------------------------------------------------------
# pbmc_small_nested_interactions <-
# pbmc_small_nested_reanalysed %>%
# # Unnest based on cell category
# unnest(data) %>%
# # Create unambiguous clusters
# mutate(integrated_clusters=first.labels %>% as.factor() %>% as.integer()) %>%
# # Nest based on sample
# nest(data=-sample) %>%
# mutate(interactions=map(data, ~ {
# # Produce variables. Yuck!
# cluster <- colData(.x)$integrated_clusters
# data <- data.frame(assay(.x) %>% as.matrix())
# # Ligand/Receptor analysis using 'SingleCellSignalR'
# data %>%
# cell_signaling(genes=rownames(data), cluster=cluster) %>%
# inter_network(data=data, signal=., genes=rownames(data), cluster=cluster) %$%
# `individual-networks` %>%
# map_dfr(~ bind_rows(as_tibble(.x)))
# }))
#
# pbmc_small_nested_interactions %>%
# select(-data) %>%
# unnest(interactions)
## -----------------------------------------------------------------------------
data(pbmc_small_nested_interactions)
pbmc_small_nested_interactions
## -----------------------------------------------------------------------------
pbmc_small_tidy %>%
aggregate_cells(groups, assays="counts")
## -----------------------------------------------------------------------------
sessionInfo()