Skip to main content
Version: 0.17.x

3D Cell visualization

In this tutorial, we will make 3D visualization of the graph components that constitute cells in MPX data using R.

Setup

library(dplyr)
library(tidyr)
library(here)
library(ggplot2)
library(pixelatorR)
library(stringr)
library(SeuratObject)

Load data

We will continue with the stimulated T cell data that we analyzed in the previous tutorial (Karlsson et al., Nature Methods 2024). Here, we read the data files and apply the same filtering steps as in Spatial Analysis.

# The folders the data is located in:
data_files <-
c(Control = here("data/Uropod_control.layout.dataset.pxl"),
Stimulated = here("data/Uropod_CD54_fixed_RANTES_stimulated.layout.dataset.pxl"))


# Read .pxl files as a list of Seurat objects
pg_data <- lapply(data_files, ReadMPX_Seurat)

# Combine data to Seurat object
pg_data_combined <-
merge(pg_data[[1]],
y = pg_data[-1],
add.cell.ids = names(pg_data))

# Add sample column to meta data
pg_data_combined <-
AddMetaData(pg_data_combined,
metadata = str_remove(colnames(pg_data_combined), "_.*"),
col.name = "sample")
rm(pg_data)
gc()
           used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells 3490397 186.5 5342907 285.4 5342907 285.4
Vcells 58727470 448.1 195331305 1490.3 194748656 1485.9
# Filter cells
pg_data_combined <-
pg_data_combined %>%
subset(molecules >= 8000 & tau_type == "normal")

Now we can select components with high colocalization for a specific marker pair to get some good example components to work with. First, we will filter our colocalization scores to keep only CD162/CD50 and then we will add the abundance data for CD162, CD37 and CD50 to the colocalization table.

colocalization_scores <- 
ColocalizationScores(pg_data_combined, meta_data_columns = "sample") %>%
unite(contrast, marker_1, marker_2, sep = "/") %>%
filter(contrast == "CD162/CD50") %>%
left_join(y = FetchData(pg_data_combined, vars = c("CD162", "CD50", "CD37"), layer = "counts") %>%
as_tibble(rownames = "component"),
by = "component")

Next, we will select components with high colocalization scores for CD162/CD50. In the histogram below, we see the distribution of Pearson Z scores for the CD162/CD50 colocalization. Since we want to narrow down the selection of components, we set a threshold for the Pearson Z scores. Here, we use a threshold of 17 to keep the components with the highest colocalization scores. Note that the threshold is arbitrary and should be adjusted based on the data.

# Set Pearson Z threshold
pearson_z_threshold <- 17

ggplot(colocalization_scores, aes(pearson_z)) +
geom_histogram(bins = 50) +
geom_vline(xintercept = pearson_z_threshold, linetype = "dashed") +
labs(title = "Pearson Z scores [CD162/CD50]")

Then we apply our threshold to get a final selection of 7 components.

plot_components <- 
colocalization_scores %>%
filter(pearson_z > pearson_z_threshold, sample == "Stimulated") %>%
arrange(-CD162) %>%
select(component, pearson_z, CD162, CD50, CD37)
plot_components
# A tibble: 7 × 5
component pearson_z CD162 CD50 CD37
<chr> <dbl> <dbl> <dbl> <dbl>
1 Stimulated_RCVCMP0000337 19.6 209 113 40
2 Stimulated_RCVCMP0000219 23.1 153 229 83
3 Stimulated_RCVCMP0001498 18.0 136 111 38
4 Stimulated_RCVCMP0000278 23.5 68 118 69
5 Stimulated_RCVCMP0001465 17.4 59 60 86
6 Stimulated_RCVCMP0001881 17.4 53 102 48
7 Stimulated_RCVCMP0002508 17.2 49 59 84

Visualize single cell graph components in 3D

As one of the main features of the MPX technology, we can create three-dimensional layouts to visualize single cells. The layout algorithms pmds and fr that we used in 2D Cell Visualization can also be used here.

Pre-computed layouts

As we learned in the 2D visualization tutorial, we can use pre-computed layouts if they are available in our PXL files. Alternatively, we use the ComputeLayout method to generate 3D layouts after we have loaded our component graphs with LoadCellGraphs. For the latter option, we need to set dim = 3 to make sure that the layout is computed in 3D:

# Run this to load graphs and compute 3D layouts
pg_data_combined <- pg_data_combined %>%
LoadCellGraphs(cells = plot_components %>%
pull(component)) %>%
ComputeLayout(dim = 3, layout_name = "pmds_3d")

In this tutorial, we load the pre-computed layouts from the PXL files.

pg_data_combined <- pg_data_combined %>% 
LoadCellGraphs(cells = plot_components %>%
pull(component),
force = TRUE,
load_layouts = TRUE)

Now when we have the 3D layouts, we can visualize component graphs as 3D scatter plots where each point corresponds to a UPI. Just as in the 2D plots, we can color the points based on marker abundance to see how it is distributed on the cell surface. However, now that we have a 3D layout, we get an interactive plot that we can rotate, pan and zoom. Since the 3D scatter plots are quite complex, Plot3DGraph only displays the layout for a single component.

From our table above, we can select any component. The “Stimulated_RCVCMP0001498” cell is a good candidate to visualize in 3D as it has a high Pearson Z score for CD162/CD50 and high abundance of CD162, CD50 and CD37. If we project the abundance of one of the uropod markers (CD162) on the nodes, we see that counts are indeed clustered to a small area of the cell surface.

Plot3DGraph(pg_data_combined, 
cell_id = "Stimulated_RCVCMP0001498",
marker = "CD162",
layout_method = "pmds_3d", node_size = 3,
colors = colorRampPalette(colors = c("gray80", "orangered"))(20)) %>%
plotly::layout(scene = list(xaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F),
yaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F),
zaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F)))

We can also normalize the 3D coordinates to a unit sphere.

Plot3DGraph(pg_data_combined, 
cell_id = "Stimulated_RCVCMP0001498",
marker = "CD162",
layout_method = "pmds_3d", node_size = 3,
project = TRUE,
colors = colorRampPalette(colors = c("gray80", "orangered"))(20)) %>%
plotly::layout(scene = list(xaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F),
yaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F),
zaxis = list(title = "", showgrid = F, nticks = 0, showticklabels = F)))

Viewing colocalization of markers

We can also visualize the colocalization of two markers in 3D. Here, we visualize CD162 and CD50 on the same MPX cell component. We use the same approach as above, but now we create two plots and combine them in one view.

# Create a plot for CD162
plotly_obj_cd162 <- Plot3DGraph(pg_data_combined,
cell_id = "Stimulated_RCVCMP0001498",
marker = "CD162",
layout_method = "pmds_3d", node_size = 3,
colors = colorRampPalette(colors = c("gray80", "orangered"))(20),
scene = "scene") %>%
plotly::colorbar(x = -0.1, y = 1) %>%
plotly::layout(title = list(title = ''))
plotly_obj_cd162$x$layout$annotations <- NULL

# Create a plot for CD50
plotly_obj_cd50 <- Plot3DGraph(pg_data_combined,
cell_id = "Stimulated_RCVCMP0001498",
marker = "CD50",
layout_method = "pmds_3d", node_size = 3,
colors = colorRampPalette(colors = c("gray80", "orangered"))(20),
scene = "scene2") %>%
plotly::colorbar(y = 1)
plotly_obj_cd50$x$layout$annotations <- NULL

# Combine plots
main_plot <- plotly::subplot(plotly_obj_cd162,
plotly_obj_cd50, nrows = 1, margin = 0.06) %>%
plotly::layout(scene = list(domain = list(x = c(0, 0.5),
aspectmode = "cube")),
scene2 = list(domain = list(x = c(0.5, 1),
aspectmode = "cube"))
) %>%
plotly::hide_legend()

main_plot$x$layout[1:2] <- NULL

# Add marker annotations
annotations = list(
list(
x = 0.2,
y = 1.0,
text = "CD162",
xref = "paper",
yref = "paper",
xanchor = "center",
yanchor = "bottom",
showarrow = FALSE
),
list(
x = 0.8,
y = 1,
text = "CD50",
xref = "paper",
yref = "paper",
xanchor = "center",
yanchor = "bottom",
showarrow = FALSE
)
)

main_plot <- main_plot %>% plotly::layout(annotations = annotations)

# We can sync the plots with some custom JavaScript code
main_plot %>%
htmlwidgets::onRender(
"function(x, el) {
x.on('plotly_relayout', function(d) {
const camera = Object.keys(d).filter((key) => /\\.camera$/.test(key));
if (camera.length) {
const scenes = Object.keys(x.layout).filter((key) => /^scene\\d*/.test(key));
const new_layout = {};
scenes.forEach(key => {
new_layout[key] = {...x.layout[key], camera: {...d[camera]}};
});
Plotly.relayout(x, new_layout);
}
});
}")