--- title: "Statial" date: "`r BiocStyle::doc_date()`" params: test: FALSE author: - name: Farhan Ameen affiliation: - School of Mathematics and Statistics, University of Sydney, Australia - name: Sourish Iyengar affiliation: - School of Mathematics and Statistics, University of Sydney, Australia - name: Alex Qin affiliation: - School of Mathematics and Statistics, University of Sydney, Australia - name: Ellis Patrick affiliation: - &WIMR Westmead Institute for Medical Research, University of Sydney, Australia - School of Mathematics and Statistics, University of Sydney, Australia vignette: > %\VignetteIndexEntry{"Introduction to Statial"} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} output: BiocStyle::html_document --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", cache = TRUE, message = FALSE, warning = FALSE ) library(BiocStyle) ``` # Installation ```{r eval = FALSE} # Install the package from Bioconductor if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("Statial") ``` # Load packages ```{r warning = FALSE, message = FALSE} # Loading required packages library(Statial) library(spicyR) library(ClassifyR) library(lisaClust) library(dplyr) library(SingleCellExperiment) library(ggplot2) library(ggsurvfit) library(survival) library(tibble) theme_set(theme_classic()) nCores <- 1 ``` # Overview There are over 37 trillion cells in the human body, each taking up different forms and functions. The behaviour of these cells can be described by canonical characteristics, but their functions can also dynamically change based on their environmental context, leading to cells with diverse states. Understanding changes in cell state and the interplay between cells is key to understanding their mechanisms of action and how they contribute to human disease. `Statial` is a suite of functions for identifying changes in cell state. This guide will provide a step-by-step overview of some key functions within `Statial`. # Loading example data To illustrate the functionality of Statial we will use a multiplexed ion beam imaging by time-of-flight (MIBI-TOF) dataset profiling tissue from triple-negative breast cancer patients$^1$ by Keren et al., 2018. This dataset simultaneously quantifies *in situ* expression of 36 proteins in 34 immune rich patients. *Note:* The data contains some "uninformative" probes and the original cohort included 41 patients. These images are stored in a `SingleCellExperiment` object called `kerenSCE`. This object contains 57811 cells across 10 images and includes information on cell type and patient survival. *Note:* The original dataset was reduced down from 41 images to 10 images for the purposes of this vignette, due to size restrictions. ```{r} # Load head and neck data data("kerenSCE") kerenSCE ``` # Kontextual: Identifying discrete changes in cell state `Kontextual` is a method to evaluate the localisation relationship between two cell types in an image. `Kontextual` builds on the L-function by contextualising the relationship between two cell types in reference to the typical spatial behaviour of a $3^{rd}$ cell type/population. By taking this approach, `Kontextual` is invariant to changes in the window of the image as well as tissue structures which may be present. The definitions of cell types and cell states are somewhat ambiguous, cell types imply well defined groups of cells that serve different roles from one another, on the other hand cell states imply that cells are a dynamic entity which cannot be discretised, and thus exist in a continuum. For the purposes of using `Kontextual` we treat cell states as identified clusters of cells, where larger clusters represent a "parent" cell population, and finer sub-clusters representing a "child" cell population. For example a CD4 T cell may be considered a child to a larger parent population of Immune cells. `Kontextual` thus aims to see how a child population of cells deviate from the spatial behaviour of their parent population, and how that influences the localisation between the child cell state and another cell state. ## Cell type hierarchy A key input for Kontextual is an annotation of cell type hierarchies. We will need these to organise all the cells present into cell state populations or clusters, e.g. all the different B cell types are put in a vector called bcells. For the purposes of this vignette, these will be manually defined. Alternatively, you can use the `treeKor` bioconductor package [treekoR](http://www.bioconductor.org/packages/release/bioc/html/treekoR.html) to define these hierarchies in a data driven way. ```{r} # Examine all cell types in image unique(kerenSCE$cellType) # Set up cell populations tumour <- c("Keratin_Tumour", "Tumour") bcells <- c("B") tcells <- c("CD3_Cell", "CD4_Cell", "CD8_Cell", "Tregs") myeloid <- c("Dc/Mono", "DC", "Mono/Neu", "Macrophages", "Neutrophils") endothelial <- c("Endothelial") mesenchymal <- c("Mesenchymal") tissue <- c(endothelial, mesenchymal) immune <- c(bcells, tcells, myeloid, "NK", "other immune") # NK = Natural Killer cells all <- c(tumour, tissue, immune, "Unidentified") ``` ## Discrete cell state changes within a single image Here we examine an image highlighted in the Keren et al. 2018 manuscript where the relationship between two cell types depends on a parent cell population. In image 6 of the Keren et al. dataset, we can see that *p53+ tumour cells* and *immune cells* are dispersed. However when the behaviour of *p53+ tumour cells* are placed in the context of the spatial behaviour of its broader parent population *tumour cells*, *p53+ tumour cells* and *immune* would appear localised.w ```{r} # Lets define a new cell type vector kerenSCE$cellTypeNew <- kerenSCE$cellType # Select for all cells that express higher than baseline level of p53 p53Pos = assay(kerenSCE)["p53",] > -0.300460 # Find p53+ tumour cells kerenSCE$cellTypeNew[kerenSCE$cellType %in% tumour] <- "Tumour" kerenSCE$cellTypeNew[p53Pos & kerenSCE$cellType %in% tumour] <- "p53_Tumour" #Group all immune cells under the name "Immune" kerenSCE$cellTypeNew[kerenSCE$cellType %in% immune] <- "Immune" # Plot image 6 kerenSCE |> colData() |> as.data.frame() |> filter(imageID == "6") |> filter(cellTypeNew %in% c("Immune", "Tumour", "p53_Tumour")) |> arrange(cellTypeNew) |> ggplot(aes(x = x, y = y, color = cellTypeNew)) + geom_point(size = 1) + scale_colour_manual(values = c("#505050", "#64BC46","#D6D6D6")) + guides(colour = guide_legend(title = "Cell types", override.aes = list(size=3))) ``` `Kontextual` accepts a `SingleCellExperiment` object, a single image, or list of images from a `SingleCellExperiment` object, which gets passed into the `cells` argument. The two cell types which will be evaluated are specified in the `to` and `from` arguments. A parent population must also be specified in the `parent` argument, note the parent cell population must include the `to` cell type. The argument `r` will specify the radius which the cell relationship will be evaluated on. `Kontextual` supports parallel processing, the number of cores can be specified using the `cores` argument. `Kontextual` can take a single value or multiple values for each argument and will test all combinations of the arguments specified. We can calculate these relationships across all images for a single radius. ```{r} p53_Kontextual <- Kontextual( cells = kerenSCE, r = 100, from = "Immune", to = "p53_Tumour", parent = c("p53_Tumour", "Tumour"), cellType = "cellTypeNew" ) p53_Kontextual ``` The `kontextCurve` function plots the L-function value and Kontextual values over a range of radii. If the points lie above the red line (expected pattern) then localisation is indicated for that radius, if the points lie below the red line then dispersion is indicated. As seen in the following plot Kontextual is able to correctly identify localisation between `Immune` and `p53` in the example image for a certain range of radii. When the radius gets too large the overall relationship between `Immune` and `p53` looks dispersed. The original L-function is not able to identify localisation at any value of radii. ```{r} curves <- kontextCurve( cells = kerenSCE, from = "Immune", to = "p53_Tumour", parent = c("p53_Tumour", "Tumour"), rs = seq(50, 510, 50), image = "6", cellType = "cellTypeNew", cores = nCores ) kontextPlot(curves) ``` Alternatively all pairwise cell relationships and their corresponding parent in the dataset can be tested. A data frame with all pairwise combinations can be creating using the `parentCombinations` function. This function takes in a vector of all the cells, as well as all the parent vectors set up earlier. As shown below the output is a data frame specifying the `to`, `from`, and `parent` arguments for `Kontextual`. ```{r} # Get all relationships between cell types and their parents parentDf <- parentCombinations( all = all, tumour, bcells, tcells, myeloid, endothelial, mesenchymal, tissue, immune ) ``` ## Discrete cell state changes across all images Rather than specifying `to`, `from`, and `parent` in Kontextual, the output from `parentCombinations` can be inputed into `Kontextual` using the `parentDf` argument, to examine all pairwise relationships in the dataset. This chunk will take a signficant amount of time to run, for demonstration the results have been saved and are loaded in. ```{r eval=FALSE} # Running Kontextual on all relationships across all images. kerenKontextual <- Kontextual( cells = kerenSCE, parentDf = parentDf, r = 100, cores = nCores ) ``` ```{r} data("kerenKontextual") bigDiff <- (kerenKontextual$original - kerenKontextual$kontextual) head(kerenKontextual[order(bigDiff),], 10) ``` ## Associate discrete state changes with survival outcomes To examine whether the features obtained from `Statial` are associated with patient outcomes or groupings, we can use the `colTest` function from `SpicyR`. To understand if survival outcomes differ significantly between 2 patient groups, specify `type = "survival"` in `colTest`. Here we examine which features are most associated with patient survival using the Kontextual values as an example. To do so, survival data is extracted from `kerenSCE` and converted into the survival object `kerenSurv`. ```{r} # Extracting survival data survData = kerenSCE |> colData() |> data.frame() |> select(imageID, Survival_days_capped, Censored) |> mutate(event = 1 - Censored) |> unique() # Creating survival vector kerenSurv = Surv(survData$Survival_days_capped, survData$event) names(kerenSurv) = survData$imageID ``` In addition to this, the Kontextual results must be converted from a `data.frame` to a wide `matrix`, this can be done using `prepMatrix`. Note, to extract the original L-function values, specify `column = "original"` in `prepMatrix`. ```{r} # Converting Kontextual result into data matrix kontextMat = prepMatrix(kerenKontextual) # Ensuring rownames of kontextMat match up with rownames of the survival vector kontextMat = kontextMat[names(kerenSurv), ] # Replace NAs with 0 kontextMat[is.na(kontextMat )] <- 0 ``` Finally, both the Kontextual matrix and survival object are passed into `colTest`, with `type = "survival"` to obtain the survival results. ```{r} # Running survival analysis survivalResults = spicyR::colTest(kontextMat, kerenSurv, type = "survival") head(survivalResults) ``` As we can see from the results `CD8_Cell__Neutrophils__myeloid` is the one of the most significant pairwise relationship which contributes to patient survival. That is the relationship between CD8 T cells and Neutrophils, relative to the parent population of all myeloid cells. We can see that there is a negative coefficient associated with this relationship, which tells us an increase in localisation of CD8 T cells and Neutrophils leads to poorer survival outcomes for patients. The association between `CD8_Cell__Neutrophils__myeloid` and survival can also be visualised on a Kaplan-Meier curve. We must first extract the Kontextual values of this relationship across all images. Next we determine if CD8 T cells and Neutrophils are relatively attracted or avoiding in each image, by comparing the Kontextual value in each image to the median Kontextual value. Finally we plot the Kaplan-Meier curve using the `ggsurvfit` package. As shown below, when Neutrophils and Dc/Mono are relatively more localised to one another, patients tend to have worse survival outcomes. ```{r, fig.width=5, fig.height=4} # Selecting most significant relationship survRelationship = kontextMat[["CD8_Cell__Neutrophils__myeloid"]] survRelationship = ifelse(survRelationship > median(survRelationship), "Localised", "Dispersed") # Plotting Kaplan-Meier curve survfit2(kerenSurv ~ survRelationship) |> ggsurvfit() + ggtitle("CD8_Cell__Neutrophils__myeloid") ``` # SpatioMark: Identifying continuous changes in cell state Changes in cell states can be analytically framed as the change in abundance of a gene or protein within a particular cell type. We can use marker expression to identify and quantify evidence of cell interactions that catalyse cell state changes. This approach measures how protein markers in a cell change with spatial proximity and abundance to other cell types. The methods utilised here will thereby provide a framework to explore how the dynamic behaviour of cells are altered by the agents they are surrounded by. ## Continuous cell state changes within a single image The first step in analysing these changes is to calculate the spatial proximity (`getDistances`) and abundance (`getAbundances`) of each cell to every cell type. These values will then be stored in the `reducedDims` slot of the `SingleCellExperiment` object under the names `distances` and `abundances` respectively. ```{r} kerenSCE <- getDistances(kerenSCE, maxDist = 200, nCores = 1) kerenSCE <- getAbundances(kerenSCE, r = 200, nCores = 1) ``` First, let's examine the same effect observed earlier with Kontextual - the localisation between p53-positive keratin/tumour cells and macrophages in the context of total keratin/tumour cells for image 6 of the Keren et al. dataset. Statial provides two main functions to assess this relationship - `calcStateChanges` and `plotStateChanges`. We can use `calcStateChanges` to examine the relationship between 2 cell types for 1 marker in a specific image. In this case, we're examining the relationship between keratin/tumour cells (`from = Keratin_Tumour`) and macrophages (`to = "Macrophages"`) for the marker p53 (`marker = "p53"`) in `image = "6"`. We can appreciate that the `fdr` statistic for this relationship is significant, with a negative tvalue, indicating that the expression of p53 in keratin/tumour cells decreases as distance from macrophages increases. ```{r} stateChanges <- calcStateChanges( cells = kerenSCE, type = "distances", image = "6", from = "Keratin_Tumour", to = "Macrophages", marker = "p53", nCores = 1) stateChanges ``` Statial also provides a convenient function for visualising this interaction - `plotStateChanges`. Here, again we can specify `image = 6` and our main cell types of interest, keratin/tumour cells and macrophages, and our marker p53, in the same format as `calcStateChanges`. Through this analysis, we can observe that keratin/tumour cells closer to a group of macrophages tend to have higher expression of p53, as observed in the first graph. This relationship is quantified with the second graph, showing an overall decrease of p53 expression in keratin/tumour cells as distance to macrophages increase. These results allow us to essentially arrive at the same result as Kontextual, which calculated a localisation between p53+ keratin/tumour cells and macrophages in the wider context of keratin/tumour cells. ```{r} p <- plotStateChanges( cells = kerenSCE, type = "distances", image = "6", from = "Keratin_Tumour", to = "Macrophages", marker = "p53", size = 1, shape = 19, interactive = FALSE, plotModelFit = FALSE, method = "lm") p ``` ## Continuous cell state changes across all images Beyond looking at single cell-to-cell interactions for a single image, we can also look at all interactions across all images. The `calcStateChanges` function provided by Statial can be expanded for this exact purpose - by not specifying cell types, a marker, or an image, `calcStateChanges` will examine the most significant correlations between distance and marker expression across the entire dataset. Here, we've filtered out the most significant interactions to only include those found within image 6 of the Keren et al. dataset. ```{r} stateChanges <- calcStateChanges( cells = kerenSCE, type = "distances", nCores = 1, minCells = 100) stateChanges |> filter(imageID == 6) |> head(n = 10) ``` In image 6, the majority of the top 10 most significant interactions occur between keratin/tumour cells and an immune population, and many of these interactions appear to involve the HLA class I ligand. We can examine some of these interactions further with the `plotStateChanges` function. Taking a closer examination of the relationship between macrophages and keratin/tumour HLA class I expression, the plot below shows us a clear visual correlation - as macrophage density increases, keratin/tumour cells increase their expression HLA class I. Biologically, HLA Class I is a ligand which exists on all nucleated cells, tasked with presenting internal cell antigens for recognition by the immune system, marking aberrant cells for destruction by either CD8+ T cells or NK cells. ```{r} p <- plotStateChanges( cells = kerenSCE, type = "distances", image = "6", from = "Keratin_Tumour", to = "Macrophages", marker = "HLA_Class_1", size = 1, shape = 19, interactive = FALSE, plotModelFit = FALSE, method = "lm") p ``` Next, let's take a look at the top 10 most significant results across all images. ```{r} stateChanges |> head(n = 10) ``` Immediately, we can appreciate that a couple of these interactions are not biologically plausible. One of the most significant interactions occurs between B cells and CD4 T cells in image 35, where CD4 T cells are found to increase in CD20 expression when in close proximity to B cells. Biologically, CD20 is a highly specific ligand for B cells, and under healthy circumstances are usually not expressed in T cells. Could this potentially be an artefact of `calcStateChanges`? We can examine the image through the `plotStateChanges` function, where we indeed observe a strong increase in CD20 expression in T cells nearby B cell populations. ```{r} p <- plotStateChanges( cells = kerenSCE, type = "distances", image = "35", from = "CD4_Cell", to = "B", marker = "CD20", size = 1, shape = 19, interactive = FALSE, plotModelFit = FALSE, method = "lm") p ``` So why are T cells expressing CD20? This brings us to a key problem of cell segmentation - contamination. ## Contamination (Lateral marker spill over) Contamination, or lateral marker spill over is an issue that results in a cell’s marker expressions being wrongly attributed to another adjacent cell. This issue arises from incorrect segmentation where components of one cell are wrongly determined as belonging to another cell. Alternatively, this issue can arise when antibodies used to tag and measure marker expressions don't latch on properly to a cell of interest, thereby resulting in residual markers being wrongly assigned as belonging to a cell near the intended target cell. It is important that we either correct or account for this incorrect attribution of markers in our modelling process. This is critical in understanding whether significant cell-cell interactions detected are an artefact of technical measurement errors driven by spill over or are real biological changes that represent a shift in a cell’s state. To circumvent this problem, Statial provides a function that predicts the probability that a cell is any particular cell type - `calcContamination`. `calcContamination` returns a dataframe of probabilities demarcating the chance of a cell being any particular cell type. This dataframe is stored under `contaminations` in the `reducedDim` slot of the `SingleCellExperiment` object. It also provides the `rfMainCellProb` column, which provides the probability that a cell is indeed the cell type it has been designated. E.g. For a cell designated as CD8, rfMainCellProb could give a 80% chance that the cell is indeed CD8, due to contamination. We can then introduce these probabilities as covariates into our linear model by setting `contamination = TRUE` as a parameter in our `calcStateChanges` function. However, this is not a perfect solution for the issue of contamination. As we can see, despite factoring in contamination into our linear model, the correlation between B cell density and CD20 expression in CD4 T cells remains one of the most significant interactions in our model. ```{r} kerenSCE <- calcContamination(kerenSCE) stateChangesCorrected <- calcStateChanges( cells = kerenSCE, type = "distances", nCores = 1, minCells = 100, contamination = TRUE) stateChangesCorrected |> head(n = 20) ``` However, this does not mean factoring in contamination into our linear model was ineffective. Whilst our correction attempts do not rectify every relationship which arises due to contamination, we show that a significant portion of these relationships are rectified. We can show this by plotting a ROC curve of true positives against false positives. In general, cell type specific markers such as CD4, CD8, and CD20 should not change in cells they are not specific to. Therefore, relationships detected to be significant involving these cell type markers are likely false positives and will be treated as such for the purposes of evaluation. Meanwhile, cell state markers are predominantly likely to be true positives. Plotting the relationship between false positives and true positives, we'd expect the contamination correction to be greatest in the relationships with the top 100 lowest p values, where we indeed see more true positives than false positives with contamination correction. ```{r} cellTypeMarkers <- c("CD3", "CD4", "CD8", "CD56", "CD11c", "CD68", "CD45", "CD20") values = c("blue", "red") names(values) <- c("None", "Corrected") df <- rbind(data.frame(TP =cumsum(stateChanges$marker %in% cellTypeMarkers), FP = cumsum(!stateChanges$marker %in% cellTypeMarkers), type = "None"), data.frame(TP =cumsum(stateChangesCorrected$marker %in% cellTypeMarkers), FP = cumsum(!stateChangesCorrected$marker %in% cellTypeMarkers), type = "Corrected")) ggplot(df, aes(x = TP, y = FP, colour = type)) + geom_line()+ labs(y = "Cell state marker", x = "Cell type marker") + scale_colour_manual(values = values) ``` Here, we zoom in on the ROC curve where the top 100 lowest p values occur, where we indeed see more true positives than false positives with contamination correction. ```{r} ggplot(df, aes(x = TP, y = FP, colour = type)) + geom_line()+ xlim(0,100) + ylim(0,1000)+ labs(y = "Cell state marker", x = "Cell type marker") + scale_colour_manual(values = values) ``` ## Associate continuous state changes with survival outcomes Similiar to `Kontextual`, we can run a similar survival analysis using our state changes results. Here, `prepMatrix` extracts the coefficients, or the `coef` column of `stateChanges` by default. To use the t values instead, specify `column = "tval"` in the `prepMatrix` function. ```{r} # Preparing features for Statial stateMat <- prepMatrix(stateChanges) # Ensuring rownames of stateMat match up with rownames of the survival vector stateMat <- stateMat[names(kerenSurv), ] # Remove some very small values stateMat <- stateMat[,colMeans(abs(stateMat)>0.0001)>.8] survivalResults <- colTest(stateMat, kerenSurv, type = "survival") head(survivalResults) ``` For our state changes results, `Keratin_Tumour__CD4_Cell__Keratin6` is the most significant pairwise relationship which contributes to patient survival. That is, the relationship between HLA class I expression in keratin/tumour cells and their spatial proximity to mesenchymal cells. As there is a negative coeffcient associated with this relationship, which tells us that higher HLA class I expression in keratin/tumour cells nearby mesenchymal cell populations lead to poorer survival outcomes for patients. ```{r, fig.width=5, fig.height=4} # Selecting the most significant relationship survRelationship = stateMat[["Keratin_Tumour__CD4_Cell__Keratin6"]] survRelationship = ifelse(survRelationship > median(survRelationship), "Higher expression in close cells", "Lower expression in close cells") # Plotting Kaplan-Meier curve survfit2(kerenSurv ~ survRelationship) |> ggsurvfit() + add_pvalue() + ggtitle("Keratin_Tumour__CD4_Cell__Keratin6") ``` # Region analysis using lisaClust Next we can cluster areas with similar spatial interactions to identify regions using lisaClust. Here we set `k = 5` to identify 5 regions. ```{r} set.seed(51773) # Preparing features for lisaClust kerenSCE <- lisaClust::lisaClust(kerenSCE, k = 5) ``` The regions identified by licaClust can be visualised using the `hatchingPlot` function. ```{r, fig.height=5, fig.width=6.5} # Use hatching to visualise regions and cell types. lisaClust::hatchingPlot(kerenSCE, useImages = "5", line.spacing = 41, # spacing of lines nbp = 100 # smoothness of lines ) ``` # Marker Means `Statial` provides functionality to identify the average marker expression of a given cell type in a given region, using the `getMarkerMeans` function. Similar to the analysis above, these features can also be used for survival analysis. ```{r lisaClust} cellTypeRegionMeans <- getMarkerMeans(kerenSCE, imageID = "imageID", cellType = "cellType", region = "region") survivalResults = colTest(cellTypeRegionMeans[names(kerenSurv),], kerenSurv, type = "survival") head(survivalResults) ``` # Patient classification Finally we demonstrate how we can use `ClassifyR` to perform patient classification with the features generated from `Statial`. In addition to the kontextual, state changes, and marker means values, we also calculate cell type proportions and region proportions using the `getProp` function in `spicyR`. Here we perform 3 fold cross validation with 10 repeats, using a CoxPH model for survival classification, feature selection is also performed by selecting the top 10 features per fold using a CoxPH model. ```{r} # Calculate cell type and region proportions cellTypeProp <- getProp(kerenSCE, feature = "cellType", imageID = "imageID") regionProp <- getProp(kerenSCE, feature = "region", imageID = "imageID") # Combine all the features into a list for classification featureList <- list(states = stateMat, kontextual = kontextMat, regionMarkerMeans = cellTypeRegionMeans, cellTypeProp = cellTypeProp, regionProp = regionProp ) # Ensure the rownames of the features match the order of the survival vector featureList <- lapply(featureList, function(x)x[names(kerenSurv),]) set.seed(51773) kerenCV = crossValidate( measurements = featureList, outcome = kerenSurv, classifier = "CoxPH", selectionMethod = "CoxPH", nFolds = 5, nFeatures = 10, nRepeats = 20, nCores = 1 ) ``` Here, we use the `performancePlot` function to assess the C-index from each repeat of the 3-fold cross-validation. We can see the resulting C-indexes are very variable due to the dataset only containing 10 images. ```{r} # Calculate AUC for each cross-validation repeat and plot. performancePlot(kerenCV, characteristicsList = list(x = "Assay Name") ) + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) ``` # References Keren, L., Bosse, M., Marquez, D., Angoshtari, R., Jain, S., Varma, S., Yang, S. R., Kurian, A., Van Valen, D., West, R., Bendall, S. C., & Angelo, M. (2018). A Structured Tumor-Immune Microenvironment in Triple Negative Breast Cancer Revealed by Multiplexed Ion Beam Imaging. Cell, 174(6), 1373-1387.e1319. ([DOI](https://doi.org/10.1016/j.cell.2018.08.039)) # sessionInfo ```{r} sessionInfo() ```