57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411 | def addROI_image (image_path, adata, subset=None,imageid='imageid', overlay=None, flip_y=True,
overlay_category=None,markers=None,channel_names='default',
x_coordinate='X_centroid',y_coordinate='Y_centroid',point_size=10,
point_color=None,seg_mask=None,
n_jobs=-1, verbose=False,
overwrite=True, label='ROI', **kwargs):
"""
Parameters:
image_path : string
Location to the image file (TIFF, OME.TIFF, ZARR supported)
adata : AnnData Object
subset : list, optional
Name of the image to which the ROI is to be added. Note if you have multiple images in the
adata object, you will need to add ROI's to each image one after the other independently.
imageid : string, optional
In the event that the adata object contains multiple images, it is
important that ROIs are added to each image seperately. Pass the name
of the column that contains the `imageid` and use it in conjunction with
the `subset` parameter to add ROI's to a specific image.
seg_mask: string
Location to the segmentation mask file.
overlay : string, optional
Name of the column with any categorical data such as phenotypes or clusters.
flip_y : bool, optional
Flip the Y-axis if needed. Some algorithms output the XY with the Y-coordinates flipped.
If the image overlays do not align to the cells, try again by setting this to `False`.
overlay_category : list, optional
If only specfic categories within the overlay column is needed, pass their names as a list.
If None, all categories will be used.
markers : list, optional
Markers to be included. If none, all markers will be displayed.
channel_names : list, optional
List of channels in the image in the exact order as image. The default is `adata.uns['all_markers']`
x_coordinate : string, optional
X axis coordinate column name in AnnData object.
y_coordinate : string, optional
Y axis coordinate column name in AnnData object.
point_size : int, optional
point size in the napari plot.
overwrite : bool, optional
In the event you have multiple images in the adata object, ROI can be added to each image
independently using the `imageid` and `subset` parameter. If you wish the results to be
all saved with in the same column set this parameter to `False`. By default, the
function will overwrite the `label` column.
n_jobs : int, optional
Number of cores to use. Default is to use all available cores.
label : string, optional
Key for the returned data, stored in `adata.obs`.
**kwargs
Other arguments that can be passed to napari viewer
Returns:
Modified Anndata object.
Example:
```python
image_path = '/Users/aj/Desktop/ptcl_tma/image.ome.tif'
adata = sm.pl.addROI_image (image_path, adata, label="Tumor Regions")
```
"""
#TODO
# - ADD Subset markers for ZARR ssection
# - Ability to use ZARR metadata if available
# create data matrix that has the co-ordinates
data = pd.DataFrame(adata.obs)[[x_coordinate, y_coordinate, imageid]]
# subset the data if needed
if subset is not None:
# convert string to list
if isinstance(subset, str):
subset = [subset]
# subset data
sub_data = data[data['imageid'].isin(subset)]
else:
sub_data = data
# adding option to load just the image without an adata object
if adata is None:
channel_names = None
else:
# All operations on the AnnData object is performed first
# Plot only the Image that is requested
if subset is not None:
adata_subset = adata[adata.obs[imageid].isin(subset)]
else:
adata_subset = adata.copy()
# Recover the channel names from adata
if channel_names == 'default':
channel_names = adata_subset.uns['all_markers']
else:
channel_names = channel_names
# Index of the marker of interest and corresponding names
if markers is None:
idx = list(range(len(channel_names)))
channel_names = channel_names
else:
idx = []
for i in markers:
idx.append(list(channel_names).index(i))
channel_names = markers
# Load the segmentation mask
if seg_mask is not None:
seg_m = tiff.imread(seg_mask)
if (len(seg_m.shape) > 2) and (seg_m.shape[0] > 1):
seg_m = seg_m[0]
# Operations on the OME TIFF image is performed next
# check the format of image
if os.path.isfile(image_path) is True:
image = tiff.TiffFile(image_path, is_ome=False) #is_ome=False
z = zarr.open(image.aszarr(), mode='r') # convert image to Zarr array
# Identify the number of pyramids
n_levels = len(image.series[0].levels) # pyramid
# If and if not pyramids are available
if n_levels > 1:
pyramid = [da.from_zarr(z[i]) for i in range(n_levels)]
multiscale = True
else:
pyramid = da.from_zarr(z)
multiscale = False
# subset channels of interest
if markers is not None:
if n_levels > 1:
for i in range(n_levels-1):
pyramid[i] = pyramid[i][idx, :, :]
n_channels = pyramid[0].shape[0] # identify the number of channels
else:
pyramid = pyramid[idx, :, :]
n_channels = pyramid.shape[0] # identify the number of channels
else:
if n_levels > 1:
n_channels = pyramid[0].shape[0]
else:
n_channels = pyramid.shape[0]
# check if channel names have been passed to all channels
if channel_names is not None:
assert n_channels == len(channel_names), (
f'number of channel names ({len(channel_names)}) must '
f'match number of channels ({n_channels})'
)
# Load the viewer
viewer = napari.view_image(
pyramid, multiscale=multiscale, channel_axis=0,
visible=False,
name = None if channel_names is None else channel_names,
**kwargs
)
# Operations on the ZARR image
# check the format of image
if os.path.isfile(image_path) is False:
#print(image_path)
viewer = napari.Viewer()
viewer.open(image_path, multiscale=True,
visible=False,
name = None if channel_names is None else channel_names)
# Add the seg mask
if seg_mask is not None:
viewer.add_labels(seg_m, name='segmentation mask', visible=False)
# Add phenotype layer function
def add_phenotype_layer (adata, overlay, phenotype_layer,x,y,viewer,point_size,point_color):
coordinates = adata[adata.obs[overlay] == phenotype_layer]
# Flip Y AXIS if needed
if flip_y is True:
coordinates = pd.DataFrame({'y': coordinates.obs[y],'x': coordinates.obs[x]})
else:
coordinates = pd.DataFrame({'x': coordinates.obs[x],'y': coordinates.obs[y]})
#points = coordinates.values.tolist()
points = coordinates.values
if point_color is None:
r = lambda: random.randint(0,255) # random color generator
point_color = '#%02X%02X%02X' % (r(),r(),r()) # random color generator
viewer.add_points(points, size=point_size,face_color=point_color,visible=False,name=phenotype_layer)
if overlay is not None:
# categories under investigation
if overlay_category is None:
available_phenotypes = list(adata_subset.obs[overlay].unique())
else:
available_phenotypes = overlay_category
# Run the function on all phenotypes
for i in available_phenotypes:
add_phenotype_layer (adata=adata_subset, overlay=overlay,
phenotype_layer=i, x=x_coordinate, y=y_coordinate, viewer=viewer,
point_size=point_size,point_color=point_color)
# Intiate an ROI layer
shape_layer = viewer.add_shapes(name=label)
shape_layer.mode = 'add_polygon'
#_ = show_info('Draw ROIs')
# helper functions
def ellipse_points_to_patch(vertex_1, vertex_2,co_vertex_1, co_vertex_2):
"""
Parameters
----------
vertex_1, vertex_2, co_vertex_1, co_vertex_2: array like, in the form of (x-coordinate, y-coordinate)
"""
v_and_co_v = np.array([
vertex_1, vertex_2,
co_vertex_1, co_vertex_2
])
centers = v_and_co_v.mean(axis=0)
d = sdistance.cdist(v_and_co_v, v_and_co_v, metric='euclidean')
width = d[0, 1]
height = d[2, 3]
vector_2 = v_and_co_v[1] - v_and_co_v[0]
vector_2 /= np.linalg.norm(vector_2)
angle = np.degrees(np.arccos([1, 0] @ vector_2))
ellipse_patch = mpatches.Ellipse(
centers, width=width, height=height, angle=angle
)
return ellipse_patch
# block the viewer until ROI is added
a = """
Opening Napari;
Add shape layers (on left) to draw ROI's.
Rename the shape layer to give a name to your ROI
Multiple shape layers are supported
ROI's should not overlap
Close Napari to save ROI's.
"""
print(a)
viewer.show(block=True)
# Find all the shape layers
my_shapes = [layer for layer in viewer.layers if isinstance(layer, Shapes)]
# loop through the layers to find their names
shape_names = []
added_rois = []
for i in my_shapes:
shape_names.append(i.name)
added_rois.append(len(viewer.layers[i.name].data))
if any(y > 0 for y in added_rois):
# Loop through all the Shape layers and extract the vertices and shape type
all_rois = pd.DataFrame()
for i in shape_names:
# return shape vertices
ver = viewer.layers[i].data
# return shape shape
structure = viewer.layers[i].shape_type
# Each layer may contain multiple ROI's with different shapes (handle that)
napari_roi_table = pd.DataFrame(dict(vertices=[np.fliplr(v) for v in ver], type=[str(s) for s in structure], ROI=i))
# convert gathered ROIs into mpatches
for i in range(len(napari_roi_table.index)):
# ellipse
if napari_roi_table['type'][i] == 'ellipse':
napari_roi_table.loc[i:, 'mpatch'] = ellipse_points_to_patch(napari_roi_table['vertices'][i][0],
napari_roi_table['vertices'][i][1],
napari_roi_table['vertices'][i][2],
napari_roi_table['vertices'][i][3])
# polygon, rectangle, line
elif napari_roi_table['type'][i] in ['rectangle', 'polygon', 'path']:
napari_roi_table.loc[i:, 'mpatch'] = mpatches.Polygon(napari_roi_table['vertices'][i], closed=True)
else:
raise ValueError
# merge the final ROI's across all shape layers
all_rois = pd.concat ([all_rois, napari_roi_table])
all_rois['id'] = range(all_rois.shape[0])
all_rois.index = all_rois['id']
# Find cells within ROI and add it to the dataframe
def add_roi_internal (roi_id):
roi_mpatch = all_rois[all_rois['id'] == roi_id]['mpatch'][roi_id]
inside = sub_data[roi_mpatch.contains_points(sub_data[[x_coordinate, y_coordinate]])]
inside['ROI_internal'] = all_rois[all_rois['id'] == roi_id]['ROI'][roi_id]
# return
return inside
print("Identifying cells within selected ROI's")
# all ROI cells
roi_list = all_rois['id'].unique()
#final_roi = list()
#for i in roi_list:
# roi_mpatch = all_rois[all_rois['id'] == i]['mpatch'][i]
# inside = sub_data[roi_mpatch.contains_points(sub_data[[x_coordinate, y_coordinate]])]
# inside['ROI_internal'] = all_rois[all_rois['id'] == i]['ROI'][i]
# final_roi.append(inside)
final_roi = Parallel(n_jobs=n_jobs, verbose=verbose)(delayed(add_roi_internal)(roi_id=i) for i in roi_list)
# Merge all into a single DF
final_roi = pd.concat(final_roi)[['ROI_internal']]
# Add the list to obs
result = pd.merge(data, final_roi, left_index=True, right_index=True, how='outer')
# Reindex
result = result.reindex(adata.obs.index)
# check if adata already has a column with the supplied label
# if avaialble overwrite or append depending on users choice
if label in adata.obs.columns:
if overwrite is False:
# Append
# retreive the ROI information
old_roi = adata.obs[label]
combined_roi = pd.merge(result, old_roi, left_index=True, right_index=True, how='outer')
combined_roi['ROI_internal'] = combined_roi['ROI_internal'].fillna(combined_roi[label])
else:
# Over write
combined_roi = result.copy()
combined_roi['ROI_internal'] = combined_roi['ROI_internal'].fillna('Other')
else:
# if label is not present just create a new one
combined_roi = result.copy()
combined_roi['ROI_internal'] = combined_roi['ROI_internal'].fillna('Other')
# Add to adata
adata.obs[label] = combined_roi['ROI_internal']
# return
return adata
|