## ----include = FALSE-------------------------------------------------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE-------------------------------------------------------- ## For links library("BiocStyle") ## Track time spent on making the vignette startTime <- Sys.time() ## Bib setup library("RefManageR") ## Write bibliography information bib <- c( R = citation(), AnnotationHub = citation("AnnotationHub")[1], benchmarkme = citation("benchmarkme")[1], BiocFileCache = citation("BiocFileCache")[1], BiocGenerics = citation("BiocGenerics")[1], BiocStyle = citation("BiocStyle")[1], cowplot = citation("cowplot")[1], DT = citation("DT")[1], edgeR = citation("edgeR")[1], ExperimentHub = citation("ExperimentHub")[1], fields = citation("fields")[1], GenomicRanges = citation("GenomicRanges")[1], ggplot2 = citation("ggplot2")[1], golem = citation("golem")[1], IRanges = citation("IRanges")[1], knitr = citation("knitr")[3], limma = citation("limma")[1], magick = citation("magick")[1], Matrix = citation("Matrix")[1], paletteer = citation("paletteer")[1], plotly = citation("plotly")[1], RColorBrewer = citation("RColorBrewer")[1], RefManageR = citation("RefManageR")[1], rmarkdown = citation("rmarkdown")[1], rtracklayer = citation("rtracklayer")[1], S4Vectors = citation("S4Vectors")[1], scater = citation("scater")[1], scuttle = citation("scuttle")[1], sessioninfo = citation("sessioninfo")[1], SingleCellExperiment = citation("SingleCellExperiment")[1], shiny = citation("shiny")[1], SpatialExperiment = citation("SpatialExperiment")[1], spatialLIBD = citation("spatialLIBD")[1], spatialLIBDpaper = citation("spatialLIBD")[2], spatialDLPFC = citation("spatialLIBD")[3], VisiumSPGAD = citation("spatialLIBD")[4], statmod = citation("statmod")[1], SummarizedExperiment = citation("SummarizedExperiment")[1], testthat = citation("testthat")[1], viridisLite = citation("viridisLite")[1] ) ## ----'install', eval = FALSE------------------------------------------------------------------------------------------ # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # # BiocManager::install("spatialLIBD") # # ## Check that you have a valid Bioconductor installation # BiocManager::valid() ## ----"install_deps", eval = FALSE------------------------------------------------------------------------------------- # BiocManager::install("spatialLIBD", dependencies = TRUE, force = TRUE) ## ----"install_devel", eval = FALSE------------------------------------------------------------------------------------ # BiocManager::install("LieberInstitute/spatialLIBD") ## ----'citation'------------------------------------------------------------------------------------------------------- ## Citation info citation("spatialLIBD") ## ----setup, message = FALSE, warning = FALSE-------------------------------------------------------------------------- library("spatialLIBD") ## ----'experiment_hub'------------------------------------------------------------------------------------------------- ## Connect to ExperimentHub ehub <- ExperimentHub::ExperimentHub() ## ----'download_data'-------------------------------------------------------------------------------------------------- ## Download the small example sce data sce <- fetch_data(type = "sce_example", eh = ehub) ## Convert to a SpatialExperiment object spe <- sce_to_spe(sce) ## If you want to download the full real data (about 2.1 GB in RAM) use: if (FALSE) { if (!exists("spe")) spe <- fetch_data(type = "spe", eh = ehub) } ## Query ExperimentHub and download the data if (!exists("sce_layer")) sce_layer <- fetch_data(type = "sce_layer", eh = ehub) modeling_results <- fetch_data("modeling_results", eh = ehub) ## ----'explore data'--------------------------------------------------------------------------------------------------- ## spot-level data spe ## layer-level data sce_layer ## list of modeling result tables sapply(modeling_results, class) sapply(modeling_results, dim) sapply(modeling_results, function(x) { head(colnames(x)) }) ## ----'create_sig_genes'----------------------------------------------------------------------------------------------- ## Convert to a long format the modeling results ## This takes a few seconds to run system.time( sig_genes <- sig_genes_extract_all( n = nrow(sce_layer), modeling_results = modeling_results, sce_layer = sce_layer ) ) ## Explore the result class(sig_genes) dim(sig_genes) ## --------------------------------------------------------------------------------------------------------------------- if (interactive()) { run_app( spe = spe, sce_layer = sce_layer, modeling_results = modeling_results, sig_genes = sig_genes ) } ## ----'vis_clus', fig.small = TRUE------------------------------------------------------------------------------------- ## View our LIBD layers for one sample vis_clus( spe = spe, clustervar = "layer_guess_reordered", sampleid = "151673", colors = libd_layer_colors, ... = " LIBD Layers" ) ## ----'vis_clus_variables'--------------------------------------------------------------------------------------------- ## This is not fully precise, but gives you a rough idea ## Some integer columns are actually continuous variables table(sapply(colData(spe), class) %in% c("factor", "integer")) ## This is more precise (one cluster has 28 unique values) table(sapply(colData(spe), function(x) length(unique(x))) < 29) ## ----'vis_clus_no_spatial', fig.small = TRUE-------------------------------------------------------------------------- ## View our LIBD layers for one sample ## without spatial information vis_clus( spe = spe, clustervar = "layer_guess_reordered", sampleid = "151673", colors = libd_layer_colors, ... = " LIBD Layers", spatial = FALSE ) ## ----'libd_layer_colors'---------------------------------------------------------------------------------------------- ## Color palette designed by Lukas M. Weber with feedback from the team. libd_layer_colors ## ----'vis_gene', fig.small = TRUE------------------------------------------------------------------------------------- ## Available gene expression assays assayNames(spe) ## Not all of these make sense to visualize ## In particular, the key is not useful to visualize. colnames(colData(spe)) ## Visualize a gene vis_gene( spe = spe, sampleid = "151673", viridis = FALSE ) ## Visualize the estimated number of cells per spot vis_gene( spe = spe, sampleid = "151673", geneid = "cell_count" ) ## Visualize the fraction of chrM expression per spot ## without the spatial layer vis_gene( spe = spe, sampleid = "151673", geneid = "expr_chrM_ratio", spatial = FALSE ) ## ----'sig_genes_detail'----------------------------------------------------------------------------------------------- head(sig_genes) ## ----'layer_boxplot'-------------------------------------------------------------------------------------------------- ## Note that we recommend setting the random seed so the jittering of the ## points will be reproducible. Given the requirements by BiocCheck this ## cannot be done inside the layer_boxplot() function. ## Create a boxplot of the first gene in `sig_genes`. set.seed(20200206) layer_boxplot(sig_genes = sig_genes, sce_layer = sce_layer) ## Viridis colors displayed in the shiny app ## showing the first pairwise model result ## (which illustrates the background colors used for the layers not ## involved in the pairwise comparison) set.seed(20200206) layer_boxplot( i = which(sig_genes$model_type == "pairwise")[1], sig_genes = sig_genes, sce_layer = sce_layer, col_low_box = viridisLite::viridis(4)[2], col_low_point = viridisLite::viridis(4)[1], col_high_box = viridisLite::viridis(4)[3], col_high_point = viridisLite::viridis(4)[4] ) ## Paper colors displayed in the shiny app set.seed(20200206) layer_boxplot( sig_genes = sig_genes, sce_layer = sce_layer, short_title = FALSE, col_low_box = "palegreen3", col_low_point = "springgreen2", col_high_box = "darkorange2", col_high_point = "orange1" ) ## ----'matt_t_stats'--------------------------------------------------------------------------------------------------- ## Explore the enrichment t-statistics derived from Tran et al's snRNA-seq ## DLPFC data dim(tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer) tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer[seq_len(3), seq_len(6)] ## ----'layer_stat_cor'------------------------------------------------------------------------------------------------- ## Compute the correlation matrix of enrichment t-statistics between our data ## and Tran et al's snRNA-seq data cor_stats_layer <- layer_stat_cor( tstats_Human_DLPFC_snRNAseq_Nguyen_topLayer, modeling_results, "enrichment" ) ## Explore the correlation matrix head(cor_stats_layer[, seq_len(3)]) ## ----'layer_stat_cor_plot', fig.wide = TRUE--------------------------------------------------------------------------- ## Visualize the correlation matrix layer_stat_cor_plot(cor_stats_layer, max = max(cor_stats_layer)) ## ----'read_sfari'----------------------------------------------------------------------------------------------------- ## Read in the SFARI gene sets included in the package asd_sfari <- utils::read.csv( system.file( "extdata", "SFARI-Gene_genes_01-03-2020release_02-04-2020export.csv", package = "spatialLIBD" ), as.is = TRUE ) ## Format them appropriately asd_sfari_geneList <- list( Gene_SFARI_all = asd_sfari$ensembl.id, Gene_SFARI_high = asd_sfari$ensembl.id[asd_sfari$gene.score < 3], Gene_SFARI_syndromic = asd_sfari$ensembl.id[asd_sfari$syndromic == 1] ) ## ----'sfari_enrichment'----------------------------------------------------------------------------------------------- ## Compute the gene set enrichment results asd_sfari_enrichment <- gene_set_enrichment( gene_list = asd_sfari_geneList, modeling_results = modeling_results, model_type = "enrichment" ) ## Explore the results head(asd_sfari_enrichment) ## ----'sfari_enrichment_plot'------------------------------------------------------------------------------------------ ## Visualize gene set enrichment results gene_set_enrichment_plot( asd_sfari_enrichment, xlabs = gsub(".*_", "", unique(asd_sfari_enrichment$ID)), layerHeights = c(0, 40, 55, 75, 85, 110, 120, 135) ) ## ----'run_check_functions'-------------------------------------------------------------------------------------------- check_spe(spe) check_sce_layer(sce_layer) ## The output here is too long to print xx <- check_modeling_results(modeling_results) identical(xx, modeling_results) ## ----createVignette, eval=FALSE--------------------------------------------------------------------------------------- # ## Create the vignette # library("rmarkdown") # system.time(render("spatialLIBD.Rmd")) # # ## Extract the R code # library("knitr") # knit("spatialLIBD.Rmd", tangle = TRUE) ## ----reproduce1, echo=FALSE------------------------------------------------------------------------------------------- ## Date the vignette was generated Sys.time() ## ----reproduce2, echo=FALSE------------------------------------------------------------------------------------------- ## Processing time in seconds totalTime <- diff(c(startTime, Sys.time())) round(totalTime, digits = 3) ## ----reproduce3, echo=FALSE------------------------------------------------------------------------------------------- ## Session info library("sessioninfo") options(width = 120) session_info() ## ----vignetteBiblio, results = 'asis', echo = FALSE, warning = FALSE, message = FALSE--------------------------------- ## Print bibliography PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))