---
title: "Introduction to PhiSpace"
subtitle: "Annotating single cells with a bulk reference"
author: "Jiadong Mao"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
vignette: >
  %\VignetteIndexEntry{Introduction to PhiSpace}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5
)
```

# Introduction

PhiSpace 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](https://mao-lab-org.github.io/PhiSpace/).

# Installation

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("PhiSpace")
```

```{r packages, message = FALSE}
library(PhiSpace)
library(SingleCellExperiment)
library(celldex)
library(scRNAseq)
```

# Data

## Reference: ImmGen stem and progenitor cells

The ImmGen reference in `r Biocpkg("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).

```{r reference, message = FALSE}
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)
```

## Query: single-cell RNA-seq of stem and progenitor cells

The Nestorowa dataset in `r Biocpkg("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
`r Biocpkg("org.Mm.eg.db")`.

```{r query, message = FALSE}
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)
```

## Normalisation

PhiSpace 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.

```{r normalise}
genes <- intersect(rownames(reference), rownames(query))
length(genes)
reference <- RankTransf(reference[genes, ], "logcounts", sparse = FALSE)
query <- logNormTransf(query[genes, ])
query <- RankTransf(query, "logcounts")
```

# Annotation

`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.

```{r phispace}
query <- PhiSpace(
    reference, query,
    phenotypes = "stage", refAssay = "rank", ncomp = 13
)
scores <- normPhiScores(reducedDim(query, "PhiSpace"))
dim(scores)
```

`normPhiScores()` 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.

# Results

## Scores by FACS gate

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.

```{r heatmap, fig.height = 6}
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.

```{r means}
gateMeans <- vapply(
    gates,
    function(g) colMeans(scores[query$broad == g, , drop = FALSE]),
    numeric(ncol(scores))
)
round(t(gateMeans), 2)
```

On 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.

## Discrete labels

A discrete label for each cell is the stage with the highest score. The table
compares these labels with the FACS gates.

```{r labels}
topStage <- colnames(scores)[max.col(scores, ties.method = "first")]
table(FACS = query$broad, PhiSpace = topStage)
```

```{r mppTop, include = FALSE}
mppTop <- sort(table(topStage[query$broad == "MPP"]), decreasing = TRUE)
```

Discrete labels hide gradual differences between stages. For example, the
`r sum(query$broad == "MPP")` cells of the MPP gate have four common top
stages, all early ones:
`r paste0(names(mppTop)[1:4], " (", mppTop[1:4], " cells)", collapse = ", ")`.
The scores keep this information, because each cell has a score for every
stage.

## The phenotype space

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.

```{r pca}
pc <- getPC(scores, ncomp = 3)
round(pc$props, 2)
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.

# Session information

```{r session}
sessionInfo()
```