scGOclust_vignette
Load packages
scGOclust is a package that leverages Gene Ontology to
analyse the functional profile of cells with scRNA-seq data.
# load required libraries
library(Seurat)
#> Loading required package: SeuratObject
#> Loading required package: sp
#> 'SeuratObject' was built under R 4.6.0 but the current version is
#> 4.6.1; it is recomended that you reinstall 'SeuratObject' as the ABI
#> for R may have changed
#> 'SeuratObject' was built with package 'Matrix' 1.7.5 but the current
#> version is 1.7.6; it is recomended that you reinstall 'SeuratObject' as
#> the ABI for 'Matrix' may have changed
#>
#> Attaching package: 'SeuratObject'
#> The following objects are masked from 'package:base':
#>
#> intersect, t
library(pheatmap)
library(httr)
## if (!require("devtools")) install.packages("devtools")
## install latest from source
## for reproducibility we do not update dependencies
# devtools::install_github("Papatheodorou-Group/scGOclust", upgrade = FALSE)
library(scGOclust)2. Load input data
# get a gene to GO BP terms mapping table
# remove electronically inferred records
# sometimes ensembl complains about ssh certificate has expired, this is a known issue, run this code
httr::set_config(httr::config(ssl_verifypeer = FALSE))
#mmu_tbl = ensemblToGo(species = 'mmusculus', GO_linkage_type = c('experimental', 'phylogenetic', 'computational', 'author', 'curator' ))
#dme_tbl = ensemblToGo(species = 'dmelanogaster', GO_linkage_type = c('experimental', 'phylogenetic', 'computational', 'author', 'curator' ))
# here we load the example data for convenience
data(mmu_tbl)
data(dme_tbl)3. Build GO BP profile
## construct a Seurat object with GO BP as features
mmu_go_obj <- makeGOSeurat(ensembl_to_GO = mmu_tbl, feature_type = 'external_gene_name', seurat_obj = mmu_subset)
#> collect data
#> compute GO to cell matrix, might take a few secs
#> Warning: Data is of class matrix. Coercing to dgCMatrix.
#> time used: 0.55 secs
#> removing 0 GO terms with all 0 in dataset
#> returning GO Seurat object
dme_go_obj <- makeGOSeurat(ensembl_to_GO = dme_tbl, feature_type = 'external_gene_name', seurat_obj = dme_subset)
#> collect data
#> compute GO to cell matrix, might take a few secs
#> Warning: Data is of class matrix. Coercing to dgCMatrix.
#> time used: 0.23 secs
#> removing 0 GO terms with all 0 in dataset
#> returning GO Seurat object4. Calculate cell type average GO BP profile
# specify the column with cell type annotation in [email protected]
mmu_ct_go <- getCellTypeGO(go_seurat_obj = mmu_go_obj, cell_type_col = 'cell_type_annotation')
#> perform normalization and log1p for mmu_go_obj
#> Normalizing layer: counts
#> Centering and scaling data matrix
#> As of Seurat v5, we recommend using AggregateExpression to perform pseudo-bulk analysis.
#> Names of identity class contain underscores ('_'), replacing with dashes ('-')
#> This message is displayed once per session.
dme_ct_go <- getCellTypeGO(go_seurat_obj = dme_go_obj, cell_type_col = 'annotation')
#> perform normalization and log1p for dme_go_obj
#> Normalizing layer: counts
#> Centering and scaling data matrix5. Calculate within-species cell type functional similariy
# heatmap of Pearson's correlation coefficient of cell type average BP profiles within species
mmu_corr = cellTypeGOCorr(cell_type_go = mmu_ct_go, corr_method = 'pearson')
pheatmap(mmu_corr)
5. Calculate cross-species cell type functional similariy
# calculate Pearson's correlation coefficient of cell type average BP profiles across species
corr = crossSpeciesCellTypeGOCorr(species_1 = 'mmusculus', species_2 = 'dmelanogaster', cell_type_go_sp1 = mmu_ct_go, cell_type_go_sp2 = dme_ct_go, corr_method = 'pearson')
# tries to put cells with higher values on the diagonal
# helpful when cross-species cell type similarity signal is less clear
plotCellTypeCorrHeatmap(corr, width = 9, height = 10)
# scale per column or row to see the relative similarity
plotCellTypeCorrHeatmap(corr, scale = 'column', width = 9, height = 10)
6. Dimensional reduction and UMAP visualization of cells with GO profile
# analyze the cell-by-GO BP profile as a count matrix
# Note that the example data has very few cells (for reducing data size), so I am using a small UMAP min_dist and small spread here. Default min_dist is 0.3 and spread is 1.
mmu_go_analyzed = analyzeGOSeurat(go_seurat_obj = mmu_go_obj, cell_type_col = 'cell_type_annotation', min.dist = 0.1, spread = 0.5)
#> Warning in analyzeGOSeurat(go_seurat_obj = mmu_go_obj, cell_type_col = "cell_type_annotation", : 'seurat_clusters' is presented in the metadata column names, this will be overwritten.
#> If you would like to keep annotations in this column, rename the column.
#> perform normalization and log1p for mmu_go_obj
#> Normalizing layer: counts
#> Warning: The following arguments are not used: spread
#> Finding variable features for layer counts
#> Warning: The following arguments are not used: spread
#> Warning: The following arguments are not used: spread
#> Warning: The following arguments are not used: spread
#> Computing nearest neighbor graph
#> Computing SNN
#> Warning: The following arguments are not used: spread
#> Warning: The following arguments are not used: spread
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#>
#> Number of nodes: 219
#> Number of edges: 9890
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.4728
#> Number of communities: 4
#> Elapsed time: 0 seconds
#> Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
#> To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
#> This message will be shown once per session
#> 08:15:58 UMAP embedding parameters a = 5.069 b = 1.003
#> 08:15:58 Read 219 rows and found 50 numeric columns
#> 08:15:58 Using Annoy for neighbor search, n_neighbors = 30
#> 08:15:58 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 08:15:58 Writing NN index file to temp file /tmp/RtmpwmVOiO/file20a7276c85fe
#> 08:15:58 Searching Annoy index using 1 thread, search_k = 3000
#> 08:15:58 Annoy recall = 100%
#> 08:15:59 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
#> 08:16:00 Initializing from normalized Laplacian + noise (using RSpectra)
#> 08:16:00 Commencing optimization for 500 epochs, with 7820 positive edges
#> 08:16:00 Using rng type: pcg
#> 08:16:01 Optimization finished# UMAP plot of the analyzed cell-by-GO BP profile
# labeled by previously specified cell annotation column in meta.data
DimPlot(mmu_go_analyzed, label = TRUE) + NoLegend()
dme_go_analyzed = analyzeGOSeurat(go_seurat_obj = dme_go_obj, cell_type_col = 'annotation', min_dist=0.1, spread=0.5)
#> Warning in analyzeGOSeurat(go_seurat_obj = dme_go_obj, cell_type_col = "annotation", : 'seurat_clusters' is presented in the metadata column names, this will be overwritten.
#> If you would like to keep annotations in this column, rename the column.
#> perform normalization and log1p for dme_go_obj
#> Normalizing layer: counts
#> Warning: The following arguments are not used: min_dist, spread
#> Finding variable features for layer counts
#> Warning: The following arguments are not used: min_dist, spread
#> Warning: The following arguments are not used: min_dist, spread
#> Warning: The following arguments are not used: min_dist, spread
#> Computing nearest neighbor graph
#> Computing SNN
#> Warning: The following arguments are not used: min_dist, spread
#> Warning: The following arguments are not used: min_dist, spread
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#>
#> Number of nodes: 180
#> Number of edges: 6418
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.5467
#> Number of communities: 4
#> Elapsed time: 0 seconds
#> Warning: The following arguments are not used: min_dist
#> 08:16:02 UMAP embedding parameters a = 3.244 b = 1.447
#> 08:16:02 Read 180 rows and found 50 numeric columns
#> 08:16:02 Using Annoy for neighbor search, n_neighbors = 30
#> 08:16:02 Building Annoy index with metric = cosine, n_trees = 50
#> 0% 10 20 30 40 50 60 70 80 90 100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 08:16:02 Writing NN index file to temp file /tmp/RtmpwmVOiO/file20a722786043
#> 08:16:02 Searching Annoy index using 1 thread, search_k = 3000
#> 08:16:02 Annoy recall = 100%
#> 08:16:03 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
#> 08:16:04 Initializing from normalized Laplacian + noise (using RSpectra)
#> 08:16:04 Commencing optimization for 500 epochs, with 6610 positive edges
#> 08:16:04 Using rng type: pcg
#> 08:16:04 Optimization finished
7. Get co-up and co-down regulated terms between pairs of cell types
## calculation takes a few minutes due to the Wilcox signed rank test
ct_shared_go = scGOclust::getCellTypeSharedGO(species_1 = 'mmusculus', species_2 = 'dmelanogaster', analyzed_go_seurat_sp1 = mmu_go_analyzed, analyzed_go_seurat_sp2 = dme_go_analyzed, cell_type_col_sp1 = 'cell_type_annotation', cell_type_col_sp2 = 'annotation', layer_use = 'data')
# query shared GO terms for specific cell type pairs
getCellTypeSharedTerms(shared_go = ct_shared_go,
cell_type_sp1 = 'intestine_Enteroendocrine cell',
cell_type_sp2 = 'enteroendocrine cell',
return_full = FALSE)
# viola
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.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
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] future_1.76.0 scGOclust_0.2.1 httr_1.4.9 pheatmap_1.0.13
#> [5] Seurat_5.5.1 SeuratObject_5.4.0 sp_2.2-3 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_2.0.0
#> [4] magrittr_2.0.5 spatstat.utils_3.2-5 farver_2.1.2
#> [7] vctrs_0.7.3 ROCR_1.0-12 memoise_2.0.1
#> [10] spatstat.explore_3.8-3 BiocBaseUtils_1.15.1 progress_1.2.3
#> [13] htmltools_0.5.9 curl_8.0.0 pracma_2.4.6
#> [16] sass_0.4.10 sctransform_0.4.3 parallelly_1.48.0
#> [19] KernSmooth_2.23-27 bslib_0.12.0 htmlwidgets_1.6.4
#> [22] ica_1.0-3 httr2_1.3.0 plyr_1.8.9
#> [25] plotly_4.12.1 zoo_1.9-1 cachem_1.1.0
#> [28] networkD3_0.4.1 buildtools_1.0.0 igraph_2.3.4
#> [31] mime_0.13 lifecycle_1.0.5 pkgconfig_2.0.3
#> [34] Matrix_1.7-6 R6_2.6.1 fastmap_1.2.0
#> [37] fitdistrplus_1.2-6 shiny_1.14.0 digest_0.6.39
#> [40] S4Vectors_0.51.10 AnnotationDbi_1.75.2 patchwork_1.3.2
#> [43] tensor_1.5.1 RSpectra_0.16-2 irlba_2.3.7
#> [46] RSQLite_3.53.3 labeling_0.4.3 filelock_1.0.3
#> [49] progressr_1.0.0 spatstat.sparse_3.2-0 slanter_0.2-0
#> [52] polyclip_1.10-7 abind_1.4-8 compiler_4.6.1
#> [55] withr_3.0.3 bit64_4.8.6 S7_0.2.2
#> [58] DBI_1.3.0 fastDummies_1.7.6 biomaRt_2.69.4
#> [61] MASS_7.3-66 tools_4.6.1 lmtest_0.9-40
#> [64] otel_0.2.0 httpuv_1.6.17 future.apply_1.20.2
#> [67] goftest_1.2-3 glue_1.8.1 nlme_3.1-171
#> [70] promises_1.5.0 grid_4.6.1 Rtsne_0.17
#> [73] cluster_2.1.8.3 reshape2_1.4.5 generics_0.1.4
#> [76] gtable_0.3.6 spatstat.data_3.1-9 tidyr_1.3.2
#> [79] hms_1.1.4 data.table_1.18.6.1 XVector_0.53.0
#> [82] BiocGenerics_0.59.12 spatstat.geom_3.8-3 RcppAnnoy_0.0.23
#> [85] ggrepel_0.9.8 RANN_2.6.3 pillar_1.11.1
#> [88] stringr_1.6.0 spam_2.11-4 RcppHNSW_0.7.0
#> [91] limma_3.99.0 later_1.4.8 splines_4.6.1
#> [94] dplyr_1.2.1 BiocFileCache_3.3.0 lattice_0.23-1
#> [97] survival_3.8-12 bit_4.6.0 deldir_2.0-4
#> [100] tidyselect_1.2.1 Biostrings_2.81.9 maketools_1.3.2
#> [103] miniUI_0.1.2 pbapply_1.7-5 knitr_1.52
#> [106] gridExtra_2.3.1 Seqinfo_1.3.2 IRanges_2.47.5
#> [109] scattermore_1.2 stats4_4.6.1 xfun_0.61
#> [112] Biobase_2.73.2 statmod_1.5.2 matrixStats_1.5.0
#> [115] stringi_1.8.9 yaml_2.3.12 evaluate_1.0.5
#> [118] codetools_0.2-20 data.tree_1.2.0 tibble_3.3.1
#> [121] cli_3.6.6 uwot_0.2.5 xtable_1.8-8
#> [124] reticulate_1.47.0 jquerylib_0.1.4 Rcpp_1.1.2
#> [127] globals_0.19.1 spatstat.random_3.5-2 dbplyr_2.6.0
#> [130] png_0.1-9 spatstat.univar_3.2-0 parallel_4.6.1
#> [133] ggplot2_4.0.3 blob_1.3.0 prettyunits_1.2.0
#> [136] dotCall64_1.2 listenv_1.0.0 viridisLite_0.4.3
#> [139] scales_1.4.0 ggridges_0.5.7 crayon_1.5.3
#> [142] purrr_1.2.2 rlang_1.3.0 KEGGREST_1.53.6
#> [145] cowplot_1.2.0- Load packages
- 2. Load input data
- 3. Build GO BP profile
- 4. Calculate cell type average GO BP profile
- 5. Calculate within-species cell type functional similariy
- 5. Calculate cross-species cell type functional similariy
- 6. Dimensional reduction and UMAP visualization of cells with GO profile
- 7. Get co-up and co-down regulated terms between pairs of cell types