vignettes/PhiSpace.Rmd
PhiSpace.RmdPhiSpace annotates the cells of a query dataset with the phenotypes of an annotated reference, such as cell types or differentiation stages. It fits a partial least squares (PLS) regression of the reference phenotypes on gene expression, and then predicts a score for every phenotype in every query cell. These PhiSpace scores are continuous: a cell can score highly for two related phenotypes, which describes a transitional state better than a single label. The scores of a cell, one per phenotype, are its coordinates in the phenotype space.
The reference and the query can come from different technologies. In this vignette the reference is bulk microarray data of sorted mouse blood stem and progenitor cells from the Immunological Genome Project (ImmGen), and the query is single-cell RNA-seq of the same compartment (Nestorowa et al., 2016). The query cells were sorted with fluorescence-activated cell sorting (FACS) into known gates, so we can compare the PhiSpace scores with the FACS labels.
The method is described in Mao et al. (2025), Genome Biology
26, 323 (https://doi.org/10.1186/s13059-025-03755-8). The spatial
vignette, vignette("PhiSpaceST", package = "PhiSpace"),
applies PhiSpace to spatial transcriptomics. Larger case studies are on
the PhiSpace
website.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("PhiSpace")The ImmGen reference in celldex
contains 830 bulk microarray samples of sorted mouse immune cells. We
keep the 36 samples of stem and progenitor cells. Their labels
(label.fine) name the sorted populations, for example
LTHSC and SC.LT34F (long-term haematopoietic
stem cells), SC.MEP (megakaryocyte-erythroid progenitors),
SC.CMP.DR (common myeloid progenitors), GMP
(granulocyte-monocyte progenitors), MLP (multi-lymphoid
progenitors) and proB.CLP (common lymphoid
progenitors).
immgen <- fetchReference("immgen", "2024-02-26")
immgen <- immgen[, immgen$label.main == "Stem cells"]
reference <- SingleCellExperiment(
assays = list(logcounts = as.matrix(assay(immgen, "logcounts"))),
colData = colData(immgen)
)
reference$stage <- sub("^Stem cells \\((.*)\\)$", "\\1", reference$label.fine)
table(reference$stage)
#>
#> GMP LTHSC MLP proB.CLP SC.CD150-CD48-
#> 3 3 7 1 1
#> SC.CDP SC.CMP.DR SC.LT34F SC.MDP SC.MEP
#> 3 2 3 3 2
#> SC.MPP34F SC.ST34F SC.STSL
#> 2 2 4The Nestorowa dataset in scRNAseq
contains 1,920 mouse stem and progenitor cells sequenced with
Smart-seq2. The column broad holds the FACS gate of each
cell; we keep the cells that have one. The genes are Ensembl
identifiers, so we map them to the gene symbols of the reference with
org.Mm.eg.db.
query <- fetchDataset("nestorowa-hsc-2016", "2024-04-18")
query <- query[, !is.na(query$broad)]
symbols <- AnnotationDbi::mapIds(
org.Mm.eg.db::org.Mm.eg.db, keys = rownames(query),
column = "SYMBOL", keytype = "ENSEMBL"
)
keep <- !is.na(symbols) & !duplicated(symbols)
query <- query[keep, ]
rownames(query) <- symbols[keep]
table(query$broad)
#>
#> CMP GMP LMPP LTHSC MEP MPP
#> 328 123 280 323 362 368PhiSpace needs the reference and the query to be normalised in the
same way. Here they are measured by different technologies, microarray
and single-cell RNA-seq, so we rank the genes within each sample with
RankTransf(). Ranks remove most differences of scale
between the two technologies. We first keep the genes measured in both
datasets.
genes <- intersect(rownames(reference), rownames(query))
length(genes)
#> [1] 18482
reference <- RankTransf(reference[genes, ], "logcounts", sparse = FALSE)
query <- logNormTransf(query[genes, ])
query <- RankTransf(query, "logcounts")PhiSpace() fits the PLS model on the reference and
stores the scores of the query cells in
reducedDim(query, "PhiSpace"). The argument
phenotypes names the reference column to predict, and
refAssay the assay to use (the query uses the same assay by
default). The number of PLS components, ncomp, controls the
complexity of the model; a good default is the number of phenotype
labels, here 13.
query <- PhiSpace(
reference, query,
phenotypes = "stage", refAssay = "rank", ncomp = 13
)
scores <- normPhiScores(reducedDim(query, "PhiSpace"))
dim(scores)
#> [1] 1784 13normPhiScores() centres the scores of each phenotype at
their median and scales them to a maximum absolute value of 1, so that
the scores of different phenotypes are comparable.
Each query cell has a score for each of the 13 reference stages. The heatmap shows the scores of all cells, grouped by their FACS gate.
gates <- c("LTHSC", "MPP", "LMPP", "CMP", "GMP", "MEP")
plotPhiSpaceHeatMap(
scores,
queryLabs = query$broad,
queryLvls = gates,
cluster_columns = TRUE,
name = "Score"
)
The mean score of each stage within each FACS gate summarises the heatmap.
gateMeans <- vapply(
gates,
function(g) colMeans(scores[query$broad == g, , drop = FALSE]),
numeric(ncol(scores))
)
round(t(gateMeans), 2)
#> SC.STSL SC.CD150-CD48- proB.CLP SC.LT34F SC.MDP SC.MEP SC.MPP34F SC.ST34F
#> LTHSC 0.21 0.08 -0.02 0.34 -0.06 0.02 0.03 -0.14
#> MPP 0.06 0.04 -0.03 0.07 -0.01 -0.03 0.24 0.03
#> LMPP 0.08 0.02 0.10 0.06 -0.04 -0.05 0.41 -0.15
#> CMP -0.14 -0.03 0.03 -0.11 0.13 0.04 -0.09 0.12
#> GMP -0.04 -0.04 0.05 -0.18 0.33 0.00 -0.26 0.12
#> MEP -0.11 -0.04 0.03 -0.17 0.01 0.46 -0.28 0.02
#> SC.CDP SC.CMP.DR GMP MLP LTHSC
#> LTHSC 0.09 -0.06 0.16 -0.32 0.14
#> MPP 0.01 0.01 0.01 -0.06 -0.03
#> LMPP 0.03 0.03 -0.07 -0.03 -0.18
#> CMP -0.02 0.12 0.03 0.09 0.08
#> GMP -0.08 0.05 0.17 0.10 0.02
#> MEP -0.03 -0.09 -0.13 0.07 0.02On average, cells sorted as MEP score highest for the ImmGen MEP
sample (SC.MEP), and cells sorted as long-term stem cells
(LTHSC) score highest for the long-term stem cell sample
SC.LT34F. The MPP and LMPP gates score highest for the
multipotent progenitor sample SC.MPP34F. The GMP gate
scores highest for SC.MDP (macrophage-dendritic cell
progenitors) and next for GMP; both are myeloid
progenitors. The CMP gate has no single top stage: its mean scores for
SC.MDP, SC.CMP.DR and SC.ST34F
are similar.
A discrete label for each cell is the stage with the highest score. The table compares these labels with the FACS gates.
topStage <- colnames(scores)[max.col(scores, ties.method = "first")]
table(FACS = query$broad, PhiSpace = topStage)
#> PhiSpace
#> FACS GMP LTHSC MLP proB.CLP SC.CD150-CD48- SC.CDP SC.CMP.DR SC.LT34F SC.MDP
#> CMP 13 30 26 12 22 2 66 2 61
#> GMP 17 7 1 0 9 0 12 1 52
#> LMPP 1 4 9 13 23 1 19 10 1
#> LTHSC 10 31 3 16 22 5 11 124 0
#> MEP 1 16 9 3 11 0 63 3 1
#> MPP 13 21 11 10 54 2 27 17 10
#> PhiSpace
#> FACS SC.MEP SC.MPP34F SC.ST34F SC.STSL
#> CMP 21 21 48 4
#> GMP 4 3 14 3
#> LMPP 1 167 9 22
#> LTHSC 28 12 12 49
#> MEP 239 2 11 3
#> MPP 2 113 47 41Discrete labels hide gradual differences between stages. For example, the 368 cells of the MPP gate have four common top stages, all early ones: SC.MPP34F (113 cells), SC.CD150-CD48- (54 cells), SC.ST34F (47 cells), SC.STSL (41 cells). The scores keep this information, because each cell has a score for every stage.
A principal component analysis (PCA) of the scores shows the cells in
the phenotype space. getPC() returns the principal
component scores and the proportion of variance that each component
explains.
pc <- getPC(scores, ncomp = 3)
round(pc$props, 2)
#> [1] 0.31 0.16 0.12
pcDat <- data.frame(pc$scores, gate = factor(query$broad, levels = gates))
ggplot2::ggplot(pcDat, ggplot2::aes(comp1, comp2, colour = gate)) +
ggplot2::geom_point(size = 0.8) +
ggplot2::labs(x = "PC1", y = "PC2", colour = "FACS gate") +
ggplot2::theme_bw()
The first principal component orders the cells from the long-term stem cell gate to the MEP gate, and the second separates the LMPP gate from the others. The FACS gates overlap in the phenotype space, as expected for cells that differentiate gradually.
sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 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.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] scRNAseq_2.24.0 celldex_1.20.0
#> [3] BiocStyle_2.38.0 PhiSpace_0.99.0
#> [5] testthat_3.3.2 SpatialExperiment_1.20.0
#> [7] SingleCellExperiment_1.32.0 SummarizedExperiment_1.40.0
#> [9] Biobase_2.70.0 GenomicRanges_1.62.1
#> [11] Seqinfo_1.0.0 IRanges_2.44.0
#> [13] S4Vectors_0.48.1 BiocGenerics_0.56.0
#> [15] generics_0.1.4 MatrixGenerics_1.22.0
#> [17] matrixStats_1.5.0
#>
#> loaded via a namespace (and not attached):
#> [1] BiocIO_1.20.0 bitops_1.1-0
#> [3] filelock_1.0.3 tibble_3.3.1
#> [5] XML_3.99-0.25 lifecycle_1.0.5
#> [7] httr2_1.3.0 rstatix_1.1.0
#> [9] edgeR_4.8.2 doParallel_1.0.17
#> [11] rprojroot_2.1.1 ensembldb_2.34.0
#> [13] lattice_0.22-6 alabaster.base_1.10.0
#> [15] backports_1.5.1 magrittr_2.0.5
#> [17] limma_3.66.0 sass_0.4.10
#> [19] rmarkdown_2.32 jquerylib_0.1.4
#> [21] yaml_2.3.12 metapod_1.18.0
#> [23] otel_0.2.0 askpass_1.2.1
#> [25] pkgbuild_1.4.8 reticulate_1.47.0
#> [27] DBI_1.3.0 RColorBrewer_1.1-3
#> [29] abind_1.4-8 pkgload_1.5.3
#> [31] purrr_1.2.2 downlit_0.4.5
#> [33] AnnotationFilter_1.34.0 RCurl_1.98-1.20
#> [35] rappdirs_0.3.4 circlize_0.4.18
#> [37] irlba_2.4.1 alabaster.sce_1.10.0
#> [39] umap_0.2.10.0 RSpectra_0.16-2
#> [41] dqrng_0.4.1 pkgdown_2.2.0
#> [43] DelayedMatrixStats_1.32.0 codetools_0.2-20
#> [45] DelayedArray_0.36.1 scuttle_1.20.0
#> [47] xml2_1.6.0 tidyselect_1.2.1
#> [49] shape_1.4.6.1 UCSC.utils_1.6.1
#> [51] farver_2.1.2 ScaledMatrix_1.18.0
#> [53] BiocFileCache_3.0.0 GenomicAlignments_1.46.0
#> [55] jsonlite_2.0.0 GetoptLong_1.1.1
#> [57] BiocNeighbors_2.4.0 Formula_1.2-6
#> [59] iterators_1.0.14 systemfonts_1.3.2
#> [61] foreach_1.5.2 tools_4.5.0
#> [63] ragg_1.5.2 Rcpp_1.1.2
#> [65] glue_1.8.1 SparseArray_1.10.10
#> [67] xfun_0.61 GenomeInfoDb_1.46.2
#> [69] dplyr_1.2.1 HDF5Array_1.38.0
#> [71] gypsum_1.6.0 withr_3.0.3
#> [73] BiocManager_1.30.27 fastmap_1.2.0
#> [75] rhdf5filters_1.22.0 bluster_1.20.0
#> [77] openssl_2.4.2 digest_0.6.39
#> [79] rsvd_1.0.5 R6_2.6.1
#> [81] textshaping_1.0.5 colorspace_2.1-4
#> [83] RSQLite_3.53.3 cigarillo_1.0.0
#> [85] h5mread_1.2.1 tidyr_1.3.2
#> [87] rtracklayer_1.70.1 FNN_1.1.4.1
#> [89] httr_1.4.9 htmlwidgets_1.6.4
#> [91] S4Arrays_1.10.1 org.Mm.eg.db_3.22.0
#> [93] whisker_0.4.1 pkgconfig_2.0.3
#> [95] gtable_0.3.6 blob_1.3.0
#> [97] ComplexHeatmap_2.26.1 S7_0.2.2
#> [99] XVector_0.50.0 brio_1.1.5
#> [101] htmltools_0.5.9 carData_3.0-6
#> [103] bookdown_0.48 ProtGenerics_1.42.0
#> [105] clue_0.3-68 scales_1.4.0
#> [107] alabaster.matrix_1.10.0 png_0.1-9
#> [109] scran_1.38.1 knitr_1.52
#> [111] rstudioapi_0.19.0 rjson_0.2.23
#> [113] curl_8.0.0 cachem_1.1.0
#> [115] rhdf5_2.54.1 GlobalOptions_0.1.4
#> [117] BiocVersion_3.22.0 parallel_4.5.0
#> [119] AnnotationDbi_1.72.0 restfulr_0.0.17
#> [121] desc_1.4.3 pillar_1.11.1
#> [123] grid_4.5.0 alabaster.schemas_1.10.0
#> [125] vctrs_0.7.3 ggpubr_1.0.0
#> [127] car_3.1-5 BiocSingular_1.26.1
#> [129] dbplyr_2.6.0 beachmat_2.26.0
#> [131] cluster_2.1.8.1 evaluate_1.0.5
#> [133] GenomicFeatures_1.62.0 magick_2.9.1
#> [135] Rsamtools_2.26.0 cli_3.6.6
#> [137] locfit_1.5-9.12 compiler_4.5.0
#> [139] rlang_1.3.0 crayon_1.5.3
#> [141] kerndwd_2.0.3 ggsignif_0.6.4
#> [143] labeling_0.4.3 plyr_1.8.9
#> [145] fs_2.1.0 alabaster.se_1.10.0
#> [147] BiocParallel_1.44.0 Biostrings_2.78.0
#> [149] lazyeval_0.2.3 Matrix_1.7-3
#> [151] ExperimentHub_3.0.0 sparseMatrixStats_1.22.0
#> [153] bit64_4.8.6 ggplot2_4.0.3
#> [155] Rhdf5lib_1.32.0 KEGGREST_1.50.0
#> [157] statmod_1.5.2 alabaster.ranges_1.10.0
#> [159] AnnotationHub_4.0.0 igraph_2.3.4
#> [161] broom_1.0.13 memoise_2.0.1
#> [163] bslib_0.12.0 bit_4.6.0