Analyze MIBI-TOF image data
This tutorial shows how to apply Squidpy to MIBI-TOF data.
The data used here comes from a recent paper from [Hartmann et al., 2020].
We provide a pre-processed subset of the data, in anndata.AnnData format.
For details on how it was pre-processed, please refer to the original paper.
See also
See Analyze Visium H&E data for additional analysis using images and Analyze seqFISH data for analysis using spatial graph functions.
Import packages & data
To run the notebook locally, create a conda environment as conda env create -f environment.yml using this environment.yml.
import squidpy as sq
import numpy as np
import matplotlib.pyplot as plt
adata = sq.datasets.mibitof()
The subset of the data we consider here comprises three biopsies colorectal carcinoma biopsies from different donors, where MIBI-TOF was used to measure single-cell metabolic profiles. As imaging information, we included three raw image channels:
145_CD45 - a immune cell marker (cyan).
174_CK - a tumor marker (magenta).
113_vimentin - a mesenchymal cell marker (yellow).
and a cell segmentation mask provided by the authors of the original paper.
The adata object contains three different libraries, one for each biopsy.
The images are contained in adata.uns['spatial'][<library_id>]['images'].
Let us visualize the cluster annotations for each library using squidpy.pl.spatial_scatter().
sq.pl.spatial_segment(adata, color="Cluster", library_key="library_id", seg_cell_id="cell_id")
Let us create an ImageContainer from the images contained in adata. As all three biopsies are already joined in adata, let us also create one ImageContainer for all three biopsies using a z-stack. For more information on how to use ImageContainer with z-stacks, also have a look at Use z-stacks with ImageContainer.
imgs = []
for library_id in adata.uns["spatial"].keys():
img = sq.im.ImageContainer(adata.uns["spatial"][library_id]["images"]["hires"], library_id=library_id)
img.add_img(adata.uns["spatial"][library_id]["images"]["segmentation"], library_id=library_id, layer="segmentation")
img["segmentation"].attrs["segmentation"] = True
imgs.append(img)
img = sq.im.ImageContainer.concat(imgs)
Note that we also added the segmentation as an additional layer to img, and set the segmentation attribute in the ImageContainer. This allows visualization of the segmentation layer as a labels layer in Napari.
img
If you have Napari installed, you can have a look at the data using the interactive viewer: Note that you can load the segmentation layer as an overlay over the image.
img.interactive(adata, library_key='library_id')
Let us also statically visualize the data in img, using squidpy.im.ImageCntainer.show():
img.show("image")
img.show("image", segmentation_layer="segmentation")
In the following we show how to use Squidpy to extract cellular mean intensity information using raw images
and a provided segmentation mask.
In the present case, adata of course already contains the post-processed cellular mean intensity
for the raw image channels.
The aim of this tutorial, however, is to showcase how the extraction of such features is possible using Squidpy.
As Squidpy is backed by dask and supports chunked image processing,
also large images can be processed in this way.
Convert image to CMYK
As already mentioned, the images contain information from three raw channels, 145_CD45, 174_CK, and 113_vimentin. As the channel information is encoded in CMYK space, we first need to convert the RGB images to CMYK.
For this, we can use squidpy.im.ImageContainer.apply().
def rgb2cmyk(arr):
"""Convert arr from RGB to CMYK color space."""
R = arr[..., 0] / 255
G = arr[..., 1] / 255
B = arr[..., 2] / 255
K = 1 - (np.max(arr, axis=-1) / 255)
C = (1 - R - K) / (1 - K + np.finfo(float).eps) # avoid division by 0
M = (1 - G - K) / (1 - K + np.finfo(float).eps)
Y = (1 - B - K) / (1 - K + np.finfo(float).eps)
return np.stack([C, M, Y, K], axis=3)
img.apply(rgb2cmyk, layer="image", new_layer="image_cmyk", copy=False)
img.show("image_cmyk", channelwise=True)
Extract per-cell mean intensity
Now that we have disentangled the individual channels, let use use the provided segmentation mask to extract per-cell mean intensities.
By default, the segmentation feature extractor extracts information using all segments (cells) in the current crop. As we would like to only get information of the segment (cell) in the center of the current crop, let us use a custom feature extractor.
Fist, define a custom feature extraction function. This function needs to get the segmentation mask
and the original image as input.
We will achieve this by passing an additional_layers argument to the custom feature extractor.
This special argument will pass the values of every layer in additional_layers
to the custom feature extraction function.
def segmentation_image_intensity(arr, image_cmyk):
"""
Calculate per-channel mean intensity of the center segment.
arr: the segmentation
image_cmyk: the raw image values
"""
import skimage.measure
# the center of the segmentation mask contains the current label
# use that to calculate the mask
s = arr.shape[0]
mask = (arr == arr[s // 2, s // 2, 0, 0]).astype(int)
# use skimage.measure.regionprops to get the intensity per channel
features = []
for c in range(image_cmyk.shape[-1]):
feature = skimage.measure.regionprops_table(
np.squeeze(mask), # skimage needs 3d or 2d images, so squeeze excess dims
intensity_image=np.squeeze(image_cmyk[:, :, :, c]),
properties=["mean_intensity"],
)["mean_intensity"][0]
features.append(feature)
return features
Now, use squidpy.im.calculate_image_features() with the custom feature extractor,
specifying the function (func) to use, and the additional layers (additional_layers)
to pass to the function.
We will use spot_scale = 10 to ensure that we also cover big segments fully by one crop.
sq.im.calculate_image_features(
adata,
img,
library_id="library_id",
features="custom",
spot_scale=10,
layer="segmentation",
features_kwargs={"custom": {"func": segmentation_image_intensity, "additional_layers": ["image_cmyk"]}},
)
Out:
0%| | 0/3309 [00:00<?, ?/s]
0%| | 1/3309 [00:05<4:55:43, 5.36s/]
1%| | 17/3309 [00:05<12:40, 4.33/s]
1%| | 32/3309 [00:05<05:43, 9.53/s]
1%|1 | 49/3309 [00:05<03:08, 17.31/s]
2%|1 | 66/3309 [00:05<01:59, 27.13/s]
3%|2 | 83/3309 [00:05<01:22, 39.00/s]
3%|3 | 100/3309 [00:05<01:00, 52.88/s]
4%|3 | 118/3309 [00:06<00:46, 69.20/s]
4%|4 | 135/3309 [00:06<00:37, 85.05/s]
5%|4 | 152/3309 [00:06<00:31, 100.52/s]
5%|5 | 169/3309 [00:06<00:27, 114.18/s]
6%|5 | 186/3309 [00:06<00:24, 126.53/s]
6%|6 | 203/3309 [00:06<00:22, 135.61/s]
7%|6 | 220/3309 [00:06<00:21, 143.07/s]
7%|7 | 237/3309 [00:06<00:20, 150.23/s]
8%|7 | 254/3309 [00:06<00:19, 154.10/s]
8%|8 | 271/3309 [00:07<00:19, 158.16/s]
9%|8 | 289/3309 [00:07<00:18, 162.97/s]
9%|9 | 306/3309 [00:07<00:18, 164.49/s]
10%|9 | 324/3309 [00:07<00:17, 167.86/s]
10%|# | 342/3309 [00:07<00:17, 170.17/s]
11%|# | 360/3309 [00:07<00:17, 172.05/s]
11%|#1 | 378/3309 [00:07<00:16, 172.98/s]
12%|#1 | 396/3309 [00:07<00:16, 174.73/s]
13%|#2 | 414/3309 [00:07<00:16, 174.55/s]
13%|#3 | 432/3309 [00:07<00:16, 175.12/s]
14%|#3 | 450/3309 [00:08<00:16, 174.17/s]
14%|#4 | 468/3309 [00:08<00:16, 174.91/s]
15%|#4 | 486/3309 [00:08<00:16, 175.64/s]
15%|#5 | 504/3309 [00:08<00:15, 175.83/s]
16%|#5 | 522/3309 [00:08<00:15, 175.65/s]
16%|#6 | 540/3309 [00:08<00:15, 174.78/s]
17%|#6 | 558/3309 [00:08<00:15, 174.77/s]
17%|#7 | 577/3309 [00:08<00:15, 176.44/s]
18%|#7 | 595/3309 [00:08<00:15, 174.93/s]
19%|#8 | 613/3309 [00:08<00:15, 169.05/s]
19%|#9 | 630/3309 [00:09<00:16, 165.18/s]
20%|#9 | 647/3309 [00:09<00:16, 165.57/s]
20%|## | 665/3309 [00:09<00:15, 167.36/s]
21%|## | 682/3309 [00:09<00:16, 154.77/s]
21%|##1 | 698/3309 [00:09<00:17, 152.13/s]
22%|##1 | 715/3309 [00:09<00:16, 155.77/s]
22%|##2 | 733/3309 [00:09<00:15, 161.17/s]
23%|##2 | 751/3309 [00:09<00:15, 164.73/s]
23%|##3 | 768/3309 [00:09<00:15, 165.33/s]
24%|##3 | 786/3309 [00:10<00:15, 167.87/s]
24%|##4 | 804/3309 [00:10<00:14, 168.99/s]
25%|##4 | 822/3309 [00:10<00:14, 169.34/s]
25%|##5 | 839/3309 [00:10<00:14, 169.03/s]
26%|##5 | 857/3309 [00:10<00:14, 171.23/s]
26%|##6 | 875/3309 [00:10<00:14, 170.72/s]
27%|##6 | 893/3309 [00:10<00:14, 171.00/s]
28%|##7 | 911/3309 [00:10<00:14, 170.58/s]
28%|##8 | 929/3309 [00:10<00:13, 170.82/s]
29%|##8 | 947/3309 [00:10<00:13, 171.13/s]
29%|##9 | 965/3309 [00:11<00:13, 171.83/s]
30%|##9 | 983/3309 [00:11<00:13, 171.83/s]
30%|### | 1001/3309 [00:11<00:13, 172.49/s]
31%|### | 1019/3309 [00:11<00:13, 172.44/s]
31%|###1 | 1037/3309 [00:11<00:13, 172.72/s]
32%|###1 | 1055/3309 [00:11<00:13, 172.86/s]
32%|###2 | 1073/3309 [00:11<00:12, 172.60/s]
33%|###2 | 1091/3309 [00:11<00:12, 173.22/s]
34%|###3 | 1109/3309 [00:11<00:12, 173.40/s]
34%|###4 | 1127/3309 [00:12<00:12, 173.79/s]
35%|###4 | 1145/3309 [00:12<00:12, 173.75/s]
35%|###5 | 1163/3309 [00:12<00:12, 173.77/s]
36%|###5 | 1181/3309 [00:12<00:12, 170.88/s]
36%|###6 | 1199/3309 [00:12<00:12, 170.00/s]
37%|###6 | 1217/3309 [00:12<00:12, 165.61/s]
37%|###7 | 1234/3309 [00:12<00:12, 166.74/s]
38%|###7 | 1251/3309 [00:12<00:12, 166.80/s]
38%|###8 | 1268/3309 [00:12<00:12, 164.63/s]
39%|###8 | 1285/3309 [00:12<00:12, 163.77/s]
39%|###9 | 1302/3309 [00:13<00:12, 164.30/s]
40%|###9 | 1320/3309 [00:13<00:11, 167.19/s]
40%|#### | 1337/3309 [00:13<00:11, 166.89/s]
41%|#### | 1354/3309 [00:13<00:11, 167.49/s]
41%|####1 | 1371/3309 [00:13<00:11, 167.57/s]
42%|####1 | 1389/3309 [00:13<00:11, 169.03/s]
43%|####2 | 1407/3309 [00:13<00:11, 170.58/s]
43%|####3 | 1425/3309 [00:13<00:10, 171.37/s]
44%|####3 | 1443/3309 [00:13<00:10, 170.19/s]
44%|####4 | 1461/3309 [00:14<00:10, 170.79/s]
45%|####4 | 1479/3309 [00:14<00:10, 169.73/s]
45%|####5 | 1497/3309 [00:14<00:10, 170.51/s]
46%|####5 | 1515/3309 [00:14<00:10, 170.69/s]
46%|####6 | 1533/3309 [00:14<00:10, 170.01/s]
47%|####6 | 1551/3309 [00:14<00:10, 169.60/s]
47%|####7 | 1568/3309 [00:14<00:10, 169.58/s]
48%|####7 | 1585/3309 [00:14<00:10, 169.68/s]
48%|####8 | 1603/3309 [00:14<00:10, 170.55/s]
49%|####8 | 1621/3309 [00:14<00:09, 170.19/s]
50%|####9 | 1639/3309 [00:15<00:09, 170.42/s]
50%|##### | 1657/3309 [00:15<00:09, 171.61/s]
51%|##### | 1675/3309 [00:15<00:09, 171.84/s]
51%|#####1 | 1693/3309 [00:15<00:09, 171.74/s]
52%|#####1 | 1711/3309 [00:15<00:09, 172.29/s]
52%|#####2 | 1729/3309 [00:15<00:09, 172.32/s]
53%|#####2 | 1747/3309 [00:15<00:09, 172.79/s]
53%|#####3 | 1765/3309 [00:15<00:08, 171.75/s]
54%|#####3 | 1783/3309 [00:15<00:08, 172.00/s]
54%|#####4 | 1801/3309 [00:15<00:08, 173.07/s]
55%|#####4 | 1819/3309 [00:16<00:08, 172.91/s]
56%|#####5 | 1837/3309 [00:16<00:08, 172.84/s]
56%|#####6 | 1855/3309 [00:16<00:08, 173.83/s]
57%|#####6 | 1873/3309 [00:16<00:08, 173.76/s]
57%|#####7 | 1891/3309 [00:16<00:08, 171.88/s]
58%|#####7 | 1909/3309 [00:16<00:08, 171.71/s]
58%|#####8 | 1927/3309 [00:16<00:07, 173.04/s]
59%|#####8 | 1945/3309 [00:16<00:07, 173.01/s]
59%|#####9 | 1963/3309 [00:16<00:07, 172.65/s]
60%|#####9 | 1981/3309 [00:17<00:07, 173.07/s]
60%|###### | 1999/3309 [00:17<00:07, 172.09/s]
61%|###### | 2017/3309 [00:17<00:07, 172.93/s]
61%|######1 | 2035/3309 [00:17<00:07, 171.61/s]
62%|######2 | 2053/3309 [00:17<00:07, 172.03/s]
63%|######2 | 2071/3309 [00:17<00:07, 172.06/s]
63%|######3 | 2089/3309 [00:17<00:07, 172.57/s]
64%|######3 | 2107/3309 [00:17<00:06, 173.57/s]
64%|######4 | 2125/3309 [00:17<00:06, 173.74/s]
65%|######4 | 2143/3309 [00:17<00:06, 171.90/s]
65%|######5 | 2161/3309 [00:18<00:06, 172.40/s]
66%|######5 | 2179/3309 [00:18<00:06, 172.68/s]
66%|######6 | 2197/3309 [00:18<00:06, 174.30/s]
67%|######6 | 2215/3309 [00:18<00:06, 173.35/s]
67%|######7 | 2233/3309 [00:18<00:06, 169.64/s]
68%|######7 | 2250/3309 [00:18<00:06, 163.37/s]
69%|######8 | 2267/3309 [00:18<00:06, 162.51/s]
69%|######9 | 2284/3309 [00:18<00:06, 163.55/s]
70%|######9 | 2302/3309 [00:18<00:06, 165.64/s]
70%|####### | 2319/3309 [00:19<00:06, 164.79/s]
71%|####### | 2336/3309 [00:19<00:05, 165.90/s]
71%|#######1 | 2353/3309 [00:19<00:05, 165.43/s]
72%|#######1 | 2370/3309 [00:19<00:05, 164.74/s]
72%|#######2 | 2388/3309 [00:19<00:05, 166.77/s]
73%|#######2 | 2406/3309 [00:19<00:05, 170.03/s]
73%|#######3 | 2424/3309 [00:19<00:05, 169.53/s]
74%|#######3 | 2441/3309 [00:19<00:05, 161.98/s]
74%|#######4 | 2458/3309 [00:19<00:05, 163.66/s]
75%|#######4 | 2475/3309 [00:19<00:05, 162.99/s]
75%|#######5 | 2493/3309 [00:20<00:04, 166.33/s]
76%|#######5 | 2510/3309 [00:20<00:04, 162.08/s]
76%|#######6 | 2527/3309 [00:20<00:04, 160.52/s]
77%|#######6 | 2545/3309 [00:20<00:04, 164.76/s]
77%|#######7 | 2563/3309 [00:20<00:04, 166.76/s]
78%|#######7 | 2580/3309 [00:20<00:04, 166.79/s]
79%|#######8 | 2598/3309 [00:20<00:04, 168.43/s]
79%|#######9 | 2616/3309 [00:20<00:04, 170.15/s]
80%|#######9 | 2634/3309 [00:20<00:03, 171.15/s]
80%|######## | 2652/3309 [00:21<00:03, 171.70/s]
81%|######## | 2670/3309 [00:21<00:03, 172.69/s]
81%|########1 | 2688/3309 [00:21<00:03, 173.70/s]
82%|########1 | 2706/3309 [00:21<00:03, 174.56/s]
82%|########2 | 2724/3309 [00:21<00:03, 173.93/s]
83%|########2 | 2742/3309 [00:21<00:03, 174.26/s]
83%|########3 | 2760/3309 [00:21<00:03, 174.57/s]
84%|########3 | 2778/3309 [00:21<00:03, 174.76/s]
84%|########4 | 2796/3309 [00:21<00:02, 174.02/s]
85%|########5 | 2814/3309 [00:21<00:02, 174.46/s]
86%|########5 | 2832/3309 [00:22<00:02, 174.68/s]
86%|########6 | 2850/3309 [00:22<00:02, 175.01/s]
87%|########6 | 2868/3309 [00:22<00:02, 176.23/s]
87%|########7 | 2886/3309 [00:22<00:02, 177.18/s]
88%|########7 | 2904/3309 [00:22<00:02, 176.51/s]
88%|########8 | 2922/3309 [00:22<00:02, 175.38/s]
89%|########8 | 2940/3309 [00:22<00:02, 175.01/s]
89%|########9 | 2958/3309 [00:22<00:02, 171.93/s]
90%|########9 | 2976/3309 [00:22<00:01, 171.50/s]
90%|######### | 2994/3309 [00:22<00:01, 172.34/s]
91%|#########1| 3012/3309 [00:23<00:01, 173.74/s]
92%|#########1| 3030/3309 [00:23<00:01, 173.12/s]
92%|#########2| 3048/3309 [00:23<00:01, 172.76/s]
93%|#########2| 3066/3309 [00:23<00:01, 172.84/s]
93%|#########3| 3084/3309 [00:23<00:01, 173.21/s]
94%|#########3| 3102/3309 [00:23<00:01, 174.17/s]
94%|#########4| 3120/3309 [00:23<00:01, 174.47/s]
95%|#########4| 3138/3309 [00:23<00:00, 174.57/s]
95%|#########5| 3156/3309 [00:23<00:00, 175.47/s]
96%|#########5| 3174/3309 [00:24<00:00, 174.17/s]
96%|#########6| 3192/3309 [00:24<00:00, 174.15/s]
97%|#########7| 3210/3309 [00:24<00:00, 174.34/s]
98%|#########7| 3228/3309 [00:24<00:00, 175.80/s]
98%|#########8| 3246/3309 [00:24<00:00, 174.66/s]
99%|#########8| 3264/3309 [00:24<00:00, 172.75/s]
99%|#########9| 3282/3309 [00:24<00:00, 170.44/s]
100%|#########9| 3300/3309 [00:24<00:00, 170.25/s]
100%|##########| 3309/3309 [00:24<00:00, 133.40/s]
The resulting features are stored in adata.obs['img_features'],
with channel 0 representing 145_CD45, channel 1 174_CK, and channel 2 113_vimentin.
adata.obsm["img_features"]
As described in [Hartmann et al., 2020], let us transformed using an inverse hyperbolic sine (arcsinh) co-factor of 0.05, to allow us to compare the computed mean intensities with the values contained in adata.
adata.obsm["img_features_transformed"] = np.arcsinh(adata.obsm["img_features"] / 0.05)
Now, let’s visualize the result:
channels = ["CD45", "CK", "vimentin"]
fig, axes = plt.subplots(1, 3, figsize=(15, 3))
for i, ax in enumerate(axes):
X = np.array(adata[:, channels[i]].X.todense())[:, 0]
Y = adata.obsm["img_features_transformed"][f"segmentation_image_intensity_{i}"]
ax.scatter(X, Y)
ax.set_xlabel("true value in adata.X")
ax.set_ylabel("computed mean intensity")
corr = np.corrcoef(X, Y)[1, 0]
ax.set_title(f"{channels[i]}, corr: {corr:.2f}")
We get high correlations between the original values and our computation using Squidpy. The remaining differences are probably due to more pre-processing applied by the authors of [Hartmann et al., 2020].
In this tutorial we have shown how to pre-process imaging data to extract per-cell
counts / mean intensities using Squidpy.
Of course it is also possible to apply spatial statistics functions provided by the
squidpy.gr module to MIBI-TOF data.
For examples of this, please see our other Analysis tutorials, e.g.
Analyze seqFISH data.
Total running time of the script: ( 1 minutes 2.671 seconds)
Estimated memory usage: 802 MB