Heatmap colours

[ncborcherding]

Heatmap colours

[NGuslund]

I fell down a rabbit hole of updating modules and so on, but managed to get it to work. Thanks for your help Nick!

[NGuslund]

This is the script I ran to produce the (bad) heatmap:
pbmc_small <- get("pbmc_small")
data("pbmc_small")
pbmc_small <- suppressMessages(UpdateSeuratObject(pbmc_small))
ES.seurat <- enrichIt(obj = pbmc_small, gene.sets = GS.hallmark, groups = 1000, cores = 2)
pbmc_small <- Seurat::AddMetaData(pbmc_small, ES.seurat)
colors <- colorRampPalette(c("#0348A6", "#7AC5FF", "#C6FDEC", "#FFB433", "#FF4B20"))
dittoHeatmap(pbmc_small, genes = NULL, metas = names(ES.seurat),
annot.by = "groups",
fontsize = 7,
cluster_cols = TRUE,
heatmap.colors = colors(50))

[ncborcherding]

Hey NGuslund,

Thanks for reaching out and the excellent summary of your issue. Would you mind giving me the output of sessioninfo(). I've actually seen this issue before from my end and if I remember correctly - it has to do with the version of dittoSeq - which is a separate package that my co-author Jared help develop.

Nick

[NGuslund]

Thank you for your quick response!

sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

Random number generation:
RNG: Mersenne-Twister
Normal: Inversion
Sample: Rounding

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base

other attached packages:
[1] tidyseurat_0.3.0 SeuratObject_4.0.2 Seurat_4.0.5 SingleCellExperiment_1.12.0 SummarizedExperiment_1.20.0
[6] Biobase_2.50.0 GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1 S4Vectors_0.28.1
[11] BiocGenerics_0.36.1 MatrixGenerics_1.2.1 matrixStats_0.61.0 dittoSeq_1.2.6 ggplot2_3.3.5
[16] escape_1.0.1 ReactomeGSA_1.4.2

loaded via a namespace (and not attached):
[1] snow_0.4-3 plyr_1.8.6 igraph_1.2.7 lazyeval_0.2.2 GSEABase_1.52.1 splines_4.0.4
[7] BiocParallel_1.24.1 listenv_0.8.0 scattermore_0.7 digest_0.6.28 htmltools_0.5.2 fansi_0.5.0
[13] magrittr_2.0.1 memoise_2.0.0 tensor_1.5 cluster_2.1.2 ROCR_1.0-11 limma_3.46.0
[19] globals_0.14.0 annotate_1.68.0 spatstat.sparse_2.0-0 colorspace_2.0-2 blob_1.2.2 ggrepel_0.9.1
[25] xfun_0.27 dplyr_1.0.7 crayon_1.4.1 RCurl_1.98-1.5 jsonlite_1.7.2 graph_1.68.0
[31] spatstat.data_2.1-0 survival_3.2-13 zoo_1.8-9 glue_1.4.2 polyclip_1.10-0 gtable_0.3.0
[37] zlibbioc_1.36.0 XVector_0.30.0 leiden_0.3.9 DelayedArray_0.16.3 future.apply_1.8.1 msigdbr_7.4.1
[43] abind_1.4-5 scales_1.1.1 pheatmap_1.0.12 edgeR_3.32.1 DBI_1.1.1 miniUI_0.1.1.1
[49] Rcpp_1.0.7 viridisLite_0.4.0 xtable_1.8-4 reticulate_1.22 spatstat.core_2.3-0 bit_4.0.4
[55] GSVA_1.38.2 htmlwidgets_1.5.4 httr_1.4.2 gplots_3.1.1 RColorBrewer_1.1-2 ellipsis_0.3.2
[61] ica_1.0-2 farver_2.1.0 pkgconfig_2.0.3 XML_3.99-0.8 uwot_0.1.10 dbplyr_2.1.1
[67] deldir_1.0-5 locfit_1.5-9.4 utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1 rlang_0.4.12
[73] reshape2_1.4.4 later_1.3.0 AnnotationDbi_1.52.0 munsell_0.5.0 tools_4.0.4 cachem_1.0.6
[79] cli_3.0.1 generics_0.1.0 RSQLite_2.2.8 ggridges_0.5.3 stringr_1.4.0 fastmap_1.1.0
[85] goftest_1.2-3 babelgene_21.4 knitr_1.36 bit64_4.0.5 fitdistrplus_1.1-6 caTools_1.18.2
[91] purrr_0.3.4 RANN_2.6.1 pbapply_1.5-0 future_1.22.1 nlme_3.1-153 mime_0.12
[97] rstudioapi_0.13 compiler_4.0.4 plotly_4.10.0 png_0.1-7 spatstat.utils_2.2-0 tibble_3.1.5
[103] stringi_1.7.5 lattice_0.20-45 Matrix_1.3-4 vctrs_0.3.8 pillar_1.6.4 lifecycle_1.0.1
[109] BiocManager_1.30.16 spatstat.geom_2.3-0 lmtest_0.9-38 RcppAnnoy_0.0.19 data.table_1.14.2 cowplot_1.1.1
[115] bitops_1.0-7 irlba_2.3.3 httpuv_1.6.3 patchwork_1.1.1 R6_2.5.1 promises_1.2.0.1
[121] KernSmooth_2.23-20 gridExtra_2.3 parallelly_1.28.1 codetools_0.2-18 MASS_7.3-54 gtools_3.9.2
[127] assertthat_0.2.1 withr_2.4.2 sctransform_0.3.2 GenomeInfoDbData_1.2.4 mgcv_1.8-38 grid_4.0.4
[133] rpart_4.1-15 tidyr_1.1.4 Rtsne_0.15 shiny_1.7.1

[ncborcherding]

Hey NGusland,

So you have the current version of dittoSeq from BioConductor, but not the development version, which you can get using the code below.

devtools::install_github("dtm2451/dittoSeq")

After restarting your R session, let me know if the problem persists and then I can help troubleshoot.

Thanks,
Nick

[NGuslund]

Thanks Nick, that worked!

I have a couple of additional questions.
In the vignette the hallmark gene set is used:
## We'll use the HALLMARK gene set collection ("H")
GS.hallmark <- getGeneSets(library = "H")

and then furthur down specific pathways are selected:
multi_dittoPlot(sce, vars = c("HALLMARK_APOPTOSIS", "HALLMARK_DNA_REPAIR", "HALLMARK_P53_PATHWAY"),
group.by = "groups", plots = c("jitter", "vlnplot", "boxplot"),
ylab = "Enrichment Scores",
theme = theme_classic() + theme(plot.title = element_text(size = 10)))

I want to use the reactome pathway database.
I selected it like this, but I think this has not selected only the reactome database but several, how can I be more specific?
getGeneSets(library = "C2")

And then I want to visualise these pathways:
multi_dittoPlot(Macs_seurat_PreImmvsVacD1, vars = c("R-HSA-983170", "R-HSA-388841", "R-HSA-2132295"),
group.by = "treatment", plots = c("jitter", "vlnplot", "boxplot"),
ylab = "Enrichment Scores",
theme = theme_classic() + theme(plot.title = element_text(size = 10)))

But I get this error:
Error in .var_OR_get_meta_or_gene(main.var, object, assay, slot, adjustment) : 'var' is not a metadata or gene nor equal in length to ncol('object')

Additionally, I am would like to carry out a similar analysis as this:

Are there any vignettes I can follow? I am not even sure of the name of this type of GSEA so it is hard to search for :)

[ncborcherding]

Hey NGuslund,

These actually are great questions and areas of development for escape. So I'd direct you to the development version which has some of the things incoporated:

devtools::install_github("ncborcherding/escape@dev")

I want to use the reactome pathway database.
I selected it like this, but I think this has not selected only the reactome database but several, how can I be more specific?
getGeneSets(library = "C2")

The new getGeneSets() has a parameter for subcategories -

getGeneSets(species = "Homo sapiens", 
                        library = "C2, 
                        subcategory = "REACTOME")

And then I want to visualise these pathways:
multi_dittoPlot(Macs_seurat_PreImmvsVacD1, vars = c("R-HSA-983170", "R-HSA-388841", "R-HSA-2132295"),
group.by = "treatment", plots = c("jitter", "vlnplot", "boxplot"),
ylab = "Enrichment Scores",
theme = theme_classic() + theme(plot.title = element_text(size = 10)))

But I get this error:
Error in .var_OR_get_meta_or_gene(main.var, object, assay, slot, adjustment) : 'var' is not a metadata or gene nor equal in length to ncol('object')

I would try this with the new version of dittoSeq and see if this problem persists. DittoSeq has wider support for single cell objects in the newer versions.

Additionally, I am would like to carry out a similar analysis as this:

Just got done adding the rank plot feature here. Still need to develop the significance testing strategy for this, but it is at least generating the enrichment density between 2 or more groups.

Nick

[NGuslund]

Thank you Nick, I have run this:
devtools::install_github("ncborcherding/escape@dev")

However, this script below brings an error message
GS.hallmark <- getGeneSets(species = "Homo sapiens",
library = "C2",
subcategory = "REACTOME")

Error in getGeneSets(species = "Homo sapiens", library = "C2", subcategory = "REACTOME") :
unused argument (subcategory = "REACTOME")

[ncborcherding]

Hey NGusland,

Did you restart your R session? The error you list indicates that there is no parameter "subcategory". Does the documentation page have the subcategory parameter?

You can see this by using:

??getGeneSets

Nick