Script error: peak_local_max() function of most recent skimage version changed

Hannah Scharff
Hannah Scharff
May 20, 2024
protocol Protocol: Automated Layer Analysis (ALAn): An Image Analysis Tool for the Unbiased Characterization of Mammalian Epithelial Architecture in Culture

Hello,


I am trying to use your code starting with the available data (http://tinyurl.com/4kv7jfav.) and it is not working properly. The problem is because of this line (in the xy_segment Read more Read less more less

answer Answer
recommend recommend recommend recommend Recommend
follow follow follow follow Follow
share share Share
1 answer
Sort by: most helpful / most recent

Christian Cammarota Author Answered May 22, 2024

Rochester Institute of Technology Rochester

Hi Hannah,


Thank you for reaching out! I just checked and you're definitely right. It looks like in the recent versions of skimage, the output of peak_local_max was forced to be a list of indices. For the XY segmentation function, we use a projection of the cell layer onto the XY plan and run the segmentation through by creating a distance map of the binarized actin signals, and identifying unique cells with peak_local_max.


I quickly wrote up a 'band-aid' for the functionality that converts skimage's current output to the pervious output we used with indices=False. I'm currently working on my postdoc and can't really test this as thoroughly as I'd like to (and as such don't think it's time to update our github just yet before it's fixed), but I'll be in contact with the Finegan-Bergstralh Lab to make sure this is updated well. I hope the 'band-aid' works for now, and thanks for bringing this to our attention! I'm happy to send a Jupyter notebook of the 'band-aid' to you if you'd like it, otherwise the instructions are listed below.


To fix, you can replace the 'markers = measure.label(local_peaks)' with the following (making sure to indent the aux[loacl_peaks[i][o]....... line to be within that for loop):


  aux = np.zeros_like(distance_map)

  for i in range(len(local_peaks)):

    aux[local_peaks[i][0], local_peaks[i][1]] = True

  markers = measure.label(aux)



Alternatively, replace the entire xy_segmentation function with the following (again making sure that all of the appropriate indents are in place):

def xy_segmentation(df, image, **kwargs): 

  image = image_shuffle(image)

  x_max, y_max, z_max, x_min, y_min, z_min, slice_heights = get_layer_position(df, image)

  scale = (x_max - x_min)/image.shape[2]

  df_cleared = clear_debris(df)

  density = len(df_cleared)/(x_max-x_min)**2*1000

  plot = kwargs.get('plot', False)

  save_name = kwargs.get('save_name', False)

  actin_channel = kwargs.get('actin_channel', 1)

  dapi_channel = 0

  invert = kwargs.get('invert', False)

   

  if invert != False:

    xy_proj_actin = np.sum(np.sum(image, axis = 3), axis = 2)[:,actin_channel].copy()[::-1] 

  else:   

    xy_proj_actin = np.sum(np.sum(image, axis = 3), axis = 2)[:,actin_channel].copy()

   

  norm_intensities_actin = (xy_proj_actin-np.min(xy_proj_actin))/np.max(xy_proj_actin-np.min(xy_proj_actin))

  df_cleared = clear_debris(df)

  density = len(df_cleared)/(x_max-x_min)**2*1000

  if density <=2:

    bot_cutoff = 0.5

    top_cutoff = 0.6

  elif density>=6:

    bot_cutoff = 0.2

    top_cutoff = 0.8

  else:

    bot_cutoff = 0.5-0.3*(density-2)/4

    top_cutoff = 0.6+0.2*(density-2)/4

   

  start = np.argwhere(norm_intensities_actin >= bot_cutoff)[0][0]

  end = np.argwhere(norm_intensities_actin >= top_cutoff)[-1][0]

   

  truncated_stack = image[start+3:end-3, :, :, :].copy()

   

  sum_project = np.sum(truncated_stack, axis = 0)

  DNA = sum_project[dapi_channel,:,:]

  Actin = sum_project[actin_channel,:,:]

   

  filtered = filters.gaussian(Actin)


  thresh = filters.threshold_local(filtered, block_size = 55)

  binary = filtered > thresh

   

  clean_image = morphology.remove_small_objects(morphology.remove_small_objects(binary))

  dilate = morphology.dilation(morphology.dilation(morphology.dilation(clean_image)))

  skeleton = morphology.skeletonize(dilate)

  dilate_skeleton = morphology.dilation(morphology.dilation(morphology.dilation(morphology.dilation(skeleton))))

  erosion_skeleton = morphology.erosion(morphology.erosion(morphology.erosion(dilate_skeleton)))

  cell_mask = np.invert(erosion_skeleton)

  clean_mask = clear_border(cell_mask)

  distance_map = ndimage.distance_transform_edt(cell_mask)

  local_peaks = peak_local_max(distance_map, min_distance = 10, threshold_abs = 2)

  aux = np.zeros_like(distance_map)

  for i in range(len(local_peaks)):

    aux[local_peaks[i][0], local_peaks[i][1]] = True

  markers = measure.label(aux)

  watershed_map = watershed(-distance_map, markers)

  labeled_image = watershed_map * clean_mask


  cell_props = regionprops(labeled_image)


  cell_ids = []

  cell_areas = []

  cell_centroid_rows = []

  cell_centroid_cols = []

  cell_perimeters = []

  cell_eccentricities = []

  cell_circularities = []

     


  for cell in cell_props:

    cell_ids.append(cell.label)

    cell_areas.append(cell.area * scale**2)

    row, col = cell.centroid

    cell_centroid_rows.append(row)

    cell_centroid_cols.append(col)

    cell_perimeters.append(cell.perimeter * scale)

    cell_eccentricities.append(cell.eccentricity)

    cell_circularities.append(4. * np.pi * cell.area / cell.perimeter**2)


   

  if save_name != False:

     

    data_dict = {'labels': cell_ids,

           'areas': cell_areas,

           'crows': cell_centroid_rows,

           'ccols': cell_centroid_cols,

           'perimeters': cell_perimeters,

           'eccentricities': cell_eccentricities,

           'circularities': cell_circularities}


    df = pd.DataFrame(data_dict)

    df.to_csv(save_name + '.csv')


  return np.mean(cell_areas), np.mean(cell_perimeters), np.mean(cell_circularities)

0 helpful
0 unhelpful
comment 0 comments down up
share Share

My answer

Write your answer...

References (optional)

add Add more

post Post an Answer