Skip to main content
Version: 0.18.x

3D Cell visualization

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

Setup

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

from pathlib import Path
from pixelator import read, simple_aggregate

import plotly.graph_objects as go
import plotly.io as pio
from plotly import subplots

pio.renderers.default = "plotly_mimetype+notebook_connected"
DATA_DIR = Path("<path to the directory to save datasets to>")

Load data

In this tutorial we will continue with the stimulated T cell data that we analyzed in the previous tutorial (Karlsson et al., Nature Methods 2024).

pg_data_combined = read(DATA_DIR/ "uropod_spatial.pxl")

polarity_scores = pg_data_combined.polarization
colocalization_scores = pg_data_combined.colocalization
edgelist = pg_data_combined.edgelist

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.

# Create a new columns with marker_1/marker_2 using iloc
colocalization_scores['contrast'] = colocalization_scores['marker_1'] + '/' + colocalization_scores['marker_2']
colocalization_scores = colocalization_scores[colocalization_scores['contrast'] == 'CD162/CD50']

# Fetch data from pg_data_combined and convert to a DataFrame
marker_counts = pg_data_combined.adata.X[
:, pg_data_combined.adata.var.index.isin(["CD162", "CD37", "CD50"])
]
marker_counts = pd.DataFrame(
marker_counts,
columns=["CD162", "CD37", "CD50"]
)
marker_counts["component"] = pg_data_combined.adata.obs.index

# Merge the colocalization scores and marker count data dataframes on 'component'
colocalization_scores = colocalization_scores.merge(marker_counts, on='component', how='left')

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 35 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 = 35

sns.histplot(colocalization_scores['pearson_z'], bins=50)
plt.axvline(pearson_z_threshold, linestyle="--")
plt.title("Pearson Z scores [CD162/CD50]")
Text(0.5, 1.0, 'Pearson Z scores [CD162/CD50]')

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

plot_components = colocalization_scores[
(colocalization_scores['pearson_z'] > pearson_z_threshold)
].sort_values(by='CD162', ascending=False)

plot_components = plot_components[["component", "pearson_z", "CD162", "CD50", "CD37"]]
plot_components
component pearson_z CD162 CD50 CD37
165 RCVCMP0000082_Stimulated 57.554885 1039 1105 598
156 RCVCMP0000045_Stimulated 39.529305 929 1205 621
179 RCVCMP0000132_Stimulated 42.860281 688 700 444
273 RCVCMP0000568_Stimulated 36.800132 566 462 185
216 RCVCMP0000296_Stimulated 38.464479 486 738 248
220 RCVCMP0000308_Stimulated 47.187758 485 557 347
265 RCVCMP0000535_Stimulated 39.316217 412 651 286
151 RCVCMP0000059_Stimulated 36.645624 302 467 115
342 RCVCMP0001141_Stimulated 37.816523 211 113 41
169 RCVCMP0000122_Stimulated 41.991363 193 241 79
248 RCVCMP0000453_Stimulated 43.902488 192 212 93
339 RCVCMP0001049_Stimulated 35.784675 108 212 75
243 RCVCMP0000431_Stimulated 41.535092 68 118 69

In addition to 2D layouts, we can also create three-dimensional layouts to visualize single cells. The pivot multidimensional scaling (pmds), Kamada–Kawai and Fruchterman-Reingold algorithms that we used in 2D Cell Visualization can also be used to create 3D layouts. With pixelator version >0.17 you can specify to pre-compute layouts which will be stored in the PXL file.

Now we can check if a pre-computed layout exists in our data set and otherwise generate a new 3D layout. Here we will use weighted pmds (wpmds_3d) to generate the 3D layout, which is similar to pmds but improves the results by leveraging edge weights based on transition probabilities in the computation.

Here, we select one of the components from our plot_components table.

component = "RCVCMP0001141_Stimulated"

if (hasattr(pg_data_combined, "precomputed_layouts") and pg_data_combined.precomputed_layouts is not None
):
graph_layout_data = pg_data_combined.precomputed_layouts.filter(
component_ids=component
).to_df()
else:
graph = pg_data_combined.graph(component,
add_node_marker_counts=True,
simplify=True,
use_full_bipartite=True)
graph_layout_data = graph.layout_coordinates("wpmds_3d",
only_keep_a_pixels=True,
get_node_marker_matrix=True)

Now when we have generated the 3D layouts, we can visualize single cells with protein abundance levels projected upon single UPIs on the reconstructed cell surface. The layout can be represented as a point cloud, and if we project the abundance of one of the uropod markers upon it, we see that counts are indeed clustered to a small area of the cell surface. Here we are displaying the log1p-transformed counts for CD162.

fig = go.Figure(
data=[
go.Scatter3d(
x=graph_layout_data["x"],
y=graph_layout_data["y"],
z=graph_layout_data["z"],
mode="markers",
marker=dict(size=3, opacity=1, colorscale=["#cccccc", "#ff4500"]),
marker_color=np.log1p(graph_layout_data["CD162"]),
),
]
)
fig.update_layout(title=component)
fig

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

fig = go.Figure(
data=[
go.Scatter3d(
x=graph_layout_data["x_norm"],
y=graph_layout_data["y_norm"],
z=graph_layout_data["z_norm"],
mode="markers",
marker=dict(size=3, opacity=1, colorscale=["#cccccc", "#ff4500"]),
marker_color=np.log1p(graph_layout_data["CD162"]),
),
]
)
fig.update_layout(title=component)
fig

Viewing colocalization of markers

marker_1 = "CD50"
marker_2 = "CD162"

trace1 = go.Scatter3d(
x=graph_layout_data["x"],
y=graph_layout_data["y"],
z=graph_layout_data["z"],
mode="markers",
marker=dict(size=3, opacity=1, colorscale=["#cccccc", "#ff4500"], color=np.log1p(graph_layout_data[marker_1])),
scene="scene1",
)


# Create the second 3D plot trace
trace2 = go.Scatter3d(
x=graph_layout_data["x"],
y=graph_layout_data["y"],
z=graph_layout_data["z"],
mode="markers",
marker=dict(size=3, opacity=1, colorscale=["#cccccc", "#ff4500"], color=np.log1p(graph_layout_data[marker_2])),
scene="scene2",
)


f = subplots.make_subplots(rows=1, cols=2, specs=[[{'is_3d': True}, {'is_3d': True}]])

f.append_trace(trace1, 1, 1)
f.append_trace(trace2, 1, 2)

fig = go.FigureWidget(f)

def cam_change1(layout, camera):
if fig.layout.scene2.camera == camera:
return
fig.layout.scene2.camera = camera

def cam_change2(layout, camera):
if fig.layout.scene1.camera == camera:
return
fig.layout.scene1.camera = camera

fig.layout.scene1.on_change(cam_change1, 'camera')
fig.layout.scene2.on_change(cam_change2, 'camera')
display(fig)