07-classification.qmd

# Finding associations between clinical variables and spatial features

When it comes to biological datasets, the end goal is either mechanistic or translational. For example, if we had a mechanistic end goal, we might want to find what genes are differentially expressed between two conditions, and further aim to characterise the pathways which lead to this differential expression. Alternatively, if the end goal is translational, we might want to use a biological dataset that can be relatively cheaply obtained (e.g. IMC) to predict whether a patient's disease will progress or not progress (e.g. metastasize in cancer).

<!-- Steps: -->

<!-- 1.  Introduction to ClassifyR -->
<!-- 2.  Classification of patients by condition -->
<!-- 3.  Classification of patients by survival -->
<!-- 4.  Easy and Hard to classify patients (samplesMetricMap) -->
<!-- 5.  Maximising accuracy during classification (parameter tuning for crossValidate) -->

```{r load libraries, echo=FALSE, results="hide", warning=FALSE}
suppressPackageStartupMessages({
  library(ClassifyR)
  library(lisaClust)
  library(Statial)
  library(spicyR)
  library(ggplot2)
})
```

```{r, eval = FALSE}
library(ClassifyR)
library(lisaClust)
library(Statial)
library(spicyR)
library(ggplot2)
```

```{r, setParam}
# set parameters
set.seed(51773)
use_mc <- TRUE
if (use_mc) {
  nCores <- max(parallel::detectCores()/2, 1)
} else {
  nCores <- 1
}
BPPARAM <- simpleSeg:::generateBPParam(nCores)
theme_set(theme_classic())
```

## Introduction to `ClassifyR`

[ClassifyR](https://www.bioconductor.org/packages/release/bioc/html/ClassifyR.html) provides a structured pipeline for cross-validated classification. Classification is viewed in terms of four stages: data transformation, feature selection, classifier training, and prediction. The driver functions crossValidate and runTests implements varieties of cross-validation. They are:

-   Permutation of the order of samples followed by k-fold cross-validation
-   Repeated $x\%$ test set cross-validation
-   leave-k-out cross-validation

We will use the Keren 2018 dataset to perform two predictions: 1) predict the patient's tumour type (cold, compartmentalised, or mixed) and 2) predict the patient's survival outcome.

```{r, message=FALSE,warning=FALSE}
# load the Keren 2018 dataset
kerenSPE <- SpatialDatasets::spe_Keren_2018()

# remove any missing data in our outcome columns
kerenSPE = kerenSPE[, complete.cases(colData(kerenSPE)[, c("Censored", "Survival_days_capped*",
                                                           "tumour_type")])]

kontextMat <- readRDS("data/kontextMat.rds")

stateMat <- readRDS("data/stateMat.rds")
```

## Classification of patients by condition

We will use the Keren 2018 dataset to classify a patient's tumour into one of three types: cold, compartmentalised, and mixed.

We will build a list of feature matrices using the features we've generated in the previous chapters:

1.  Cell type proportions (FuseSOM)
2.  Cell type co-localistion using the L-function (spicyR)
3.  Cell type co-localisation with respect to a parent (Kontextual)
4.  Spatial domains (lisaClust)
5.  Marker means in each cell type (Statial)
6.  Marker means in each cell type in each region (Statial)
7.  Proximity-associated changes in marker expression (SpatioMark)

```{r}
kontextMat <- readRDS("data/kontextMat.rds")
stateMat <- readRDS("data/stateMat.rds")
```

```{r}
data <- list()

# Cell type proportions (FuseSOM)
data[["proportions"]] <- getProp(kerenSPE, "cellType")

# Cell localisation (spicyR)
data[["spicyR"]] <- getPairwise(kerenSPE,
                                BPPARAM = BPPARAM)

# Cell localisation with respect to a parent (Kontextual)
data[["Kontextual"]] <- kontextMat

# Spatial Domains (lisaClust)
data[["lisaClust"]] <- getProp(lisaClust(kerenSPE, k = 5), "region")

# Marker means in each cell type (Statial)
data[["celltypeMarkerMeans"]] <- getMarkerMeans(kerenSPE, imageID = "imageID",
                                        cellType = "cellType",
                                        region = "cellType")

# Marker means in each cell type in each region (Statial)
data[["regionMarkerMeans"]] <- getMarkerMeans(lisaClust(kerenSPE, k = 5), 
                                        imageID = "imageID",
                                        cellType = "cellType",
                                        region = "region")

# Proximity-associated changes in marker expression (SpatioMark)
data[["SpatioMark"]] <- stateMat
```

We will then generate a factor vector of our outcome variable.

```{r}
# outcome vector
outcome = kerenSPE$tumour_type[!duplicated(kerenSPE$imageID)]
names(outcome) = kerenSPE$imageID[!duplicated(kerenSPE$imageID)]

head(outcome, 5)
```

`ClassifyR` provides a convenient function, `crossValidate`, to build and test models. `crossValidate` must be supplied with `measurements`, a simple tabular data container or a list-like structure of such related tabular data on common samples. It can be in the form of a `matrix`, `data.frame`, `DataFrame`, `MultiAssayExperiment` or a list of `data.frames`.

`crossValidate` must also be supplied with `outcome`, which represents the prediction to be made. `outcome` can be either a factor containing the class labels for each observation, or a `character` of length 1 that matches a column name in a data frame which holds the classes. If a `character` is provided, `crossValidate` will automatically remove the classes before training.

By default, `crossValidate` will build and train a random forest. Alternative classifiers can be specified using the `classifier` argument. To view all available feature selection and classification approaches, use the `available()` function.

```{r}
# perform 50 repeats of 5-fold cross-validation
cv = crossValidate(measurements = data,
                   outcome = outcome,
                   nFolds = 5,
                   nRepeats = 50,
                   nCores = 10)
```

We can use `performancePlot` to visualise performance metrics for all our features. Here, we visualise the AUC for each of the six feature matrices we tested. Additional performance metrics can be specified in the `metric` argument.

```{r fig.width = 12, fig.height = 6}
performancePlot(
  cv,
  metric = "AUC",
  characteristicsList = list(x = "Assay Name"),
  orderingList = list("Assay Name" = c("proportions", "spicyR", "lisaClust", "Kontextual", "celltypeMarkerMeans", "regionMarkerMeans", "SpatioMark"))
)
```

From the graph, both `lisaClust` and `proportions` appear to capture information which is predictive of the tumour type of patients.

## Classification of patients by survival

`crossValidate` also has the capacity to test classification performance for a survival outcome. In this case, `outcome` must be a `Surv` object of the same length as the number of samples in the feature matrix and should information about the time and censoring of the samples. Alternatively, we can specify `outcome` to be a `character` of length 2 or 3 that each match a column name in a data frame which holds information about the time and censoring of the samples. The time-to-event columns will automatically be removed before training is done.

We will first add a `survival` column to the `kerenSPE` object, and remove any rows that contain missing survival data.

```{r}
# create a Surv object named "survival"
kerenSPE$event = 1 - kerenSPE$Censored
kerenSPE$survival = Surv(kerenSPE$`Survival_days_capped*`, kerenSPE$event)

# outcome vector
surv_outcome = kerenSPE$survival[!duplicated(kerenSPE$imageID)]
names(surv_outcome) = kerenSPE$imageID[!duplicated(kerenSPE$imageID)]

surv_outcome
```

Next, we filter `data` to ensure the image IDs it contains are consistent with those in `surv_outcome`.

We can then run `crossValidate` and specify the outcome to be `surv_outcome`, and use `performancePlot` to visualise the performance of the cross-validation. Since we are performing survival analysis, we will specify `metric = "C-index"`.

```{r fig.width = 12, fig.height = 6}
# perform 50 repeats of 5-fold cross-validation
surv_cv = crossValidate(measurements = data,
                   outcome = surv_outcome,
                   nFolds = 5,
                   nRepeats = 50,
                   nCores = 10)

performancePlot(surv_cv,
  metric = "C-index",
  characteristicsList = list(x = "Assay Name"),
  orderingList = list("Assay Name" = c("proportions", "spicyR", "lisaClust", "Kontextual", "celltypeMarkerMeans", "regionMarkerMeans", "SpatioMark"))
)
```

From the graph, we can see that `lisaClust` appears to capture information that is predictive of survival outcomes comparatively well.

## Easy and hard to classify patients

The `samplesMetricMap` function allows the visual comparison of sample-wise error rate or accuracy measures from the cross-validation process.

### Predicting tumour type

```{r, fig.width=15}
samplesMetricMap(cv,  
                 classColours = c("#3F48CC", "#880015", "#00b409"),
                 metricColours = list(c("#FFFFFF", "#CFD1F2", "#9FA3E5", "#6F75D8", "#3F48CC"),
                                      c("#FFFFFF", "#E1BFC4", "#C37F8A", "#A53F4F", "#880015"),
                                      c("#FFFFFF", "#c1ffc4", "#88ff8e", "#5aff62", "#00b409")))
```

### Predicting survival outcomes

```{r}
samplesMetricMap(surv_cv, showXtickLabels = FALSE)
```

## Maximising accuracy during classification