From a0b1e1ae7691b7fb9398dc55f61a104b25f6683a Mon Sep 17 00:00:00 2001 From: geertvangeest Date: Mon, 29 Apr 2024 15:09:31 +0200 Subject: [PATCH] adds quarto website --- .Rprofile | 1 + .github/workflows/publish.yml | 47 + .gitignore | 14 + _quarto.yml | 48 + assets/SIB_LogoQ_GBv.svg | 1 + assets/SIB_logo.svg | 138 + custom.scss | 13 + index.qmd | 3 + normalization_scaling.qmd | 0 ...s_tutorial_-_spatial_transcriptomics.Rproj | 13 + quality_control.qmd | 28 + renv.lock | 2513 +++++++++++++++++ renv/.gitignore | 7 + renv/activate.R | 1180 ++++++++ renv/settings.json | 19 + scripts/spatial_RNAseq-clustering.R | 763 +++++ .../spatial_RNAseq_integration-clustering.R | 763 +++++ 17 files changed, 5551 insertions(+) create mode 100644 .Rprofile create mode 100644 .github/workflows/publish.yml create mode 100644 _quarto.yml create mode 100644 assets/SIB_LogoQ_GBv.svg create mode 100644 assets/SIB_logo.svg create mode 100644 custom.scss create mode 100644 index.qmd create mode 100644 normalization_scaling.qmd create mode 100644 p907_SIBdays_tutorial_-_spatial_transcriptomics.Rproj create mode 100644 quality_control.qmd create mode 100644 renv.lock create mode 100644 renv/.gitignore create mode 100644 renv/activate.R create mode 100644 renv/settings.json create mode 100644 scripts/spatial_RNAseq-clustering.R create mode 100644 scripts/spatial_RNAseq_integration-clustering.R diff --git a/.Rprofile b/.Rprofile new file mode 100644 index 0000000..81b960f --- /dev/null +++ b/.Rprofile @@ -0,0 +1 @@ +source("renv/activate.R") diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100644 index 0000000..1367237 --- /dev/null +++ b/.github/workflows/publish.yml @@ -0,0 +1,47 @@ +# on: +# # manually trigger +# workflow_dispatch: +# # pushing to main +# push: +# branches: main +# paths: +# - *.qmd +# - custom.scss +# - renv.lock +# - _quarto.yml +# # every week on sunday +# schedule: +# - cron: "0 0 * * SUN" +# +# name: Quarto Publish +# +# jobs: +# build-deploy: +# runs-on: ubuntu-latest +# permissions: +# contents: write +# steps: +# - name: Check out repository +# uses: actions/checkout@v4 +# +# - name: Set up Quarto +# uses: quarto-dev/quarto-actions/setup@v2 +# +# - name: Install R +# uses: r-lib/actions/setup-r@v2 +# with: +# r-version: '4.3.1' +# +# - name: Install R Dependencies +# uses: r-lib/actions/setup-renv@v2 +# with: +# cache-version: 1 +# +# - name: Render and Publish +# uses: quarto-dev/quarto-actions/publish@v2 +# with: +# target: gh-pages +# env: +# GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} +# PAT: ${{ secrets.PAT }} +# GOOGLE_API_KEY: ${{ secrets.GOOGLE_API_KEY }} \ No newline at end of file diff --git a/.gitignore b/.gitignore index e75435c..db63e81 100644 --- a/.gitignore +++ b/.gitignore @@ -47,3 +47,17 @@ po/*~ # RStudio Connect folder rsconnect/ + +/.quarto/ + +# rstudio server files +var-lib-rstudio-server/ +run/ +database.conf + +# quarto large dirs +_freeze/ +_site/ + +# analysis dir +analysis/ diff --git a/_quarto.yml b/_quarto.yml new file mode 100644 index 0000000..519d639 --- /dev/null +++ b/_quarto.yml @@ -0,0 +1,48 @@ +project: + type: website + execute-dir: project + render: + - "*.qmd" + - "!CONTRIBUTIONG.md" + - "!LICENSE.md" + +website: + title: "Spatial transcriptomics" + favicon: assets/SIB_logo.svg + search: true + page-navigation: true + navbar: + background: "#003eaa" + left: + - href: index.qmd + text: Home + - href: precourse_preparations.qmd + - href: course_schedule.qmd + - text: "Course material" + menu: + - href: "quality_control.qmd" + - href: "normalization_scaling.qmd" + right: + - icon: github + href: https://github.com/sib-swiss/spatial-transcriptomics-training/ + aria-label: GitHub + sidebar: + logo: assets/SIB_LogoQ_GBv.svg + +format: + html: + code-link: true + theme: + - default + - custom.scss + toc: true + header-includes: | + + +execute: + freeze: auto + cache: true + tidy: true + + + diff --git a/assets/SIB_LogoQ_GBv.svg b/assets/SIB_LogoQ_GBv.svg new file mode 100644 index 0000000..14bac8b --- /dev/null +++ b/assets/SIB_LogoQ_GBv.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/assets/SIB_logo.svg b/assets/SIB_logo.svg new file mode 100644 index 0000000..b5f2113 --- /dev/null +++ b/assets/SIB_logo.svg @@ -0,0 +1,138 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/custom.scss b/custom.scss new file mode 100644 index 0000000..39787fd --- /dev/null +++ b/custom.scss @@ -0,0 +1,13 @@ +@import url('https://fonts.cdnfonts.com/css/cabinet-grotesk'); +@import url('https://fonts.googleapis.com/css?family=Inter'); + +/*-- scss:defaults --*/ +// Base document colors +$headings-font-family: 'Cabinet grotesk'; +$font-family-sans-serif: 'Inter'; +$code-color: #ae191a; +$link-color: #009ee3; +$navbar-bg: #ae191a; +$kbd-color: #ae191a; +$kbd-bg: #ae191a; +$pre-color: #747473; diff --git a/index.qmd b/index.qmd new file mode 100644 index 0000000..250972f --- /dev/null +++ b/index.qmd @@ -0,0 +1,3 @@ +# SIB days tutorial: introduction to spatial transcriptomics + +Lorem ipsum \ No newline at end of file diff --git a/normalization_scaling.qmd b/normalization_scaling.qmd new file mode 100644 index 0000000..e69de29 diff --git a/p907_SIBdays_tutorial_-_spatial_transcriptomics.Rproj b/p907_SIBdays_tutorial_-_spatial_transcriptomics.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/p907_SIBdays_tutorial_-_spatial_transcriptomics.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/quality_control.qmd b/quality_control.qmd new file mode 100644 index 0000000..356d698 --- /dev/null +++ b/quality_control.qmd @@ -0,0 +1,28 @@ + +# Mapping + counting (spaceranger) + +we downloaded the data from https://www.10xgenomics.com/datasets/fresh-frozen-visium-on-cytassist-mouse-brain-probe-based-whole-transcriptome-profiling-2-standard + +Fresh Frozen Visium on CytAssist: Mouse Brain, Probe-Based Whole Transcriptome Profiling + +```{r} +#| output: false +library(Seurat) +library(umap) +library(ggpubr) +library(hdf5r) +``` + +```{r} +countFolder <- "raw_data/Posterior" +analysisFolder <- "analysis" +dir.create(analysisFolder, showWarnings = FALSE) +``` + +```{r} +seu <- Load10X_Spatial(countFolder, + filename = "filtered_feature_bc_matrix.h5", + slice = "Posterior") +``` + + diff --git a/renv.lock b/renv.lock new file mode 100644 index 0000000..e34b2e4 --- /dev/null +++ b/renv.lock @@ -0,0 +1,2513 @@ +{ + "R": { + "Version": "4.3.1", + "Repositories": [ + { + "Name": "CRAN", + "URL": "https://p3m.dev/cran/2023-10-30" + } + ] + }, + "Packages": { + "BH": { + "Package": "BH", + "Version": "1.84.0-0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a8235afbcd6316e6e91433ea47661013" + }, + "FNN": { + "Package": "FNN", + "Version": "1.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "eaabdc7938aa3632a28273f53a0d226d" + }, + "KernSmooth": { + "Package": "KernSmooth", + "Version": "2.23-21", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats" + ], + "Hash": "6314fc110e09548ba889491db6ae67fb" + }, + "MASS": { + "Package": "MASS", + "Version": "7.3-60", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "methods", + "stats", + "utils" + ], + "Hash": "a56a6365b3fa73293ea8d084be0d9bb0" + }, + "Matrix": { + "Package": "Matrix", + "Version": "1.6-5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "grid", + "lattice", + "methods", + "stats", + "utils" + ], + "Hash": "8c7115cd3a0e048bda2a7cd110549f7a" + }, + "MatrixModels": { + "Package": "MatrixModels", + "Version": "0.5-3", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "Matrix", + "R", + "methods", + "stats" + ], + "Hash": "0776bf7526869e0286b0463cb72fb211" + }, + "R6": { + "Package": "R6", + "Version": "2.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "470851b6d5d0ac559e9d01bb352b4021" + }, + "RANN": { + "Package": "RANN", + "Version": "2.6.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d128ea05a972d3e67c6f39de52c72bd7" + }, + "RColorBrewer": { + "Package": "RColorBrewer", + "Version": "1.1-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "45f0398006e83a5b10b72a90663d8d8c" + }, + "ROCR": { + "Package": "ROCR", + "Version": "1.0-11", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "gplots", + "grDevices", + "graphics", + "methods", + "stats" + ], + "Hash": "cc151930e20e16427bc3d0daec62b4a9" + }, + "RSpectra": { + "Package": "RSpectra", + "Version": "0.16-1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "Rcpp", + "RcppEigen" + ], + "Hash": "6b5ab997fd5ff6d46a5f1d9f8b76961c" + }, + "Rcpp": { + "Package": "Rcpp", + "Version": "1.0.12", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "methods", + "utils" + ], + "Hash": "5ea2700d21e038ace58269ecdbeb9ec0" + }, + "RcppAnnoy": { + "Package": "RcppAnnoy", + "Version": "0.0.22", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "methods" + ], + "Hash": "f6baa1e06fb6c3724f601a764266cb0d" + }, + "RcppArmadillo": { + "Package": "RcppArmadillo", + "Version": "0.12.8.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "methods", + "stats", + "utils" + ], + "Hash": "85c0c4491c06db4814ed7e39b960973c" + }, + "RcppEigen": { + "Package": "RcppEigen", + "Version": "0.3.4.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "stats", + "utils" + ], + "Hash": "df49e3306f232ec28f1604e36a202847" + }, + "RcppHNSW": { + "Package": "RcppHNSW", + "Version": "0.6.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Rcpp", + "methods" + ], + "Hash": "1f2dc32c27746a35196aaf95adb357be" + }, + "RcppProgress": { + "Package": "RcppProgress", + "Version": "0.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "1c0aa18b97e6aaa17f93b8b866c0ace5" + }, + "RcppTOML": { + "Package": "RcppTOML", + "Version": "0.2.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp" + ], + "Hash": "c232938949fcd8126034419cc529333a" + }, + "Rtsne": { + "Package": "Rtsne", + "Version": "0.17", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Rcpp", + "stats" + ], + "Hash": "f81f7764a3c3e310b1d40e1a8acee19e" + }, + "Seurat": { + "Package": "Seurat", + "Version": "5.0.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "KernSmooth", + "MASS", + "Matrix", + "R", + "RANN", + "RColorBrewer", + "ROCR", + "RSpectra", + "Rcpp", + "RcppAnnoy", + "RcppEigen", + "RcppHNSW", + "RcppProgress", + "Rtsne", + "SeuratObject", + "cluster", + "cowplot", + "fastDummies", + "fitdistrplus", + "future", + "future.apply", + "generics", + "ggplot2", + "ggrepel", + "ggridges", + "grDevices", + "graphics", + "grid", + "httr", + "ica", + "igraph", + "irlba", + "jsonlite", + "leiden", + "lifecycle", + "lmtest", + "matrixStats", + "methods", + "miniUI", + "patchwork", + "pbapply", + "plotly", + "png", + "progressr", + "purrr", + "reticulate", + "rlang", + "scales", + "scattermore", + "sctransform", + "shiny", + "spatstat.explore", + "spatstat.geom", + "stats", + "tibble", + "tools", + "utils", + "uwot" + ], + "Hash": "f5325a968489ef65254a51d74dc0c520" + }, + "SeuratObject": { + "Package": "SeuratObject", + "Version": "5.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "Rcpp", + "RcppEigen", + "future", + "future.apply", + "generics", + "grDevices", + "grid", + "lifecycle", + "methods", + "progressr", + "rlang", + "sp", + "spam", + "stats", + "tools", + "utils" + ], + "Hash": "424c2b9e50b053c86fc38f17e09917da" + }, + "SparseM": { + "Package": "SparseM", + "Version": "1.81", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "methods", + "stats", + "utils" + ], + "Hash": "2042cd9759cc89a453c4aefef0ce9aae" + }, + "abind": { + "Package": "abind", + "Version": "1.4-5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods", + "utils" + ], + "Hash": "4f57884290cc75ab22f4af9e9d4ca862" + }, + "askpass": { + "Package": "askpass", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "sys" + ], + "Hash": "cad6cf7f1d5f6e906700b9d3e718c796" + }, + "backports": { + "Package": "backports", + "Version": "1.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "c39fbec8a30d23e721980b8afb31984c" + }, + "base64enc": { + "Package": "base64enc", + "Version": "0.1-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "543776ae6848fde2f48ff3816d0628bc" + }, + "bitops": { + "Package": "bitops", + "Version": "1.0-7", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b7d8d8ee39869c18d8846a184dd8a1af" + }, + "boot": { + "Package": "boot", + "Version": "1.3-28.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "stats" + ], + "Hash": "9a052fbcbe97a98ceb18dbfd30ebd96e" + }, + "brio": { + "Package": "brio", + "Version": "1.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "68bd2b066e1fe780bbf62fc8bcc36de3" + }, + "broom": { + "Package": "broom", + "Version": "1.0.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "backports", + "dplyr", + "ellipsis", + "generics", + "glue", + "lifecycle", + "purrr", + "rlang", + "stringr", + "tibble", + "tidyr" + ], + "Hash": "fd25391c3c4f6ecf0fa95f1e6d15378c" + }, + "bslib": { + "Package": "bslib", + "Version": "0.7.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "base64enc", + "cachem", + "fastmap", + "grDevices", + "htmltools", + "jquerylib", + "jsonlite", + "lifecycle", + "memoise", + "mime", + "rlang", + "sass" + ], + "Hash": "8644cc53f43828f19133548195d7e59e" + }, + "caTools": { + "Package": "caTools", + "Version": "1.18.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "bitops" + ], + "Hash": "34d90fa5845004236b9eacafc51d07b2" + }, + "cachem": { + "Package": "cachem", + "Version": "1.0.8", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "fastmap", + "rlang" + ], + "Hash": "c35768291560ce302c0a6589f92e837d" + }, + "callr": { + "Package": "callr", + "Version": "3.7.6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "processx", + "utils" + ], + "Hash": "d7e13f49c19103ece9e58ad2d83a7354" + }, + "car": { + "Package": "car", + "Version": "3.1-2", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "MASS", + "R", + "abind", + "carData", + "grDevices", + "graphics", + "lme4", + "mgcv", + "nlme", + "nnet", + "pbkrtest", + "quantreg", + "scales", + "stats", + "utils" + ], + "Hash": "839b351f31d56e0147439eb22c00a66a" + }, + "carData": { + "Package": "carData", + "Version": "3.0-5", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R" + ], + "Hash": "ac6cdb8552c61bd36b0e54d07cf2aab7" + }, + "cli": { + "Package": "cli", + "Version": "3.6.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "utils" + ], + "Hash": "1216ac65ac55ec0058a6f75d7ca0fd52" + }, + "cluster": { + "Package": "cluster", + "Version": "2.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "stats", + "utils" + ], + "Hash": "5edbbabab6ce0bf7900a74fd4358628e" + }, + "codetools": { + "Package": "codetools", + "Version": "0.2-19", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "c089a619a7fae175d149d89164f8c7d8" + }, + "colorspace": { + "Package": "colorspace", + "Version": "2.1-0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "methods", + "stats" + ], + "Hash": "f20c47fd52fae58b4e377c37bb8c335b" + }, + "commonmark": { + "Package": "commonmark", + "Version": "1.9.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "5d8225445acb167abf7797de48b2ee3c" + }, + "corrplot": { + "Package": "corrplot", + "Version": "0.92", + "Source": "Repository", + "Repository": "RSPM", + "Hash": "fcf11a91936fd5047b2ee9bc00595e36" + }, + "cowplot": { + "Package": "cowplot", + "Version": "1.1.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "ggplot2", + "grDevices", + "grid", + "gtable", + "methods", + "rlang", + "scales" + ], + "Hash": "8ef2084dd7d28847b374e55440e4f8cb" + }, + "cpp11": { + "Package": "cpp11", + "Version": "0.4.7", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "5a295d7d963cc5035284dcdbaf334f4e" + }, + "crayon": { + "Package": "crayon", + "Version": "1.5.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grDevices", + "methods", + "utils" + ], + "Hash": "e8a1e41acf02548751f45c718d55aa6a" + }, + "crosstalk": { + "Package": "crosstalk", + "Version": "1.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R6", + "htmltools", + "jsonlite", + "lazyeval" + ], + "Hash": "ab12c7b080a57475248a30f4db6298c0" + }, + "curl": { + "Package": "curl", + "Version": "5.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "411ca2c03b1ce5f548345d2fc2685f7a" + }, + "data.table": { + "Package": "data.table", + "Version": "1.15.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "8ee9ac56ef633d0c7cab8b2ca87d683e" + }, + "deldir": { + "Package": "deldir", + "Version": "2.0-4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics" + ], + "Hash": "24754fce82729ff85317dd195b6646a8" + }, + "desc": { + "Package": "desc", + "Version": "1.4.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "cli", + "utils" + ], + "Hash": "99b79fcbd6c4d1ce087f5c5c758b384f" + }, + "diffobj": { + "Package": "diffobj", + "Version": "0.3.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "crayon", + "methods", + "stats", + "tools", + "utils" + ], + "Hash": "bcaa8b95f8d7d01a5dedfd959ce88ab8" + }, + "digest": { + "Package": "digest", + "Version": "0.6.35", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "utils" + ], + "Hash": "698ece7ba5a4fa4559e3d537e7ec3d31" + }, + "dotCall64": { + "Package": "dotCall64", + "Version": "1.1-1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "80f374ef8500fcdc5d84a0345b837227" + }, + "dplyr": { + "Package": "dplyr", + "Version": "1.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "cli", + "generics", + "glue", + "lifecycle", + "magrittr", + "methods", + "pillar", + "rlang", + "tibble", + "tidyselect", + "utils", + "vctrs" + ], + "Hash": "fedd9d00c2944ff00a0e2696ccf048ec" + }, + "dqrng": { + "Package": "dqrng", + "Version": "0.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "BH", + "R", + "Rcpp", + "sitmo" + ], + "Hash": "824df2aeba88d701df5e79018b35b815" + }, + "ellipsis": { + "Package": "ellipsis", + "Version": "0.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "rlang" + ], + "Hash": "bb0eec2fe32e88d9e2836c2f73ea2077" + }, + "evaluate": { + "Package": "evaluate", + "Version": "0.23", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "daf4a1246be12c1fa8c7705a0935c1a0" + }, + "fansi": { + "Package": "fansi", + "Version": "1.0.6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "utils" + ], + "Hash": "962174cf2aeb5b9eea581522286a911f" + }, + "farver": { + "Package": "farver", + "Version": "2.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "8106d78941f34855c440ddb946b8f7a5" + }, + "fastDummies": { + "Package": "fastDummies", + "Version": "1.7.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "data.table", + "stringr", + "tibble" + ], + "Hash": "e0f9c0c051e0e8d89996d7f0c400539f" + }, + "fastmap": { + "Package": "fastmap", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "f7736a18de97dea803bde0a2daaafb27" + }, + "fitdistrplus": { + "Package": "fitdistrplus", + "Version": "1.1-11", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "R", + "grDevices", + "methods", + "stats", + "survival" + ], + "Hash": "f40ef9686e85681a1ccbf33d9236aeb9" + }, + "fontawesome": { + "Package": "fontawesome", + "Version": "0.5.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "htmltools", + "rlang" + ], + "Hash": "c2efdd5f0bcd1ea861c2d4e2a883a67d" + }, + "fs": { + "Package": "fs", + "Version": "1.6.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "47b5f30c720c23999b913a1a635cf0bb" + }, + "future": { + "Package": "future", + "Version": "1.33.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "digest", + "globals", + "listenv", + "parallel", + "parallelly", + "utils" + ], + "Hash": "fd7b1d69d16d0d114e4fa82db68f184c" + }, + "future.apply": { + "Package": "future.apply", + "Version": "1.11.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "future", + "globals", + "parallel", + "utils" + ], + "Hash": "afe1507511629f44572e6c53b9baeb7c" + }, + "generics": { + "Package": "generics", + "Version": "0.1.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "15e9634c0fcd294799e9b2e929ed1b86" + }, + "ggplot2": { + "Package": "ggplot2", + "Version": "3.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "R", + "cli", + "glue", + "grDevices", + "grid", + "gtable", + "isoband", + "lifecycle", + "mgcv", + "rlang", + "scales", + "stats", + "tibble", + "vctrs", + "withr" + ], + "Hash": "52ef83f93f74833007f193b2d4c159a2" + }, + "ggpubr": { + "Package": "ggpubr", + "Version": "0.6.0", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R", + "cowplot", + "dplyr", + "ggplot2", + "ggrepel", + "ggsci", + "ggsignif", + "glue", + "grid", + "gridExtra", + "magrittr", + "polynom", + "purrr", + "rlang", + "rstatix", + "scales", + "stats", + "tibble", + "tidyr", + "utils" + ], + "Hash": "c957612b8bb1ee9ab7b2450d26663e7e" + }, + "ggrepel": { + "Package": "ggrepel", + "Version": "0.9.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "ggplot2", + "grid", + "rlang", + "scales", + "withr" + ], + "Hash": "cc3361e234c4a5050e29697d675764aa" + }, + "ggridges": { + "Package": "ggridges", + "Version": "0.5.6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "ggplot2", + "grid", + "scales", + "withr" + ], + "Hash": "66488692cb8621bc78df1b9b819497a6" + }, + "ggsci": { + "Package": "ggsci", + "Version": "3.0.0", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R", + "ggplot2", + "grDevices", + "scales" + ], + "Hash": "93664e03010c3f4b570c890dda99ade5" + }, + "ggsignif": { + "Package": "ggsignif", + "Version": "0.6.4", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "ggplot2" + ], + "Hash": "a57f0f5dbcfd0d77ad4ff33032f5dc79" + }, + "globals": { + "Package": "globals", + "Version": "0.16.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "codetools" + ], + "Hash": "2580567908cafd4f187c1e5a91e98b7f" + }, + "glue": { + "Package": "glue", + "Version": "1.7.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "e0b3a53876554bd45879e596cdb10a52" + }, + "goftest": { + "Package": "goftest", + "Version": "1.2-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats" + ], + "Hash": "dbe0201f91eeb15918dd3fbf01ee689a" + }, + "gplots": { + "Package": "gplots", + "Version": "3.1.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "KernSmooth", + "R", + "caTools", + "gtools", + "methods", + "stats" + ], + "Hash": "f72b5d1ed587f8905e38ee7295e88d80" + }, + "gridExtra": { + "Package": "gridExtra", + "Version": "2.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grDevices", + "graphics", + "grid", + "gtable", + "utils" + ], + "Hash": "7d7f283939f563670a697165b2cf5560" + }, + "gtable": { + "Package": "gtable", + "Version": "0.3.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "glue", + "grid", + "lifecycle", + "rlang" + ], + "Hash": "b29cf3031f49b04ab9c852c912547eef" + }, + "gtools": { + "Package": "gtools", + "Version": "3.9.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "methods", + "stats", + "utils" + ], + "Hash": "588d091c35389f1f4a9d533c8d709b35" + }, + "here": { + "Package": "here", + "Version": "1.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "rprojroot" + ], + "Hash": "24b224366f9c2e7534d2344d10d59211" + }, + "highr": { + "Package": "highr", + "Version": "0.10", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "xfun" + ], + "Hash": "06230136b2d2b9ba5805e1963fa6e890" + }, + "htmltools": { + "Package": "htmltools", + "Version": "0.5.8.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "base64enc", + "digest", + "fastmap", + "grDevices", + "rlang", + "utils" + ], + "Hash": "81d371a9cc60640e74e4ab6ac46dcedc" + }, + "htmlwidgets": { + "Package": "htmlwidgets", + "Version": "1.6.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grDevices", + "htmltools", + "jsonlite", + "knitr", + "rmarkdown", + "yaml" + ], + "Hash": "04291cc45198225444a397606810ac37" + }, + "httpuv": { + "Package": "httpuv", + "Version": "1.6.15", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "Rcpp", + "later", + "promises", + "utils" + ], + "Hash": "d55aa087c47a63ead0f6fc10f8fa1ee0" + }, + "httr": { + "Package": "httr", + "Version": "1.4.7", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "curl", + "jsonlite", + "mime", + "openssl" + ], + "Hash": "ac107251d9d9fd72f0ca8049988f1d7f" + }, + "ica": { + "Package": "ica", + "Version": "1.0-3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d9b52ced14e24a0e69e228c20eb5eb27" + }, + "igraph": { + "Package": "igraph", + "Version": "2.0.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "cli", + "cpp11", + "grDevices", + "graphics", + "lifecycle", + "magrittr", + "methods", + "pkgconfig", + "rlang", + "stats", + "utils", + "vctrs" + ], + "Hash": "c3b7d801d722e26e4cd888e042bf9af5" + }, + "irlba": { + "Package": "irlba", + "Version": "2.3.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "methods", + "stats" + ], + "Hash": "acb06a47b732c6251afd16e19c3201ff" + }, + "isoband": { + "Package": "isoband", + "Version": "0.2.7", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grid", + "utils" + ], + "Hash": "0080607b4a1a7b28979aecef976d8bc2" + }, + "jquerylib": { + "Package": "jquerylib", + "Version": "0.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "htmltools" + ], + "Hash": "5aab57a3bd297eee1c1d862735972182" + }, + "jsonlite": { + "Package": "jsonlite", + "Version": "1.8.8", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "methods" + ], + "Hash": "e1b9c55281c5adc4dd113652d9e26768" + }, + "knitr": { + "Package": "knitr", + "Version": "1.46", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "evaluate", + "highr", + "methods", + "tools", + "xfun", + "yaml" + ], + "Hash": "6e008ab1d696a5283c79765fa7b56b47" + }, + "labeling": { + "Package": "labeling", + "Version": "0.4.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "graphics", + "stats" + ], + "Hash": "b64ec208ac5bc1852b285f665d6368b3" + }, + "later": { + "Package": "later", + "Version": "1.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Rcpp", + "rlang" + ], + "Hash": "a3e051d405326b8b0012377434c62b37" + }, + "lattice": { + "Package": "lattice", + "Version": "0.21-8", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "grid", + "stats", + "utils" + ], + "Hash": "0b8a6d63c8770f02a8b5635f3c431e6b" + }, + "lazyeval": { + "Package": "lazyeval", + "Version": "0.2.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "d908914ae53b04d4c0c0fd72ecc35370" + }, + "leiden": { + "Package": "leiden", + "Version": "0.4.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "igraph", + "methods", + "reticulate" + ], + "Hash": "b21c4ae2bb7935504c42bcdf749c04e6" + }, + "lifecycle": { + "Package": "lifecycle", + "Version": "1.0.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "glue", + "rlang" + ], + "Hash": "b8552d117e1b808b09a832f589b79035" + }, + "listenv": { + "Package": "listenv", + "Version": "0.9.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "e2fca3e12e4db979dccc6e519b10a7ee" + }, + "lme4": { + "Package": "lme4", + "Version": "1.1-35.1", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "MASS", + "Matrix", + "R", + "Rcpp", + "RcppEigen", + "boot", + "graphics", + "grid", + "lattice", + "methods", + "minqa", + "nlme", + "nloptr", + "parallel", + "splines", + "stats", + "utils" + ], + "Hash": "07fb0c5b727b15b0ce40feb641498e4e" + }, + "lmtest": { + "Package": "lmtest", + "Version": "0.9-40", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "stats", + "zoo" + ], + "Hash": "c6fafa6cccb1e1dfe7f7d122efd6e6a7" + }, + "magrittr": { + "Package": "magrittr", + "Version": "2.0.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "7ce2733a9826b3aeb1775d56fd305472" + }, + "matrixStats": { + "Package": "matrixStats", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "33a3ca9e732b57244d14f5d732ffc9eb" + }, + "memoise": { + "Package": "memoise", + "Version": "2.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "cachem", + "rlang" + ], + "Hash": "e2817ccf4a065c5d9d7f2cfbe7c1d78c" + }, + "mgcv": { + "Package": "mgcv", + "Version": "1.8-42", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "graphics", + "methods", + "nlme", + "splines", + "stats", + "utils" + ], + "Hash": "3460beba7ccc8946249ba35327ba902a" + }, + "mime": { + "Package": "mime", + "Version": "0.12", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "tools" + ], + "Hash": "18e9c28c1d3ca1560ce30658b22ce104" + }, + "miniUI": { + "Package": "miniUI", + "Version": "0.1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "htmltools", + "shiny", + "utils" + ], + "Hash": "fec5f52652d60615fdb3957b3d74324a" + }, + "minqa": { + "Package": "minqa", + "Version": "1.2.6", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "Rcpp" + ], + "Hash": "f48238f8d4740426ca12f53f27d004dd" + }, + "munsell": { + "Package": "munsell", + "Version": "0.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "colorspace", + "methods" + ], + "Hash": "4fd8900853b746af55b81fda99da7695" + }, + "nlme": { + "Package": "nlme", + "Version": "3.1-162", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "lattice", + "stats", + "utils" + ], + "Hash": "0984ce8da8da9ead8643c5cbbb60f83e" + }, + "nloptr": { + "Package": "nloptr", + "Version": "2.0.3", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "testthat" + ], + "Hash": "277c67a08f358f42b6a77826e4492f79" + }, + "nnet": { + "Package": "nnet", + "Version": "7.3-19", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats", + "utils" + ], + "Hash": "2c797b46eea7fb58ede195bc0b1f1138" + }, + "numDeriv": { + "Package": "numDeriv", + "Version": "2016.8-1.1", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R" + ], + "Hash": "df58958f293b166e4ab885ebcad90e02" + }, + "openssl": { + "Package": "openssl", + "Version": "2.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "askpass" + ], + "Hash": "2a0dc8c6adfb6f032e4d4af82d258ab5" + }, + "parallelly": { + "Package": "parallelly", + "Version": "1.37.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "parallel", + "tools", + "utils" + ], + "Hash": "5410df8d22bd36e616f2a2343dbb328c" + }, + "patchwork": { + "Package": "patchwork", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "cli", + "ggplot2", + "grDevices", + "graphics", + "grid", + "gtable", + "rlang", + "stats", + "utils" + ], + "Hash": "9c8ab14c00ac07e9e04d1664c0b74486" + }, + "pbapply": { + "Package": "pbapply", + "Version": "1.7-2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "parallel" + ], + "Hash": "68a2d681e10cf72f0afa1d84d45380e5" + }, + "pbkrtest": { + "Package": "pbkrtest", + "Version": "0.5.2", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "MASS", + "Matrix", + "R", + "broom", + "dplyr", + "lme4", + "methods", + "numDeriv", + "parallel" + ], + "Hash": "3b5b99f4d3f067bb9c1d59317d071370" + }, + "pillar": { + "Package": "pillar", + "Version": "1.9.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "cli", + "fansi", + "glue", + "lifecycle", + "rlang", + "utf8", + "utils", + "vctrs" + ], + "Hash": "15da5a8412f317beeee6175fbc76f4bb" + }, + "pkgbuild": { + "Package": "pkgbuild", + "Version": "1.4.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "callr", + "cli", + "desc", + "processx" + ], + "Hash": "a29e8e134a460a01e0ca67a4763c595b" + }, + "pkgconfig": { + "Package": "pkgconfig", + "Version": "2.0.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "utils" + ], + "Hash": "01f28d4278f15c76cddbea05899c5d6f" + }, + "pkgload": { + "Package": "pkgload", + "Version": "1.3.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "crayon", + "desc", + "fs", + "glue", + "methods", + "pkgbuild", + "rlang", + "rprojroot", + "utils", + "withr" + ], + "Hash": "876c618df5ae610be84356d5d7a5d124" + }, + "plotly": { + "Package": "plotly", + "Version": "4.10.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "RColorBrewer", + "base64enc", + "crosstalk", + "data.table", + "digest", + "dplyr", + "ggplot2", + "htmltools", + "htmlwidgets", + "httr", + "jsonlite", + "lazyeval", + "magrittr", + "promises", + "purrr", + "rlang", + "scales", + "tibble", + "tidyr", + "tools", + "vctrs", + "viridisLite" + ], + "Hash": "a1ac5c03ad5ad12b9d1597e00e23c3dd" + }, + "plyr": { + "Package": "plyr", + "Version": "1.8.9", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp" + ], + "Hash": "6b8177fd19982f0020743fadbfdbd933" + }, + "png": { + "Package": "png", + "Version": "0.1-8", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "bd54ba8a0a5faded999a7aab6e46b374" + }, + "polyclip": { + "Package": "polyclip", + "Version": "1.10-6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "436542aadb70675e361cf359285af7c7" + }, + "polynom": { + "Package": "polynom", + "Version": "1.4-1", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "graphics", + "stats" + ], + "Hash": "ceb5c2a59ba33d42d051285a3e8a5118" + }, + "praise": { + "Package": "praise", + "Version": "1.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a555924add98c99d2f411e37e7d25e9f" + }, + "processx": { + "Package": "processx", + "Version": "3.8.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "ps", + "utils" + ], + "Hash": "0c90a7d71988856bad2a2a45dd871bb9" + }, + "progressr": { + "Package": "progressr", + "Version": "0.14.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "digest", + "utils" + ], + "Hash": "ac50c4ffa8f6a46580dd4d7813add3c4" + }, + "promises": { + "Package": "promises", + "Version": "1.3.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R6", + "Rcpp", + "fastmap", + "later", + "magrittr", + "rlang", + "stats" + ], + "Hash": "434cd5388a3979e74be5c219bcd6e77d" + }, + "ps": { + "Package": "ps", + "Version": "1.7.6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "utils" + ], + "Hash": "dd2b9319ee0656c8acf45c7f40c59de7" + }, + "purrr": { + "Package": "purrr", + "Version": "1.0.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "lifecycle", + "magrittr", + "rlang", + "vctrs" + ], + "Hash": "1cba04a4e9414bdefc9dcaa99649a8dc" + }, + "quantreg": { + "Package": "quantreg", + "Version": "5.97", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "MASS", + "Matrix", + "MatrixModels", + "R", + "SparseM", + "graphics", + "methods", + "stats", + "survival" + ], + "Hash": "1bbc97f7d637ab3917c514a69047b2c1" + }, + "rappdirs": { + "Package": "rappdirs", + "Version": "0.3.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "5e3c5dc0b071b21fa128676560dbe94d" + }, + "rematch2": { + "Package": "rematch2", + "Version": "2.1.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "tibble" + ], + "Hash": "76c9e04c712a05848ae7a23d2f170a40" + }, + "renv": { + "Package": "renv", + "Version": "1.0.3", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "utils" + ], + "Hash": "41b847654f567341725473431dd0d5ab" + }, + "reshape2": { + "Package": "reshape2", + "Version": "1.4.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "plyr", + "stringr" + ], + "Hash": "bb5996d0bd962d214a11140d77589917" + }, + "reticulate": { + "Package": "reticulate", + "Version": "1.35.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "Rcpp", + "RcppTOML", + "graphics", + "here", + "jsonlite", + "methods", + "png", + "rappdirs", + "rlang", + "utils", + "withr" + ], + "Hash": "90be16b53b955990db4aa355c03d85eb" + }, + "rlang": { + "Package": "rlang", + "Version": "1.1.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "utils" + ], + "Hash": "42548638fae05fd9a9b5f3f437fbbbe2" + }, + "rmarkdown": { + "Package": "rmarkdown", + "Version": "2.26", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "bslib", + "evaluate", + "fontawesome", + "htmltools", + "jquerylib", + "jsonlite", + "knitr", + "methods", + "tinytex", + "tools", + "utils", + "xfun", + "yaml" + ], + "Hash": "9b148e7f95d33aac01f31282d49e4f44" + }, + "rprojroot": { + "Package": "rprojroot", + "Version": "2.0.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "4c8415e0ec1e29f3f4f6fc108bef0144" + }, + "rstatix": { + "Package": "rstatix", + "Version": "0.7.2", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R", + "broom", + "car", + "corrplot", + "dplyr", + "generics", + "magrittr", + "purrr", + "rlang", + "stats", + "tibble", + "tidyr", + "tidyselect", + "utils" + ], + "Hash": "5045fbb71b143878d8c51975d1d7d56d" + }, + "sass": { + "Package": "sass", + "Version": "0.4.9", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R6", + "fs", + "htmltools", + "rappdirs", + "rlang" + ], + "Hash": "d53dbfddf695303ea4ad66f86e99b95d" + }, + "scales": { + "Package": "scales", + "Version": "1.3.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "RColorBrewer", + "cli", + "farver", + "glue", + "labeling", + "lifecycle", + "munsell", + "rlang", + "viridisLite" + ], + "Hash": "c19df082ba346b0ffa6f833e92de34d1" + }, + "scattermore": { + "Package": "scattermore", + "Version": "1.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "ggplot2", + "grDevices", + "graphics", + "grid", + "scales" + ], + "Hash": "d316e4abb854dd1677f7bd3ad08bc4e8" + }, + "sctransform": { + "Package": "sctransform", + "Version": "0.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "Matrix", + "R", + "Rcpp", + "RcppArmadillo", + "dplyr", + "future", + "future.apply", + "ggplot2", + "gridExtra", + "magrittr", + "matrixStats", + "methods", + "reshape2", + "rlang" + ], + "Hash": "0242402f321be0246fb67cf8c63b3572" + }, + "shiny": { + "Package": "shiny", + "Version": "1.8.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "bslib", + "cachem", + "commonmark", + "crayon", + "fastmap", + "fontawesome", + "glue", + "grDevices", + "htmltools", + "httpuv", + "jsonlite", + "later", + "lifecycle", + "methods", + "mime", + "promises", + "rlang", + "sourcetools", + "tools", + "utils", + "withr", + "xtable" + ], + "Hash": "54b26646816af9960a4c64d8ceec75d6" + }, + "sitmo": { + "Package": "sitmo", + "Version": "2.0.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp" + ], + "Hash": "c956d93f6768a9789edbc13072b70c78" + }, + "sourcetools": { + "Package": "sourcetools", + "Version": "0.1.7-1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "5f5a7629f956619d519205ec475fe647" + }, + "sp": { + "Package": "sp", + "Version": "2.1-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "grid", + "lattice", + "methods", + "stats", + "utils" + ], + "Hash": "1a0cc0cec2915700e63fd0921085cf6a" + }, + "spam": { + "Package": "spam", + "Version": "2.10-0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp", + "dotCall64", + "grid", + "methods" + ], + "Hash": "ffe1f9e95a4375530747b268f82b5086" + }, + "spatstat.data": { + "Package": "spatstat.data", + "Version": "3.0-4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "spatstat.utils" + ], + "Hash": "114a35341cb4955e1b8d30e28ec356c6" + }, + "spatstat.explore": { + "Package": "spatstat.explore", + "Version": "3.2-7", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "abind", + "goftest", + "grDevices", + "graphics", + "methods", + "nlme", + "spatstat.data", + "spatstat.geom", + "spatstat.random", + "spatstat.sparse", + "spatstat.utils", + "stats", + "utils" + ], + "Hash": "b590c5d278d0606d53278afac666ad60" + }, + "spatstat.geom": { + "Package": "spatstat.geom", + "Version": "3.2-9", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "deldir", + "grDevices", + "graphics", + "methods", + "polyclip", + "spatstat.data", + "spatstat.utils", + "stats", + "utils" + ], + "Hash": "96b9f23da42d16660aa6d136ad99cf5f" + }, + "spatstat.random": { + "Package": "spatstat.random", + "Version": "3.2-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "methods", + "spatstat.data", + "spatstat.geom", + "spatstat.utils", + "stats", + "utils" + ], + "Hash": "01aff260c49550e820a47d97e566af6a" + }, + "spatstat.sparse": { + "Package": "spatstat.sparse", + "Version": "3.0-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "abind", + "methods", + "spatstat.utils", + "stats", + "tensor", + "utils" + ], + "Hash": "1daaecefd754bb259a5ad5ce95a2cdcc" + }, + "spatstat.utils": { + "Package": "spatstat.utils", + "Version": "3.0-4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "methods", + "stats", + "utils" + ], + "Hash": "81d929f43f8532c1f8180a105449d414" + }, + "stringi": { + "Package": "stringi", + "Version": "1.8.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats", + "tools", + "utils" + ], + "Hash": "058aebddea264f4c99401515182e656a" + }, + "stringr": { + "Package": "stringr", + "Version": "1.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "glue", + "lifecycle", + "magrittr", + "rlang", + "stringi", + "vctrs" + ], + "Hash": "960e2ae9e09656611e0b8214ad543207" + }, + "survival": { + "Package": "survival", + "Version": "3.5-5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "graphics", + "methods", + "splines", + "stats", + "utils" + ], + "Hash": "d683341b1fa2e8d817efde27d6e6d35b" + }, + "sys": { + "Package": "sys", + "Version": "3.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "3a1be13d68d47a8cd0bfd74739ca1555" + }, + "tensor": { + "Package": "tensor", + "Version": "1.5", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "25cfab6cf405c15bccf7e69ec39df090" + }, + "testthat": { + "Package": "testthat", + "Version": "3.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "R6", + "brio", + "callr", + "cli", + "desc", + "digest", + "evaluate", + "jsonlite", + "lifecycle", + "magrittr", + "methods", + "pkgload", + "praise", + "processx", + "ps", + "rlang", + "utils", + "waldo", + "withr" + ], + "Hash": "4767a686ebe986e6cb01d075b3f09729" + }, + "tibble": { + "Package": "tibble", + "Version": "3.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "fansi", + "lifecycle", + "magrittr", + "methods", + "pillar", + "pkgconfig", + "rlang", + "utils", + "vctrs" + ], + "Hash": "a84e2cc86d07289b3b6f5069df7a004c" + }, + "tidyr": { + "Package": "tidyr", + "Version": "1.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "cpp11", + "dplyr", + "glue", + "lifecycle", + "magrittr", + "purrr", + "rlang", + "stringr", + "tibble", + "tidyselect", + "utils", + "vctrs" + ], + "Hash": "915fb7ce036c22a6a33b5a8adb712eb1" + }, + "tidyselect": { + "Package": "tidyselect", + "Version": "1.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "glue", + "lifecycle", + "rlang", + "vctrs", + "withr" + ], + "Hash": "829f27b9c4919c16b593794a6344d6c0" + }, + "tinytex": { + "Package": "tinytex", + "Version": "0.50", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "xfun" + ], + "Hash": "be7a76845222ad20adb761f462eed3ea" + }, + "umap": { + "Package": "umap", + "Version": "0.2.10.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "RSpectra", + "Rcpp", + "methods", + "openssl", + "reticulate", + "stats" + ], + "Hash": "249dae2d91909b98bc93086d82cd7c7f" + }, + "utf8": { + "Package": "utf8", + "Version": "1.2.4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "62b65c52671e6665f803ff02954446e9" + }, + "uwot": { + "Package": "uwot", + "Version": "0.1.16", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "FNN", + "Matrix", + "Rcpp", + "RcppAnnoy", + "RcppProgress", + "dqrng", + "irlba", + "methods" + ], + "Hash": "252deaa1c1d6d3da6946694243781ea3" + }, + "vctrs": { + "Package": "vctrs", + "Version": "0.6.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "glue", + "lifecycle", + "rlang" + ], + "Hash": "c03fa420630029418f7e6da3667aac4a" + }, + "viridisLite": { + "Package": "viridisLite", + "Version": "0.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "c826c7c4241b6fc89ff55aaea3fa7491" + }, + "waldo": { + "Package": "waldo", + "Version": "0.5.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "cli", + "diffobj", + "fansi", + "glue", + "methods", + "rematch2", + "rlang", + "tibble" + ], + "Hash": "c7d3fd6d29ab077cbac8f0e2751449e6" + }, + "withr": { + "Package": "withr", + "Version": "3.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics" + ], + "Hash": "d31b6c62c10dcf11ec530ca6b0dd5d35" + }, + "xfun": { + "Package": "xfun", + "Version": "0.43", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grDevices", + "stats", + "tools" + ], + "Hash": "ab6371d8653ce5f2f9290f4ec7b42a8e" + }, + "xtable": { + "Package": "xtable", + "Version": "1.8-4", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats", + "utils" + ], + "Hash": "b8acdf8af494d9ec19ccb2481a9b11c2" + }, + "yaml": { + "Package": "yaml", + "Version": "2.3.8", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "29240487a071f535f5e5d5a323b7afbd" + }, + "zoo": { + "Package": "zoo", + "Version": "1.8-12", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "lattice", + "stats", + "utils" + ], + "Hash": "5c715954112b45499fb1dadc6ee6ee3e" + } + } +} diff --git a/renv/.gitignore b/renv/.gitignore new file mode 100644 index 0000000..0ec0cbb --- /dev/null +++ b/renv/.gitignore @@ -0,0 +1,7 @@ +library/ +local/ +cellar/ +lock/ +python/ +sandbox/ +staging/ diff --git a/renv/activate.R b/renv/activate.R new file mode 100644 index 0000000..cb5401f --- /dev/null +++ b/renv/activate.R @@ -0,0 +1,1180 @@ + +local({ + + # the requested version of renv + version <- "1.0.3" + attr(version, "sha") <- NULL + + # the project directory + project <- getwd() + + # use start-up diagnostics if enabled + diagnostics <- Sys.getenv("RENV_STARTUP_DIAGNOSTICS", unset = "FALSE") + if (diagnostics) { + start <- Sys.time() + profile <- tempfile("renv-startup-", fileext = ".Rprof") + utils::Rprof(profile) + on.exit({ + utils::Rprof(NULL) + elapsed <- signif(difftime(Sys.time(), start, units = "auto"), digits = 2L) + writeLines(sprintf("- renv took %s to run the autoloader.", format(elapsed))) + writeLines(sprintf("- Profile: %s", profile)) + print(utils::summaryRprof(profile)) + }, add = TRUE) + } + + # figure out whether the autoloader is enabled + enabled <- local({ + + # first, check config option + override <- getOption("renv.config.autoloader.enabled") + if (!is.null(override)) + return(override) + + # next, check environment variables + # TODO: prefer using the configuration one in the future + envvars <- c( + "RENV_CONFIG_AUTOLOADER_ENABLED", + "RENV_AUTOLOADER_ENABLED", + "RENV_ACTIVATE_PROJECT" + ) + + for (envvar in envvars) { + envval <- Sys.getenv(envvar, unset = NA) + if (!is.na(envval)) + return(tolower(envval) %in% c("true", "t", "1")) + } + + # enable by default + TRUE + + }) + + if (!enabled) + return(FALSE) + + # avoid recursion + if (identical(getOption("renv.autoloader.running"), TRUE)) { + warning("ignoring recursive attempt to run renv autoloader") + return(invisible(TRUE)) + } + + # signal that we're loading renv during R startup + options(renv.autoloader.running = TRUE) + on.exit(options(renv.autoloader.running = NULL), add = TRUE) + + # signal that we've consented to use renv + options(renv.consent = TRUE) + + # load the 'utils' package eagerly -- this ensures that renv shims, which + # mask 'utils' packages, will come first on the search path + library(utils, lib.loc = .Library) + + # unload renv if it's already been loaded + if ("renv" %in% loadedNamespaces()) + unloadNamespace("renv") + + # load bootstrap tools + `%||%` <- function(x, y) { + if (is.null(x)) y else x + } + + catf <- function(fmt, ..., appendLF = TRUE) { + + quiet <- getOption("renv.bootstrap.quiet", default = FALSE) + if (quiet) + return(invisible()) + + msg <- sprintf(fmt, ...) + cat(msg, file = stdout(), sep = if (appendLF) "\n" else "") + + invisible(msg) + + } + + header <- function(label, + ..., + prefix = "#", + suffix = "-", + n = min(getOption("width"), 78)) + { + label <- sprintf(label, ...) + n <- max(n - nchar(label) - nchar(prefix) - 2L, 8L) + if (n <= 0) + return(paste(prefix, label)) + + tail <- paste(rep.int(suffix, n), collapse = "") + paste0(prefix, " ", label, " ", tail) + + } + + startswith <- function(string, prefix) { + substring(string, 1, nchar(prefix)) == prefix + } + + bootstrap <- function(version, library) { + + friendly <- renv_bootstrap_version_friendly(version) + section <- header(sprintf("Bootstrapping renv %s", friendly)) + catf(section) + + # attempt to download renv + catf("- Downloading renv ... ", appendLF = FALSE) + withCallingHandlers( + tarball <- renv_bootstrap_download(version), + error = function(err) { + catf("FAILED") + stop("failed to download:\n", conditionMessage(err)) + } + ) + catf("OK") + on.exit(unlink(tarball), add = TRUE) + + # now attempt to install + catf("- Installing renv ... ", appendLF = FALSE) + withCallingHandlers( + status <- renv_bootstrap_install(version, tarball, library), + error = function(err) { + catf("FAILED") + stop("failed to install:\n", conditionMessage(err)) + } + ) + catf("OK") + + # add empty line to break up bootstrapping from normal output + catf("") + + return(invisible()) + } + + renv_bootstrap_tests_running <- function() { + getOption("renv.tests.running", default = FALSE) + } + + renv_bootstrap_repos <- function() { + + # get CRAN repository + cran <- getOption("renv.repos.cran", "https://cloud.r-project.org") + + # check for repos override + repos <- Sys.getenv("RENV_CONFIG_REPOS_OVERRIDE", unset = NA) + if (!is.na(repos)) { + + # check for RSPM; if set, use a fallback repository for renv + rspm <- Sys.getenv("RSPM", unset = NA) + if (identical(rspm, repos)) + repos <- c(RSPM = rspm, CRAN = cran) + + return(repos) + + } + + # check for lockfile repositories + repos <- tryCatch(renv_bootstrap_repos_lockfile(), error = identity) + if (!inherits(repos, "error") && length(repos)) + return(repos) + + # retrieve current repos + repos <- getOption("repos") + + # ensure @CRAN@ entries are resolved + repos[repos == "@CRAN@"] <- cran + + # add in renv.bootstrap.repos if set + default <- c(FALLBACK = "https://cloud.r-project.org") + extra <- getOption("renv.bootstrap.repos", default = default) + repos <- c(repos, extra) + + # remove duplicates that might've snuck in + dupes <- duplicated(repos) | duplicated(names(repos)) + repos[!dupes] + + } + + renv_bootstrap_repos_lockfile <- function() { + + lockpath <- Sys.getenv("RENV_PATHS_LOCKFILE", unset = "renv.lock") + if (!file.exists(lockpath)) + return(NULL) + + lockfile <- tryCatch(renv_json_read(lockpath), error = identity) + if (inherits(lockfile, "error")) { + warning(lockfile) + return(NULL) + } + + repos <- lockfile$R$Repositories + if (length(repos) == 0) + return(NULL) + + keys <- vapply(repos, `[[`, "Name", FUN.VALUE = character(1)) + vals <- vapply(repos, `[[`, "URL", FUN.VALUE = character(1)) + names(vals) <- keys + + return(vals) + + } + + renv_bootstrap_download <- function(version) { + + sha <- attr(version, "sha", exact = TRUE) + + methods <- if (!is.null(sha)) { + + # attempting to bootstrap a development version of renv + c( + function() renv_bootstrap_download_tarball(sha), + function() renv_bootstrap_download_github(sha) + ) + + } else { + + # attempting to bootstrap a release version of renv + c( + function() renv_bootstrap_download_tarball(version), + function() renv_bootstrap_download_cran_latest(version), + function() renv_bootstrap_download_cran_archive(version) + ) + + } + + for (method in methods) { + path <- tryCatch(method(), error = identity) + if (is.character(path) && file.exists(path)) + return(path) + } + + stop("All download methods failed") + + } + + renv_bootstrap_download_impl <- function(url, destfile) { + + mode <- "wb" + + # https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17715 + fixup <- + Sys.info()[["sysname"]] == "Windows" && + substring(url, 1L, 5L) == "file:" + + if (fixup) + mode <- "w+b" + + args <- list( + url = url, + destfile = destfile, + mode = mode, + quiet = TRUE + ) + + if ("headers" %in% names(formals(utils::download.file))) + args$headers <- renv_bootstrap_download_custom_headers(url) + + do.call(utils::download.file, args) + + } + + renv_bootstrap_download_custom_headers <- function(url) { + + headers <- getOption("renv.download.headers") + if (is.null(headers)) + return(character()) + + if (!is.function(headers)) + stopf("'renv.download.headers' is not a function") + + headers <- headers(url) + if (length(headers) == 0L) + return(character()) + + if (is.list(headers)) + headers <- unlist(headers, recursive = FALSE, use.names = TRUE) + + ok <- + is.character(headers) && + is.character(names(headers)) && + all(nzchar(names(headers))) + + if (!ok) + stop("invocation of 'renv.download.headers' did not return a named character vector") + + headers + + } + + renv_bootstrap_download_cran_latest <- function(version) { + + spec <- renv_bootstrap_download_cran_latest_find(version) + type <- spec$type + repos <- spec$repos + + baseurl <- utils::contrib.url(repos = repos, type = type) + ext <- if (identical(type, "source")) + ".tar.gz" + else if (Sys.info()[["sysname"]] == "Windows") + ".zip" + else + ".tgz" + name <- sprintf("renv_%s%s", version, ext) + url <- paste(baseurl, name, sep = "/") + + destfile <- file.path(tempdir(), name) + status <- tryCatch( + renv_bootstrap_download_impl(url, destfile), + condition = identity + ) + + if (inherits(status, "condition")) + return(FALSE) + + # report success and return + destfile + + } + + renv_bootstrap_download_cran_latest_find <- function(version) { + + # check whether binaries are supported on this system + binary <- + getOption("renv.bootstrap.binary", default = TRUE) && + !identical(.Platform$pkgType, "source") && + !identical(getOption("pkgType"), "source") && + Sys.info()[["sysname"]] %in% c("Darwin", "Windows") + + types <- c(if (binary) "binary", "source") + + # iterate over types + repositories + for (type in types) { + for (repos in renv_bootstrap_repos()) { + + # retrieve package database + db <- tryCatch( + as.data.frame( + utils::available.packages(type = type, repos = repos), + stringsAsFactors = FALSE + ), + error = identity + ) + + if (inherits(db, "error")) + next + + # check for compatible entry + entry <- db[db$Package %in% "renv" & db$Version %in% version, ] + if (nrow(entry) == 0) + next + + # found it; return spec to caller + spec <- list(entry = entry, type = type, repos = repos) + return(spec) + + } + } + + # if we got here, we failed to find renv + fmt <- "renv %s is not available from your declared package repositories" + stop(sprintf(fmt, version)) + + } + + renv_bootstrap_download_cran_archive <- function(version) { + + name <- sprintf("renv_%s.tar.gz", version) + repos <- renv_bootstrap_repos() + urls <- file.path(repos, "src/contrib/Archive/renv", name) + destfile <- file.path(tempdir(), name) + + for (url in urls) { + + status <- tryCatch( + renv_bootstrap_download_impl(url, destfile), + condition = identity + ) + + if (identical(status, 0L)) + return(destfile) + + } + + return(FALSE) + + } + + renv_bootstrap_download_tarball <- function(version) { + + # if the user has provided the path to a tarball via + # an environment variable, then use it + tarball <- Sys.getenv("RENV_BOOTSTRAP_TARBALL", unset = NA) + if (is.na(tarball)) + return() + + # allow directories + if (dir.exists(tarball)) { + name <- sprintf("renv_%s.tar.gz", version) + tarball <- file.path(tarball, name) + } + + # bail if it doesn't exist + if (!file.exists(tarball)) { + + # let the user know we weren't able to honour their request + fmt <- "- RENV_BOOTSTRAP_TARBALL is set (%s) but does not exist." + msg <- sprintf(fmt, tarball) + warning(msg) + + # bail + return() + + } + + catf("- Using local tarball '%s'.", tarball) + tarball + + } + + renv_bootstrap_download_github <- function(version) { + + enabled <- Sys.getenv("RENV_BOOTSTRAP_FROM_GITHUB", unset = "TRUE") + if (!identical(enabled, "TRUE")) + return(FALSE) + + # prepare download options + pat <- Sys.getenv("GITHUB_PAT") + if (nzchar(Sys.which("curl")) && nzchar(pat)) { + fmt <- "--location --fail --header \"Authorization: token %s\"" + extra <- sprintf(fmt, pat) + saved <- options("download.file.method", "download.file.extra") + options(download.file.method = "curl", download.file.extra = extra) + on.exit(do.call(base::options, saved), add = TRUE) + } else if (nzchar(Sys.which("wget")) && nzchar(pat)) { + fmt <- "--header=\"Authorization: token %s\"" + extra <- sprintf(fmt, pat) + saved <- options("download.file.method", "download.file.extra") + options(download.file.method = "wget", download.file.extra = extra) + on.exit(do.call(base::options, saved), add = TRUE) + } + + url <- file.path("https://api.github.com/repos/rstudio/renv/tarball", version) + name <- sprintf("renv_%s.tar.gz", version) + destfile <- file.path(tempdir(), name) + + status <- tryCatch( + renv_bootstrap_download_impl(url, destfile), + condition = identity + ) + + if (!identical(status, 0L)) + return(FALSE) + + renv_bootstrap_download_augment(destfile) + + return(destfile) + + } + + # Add Sha to DESCRIPTION. This is stop gap until #890, after which we + # can use renv::install() to fully capture metadata. + renv_bootstrap_download_augment <- function(destfile) { + sha <- renv_bootstrap_git_extract_sha1_tar(destfile) + if (is.null(sha)) { + return() + } + + # Untar + tempdir <- tempfile("renv-github-") + on.exit(unlink(tempdir, recursive = TRUE), add = TRUE) + untar(destfile, exdir = tempdir) + pkgdir <- dir(tempdir, full.names = TRUE)[[1]] + + # Modify description + desc_path <- file.path(pkgdir, "DESCRIPTION") + desc_lines <- readLines(desc_path) + remotes_fields <- c( + "RemoteType: github", + "RemoteHost: api.github.com", + "RemoteRepo: renv", + "RemoteUsername: rstudio", + "RemotePkgRef: rstudio/renv", + paste("RemoteRef: ", sha), + paste("RemoteSha: ", sha) + ) + writeLines(c(desc_lines[desc_lines != ""], remotes_fields), con = desc_path) + + # Re-tar + local({ + old <- setwd(tempdir) + on.exit(setwd(old), add = TRUE) + + tar(destfile, compression = "gzip") + }) + invisible() + } + + # Extract the commit hash from a git archive. Git archives include the SHA1 + # hash as the comment field of the tarball pax extended header + # (see https://www.kernel.org/pub/software/scm/git/docs/git-archive.html) + # For GitHub archives this should be the first header after the default one + # (512 byte) header. + renv_bootstrap_git_extract_sha1_tar <- function(bundle) { + + # open the bundle for reading + # We use gzcon for everything because (from ?gzcon) + # > Reading from a connection which does not supply a 'gzip' magic + # > header is equivalent to reading from the original connection + conn <- gzcon(file(bundle, open = "rb", raw = TRUE)) + on.exit(close(conn)) + + # The default pax header is 512 bytes long and the first pax extended header + # with the comment should be 51 bytes long + # `52 comment=` (11 chars) + 40 byte SHA1 hash + len <- 0x200 + 0x33 + res <- rawToChar(readBin(conn, "raw", n = len)[0x201:len]) + + if (grepl("^52 comment=", res)) { + sub("52 comment=", "", res) + } else { + NULL + } + } + + renv_bootstrap_install <- function(version, tarball, library) { + + # attempt to install it into project library + dir.create(library, showWarnings = FALSE, recursive = TRUE) + output <- renv_bootstrap_install_impl(library, tarball) + + # check for successful install + status <- attr(output, "status") + if (is.null(status) || identical(status, 0L)) + return(status) + + # an error occurred; report it + header <- "installation of renv failed" + lines <- paste(rep.int("=", nchar(header)), collapse = "") + text <- paste(c(header, lines, output), collapse = "\n") + stop(text) + + } + + renv_bootstrap_install_impl <- function(library, tarball) { + + # invoke using system2 so we can capture and report output + bin <- R.home("bin") + exe <- if (Sys.info()[["sysname"]] == "Windows") "R.exe" else "R" + R <- file.path(bin, exe) + + args <- c( + "--vanilla", "CMD", "INSTALL", "--no-multiarch", + "-l", shQuote(path.expand(library)), + shQuote(path.expand(tarball)) + ) + + system2(R, args, stdout = TRUE, stderr = TRUE) + + } + + renv_bootstrap_platform_prefix <- function() { + + # construct version prefix + version <- paste(R.version$major, R.version$minor, sep = ".") + prefix <- paste("R", numeric_version(version)[1, 1:2], sep = "-") + + # include SVN revision for development versions of R + # (to avoid sharing platform-specific artefacts with released versions of R) + devel <- + identical(R.version[["status"]], "Under development (unstable)") || + identical(R.version[["nickname"]], "Unsuffered Consequences") + + if (devel) + prefix <- paste(prefix, R.version[["svn rev"]], sep = "-r") + + # build list of path components + components <- c(prefix, R.version$platform) + + # include prefix if provided by user + prefix <- renv_bootstrap_platform_prefix_impl() + if (!is.na(prefix) && nzchar(prefix)) + components <- c(prefix, components) + + # build prefix + paste(components, collapse = "/") + + } + + renv_bootstrap_platform_prefix_impl <- function() { + + # if an explicit prefix has been supplied, use it + prefix <- Sys.getenv("RENV_PATHS_PREFIX", unset = NA) + if (!is.na(prefix)) + return(prefix) + + # if the user has requested an automatic prefix, generate it + auto <- Sys.getenv("RENV_PATHS_PREFIX_AUTO", unset = NA) + if (auto %in% c("TRUE", "True", "true", "1")) + return(renv_bootstrap_platform_prefix_auto()) + + # empty string on failure + "" + + } + + renv_bootstrap_platform_prefix_auto <- function() { + + prefix <- tryCatch(renv_bootstrap_platform_os(), error = identity) + if (inherits(prefix, "error") || prefix %in% "unknown") { + + msg <- paste( + "failed to infer current operating system", + "please file a bug report at https://github.com/rstudio/renv/issues", + sep = "; " + ) + + warning(msg) + + } + + prefix + + } + + renv_bootstrap_platform_os <- function() { + + sysinfo <- Sys.info() + sysname <- sysinfo[["sysname"]] + + # handle Windows + macOS up front + if (sysname == "Windows") + return("windows") + else if (sysname == "Darwin") + return("macos") + + # check for os-release files + for (file in c("/etc/os-release", "/usr/lib/os-release")) + if (file.exists(file)) + return(renv_bootstrap_platform_os_via_os_release(file, sysinfo)) + + # check for redhat-release files + if (file.exists("/etc/redhat-release")) + return(renv_bootstrap_platform_os_via_redhat_release()) + + "unknown" + + } + + renv_bootstrap_platform_os_via_os_release <- function(file, sysinfo) { + + # read /etc/os-release + release <- utils::read.table( + file = file, + sep = "=", + quote = c("\"", "'"), + col.names = c("Key", "Value"), + comment.char = "#", + stringsAsFactors = FALSE + ) + + vars <- as.list(release$Value) + names(vars) <- release$Key + + # get os name + os <- tolower(sysinfo[["sysname"]]) + + # read id + id <- "unknown" + for (field in c("ID", "ID_LIKE")) { + if (field %in% names(vars) && nzchar(vars[[field]])) { + id <- vars[[field]] + break + } + } + + # read version + version <- "unknown" + for (field in c("UBUNTU_CODENAME", "VERSION_CODENAME", "VERSION_ID", "BUILD_ID")) { + if (field %in% names(vars) && nzchar(vars[[field]])) { + version <- vars[[field]] + break + } + } + + # join together + paste(c(os, id, version), collapse = "-") + + } + + renv_bootstrap_platform_os_via_redhat_release <- function() { + + # read /etc/redhat-release + contents <- readLines("/etc/redhat-release", warn = FALSE) + + # infer id + id <- if (grepl("centos", contents, ignore.case = TRUE)) + "centos" + else if (grepl("redhat", contents, ignore.case = TRUE)) + "redhat" + else + "unknown" + + # try to find a version component (very hacky) + version <- "unknown" + + parts <- strsplit(contents, "[[:space:]]")[[1L]] + for (part in parts) { + + nv <- tryCatch(numeric_version(part), error = identity) + if (inherits(nv, "error")) + next + + version <- nv[1, 1] + break + + } + + paste(c("linux", id, version), collapse = "-") + + } + + renv_bootstrap_library_root_name <- function(project) { + + # use project name as-is if requested + asis <- Sys.getenv("RENV_PATHS_LIBRARY_ROOT_ASIS", unset = "FALSE") + if (asis) + return(basename(project)) + + # otherwise, disambiguate based on project's path + id <- substring(renv_bootstrap_hash_text(project), 1L, 8L) + paste(basename(project), id, sep = "-") + + } + + renv_bootstrap_library_root <- function(project) { + + prefix <- renv_bootstrap_profile_prefix() + + path <- Sys.getenv("RENV_PATHS_LIBRARY", unset = NA) + if (!is.na(path)) + return(paste(c(path, prefix), collapse = "/")) + + path <- renv_bootstrap_library_root_impl(project) + if (!is.null(path)) { + name <- renv_bootstrap_library_root_name(project) + return(paste(c(path, prefix, name), collapse = "/")) + } + + renv_bootstrap_paths_renv("library", project = project) + + } + + renv_bootstrap_library_root_impl <- function(project) { + + root <- Sys.getenv("RENV_PATHS_LIBRARY_ROOT", unset = NA) + if (!is.na(root)) + return(root) + + type <- renv_bootstrap_project_type(project) + if (identical(type, "package")) { + userdir <- renv_bootstrap_user_dir() + return(file.path(userdir, "library")) + } + + } + + renv_bootstrap_validate_version <- function(version, description = NULL) { + + # resolve description file + # + # avoid passing lib.loc to `packageDescription()` below, since R will + # use the loaded version of the package by default anyhow. note that + # this function should only be called after 'renv' is loaded + # https://github.com/rstudio/renv/issues/1625 + description <- description %||% packageDescription("renv") + + # check whether requested version 'version' matches loaded version of renv + sha <- attr(version, "sha", exact = TRUE) + valid <- if (!is.null(sha)) + renv_bootstrap_validate_version_dev(sha, description) + else + renv_bootstrap_validate_version_release(version, description) + + if (valid) + return(TRUE) + + # the loaded version of renv doesn't match the requested version; + # give the user instructions on how to proceed + remote <- if (!is.null(description[["RemoteSha"]])) { + paste("rstudio/renv", description[["RemoteSha"]], sep = "@") + } else { + paste("renv", description[["Version"]], sep = "@") + } + + # display both loaded version + sha if available + friendly <- renv_bootstrap_version_friendly( + version = description[["Version"]], + sha = description[["RemoteSha"]] + ) + + fmt <- paste( + "renv %1$s was loaded from project library, but this project is configured to use renv %2$s.", + "- Use `renv::record(\"%3$s\")` to record renv %1$s in the lockfile.", + "- Use `renv::restore(packages = \"renv\")` to install renv %2$s into the project library.", + sep = "\n" + ) + catf(fmt, friendly, renv_bootstrap_version_friendly(version), remote) + + FALSE + + } + + renv_bootstrap_validate_version_dev <- function(version, description) { + expected <- description[["RemoteSha"]] + is.character(expected) && startswith(expected, version) + } + + renv_bootstrap_validate_version_release <- function(version, description) { + expected <- description[["Version"]] + is.character(expected) && identical(expected, version) + } + + renv_bootstrap_hash_text <- function(text) { + + hashfile <- tempfile("renv-hash-") + on.exit(unlink(hashfile), add = TRUE) + + writeLines(text, con = hashfile) + tools::md5sum(hashfile) + + } + + renv_bootstrap_load <- function(project, libpath, version) { + + # try to load renv from the project library + if (!requireNamespace("renv", lib.loc = libpath, quietly = TRUE)) + return(FALSE) + + # warn if the version of renv loaded does not match + renv_bootstrap_validate_version(version) + + # execute renv load hooks, if any + hooks <- getHook("renv::autoload") + for (hook in hooks) + if (is.function(hook)) + tryCatch(hook(), error = warnify) + + # load the project + renv::load(project) + + TRUE + + } + + renv_bootstrap_profile_load <- function(project) { + + # if RENV_PROFILE is already set, just use that + profile <- Sys.getenv("RENV_PROFILE", unset = NA) + if (!is.na(profile) && nzchar(profile)) + return(profile) + + # check for a profile file (nothing to do if it doesn't exist) + path <- renv_bootstrap_paths_renv("profile", profile = FALSE, project = project) + if (!file.exists(path)) + return(NULL) + + # read the profile, and set it if it exists + contents <- readLines(path, warn = FALSE) + if (length(contents) == 0L) + return(NULL) + + # set RENV_PROFILE + profile <- contents[[1L]] + if (!profile %in% c("", "default")) + Sys.setenv(RENV_PROFILE = profile) + + profile + + } + + renv_bootstrap_profile_prefix <- function() { + profile <- renv_bootstrap_profile_get() + if (!is.null(profile)) + return(file.path("profiles", profile, "renv")) + } + + renv_bootstrap_profile_get <- function() { + profile <- Sys.getenv("RENV_PROFILE", unset = "") + renv_bootstrap_profile_normalize(profile) + } + + renv_bootstrap_profile_set <- function(profile) { + profile <- renv_bootstrap_profile_normalize(profile) + if (is.null(profile)) + Sys.unsetenv("RENV_PROFILE") + else + Sys.setenv(RENV_PROFILE = profile) + } + + renv_bootstrap_profile_normalize <- function(profile) { + + if (is.null(profile) || profile %in% c("", "default")) + return(NULL) + + profile + + } + + renv_bootstrap_path_absolute <- function(path) { + + substr(path, 1L, 1L) %in% c("~", "/", "\\") || ( + substr(path, 1L, 1L) %in% c(letters, LETTERS) && + substr(path, 2L, 3L) %in% c(":/", ":\\") + ) + + } + + renv_bootstrap_paths_renv <- function(..., profile = TRUE, project = NULL) { + renv <- Sys.getenv("RENV_PATHS_RENV", unset = "renv") + root <- if (renv_bootstrap_path_absolute(renv)) NULL else project + prefix <- if (profile) renv_bootstrap_profile_prefix() + components <- c(root, renv, prefix, ...) + paste(components, collapse = "/") + } + + renv_bootstrap_project_type <- function(path) { + + descpath <- file.path(path, "DESCRIPTION") + if (!file.exists(descpath)) + return("unknown") + + desc <- tryCatch( + read.dcf(descpath, all = TRUE), + error = identity + ) + + if (inherits(desc, "error")) + return("unknown") + + type <- desc$Type + if (!is.null(type)) + return(tolower(type)) + + package <- desc$Package + if (!is.null(package)) + return("package") + + "unknown" + + } + + renv_bootstrap_user_dir <- function() { + dir <- renv_bootstrap_user_dir_impl() + path.expand(chartr("\\", "/", dir)) + } + + renv_bootstrap_user_dir_impl <- function() { + + # use local override if set + override <- getOption("renv.userdir.override") + if (!is.null(override)) + return(override) + + # use R_user_dir if available + tools <- asNamespace("tools") + if (is.function(tools$R_user_dir)) + return(tools$R_user_dir("renv", "cache")) + + # try using our own backfill for older versions of R + envvars <- c("R_USER_CACHE_DIR", "XDG_CACHE_HOME") + for (envvar in envvars) { + root <- Sys.getenv(envvar, unset = NA) + if (!is.na(root)) + return(file.path(root, "R/renv")) + } + + # use platform-specific default fallbacks + if (Sys.info()[["sysname"]] == "Windows") + file.path(Sys.getenv("LOCALAPPDATA"), "R/cache/R/renv") + else if (Sys.info()[["sysname"]] == "Darwin") + "~/Library/Caches/org.R-project.R/R/renv" + else + "~/.cache/R/renv" + + } + + renv_bootstrap_version_friendly <- function(version, shafmt = NULL, sha = NULL) { + sha <- sha %||% attr(version, "sha", exact = TRUE) + parts <- c(version, sprintf(shafmt %||% " [sha: %s]", substring(sha, 1L, 7L))) + paste(parts, collapse = "") + } + + renv_bootstrap_exec <- function(project, libpath, version) { + if (!renv_bootstrap_load(project, libpath, version)) + renv_bootstrap_run(version, libpath) + } + + renv_bootstrap_run <- function(version, libpath) { + + # perform bootstrap + bootstrap(version, libpath) + + # exit early if we're just testing bootstrap + if (!is.na(Sys.getenv("RENV_BOOTSTRAP_INSTALL_ONLY", unset = NA))) + return(TRUE) + + # try again to load + if (requireNamespace("renv", lib.loc = libpath, quietly = TRUE)) { + return(renv::load(project = getwd())) + } + + # failed to download or load renv; warn the user + msg <- c( + "Failed to find an renv installation: the project will not be loaded.", + "Use `renv::activate()` to re-initialize the project." + ) + + warning(paste(msg, collapse = "\n"), call. = FALSE) + + } + + renv_json_read <- function(file = NULL, text = NULL) { + + jlerr <- NULL + + # if jsonlite is loaded, use that instead + if ("jsonlite" %in% loadedNamespaces()) { + + json <- catch(renv_json_read_jsonlite(file, text)) + if (!inherits(json, "error")) + return(json) + + jlerr <- json + + } + + # otherwise, fall back to the default JSON reader + json <- catch(renv_json_read_default(file, text)) + if (!inherits(json, "error")) + return(json) + + # report an error + if (!is.null(jlerr)) + stop(jlerr) + else + stop(json) + + } + + renv_json_read_jsonlite <- function(file = NULL, text = NULL) { + text <- paste(text %||% read(file), collapse = "\n") + jsonlite::fromJSON(txt = text, simplifyVector = FALSE) + } + + renv_json_read_default <- function(file = NULL, text = NULL) { + + # find strings in the JSON + text <- paste(text %||% read(file), collapse = "\n") + pattern <- '["](?:(?:\\\\.)|(?:[^"\\\\]))*?["]' + locs <- gregexpr(pattern, text, perl = TRUE)[[1]] + + # if any are found, replace them with placeholders + replaced <- text + strings <- character() + replacements <- character() + + if (!identical(c(locs), -1L)) { + + # get the string values + starts <- locs + ends <- locs + attr(locs, "match.length") - 1L + strings <- substring(text, starts, ends) + + # only keep those requiring escaping + strings <- grep("[[\\]{}:]", strings, perl = TRUE, value = TRUE) + + # compute replacements + replacements <- sprintf('"\032%i\032"', seq_along(strings)) + + # replace the strings + mapply(function(string, replacement) { + replaced <<- sub(string, replacement, replaced, fixed = TRUE) + }, strings, replacements) + + } + + # transform the JSON into something the R parser understands + transformed <- replaced + transformed <- gsub("{}", "`names<-`(list(), character())", transformed, fixed = TRUE) + transformed <- gsub("[[{]", "list(", transformed, perl = TRUE) + transformed <- gsub("[]}]", ")", transformed, perl = TRUE) + transformed <- gsub(":", "=", transformed, fixed = TRUE) + text <- paste(transformed, collapse = "\n") + + # parse it + json <- parse(text = text, keep.source = FALSE, srcfile = NULL)[[1L]] + + # construct map between source strings, replaced strings + map <- as.character(parse(text = strings)) + names(map) <- as.character(parse(text = replacements)) + + # convert to list + map <- as.list(map) + + # remap strings in object + remapped <- renv_json_remap(json, map) + + # evaluate + eval(remapped, envir = baseenv()) + + } + + renv_json_remap <- function(json, map) { + + # fix names + if (!is.null(names(json))) { + lhs <- match(names(json), names(map), nomatch = 0L) + rhs <- match(names(map), names(json), nomatch = 0L) + names(json)[rhs] <- map[lhs] + } + + # fix values + if (is.character(json)) + return(map[[json]] %||% json) + + # handle true, false, null + if (is.name(json)) { + text <- as.character(json) + if (text == "true") + return(TRUE) + else if (text == "false") + return(FALSE) + else if (text == "null") + return(NULL) + } + + # recurse + if (is.recursive(json)) { + for (i in seq_along(json)) { + json[i] <- list(renv_json_remap(json[[i]], map)) + } + } + + json + + } + + # load the renv profile, if any + renv_bootstrap_profile_load(project) + + # construct path to library root + root <- renv_bootstrap_library_root(project) + + # construct library prefix for platform + prefix <- renv_bootstrap_platform_prefix() + + # construct full libpath + libpath <- file.path(root, prefix) + + # run bootstrap code + renv_bootstrap_exec(project, libpath, version) + + invisible() + +}) diff --git a/renv/settings.json b/renv/settings.json new file mode 100644 index 0000000..ffdbb32 --- /dev/null +++ b/renv/settings.json @@ -0,0 +1,19 @@ +{ + "bioconductor.version": null, + "external.libraries": [], + "ignored.packages": [], + "package.dependency.fields": [ + "Imports", + "Depends", + "LinkingTo" + ], + "ppm.enabled": null, + "ppm.ignored.urls": [], + "r.version": null, + "snapshot.type": "implicit", + "use.cache": true, + "vcs.ignore.cellar": true, + "vcs.ignore.library": true, + "vcs.ignore.local": true, + "vcs.manage.ignores": true +} diff --git a/scripts/spatial_RNAseq-clustering.R b/scripts/spatial_RNAseq-clustering.R new file mode 100644 index 0000000..fc8bfa8 --- /dev/null +++ b/scripts/spatial_RNAseq-clustering.R @@ -0,0 +1,763 @@ +# spatial RNAseq clustering +# Interfaculty Bioinformatics Unit (IBU) (https://www.bioinformatics.unibe.ch + +library(Seurat) #v4.1 +library(ggplot2) +library(patchwork) +library(dplyr) +library(ggpubr) +library(clusterProfiler) +library(clustree) +library(future) +library(openxlsx) +library(RColorBrewer) + + + + +#library(STutility) +#library(topGO) +#library(limma) +#library(biomaRt) +#library(msigdbr) +#library(scales) + + +# setup -------------------------- + +# Specify input directories and files +workPath <- "/data/projects/p907_SIBdays_tutorial_-_spatial_transcriptomics/" ### ADJUST + +analysisFolder <- paste0(workPath, "scripts/analysis") + +name <- "Posterior" ### ADJUST +#--------------------------------- + + +seu.norm <- readRDS(paste0(analysisFolder, "/normalized.rds")) + +#parallelization +plan("multisession", workers = availableCores()) +options(future.globals.maxSize = 8000 * 1024^2) + + +seu.norm <- RunPCA(seu.norm, assay = "SCT", npcs = 50, verbose = FALSE) +seu.norm <- RunUMAP(seu.norm, reduction = "pca", dims = 1:50) + +DimPlot(seu.norm, reduction = "pca", group.by = "orig.ident") + NoLegend() +DimPlot(seu.norm, reduction = "umap", group.by = "orig.ident")+ ggtitle("") + NoLegend() + +#seu.norm <- ScaleData(seu.norm, features=rownames(seu.norm), verbose = FALSE) + + +# 1. Dimensionality reduction ######################### + Tool: [Seurat](https://satijalab.org/seurat/articles/spatial_vignette.html) + +# Selecting which PCs to use: --------- +# To overcome the extensive technical noise in any single gene, Seurat clusters spots based +# on their PCA scores, with each PC essentially representing a metagene that combines information +# across a correlated group of genes. + +# ElbowPlot ranks the principal components based on the variance explained by each. This plot typically +# shows an "elbow", which can be used to assess how many PCs are needed to capture most of the signal in the data. + +use.pcs <- 50 +ElbowPlot(seu.norm, ndims=use.pcs) + +# The sctransform workflow performs more effective normalization, strongly removing technical effects from the data, +# this means that higher PCs are more likely to represent subtle, but biologically relevant, sources of heterogeneity +# - so including them may improve downstream analysis. Therefore, higher number of PC can be used. +# By default we are using the first 50 PCs + + + +# 2. Identifying clusters ######################### +# Seurat implements a graph-based clustering approach. Distances between the spots are calculated based on +# previously identified PCs. Briefly, Seurat identifies clusters of spots by a shared nearest neighbor (SNN) +# modularity optimization based clustering algorithm. First, it identifies k-nearest neighbors (KNN) and constructs +# the SNN graph. Then it optimizes the modularity function to determine clusters. For a full description of the +# algorithms, see Waltman and van Eck (2013) The European Physical Journal B. + +# The FindClusters function implements the procedure, and contains a resolution parameter that sets the granularity +# of the downstream clustering, with increased values leading to a greater number of clusters. + +# Selecting which resolution to use: ------- +resolution_vector <- seq(0.2,2,0.2) +seu.norm <- FindNeighbors(seu.norm, reduction="pca", dims=1:use.pcs) +seu.norm <- FindClusters(object=seu.norm, resolution=resolution_vector, verbose=FALSE) + +resTable <- as.data.frame(sapply(grep("res",colnames(seu.norm@meta.data),value = TRUE), function(x) length(unique(seu.norm@meta.data[,x])))) +colnames(resTable) <- "number of clusters" + +# how many clusters in each resolution +resTable + +for(i in seq(resolution_vector)){ + print(DimPlot(seu.norm, reduction = "umap", label=TRUE, group.by=paste0("SCT_snn_res.", resolution_vector[i])) + labs(title=paste0("res_", resolution_vector[i]))) + print(SpatialDimPlot(seu.norm, label = TRUE, combine=FALSE, label.size = 3, group.by=paste0("SCT_snn_res.", resolution_vector[i]))) +} + +# Plotting clustering trees +# To build a clustering tree we need to look at how spots move as the clustering resolution is increased. +# Each cluster forms a node in the tree and edges are constructed by considering the spots in a cluster +# at a lower resolution that end up in a cluster at the next higher resolution. By connecting clusters +# in this way we can see how clusters are related to each other, which are clearly distinct and which are unstable. + +# If there are nodes with multiple incoming edges, it is a good indication that we have over-clustered the data. + +clustree(seu.norm, prefis="res.") + + + +# -> Set the resolution to 0.4 ############ +resolution <- "0.4" + +# Final clustering: +Idents(seu.norm) <- paste0("SCT_snn_res.", resolution) +seu.norm$finalCluster <- Idents(seu.norm) + +DimPlot(seu.norm, reduction = "umap", label=TRUE, group.by = "ident") +SpatialDimPlot(seu.norm, label = TRUE, combine=FALSE, label.size = 3) + + +# As there are many colors, it can be challenging to visualize which spot belongs to which cluster. +# Therefore, we highlight spots of each cluster seperately: +SpatialDimPlot(seu.norm, combine=TRUE, images=name, cells.highlight = CellsByIdentities(object = seu.norm), cols.highlight = c("red", "NA"), facet.highlight = TRUE, ncol = 4) + + +# Total number of cells per cluster per slice +seu.norm@meta.data %>% dplyr::select(orig.ident, finalCluster) %>% table() + + + +# 3. Identify marker genes for each cluster ####################### +# Identify unique upregulated marker genes for each clusters (based on normalized data) + +## gvg: clustering done, set default assay back to SCT for DGE analyses +DefaultAssay(seu.norm) <- "SCT" + +#prepare object to run differential expression on SCT assay with multiple models +integrated <- PrepSCTFindMarkers(seu.norm) + +saveRDS(seu.norm, paste0(analysisFolder, "/clustered_res", resolution, ".rds")) + +# find marker genes of all clusters (only overexpressed genes -> only.pos=TRUE) +markers_all <- FindAllMarkers(object=seu.norm, only.pos=TRUE, min.pct = 0.25) +markers_all <- markers_all[markers_all$p_val_adj < 0.05,] + +#average log2 fold change is calculated: +# data1 = log2(mean(expm1(ExpressionData1)) + 1) +# data2 = log2(mean(expm1(ExpressionData2)) + 1) +# avg_log2FC = data1-data2 + +#average log2 fold change is calculated (for "negbinom", "poisson", or "DESeq2" DE tests): +# data1 = log2(mean(ExpressionData1) + 1) +# data2 = log2(mean(ExpressionData2) + 1) +# avg_log2FC = data1-data2 + + +# only genes unique between clusters +markers_all_single <- markers_all[markers_all$gene %in% names(table(markers_all$gene))[table(markers_all$gene) == 1],] + +# nb of unique marker genes per cluster +summary <- as.data.frame(table(markers_all_single$cluster)) +colnames(summary) <- c("cluster", "Nb unique marker genes") +summary + +write.xlsx(markers_all_single %>% dplyr::select(cluster, everything()), paste0(analysisFolder, "/uniqueUpregulatedGenes_perCluster_res", resolution, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) +# *Note: p_val: p-value, avg_log2FC: log2 fold-change in the average expression between the two groups. +# Positive values indicate that the gene is more highly expressed within the cluster, pct.1: percentage +# of spots where the gene is detected within the cluster, pct.2: percentage of spots where the gene is +# detected outside the cluster, p_val_adj: adjusted p-value based on Bonferroni correction using all +# genes in the dataset, cluster: cluster ID, gene: gene name + +# Heatmap of top 8 unique marker genes between clusters: +# heatmap of genes by cluster for the top 8 marker genes per cluster +top <- markers_all_single %>% group_by(cluster) %>% top_n(8, avg_log2FC) +DoHeatmap(object = seu.norm, features = top$gene) + + +# Dotplot of top 8 **unique marker genes** between clusters: +DotPlot(seu.norm, features = top$gene, dot.scale = 8) + coord_flip() + RotatedAxis() + + +# Identify all upregulated genes for each cluster ---------------- +# For each cluster, we identify all genes that are upregulated in spots in the cluster compared to spots outside the cluster. +# In contrast to the analysis of unique upregulated marker genes above, we do not require that a gene is specific to a single +# cluster, i.e. a given gene can be detected for more than one cluster. + +for(cluster in sort(unique(markers_all$cluster))){ + + markers_all_cluster <- markers_all[markers_all$cluster == cluster, ] + rownames(markers_all_cluster) <- markers_all_cluster$gene + + print(paste0("### Cluster: ",cluster)) + print(paste0("Markers genes for cluster ",cluster, ":")) + write.xlsx(markers_all_cluster, paste0(analysisFolder, "/upregulatedGenes_res", resolution, "_cluster", cluster, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + # *Note: p_val: p-value, avg_log2FC: log2 fold-change in the average expression between the two groups. Positive values + # indicate that the gene is more highly expressed within the cluster, pct.1: percentage of spots where the gene is detected + # within the cluster, pct.2: percentage of spots where the gene is detected outside the cluster, p_val_adj: adjusted p-value + # based on Bonferroni correction using all genes in the dataset, cluster: cluster ID, gene: gene name + + print(FeaturePlot(seu.norm, head(markers_all_cluster$gene, n=4), cols = c("lightgrey", "blue"), ncol = 2)) + print(SpatialFeaturePlot(seu.norm, features = head(markers_all_cluster$gene, n=4), alpha = c(0.1, 1))) +} + + + + +# 4. Identify spatial variable features ########## +# Search for features exhibiting spatial patterning in the absence of pre-annotation. The method, is inspired +# by the Trendsceek, which models spatial transcriptomics data as a mark point process and computes a variogram, +# which identifies genes whose expression level is dependent on their spatial location. More specifically, this +# process calculates gamma(r) values measuring the dependence between two spots a certain "r" distance apart. +# By default, we use an r-value of '5' in these analyses, and only compute these values for variable genes +# (where variation is calculated independently of spatial location) to save time.\n\nTop 9 of spatial variable genes: + +seu.norm <- FindSpatiallyVariableFeatures(seu.norm, assay = "SCT", features = rownames(GetAssayData(seu.norm, slot = "scale.data")), selection.method = "markvariogram") + +spatialFeatures <- SVFInfo(seu.norm, selection.method="markvariogram", status=TRUE) +spatialFeatures <- na.omit(spatialFeatures) +spatialFeatures <- spatialFeatures %>% filter(variable == TRUE) %>% arrange(rank) +spatialFeatures <- spatialFeatures[,c("r.metric.5"), drop=FALSE] + +write.xlsx(spatialFeatures, paste0(analysisFolder, "/spatialVariableFeatures_res", resolution, "_", names, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) +top.features <- head(SpatiallyVariableFeatures(integrated, selection.method = "markvariogram"), 9) +SpatialFeaturePlot(seu.norm, features=top.features, ncol=3, alpha = c(0.1, 1)) + + + + + + + +############ Heidi: until here ##################### + + + + +# 5. Cell type identification (only human or mouse) +# ScType a computational method for automated selection of marker genes based merely on scRNA-seq data. +# Reference from CellMarker 2.0: One of the largest available databases for cell-type marker. Cell markers +# of different cell types from different tissues with more than one evidence were used. + +# Be aware this is an automatic annotation and requires manual curation! + + +# DB file +doAnnotation <- TRUE +species <- "" +if(msigdb_species == "Homo sapiens"){ + species <- "Human" +} else if(msigdb_species == "Mus musculus"){ + species <- "Mouse" +} else { + cat("**--> nothing done, as species is neither human nor mouse!") + doAnnotation <- FALSE +} +``` + +```{r databaseSetup , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} +# load gene set preparation function +library("HGNChelper") +source("utils/sc-type/gene_sets_prepare.R") +# load cell type annotation function +source("utils/sc-type/sctype_score_.R") + +# create own db (downloaded from CellMarker 2.0) +tryCatch({ + own_db_file <- read.xlsx("http://bio-bigdata.hrbmu.edu.cn/CellMarker/CellMarker_download_files/file/Cell_marker_All2.xlsx") +},error =function(e){ + own_db_file <<- read.xlsx("utils/Cell_marker_All.xlsx") +} +) +own_db_file <- own_db_file[own_db_file$species == species, c("tissue_class", "cell_name", "Symbol")] + +# remove empty symbol lines +own_db_file <- own_db_file[!is.na(own_db_file$Symbol), ] + +# merge marker names into one line and remove entries with single support +own_db_file$tissue_cell_name <- paste0(own_db_file$tissue_class, own_db_file$cell_name) +own_db <- own_db_file %>% group_by(tissue_cell_name) %>% + mutate(Symbol = paste(unique(Symbol[duplicated(Symbol)]), collapse=",")) %>% + distinct() %>% arrange(tissue_class) + +# remove empty symbol lines +own_db <- own_db[own_db$Symbol != "", ] + + +own_db <- own_db[,c("tissue_class", "cell_name", "Symbol")] +colnames(own_db) <- c("tissueType", "cellName", "geneSymbolmore1") +own_db$geneSymbolmore2 <- "" +own_db$shortName <- "" + +#remove entries with < 3 cellName +own_db <- own_db[own_db$tissueType %in% names(table(own_db$tissueType)[table(own_db$tissueType) > 2]),] + +write.xlsx(own_db, paste0(analysisFolder, "/referenceDB.xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=FALSE, append=FALSE) + +# automatically detect tissue type of the dataset if not set +if(tissue == "" || length(grep(tissue, own_db$tissueType, ignore.case=TRUE, value=TRUE)) == 0){ + source("utils/sc-type/auto_detect_tissue_type.R") + # guess a tissue type + # if saled = TRUE, make sure the data is scaled, as seuratObject[[assay]]@scale.data is used. + # If you just created a Seurat object, without any scaling and normalization, set scaled = FALSE, seuratObject[[assay]]@counts will be used + + tissue_guess <- auto_detect_tissue_type(path_to_db_file = paste0(analysisFolder, "/referenceDB.xlsx"), seuratObject = integrated, scaled = TRUE, assay = "SCT") + write.table(tissue_guess, paste0(analysisFolder, "/tissue_guess.txt"), quote=FALSE, row.names=FALSE) + + tissue <- tissue_guess$tissue[1] +} else { + #make sure its written the correct way + tissue <- grep(tissue, own_db$tissueType, ignore.case=TRUE, value=TRUE)[1] +} +``` + +```{r autodetectedTissue , eval=doAnnotation,echo=FALSE, warning=FALSE, message=FALSE,results="asis"} +cat("\n") +cat(paste0("**-> Tissue autodetection: ", tissue, "**")) + +``` + +`r if(doAnnotation){"
"}` + +`r if(doAnnotation){"### Top 5 cell type annotation for each cluster"}` + +```{r cellTypeIdentification , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE} +# prepare gene sets +gs_list <- gene_sets_prepare(paste0(analysisFolder, "/referenceDB.xlsx"), tissue) + +# get cell-type by cell matrix +es.max <- sctype_score(scRNAseqData = integrated[["SCT"]]@scale.data, scaled = TRUE, + gs = gs_list$gs_positive, gs2 = gs_list$gs_negative) + +# merge by cluster +cL_resutls <- do.call("rbind", lapply(unique(integrated@meta.data$finalCluster), function(cl){ + es.max.cl <- sort(rowSums(es.max[ ,rownames(integrated@meta.data[integrated@meta.data$finalCluster==cl, ])]), decreasing = !0) + head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(integrated@meta.data$finalCluster==cl)), 10) +})) +sctype_scores <- cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores) +sctype_scores_top5 <- cL_resutls %>% group_by(cluster) %>% top_n(n = 5, wt = scores) + +sctype_scores_top5 %>% kable() %>% kable_styling() %>% + scroll_box(width = "100%", height = "300px") +``` +`r if(doAnnotation){"*Note: cluster: cluster id, type: cell type annotation, scores: Cell-type specificity score provides a quantitative measure of how uniquely markers identify a specific cell-type of the given tissue. A cell-type with the highest ScType score is used for assignment to the cluster. Low ScType score (less than quarter the number of cells in a cluster), or a negative ScType score indicates a low-confidence cell-type annotation, which will be assigned as “unknown” cell type, ncells: number of cells in this cluster*"}` + +`r if(doAnnotation){"
"}` + +```{r cellTypeIdentification_dimPlot , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=5, fig.width = 9} +# set low-confident (low ScType score) clusters to "unknown" +sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown" + +# UMAP +integrated@meta.data$clusterAnnotation = "" +for(j in unique(sctype_scores$cluster)){ + cl_type = sctype_scores[sctype_scores$cluster==j,]; + integrated@meta.data$clusterAnnotation[integrated@meta.data$finalCluster == j] = as.character(cl_type$type[1]) +} + +DimPlot(integrated, reduction = "umap", label = TRUE, repel = TRUE, group.by = 'clusterAnnotation') +``` + + +`r if(doAnnotation){"
"}` + +`r if(doAnnotation){"Bubble plot showing all the cell types that were considered by ScType for cluster annotation. The outter (grey) bubbles correspond to each cluster (the bigger bubble, the more cells in the cluster), while the inner bubbles correspond to considered cell types for each cluster, with the biggest bubble corresponding to assigned cell type."}` + +```{r cellTypeIdentification_bubblePlot , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=7, fig.width = 9} +# load libraries +library(ggraph) +library(igraph) + +# prepare edges +cL_resutls <- cL_resutls[order(cL_resutls$cluster),] +edges <- cL_resutls +edges$type <- paste0(edges$type,"_",edges$cluster) +edges$cluster <- paste0("cluster ", edges$cluster) +edges <- edges[,c("cluster", "type")] +colnames(edges) <- c("from", "to") +rownames(edges) <- NULL + +# prepare nodes +nodes_lvl1 <- sctype_scores[,c("cluster", "ncells")] +nodes_lvl1$cluster <- paste0("cluster ", nodes_lvl1$cluster) +nodes_lvl1$Colour <- "#f1f1ef" +nodes_lvl1$ord <- 1 +nodes_lvl1$realname <- nodes_lvl1$cluster +nodes_lvl1 <- as.data.frame(nodes_lvl1) +nodes_lvl2 <- c() +ccolss <- hue_pal()(length(unique(cL_resutls$cluster))) +for (i in 1:length(unique(cL_resutls$cluster))){ + dt_tmp <- cL_resutls[cL_resutls$cluster == unique(cL_resutls$cluster)[i], ] + nodes_lvl2 <- rbind(nodes_lvl2, data.frame(cluster = paste0(dt_tmp$type,"_",dt_tmp$cluster), ncells = dt_tmp$scores, Colour = ccolss[i], ord = 2, realname = dt_tmp$type)) +} +nodes <- rbind(nodes_lvl1, nodes_lvl2) +nodes$ncells[nodes$ncells<1] <- 1 +files_db <- own_db[,c("cellName","shortName")] +#files_db <- openxlsx::read.xlsx(db_)[,c("cellName","shortName")] +files_db <- unique(files_db) +nodes <- merge(nodes, files_db, all.x = T, all.y = F, by.x = "realname", by.y = "cellName", sort = F) +nodes$shortName[is.na(nodes$shortName) | nodes$shortName == ""] <- nodes$realname[is.na(nodes$shortName) | nodes$shortName == ""] +nodes <- nodes[,c("cluster", "ncells", "Colour", "ord", "shortName", "realname")] + +mygraph <- graph_from_data_frame(edges, vertices=nodes) + +# Make the graph +plotList <- list() +plotList[[1]] <- ggraph(mygraph, layout = 'circlepack', weight=I(ncells)) + + geom_node_circle(aes(filter=ord==1,fill=I("#F5F5F5"), colour=I("#D3D3D3")), alpha=0.9) + + geom_node_circle(aes(filter=ord==2,fill=I(Colour), colour=I("#D3D3D3")), alpha=0.9) + + theme_void() + + geom_node_text(aes(filter=ord==2, label=shortName, colour=I("#000000"), fill="white", repel = !1, parse = T, size = I(log(ncells,25)*1.2))) + + geom_node_label(aes(filter=ord==1, label=shortName, colour=I("#000000"), size = I(3), fill="white", parse = T), repel = !0, segment.linetype="dotted") + +plotList[[2]] <- DimPlot(integrated, reduction = "umap", label = TRUE, repel = TRUE, cols = ccolss) +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) + + + +# reset finalCluster +integrated$finalCluster_number <- integrated$finalCluster + +Idents(integrated) <- integrated@meta.data$clusterAnnotation +integrated$finalCluster <- Idents(integrated) +saveRDS(integrated, paste0(analysisFolder, "/clustered_cellTypeAnnotated_res", resolution, "_", comparison, ".rds")) +``` + + +
+ + `r if(length(uniqueConditions) <= 1){"# 8. Cell-type deconvolution"}` +`r if(length(uniqueConditions) > 1){"# 10. Cell-type deconvolution"}` +Tools: + + * [STdeconvolve](https://bioconductor.org/packages/release/bioc/html/STdeconvolve.html) +* [CellMarker](http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html) + +STdeconvolve is an unsupervised, reference-free approach to infer latent cell-type proportions and transcriptional profiles within multi-cellular spatially-resolved pixels from spatial transcriptomics (ST) datasets. STdeconvolve builds on latent Dirichlet allocation (LDA), a generative statistical model commonly used in natural language processing. In the context of ST data, given a count matrix of gene expression in multi-cellular ST pixels, STdeconvolve applies LDA to infer the putative transcriptional profile for each cell-type and the proportional representation of each cell-type in each multi-cellular ST pixel. + +
+ + ```{r STdeconvolve_selectK , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} + +library(STdeconvolve) + +# this is the genes x barcode sparse count matrix +cd <- integrated[["SCT"]]@counts + +# xy coordinates +addition <- 0 +pos <- data.frame() +posList <- list() +for(name in names){ + posImage <- GetTissueCoordinates(integrated, image = name) + colnames(posImage) <- c("y", "x") + #reverse y axis + posImage$y <- 2*tail(posImage$y,1)-posImage$y + posList[[name]] <- posImage + posImage$y <- posImage$y + addition + pos <- rbind(pos, posImage) + addition <- addition + 1000 +} + + +## feature select for genes +corpus <- restrictCorpus(cd, + removeAbove = 1.0, + removeBelow = 0.05, + alpha = 0.05, + plot = FALSE, + verbose = TRUE) + +#choose optimal K +ldas <- fitLDA(t(as.matrix(corpus)), Ks = seq(2, 15, by = 1), + perc.rare.thresh = 0.05, + plot=TRUE, + verbose=TRUE) + +## select model +optLDA <- optimalModel(models = ldas, opt = "kneed") +optLDA2 <- optimalModel(models = ldas, opt = "min") + +if(optLDA2@k > optLDA@k){ + optLDA <- optimalModel(models = ldas, opt = round(mean(c(optLDA2@k, optLDA@k)))) +} +``` + +**-> K of `r optLDA@k` chosen** + +
+ + #### Scatterpies of proportions of all deconvolved cell-types: + Cell-types from spots that contribute less than 5% of the spot proportions are removed and the deconvolved transcriptional profiles are scaled by 1000 + +```{r STdeconvolve_proportionPlot , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=9, fig.width = 9} +## Extract pixel cell-type proportions (theta) and cell-type gene expression +## profiles (beta) for the given dataset. +## We can also remove cell-types from pixels that contribute less than 5% of the +## pixel proportion and scale the deconvolved transcriptional profiles by 1000 +results <- getBetaTheta(optLDA, + perc.filt = 0.05, + betaScale = 1000) + +deconProp <- results$theta +deconGexp <- results$beta + +plotList <- list() +for(name in names){ + cells <- rownames(posList[[name]]) + #visualize the proportion of each deconvolved cell-type across the original spatially resolved pixels. (long runtime!) + plotList[[name]] <- vizAllTopics(deconProp[cells,], posList[[name]], r=1.8, lwd = 0, showLegend = TRUE, plotTitle = name) + + ggplot2::guides(colour = "none") +} +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) + +#plot seperately for each celltype +#vizTopic(theta = deconProp, pos = pos, topic = "5", plotTitle = "X5", +# size = 5, stroke = 1, alpha = 0.5, +# low = "white", +# high = "red") +``` + + +#### Compare deconvolved cell-types to clustering +Correlation matrix of the correlations between the proportions of each cell-type and the clustering + +```{r STdeconvolve_correlation , eval=TRUE,echo=FALSE, warning=FALSE, message=FALSE, results="asis"} +cluster_proxyTheta <- model.matrix(~ 0 + integrated@meta.data$finalCluster_number) +rownames(cluster_proxyTheta) <- rownames(integrated@meta.data) + +corMat_prop <- getCorrMtx(m1 = as.matrix(cluster_proxyTheta), + m2 = deconProp, + type = "t") + +rownames(corMat_prop) <- paste0("cluster_", levels(integrated@meta.data$finalCluster_number)) +colnames(corMat_prop) <- paste0("decon_", seq(ncol(corMat_prop))) + +correlationPlot(mat = corMat_prop, colLabs = "Transcriptional clusters", rowLabs = "STdeconvolve") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90) + ) +``` + +
+ + ## Marker genes for each deconvoluted cell-type + Marker genes are defined here as genes highly expressed in the deconvolved cell-type (count > 2) that also have the highest log2(fold change) (>1) when comparing the deconvolved cell-type’s expression profile to the average of all other deconvolved cell-types expression profiles. + +```{r STdeconvolve_markerGenes , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} + +#Now, let’s get the differentially expressed genes for each deconvolved cell-type transcriptional profile and label the top expressed genes for each cell-type: + +markerList <- lapply(colnames(deconProp), function(celltype) { + celltype <- as.numeric(celltype) + ## highly expressed in cell-type of interest + highgexp <- names(which(deconGexp[celltype,] > 2)) + + ## high log2(fold-change) compared to other deconvolved cell-types + markerTable <- data.frame(celltype=celltype, gene=highgexp, log2FC=log2(deconGexp[celltype,highgexp]/colMeans(deconGexp[-celltype,highgexp])), + exp.profile.celltype=deconGexp[celltype,highgexp], + mean.exp.profile.other.celltypes=colMeans(deconGexp[-celltype,highgexp])) + + markerTable <- markerTable[markerTable$log2FC > 1,] + + markerTable <- markerTable %>% arrange(desc(log2FC)) + return(markerTable) +}) + + +``` + + +```{r STdeconvolve_markerGenesPlot , eval=TRUE,echo=FALSE, warning=FALSE, message=FALSE, fig.height=10, fig.width=13, results="asis"} +celltype <- 0 +for(markerTable in markerList){ + celltype <- celltype + 1; + cat("\n") + cat("\n") + cat(paste0("### Cell-tpye: ",celltype)) + cat("\n") + cat(paste0("Markers genes for cell-type ",celltype, ":")) + cat("\n") + cat("\n") + write.xlsx(markerTable, paste0(analysisFolder, "/STdeconvolve_markerGenes_cellType", i, "_", comparison, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + + cat(markerTable %>% + kable() %>% + kable_styling() %>% + scroll_box(width = "100%", height = "300px")) + cat("\n") + cat("*Note: celltype: deconvolved cell-type, gene: gene name, log2FC: log2 fold-change in comparing deconvolved cell-type expression profile to the average of all other deconvolved cell-types expression profiles. Positive values indicate that the gene is more highly expressed within the deconvolved cell-type, exp.profile.celltype: expression profile of deconvolved cell-type, mean.exp.profile.other.celltypes: mean expression profile of all other deconvolved cell-types*\n") + cat("\n") + cat("\n") + + if(nrow(markerTable) > 0){ + print(FeaturePlot(integrated, head(markerTable$gene, n=4), cols = c("lightgrey", "blue"), ncol = 2)) + } + + if(nrow(markerTable) > 1){ + for(i in c(1:2)){ + if(i <= length(markerTable$gene)){ + #split into images of 2 + imageNames <- names(integrated@images) + j <- 1 + while(j <= length(imageNames)){ + if(j+1 <= length(imageNames)){ + print(SpatialFeaturePlot(integrated, features = markerTable$gene[i], pt.size.factor=point_size, alpha = c(0.1, 1), images=imageNames[j:(j+1)])) + } else { + print(SpatialFeaturePlot(integrated, features = markerTable$gene[i], pt.size.factor=point_size, alpha = c(0.1, 1), images=imageNames[j]) + ggtitle(imageNames[j]) + theme(plot.title = element_text(hjust = 0.5))) + } + j <- j + 2 + } + + } + } + cat("\n") + cat("\n") + } + + + # GO term enrichment of marker genes in a cell-type + all.genes <- colnames(deconGexp) + + marker.genes <- markerTable$gene + + # define geneList as 1 if gene is in marker.genes, 0 otherwise + geneList <- ifelse(all.genes %in% marker.genes, 1, 0) + names(geneList) <- all.genes + + onts <- c("BP", "MF", "CC") + gotable <- data.frame() + for(ont in onts){ + tryCatch({ + # Create topGOdata object + GOdata <- new("topGOdata", ontology = ont, # use biological process ontology + allGenes= geneList, geneSelectionFun = function(x)(x == 1), + annot=annFUN.org, mapping=topgo_db, ID="symbol") + + # Test for enrichment using Fisher's Exact Test + resultFisher <- runTest(GOdata, algorithm = "elim", statistic = "fisher", cutOff=0.05) + gotable.ont <- GenTable(GOdata, Fisher = resultFisher, numChar = 60, topNodes=20) + gotable.ont <- gotable.ont[as.numeric(gotable.ont$Fisher) <= 0.05,] + if(nrow(gotable.ont) > 0){ + gotable <- rbind(gotable, cbind(ont, gotable.ont)) + } + }, + error = function(e) { + cat(paste0("**-> GO enrichment with ", onts, " failed!**")) + cat("\n") + cat("\n") + } + ) + } + + if(nrow(gotable) > 0){ + write.xlsx(gotable, paste0(analysisFolder, "/STdeconvolve_markerGenes_cellType", i,"_enrichedGOterms_", comparison, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + + cat("Enriched GO-terms (top 20 of each ontology domain):") + cat("\n") + cat("\n") + cat(gotable %>% + kable() %>% + kable_styling() %>% + scroll_box(width = "100%", height = "300px")) + cat("\n") + cat("*Note: ont: gene ontology domain, GO.ID: Gene Ontology ID, Term: GO term, Annotated: Number of genes annotated within this GO term, Significant: Number of upregulated genes within this GO term, Expected: Expected number of genes (under null hypothesis) within this GO term, Fisher: p-value of Fisher's exact test*\n") + cat("\n") + cat("\n") + } else{ + cat("\n") + cat("\n") + cat("**-> no GO-terms enriched**") + cat("\n") + cat("
") + cat("
") + } +} +``` + +
+ + + `r if(doAnnotation){"## Cell-type annotation"}` +`r if(doAnnotation){"Given a list of reference gene sets for different cell types, we can performed gene set enrichment analysis (GSEA) on the deconvolved transcriptional profiles to test for significant enrichment of any known ground truth cell-types."}` + +`r if(doAnnotation){"Reference from CellMarker 2.0: One of the largest available databases for cell-type marker. Human cell markers of different cell types from different tissues with more than one evidence were used."}` + +`r if(doAnnotation){"**Be aware this is an automatic annotation and requires manual curation!**"}` + +`r if(doAnnotation){"
"}` +```{r STdeconvolve_annotation , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} +#GSEA annotation +# prepare gene sets +gs_list <- gene_sets_prepare(paste0(analysisFolder, "/referenceDB.xlsx"), tissue) + +celltype_annotations <- annotateCellTypesGSEA(beta = deconGexp , gset = gs_list$gs_positive, qval = 0.05) + +celltype_annotation_table <- data.frame() +i <- 0 +for(result in celltype_annotations$results){ + i <- i + 1 + result <- cbind(data.frame(celltype=i, annotation=rownames(result)), result) + celltype_annotation_table <- rbind(celltype_annotation_table, result) +} +``` + +```{r STdeconvolve_annotation_table , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="asis"} +celltype_annotation_table %>% kable(row.names=FALSE) %>% kable_styling() %>% + scroll_box(width = "100%", height = "300px") +``` +`r if(doAnnotation){"*Note: celltype: deconvolved cell-type, annotation: annotation, p.val: p-value of the enrichment score, q.val: q-value estimation for false discovery rate. This is the significance threshold that should be considered, sscore: enrichment score represent the degree to which a set is over-represented at the top or bottom of the ranked list, edge: edge score*"}` + + +```{r STdeconvolve_proportionPlot2 , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=9, fig.width = 9} + +annotationNames <- as.character(celltype_annotations$predictions) +annotationNames[is.na(annotationNames)] <- "unknown" +colnames(deconProp) <- annotationNames + +plotList <- list() +for(name in names){ + cells <- rownames(posList[[name]]) + #visualize the proportion of each deconvolved cell-type across the original spatially resolved pixels. (long runtime!) + plotList[[name]] <- vizAllTopics(deconProp[cells,], posList[[name]], r=1.8, lwd = 0, showLegend = TRUE, plotTitle = name) + + ggplot2::guides(colour = "none") +} +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) +``` + + +`r if(doAnnotation){"#### Compare deconvolved cell-types annotation to clustering annotation"}` +`r if(doAnnotation){"Correlation matrix of the correlations between the proportions of each cell-type annotation and the clustering annotation"}` + +```{r STdeconvolve_correlationAnnotation , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="asis"} + +deconPropMerged <- t(apply(deconProp, 1, function(x) tapply(x,colnames(deconProp),sum))) + +# correlation with clustering +cluster_proxyTheta <- model.matrix(~ 0 + integrated@meta.data$finalCluster) +rownames(cluster_proxyTheta) <- rownames(integrated@meta.data) + +corMat_prop <- getCorrMtx(m1 = as.matrix(cluster_proxyTheta), + m2 = deconPropMerged, + type = "t") + +rownames(corMat_prop) <- levels(integrated@meta.data$finalCluster) + +correlationPlot(mat = corMat_prop, colLabs = "Transcriptional clusters", rowLabs = "STdeconvolve") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90) + ) +``` + + + + + + + + +sessionInfo() + diff --git a/scripts/spatial_RNAseq_integration-clustering.R b/scripts/spatial_RNAseq_integration-clustering.R new file mode 100644 index 0000000..d244366 --- /dev/null +++ b/scripts/spatial_RNAseq_integration-clustering.R @@ -0,0 +1,763 @@ +# spatial RNAseq clustering +# Interfaculty Bioinformatics Unit (IBU) (https://www.bioinformatics.unibe.ch + +library(Seurat) #v4.1 +library(ggplot2) +library(patchwork) +library(dplyr) +library(ggpubr) +library(clusterProfiler) +library(clustree) +library(future) +library(openxlsx) +library(RColorBrewer) + + + + +#library(STutility) +#library(topGO) +#library(limma) +#library(biomaRt) +#library(msigdbr) +#library(scales) + + +# setup -------------------------- + +# Specify input directories and files +workPath <- "/data/projects/p907_SIBdays_tutorial_-_spatial_transcriptomics/" ### ADJUST + +analysisFolder <- paste0(workPath, "scripts/analysis") + +name <- "Posterior" ### ADJUST +#--------------------------------- + + +seu.norm <- readRDS(paste0(analysisFolder, "/normalized.rds")) + +#parallelization +plan("multisession", workers = availableCores()) +options(future.globals.maxSize = 8000 * 1024^2) + + +seu.norm <- RunPCA(seu.norm, assay = "SCT", npcs = 50, verbose = FALSE) +seu.norm <- RunUMAP(seu.norm, reduction = "pca", dims = 1:50) + +DimPlot(seu.norm, reduction = "pca", group.by = "orig.ident") + NoLegend() +DimPlot(seu.norm, reduction = "umap", group.by = "orig.ident")+ ggtitle("") + NoLegend() + +#seu.norm <- ScaleData(seu.norm, features=rownames(seu.norm), verbose = FALSE) + + +# 1. Dimensionality reduction ######################### + # Tool: [Seurat](https://satijalab.org/seurat/articles/spatial_vignette.html) + +# Selecting which PCs to use: --------- +# To overcome the extensive technical noise in any single gene, Seurat clusters spots based +# on their PCA scores, with each PC essentially representing a metagene that combines information +# across a correlated group of genes. + +# ElbowPlot ranks the principal components based on the variance explained by each. This plot typically +# shows an "elbow", which can be used to assess how many PCs are needed to capture most of the signal in the data. + +use.pcs <- 50 +ElbowPlot(seu.norm, ndims=use.pcs) + +# The sctransform workflow performs more effective normalization, strongly removing technical effects from the data, +# this means that higher PCs are more likely to represent subtle, but biologically relevant, sources of heterogeneity +# - so including them may improve downstream analysis. Therefore, higher number of PC can be used. +# By default we are using the first 50 PCs + + + +# 2. Identifying clusters ######################### +# Seurat implements a graph-based clustering approach. Distances between the spots are calculated based on +# previously identified PCs. Briefly, Seurat identifies clusters of spots by a shared nearest neighbor (SNN) +# modularity optimization based clustering algorithm. First, it identifies k-nearest neighbors (KNN) and constructs +# the SNN graph. Then it optimizes the modularity function to determine clusters. For a full description of the +# algorithms, see Waltman and van Eck (2013) The European Physical Journal B. + +# The FindClusters function implements the procedure, and contains a resolution parameter that sets the granularity +# of the downstream clustering, with increased values leading to a greater number of clusters. + +# Selecting which resolution to use: ------- +resolution_vector <- seq(0.2,2,0.2) +seu.norm <- FindNeighbors(seu.norm, reduction="pca", dims=1:use.pcs) +seu.norm <- FindClusters(object=seu.norm, resolution=resolution_vector, verbose=FALSE) + +resTable <- as.data.frame(sapply(grep("res",colnames(seu.norm@meta.data),value = TRUE), function(x) length(unique(seu.norm@meta.data[,x])))) +colnames(resTable) <- "number of clusters" + +# how many clusters in each resolution +resTable + +for(i in seq(resolution_vector)){ + print(DimPlot(seu.norm, reduction = "umap", label=TRUE, group.by=paste0("SCT_snn_res.", resolution_vector[i])) + labs(title=paste0("res_", resolution_vector[i]))) + print(SpatialDimPlot(seu.norm, label = TRUE, combine=FALSE, label.size = 3, group.by=paste0("SCT_snn_res.", resolution_vector[i]))) +} + +# Plotting clustering trees +# To build a clustering tree we need to look at how spots move as the clustering resolution is increased. +# Each cluster forms a node in the tree and edges are constructed by considering the spots in a cluster +# at a lower resolution that end up in a cluster at the next higher resolution. By connecting clusters +# in this way we can see how clusters are related to each other, which are clearly distinct and which are unstable. + +# If there are nodes with multiple incoming edges, it is a good indication that we have over-clustered the data. + +clustree(seu.norm, prefis="res.") + + + +# -> Set the resolution to 0.4 ############ +resolution <- "0.4" + +# Final clustering: +Idents(seu.norm) <- paste0("SCT_snn_res.", resolution) +seu.norm$finalCluster <- Idents(seu.norm) + +DimPlot(seu.norm, reduction = "umap", label=TRUE, group.by = "ident") +SpatialDimPlot(seu.norm, label = TRUE, combine=FALSE, label.size = 3) + + +# As there are many colors, it can be challenging to visualize which spot belongs to which cluster. +# Therefore, we highlight spots of each cluster seperately: +SpatialDimPlot(seu.norm, combine=TRUE, images=name, cells.highlight = CellsByIdentities(object = seu.norm), cols.highlight = c("red", "NA"), facet.highlight = TRUE, ncol = 4) + + +# Total number of cells per cluster per slice +seu.norm@meta.data %>% dplyr::select(orig.ident, finalCluster) %>% table() + + + +# 3. Identify marker genes for each cluster ####################### +# Identify unique upregulated marker genes for each clusters (based on normalized data) + +## gvg: clustering done, set default assay back to SCT for DGE analyses +DefaultAssay(seu.norm) <- "SCT" + +#prepare object to run differential expression on SCT assay with multiple models +integrated <- PrepSCTFindMarkers(seu.norm) + +saveRDS(seu.norm, paste0(analysisFolder, "/clustered_res", resolution, ".rds")) + +# find marker genes of all clusters (only overexpressed genes -> only.pos=TRUE) +markers_all <- FindAllMarkers(object=seu.norm, only.pos=TRUE, min.pct = 0.25) +markers_all <- markers_all[markers_all$p_val_adj < 0.05,] + +#average log2 fold change is calculated: +# data1 = log2(mean(expm1(ExpressionData1)) + 1) +# data2 = log2(mean(expm1(ExpressionData2)) + 1) +# avg_log2FC = data1-data2 + +#average log2 fold change is calculated (for "negbinom", "poisson", or "DESeq2" DE tests): +# data1 = log2(mean(ExpressionData1) + 1) +# data2 = log2(mean(ExpressionData2) + 1) +# avg_log2FC = data1-data2 + + +# only genes unique between clusters +markers_all_single <- markers_all[markers_all$gene %in% names(table(markers_all$gene))[table(markers_all$gene) == 1],] + +# nb of unique marker genes per cluster +summary <- as.data.frame(table(markers_all_single$cluster)) +colnames(summary) <- c("cluster", "Nb unique marker genes") +summary + +write.xlsx(markers_all_single %>% dplyr::select(cluster, everything()), paste0(analysisFolder, "/uniqueUpregulatedGenes_perCluster_res", resolution, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) +# *Note: p_val: p-value, avg_log2FC: log2 fold-change in the average expression between the two groups. +# Positive values indicate that the gene is more highly expressed within the cluster, pct.1: percentage +# of spots where the gene is detected within the cluster, pct.2: percentage of spots where the gene is +# detected outside the cluster, p_val_adj: adjusted p-value based on Bonferroni correction using all +# genes in the dataset, cluster: cluster ID, gene: gene name + +# Heatmap of top 8 unique marker genes between clusters: +# heatmap of genes by cluster for the top 8 marker genes per cluster +top <- markers_all_single %>% group_by(cluster) %>% top_n(8, avg_log2FC) +DoHeatmap(object = seu.norm, features = top$gene) + + +# Dotplot of top 8 **unique marker genes** between clusters: +DotPlot(seu.norm, features = top$gene, dot.scale = 8) + coord_flip() + RotatedAxis() + + +# Identify all upregulated genes for each cluster ---------------- +# For each cluster, we identify all genes that are upregulated in spots in the cluster compared to spots outside the cluster. +# In contrast to the analysis of unique upregulated marker genes above, we do not require that a gene is specific to a single +# cluster, i.e. a given gene can be detected for more than one cluster. + +for(cluster in sort(unique(markers_all$cluster))){ + + markers_all_cluster <- markers_all[markers_all$cluster == cluster, ] + rownames(markers_all_cluster) <- markers_all_cluster$gene + + print(paste0("### Cluster: ",cluster)) + print(paste0("Markers genes for cluster ",cluster, ":")) + write.xlsx(markers_all_cluster, paste0(analysisFolder, "/upregulatedGenes_res", resolution, "_cluster", cluster, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + # *Note: p_val: p-value, avg_log2FC: log2 fold-change in the average expression between the two groups. Positive values + # indicate that the gene is more highly expressed within the cluster, pct.1: percentage of spots where the gene is detected + # within the cluster, pct.2: percentage of spots where the gene is detected outside the cluster, p_val_adj: adjusted p-value + # based on Bonferroni correction using all genes in the dataset, cluster: cluster ID, gene: gene name + + print(FeaturePlot(seu.norm, head(markers_all_cluster$gene, n=4), cols = c("lightgrey", "blue"), ncol = 2)) + print(SpatialFeaturePlot(seu.norm, features = head(markers_all_cluster$gene, n=4), alpha = c(0.1, 1))) +} + + + + +# 4. Identify spatial variable features ########## +# Search for features exhibiting spatial patterning in the absence of pre-annotation. The method, is inspired +# by the Trendsceek, which models spatial transcriptomics data as a mark point process and computes a variogram, +# which identifies genes whose expression level is dependent on their spatial location. More specifically, this +# process calculates gamma(r) values measuring the dependence between two spots a certain "r" distance apart. +# By default, we use an r-value of '5' in these analyses, and only compute these values for variable genes +# (where variation is calculated independently of spatial location) to save time.\n\nTop 9 of spatial variable genes: + +seu.norm <- FindSpatiallyVariableFeatures(seu.norm, assay = "SCT", features = rownames(GetAssayData(seu.norm, slot = "scale.data")), selection.method = "markvariogram") + +spatialFeatures <- SVFInfo(seu.norm, selection.method="markvariogram", status=TRUE) +spatialFeatures <- na.omit(spatialFeatures) +spatialFeatures <- spatialFeatures %>% filter(variable == TRUE) %>% arrange(rank) +spatialFeatures <- spatialFeatures[,c("r.metric.5"), drop=FALSE] + +write.xlsx(spatialFeatures, paste0(analysisFolder, "/spatialVariableFeatures_res", resolution, "_", names, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) +top.features <- head(SpatiallyVariableFeatures(integrated, selection.method = "markvariogram"), 9) +SpatialFeaturePlot(seu.norm, features=top.features, ncol=3, alpha = c(0.1, 1)) + + + + + + + +############ Heidi: until here ##################### + + + + +# 5. Cell type identification (only human or mouse) +# ScType a computational method for automated selection of marker genes based merely on scRNA-seq data. +# Reference from CellMarker 2.0: One of the largest available databases for cell-type marker. Cell markers +# of different cell types from different tissues with more than one evidence were used. + +# Be aware this is an automatic annotation and requires manual curation! + + +# DB file +doAnnotation <- TRUE +species <- "" +if(msigdb_species == "Homo sapiens"){ + species <- "Human" +} else if(msigdb_species == "Mus musculus"){ + species <- "Mouse" +} else { + cat("**--> nothing done, as species is neither human nor mouse!") + doAnnotation <- FALSE +} +# ``` + +# ```{r databaseSetup , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} +# load gene set preparation function +library("HGNChelper") +source("utils/sc-type/gene_sets_prepare.R") +# load cell type annotation function +source("utils/sc-type/sctype_score_.R") + +# create own db (downloaded from CellMarker 2.0) +tryCatch({ + own_db_file <- read.xlsx("http://bio-bigdata.hrbmu.edu.cn/CellMarker/CellMarker_download_files/file/Cell_marker_All2.xlsx") +},error =function(e){ + own_db_file <<- read.xlsx("utils/Cell_marker_All.xlsx") +} +) +own_db_file <- own_db_file[own_db_file$species == species, c("tissue_class", "cell_name", "Symbol")] + +# remove empty symbol lines +own_db_file <- own_db_file[!is.na(own_db_file$Symbol), ] + +# merge marker names into one line and remove entries with single support +own_db_file$tissue_cell_name <- paste0(own_db_file$tissue_class, own_db_file$cell_name) +own_db <- own_db_file %>% group_by(tissue_cell_name) %>% + mutate(Symbol = paste(unique(Symbol[duplicated(Symbol)]), collapse=",")) %>% + distinct() %>% arrange(tissue_class) + +# remove empty symbol lines +own_db <- own_db[own_db$Symbol != "", ] + + +own_db <- own_db[,c("tissue_class", "cell_name", "Symbol")] +colnames(own_db) <- c("tissueType", "cellName", "geneSymbolmore1") +own_db$geneSymbolmore2 <- "" +own_db$shortName <- "" + +#remove entries with < 3 cellName +own_db <- own_db[own_db$tissueType %in% names(table(own_db$tissueType)[table(own_db$tissueType) > 2]),] + +write.xlsx(own_db, paste0(analysisFolder, "/referenceDB.xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=FALSE, append=FALSE) + +# automatically detect tissue type of the dataset if not set +if(tissue == "" || length(grep(tissue, own_db$tissueType, ignore.case=TRUE, value=TRUE)) == 0){ + source("utils/sc-type/auto_detect_tissue_type.R") + # guess a tissue type + # if saled = TRUE, make sure the data is scaled, as seuratObject[[assay]]@scale.data is used. + # If you just created a Seurat object, without any scaling and normalization, set scaled = FALSE, seuratObject[[assay]]@counts will be used + + tissue_guess <- auto_detect_tissue_type(path_to_db_file = paste0(analysisFolder, "/referenceDB.xlsx"), seuratObject = integrated, scaled = TRUE, assay = "SCT") + write.table(tissue_guess, paste0(analysisFolder, "/tissue_guess.txt"), quote=FALSE, row.names=FALSE) + + tissue <- tissue_guess$tissue[1] +} else { + #make sure its written the correct way + tissue <- grep(tissue, own_db$tissueType, ignore.case=TRUE, value=TRUE)[1] +} +# ``` + +# ```{r autodetectedTissue , eval=doAnnotation,echo=FALSE, warning=FALSE, message=FALSE,results="asis"} +cat("\n") +cat(paste0("**-> Tissue autodetection: ", tissue, "**")) +# +# ``` + +# `r if(doAnnotation){"
"}` + +# `r if(doAnnotation){"### Top 5 cell type annotation for each cluster"}` + +# ```{r cellTypeIdentification , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE} +# prepare gene sets +gs_list <- gene_sets_prepare(paste0(analysisFolder, "/referenceDB.xlsx"), tissue) + +# get cell-type by cell matrix +es.max <- sctype_score(scRNAseqData = integrated[["SCT"]]@scale.data, scaled = TRUE, + gs = gs_list$gs_positive, gs2 = gs_list$gs_negative) + +# merge by cluster +cL_resutls <- do.call("rbind", lapply(unique(integrated@meta.data$finalCluster), function(cl){ + es.max.cl <- sort(rowSums(es.max[ ,rownames(integrated@meta.data[integrated@meta.data$finalCluster==cl, ])]), decreasing = !0) + head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(integrated@meta.data$finalCluster==cl)), 10) +})) +sctype_scores <- cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores) +sctype_scores_top5 <- cL_resutls %>% group_by(cluster) %>% top_n(n = 5, wt = scores) + +sctype_scores_top5 %>% kable() %>% kable_styling() %>% + scroll_box(width = "100%", height = "300px") +# ``` +# `r if(doAnnotation){"*Note: cluster: cluster id, type: cell type annotation, scores: Cell-type specificity score provides a quantitative measure of how uniquely markers identify a specific cell-type of the given tissue. A cell-type with the highest ScType score is used for assignment to the cluster. Low ScType score (less than quarter the number of cells in a cluster), or a negative ScType score indicates a low-confidence cell-type annotation, which will be assigned as “unknown” cell type, ncells: number of cells in this cluster*"}` +# +# `r if(doAnnotation){"
"}` +# +# ```{r cellTypeIdentification_dimPlot , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=5, fig.width = 9} +# set low-confident (low ScType score) clusters to "unknown" +sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown" + +# UMAP +integrated@meta.data$clusterAnnotation = "" +for(j in unique(sctype_scores$cluster)){ + cl_type = sctype_scores[sctype_scores$cluster==j,]; + integrated@meta.data$clusterAnnotation[integrated@meta.data$finalCluster == j] = as.character(cl_type$type[1]) +} + +DimPlot(integrated, reduction = "umap", label = TRUE, repel = TRUE, group.by = 'clusterAnnotation') +``` + + +# `r if(doAnnotation){"
"}` + +# `r if(doAnnotation){"Bubble plot showing all the cell types that were considered by ScType for cluster annotation. The outter (grey) bubbles correspond to each cluster (the bigger bubble, the more cells in the cluster), while the inner bubbles correspond to considered cell types for each cluster, with the biggest bubble corresponding to assigned cell type."}` + +# ```{r cellTypeIdentification_bubblePlot , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=7, fig.width = 9} +# load libraries +library(ggraph) +library(igraph) + +# prepare edges +cL_resutls <- cL_resutls[order(cL_resutls$cluster),] +edges <- cL_resutls +edges$type <- paste0(edges$type,"_",edges$cluster) +edges$cluster <- paste0("cluster ", edges$cluster) +edges <- edges[,c("cluster", "type")] +colnames(edges) <- c("from", "to") +rownames(edges) <- NULL + +# prepare nodes +nodes_lvl1 <- sctype_scores[,c("cluster", "ncells")] +nodes_lvl1$cluster <- paste0("cluster ", nodes_lvl1$cluster) +nodes_lvl1$Colour <- "#f1f1ef" +nodes_lvl1$ord <- 1 +nodes_lvl1$realname <- nodes_lvl1$cluster +nodes_lvl1 <- as.data.frame(nodes_lvl1) +nodes_lvl2 <- c() +ccolss <- hue_pal()(length(unique(cL_resutls$cluster))) +for (i in 1:length(unique(cL_resutls$cluster))){ + dt_tmp <- cL_resutls[cL_resutls$cluster == unique(cL_resutls$cluster)[i], ] + nodes_lvl2 <- rbind(nodes_lvl2, data.frame(cluster = paste0(dt_tmp$type,"_",dt_tmp$cluster), ncells = dt_tmp$scores, Colour = ccolss[i], ord = 2, realname = dt_tmp$type)) +} +nodes <- rbind(nodes_lvl1, nodes_lvl2) +nodes$ncells[nodes$ncells<1] <- 1 +files_db <- own_db[,c("cellName","shortName")] +#files_db <- openxlsx::read.xlsx(db_)[,c("cellName","shortName")] +files_db <- unique(files_db) +nodes <- merge(nodes, files_db, all.x = T, all.y = F, by.x = "realname", by.y = "cellName", sort = F) +nodes$shortName[is.na(nodes$shortName) | nodes$shortName == ""] <- nodes$realname[is.na(nodes$shortName) | nodes$shortName == ""] +nodes <- nodes[,c("cluster", "ncells", "Colour", "ord", "shortName", "realname")] + +mygraph <- graph_from_data_frame(edges, vertices=nodes) + +# Make the graph +plotList <- list() +plotList[[1]] <- ggraph(mygraph, layout = 'circlepack', weight=I(ncells)) + + geom_node_circle(aes(filter=ord==1,fill=I("#F5F5F5"), colour=I("#D3D3D3")), alpha=0.9) + + geom_node_circle(aes(filter=ord==2,fill=I(Colour), colour=I("#D3D3D3")), alpha=0.9) + + theme_void() + + geom_node_text(aes(filter=ord==2, label=shortName, colour=I("#000000"), fill="white", repel = !1, parse = T, size = I(log(ncells,25)*1.2))) + + geom_node_label(aes(filter=ord==1, label=shortName, colour=I("#000000"), size = I(3), fill="white", parse = T), repel = !0, segment.linetype="dotted") + +plotList[[2]] <- DimPlot(integrated, reduction = "umap", label = TRUE, repel = TRUE, cols = ccolss) +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) + + + +# reset finalCluster +integrated$finalCluster_number <- integrated$finalCluster + +Idents(integrated) <- integrated@meta.data$clusterAnnotation +integrated$finalCluster <- Idents(integrated) +saveRDS(integrated, paste0(analysisFolder, "/clustered_cellTypeAnnotated_res", resolution, "_", comparison, ".rds")) +# ``` + + +#
+# +# `r if(length(uniqueConditions) <= 1){"# 8. Cell-type deconvolution"}` +# `r if(length(uniqueConditions) > 1){"# 10. Cell-type deconvolution"}` +Tools: + + # * [STdeconvolve](https://bioconductor.org/packages/release/bioc/html/STdeconvolve.html) +# * [CellMarker](http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html) + +# STdeconvolve is an unsupervised, reference-free approach to infer latent cell-type proportions and transcriptional profiles within multi-cellular spatially-resolved pixels from spatial transcriptomics (ST) datasets. STdeconvolve builds on latent Dirichlet allocation (LDA), a generative statistical model commonly used in natural language processing. In the context of ST data, given a count matrix of gene expression in multi-cellular ST pixels, STdeconvolve applies LDA to infer the putative transcriptional profile for each cell-type and the proportional representation of each cell-type in each multi-cellular ST pixel. + +#
+ + # ```{r STdeconvolve_selectK , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} + +library(STdeconvolve) + +# this is the genes x barcode sparse count matrix +cd <- integrated[["SCT"]]@counts + +# xy coordinates +addition <- 0 +pos <- data.frame() +posList <- list() +for(name in names){ + posImage <- GetTissueCoordinates(integrated, image = name) + colnames(posImage) <- c("y", "x") + #reverse y axis + posImage$y <- 2*tail(posImage$y,1)-posImage$y + posList[[name]] <- posImage + posImage$y <- posImage$y + addition + pos <- rbind(pos, posImage) + addition <- addition + 1000 +} + + +## feature select for genes +corpus <- restrictCorpus(cd, + removeAbove = 1.0, + removeBelow = 0.05, + alpha = 0.05, + plot = FALSE, + verbose = TRUE) + +#choose optimal K +ldas <- fitLDA(t(as.matrix(corpus)), Ks = seq(2, 15, by = 1), + perc.rare.thresh = 0.05, + plot=TRUE, + verbose=TRUE) + +## select model +optLDA <- optimalModel(models = ldas, opt = "kneed") +optLDA2 <- optimalModel(models = ldas, opt = "min") + +if(optLDA2@k > optLDA@k){ + optLDA <- optimalModel(models = ldas, opt = round(mean(c(optLDA2@k, optLDA@k)))) +} +# ``` + +# **-> K of `r optLDA@k` chosen** + + #
+ + #### Scatterpies of proportions of all deconvolved cell-types: + # Cell-types from spots that contribute less than 5% of the spot proportions are removed and the deconvolved transcriptional profiles are scaled by 1000 + +# ```{r STdeconvolve_proportionPlot , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=9, fig.width = 9} +## Extract pixel cell-type proportions (theta) and cell-type gene expression +## profiles (beta) for the given dataset. +## We can also remove cell-types from pixels that contribute less than 5% of the +## pixel proportion and scale the deconvolved transcriptional profiles by 1000 +results <- getBetaTheta(optLDA, + perc.filt = 0.05, + betaScale = 1000) + +deconProp <- results$theta +deconGexp <- results$beta + +plotList <- list() +for(name in names){ + cells <- rownames(posList[[name]]) + #visualize the proportion of each deconvolved cell-type across the original spatially resolved pixels. (long runtime!) + plotList[[name]] <- vizAllTopics(deconProp[cells,], posList[[name]], r=1.8, lwd = 0, showLegend = TRUE, plotTitle = name) + + ggplot2::guides(colour = "none") +} +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) + +#plot seperately for each celltype +#vizTopic(theta = deconProp, pos = pos, topic = "5", plotTitle = "X5", +# size = 5, stroke = 1, alpha = 0.5, +# low = "white", +# high = "red") +# ``` + + +#### Compare deconvolved cell-types to clustering +Correlation matrix of the correlations between the proportions of each cell-type and the clustering + +# ```{r STdeconvolve_correlation , eval=TRUE,echo=FALSE, warning=FALSE, message=FALSE, results="asis"} +cluster_proxyTheta <- model.matrix(~ 0 + integrated@meta.data$finalCluster_number) +rownames(cluster_proxyTheta) <- rownames(integrated@meta.data) + +corMat_prop <- getCorrMtx(m1 = as.matrix(cluster_proxyTheta), + m2 = deconProp, + type = "t") + +rownames(corMat_prop) <- paste0("cluster_", levels(integrated@meta.data$finalCluster_number)) +colnames(corMat_prop) <- paste0("decon_", seq(ncol(corMat_prop))) + +correlationPlot(mat = corMat_prop, colLabs = "Transcriptional clusters", rowLabs = "STdeconvolve") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90) + ) +# ``` + +#
+ + ## Marker genes for each deconvoluted cell-type + Marker genes are defined here as genes highly expressed in the deconvolved cell-type (count > 2) that also have the highest log2(fold change) (>1) when comparing the deconvolved cell-type’s expression profile to the average of all other deconvolved cell-types expression profiles. + +# ```{r STdeconvolve_markerGenes , eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} + +#Now, let’s get the differentially expressed genes for each deconvolved cell-type transcriptional profile and label the top expressed genes for each cell-type: + +markerList <- lapply(colnames(deconProp), function(celltype) { + celltype <- as.numeric(celltype) + ## highly expressed in cell-type of interest + highgexp <- names(which(deconGexp[celltype,] > 2)) + + ## high log2(fold-change) compared to other deconvolved cell-types + markerTable <- data.frame(celltype=celltype, gene=highgexp, log2FC=log2(deconGexp[celltype,highgexp]/colMeans(deconGexp[-celltype,highgexp])), + exp.profile.celltype=deconGexp[celltype,highgexp], + mean.exp.profile.other.celltypes=colMeans(deconGexp[-celltype,highgexp])) + + markerTable <- markerTable[markerTable$log2FC > 1,] + + markerTable <- markerTable %>% arrange(desc(log2FC)) + return(markerTable) +}) + + +# ``` + + +# ```{r STdeconvolve_markerGenesPlot , eval=TRUE,echo=FALSE, warning=FALSE, message=FALSE, fig.height=10, fig.width=13, results="asis"} +celltype <- 0 +for(markerTable in markerList){ + celltype <- celltype + 1; + cat("\n") + cat("\n") + cat(paste0("### Cell-tpye: ",celltype)) + cat("\n") + cat(paste0("Markers genes for cell-type ",celltype, ":")) + cat("\n") + cat("\n") + write.xlsx(markerTable, paste0(analysisFolder, "/STdeconvolve_markerGenes_cellType", i, "_", comparison, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + + cat(markerTable %>% + kable() %>% + kable_styling() %>% + scroll_box(width = "100%", height = "300px")) + cat("\n") + cat("*Note: celltype: deconvolved cell-type, gene: gene name, log2FC: log2 fold-change in comparing deconvolved cell-type expression profile to the average of all other deconvolved cell-types expression profiles. Positive values indicate that the gene is more highly expressed within the deconvolved cell-type, exp.profile.celltype: expression profile of deconvolved cell-type, mean.exp.profile.other.celltypes: mean expression profile of all other deconvolved cell-types*\n") + cat("\n") + cat("\n") + + if(nrow(markerTable) > 0){ + print(FeaturePlot(integrated, head(markerTable$gene, n=4), cols = c("lightgrey", "blue"), ncol = 2)) + } + + if(nrow(markerTable) > 1){ + for(i in c(1:2)){ + if(i <= length(markerTable$gene)){ + #split into images of 2 + imageNames <- names(integrated@images) + j <- 1 + while(j <= length(imageNames)){ + if(j+1 <= length(imageNames)){ + print(SpatialFeaturePlot(integrated, features = markerTable$gene[i], pt.size.factor=point_size, alpha = c(0.1, 1), images=imageNames[j:(j+1)])) + } else { + print(SpatialFeaturePlot(integrated, features = markerTable$gene[i], pt.size.factor=point_size, alpha = c(0.1, 1), images=imageNames[j]) + ggtitle(imageNames[j]) + theme(plot.title = element_text(hjust = 0.5))) + } + j <- j + 2 + } + + } + } + cat("\n") + cat("\n") + } + + + # GO term enrichment of marker genes in a cell-type + all.genes <- colnames(deconGexp) + + marker.genes <- markerTable$gene + + # define geneList as 1 if gene is in marker.genes, 0 otherwise + geneList <- ifelse(all.genes %in% marker.genes, 1, 0) + names(geneList) <- all.genes + + onts <- c("BP", "MF", "CC") + gotable <- data.frame() + for(ont in onts){ + tryCatch({ + # Create topGOdata object + GOdata <- new("topGOdata", ontology = ont, # use biological process ontology + allGenes= geneList, geneSelectionFun = function(x)(x == 1), + annot=annFUN.org, mapping=topgo_db, ID="symbol") + + # Test for enrichment using Fisher's Exact Test + resultFisher <- runTest(GOdata, algorithm = "elim", statistic = "fisher", cutOff=0.05) + gotable.ont <- GenTable(GOdata, Fisher = resultFisher, numChar = 60, topNodes=20) + gotable.ont <- gotable.ont[as.numeric(gotable.ont$Fisher) <= 0.05,] + if(nrow(gotable.ont) > 0){ + gotable <- rbind(gotable, cbind(ont, gotable.ont)) + } + }, + error = function(e) { + cat(paste0("**-> GO enrichment with ", onts, " failed!**")) + cat("\n") + cat("\n") + } + ) + } + + if(nrow(gotable) > 0){ + write.xlsx(gotable, paste0(analysisFolder, "/STdeconvolve_markerGenes_cellType", i,"_enrichedGOterms_", comparison, ".xlsx"), sheetName="Sheet1", colNames=TRUE, rowNames=TRUE, append=FALSE) + + cat("Enriched GO-terms (top 20 of each ontology domain):") + cat("\n") + cat("\n") + cat(gotable %>% + kable() %>% + kable_styling() %>% + scroll_box(width = "100%", height = "300px")) + cat("\n") + cat("*Note: ont: gene ontology domain, GO.ID: Gene Ontology ID, Term: GO term, Annotated: Number of genes annotated within this GO term, Significant: Number of upregulated genes within this GO term, Expected: Expected number of genes (under null hypothesis) within this GO term, Fisher: p-value of Fisher's exact test*\n") + cat("\n") + cat("\n") + } else{ + cat("\n") + cat("\n") + cat("**-> no GO-terms enriched**") + cat("\n") + cat("
") + cat("
") + } +} +# ``` + +#
+# +# +# `r if(doAnnotation){"## Cell-type annotation"}` +# `r if(doAnnotation){"Given a list of reference gene sets for different cell types, we can performed gene set enrichment analysis (GSEA) on the deconvolved transcriptional profiles to test for significant enrichment of any known ground truth cell-types."}` +# +# `r if(doAnnotation){"Reference from CellMarker 2.0: One of the largest available databases for cell-type marker. Human cell markers of different cell types from different tissues with more than one evidence were used."}` +# +# `r if(doAnnotation){"**Be aware this is an automatic annotation and requires manual curation!**"}` +# +# `r if(doAnnotation){"
"}` +# ```{r STdeconvolve_annotation , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide"} +#GSEA annotation +# prepare gene sets +gs_list <- gene_sets_prepare(paste0(analysisFolder, "/referenceDB.xlsx"), tissue) + +celltype_annotations <- annotateCellTypesGSEA(beta = deconGexp , gset = gs_list$gs_positive, qval = 0.05) + +celltype_annotation_table <- data.frame() +i <- 0 +for(result in celltype_annotations$results){ + i <- i + 1 + result <- cbind(data.frame(celltype=i, annotation=rownames(result)), result) + celltype_annotation_table <- rbind(celltype_annotation_table, result) +} +``` + +# ```{r STdeconvolve_annotation_table , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="asis"} +celltype_annotation_table %>% kable(row.names=FALSE) %>% kable_styling() %>% + scroll_box(width = "100%", height = "300px") +# ``` +# `r if(doAnnotation){"*Note: celltype: deconvolved cell-type, annotation: annotation, p.val: p-value of the enrichment score, q.val: q-value estimation for false discovery rate. This is the significance threshold that should be considered, sscore: enrichment score represent the degree to which a set is over-represented at the top or bottom of the ranked list, edge: edge score*"}` + + +# ```{r STdeconvolve_proportionPlot2 , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="hide", fig.height=9, fig.width = 9} + +annotationNames <- as.character(celltype_annotations$predictions) +annotationNames[is.na(annotationNames)] <- "unknown" +colnames(deconProp) <- annotationNames + +plotList <- list() +for(name in names){ + cells <- rownames(posList[[name]]) + #visualize the proportion of each deconvolved cell-type across the original spatially resolved pixels. (long runtime!) + plotList[[name]] <- vizAllTopics(deconProp[cells,], posList[[name]], r=1.8, lwd = 0, showLegend = TRUE, plotTitle = name) + + ggplot2::guides(colour = "none") +} +ggarrange(plotlist=plotList, ncol=1, nrow=1, hjust=-1) +# ``` + + +# `r if(doAnnotation){"#### Compare deconvolved cell-types annotation to clustering annotation"}` +# `r if(doAnnotation){"Correlation matrix of the correlations between the proportions of each cell-type annotation and the clustering annotation"}` + +# ```{r STdeconvolve_correlationAnnotation , eval=doAnnotation, echo=FALSE, warning=FALSE, message=FALSE, results="asis"} + +deconPropMerged <- t(apply(deconProp, 1, function(x) tapply(x,colnames(deconProp),sum))) + +# correlation with clustering +cluster_proxyTheta <- model.matrix(~ 0 + integrated@meta.data$finalCluster) +rownames(cluster_proxyTheta) <- rownames(integrated@meta.data) + +corMat_prop <- getCorrMtx(m1 = as.matrix(cluster_proxyTheta), + m2 = deconPropMerged, + type = "t") + +rownames(corMat_prop) <- levels(integrated@meta.data$finalCluster) + +correlationPlot(mat = corMat_prop, colLabs = "Transcriptional clusters", rowLabs = "STdeconvolve") + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90) + ) +# ``` + + +# + + + + + +sessionInfo() +