Load packages, functions and data
library(Seurat)
library(tidyverse)
library(readxl)
source("R/ReadPackages.R")
source("R/Functions.R")
# Read CODEX data (only meta data)
# Available at BioStudies database (https://www.ebi.ac.uk/biostudies/) under accession number S-BIAD565
codex_annotation <- data.table::fread("data/cells_annotation.csv") %>% tibble() %>%
filter(Merged_final!="na")
Function
# Run nearest neighbor analysis
run_NNanalysis <- function(data=NULL, plan_session="sequential", n_workers=1, nn=25, regions=NULL,
add.prop=TRUE){
require(future)
require(future.apply)
require(FNN)
if(plan_session=="sequential"){
plan(strategy = plan_session)
}
if(plan_session=="multisession"){
plan(strategy = plan_session, workers=n_workers)
}
options(future.globals.maxSize= 8000*1024^2)
output <- future_lapply(regions, function(r){
df <- data %>%
filter(unique_region==r)
df_nn <- get.knn(select(df,x,y), nn)
df_nn <- df_nn$nn.index %>%
`rownames<-`(df$unique_cell_id) %>%
data.frame() %>%
rownames_to_column("unique_cell_id") %>%
pivot_longer(cols = 2:ncol(.)) %>%
left_join(., df %>% select(unique_cell_id_nn=unique_cell_id, Merged_final) %>% mutate(value=1:nrow(.)), by="value") %>%
select(-value)
if(add.prop==TRUE){
mat <- df_nn %>%
add_prop(vars = c("unique_cell_id", "Merged_final"), group.vars = 1) %>%
pivot_wider(names_from = "Merged_final", values_from = "Prop", values_fill = 0) %>%
column_to_rownames("unique_cell_id")
mat <- mat[names(which(rowSums(mat)>0)), ]
} else {
mat <- df_nn
}
return(mat)
})
plan(strategy = "sequential")
mat <- do.call(rbind, output)
return(mat)
}
Run Analysis
# Read and process codex data (uses future package)
nn <- run_NNanalysis(data = codex_annotation, regions = unique(codex_annotation$unique_region),
plan_session = "multisession",
add.prop=TRUE,
n_workers = 10,
nn = 25)
# Use k-means clustering to determine neighborhood classes
set.seed(1)
hclust <- kmeans(nn, centers = 10, iter.max = 50, nstart = 1)
nn_classes <- data.frame(Region=hclust$cluster) %>%
mutate(Region=as.character(Region)) %>%
rownames_to_column("unique_cell_id")
# Use PCA on nearst neighbour matrix
pca_codex <- prcomp(nn, scale. = T, center = T)
# Remove objects
rm(list = setdiff(ls(), c("nn", "nn_classes", "pca_codex", "run_NNanalysis")))
Save
save.image("output/Neighborhood_results.RData")
Session Info
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /g/easybuild/x86_64/CentOS/7/haswell/software/FlexiBLAS/3.0.4-GCC-11.2.0/lib64/libflexiblas.so.3.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] FNN_1.1.3 future.apply_1.8.1 future_1.23.0 pamr_1.56.1
## [5] cluster_2.1.2 glmnet_4.1-2 Matrix_1.5-1 immunarch_0.7.0
## [9] data.table_1.14.2 dtplyr_1.2.2 rmdformats_1.0.4 ggplotify_0.1.0
## [13] ggraph_2.0.6 igraph_1.3.5 ggrastr_1.0.1 ggtext_0.1.1
## [17] ggalluvial_0.12.3 maxstat_0.7-25 survival_3.2-13 survminer_0.4.9
## [21] ggridges_0.5.3 cowplot_1.1.1 R.utils_2.11.0 R.oo_1.24.0
## [25] R.methodsS3_1.8.1 caret_6.0-90 lattice_0.20-45 patchwork_1.1.2
## [29] rstatix_0.7.0 ggpubr_0.4.0 ggrepel_0.9.1 matrixStats_0.61.0
## [33] scales_1.2.1 RColorBrewer_1.1-3 viridis_0.6.2 viridisLite_0.4.1
## [37] readxl_1.4.1 forcats_0.5.1 stringr_1.4.1 dplyr_1.0.10
## [41] purrr_0.3.4 readr_2.1.2 tidyr_1.2.1 tibble_3.1.8
## [45] ggplot2_3.3.6 tidyverse_1.3.1 SeuratObject_4.0.4 Seurat_4.1.0
## [49] knitr_1.40
##
## loaded via a namespace (and not attached):
## [1] scattermore_0.8 prabclus_2.3-2 ModelMetrics_1.2.2.2
## [4] exactRankTests_0.8-34 irlba_2.3.5 rpart_4.1-15
## [7] doParallel_1.0.17 generics_0.1.3 RANN_2.6.1
## [10] tzdb_0.3.0 rlist_0.4.6.2 spatstat.data_2.1-2
## [13] xml2_1.3.2 lubridate_1.8.0 httpuv_1.6.6
## [16] assertthat_0.2.1 gower_0.2.2 xfun_0.33
## [19] hms_1.1.2 jquerylib_0.1.4 evaluate_0.16
## [22] promises_1.2.0.1 DEoptimR_1.0-11 fansi_1.0.3
## [25] dbplyr_2.1.1 km.ci_0.5-2 DBI_1.1.2
## [28] htmlwidgets_1.5.4 spatstat.geom_2.3-2 stringdist_0.9.8
## [31] stats4_4.1.2 ellipsis_0.3.2 backports_1.4.1
## [34] bookdown_0.29 deldir_1.0-6 vctrs_0.4.2
## [37] ROCR_1.0-11 abind_1.4-5 cachem_1.0.6
## [40] withr_2.5.0 ggforce_0.4.0 robustbase_0.95-0
## [43] sctransform_0.3.3 mclust_5.4.10 goftest_1.2-3
## [46] ape_5.6-2 lazyeval_0.2.2 crayon_1.5.2
## [49] recipes_0.1.17 pkgconfig_2.0.3 tweenr_2.0.2
## [52] nlme_3.1-153 vipor_0.4.5 nnet_7.3-16
## [55] rlang_1.0.6 globals_0.14.0 diptest_0.76-0
## [58] lifecycle_1.0.2 miniUI_0.1.1.1 modelr_0.1.8
## [61] cellranger_1.1.0 polyclip_1.10-0 lmtest_0.9-39
## [64] phangorn_2.10.0 ggseqlogo_0.1 KMsurv_0.1-5
## [67] carData_3.0-5 zoo_1.8-9 reprex_2.0.1
## [70] beeswarm_0.4.0 GlobalOptions_0.1.2 pheatmap_1.0.12
## [73] png_0.1-7 KernSmooth_2.23-20 pROC_1.18.0
## [76] shape_1.4.6 parallelly_1.30.0 spatstat.random_2.1-0
## [79] gridGraphics_0.5-1 ggsignif_0.6.3 magrittr_2.0.3
## [82] plyr_1.8.7 ica_1.0-2 compiler_4.1.2
## [85] factoextra_1.0.7 fitdistrplus_1.1-6 cli_3.4.1
## [88] listenv_0.8.0 pbapply_1.5-0 MASS_7.3-54
## [91] mgcv_1.8-38 tidyselect_1.1.2 stringi_1.7.8
## [94] yaml_2.3.5 survMisc_0.5.5 sass_0.4.2
## [97] fastmatch_1.1-3 tools_4.1.2 parallel_4.1.2
## [100] circlize_0.4.15 rstudioapi_0.13 uuid_1.1-0
## [103] foreach_1.5.2 gridExtra_2.3 prodlim_2019.11.13
## [106] farver_2.1.1 Rtsne_0.16 digest_0.6.29
## [109] shiny_1.7.2 lava_1.6.10 quadprog_1.5-8
## [112] fpc_2.2-9 Rcpp_1.0.9 gridtext_0.1.4
## [115] car_3.1-0 broom_1.0.1 later_1.3.0
## [118] RcppAnnoy_0.0.19 httr_1.4.2 kernlab_0.9-31
## [121] colorspace_2.0-3 rvest_1.0.2 fs_1.5.2
## [124] tensor_1.5 reticulate_1.24 splines_4.1.2
## [127] uwot_0.1.11 yulab.utils_0.0.4 spatstat.utils_2.3-0
## [130] graphlayouts_0.8.2 shinythemes_1.2.0 flexmix_2.3-18
## [133] plotly_4.10.0 xtable_1.8-4 jsonlite_1.8.0
## [136] tidygraph_1.2.2 timeDate_3043.102 UpSetR_1.4.0
## [139] modeltools_0.2-23 ipred_0.9-12 R6_2.5.1
## [142] pillar_1.8.1 htmltools_0.5.3 mime_0.12
## [145] glue_1.6.1 fastmap_1.1.0 class_7.3-19
## [148] codetools_0.2-18 mvtnorm_1.1-3 utf8_1.2.2
## [151] bslib_0.4.0 spatstat.sparse_2.1-0 ggbeeswarm_0.6.0
## [154] leiden_0.3.9 rmarkdown_2.17 munsell_0.5.0
## [157] iterators_1.0.14 haven_2.4.3 reshape2_1.4.4
## [160] gtable_0.3.1 spatstat.core_2.4-0