Unlocking single cell spatial omics analyses with scdney
Published
June 13, 2025
Presenting authors
Harry Robertson\(^{1,2,3}\), Lijia Yu\(^{1,2}\), Beilei Bian\(^{1,2}\), Andrew Zhang\(^{1,4}\), Jean Yang\(^{1,2,3}\).
Contributors
Yue Cao\(^{1,2}\), Lijia Yu\(^{1,2}\), Andy Tran\(^{1,2}\), Daniel Kim\(^{1,3}\), Dario Strbenac\(^{1,2}\), Nicholas Robertson\(^{1}\), Helen Fu\(^{3,4}\), Jean Yang\(^{1,2,4}\).
\(^1\) Sydney Precision Data Science Centre, University of Sydney, Australia \(^2\) School of Mathematics and Statistics, University of Sydney, Australia \(^3\) Charles Perkins Centre, University of Sydney, Australia \(^4\) School of Computer Science, University of Sydney, Australia
Contact: jean.yang@sydney.edu.au
Overview
The emergence of high-resolution spatial omics technologies has revolutionized our ability to map cellular ecosystems in situ. This workshop explores approaches in multi-sample spatial data analysis. We cover end-to-end workflow considerations—from experimental design and QC to spatial feature interpretation—with case studies in disease prediction.
Understand approaches for quality control of spatial proteomics data.
Calculate spatial statistics at the cell-type level using scFeatures.
Perform multi-view disease outcome prediction with the package ClassifyR.
Develop an understanding of:
evaluation of classification and survival models .
evaluate cohort heterogeneity given a survival model.
Explore various strategies for disease outcome prognosis using spatial omics data.
Our question of interest
We want to know if the risk of recurrence in the METABRIC breast cancer cohort can be accurately estimated to inform how aggressively they need to be treated.
Part 0: Loading an example image in R
To begin, we will start with loading an example image in R. It is not neccessary for you to be familiar with this part, as we will provide you with pre-processed data for the rest of the workshop.
suppressPackageStartupMessages({library(cytomapper)library(EBImage)library(S4Vectors)})img_path <-"../data/raw_image_example/MB0002_1_345_fullstack.tiff"mask_path <-"../data/raw_image_example/MB0002_1_345_cellmask.tiff"sample_id <-"../data/raw_image_example/MB0002_1_345"#Load image (multi-channel)images <-loadImages(img_path, single_channel =FALSE)names(images) <- sample_id#Load maskmask_img <-readImage(mask_path,as.is =T)masks <-CytoImageList(list(mask_img))names(masks) <- sample_id# Add metadata rows so img_id can be matched by COLUMN NAMEmcols(images) <-DataFrame(sample_id =names(images))mcols(masks) <-DataFrame(sample_id =names(masks))# Build channel names from CSVs and assign themmetal_order <-read.csv("../data/raw_image_example/channel_to_metal_order.csv",header =FALSE, stringsAsFactors =FALSE)names(metal_order) <-"metal"metal_order$Channel <-seq_len(nrow(metal_order))metal_map <-read.csv("../data/raw_image_example/metal_to_marker.csv",header =TRUE, stringsAsFactors =FALSE)merged <-merge(metal_order, metal_map, by ="metal", all.x =TRUE)merged <- merged[order(merged$Channel), ]labels <-ifelse(is.na(merged$marker), merged$metal, merged$marker)# Assign channel names to the multi-channel imagechannelNames(images) <- labelsmarker =plotPixels(image = images,colour_by ="Total HH3")
Code
# Plot the maskplotCells(masks)
Code
# Plot overlay combine =plotPixels(image = images,mask = masks,img_id ="sample_id", # column name in mcols(images)/mcols(masks)colour_by ="Total HH3")
Part 1: Exploring the data
1.1: Initial data exploration
At the start of any analysis pipeline, it is often good to explore the data to get a sense of the structure and its complexity. Let’s explore the data to answer the questions below:
Questions
How many features and observations are there in the data and what do the features represent?
What covariates are in our data?
Given our question of interest, what variable would be our outcome variable?
# How many features and observations are in our dataset?dim(data_sce)
[1] 38 76307
In addition to the proteomics data, it’s important to understand what other covariates, such as clinical variables, are in our dataset. These can help us in answering our question(s) of interest or formulate new questions. For example, we can’t explore the association between smoking status and breast cancer severity if the variable for smoking status doesn’t exist in our data.
Now that we have a basic idea of what our data looks like, we can look at it in more detail. While initial data exploration reveals fundamental patterns, deeper examination is very helpful. As it serves two critical purposes: first, to detect anomalies or biases requiring remediation, and second, to inform our choice of analytical methods tailored to the biological questions at hand.
Examining the proportion of missing values in a dataset is crucial for ensuring the accuracy and validity of data analysis. Missing values can significantly impact the reliability of statistical results, potentially leading to biased conclusions and reduced statistical power. Here, we plot the proportion of missing values for each of our clinical variables.
Code
# What is the proportion of missing values in our clinical data?gg_miss_var(clinical, show_pct =TRUE) +theme_bw()
Looking at imbalance in data is crucial because it can lead to biased models and inaccurate predictions, especially in classification tasks. Imbalance occurs when one class is significantly underrepresented compared to others. This can cause models to be overly influenced by the majority class, leading to poor performance on the minority class. Here, we tabularise some of the relevant clinical variables and plot the distribution of the Time to Recurrence-Free Survival variable.
Code
# Estrogen receptor statusggplot(clinical, aes(x = ER.Status)) +geom_bar(fill ="#0099B4", color ="white", alpha =0.8) +labs(y ="Count\n", x ="\nER status") +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold"))
Code
# Histological typeggplot(clinical, aes(x = Histological.Type)) +geom_bar(fill ="#fc6203", color ="white", alpha =0.8) +labs(y ="Count\n", x ="\nHistological Type") +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold"))
Code
# eventRFS: "Event in Recurrence-Free Survival."It indicates whether the event has occurred.#ggplot(clinical, aes(x =factor(eventRFS))) +geom_bar(fill ="#40dbb2", color ="white", alpha =0.8) +labs(y ="Count\n", x ="\nEvent in Recurrence-Free Survival") +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold"))
Code
# timeRFS: "Time to Recurrence-Free Survival." It is the time period until recurrence occurs. ggplot(clinical, aes(x = timeRFS)) +geom_histogram(fill ="#de9921", color ="white", alpha =0.8, bins=20) +labs(y ="Frequency\n", x ="\ntimeRFS") +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold"))
Here, we explore the distribution of the outcomes and variables in the meta-data. We use cross-tabulation to examine the following variables: Surgery vs death, ER status, and Grade.
Code
# Stage and deathtable(clinical$Breast.Surgery, clinical$Death, useNA ="ifany")
# ER status and gradetable(clinical$ER.Status, clinical$Grade)
1 2 3
neg 0 0 11
pos 13 30 17
Code
# "Number of individuals based on Gradetable(clinical$Grade, clinical$Death)
0 1
1 12 2
2 25 5
3 17 11
To assess potential batch effects and sample-specific clustering patterns, we visualize the cells in our data colored by the sample of origin. This qualitative inspection provides an intuitive first assessment of any batch effects in our data integration. If there is batch effects then we’d need to apply batch correction to resolve this issue.
Code
#data_sce <- runUMAP(data_sce, scale=TRUE)# With the UMAP function we can highlight by meta data of interest# Here we highlight the UMAP by sample IDplotUMAP(data_sce, colour_by ="metabricId") +theme(legend.position ="none") +coord_equal()
Discussion
Are there any anomalies or biases that might need to be corrected to ensure our analysis is robust?
Given our question of interest and the characteristics of our data, are there any particular analytic techniques that would be appropriate? Are there any that would not be appropriate?
Part 2: Quality control
Here, we explore key approaches for evaluating the quality of IMC data. A robust assessment should consider multiple factors, including the density of marker expression, marker correlations, and their co-expression. Something we might want to think about is what characeristics might indicate whether a sample(s) is low/high quality?
Here, we want to assess the quality of the markers: it is desirable to see at least two peaks, which would indicate the existence of cell type specific markers in our data.
Questions
Can we see any cell type specific markers, which are they?
Are there any non-cell type specific markers?
What can you notice about the distribution of markers across images?
We assess whether the marker intensities require transformation or normalisation. This step is important for two main reasons:
Skewed distributions: Marker intensities are often right-skewed, which can distort downstream analyses such as clustering or dimensionality reduction.
Inconsistent scales across images: Due to technical variation, the same marker may show very different intensity ranges across images. This can shift what’s considered “positive” or “negative” expression, making it difficult to label cells accurately.
By applying transformation and normalisation, we aim to stabilise variance and bring the data onto a more comparable scale across images.
Do the CD31+ and CD31- peaks clearly separate out in the density plot? To ensure that downstream analysis goes smoothly, we want our cell type specific markers to show 2 distinct peaks representing our CD31+ and CD31- cells. If we see 3 or more peaks where we don’t expect, this might be an indicator that further normalisation is required.
Are our CD31+ and CD31- peaks consistent across our images? We want to make sure that our density plots for CD3 are largely the same across images so that a CD3+ cell in one image is equivalent to a CD3+ cell in another image.
Code
# leave out the nuclei markers from our normalisation processuseMarkers <-rownames(data_sce)[!rownames(data_sce) %in%c("DNA1", "DNA2", "HH3", "HH3_total", "HH3_ph")]# transform and normalise the marker expression of each cell typedata_sce <-normalizeCells(data_sce,markers = useMarkers,transformation =NULL,method =c("trim99", "minMax", "PC1"),assayIn ="counts",imageID ="metabricId")selected_patients =which(colData(data_sce)$metabricId %in% cell_counts[1:10,]$metabricId)norm =plot_marker_densities(data_sce[,selected_patients], sample_id_col ="metabricId", markers_to_plot ="vWF_CD31", assay_name ="norm")$vWF_CD31 +theme(legend.position ="none") +ggtitle("Normalized Counts - CD31")log =plot_marker_densities(data_sce[,selected_patients], sample_id_col ="metabricId", markers_to_plot =c("vWF_CD31"), assay_name ="logcounts")$vWF_CD31 +theme(legend.position ="none") +ggtitle("Log Counts - CD31")ggarrange(log, norm, ncol=2)
In the plot above, the normalised data appeas more bimodal. We can observe one clear CD31- peak at around 0.50, and a CD31+ peak at 1.00. Image-level batch effects also appear to have been mitigated, since most peaks occur at around the same CD31 intensity.
Discussion
Do we have a sufficient amount of cell type specific markers?
Is the panel of genes sufficient for our study?
2.2: QC of individual samples
Here we plot the correlation of different markers for a subset of samples. This helps us determine whether a sample has low/high contamination of markers.
Questions
Which pairs of markers do we expect to be correlated/not correlated.
Is there any evidence of marker contamination for a given sample?
Are there any samples you would like to remove from the data?
Which samples would you like to keep?
Code
cell_counts <- IMC@colData |>as.data.frame() |> dplyr::count(metabricId, name ="cell_count") |># Count rows per metabricId dplyr::arrange(desc(cell_count)) mean_coexp_per_sample <- coexp_df %>%select(-cell_id, -cell_type) %>%# remove columns not related to prob valuesgroup_by(metabricId) %>%summarise(mean_coexp_prob =mean(c_across(where(is.numeric)), na.rm =TRUE)) %>%ungroup()merged_df <-left_join(mean_coexp_per_sample, cell_counts, by ="metabricId")# Scatter plotggplot(merged_df, aes(x = cell_count, y = mean_coexp_prob)) +geom_point(size =3, color ="steelblue") +labs(x ="Cell Count per Sample",y ="Mean Co-expression",title ="Mean Co-expression vs. Cell Count" ) +theme_minimal(base_size =14)
2.4: Does my data agree with literature?
Another way to QC our data is to explore whether the DE genes or associations we find matches current literature. For example, we have a breast cancer cohort. Breast cancer is very well studied and so it should be very easy to perform a DE analysis after pseudobulking our data to see if it makes findings in previous bulk RNA-seq studies of breast cancer cohorts. We can also do this at the single-cell level and compare with previous single-cell proteomic/RNA-seq studies of breast cancer cohorts. We won’t be running the code below. But here is an example of how one might perform a basic DE analysis using package limma.
Code
# Note: The same approach can be used for pseudobulk or single-cell level data. # Just change `logcounts(data_sce)`# We want low ki67 to be the reference level (baseline level) data_sce$ki67 <-factor(data_sce$ki67, levels=c("low_ki67", "high_ki67"))# Here we specify the design matrix. We are specifying that Y are the expression values and X is ki67 status (but it could be any other variable - such as "good" or "bad" prognosis)# Y~X: Expression as a function of ki67 status. design_matrix <-model.matrix(~data_sce$ki67)# Fit modelfit <-lmFit(assay(data_sce, "norm"), design = design_matrix)# Estimate variance and SE of coefficientsfit <-eBayes(fit)tt <-topTable(fit, coef =2, n =5, adjust.method ="BH", sort.by ="p")tt
Here, we can also check whether the cell type annotations we have make sense. One way to do this is to check whether the marker expression levels agree with the cell type annotations. For example, we expect CD3 to be highly expressed in T cells, CD20 to be highly expressed in B cells, CD68 to be highly expressed in Macrophages, vWF_CD31 to be highly expressed in Endothelial cells, Vimentin to be highly expressed in Fibroblasts, and HER2 to be highly expressed in Tumor HER2+ cells.
Code
# Calculate average expression per cell typeavg_expr <-aggregate(t(assay(data_sce, "norm")), by =list(CellType = IMC$high_level_category), FUN = mean) %>% tibble::column_to_rownames(var ="CellType")avg_expr <- avg_expr[, c("CD20", "CD3", "CD68", "CD45", "vWF_CD31", "Vimentin", "HER2")] # select only our mark# Plot heatmap of marker by celltypepheatmap(avg_expr, scale ="column", main ="Average Marker Expression per Cell Type\n", cluster_rows =FALSE,color =colorRampPalette(c("blue", "white", "red"))(10000))
Discussion
How might be check whether these results make sense? What potential issues might we find?
Part 3: Exploring spatial data
This section assumes pre-existing cell type annotations. For unannotated data, two primary annotation strategies are available:
Unsupervised approach: cluster cells and identify marker genes for manual annotation.
Supervised approach: transfer labels from reference datasets using supervised learning/classification tools.
For supervised annotation, we recommend scClassify, a robust framework for cell-type classification. Below we visualise the data with and without spatial information before exploring the data further.
Cell types
Following initial assessment of potential batch effects through sample-origin visualization, we now examine cell-type-specific clustering patterns. Distinct, biologically meaningful clusters should emerge for each annotated cell type. The presence of heterogeneous clusters containing unrelated cell types may indicate incomplete or inaccurate cell-type annotations.
Code
plotUMAP(data_sce, colour_by ="description")
Spatial structure
The advantage with spatial omics is that we can examine the organisation of the cell-types as it occurs on the tissue slide. Here, we visualise one of the slides from a patient. We select a particular patient “MB-0002” and visualise the tissue sample from this patient using ggplot. Do we see any spatial patterning or does it look randomly distributed?
Code
# obtaining the meta data for this patient one_sample <- data_sce[, data_sce$metabricId =="MB-0002"]one_sample <-data.frame(colData(one_sample))ggplot(one_sample, aes(x = Location_Center_X, y = Location_Center_Y, colour = description)) +geom_point(alpha=0.7) +theme(panel.background=element_blank(),axis.line=element_line(color="black")) +ggtitle("Original slide")
Questions
What kinds of spatial information might be of interest given the question we’d like to answer?
How might we try to capture these spatial relationships?
Background: One of the most common questions or analyses in spatial data is spatial domain detection. “Spatial domains” are regions within a tissue where cells share similar gene expression profiles and are physically clustered together. Here, we will use the terminology “spatial domain” and “regions” interchangeably. Most common analytical strategies involve spatial clustering, where different methods use different levels of information, such as gene expression data only, cell type information, and spatial coordinates.
The strategy we demonstrate here uses the lisaClust function, which uses cell type information to cluster cells into five distinct spatial regions. As a case study, we will compare individuals with good or poor prognosis and examine them graphically to see if any regions appear to be different between good or poor prognosis. We define: - Good prognosis as individuals with > 10 years recurrence-free survival and - Poor prognosis as individuals with < 5 years recurrence-free survival.
Below we visualise the spatial domain (regions) detection result based on one individual. Here we will use the terminology “spatial domain” and “regions” interchangeably.
Depending on the number of regions, it may be more useful to visualise the spatial regions either collectively in a single graph or separately in multiple graphs. To visualise it in a single graph, the hatchingPlot() function is used to create hatching patterns for representating spatial regions and cell-types.
Code
## Extract time to recurrence-free survivalclinical$survivalgroup <-"neither"## Define poor and good prognosisclinical$survivalgroup[which( clinical$timeRFS <5*365) ] <-"poor"clinical$survivalgroup[which( clinical$timeRFS >10*365) ] <-"good"## To visualise it in a single graphhatchingPlot( data_sce,region ="region",imageID ="metabricId",cellType ="description",spatialCoords =c("Location_Center_X", "Location_Center_Y") )
We have written a small function draw_region_clustering_result to visualise the data separately in multiple graphs for the individual MB-0002.
If we’d like to examine multiple individuals, is there a better visualise this information?
To perform spatial region detection across multiple individuals, we will use the Banksy package. The Banksy method identifies spatial regions by integrating both gene expression and spatial information. Below is the code to perform spatial region detection using Banksy and visualise the results.
Code
library(Banksy)library(SpatialExperiment)library(cowplot)### The code below takes a while to run, so we have saved the output as an RDS file that can be loaded directly.# spe_base <- SpatialExperiment(# assays = assays(data_sce),# rowData = rowData(data_sce),# colData = colData(data_sce),# spatialCoords = as.matrix(colData(data_sce)[, c("Location_Center_X", "Location_Center_Y")])# )# assay_to_use <- if ("counts" %in% assayNames(spe_base)) "counts" else assayNames(spe_base)[1]# sample_ids <- unique(spe_base$metabricId)# spe_list <- lapply(sample_ids, function(id) {# spe_subset <- spe_base[, spe_base$metabricId == id]# Banksy::computeBanksy(# spe_subset,# assay_name = assay_to_use,# k_geom = 18# )# })# spe_joint <- do.call(cbind, spe_list)# spe_joint <- Banksy::runBanksyPCA(# spe_joint,# lambda = 0.8,# npcs = 30,# group = "metabricId"# )# spe_joint <- Banksy::clusterBanksy(# spe_joint,# lambda = 0.8,# npcs = 30,# algo = "kmeans",# kmeans.centers = 5# )# clust_col <- grep("^clust_", colnames(colData(spe_joint)), value = TRUE)[1]# data_sce$banksy_region <- colData(spe_joint)[colnames(data_sce), clust_col]banksy_ouput_spe_joint <-readRDS("../data/banksy_ouput_spe_joint.rds")hatchingPlot( banksy_ouput_spe_joint,region ="banksy_region",imageID ="metabricId",cellType ="description",spatialCoords =c("Location_Center_X", "Location_Center_Y"))
We can also visualise the spatial region detection results across multiple individuals side by side for comparison. Below is an example of visualising two individuals, MB-0002 and MB-0064. Notice the different cell type compositions in the different regions between the two individuals.
We can better interpret the region output by summarising the proportion of each cell-type in a region across the individuals. We look at the composition of cell-types in each region and compare between prognostic outcomes.
`summarise()` has grouped output by 'Var1', 'Var2'. You can override using the
`.groups` argument.
The number of sub-cell types increase considerably when we want to add spatial domain (region) information. To enhance clarity and facilitate understanding, it may be helpful to choose a predetermined region. The code generates a set of boxplots that enable the comparison of cell-type proportions between individuals with good and poor prognosis in region_5.
Here, we use scFeatures to generate molecular features for each individual using the features x cells matrices. These features are interpretable and can be used for downstream analyses. In general, scFeatures generates features across six categories:
Cell-type proportions.
Cell-type specific gene expressions.
Cell-type specific pathway expressions.
Cell-type interaction scores.
Aggregated gene expressions.
Spatial metrics: Nearest neighbour’s correlation, L statistics, and Moran’s I.
The different types of features constructed enable a more comprehensive multi-view understanding of each individual. By default, the function will generate a total of 13 different types of features and are stored as 13 samples x features matrices, one for each feature type.
In this section, we will examine spatial information from two perspectives. Utilising spatial domain detection described in part 3, we select a specific region of interest and create molecular representations of that region for each individual (Section 4.1). Second, we will utilise spatial statistics to capture spatial relationships within the region of interest, such as Moran’s I (Section 4.2)
Code
region <- data_sce$region# Define a series of sub-cell-types that is regional specificdata_sce$celltype <-paste0( data_sce$description , "-" , region)
4.1 How to create molecular representations of individuals for an ROI?
Here, we can consider regional information and the following code allows us to create cell-type specific features for each region. We use the function “paste0” to construct region-specific sub cell-types and name it as celltype in the R object data_sce. For simplicity, in this workshop, the variable celltype in the R object data_sce refers to region-specific sub-cell-types.
There are a total of 13 different types of features (feature_types) that you can choose to generate. The argument type refers the type of input data we have. This is either scrna (single-cell RNA-sequencing data), spatial_p (spatial proteomics data), or spatial_t (single-cell spatial data). In this section, we will ignore spatial information and generate non-spatial features, such as pseudobulking at the sample /cell type levels, overall expression, cell type proportions etc…
Suppose that we are interested in determining the proportion of cell-types within each region for each individual. It is necessary to specify type = spatial_p to reflect that we have spatial proteomics data and feature_types = proportion_raw to indicate we intend to calculate cell-type proportions for each of the region-specific cell-types. Suppose we are only interested in the molecular representation of HR- CK7+ within individuals within region 5.
Code
## [A] The next few lines extract specific information from data_sce as input to scFeatures. ## Extract the expression matrix from data_sceIMCmatrix <-assay(data_sce, "norm")## Extract the sample information ## append the condition to the individuals so we can easily retrieve the individuals condition sample <- data_sce$metabricId cond <- clinical[match(sample, clinical$metabricId), ]$survivalgroupsample <-paste0(sample, "_cond_", cond ) ## Extract the region-specific sub-cell-typescelltype <- data_sce$celltype## Extract the spatial coordinatesspatialCoords <-list(colData(data_sce)$Location_Center_X, colData(data_sce)$Location_Center_Y)### [B] Running scFeaturesscfeatures_result <-scFeatures(IMCmatrix, sample = sample, celltype = celltype, spatialCoords = spatialCoords,feature_types ="proportion_raw", type ="spatial_p" )
Question
Are there any regions that are associated with “good” and “poor” prognosis?
Is this the right way to visualise the results?
Code
## [A] The next few lines extract specific information from data_sce as input to scFeatures. ## Select the HR- CK7+-region sub-cell-type # There are different ways you can use `scFeatures` to generate molecular representations for individuals and it requires the following information for spatial data.# # data,\# sample,# X coordinates,# Y coordinates,# feature_types, and# typeindex <-grep("HR- CK7+-region" , data_sce$celltype, fixed=TRUE)selected_data <- IMCmatrix[, index]selected_sample <- sample[index]selected_celltype <- data_sce$celltype[ index] selected_spatialCoords <-list(colData(data_sce)$Location_Center_X[index], colData(data_sce)$Location_Center_Y[index])### [B] Running scFeaturesscfeatures_result <-scFeatures( selected_data, sample = selected_sample, celltype = selected_celltype,spatialCoords = selected_spatialCoords,feature_types ="proportion_raw", type ="spatial_p" )### [C] Visualize the regional composition makeup for each individual for HR- CK7+ and HR- CK7-feature <- scfeatures_result$proportion_rawfeature <- feature[ grep("poor|good", rownames(feature)), ]plot_barplot(feature ) +ggtitle("Proportion raw feature") +labs(y ="proportion\n")
Can we identify “differential expression” for a feature of interest? The R object scfeatures_result contains a variety of features. A important question focuses on the identification of features that reflect an association with the prognostic outcome, specifically distinguishing between good and poor outcomes. The code provided below demonstrates the use of the limma() function to fit a linear model for the purpose of analysing gene_mean_celltype as an illustration feature type. The feature type known as gene_mean_celltype represents the mean protein expression for each sub-cell-type specific to a spatial region. It is a matrix consisting of 77 individuals and 4180 features. It is important to acknowledge that within the context of our IMC data, the term “gene” is used to refer to “protein”.
Code
scfeatures_result <-readRDS("../data/scfeatures_result.RDS")# Extract cell-type specific gene expression for all regions. gene_mean_celltype <- scfeatures_result$gene_mean_celltype# Extract HR+ CK7 cell-type specific gene expression for Region5index <-grep("HR+ CK7--region5", colnames(gene_mean_celltype) , fixed= T)gene_mean_celltype <- gene_mean_celltype [, index] # transpose to ensure we have gene by sample matrixgene_mean_celltype <-t(gene_mean_celltype)# Extract the two conditions of interest - poor prognosis vs good prognosiscondition <-unlist( lapply( strsplit( colnames(gene_mean_celltype) , "_cond_"), `[`, 2))condition <-data.frame(sample =colnames(gene_mean_celltype), condition = condition )select_index <-which( condition$condition %in%c("poor", "good" ))condition <- condition[ select_index, ]gene_mean_celltype<- gene_mean_celltype [ , select_index]# Calculate log fold change each protein using limmadesign <-model.matrix(~condition, data = condition)fit <-lmFit(gene_mean_celltype, design)fit <-eBayes(fit)tT <-topTable(fit, n =Inf)tT$gene <-rownames(tT)tT[1:6] <-signif(tT[1:6], digits=2)DT::datatable(tT)
We visualise the comparison using a volcano plot and a dotplot for the cell-type specific expression feature. This is a type of scatter-plot that is used to quickly identify changes in large datasets and represent the significance (y-axis) versus effect size or fold-change (x-axis).
Code
# order the proteins by log fold change tT <- tT[ order(tT$logFC, decreasing = T), ]tT <- tT[1:20, ]ggplot( tT , aes( y =reorder(gene, logFC) , x = logFC ) )+geom_point(aes(colour=-log10(P.Value)), alpha=2/3, size=4) +scale_colour_gradient(low="blue",high="red") +theme_bw() +xlab("logFC") +ylab("region specific cell type specfic features" )
The code below enable us to generate all feature types for all cell-types in a line. Due to limitations with today’s computational capacity, Please DO NOT run it in today’s workshop, it will crash your system.
Code
# here, we specify that this is a spatial proteomics data# scFeatures support parallel computation to speed up the process # scfeatures_result <- scFeatures(IMCmatrix, # type = "spatial_p",# sample = sample, # celltype = celltype, # spatialCoords = spatialCoords,# ncores = 32)
Assuming you have already generated a collection of molecular representation for individuals, please load the prepared RDS file scfeatures_result.RDS. Again, you can remind yourself that all generated feature types are stored in a matrix of samples x features.
Code
# Upload pre-generated RDS filescfeatures_result <-readRDS("../data/scfeatures_result.RDS")# What are the features and the dimensions of features matrices that we have generated?lapply(scfeatures_result, dim)
4.2 How can we represent spatial information and relationships for a given individual?
We will now look at the spatial statistic output by scFeatures and qualitively assess whether there is any association between between these statistics and the “good” and “bad” prognosis groups.
Questions
What kind of information would spatial statistics provide?
Are there any differences in the distribution of spatial statistics between the “good” and “bad” prognosis groups?
The L function is a spatial statistic used to assess the spatial distribution of cell-types. It assesses the significance of cell-cell interactions, by calculating the density of a cell-type with other cell-types within a certain radius. High values indicate spatial association (co-localisation), low values indicate spatial avoidance. To demonstrate the L-function, we will plot a specific patient “MB-0128” who has a high L value for B cells interacting with Fibroblasts and a low L value for B cells interacting with HR- CK7- cells.
Code
tableau_palette <-scale_colour_tableau( palette ="Tableau 20")color_codes <- tableau_palette$palette(10)# Create a named color vectorcell_colors <-c("B cells"= color_codes[1],"Fibroblasts"= color_codes[2],"HR- CK7-"= color_codes[9],"others"= color_codes[4])# select one patient one_sample <- data_sce[ , data_sce$metabricId =="MB-0128" ]one_sample <-data.frame( colData(one_sample) )one_sample$celltype <- one_sample$description# High L-function value plot.one_sample <-data.frame(colData(data_sce[, data_sce$metabricId =="MB-0128"]))one_sample$celltype <- one_sample$descriptionindex <- one_sample$celltype %in%c("B cells", "Fibroblasts")one_sample$celltype[!index] <-"others"a <-ggplot(one_sample, aes(x = Location_Center_X, y = Location_Center_Y, colour = celltype)) +geom_point() +theme(panel.background=element_blank(),axis.line=element_line(color="black")) +scale_colour_manual(values = cell_colors) +# Use named vectorggtitle("Patient MB-0128 - high L value with \n B cells interacting Fibroblasts")# Low L-function value plot.one_sample$celltype <- one_sample$descriptionindex <- one_sample$celltype %in%c("Fibroblasts", "HR- CK7-")one_sample$celltype[!index] <-"others"b <-ggplot(one_sample, aes(x = Location_Center_X, y = Location_Center_Y, colour = celltype)) +geom_point() +theme(panel.background=element_blank(),axis.line=element_line(color="black")) +scale_colour_manual(values = cell_colors) +# Use named vectorggtitle("Patient MB-0128 - low L value with \n B cells interacting HR_ CK7")ggarrange(plotlist =list(a, b))
Moran’s I is a spatial autocorrelation statistic based on both location and values. It quantifies whether similar values tend to occur near each other or are dispersed.
Part 5: Disease classification with ClassifyR [Presentation]
Recurrence risk estimation is a fundamental concern in medical research, particularly in the context of patient survival analysis. In this section, we will estimate recurrence risk using the molecular representation of individuals generated from scFeatures to build a survival model. We will use classifyR to build the survival model. The patient outcome is time-to-event, so, by default, ClassifyR will use Cox proportional hazards ranking to choose a set of features and also Cox proportional hazards to predict risk scores. We will also demonstrate other available models in ClassifyR.
Questions
Are spatial features globally informative in predicting survival?
If not, for which individuals is it important in predicting survival?
Recall in the previous section that we have stored the 13 matrices of different feature types in a list. Instead of manually retrieving each matrix from the list to build separate models, classifyR can directly take a list of matrices as an input and run a repeated cross-validation model on each matrix individually. Below, we run 5 repeats of 5-fold cross-validation. A high score indicates prognosis of a worse outcome than a lower risk score. Although we have provided the code below, to save time, just load the prepared RDS file classifyr_result_IMC.rds and we will focus on the interpretation in this workshop.
# We use the following variables:
# timeRFS: "Time to Recurrence-Free Survival." It is the time period until recurrence occurs.
# eventRFS: "Event in Recurrence-Free Survival."It indicates whether the event has occurred.
# Breast.Tumour.Laterality: Laterality of tumors, eg, whether the tumor is located in left or right.
# ER.Status: Whether the tumor is ER positive or ER negative.
# Inferred.Menopausal.State: of the patient.
# Grade: of the tumor.
# Size: of the tumor.
usefulFeatures <- c("Breast.Tumour.Laterality", "ER.Status", "Inferred.Menopausal.State", "Grade", "Size")
nFeatures <- append(list(clinical = 1:3), lapply(scfeatures_result, function(metaFeature) 1:5))
clinicalAndOmics <- append(list(clinical = clinical), scfeatures_result)
### generate classfyr result
classifyr_result_IMC <- crossValidate(clinicalAndOmics, c("timeRFS", "eventRFS"),
extraParams = list(prepare = list(useFeatures = list(clinical = usefulFeatures))),
nFeatures = nFeatures, nFolds = 5, nRepeats = 5, nCores = 5)
Cox proportional hazards is a classical statistical method, as opposed to machine learning methods like Random survival forest. These machine learning methods can build remarkably complex relationships between features, however their running time can be much longer than Cox proportional hazards. We use feature selection to limit the number of features considered to at most 100 per metafeature and to save time, you can just load the prepared RDS file. We will compare the predictive performance between these methods.
To examine the distribution of prognostic performance, use performancePlot. Currently, the only metric for time-to-event data is C-index and that will automatically be used because the predictive model type is tracked inside of the result objects.
Code
## Make axis label 45 degree to improve readiability tilt <-theme(axis.text.x =element_text(angle =45, vjust =1, hjust =1))## Putting two sets of cross-validation results togethermultiresults <-append(classifyr_result_IMC, survForestCV)ordering <-c("clinical", names(scfeatures_result))performancePlot(multiresults,characteristicsList =list(x ="Assay Name", row ="Classifier Name"),orderingList =list("Assay Name"= ordering)) + tilt
Note how the resultant plot is a ggplot2 object and can be further modified. The same code could be used for a categorical classifier because the random forest implementation provided by the ranger package has the same interface for both. We will examine feature selection stability with selectionPlot.
Does each individual require a different collection of features? Using samplesMetricMap compare the per-sample C-index for Cox models for all kinds of metafeatures.
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.8769]
Appendix
Explanation of spatial features
L function:
The L function is a spatial statistic used to assess the spatial distribution of cell-types. It assesses the significance of cell-cell interactions, by calculating the density of a cell-type with other cell-types within a certain radius. High values indicate spatial association (co-localisation), low values indicate spatial avoidance.
Code
tableau_palette <-scale_colour_tableau( palette ="Tableau 20")color_codes <- tableau_palette$palette(10)# select one patient one_sample <- data_sce[ , data_sce$metabricId =="MB-0128" ]one_sample <-data.frame( colData(one_sample) )one_sample$celltype <- one_sample$description# select certain cell types to examine the interaction index <- one_sample$celltype %in%c("B cells", "Fibroblasts")one_sample$celltype[!index] <-"others"a <-ggplot( one_sample, aes(x = Location_Center_X , y = Location_Center_Y, colour = celltype )) +geom_point() +scale_colour_manual(values = color_codes) +ggtitle( "Patient MB-0128 - high L value with \n B cells interacting Fibroblasts")one_sample$celltype <- one_sample$descriptionindex <- one_sample$celltype %in%c("Fibroblasts", "HR- CK7-")one_sample$celltype[!index] <-"others"b <-ggplot( one_sample, aes(x = Location_Center_X , y = Location_Center_Y, colour = celltype )) +geom_point() +scale_colour_manual(values = color_codes) +ggtitle( "Patient MB-0128 - low L value with \n B cells interacting HR_ CK7")ggarrange(plotlist =list(a,b))
Cell type interaction composition:
We calculate the nearest neighbours of each cell and then calculate the pairs of cell-types based on the nearest neighbour. This allows us to summarise it into a cell-type interaction composition.
Moran’s I is a spatial autocorrelation statistic based on both location and values. It quantifies whether similar values tend to occur near each other or are dispersed.
Code
high <- data_sce["Ki67", data_sce$metabricId =="MB-0132" ]high_meta <-data.frame( colData(high) ) high_meta$expression <-as.vector(logcounts( high)) low <- data_sce["Ki67", data_sce$metabricId =="MB-0249" ]low_meta <-data.frame( colData(low) )low_meta$expression <-as.vector(logcounts(low))a <-ggplot(high_meta, aes(x = Location_Center_X , y = Location_Center_Y, colour =expression)) +geom_point(alpha=0.5) +scale_colour_viridis_c() +ggtitle("Patient MB-0132 - high Moran's I in Ki67")b <-ggplot(low_meta, aes(x = Location_Center_X , y = Location_Center_Y, colour =expression)) +geom_point(alpha=0.5) +scale_colour_viridis_c() +ggtitle("Patient MB-0249 - low Moran's I in Ki67")ggarrange(plotlist =list(a,b))
Nearest Neighbor Correlation:
This metric measures the correlation of proteins/genes between a cell and its nearest neighbour cell.
Code
plot_nncorrelation <-function(thissample , thisprotein){ sample_name <- thissample thissample <- data_sce[, data_sce$metabricId == sample_name] exprsMat <-logcounts(thissample)# calculate NN correlation cell_points_cts <- spatstat.geom::ppp(x =as.numeric(thissample$Location_Center_X ), y =as.numeric(thissample$Location_Center_Y),check =FALSE,xrange =c(min(as.numeric(thissample$Location_Center_X)),max(as.numeric(thissample$Location_Center_X)) ),yrange =c(min(as.numeric(thissample$Location_Center_Y)),max(as.numeric(thissample$Location_Center_Y)) ),marks =t(as.matrix(exprsMat)) ) value <- spatstat.explore::nncorr(cell_points_cts)["correlation", ] value <- value[ thisprotein]# Find the indices of the two nearest neighbors for each cell nn_indices <-nnwhich(cell_points_cts, k =1) protein <- thisprotein df <-data.frame(thiscell_exprs = exprsMat[protein, ] , exprs = exprsMat[protein,nn_indices ]) p <-ggplot(df, aes( x =thiscell_exprs , y = exprs , colour = exprs )) +geom_point(alpha =0.3) +ggtitle(paste0( "Patient ", sample_name , " nn_corr = " , round(value, 2) )) +scale_colour_viridis_c() +xlab("This cell expression") +ylab("Neighbouring cell expression")return (p ) }p1 <-plot_nncorrelation("MB-0605", "HER2")p2 <-plot_nncorrelation("MB-0258", "HER2")ggarrange(plotlist =list(p1, p2))
The authors thank all their colleagues, particularly at The University of Sydney, Sydney Precision Data Science and Charles Perkins Centre for their support and intellectual engagement. Special thanks to Ellis Patrick, Shila Ghazanfar, Andy Tran, Helen, and Daniel for guiding and supporting the building of this workshop.
