Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add textbook-like tutorial on measuring fluorescence at nuclear envelope. #5262

Merged
merged 36 commits into from
Jan 26, 2022

Conversation

mkcor
Copy link
Member

@mkcor mkcor commented Mar 9, 2021

Description

This draft PR fits the overall goal of offering more gallery examples for the life sciences (see #4601).

Here I'm following this amazing suggestion by @haesleinhuepf: https://forum.image.sc/t/looking-for-life-scientists-to-collaborate-on-scikit-image-tutorials/49073/15

Checklist

For reviewers

  • Check that the PR title is short, concise, and will make sense 1 year
    later.
  • Check that new functions are imported in corresponding __init__.py.
  • Check that new features, API changes, and deprecations are mentioned in
    doc/release/release_dev.rst.

@jni
Copy link
Member

jni commented Mar 11, 2021

👋 @mkcor this is very exciting. =) I don't think for this tutorial it's needed because the nuclei are not close together (?), but I just want to give a shout-out to @VolkerH's segmentation.expand_labels which is like dilation but smarter. =)

@grlee77 grlee77 added the 📄 type: Documentation Updates, fixes and additions to documentation label Mar 11, 2021
@mkcor
Copy link
Member Author

mkcor commented Jul 18, 2021

@mkcor this is very exciting. =) I don't think for this tutorial it's needed because the nuclei are not close together (?), but I just want to give a shout-out to @VolkerH's segmentation.expand_labels which is like dilation but smarter. =)

Thank you for the great suggestion, @jni! I decided to go for @VolkerH's segmentation.expand_labels even though it's not necessary with this dataset. Actually (reassuringly, maybe), these two results are equal:

expand = segmentation.expand_labels(label)
dilate = morphology.dilation(label)

(where label is the labelled image after blurring, thresholding, filling in holes, and keeping only the one nucleus of interest) and yield this difference, when subtracting morphology.erosion(label):
difference_1px

The rim is then about 2px-thick (I zoomed in with Plotly to check). But I tried reproducing the original processing faithfully, and it's straightforward with segmentation.expand_labels's keyword argument:

expand = segmentation.expand_labels(label, distance=4)

where I can just expand more, as much as I like! 😄
The resulting rim is about 5px-thick.

@pep8speaks
Copy link

pep8speaks commented Sep 21, 2021

Hello @mkcor! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2022-01-26 13:31:56 UTC

@mkcor mkcor marked this pull request as ready for review September 21, 2021 17:34
Copy link
Contributor

@haesleinhuepf haesleinhuepf left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @mkcor ,

super cool to see this evolving, thanks for working on this! As I translated the workflow at some point as well, I was curious and took a look at your implementation. I found some differences compared to the original implementation. I don't want to be picky; I was just wondering if you made those changes intentionally and if so, why. :-) Also did you compare quantitative measurements with the measurements generated with your script? You find quantitative measurements from the original workflow in ImageJ for example here and from an alternative implementation here.

Also one more minor comment: I personally prefer variable names such as "input_image" over "dat". Also "channel_0_image" could be more self-explanatory than "ch0t0". But maybe I'm just too used to writing longVariableNamesInJava ;-)

Great work otherwise! Many biologists will appreciate this example in scikit-image and the detailed explanation of the used commands.

Thanks again!

Best,
Robert


smooth = filters.gaussian(ch0t0, sigma=2)

thresh = smooth > 0.1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the original workflow Otsu's threshold method was used. Is there a reason why we're using a fixed constant threshold here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, it's much better to infer the threshold:
otsu

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed in 02c2823.

"Reassuringly, we find the expected result: The total fluorescence intensity at the nuclear envelope increases 1.3-fold in the initial five time points, and then becomes roughly constant."

Now that I'm more consistent with the original workflow, the overall result is preserved but, comparing with the textbook plot, we can see sharper dips (below 1.2) for time points 7 and 12:
original_wkfl


thresh = smooth > 0.1

fill = ndi.binary_fill_holes(thresh)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also this step is not in the original workflow. I could imagine that it's not necessary if one uses Otsu's thresholding method...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, well, this was always a bit of a mystery to me... Looking at steps c) and d), I've wondered how the small holes were handled: Does dilation automatically take care of filling them?

I don't see what difference using Otsu's thresholding method would make (see above #5262 (comment))... Even with multi-Otsu thresholding (defaulting to 3 classes), and keeping the lower threshold, there are still small holes inside the main region:
multi_otsu


#####################################################################
# Although it is not strictly necessary, let us discard the other nucleus part
# visible in the bottom right-hand corner.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In ImageJ macro, we're using the function of the Particle Analyzer for removing objects that touch the image border. How about using clear_border? I could imagine that the script becomes then more applicable also to other images and use-cases.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great point! Your suggestion is much more readable and much more elegant than what I had; done 3c99e28.

label_small = regions['label'][np.argmin(regions['area'])]
label[label == label_small] = 0 # background value

expand = segmentation.expand_labels(label, distance=4)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the original workflow, we erode the mask and dilate it, which means in ImageJ to take one pixel as radius.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, fixed in 02c2823. Does the 5ish-px-looking thickness in the final processing step only come from the rendering/printing then?

# Segment the nucleus rim
# =======================

smooth = filters.gaussian(ch0t0, sigma=2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for pointing me to the ImageJ Macro source! I admit I didn't look it up; I followed the book chapter mostly qualitatively, but let's try and get as close as possible quantitatively...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed in 02c2823.

assert props['label'] == 1
intensity_total = props['area'] * props['intensity_mean']
fluorescence_change.append(intensity_total)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could mention here that the following part is not in the original workflow. But as we are in python, we can do plotting and image analysis together in one script. :-)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, the conclusion reads "Here, to summarize the analysis in this chapter, we plot the changes in the total fluorescence intensity over time using ImageJ Macro code_plotResults.ijm (◘ Fig. 2.8). The code appears after the paragraph below."

So, I get that it's not part of the 'original workflow' but that's not because ImageJ would be limiting; is it?

Alternatively, I could change my introduction to make it explicit that I'm following this book chapter in particular, which itself is based on the original workflow.

@mkcor
Copy link
Member Author

mkcor commented Sep 21, 2021

@haesleinhuepf @jni please review: https://2809-2014929-gh.circle-artifacts.com/0/doc/build/html/auto_examples/applications/plot_fluorescence_nuclear_envelope.html

The CI failures come from

dat = imageio.volread('http://cmci.embl.de/sampleimages/NPCsingleNucleus.tif')

where urllib.error.URLError: <urlopen error timed out> -- I'll try and reach out to NeuBIAS to see if we could include this dataset in scikit-image's data registry.

@haesleinhuepf
Copy link
Contributor

haesleinhuepf commented Sep 21, 2021

Hey @mkcor ,

sorry, I was too fast with the review :-D

I'll try and reach out to NeuBIAS to see if we could include this dataset in scikit-image's data registry.

The image was made by Andrea Boni (I can give you his email, if you tell me yours, or via zulip or twitter or something) . You can ask him for permission and then also download it here if you like.

@mkcor
Copy link
Member Author

mkcor commented Sep 21, 2021

Hi @haesleinhuepf,

Wow! You've been lightning fast and our comments have crossed! Thank you so much for your review. Please be picky! 😉

Also one more minor comment: I personally prefer variable names such as "input_image" over "dat".

I hear you. In this case, would you agree with either image_stack or image_sequence?

@haesleinhuepf
Copy link
Contributor

image_sequence

This sounds great! :-)

@mkcor
Copy link
Member Author

mkcor commented Sep 21, 2021

The image was made by Andrea Boni (I can give you his email, if you tell me yours, or via zulip or twitter or something) . You can ask him for permission and then also download it here if you like.

Thank you, @haesleinhuepf! I can see there is data provenance info here, but there is no explicit licensing info... We prefer using datasets which are in the public domain or under CC0: https://gitlab.com/scikit-image/data/#data

Can you read my email address on my GitHub profile?

@mkcor
Copy link
Member Author

mkcor commented Sep 21, 2021

Also "channel_0_image" could be more self-explanatory than "ch0t0".

@haesleinhuepf how about t_0_channel_0_img?

@haesleinhuepf
Copy link
Contributor

@haesleinhuepf how about t_0_channel_0_img?

Sure. Sounds good :-)

@mkcor
Copy link
Member Author

mkcor commented Sep 21, 2021

@haesleinhuepf how about t_0_channel_0_img?

Sure. Sounds good :-)

Or maybe image_t_0_channel_0? 🤔 We use names such as 'intensity_mean' after all... (less English but more hierarchical)

@mkcor
Copy link
Member Author

mkcor commented Dec 31, 2021

This is ready for final review, thanks!

/cc @hmaarrfk 😉

@mkcor mkcor added this to the 0.19.2 milestone Jan 13, 2022
@mkcor
Copy link
Member Author

mkcor commented Jan 13, 2022

@grlee77 I've added the 0.19.2 milestone so this new tutorial ships along with #6087 but, actually, I'm not exactly sure this kind of addition belongs in a patch release... 🤔

@grlee77
Copy link
Contributor

grlee77 commented Jan 24, 2022

@grlee77 I've added the 0.19.2 milestone so this new tutorial ships along with #6087 but, actually, I'm not exactly sure this kind of addition belongs in a patch release... thinking

Okay, I will review this soon. I think it should be fine to add in 0.19.2 since it is just supplementing docs and not changing behavior of existing functions.

Copy link
Member

@rfezzani rfezzani left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @mkcor, I added a small comment concerning the step labels in the first figure's titles, suggested to tighten figure's layout and proposed few minor code style modifications.

@mkcor mkcor changed the title Start tutorial for measuring fluorescence at nuclear envelope. Add textbook-like tutorial on measuring fluorescence at nuclear envelope. Jan 24, 2022
Copy link
Contributor

@grlee77 grlee77 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks pretty nice @mkcor. I suggested some minor efficiency improvements that hopefully don't detract from readability. I will go ahead and approve as they are suggestions and don't change the final result.

mkcor and others added 2 commits January 26, 2022 14:13
Co-authored-by: Gregory Lee <grlee77@gmail.com>
@mkcor
Copy link
Member Author

mkcor commented Jan 26, 2022

The CI failure is unrelated to this PR. It comes from another gallery example, namely the one on image stitching:

Warning, treated as error:
/home/runner/work/scikit-image/scikit-image/doc/examples/registration/plot_stitching.py failed to execute correctly: Traceback (most recent call last):
[...]
  File "/home/runner/work/scikit-image/scikit-image/doc/examples/registration/plot_stitching.py", line 110, in <listcomp>
    trfm_list = [measure.ransac((dst, src),
AttributeError: 'NoneType' object has no attribute 'params'

where the full code is:

trfm_list = [measure.ransac((dst, src),
                            transform.EuclideanTransform, min_samples=3,
                            residual_threshold=2, max_trials=100)[0].params
             for dst in matching_corners]

But it's weird that the list of measure.ransac calls would return None just for Python 3.9 (3.8 and 3.10 being fine).

/cc @rfezzani

@grlee77
Copy link
Contributor

grlee77 commented Jan 26, 2022

But it's weird that the list of measure.ransac calls would return None just for Python 3.9 (3.8 and 3.10 being fine).

Let me repeat the failing case and see if it resolves itself. I'm pretty sure I've seen that same failure happen at least once before. Maybe we need to fix a random seed or something? (not in this PR)

Copy link
Contributor

@grlee77 grlee77 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM, thanks!

@grlee77 grlee77 merged commit f3ef598 into scikit-image:main Jan 26, 2022
@grlee77
Copy link
Contributor

grlee77 commented Jan 26, 2022

@meeseeksdev backport to v0.19.x

meeseeksmachine pushed a commit to meeseeksmachine/scikit-image that referenced this pull request Jan 26, 2022
grlee77 added a commit that referenced this pull request Jan 26, 2022
…2-on-v0.19.x

Backport PR #5262 on branch v0.19.x (Add textbook-like tutorial on measuring fluorescence at nuclear envelope.)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
📄 type: Documentation Updates, fixes and additions to documentation
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants