Tutorial¶
Here we provide some basic usage examples for the SIMA Python package. These examples can be run in a standard Python shell, or with IPython. Users new to Python may wish to consult the Python documentation.
For more details on the classes and methods that comprise the SIMA package, please consult the SIMA API.
Contents
Importing SIMA¶
Like all Python packages, the SIMA package must be imported prior to use. Here we show a simple example of importing the SIMA package, which results in printing the docstring containing view basic information about the package.
>>> import sima
In all future examples, we assume that the SIMA package has been imported as shown above.
Submodules of the SIMA package also need to be imported before use. For example, the motion correction module can be imported as follows.
>>> import sima.motion
Individual classes or functions can be imported from submodules. For example, we can import the iterable object for use with multi-page TIFF files with the following command:
>>> from sima.motion import HiddenMarkov2D
For more details on importing, consult the Python documentation.
Example data¶
The SIMA package comes with a small amount of example data that, although insufficient to allow for reasonable results from the motion correction and segmentation algorithms, can at least be used to run the functions. For the purposes of this tutorial, we copy the example data into the current working directory.
>>> from shutil import copy, copytree
>>> import sima.misc
>>> copytree(sima.misc.example_data(), 'example.sima')
>>> copy(sima.misc.example_tiff(), 'example.tif')
>>> copy(sima.misc.example_tiff(), 'example_Ch1.tif')
>>> copy(sima.misc.example_tiff(), 'example_Ch2.tif')
>>> copy(sima.misc.example_hdf5(), 'example.h5')
>>> copy(sima.misc.example_imagej_rois(), 'imageJ_ROIs.zip')
Creating an ImagingDataset object¶
The SIMA package is centers around the ImagingDataset object class. A single ImagingDataset object can contain imaging data from multiple sequences (i.e. continuous imaging epochs/trials) acquired at the same imaging location during the same imaging session. The subsections below provide examples of how to initialize ImagingDataset objects using raw data in a variety of formats, including Numpy arrays, TIFF files, and HDF5 files.
The ImagingDataset object is permanently stored in the location (ending with extension .sima) specified during initialization. Results of segmentation, signal extraction, and other alterations to the ImagingDataset object are automatically saved to this location.
Numpy arrays¶
To begin with, we create some Numpy arrays containing random data. The shape of these arrays is (num_frames, num_planes, num_rows, num_columns, num_channels).
>>> import numpy as np
>>> cycle1_channel1 = np.random.rand(100, 1, 128, 128, 1)
>>> cycle1_channel2 = np.random.rand(100, 1, 128, 128, 1)
>>> cycle2_channel1 = np.random.rand(100, 1, 128, 128, 1)
>>> cycle2_channel2 = np.random.rand(100, 1, 128, 128, 1)
Once we have the Numpy arrays containing the imaging data, we create the ImagingDataset object as follows.
>>> sequences = [sima.Sequence.join(
... sima.Sequence.create('ndarray', cycle1_channel1),
... sima.Sequence.create('ndarray', cycle1_channel2)),
... sima.Sequence.join(
... sima.Sequence.create('ndarray', cycle2_channel1),
... sima.Sequence.create('ndarray', cycle2_channel2))]
>>> dataset = sima.ImagingDataset(
... sequences, 'example_np.sima', channel_names=['green', 'red'])
Multipage TIFF files¶
For simplicity, we consider the case of only a single cycle and channel.
>>> sequences = [sima.Sequence.create('TIFF', 'example_Ch1.tif')]
>>> dataset = sima.ImagingDataset(sequences, 'example_TIFF.sima')
HDF5 files¶
The argument ‘yxt’ specifies that the first index of the HDF5 array corresponds to the row, the second to the column, and the third to the time.
>>> sequences = [sima.Sequence.create('HDF5', 'example.h5', 'yxt')]
>>> dataset = sima.ImagingDataset(sequences, 'example_HDF5.sima')
Loading ImagingDataset objects¶
A dataset object can also be loaded from a saved path with the .sima extension.
>>> dataset = sima.ImagingDataset.load('example.sima')
Motion correction¶
The SIMA package implements a variety of approaches for motion correction. To use one of these approaches, the user first creates an object encapsulating the approach and a choice of parameters. For example, an object encapsulating the approach of translating imaging planes in two dimensions with a maximum displacement of 15 rows and 30 columns can be created as follows:
>>> import sima.motion
>>> mc_approach = sima.motion.PlaneTranslation2D(max_displacement=[15, 30])
Once this object encapsulating the approach has been created, its correct() method can then be applied to create a corrected dataset from a list of sequences.
>>> sequences = [sima.Sequence.create('TIFF', 'example_Ch1.tif')]
>>> dataset = mc_approach.correct(sequences, 'example_translation2D.sima')
In the following example, the SIMA package is used for motion correction based on a hidden Markov model (HMM). The sima.motion.HiddenMarkov2D class takes the initialize arguments to specify its parameters (e.g.an optional argument is used to indicate that the maximum possible displacement is 20 rows and 30 columns). The correct() method takes the same arguments as are used to initialize an imaging dataset object, as well as some additional optional arguments.
>>> mc_approach = sima.motion.HiddenMarkov2D(
... granularity='row', max_displacement=[20, 30], verbose=False)
>>> dataset = mc_approach.correct(sequences, 'example_mc.sima')
When the signal is of interest is very sparse or highly dynamic, it is sometimes helpful to use a second static channel to estimate the displacements for motion correction. The example below is for the case where the first channel contains a dynamic GCaMP signal whose large variations would confuse the motion correction alogorithm, and the second channel contains a static tdTomato signal that provides a stable reference.
>>> sequences = [
... sima.Sequence.join(sima.Sequence.create('TIFF', 'example_Ch1.tif'),
... sima.Sequence.create('TIFF', 'example_Ch1.tif'))
... ]
>>> dataset = mc_approach.correct(
... sequences, 'example_mc2.sima', channel_names=['GCaMP', 'tdTomato'],
... correction_channels=['tdTomato'])
When motion correction is invoked as above, only the tdTomato channel is used for estimating the displacements, which are then applied to both channels.
Segmentation and ROIs¶
Automated segmentation¶
The SIMA package implements a number of approaches for automated segmentation. To use one of the approaches, the first step is to create an object encapsulating the approach and a choice of parameters. For example, an object representing the approach of segmenting a single-plane dataset with spatiotemporal independent component analysis (STICA) can be created as follows:
>>> import sima.segment
>>> stica_approach = sima.segment.STICA(components=5)
We can also add post-processing steps to this approach, for example to convert the STICA masks into sparse regions of interest (ROIs), smooth their boundaries, and merge overlapping ROIs.
>>> stica_approach.append(sima.segment.SparseROIsFromMasks())
>>> stica_approach.append(sima.segment.SmoothROIBoundaries())
>>> stica_approach.append(sima.segment.MergeOverlapping(threshold=0.5))
Once the approch has been created, it can be passed as an argument to the segment() method of an An ImagingDataset. The segment() method can also take an optional label argument for the resulting set of ROIs.
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> rois = dataset.segment(stica_approach, 'auto_ROIs')
Editing, creating, and registering ROIs with ROI Buddy¶
Note that the an ImagingDataset object can be loaded with the ROI Buddy graphical user interface (GUI) for manual editing of existing the ROI lists, creation of new ROI lists, or registration of ROI lists across multiple experiments in which the same field of view is imaged. For more details, consult the ROI Buddy documentation.
Importing ROIs from ImageJ¶
ROIS can also be imported from ImageJ, as shown in the following example.
>>> from sima.ROI import ROIList
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> rois = ROIList.load('imageJ_ROIs.zip', fmt='ImageJ')
>>> dataset.add_ROIs(rois, 'from_ImageJ')
>>> dataset.ROIs.keys() # view the labels of the existing ROILists
['from_ImageJ', 'auto_ROIs']
Mapping ROIs between datasets¶
Sometimes, for example when imaging the same field of view over multiple days, one wishes to segment the same structures in separate ImagingDataset objects. If all of the ImagingDataset objects have been segmented, then the results of the segmentations can be registered with the ROI Buddy GUI as mentioned previously. If, however, only one of the datasets has been segmented, the results of the segmentation can be applied to the other datasets by applying to each ROI the affine transformation necessary to map one imaged field of view onto the other. This can be done either with the ROI Buddy GUI or with a call to the import_transformed_ROIs() method, whose arguments allow for specification of the channels used to align the two datasets, the label of the :obj:ROIList to be transformed from one dataset to the other, the label that will be applied to the new ROIList, and whether to copy the properties of the ROIs as well as their shapes.
>>> source_dataset = sima.ImagingDataset.load('example.sima')
>>> target_dataset = sima.ImagingDataset.load('example_mc2.sima')
>>> target_dataset.ROIs.keys()
[]
>>> target_dataset.import_transformed_ROIs(
... source_dataset, source_channel='green', target_channel='GCaMP',
... source_label='from_ImageJ', target_label='transformed',
... copy_properties='True')
>>> target_dataset.ROIs.keys()
['transformed']
This approach allows the user to focus on careful manual curation of the segmentation for a single ImagingDataset, with the results of this segmentation then applied to all datasets acquired at the same field of view.
Accessing stored ROIs¶
Whenever ROIs are created or imported, they are permanently stored as part of the ImagingDataset object. The ROIs can be recovered at any time using the label specified at the time when the ROIs were created.
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> dataset.ROIs.keys() # view the labels of the existing ROILists
['from_ImageJ', 'auto_ROIs']
>>> rois = dataset.ROIs['auto_ROIs']
Extraction¶
Once the ROIs have been edited and registered, the dataset can be loaded, and then dynamic fluorescence signals can be extracted from the ROIs with the extract() method.
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> dataset.ROIs.keys()
['from_ImageJ', 'auto_ROIs']
>>> rois = dataset.ROIs['from_ImageJ']
>>> dataset.channel_names
['red', 'green']
>>> signals = dataset.extract(
... rois, signal_channel='green', label='green_signal')
The extracted signals are permanently saved with the ImagingDataset object and can be accessed at any time with the command signals() method.
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> signals = dataset.signals(channel='green')['green_signal']
Exporting data¶
Data can be exported from the SIMA ImagingDataset objects at various stages of the analysis. This allows SIMA to be used for early stages of data analysis, and then for the exported data to be analyzed with separate software. If, however, further analysis is to be performed with Python, such exporting may not be necessary. The subsections below contain examples showing how to export image data and signal data.
Image data¶
The ImagingDataset class has two methods for exporting image data, export_frames() and export_averages(), which export either all the frames or the time averages of each channel, respectively. These methods can be used to view the results of motion correction, as shown in the following example.
>>> import sima.motion
>>> sequences = [sima.Sequence.create('TIFF', 'example_Ch1.tif')]
>>> dataset = mc_approach.correct(sequences, 'example_mc3.sima')
>>> dataset.export_averages(['avgs.tif'], fmt='TIFF16')
>>> dataset.export_frames([[['frames.tif']]], fmt='TIFF16')
The paths to which the exported data are saved are organized as a list with one filename per channel for the export_averages() method, or as a list of lists (organized analogously to the sequence used to initialize an ImagingDataset object) for the export_frames() method. If however, the export format is specified to HDF5, then the filenames for export_frames() should be organized into a list with one filename per cycle, since both channels are combined into a single HDF5 file.
>>> dataset.export_frames(['exported_frames.h5'], fmt='HDF5')
Signal data¶
For users wishing to analyze the extracted signals with an external program, these signals can be exported to a CSV file.
>>> dataset = sima.ImagingDataset.load('example.sima')
>>> dataset.export_signals('example_signals.csv', channel='green')
The resulting CSV file contains the id, label, and tags for each ROI, and the extracted signal from each ROI at each frame time.
Complete example¶
Below are the contents of workflow.py in the examples directory provided with the SIMA source code.
#! /usr/bin/env python
"""
This file provides a demonstration of how to use the SIMA package.
In order to run this file, you will need to download the file
http://www.losonczylab.org/workflow_data.zip and extract it in your
current working directory.
"""
from __future__ import print_function
from builtins import input
from builtins import range
##############################################################################
# #
# PART 0: Import SIMA and necessary submodules. #
# #
##############################################################################
import sima
import sima.motion
import sima.segment
##############################################################################
# #
# PART 1: Preparing the iterables. #
# #
##############################################################################
# Generate the filenames with Python list comprehensions.
tiff_filenames = [
['workflow_data/Cycle{n1:02d}_Ch{n2}.tif'.format(n1=cycle, n2=channel)
for channel in range(1, 3)] for cycle in range(1, 16)
]
# The resulting filenames are printed for clarification.
print("TIFF filenames:\n", tiff_filenames)
# Finally, we construct a MultiPageTIFF iterable using each of the filenames.
sequences = [
sima.Sequence.join(*[sima.Sequence.create('TIFF', chan) for chan in cycle])
for cycle in tiff_filenames]
##############################################################################
# #
# PART 2: Running motion correction to create the dataset, and exporting #
# the corrected image data. #
# #
##############################################################################
dataset_path = 'workflow_data/dataset.sima'
correction_approach = sima.motion.HiddenMarkov2D(num_states_retained=30,
max_displacement=[20, 30])
dataset = correction_approach.correct(
sequences, dataset_path, channel_names=['tdTomato', 'GCaMP'],
trim_criterion=0.95)
# Export the time averages for a manuscript figure.
dataset.export_averages(['workflow_data/tdTomato.tif',
'workflow_data/GCaMP.tif'])
# Generate the output filenames with Python list comprehensions.
output_filenames = [
[[channel.replace('.tif', '_corrected.tif') for channel in cycle]]
for cycle in tiff_filenames
]
# The resulting filenames are printed for clarification.
print("Output filenames:\n", output_filenames)
# Export the corrected frames for a presentation.
dataset.export_frames(output_filenames, fill_gaps=True)
# At this point, one may wish to inspect the exported image data to evaluate
# the quality of the motion correction before continuing.
while True:
input_ = input("Continue? (y/n): ")
if input_ == 'n':
exit()
elif input_ == 'y':
break
##############################################################################
# #
# PART 3: Running automated segmentation and editing results with the ROI #
# Buddy GUI. #
# #
##############################################################################
# Segment the field of view into ROIs using the method for CA1 pyramidal cells
# and parameters that were determined based on the imaging magnification.
segmentation_approach = sima.segment.PlaneCA1PC(
channel='GCaMP',
num_pcs=30,
max_dist=(3, 6),
spatial_decay=(3, 6),
cut_max_pen=0.10,
cut_min_size=50,
cut_max_size=150,
x_diameter=14,
y_diameter=7,
circularity_threhold=.5,
min_roi_size=20,
min_cut_size=40
)
dataset.segment(segmentation_approach, 'auto_ROIs')
# At this point, one may wish to edit the automatically segmented ROIs using
# the ROI Buddy GUI before performing signal extraction.
while True:
input_ = input("Continue? (y/n): ")
if input_ == 'n':
exit()
elif input_ == 'y':
break
##############################################################################
# #
# PART 4: Extracting fluorescence signals from the ROIs. #
# #
##############################################################################
# Reload the dataset in case any changes have been made with ROI Buddy
dataset = sima.ImagingDataset.load(dataset_path)
# Extract the signals. By default, the most recently created ROIs are used.
dataset.extract(signal_channel='GCaMP', label='GCaMP_signals')
# Export the extracted signals to a CSV file.
dataset.export_signals('example_signals.csv', channel='GCaMP',
signals_label='GCaMP_signals')
##############################################################################
# #
# PART 5: Visualizing data using Python. #
# #
##############################################################################
# import necessary functions from matplotlib
from matplotlib.pyplot import plot, show
# plot the signal from an ROI object, with a different color for each cycle
raw_signals = dataset.signals('GCaMP')['GCaMP_signals']['raw']
for sequence in range(3): # plot data from the first 3 cycles
plot(raw_signals[sequence][3]) # plot the data from ROI #3
show(block=True)