Preliminaries
Introduction
Tools: scater, SingleCellExperiment, igraph, DESeq2, etc.
Setup
suppressPackageStartupMessages({
# install.packages("kableExtra")
requireNamespace("kableExtra")
# install.packages("crunch")
requireNamespace("crunch")
# https://rstudio.github.io/DT/
# install.packages("DT")
requireNamespace("DT")
# may need to say "don't restart session"
# install.packages("RCurl")
requireNamespace("RCurl")
# install.packages("dplyr")
library(dplyr)
# install.packages("ggplot2")
library(ggplot2)
#
requireNamespace("tidyr")
# BiocManager::install("org.Mm.eg.db")
requireNamespace("org.Mm.eg.db")
# https://www.bioconductor.org/packages/release/bioc/vignettes/goseq/inst/doc/goseq.pdf
# BiocManager::install("goseq")
requireNamespace("goseq")
# https://vroom.r-lib.org/articles/vroom.html
# install.packages("vroom")
requireNamespace("vroom")
# BiocManager::install("DESeq2")
requireNamespace("DESeq2")
# BiocManager::install("scater")
requireNamespace("scater")
# BiocManager::install("scran")
requireNamespace("scran")
# BiocManager::install("igraph")
requireNamespace("igraph")
# install.packages("pheatmap")
requireNamespace("pheatmap")
# install.packages("stringr")
requireNamespace("stringr")
# install.packages("reshape2")
requireNamespace("reshape2")
# install.packages("pathlibr")
requireNamespace("pathlibr")
# install.packages("pbapply")
requireNamespace("pbapply")
# install.packages("Rtsne")
requireNamespace("Rtsne")
# install.packages("uwot")
requireNamespace("uwot")
# avoid importing `exprs` that leads to clashes
requireNamespace("rlang")
})## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/glue/libs/glue.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/xml2/libs/xml2.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/colorspace/libs/colorspace.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/systemfonts/libs/systemfonts.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/svglite/libs/svglite.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/curl/libs/curl.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/bitops/libs/bitops.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/RCurl/libs/RCurl.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/vctrs/libs/vctrs.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/ellipsis/libs/ellipsis.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/fansi/libs/fansi.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/utf8/libs/utf8.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/tibble/libs/tibble.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/purrr/libs/purrr.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/dplyr/libs/dplyr.so") ...## now dyn.load("/usr/lib/R/library/grid/libs/grid.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/tidyr/libs/tidyr.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/Rcpp/libs/Rcpp.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/bit/libs/bit.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/bit64/libs/bit64.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/fastmap/libs/fastmap.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/cachem/libs/cachem.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/RSQLite/libs/RSQLite.so") ...## now dyn.load("/usr/lib/R/library/parallel/libs/parallel.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/Biobase/libs/Biobase.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/S4Vectors/libs/S4Vectors.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/IRanges/libs/IRanges.so") ...## Warning in find.package(package, fp.lib.loc, quiet = TRUE): package 'mgcv' found more than once, using the first from
## "/home/ra/R/x86_64-pc-linux-gnu-library/4.0/mgcv",
## "/usr/lib/R/library/mgcv"## now dyn.load("/usr/lib/R/library/splines/libs/splines.so") ...## now dyn.load("/usr/lib/R/library/lattice/libs/lattice.so") ...## now dyn.load("/usr/lib/R/library/nlme/libs/nlme.so") ...## now dyn.load("/usr/lib/R/library/Matrix/libs/Matrix.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/mgcv/libs/mgcv.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/XML/libs/XML.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/XVector/libs/XVector.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/GenomicRanges/libs/GenomicRanges.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/Biostrings/libs/Biostrings.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/BiocParallel/libs/BiocParallel.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/Rsamtools/libs/Rsamtools.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/matrixStats/libs/matrixStats.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/DelayedArray/libs/DelayedArray.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/GenomicAlignments/libs/GenomicAlignments.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/rtracklayer/libs/rtracklayer.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/openssl/libs/openssl.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/rappdirs/libs/rappdirs.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/BiasedUrn/libs/BiasedUrn.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/vroom/libs/vroom.so") ...## now dyn.load("/usr/lib/R/library/survival/libs/survival.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/genefilter/libs/genefilter.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/locfit/libs/locfit.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/DESeq2/libs/DESeq2.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/sparseMatrixStats/libs/sparseMatrixStats.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/scuttle/libs/scuttle.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/BiocNeighbors/libs/BiocNeighbors.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/irlba/libs/irlba.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/BiocSingular/libs/BiocSingular.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/igraph/libs/igraph.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/bluster/libs/bluster.so") ...## Note: the specification for class "integer_OR_NULL" in package 'bluster' seems equivalent to one from package 'DelayedArray': not turning on duplicate class definitions for this class.## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/dqrng/libs/dqrng.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/limma/libs/limma.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/edgeR/libs/edgeR.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/statmod/libs/statmod.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/scran/libs/scran.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/plyr/libs/plyr.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/reshape2/libs/reshape2.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/Rtsne/libs/Rtsne.so") ...## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/uwot/libs/uwot.so") ...SessionInfo
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=lzh_TW
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=lzh_TW LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=lzh_TW LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=lzh_TW LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.3 dplyr_1.0.5
##
## loaded via a namespace (and not attached):
## [1] BiocFileCache_1.14.0 systemfonts_1.0.1
## [3] plyr_1.8.6 igraph_1.2.6
## [5] splines_4.0.4 BiocParallel_1.24.1
## [7] GenomeInfoDb_1.26.2 scater_1.18.6
## [9] digest_0.6.27 httpcache_1.2.0
## [11] htmltools_0.5.1.1 viridis_0.5.1
## [13] GO.db_3.12.1 fansi_0.4.2
## [15] magrittr_2.0.1 memoise_2.0.0
## [17] limma_3.46.0 Biostrings_2.58.0
## [19] annotate_1.68.0 matrixStats_0.58.0
## [21] vroom_1.4.0 svglite_2.0.0
## [23] askpass_1.1 prettyunits_1.1.1
## [25] colorspace_2.0-0 blob_1.2.1
## [27] rvest_1.0.0 rappdirs_0.3.3
## [29] BiasedUrn_1.07 xfun_0.21
## [31] crayon_1.4.1 RCurl_1.98-1.2
## [33] jsonlite_1.7.2 org.Mm.eg.db_3.12.0
## [35] genefilter_1.72.1 survival_3.2-10
## [37] glue_1.4.2 kableExtra_1.3.4
## [39] gtable_0.3.0 zlibbioc_1.36.0
## [41] XVector_0.30.0 webshot_0.5.2
## [43] DelayedArray_0.16.2 BiocSingular_1.6.0
## [45] SingleCellExperiment_1.12.0 BiocGenerics_0.36.0
## [47] scales_1.1.1 pheatmap_1.0.12
## [49] DBI_1.1.1 edgeR_3.32.1
## [51] Rcpp_1.0.6 viridisLite_0.3.0
## [53] xtable_1.8-4 progress_1.2.2
## [55] dqrng_0.2.1 bit_4.0.4
## [57] rsvd_1.0.3 stats4_4.0.4
## [59] DT_0.17 htmlwidgets_1.5.3
## [61] httr_1.4.2 RColorBrewer_1.1-2
## [63] ellipsis_0.3.1 pkgconfig_2.0.3
## [65] XML_3.99-0.5 scuttle_1.0.4
## [67] uwot_0.1.10 sass_0.3.1
## [69] dbplyr_2.1.0 locfit_1.5-9.4
## [71] utf8_1.1.4 reshape2_1.4.4
## [73] tidyselect_1.1.0 rlang_0.4.10
## [75] AnnotationDbi_1.52.0 munsell_0.5.0
## [77] tools_4.0.4 cachem_1.0.4
## [79] generics_0.1.0 RSQLite_2.2.3
## [81] evaluate_0.14 stringr_1.4.0
## [83] fastmap_1.1.0 yaml_2.2.1
## [85] knitr_1.31 bit64_4.0.5
## [87] geneLenDataBase_1.26.0 purrr_0.3.4
## [89] pbapply_1.4-3 nlme_3.1-152
## [91] sparseMatrixStats_1.2.1 scran_1.18.5
## [93] xml2_1.3.2 biomaRt_2.46.3
## [95] compiler_4.0.4 rstudioapi_0.13
## [97] beeswarm_0.3.1 curl_4.3
## [99] tibble_3.1.0 statmod_1.4.35
## [101] geneplotter_1.68.0 bslib_0.2.4
## [103] stringi_1.5.3 GenomicFeatures_1.42.1
## [105] lattice_0.20-41 bluster_1.0.0
## [107] Matrix_1.3-2 vctrs_0.3.6
## [109] pillar_1.5.1 lifecycle_1.0.0
## [111] goseq_1.42.0 jquerylib_0.1.3
## [113] BiocNeighbors_1.8.2 bitops_1.0-6
## [115] irlba_2.3.3 rtracklayer_1.50.0
## [117] GenomicRanges_1.42.0 R6_2.5.0
## [119] gridExtra_2.3 vipor_0.4.5
## [121] IRanges_2.24.1 assertthat_0.2.1
## [123] SummarizedExperiment_1.20.0 openssl_1.4.3
## [125] DESeq2_1.30.1 withr_2.4.1
## [127] GenomicAlignments_1.26.0 Rsamtools_2.6.0
## [129] S4Vectors_0.28.1 GenomeInfoDbData_1.2.4
## [131] crunch_1.27.7 mgcv_1.8-34
## [133] parallel_4.0.4 hms_1.0.0
## [135] pathlibr_0.1.0 grid_4.0.4
## [137] beachmat_2.6.4 tidyr_1.1.3
## [139] rmarkdown_2.7 DelayedMatrixStats_1.12.3
## [141] MatrixGenerics_1.2.1 Rtsne_0.15
## [143] Biobase_2.50.0 ggbeeswarm_0.6.0Paths
BASEPATH <- pathlibr::Path$new(".")
stopifnot("main.Rmd" %in% names(BASEPATH$dir))path_to <- (function(.) Sys.glob(BASEPATH$join("../..")$join(.)$show))Plotting
ggplot2::theme_set(theme_light(base_size = 15))kable <- function(.) {
kableExtra::kbl(., align = "c") %>%
kableExtra::kable_paper("hover", full_width = F)
}options(DT.options = list(
pageLength = 16,
language = list(search = "Filter:")
))link <-
function(x) {
dplyr::case_when(
# If x is like GO:003453
startsWith(x, "GO:") ~ paste0(
"<a target='_blank'",
"href='", "http://amigo.geneontology.org/amigo/term/", x, "'",
">", x, "</a>"
),
# Otherwise try this
TRUE ~ paste0(
"<a target='_blank'",
"href='", "http://www.informatics.jax.org/marker/summary?nomen=", x, "'",
">", x, "</a>"
)
)
}Dataframes
# Use first column as index
by_col1 <- (function(.) tibble::column_to_rownames(., colnames(.)[1]))
# Use index as new column `name`
ind2col <- (function(., name) tibble::rownames_to_column(., var = name))arrange.like <- (function(X, T) X[rownames(T), , drop = FALSE])
rows.of <- (function(rows, X) X[rows, , drop = FALSE])Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 4)
from_csv <- function(.) {
vroom::vroom(., del = "\t", progress = FALSE) %>%
by_col1() %>%
suppressMessages()
}Misc
# Wrap functions that do `library(...)`
# and force-detach all attached packages
unattach_packages <- function(f) {
# install.packages("crunch")
requireNamespace("crunch")
# wrapper for f
function(...) {
base::with(
crunch::ContextManager(
enter = base::search,
exit = function() {
for (pkg in setdiff(base::search(), enter)) {
base::detach(pkg, character.only = TRUE, force = TRUE)
}
},
# this overwrites the variable `enter`, if existent
as = "enter"
),
{
f(...)
}
)
}
}
# Usage:
# unattach_packages(function(x, y) { library(goseq); x + y })(1, 2)Datasets
Notes
–
Transcript info
tx <- list(
# Mouse
mm = merge(
(
# https://doi.org/doi:10.18129/B9.bioc.org.Mm.eg.db
org.Mm.eg.db::org.Mm.egSYMBOL2EG
),
(
TxDb.Mmusculus.UCSC.mm39.refGene::TxDb.Mmusculus.UCSC.mm39.refGene %>%
GenomicFeatures::transcriptLengths() %>%
dplyr::group_by(gene_id) %>%
dplyr::summarize(tx_len = median(tx_len))
),
by = "gene_id"
)
# data.frame with transcript lengths
# (gene_id is unique, symbol is not)
# gene_id symbol avg_tx_len
# 1 100009600 Zglp1 1010
# 2 100009609 Vmn2r65 2538
# 3 100009614 Gm10024 564
# etc.
) %>% lapply(
# Reindex by symbol
function(.) {
(.) %>%
dplyr::group_by(symbol) %>%
dplyr::summarize(tx_len = median(tx_len)) %>%
by_col1()
}
)Gene sets
# C7: immunologic signature gene sets
# http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp
# ported to Mus musculus:
# http://bioinf.wehi.edu.au/MSigDB/index.htmlAB – Mouse ctx & hip (10x)
Load
The whole mouse wc&h dataset has 1M cells x 50k genes. We have restricted it to certain cell-types.
sce_abm_ch <-
SingleCellExperiment::SingleCellExperiment(
assays = list(
counts = (
path_to("*/AllenBrain/Mouse-WCH*/*fewer_cells/data.*.gz") %>%
from_csv() %>%
`[`(, sort(colnames(.))) %>%
# conversion to matrix takes longest
# but is necessary for some operations below
as.matrix()
)
),
colData = (
path_to("*/AllenBrain/Mouse-WCH*/f*cells/meta.*.gz") %>%
from_csv() %>%
`[`(sort(rownames(.)), )
)
) %>%
# Order by cell type: sometimes convenient below
`[`(, order((.)$cell_type_alias_label))To reduce the memory footprint, only keep recognized genes.
sce_abm_ch <-
sce_abm_ch %>%
`[`(sort(intersect(rownames(.), rownames(tx$mm))), )Attach a “Is VLMC” flag to samples.
sce_abm_ch$subclass_label <-
sce_abm_ch$subclass_label %>%
as.factor() %>%
relevel(ref = "VLMC")
sce_abm_ch$is_vlmc <-
sce_abm_ch$subclass_label %>%
magrittr::equals("VLMC")Cell types
exclude_celltypes <- c("Astro", "Endo")sce_abm_ch$subclass_label %>%
data.frame(x = .) %>%
ggplot() +
aes(x = x, fill = if_else(x %in% exclude_celltypes, "drop", "keep")) +
geom_bar() +
scale_y_log10() +
# hist_theme +
theme(
axis.title.x = element_blank(),
axis.text.x = element_text(a = 0, hjust = 0.5),
legend.title = element_blank(),
)## now dyn.load("/home/ra/R/x86_64-pc-linux-gnu-library/4.0/farver/libs/farver.so") ...
Diagnostics
Attach diagnostics by scater.
sce_abm_ch <-
sce_abm_ch %>%
{
scater::addPerCellQC(
x = (.),
subsets = list(Mito = grep("mt-", rownames(.)))
)
} %>%
# Used for the heatmap
scater::logNormCounts() %>%
# Not all are used further below:
scater::runPCA(name = "PCA", ncomponents = 20) %>%
# scater::runTSNE(perplexity = 10, dimred = "PCA") %>%
scater::runUMAP()Samples with fewer than 1000 counts have been discarded.
sce_abm_ch %>%
`[`(, order((.)$cell_type_alias_label)) %>%
scater::plotColData(
x = "sum",
y = "detected",
colour_by = "subclass_label"
) +
labs(x = "Sum of counts", y = "Genes detected") +
scale_color_brewer(palette = "Set1", name = "Cell type")
Each bar corresponds to a cell:
sce_abm_ch %>%
# Subsample for expediency
`[`(, union((1:ncol(.))[(.)$is_vlmc], sample(1:ncol(.), 777))) %>%
# Plot in the right order for (some) visibility
`[`(, order((.)$cell_type_alias_label)) %>%
scater::plotHighestExprs(
exprs_values = "counts",
colour_cells_by = "subclass_label",
as_percentage = TRUE
) +
scale_color_discrete(name = "Cell type") +
labs(x = "Counts, %", y = "Gene")
Analysis – Mouse
Overview
The cells cluster by cell type.
sce_abm_ch %>%
# Plot in the right order
`[`(, order((.)$cell_type_alias_label)) %>%
{
function(which_dimred) {
scater::plotReducedDim(
dimred = which_dimred,
object = (.),
colour_by = "subclass_label"
# colour_by = "cell_type_alias_label"
) +
# theme(legend.position = "none") +
scale_color_brewer(palette = "Set1", name = "Cell type")
# The colors that come with the dataset are all grey:
# scale_color_manual(values = c("#3B3F32", "#697255", "#535944"))
}
} %>%
{
gridExtra::grid.arrange((.)("PCA"), (.)("UMAP"), nrow = 2)
}
VLMC vs Other
Before focusing on the VLMC cells, we check for differential genes of VLMC vs the rest. We take a subsample of cells for expedience.
# This does not modify the original:
# sce_abm_ch %>% list(sce = .) %>% within({sce$dummy <- 2})
de_results.ctrl_vs_vlmc <-
sce_abm_ch %>%
list(sce = .) %>%
within({
# SUBSAMPLE:
# Take VLMC and sample some from the rest
jj <- 1:ncol(sce)
set.seed(15)
jj_vlmc <- jj[sce$is_vlmc]
jj_ctrl <- jj[!sce$is_vlmc] %>% sample(3 * length(jj_vlmc))
sce <- sce[, union(jj_vlmc, jj_ctrl)]
# Intermediate report (facultative)
sce$subclass_label %>%
table() %>%
print()
}) %>%
with({
sce %>%
# Package into DESeq object; DE between VLMC and Other
DESeq2::DESeqDataSet(design = ~is_vlmc) %>%
DESeq2::estimateSizeFactors() %>%
DESeq2::DESeq(quiet = TRUE, fitType = "local") %>%
DESeq2::results()
}) %>%
as.data.frame() %>%
slice_min(pvalue, n = 500)## .
## VLMC Astro Endo
## 129 317 70de_results.ctrl_vs_vlmc %>%
select(mean = baseMean, l2fc = log2FoldChange, p = pvalue, padj) %>%
ind2col("symbol") %>%
mutate(symbol = link(symbol)) %>%
DT::datatable(escape = FALSE) %>%
DT::formatSignif(c("mean", "l2fc", "p", "padj"), digits = 3)Establish a color scheme for the following.
annotation_colors <-
sce_abm_ch %>% {
# This defines annotation colors
list(
# Create colors for the given cell types
subclass_label = (
(.)$subclass_label %>%
unique() %>%
sort(decreasing = FALSE) %>%
list(
x = .,
c = RColorBrewer::brewer.pal(n = length(.), name = "Set1")
)
),
# Create colors for the given cell sub types
cell_type_alias_label = (
(.)$cell_type_alias_label %>%
unique() %>%
list(
x = .,
c = RColorBrewer::brewer.pal(n = length(.), name = "Spectral")
)
),
# Use the given colors for the brain anatomical region
region_label = data.frame(
x = (.)$region_label,
c = (.)$region_color
) %>% unique()
) %>%
lapply(function(.) {
(.) %>%
as.data.frame() %>%
pull(c, x) %>%
unlist()
})
}Heatmap of those DE genes for all cells. Clustering of samples using topmost DE genes does collect the VLMC cells together:
de_results.ctrl_vs_vlmc %>%
slice_min(pvalue, n = 55, with_ties = TRUE) %>%
rownames() %>%
rows.of(sce_abm_ch) %>%
scater::logNormCounts() %>%
# Order columns by cell type
`[`(, order((.)$cell_type_alias_label)) %>%
{
pheatmap::pheatmap(
mat = (
(.)@assays@data$logcounts
),
annotation_col = (
(.)@colData %>%
as.data.frame() %>%
dplyr::select(subclass_label, cell_type_alias_label)
),
cluster_cols = TRUE, # (!)
show_colnames = FALSE,
fontsize_row = 6,
treeheight_row = 0,
treeheight_col = 0,
color = grDevices::colorRampPalette(c("white", "blue"))(10),
annotation_colors = annotation_colors
)
}
Let’s look at the genes that are identified as markers in the “betsholtzlab dataset”. Some astrocytes invade the VLMC cluster.
sce_abm_ch %>%
`[`(
intersect(
rownames(.),
rownames(from_csv(path_to("*/Betsholtz-2018/*/mrkr.*")))
),
) %>%
{
pheatmap::pheatmap(
mat = (
(.)@assays@data$logcounts
),
annotation_col = (
(.)@colData %>%
as.data.frame() %>%
dplyr::select(subclass_label, cell_type_alias_label)
),
cluster_cols = TRUE, # (!)
show_colnames = FALSE,
fontsize_row = 6,
treeheight_row = 0,
treeheight_col = 0,
color = grDevices::colorRampPalette(c("white", "blue"))(10),
annotation_colors = annotation_colors
)
}
Filter for VLMC
sce_abm_ch.vlmc <-
sce_abm_ch %>%
`[`(, (.)$is_vlmc)Cluster within VLMC
UMAP roughly reflects the cell subtypes provided with the data.
sce_abm_ch.vlmc %>%
scater::plotReducedDim(dimred = "UMAP", colour_by = "cell_type_alias_label")
We wish to capture exactly two clusters for DE analysis. The subsequent analysis assumes this.
sce_abm_ch.vlmc$my_cluster_label <-
sce_abm_ch.vlmc %>%
list(sce = ., labels = (.)$subclass_label) %>%
with({
# Cluster membership (integers)
kk <- sce %>%
# k is hand-picked to yield two clusters
scran::buildSNNGraph(k = 25, use.dimred = "UMAP", type = "rank") %>% {
set.seed(43)
(.) %>%
# connected components
# igraph::components() %>%
# propagating labels by iterative majority voting
igraph::cluster_label_prop(.) %>%
igraph::membership()
}
# Letters instead of integers
LETTERS[kk] %>%
# make labels like A_VLMC, B_VLMC
paste(labels, sep = "_") %>%
as.factor() %>%
relevel(ref = "A_VLMC")
})
# Expect two clusters; otherwise visualize
sce_abm_ch.vlmc %>% {
print(levels((.)$my_cluster_label))
if (!length(levels((.)$my_cluster_label))) {
(.) %>%
scater::plotReducedDim(dimred = "UMAP", colour_by = "my_cluster_label")
}
}## [1] "A_VLMC" "B_VLMC"Henceforth, use consistent coloring.
my_cluster_colors <-
sce_abm_ch.vlmc %>%
`$`("my_cluster_label") %>%
levels() %>%
list(x = ., c = c("DarkRed", "DarkOrange")) %>%
as.data.frame()
annotation_colors$my_cluster_label <-
my_cluster_colors %>%
pull(c, x)See the proposed clustering.
sce_abm_ch.vlmc %>%
scater::plotReducedDim(
dimred = "UMAP",
colour_by = "my_cluster_label"
) +
scale_color_manual(
values = my_cluster_colors %>% pull(c, x),
name = "Cluster ID"
)
stopifnot(assertthat::are_equal(levels(sce_abm_ch.vlmc$my_cluster_label), c("A_VLMC", "B_VLMC")))DE genes
Identify DE genes between the two clusters within VLMC.
de_results.within_vlmc <-
sce_abm_ch.vlmc %>%
DESeq2::DESeqDataSet(design = ~my_cluster_label) %>%
DESeq2::estimateSizeFactors() %>%
DESeq2::DESeq(quiet = TRUE, fitType = "local", parallel = TRUE) %>%
DESeq2::results() %>%
as.data.frame() %>%
slice_min(pvalue, n = 500)## converting counts to integer modeSee as table.
de_results.within_vlmc %>%
select(mean = baseMean, l2fc = log2FoldChange, p = pvalue, padj = padj) %>%
ind2col("symbol") %>%
mutate(symbol = link(symbol)) %>%
DT::datatable(escape = FALSE) %>%
DT::formatSignif(c("mean", "l2fc", "p", "padj"), digits = 3)Heatmap of a few top DE genes, annotated by the inferred cluster and the available cell type info.
de_results.within_vlmc %>%
slice_min(pvalue, n = 77) %>%
rownames() %>%
rows.of(sce_abm_ch.vlmc) %>%
`[`(, order((.)$cell_type_alias_label)) %>%
{
pheatmap::pheatmap(
mat = (
(.)@assays@data$logcounts
),
annotation_col = (
(.)@colData %>%
as.data.frame() %>%
# of possible interest also: region_label
dplyr::select(region_label, cell_type_alias_label, my_cluster_label)
),
cluster_cols = FALSE,
cluster_rows = FALSE,
show_colnames = FALSE,
fontsize_row = 6,
treeheight_row = 0,
treeheight_col = 0,
color = grDevices::colorRampPalette(c("white", "blue"))(10),
annotation_colors = annotation_colors
)
}
In context: show all cells along the within-VLMC-DE genes.
sce_abm_ch$my_cluster_label <-
sce_abm_ch.vlmc@colData %>%
as.data.frame() %>%
select(my_cluster_label) %>%
`[`(colnames(sce_abm_ch), )
de_results.within_vlmc %>%
slice_min(pvalue, n = 77) %>%
rownames() %>%
rows.of(sce_abm_ch) %>%
`[`(, order((.)$cell_type_alias_label)) %>%
{
pheatmap::pheatmap(
mat = (
(.)@assays@data$logcounts
),
annotation_col = (
(.)@colData %>%
as.data.frame() %>%
# of possible interest also: region_label
dplyr::select(region_label, cell_type_alias_label, my_cluster_label)
),
cluster_cols = FALSE,
cluster_rows = FALSE,
show_colnames = FALSE,
fontsize_row = 6,
treeheight_row = 0,
treeheight_col = 0,
color = grDevices::colorRampPalette(c("white", "blue"))(10),
annotation_colors = annotation_colors
)
}
Pathway / geneset enrichment: preliminaries
Formatting of the GO-terms DE table:
show_go <- function(tbl) {
tbl %>%
dplyr::select("category", "overrep.p", "#DE", "#cat", "term", "GO") %>%
DT::datatable(escape = FALSE) %>%
DT::formatSignif("overrep.p", digits = 3)
}We take some of the top DE (up or down) genes, which requires to set a cut-off. This impacts on the size of GO categories that will be found below.
# Number of top DE genes to keep
# for pathway enrichment analysis
n <- 77
# Note: `lfc` contains NA where no DE
lfc <-
de_results.within_vlmc %>%
# Select top DE genes
slice_min(padj, n = n, with_ties = TRUE) %>%
# Focus on this column
select(lfc = log2FoldChange) %>%
# Rownames as `symbol` column
ind2col("symbol") %>%
# Attach all the other symbols
add_row(symbol = rownames(sce_abm_ch)) %>%
# Keeps only the first one
distinct(symbol, .keep_all = TRUE) %>%
# As named vector (for goseq)
pull(lfc, symbol)
stopifnot(assertthat::are_equal(length(lfc), nrow(sce_abm_ch)))
stopifnot(assertthat::are_equal(sum(na.omit(lfc < 0)) + sum(na.omit(0 < lfc)), n))GO enrichment with goseq::goseq
Compute
To find enriched GO terms with goseq we need to provide it with a list of genes as a binary vector (the fold-change is ignored).
enriched_terms <-
list(
up = (lfc > 0),
dn = (lfc < 0)
) %>% lapply(
unattach_packages(function(.) {
(.) %>%
# goseq-friendly:
tidyr::replace_na(replace = 0) %>%
# goseq::nullp takes a binary vector (DE=1, else=0) named by gene IDs
goseq::nullp(
id = "geneSymbol",
genome = "mm39",
bias.data = (
# Gene -to- transcript length
# Make a vector of same shape as the named indicator vector
# Note: this could return NULL if the column does not exist,
# so goseq will attempt to download the info from UCSC
tx$mm[names(.), "tx_len"]
),
plot.fit = FALSE
) %>%
# forward the resulting "Probability Weighting Function"
goseq::goseq(id = "geneSymbol", genome = "mm39") %>%
# need to be explicit because goseq does `library(...)`
dplyr::slice_min(over_represented_pvalue, n = 500) %>%
dplyr::select(-under_represented_pvalue) %>%
#
dplyr::mutate(category = link(category)) %>%
dplyr::rename(
# new name = old name
"overrep.p" = over_represented_pvalue,
"#cat" = numInCat,
"#DE" = numDEInCat,
"GO" = ontology
)
})
)Up
enriched_terms$up %>% show_go()Down
enriched_terms$dn %>% show_go()GO enrichment with limma::goana
Compute
enriched_terms <-
list(
up = (lfc > 0),
dn = (lfc < 0)
) %>% lapply(
function(.) {
data.frame(take = .) %>%
filter(take == TRUE) %>%
ind2col("symbol") %>%
# Get the Entrez IDs of our genes
merge(
org.Mm.eg.db::org.Mm.egSYMBOL2EG %>% as.data.frame(),
by = "symbol"
) %>%
pull(gene_id) %>%
# Enrichment analysis
limma::goana(species = "Mm") %>%
# Select topmost
dplyr::slice_min(P.DE, n = 500) %>%
# Format for display
ind2col("category") %>%
dplyr::mutate(category = link(category)) %>%
dplyr::rename(
# new name = old name
"overrep.p" = P.DE,
"#DE" = DE,
"#cat" = N,
"GO" = Ont,
"term" = Term
)
}
)Up
enriched_terms$up %>% show_go()Down
enriched_terms$dn %>% show_go()KEGG enrichment with limma::kegga
Compute
enriched_terms <-
list(
up = (lfc > 0),
dn = (lfc < 0)
) %>% lapply(
function(.) {
data.frame(take = .) %>%
filter(take == TRUE) %>%
ind2col("symbol") %>%
# Get the Entrez IDs of our genes
merge(
org.Mm.eg.db::org.Mm.egSYMBOL2EG %>% as.data.frame(),
by = "symbol"
) %>%
pull(gene_id) %>%
# Enrichment analysis
limma::kegga(species = "Mm") %>%
# https://www.bioconductor.org/packages/devel/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html#pathway-analysis
limma::topKEGG(n = 30, truncate = 34)
}
)Up
enriched_terms$up## Pathway N DE P.DE
## path:mmu03010 Ribosome 179 23 1.479474e-28
## path:mmu05171 Coronavirus disease - COVID-19 247 23 3.190087e-25
## path:mmu04914 Progesterone-mediated oocyte ma... 90 3 9.178243e-03
## path:mmu04066 HIF-1 signaling pathway 114 3 1.734642e-02
## path:mmu04015 Rap1 signaling pathway 214 4 1.918736e-02
## path:mmu04730 Long-term depression 60 2 3.372092e-02
## path:mmu04520 Adherens junction 71 2 4.580017e-02
## path:mmu01521 EGFR tyrosine kinase inhibitor ... 79 2 5.543745e-02
## path:mmu04612 Antigen processing and presenta... 90 2 6.973006e-02
## path:mmu05215 Prostate cancer 99 2 8.222306e-02
## path:mmu04151 PI3K-Akt signaling pathway 359 4 9.347831e-02
## path:mmu04014 Ras signaling pathway 232 3 1.006163e-01
## path:mmu04670 Leukocyte transendothelial migr... 118 2 1.105897e-01
## path:mmu04071 Sphingolipid signaling pathway 124 2 1.200261e-01
## path:mmu04726 Serotonergic synapse 131 2 1.312814e-01
## path:mmu04915 Estrogen signaling pathway 134 2 1.361796e-01
## path:mmu04371 Apelin signaling pathway 137 2 1.411192e-01
## path:mmu04140 Autophagy - animal 138 2 1.427746e-01
## path:mmu04550 Signaling pathways regulating p... 140 2 1.460982e-01
## path:mmu04150 mTOR signaling pathway 156 2 1.732332e-01
## path:mmu04514 Cell adhesion molecules 174 2 2.046671e-01
## path:mmu05030 Cocaine addiction 48 1 2.076260e-01
## path:mmu00600 Sphingolipid metabolism 50 1 2.152938e-01
## path:mmu04360 Axon guidance 181 2 2.170748e-01
## path:mmu04978 Mineral absorption 53 1 2.266598e-01
## path:mmu04923 Regulation of lipolysis in adip... 56 1 2.378649e-01
## path:mmu04510 Focal adhesion 201 2 2.528657e-01
## path:mmu05134 Legionellosis 61 1 2.561887e-01
## path:mmu04213 Longevity regulating pathway - ... 62 1 2.598015e-01
## path:mmu05205 Proteoglycans in cancer 205 2 2.600595e-01Down
enriched_terms$dn## Pathway N DE P.DE
## path:mmu04721 Synaptic vesicle cycle 77 2 0.0004369922
## path:mmu04724 Glutamatergic synapse 113 2 0.0009399664
## path:mmu05016 Huntington disease 302 2 0.0065605660
## path:mmu04964 Proximal tubule bicarbonate rec... 22 1 0.0098372738
## path:mmu04960 Aldosterone-regulated sodium re... 38 1 0.0169459349
## path:mmu04973 Carbohydrate digestion and abso... 48 1 0.0213693671
## path:mmu04978 Mineral absorption 53 1 0.0235754768
## path:mmu04961 Endocrine and other factor-regu... 61 1 0.0270974915
## path:mmu04918 Thyroid hormone synthesis 74 1 0.0328004331
## path:mmu04971 Gastric acid secretion 75 1 0.0332380800
## path:mmu04911 Insulin secretion 86 1 0.0380423993
## path:mmu04970 Salivary secretion 86 1 0.0380423993
## path:mmu04260 Cardiac muscle contraction 87 1 0.0384782659
## path:mmu04976 Bile secretion 100 1 0.0441310636
## path:mmu04925 Aldosterone synthesis and secre... 102 1 0.0449985074
## path:mmu04974 Protein digestion and absorptio... 108 1 0.0475972968
## path:mmu04972 Pancreatic secretion 114 1 0.0501907786
## path:mmu04919 Thyroid hormone signaling pathw... 120 1 0.0527789601
## path:mmu04261 Adrenergic signaling in cardiom... 152 1 0.0664933761
## path:mmu04022 cGMP-PKG signaling pathway 173 1 0.0754121655
## path:mmu04514 Cell adhesion molecules 174 1 0.0758352694
## path:mmu04024 cAMP signaling pathway 215 1 0.0930578679
## path:mmu05014 Amyotrophic lateral sclerosis 369 1 0.1556060316
## path:mmu00472 D-Arginine and D-ornithine meta... 1 0 1.0000000000
## path:mmu00780 Biotin metabolism 3 0 1.0000000000
## path:mmu00471 D-Glutamine and D-glutamate met... 4 0 1.0000000000
## path:mmu00785 Lipoic acid metabolism 4 0 1.0000000000
## path:mmu00290 Valine, leucine and isoleucine ... 4 0 1.0000000000
## path:mmu00524 Neomycin, kanamycin and gentami... 5 0 1.0000000000
## path:mmu00232 Caffeine metabolism 6 0 1.0000000000