brain_atlas.R
In cellGeometry: Geometric Single Cell Deconvolution

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE
)

## ----eval=FALSE---------------------------------------------------------------
# library(zellkonverter)
# library(SingleCellExperiment)
# library(cellGeometry)

## ----eval=FALSE---------------------------------------------------------------
# # 33 Gb
# # 2,480,956 cells
# brain <- readH5AD("../c2f66cd5-4ff4-4578-876c-55783a57cf8f.h5ad",
#                   use_hdf5 = TRUE, reader = "R")
# 
# mat <- brain@assays@data$X
# rownames(mat) <- rownames(brain)  # need to add rownames (genes)
# meta <- brain@colData@listData

## ----eval=FALSE---------------------------------------------------------------
# # shows possible cell cluster subclasses and groups
# sort(table(meta$roi))
# sort(table(meta$supercluster_term))
# 
# mk <- cellMarkers(mat, subclass = meta$roi,
#                   cellgroup = meta$supercluster_term,
#                   dual_mean = TRUE, cores = 8)
# # 27 mins (intel 8 cores)

## ----eval=FALSE---------------------------------------------------------------
# library(AnnotationHub)
# ah <- AnnotationHub()
# ensDb_v110 <- ah[["AH113665"]]
# 
# mk <- gene2symbol(mk, ensDb_v110)

## ----eval=FALSE---------------------------------------------------------------
# signature_heatmap(mk, text = FALSE, show_row_names = FALSE,
#                   row_title_rot = 0, column_title_rot = 45)

## ----eval=FALSE---------------------------------------------------------------
# # non-neuronal cells
# # 888,263 cells
# # 4 Gb
# brainNN <- readH5AD("../99f27be8-9fac-451e-9723-9e4c7191589e.h5ad",
#                     use_hdf5 = TRUE, reader = "R")
# 
# mat2 <- brainNN@assays@data$X
# rownames(mat2) <- rownames(brainNN)  # add rownames (genes)
# meta2 <- brainNN@colData@listData
# 
# sort(table(meta2$supercluster_term))
# sort(table(meta2$cell_type))
# 
# mkNN <- cellMarkers(mat2, subclass = meta2$cell_type,
#                     cellgroup = meta2$supercluster_term,
#                     dual_mean = TRUE, cores = 8)
# # 9 mins (intel 8 cores)
# 
# mkNN <- gene2symbol(mkNN, ensDb_v110)
# 
# signature_heatmap(mkNN)

## ----eval=FALSE---------------------------------------------------------------
# mkm <- mergeMarkers(mk, mkNN, transform = "none")
# mkm <- updateMarkers(mkm, expfilter = 0.2)

## ----eval=FALSE---------------------------------------------------------------
# signature_heatmap(mkm, top = 5,
#                   text = FALSE, show_row_names = FALSE,
#                   row_title_rot = 0, column_title_rot = 75,
#                   scale = "sphere")

## ----out.width='100%', echo=FALSE---------------------------------------------
knitr::include_graphics("brain_sig.png")

## ----eval=FALSE---------------------------------------------------------------
# # view signature for just the Amygdata excitatory group
# signature_heatmap(mkm, subset = "Amygdala excitatory")

## ----out.width='50%', fig.align = 'center', echo=FALSE------------------------
knitr::include_graphics("brain_sig_amygdala.png")

## ----eval=FALSE---------------------------------------------------------------
# # generate neuronal cell count samples
# set.seed(3)
# sim_counts <- generate_samples(mk, 30)
# sim_percent <- sim_counts / rowSums(sim_counts) * 100
# 
# # simulate neuronal bulk
# # sample from neuronal count matrix
# # 35 mins (Intel)
# sim_sampled <- simulate_bulk(mat, sim_counts, meta$roi,
#                              times = 1)
# 
# rownames(sim_sampled) <- gene2symbol(rownames(sim_sampled), ensDb_v110)

## ----eval=FALSE---------------------------------------------------------------
# fit <- deconvolute(mk, sim_sampled,
#                    arith_mean = TRUE,
#                    use_filter = FALSE)
# mset <- metric_set(sim_percent, fit$subclass$percent)
# summary(mset)
# 
# pdf("../brain_sim_neuron.pdf",
#     width = 12, height = 12.5)
# plot_set(sim_counts, fit$subclass$output, show_zero = T,
#          mfrow = c(8, 8))
# plot_set(sim_percent, fit$subclass$percent, show_zero = T,
#          mfrow = c(8, 8))
# dev.off()

## ----eval=FALSE---------------------------------------------------------------
# # generate non-neuronal cell count samples
# set.seed(3)
# sim_countsNN <- generate_samples(mkNN, 30)
# sim_percentNN <- sim_countsNN / rowSums(sim_countsNN) * 100
# 
# # simulate non-neuronal bulk
# # sample from non-neuronal count matrix
# # 6.32 mins (Intel)
# sim_sampledNN <- simulate_bulk(mat2, sim_countsNN, meta2$cell_type,
#                                times = 1)
# 
# rownames(sim_sampledNN) <- gene2symbol(rownames(sim_sampledNN), ensDb_v110)

## ----eval=FALSE---------------------------------------------------------------
# # check genenames are the same in both datasets
# identical(rownames(sim_sampled), rownames(sim_sampledNN))
# ## TRUE
# 
# # merge pseudobulk counts
# sim_sampled_merge <- sim_sampled + sim_sampledNN
# 
# # merge sample cell counts (ground truth)
# sim_counts_merge <- cbind(sim_counts, sim_countsNN)
# sim_percent_merge <- sim_counts_merge / rowSums(sim_counts_merge) * 100

## ----eval=FALSE---------------------------------------------------------------
# fitm <- deconvolute(mkm, sim_sampled_merge,
#                     arith_mean = TRUE, use_filter = FALSE)
# # 0.037 secs (ARM)

## ----eval=FALSE---------------------------------------------------------------
# mset <- metric_set(sim_percent_merge, fitm$subclass$percent)
# summary(mset)
# ##  pearson.rsq          Rsq               RMSE
# ## Min.   :0.6634   Min.   :-0.5747   Min.   :0.02479
# ## 1st Qu.:0.9298   1st Qu.: 0.9146   1st Qu.:0.05244
# ## Median :0.9789   Median : 0.9718   Median :0.07411
# ## Mean   :0.9452   Mean   : 0.8969   Mean   :0.10484
# ## 3rd Qu.:0.9917   3rd Qu.: 0.9855   3rd Qu.:0.11979
# ## Max.   :0.9999   Max.   : 0.9958   Max.   :0.51257

## ----eval=FALSE---------------------------------------------------------------
# # scatter plots
# pdf("../brain_sim_merge.pdf",
#     width = 12, height = 12.5)
# plot_set(sim_counts_merge, fitm$subclass$output, show_zero = TRUE,
#          mfrow = c(8, 8))
# plot_set(sim_percent_merge, fitm$subclass$percent, show_zero = TRUE,
#          mfrow = c(8, 8))
# dev.off()

## ----out.width='100%', echo=FALSE---------------------------------------------
knitr::include_graphics("brain_plotset.png")