Tip
Download this tutorial as a Jupyter notebook, or as a python script with code cells. We highly recommend using Visual Studio Code to execute this tutorial. Alternatively, you could run the Python script in a terminal with python analysis_single_illumination.py from the folder where the file is located.
Full analysis APD maps: single illumination (RED)¶
The goal on this jupyter notebook is to try to create APD maps from videos of panaromaic heart that were recorded using single ilumination (RED light) in a automatic or semi-automatic way, that is programatically.
this notebook is separated by sections that are summarized as follow:
First we will load the libraries required for the analysis.
Then we will base on the experiment conditions we create folders where the results will be save it at the end of the analysis.
we will then try to launch the napari viewer and attached the napari-omaas plugin to open a file.
for this example we will use only the anterior view, so we will crop & rotate the image as needed.
We will then try to visualize the plot profile of the fluorescence signal of the this example image along the whole time series using a predefined ROI.
We will clip the image to only 10 beats.
We will then on the resulting image apply motion correction
then we will preprocess the image by invert and normalization methods an thereafter we will apply spatial and temporal filteres to improve the image quality and finally segment the shape of the heart.
with this preporcessed image we will proceed on averaging temporally the 10 beats to a single one.
after this step we are ready to compute APD from a single beat image
post-processing of the APD could be done manually and after getting nicely visualy APD maps we can then save the resulting APD maps to the results folder initially created.
Load libraries¶
%load_ext autoreload
%autoreload 2
import napari
import napari_omaas
from napari_omaas import utils
from napari.utils import nbscreenshot
import numpy as np
from napari_omaas.custom_exceptions import CustomException
import sys
import os
from pathlib import Path
import tempfile
Define helper functions¶
def show_last_layer(viewer):
"""
Helper function to hide all layers but last one.
Parameters
----------
viewer : a viewer instance
"""
for index, layer in enumerate(viewer.layers):
if index!=(len(viewer.layers)-1):
layer.visible = False
# viewer = napari.Viewer(show=False)
viewer = napari.Viewer(show=True)
o = napari_omaas.OMAAS(viewer)
viewer.window.add_dock_widget(o, area='right')
##################
## SupRep-SQT1 ###
##################
# my_file = r"\\PC160\GroupOdening\OM_data_Bern\raw_data\2024\20241110\Videosdata\20241110_12h-06m-25" #control 4Hz 11.10.24 #7026 SupRep-SQT1 *red only* *done*
# my_file = r"\\PC160\GroupOdening\OM_data_Bern\raw_data\2024\20241110\Videosdata\20241110_12h-20m-25" #carbachol 2.5Hz 11.10.24 #7026 SupRep-SQT1 *red only* *done*
NOTE:¶
For this example, we will not be working directly with files from our institutional storage server (GroupOdenning). Such files are, for example, located at the path addreses listed above for experiments conducted by PhD Saranda Nimani.
Instead, we will be using one of the files listed above labeled with the following timestanmp 20241110_12h-20m-25, which is accesibale online the website of the Insititute of Physiology as single_illumination_spool_data_sample (772 MB).
The resulting ouput from this analysis pipeline will be sotrage in a temporal folder, as defined further bellow. You can change this folder to your prefered folder.
my_file = "20241110_12h-20m-25"
Define the experiment conditions¶
The conditions here defined are based on specific experiment. Details on that can be found in the lab sotrage driver wher this folder resides together with scan of the labbook with its respective annotations.
In this especific example we were using two different stimulation condition (2.5 and 4 Hz). We use the following drugs conditions: (carbachol and control). We used four different genotypes conditions: (WT, Sham-SQT1, UT-SQT1 and SupRep-SQT1) and finally we use the animal ID.
freq_estim_condition = "2.5Hz"
# freq_estim_condition = "4Hz"
condition = "carbachol"
# condition = "control"
# illumination_type = "dual_illumination"
illumination_type = "red_only"
# genotype = "WT"
# genotype = "Sham-SQT1"
# genotype = "UT-SQT1"
genotype = "SupRep-SQT1"
# animal_id = str(7028)
# animal_id = str(7004)
# animal_id = str(7009)
animal_id = str(7026)
Create result folders¶
For the pourpose of this tutorial, in this example we are using a temporal folder.
You may want to change this folder to your prefered folder to store your results.
############# here we create the saving results folder #############
# results_folder_name = fr"M:\PhD students\Saranda Nimani\Optical Mapping\SQT-SupRep\APD maps\{genotype}\{freq_estim_condition}\{condition}\{illumination_type}\{animal_id}\{os.path.basename(my_file)}"
# results_folder_name = fr"APD maps\{genotype}\{freq_estim_condition}\{condition}\{illumination_type}\{animal_id}\{os.path.basename(my_file)}"
results_folder_path = Path(tempfile.gettempdir()) / "napari_omaas_sample_data" / "APD_maps_tutorial_results" / genotype / freq_estim_condition / condition / illumination_type / animal_id / os.path.basename(my_file)
# results_folder_path.mkdir(exist_ok=True)
# results_folder_path = os.path.normpath(results_folder_name)
if not results_folder_path.exists():
print(f"Creating Folder: \n'{results_folder_path}'\n*Done*.")
results_folder_path.mkdir(parents=True, exist_ok=True)
else:
print(f"Folder: \n'{results_folder_path}'\nAlready exists")
############# here we create a folder to save the mask for segmentation #############
mask_folder = results_folder_path / "mask"
if not os.path.exists(mask_folder):
print(f"Creating Folder: \n'{mask_folder}'\n*Done*.")
mask_folder.mkdir(parents=True, exist_ok=True)
else:
print(f"Folder: \n'{mask_folder}'\nAlready exists")
Open File¶
# %%capture
# the %%capture command above will hide the ouptut of this cell bc can be very long (only used for the documentation)
try:
# viewer.open(path=my_file, plugin= "napari-omaas")
viewer.open_sample('napari-omaas', 'heart_sample_single_illumination')
except Exception as e:
raise CustomException(e, sys)
Let’s now see what we get in the viewer
nbscreenshot(viewer)
some info from this file¶
Let’s explore the content of the metadata of the recently downloaded file as indicated bellow.
viewer.layers[-1].metadata
Crop & rotate shape¶
In this section we create a rectangular shape that we will use to crop the anterior view of the (most-left) of the panaromic image.
my_shape = [np.array([[-27.10448098, 376.83908381],
[-27.10448098, 67.77374194],
[272.08514108, 67.77374194],
[272.08514108, 376.83908381]])]
viewer.add_shapes(my_shape)
o.rotate_l_crop.setChecked(True)
o.crop_from_shape_btn.click() # done
show_last_layer(viewer=viewer)
# this is just to manipulate the viewer by changing the position of the image, zooming it and changing to the first frame
viewer.camera.center = 0.0, 150, 127
viewer.camera.zoom = 2.3
viewer.dims.current_step = (0, 156, 749)
Alternatively, you can use the function to rotate automatically rotate the views in correct position amd work with the four views as shown bellow.
# # change to tab called "Shapes"
# o.tabs.setCurrentIndex(1)
# o.crop_all_views_and_rotate_btn.click()
# show_last_layer(viewer=viewer)
plot profile¶
We create a new shape that we will use to plot the data along the time axis in the image stack and visualize this way the fluorescence signal.
my_shape =[np.array([[128.03402574, 142.83315215],
[128.03402574, 187.82476519],
[178.56426079, 187.82476519],
[178.56426079, 142.83315215]])]
viewer.add_shapes(my_shape, name="ROI_1")
o.plot_last_generated_img(shape_indx=1)
nbscreenshot(viewer)
What we observe in the screenshot above is the plot profile of the croped image (Anterior view) using the ROI from the shape layer ROI_1.
We also observe an artefact at approximateley time ~ 4800 ms. In the next section we will clip the file to only 10 beat cycles.
Clip trace to 10 APs¶
# you can call the current clipping values with the following command:
o.double_slider_clip_trace.value()
clipping_values = (569.5483417085427, 4576.026331658292) #carbachol 2.5Hz 11.10.24 #7026 SupRep-SQT1 *red only*
o.plot_last_generated_img(shape_indx=1)
o.is_range_clicked_checkbox.setChecked(True)
o.double_slider_clip_trace.setValue(clipping_values)
o.clip_trace_btn.click()
we can visualize again the fluorescence intensity of the resulting image and check that only contain 10 cycles.
First le’s hide the layers we dont use anymore and keep visible only the last one from now and center the image in the viewer.
show_last_layer(viewer)
nbscreenshot(viewer)
Apply motion correction¶
We will apply motion correction to the resulting image, using borowwed method from the package optimap from Jan Lebber and Jan Chritoph.
You can explore more on this methods at their library documentation
For this example, we add known setting that work ok for most of the cases.
%%capture
# the %%capture command above will hide the ouptut of this cell bc can be very long (only used for the documentation)
curr_imag = viewer.layers[-1]
viewer.layers.selection.active = curr_imag
contrast_kernel = 3
ref_frame = 10
o.pre_smooth_temp.setValue(1)
o.pre_smooth_spat.setValue(1)
o.c_kernels.setValue(contrast_kernel)
o.ref_frame_val.setText(str(ref_frame))
o.apply_optimap_mot_corr_btn.click()
You can explore interactively the resulting image by using the viewer and compare with the original to explore anc check the quality of the motion correction.
Invert and normalize¶
# Using slide window method for normalization to detrend the signal
o.data_normalization_options.setCurrentText("Slide window")
o.slide_wind_n.setValue(200)
# o.data_normalization_options.setCurrentText("Local max") # Note: Use when stable traces or when "Slide window" is creating artifacts
o.inv_and_norm_data_btn.click()
o.plot_last_generated_img(shape_indx=1)
nbscreenshot(viewer)
Apply spatial and temporal filters¶
o.spat_filter_types.setCurrentText("Median")
o.apply_spat_filt_btn.click()
o.apply_temp_filt_btn.click()
o.plot_last_generated_img(shape_indx=1)
nbscreenshot(viewer)
Segment Image¶
show_last_layer(viewer)
mask_file_path = os.path.join(mask_folder, "Heart_labels_NullBckgrnd.tif")
# use another image layer to create mask, then use that mask to to apply segmentation to last layer
viewer.layers.selection.active = viewer.layers[5] # take the ratio image normalized
o.return_img_no_backg_btn.setChecked(False)
# o.return_img_no_backg_btn.setChecked(True)
o.apply_auto_segmentation_btn.click()
viewer.layers.select_previous()
o.apply_manual_segmentation_btn.click() # done
show_last_layer(viewer)
# save Mask if newly created
try:
print(f"Saving mask to: \n'{mask_file_path}'\n*Done*.")
viewer.layers[-1].save(mask_file_path)
except Exception as e:
CustomException(e, sys)
If you want to reuse the same mask for another similar image, uncoment the next code cell
# curr_imag = viewer.layers[-1]
# try:
# viewer.open(mask_file_path)
# except Exception as e:
# CustomException(e, sys)
# viewer.layers.selection.active = curr_imag
# o.apply_manual_segmentation_btn.click() # done
# dont forget to change mask folder!
Average APs¶
We create a new ROI from wich we will use their averaged values profile to average the 10 beat (cardiac cycles) in our trace.
# del viewer.layers[0]
my_shape = [np.array([[136.79882203, 101.27454257],
[136.79882203, 168.68341024],
[205.2137922 , 168.68341024],
[205.2137922 , 101.27454257]])]
viewer.add_shapes(my_shape)
o.plot_last_generated_img(shape_indx=2)
# preview AP splitting results
o.tabs.setCurrentIndex(3)
o.preview_AP_splitted_btn.click()
nbscreenshot(viewer)
# create average from splitted APs
o.create_average_AP_btn.click()
o.plot_last_generated_img(shape_indx=2)
o.preview_AP_splitted_btn.click()
nbscreenshot(viewer)
# Save averaged image
# o.save_img_dir_box_text.setText(results_folder_path)
# # for value in [-1]:
# for value in [-1, -2]:
# viewer.layers.selection.active = viewer.layers[value]
# o.export_image_btn.click()
Compute APD maps¶
Let’s move the current viewer to the first image frame.
viewer.dims.current_step = (0, 156, 749)
o.tabs.setCurrentIndex(3)
thresh_value = o.slider_APD_detection_threshold.value() * 0.0001
thresh_value
# You could also adjust the APD detection trheshold as shown bellow.
# thresh_value = 0.0146 #control 4Hz 30.07.24
o.slider_APD_detection_threshold.setValue(int(thresh_value * 10000))
%%capture
# the %%capture command above will hide the ouptut of this cell bc can be very long (only used for the documentation)
target_image = viewer.layers[-1]
# if you wish to compute multiples APD maps with different % values, replace
# the value inside the square brakets with your desired values
# like for example: [25, 75, 90]
for value in [90]:
viewer.layers.selection.active = target_image
o.slider_APD_map_percentage.setValue(value)
o.toggle_map_type.setChecked(True)
o.make_maps_btn.click()
Visualize the resulting map (APD-90)¶
show_last_layer(viewer)
nbscreenshot(viewer)
In the Post-processing Maps tab, the resulting map can be better visualized as a contour plot. For that, we will change the current tab to the Post-processing Maps tab and we select the map generated to plot it.
Let’s have a look:
o.tabs.setCurrentIndex(3)
o.mapping_tabs.setCurrentIndex(1)
# select the first item in the list of images maps (only one present at the moment)
o.map_imgs_selector.item(0).setSelected(True)
o.plot_curr_map_btn.click()
# o.clear_curr_map_btn.click()
nbscreenshot(viewer)
Post-process the resulting map¶
You can retouch and smooth the resulting APD map using the same Post-processing map tab for visualization.
The bellow commands would be the equivalent to manually click on preview and adjust process the image via eroding or filtering using a gaussian filer. Let’s see how:
o.tabs.setCurrentIndex(3)
o.mapping_tabs.setCurrentIndex(1)
o.preview_postProcessingMAP_btn.click()
# set the sigam value for the gaussian filter
o.InterctiveWindod_edit_map.gaussian_sigma.setValue(1.5)
# set the radius of the filter
o.InterctiveWindod_edit_map.gaussian_radius.setValue(4)
# set the numebr of pixels to erode (reduce the edge) of the image
o.InterctiveWindod_edit_map.n_pixels_erode_slider.setValue(4)
If you are happy wiht the results, you can accept the changes and close the postproceeesing window
# Accept changes
o.InterctiveWindod_edit_map.accept_post_processing_changes_btn.click()
# Close the tool
o.InterctiveWindod_edit_map.close_postprocessing_map_window_btn.click()
Let’s now plot the two MAPs images side by side to check the differences visually.
# select the first item in the list of images maps (only one present at the moment)
o.map_imgs_selector.item(0).setSelected(True)
o.map_imgs_selector.item(1).setSelected(True)
o.plot_curr_map_btn.click()
nbscreenshot(viewer=viewer)
Save maps¶
We save the map generated (APD-90) and the last image (average and preprocessed stack)
# Here we export last 3 images
o.save_img_dir_box_text.setText(str(results_folder_path))
# for value in [-1]:
for value in range(-1, -4, -1):
viewer.layers.selection.active = viewer.layers[value]
o.export_image_btn.click()