12 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).
# 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())
12.1 Introduction to ClassifyR
ClassifyR 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 (compartmentalised or mixed) and 2) predict the patient’s survival outcome.
# 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")])]
# load pre-computed data for efficiency
kontextMat <- readRDS("data/kontextMat.rds")
stateMat <- readRDS("data/stateMat.rds")12.2 Classification of patients by condition
We will use the Keren 2018 dataset to classify a patient’s tumour into one of two types: compartmentalised and mixed. First, we will filter out all cases which had cold tumours.
# filter out cold tumours
kerenSPE = kerenSPE[, colData(kerenSPE)$tumour_type != "cold"]
kerenSPE$tumour_type = droplevels(kerenSPE$tumour_type)
levels(kerenSPE$tumour_type) = c("compartmentalised", "mixed")
kontextMat = kontextMat[rownames(kontextMat) %in% unique(kerenSPE$imageID), ]
stateMat = stateMat[rownames(stateMat) %in% unique(kerenSPE$imageID), ]We will build a list of feature matrices using the features we’ve generated in the previous chapters:
- Cell type proportions (FuseSOM)
- Co-localisation between pairs of cell types using the L-function (spicyR)
- Cell type co-localisation with respect to a parent population using
Kontextual(Statial) - Regions of co-localisation, or spatial domains (lisaClust)
- Marker means in each cell type (Statial)
- Marker means in each cell type in each region (Statial)
- Proximity-associated changes in marker expression using SpatioMark (Statial)
data <- list()
# Cell type proportions (FuseSOM)
data[["proportions"]] <- getProp(kerenSPE, "cellType")
# Cell localisation (spicyR)
data[["spicyR"]] <- getPairwise(kerenSPE,
BPPARAM = BPPARAM) |> as.data.frame()
# Cell localisation with respect to a parent (Kontextual)
data[["Kontextual"]] <- kontextMat
# Spatial Domains (lisaClust)
data[["lisaClust"]] <- getProp(lisaClust(kerenSPE, k = 5,
BPPARAM = BPPARAM), "region")Generating local L-curves.# 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,
BPPARAM = BPPARAM),
imageID = "imageID",
cellType = "cellType",
region = "region")Generating local L-curves.# Proximity-associated changes in marker expression (SpatioMark)
data[["SpatioMark"]] <- stateMatTime for this code chunk to run with 5.5 cores: 33.54 seconds
We will then generate a factor vector of our outcome variable.
# outcome vector
outcome = kerenSPE$tumour_type[!duplicated(kerenSPE$imageID)]
names(outcome) = kerenSPE$imageID[!duplicated(kerenSPE$imageID)]
head(outcome, 5) 1 2 3 4
mixed mixed compartmentalised compartmentalised
5
compartmentalised
Levels: compartmentalised mixedClassifyR 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 measurements 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.
# perform 50 repeats of 5-fold cross-validation
cv = crossValidate(measurements = data,
outcome = outcome,
nFolds = 5,
nRepeats = 50,
nCores = nCores)Time for this code chunk to run with 5.5 cores: 25.54 seconds
We can use performancePlot to visualise performance metrics for all our features. Here, we visualise the AUC for each of the seven feature matrices we tested. Additional performance metrics can be specified in the metric argument.
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.
12.3 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 contain 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.
# 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 1 2 3 4 5 6 7 8 9 10 11 12 13
2612 745 3130+ 2523+ 1683+ 2275+ 584 946 3767+ 3822+ 3774+ 4353+ 1072
14 16 17 18 20 21 23 27 28 29 31 32 33
4145+ 530 2842 5063+ 4761+ 635 91 3658 3767+ 1319 1009 1568+ 1738+
34 35 36 37 39 40 41
2832+ 2759+ 3063+ 2853+ 2096+ 3573 3355+ 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".
# perform 50 repeats of 5-fold cross-validation
surv_cv = crossValidate(measurements = data,
outcome = surv_outcome,
nFolds = 5,
nRepeats = 50,
nCores = nCores)
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.
12.4 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.
12.4.1 Predicting tumour type
samplesMetricMap(cv,
classColours = c("#3F48CC", "#880015"),
metricColours = list(c("#FFFFFF", "#CFD1F2", "#9FA3E5", "#6F75D8", "#3F48CC"),
c("#FFFFFF", "#E1BFC4", "#C37F8A", "#A53F4F", "#880015")))
TableGrob (2 x 1) "arrange": 2 grobs
z cells name grob
1 1 (2-2,1-1) arrange gtable[layout]
2 2 (1-1,1-1) arrange text[GRID.text.409]The benefit of this plot is that it allows the easy identification of samples which are hard to classify and could be explained by considering additional information about them. For example, we can see that patients 36 and 13 in particular were difficult to classify. Overall, both cell type proportions (FuseSOM) and cell type co-localisation (spicyR and Kontextual) performed comparitively well at classifying patients, and compartmentalised tumours were easier to classify compared to mixed tumours.
12.4.2 Predicting survival outcomes
We can also use sampleMetricMap to identify samples that were difficult to classify with respect to a survival outcome.
samplesMetricMap(surv_cv)
TableGrob (2 x 1) "arrange": 2 grobs
z cells name grob
1 1 (2-2,1-1) arrange gtable[layout]
2 2 (1-1,1-1) arrange text[GRID.text.561]We can see that there was some difficulty in classying patient 35, and that overall, features showed similar performance across samples. lisaClust shows the overall highest performance, especially for samples 31-20.
DataFrame(), outcome, crossValParams, SVMparams) ``` ::: The index of chosen of the parameters, as well as all combinations of parameters and their associated performance metric, are stored for every validation, and can be accessed with the `tunedParameters` function. ::: {.cell} ```{.r .cell-code} tunedParameters(SVMresults)[1:5] ``` ::: -->In this way, we can predict clinical outcomes and identify samples that are hard to classify.
In subsequent sections, we will demonstrate how the end-to-end workflow can be applied to a single dataset in the form of case studies.
13 sessionInfox
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Australia/Sydney
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] SpatialDatasets_1.4.0 SpatialExperiment_1.16.0
[3] SingleCellExperiment_1.28.1 ExperimentHub_2.14.0
[5] AnnotationHub_3.14.0 BiocFileCache_2.14.0
[7] dbplyr_2.5.0 ggplot2_3.5.1
[9] spicyR_1.18.0 Statial_1.8.0
[11] lisaClust_1.14.4 ClassifyR_3.10.5
[13] survival_3.7-0 BiocParallel_1.40.0
[15] MultiAssayExperiment_1.32.0 SummarizedExperiment_1.36.0
[17] Biobase_2.66.0 GenomicRanges_1.58.0
[19] GenomeInfoDb_1.42.0 IRanges_2.40.0
[21] MatrixGenerics_1.18.0 matrixStats_1.4.1
[23] S4Vectors_0.44.0 BiocGenerics_0.52.0
[25] generics_0.1.3
loaded via a namespace (and not attached):
[1] splines_4.4.1 later_1.4.1 filelock_1.0.3
[4] bitops_1.0-9 svgPanZoom_0.3.4 tibble_3.2.1
[7] polyclip_1.10-7 XML_3.99-0.17 lifecycle_1.0.4
[10] rstatix_0.7.2 lattice_0.22-6 MASS_7.3-61
[13] backports_1.5.0 magrittr_2.0.3 limma_3.62.1
[16] plotly_4.10.4 rmarkdown_2.29 yaml_2.3.10
[19] httpuv_1.6.15 doRNG_1.8.6.1 sp_2.1-4
[22] dcanr_1.22.0 spatstat.sparse_3.1-0 DBI_1.2.3
[25] minqa_1.2.8 RColorBrewer_1.1-3 abind_1.4-8
[28] zlibbioc_1.52.0 purrr_1.0.2 RCurl_1.98-1.16
[31] rappdirs_0.3.3 tweenr_2.0.3 GenomeInfoDbData_1.2.13
[34] spatstat.utils_3.1-1 terra_1.7-78 genefilter_1.88.0
[37] pheatmap_1.0.12 goftest_1.2-3 simpleSeg_1.8.0
[40] annotate_1.84.0 spatstat.random_3.3-2 svglite_2.1.3
[43] codetools_0.2-20 DelayedArray_0.32.0 ggforce_0.4.2
[46] tidyselect_1.2.1 raster_3.6-30 UCSC.utils_1.2.0
[49] farver_2.1.2 viridis_0.6.5 lme4_1.1-35.5
[52] spatstat.explore_3.3-3 jsonlite_1.8.9 Formula_1.2-5
[55] iterators_1.0.14 systemfonts_1.1.0 foreach_1.5.2
[58] tools_4.4.1 Rcpp_1.0.13-1 glue_1.8.0
[61] gridExtra_2.3 SparseArray_1.6.0 xfun_0.49
[64] ranger_0.17.0 mgcv_1.9-1 ggthemes_5.1.0
[67] EBImage_4.48.0 dplyr_1.1.4 HDF5Array_1.34.0
[70] scam_1.2-17 shinydashboard_0.7.2 withr_3.0.2
[73] numDeriv_2016.8-1.1 BiocManager_1.30.25 fastmap_1.2.0
[76] ggh4x_0.2.8 boot_1.3-31 rhdf5filters_1.18.0
[79] fansi_1.0.6 digest_0.6.37 R6_2.5.1
[82] mime_0.12 colorspace_2.1-1 tensor_1.5
[85] RSQLite_2.3.8 jpeg_0.1-10 spatstat.data_3.1-4
[88] utf8_1.2.4 tidyr_1.3.1 data.table_1.16.2
[91] class_7.3-22 httr_1.4.7 htmlwidgets_1.6.4
[94] S4Arrays_1.6.0 pkgconfig_2.0.3 gtable_0.3.6
[97] blob_1.2.4 XVector_0.46.0 htmltools_0.5.8.1
[100] carData_3.0-5 fftwtools_0.9-11 scales_1.3.0
[103] ggupset_0.4.0 png_0.1-8 spatstat.univar_3.1-1
[106] knitr_1.49 rstudioapi_0.17.1 reshape2_1.4.4
[109] rjson_0.2.23 curl_6.0.1 nlme_3.1-166
[112] nloptr_2.1.1 bdsmatrix_1.3-7 cachem_1.1.0
[115] rhdf5_2.50.0 stringr_1.5.1 BiocVersion_3.20.0
[118] vipor_0.4.7 parallel_4.4.1 concaveman_1.1.0
[121] AnnotationDbi_1.68.0 pillar_1.9.0 grid_4.4.1
[124] vctrs_0.6.5 coxme_2.2-22 promises_1.3.2
[127] ggpubr_0.6.0 car_3.1-3 xtable_1.8-4
[130] beeswarm_0.4.0 evaluate_1.0.1 magick_2.8.5
[133] cli_3.6.3 locfit_1.5-9.10 compiler_4.4.1
[136] rlang_1.1.4 crayon_1.5.3 rngtools_1.5.2
[139] ggsignif_0.6.4 labeling_0.4.3 plyr_1.8.9
[142] ggbeeswarm_0.7.2 stringi_1.8.4 nnls_1.6
[145] viridisLite_0.4.2 deldir_2.0-4 cytomapper_1.18.0
[148] Biostrings_2.74.0 lmerTest_3.1-3 munsell_0.5.1
[151] lazyeval_0.2.2 tiff_0.1-12 spatstat.geom_3.3-4
[154] Matrix_1.7-1 bit64_4.5.2 Rhdf5lib_1.28.0
[157] KEGGREST_1.46.0 statmod_1.5.0 shiny_1.9.1
[160] memoise_2.0.1 igraph_2.1.1 broom_1.0.7
[163] bit_4.5.0