<jupyter_start><jupyter_text>Fly-QMA TutorialThis notebook is designed to walk the user through some of the basic functionalities offered by Fly-QMA. For more detailed explanations of all available functionalities please refer to the Fly-QMA documentation.Throughout this tutorial we'll use an example dataset; a single 3D confocal image stack containing three fluorescence reporter channels:- Channel 0: DAPI nuclear stain- Channel 1: HisRFP clonal marker- Channel 2: PntGFP control reporterNote that the control reporter should be expressed equally within a strip of cells near the morphogenetic furrow, but due to fluorescence bleedthrough from the clonal marker it shows variegated expression. We will correct for that later in the tutorial. To get started, let's import a handful of settings and packages to help out with some basic operations and visualization.<jupyter_code>from os.path import join
import numpy as np
from flyqma.visualization.settings import *
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load raw microscopy dataTo load our image stack, we need to define a path to an otherwise empty directory containing a multi-channel tif file.<jupyter_code>stack_path = '../data/replication_data/tutorial/stack'<jupyter_output><empty_output><jupyter_text>We may then instantiate a ```Stack``` object. This instance represents the entire 3D image stack. To limit unnecessary memory usage, the stack instance does not load the actual image data unless we explicitly tell it to using the ```load_image``` method. Note that we must correctly define the bit depth upon instantiation, e.g. whether our tif uses an 8-, 12-, or 16-bit container.<jupyter_code>from flyqma import Stack

stack = Stack(stack_path, bit_depth=12)
stack.load_image()<jupyter_output><empty_output><jupyter_text>Visualize layers within an image stack At this point we can readily select individual layers, either using the ```load_layer``` method or simply by indexing the ```Stack``` instance. Both of these operations return a ```Layer``` instance, which represents an individual layer within the image stack. The former is more verbose but offers a bit more functionality, such as a ```use_cache=True``` option that loads saved measurements directly from a file rather than re-processing (annotation, bleedthrough correction, etc.) them as they are loaded.<jupyter_code>layer_id = 0
layer = stack[layer_id]

# alternatively, with a bit more functionality
layer = stack.load_layer(layer_id)<jupyter_output><empty_output><jupyter_text>We can visualize a layer using the ```Layer.show``` method. See the documentation for details on all the various options.<jupyter_code>_ = layer.show(figsize=(8, 8))<jupyter_output><empty_output><jupyter_text>We may also visualize individual fluorescence channels within each layer by using the ```Layer.get_channel``` method, which returns an ```ImageScalar``` instance representing the individual channel. This instance provides all the same visualization options as it's parent multichannel layer.<jupyter_code>channel_id = 0
dapi = layer.get_channel(channel_id)

# visualize the DAPI stain
_ = dapi.show(figsize=(8, 8))<jupyter_output><empty_output><jupyter_text>Segment each layer to identify cell nuclei We can now segment each layer in the image stack. We can do this on a layer-by-layer basis using the ```Layer.segment``` method, or just use the ```Stack.segment``` method to perform the segmentation for all layers in the stack.Either way, we need to define which channel represents the nuclear stain, which in this case is "Channel 0" or DAPI. The segmentation algorithm used relies upon several seed detection strategies followed by a watershed segmentation. Parameters for each of these operations may be supplied as keyword arguments; see the Fly-QMA documentation for further details.<jupyter_code>segmented_channel_id = 0

stack.segment(channel=segmented_channel_id)<jupyter_output><empty_output><jupyter_text>Once the segmentation is complete, we can visualize the identified contours in the same manner as before. The ```segments``` argument determines whether or not individual contours are visualized.<jupyter_code>image = stack[layer_id].get_channel(segmented_channel_id)

_ = image.show(segments=True, figsize=(8, 8))<jupyter_output><empty_output><jupyter_text>Measure nuclear fluorescence intensities Following segmentation, we can aggregate the resultant fluorescence measurements from all layers using the ```Stack.aggregate_measurements``` method. The measurements are returned in a pandas dataframe.<jupyter_code>measurements = stack.aggregate_measurements()<jupyter_output><empty_output><jupyter_text>Save or Load an image stack We can always save our progress by using the ```Stack.save``` method to save the individual stack instance to file. This operation will store a number of new directories and files within the existing stack directory.<jupyter_code>stack.save()<jupyter_output><empty_output><jupyter_text>These items will be incorporated whenever we load the image stack from file using the same approach as before.<jupyter_code>from flyqma import Stack

stack_path = '../data/replication_data/tutorial/stack'
stack = Stack(stack_path)
stack.load_image()<jupyter_output><empty_output><jupyter_text>Save or load a collection of stacksOnce initialized, we can instantiate a collection of related image stacks (an "experiment") using the ```Experiment``` object. This object allow us to access child ```Stack``` instances much in the same way we accessed layers within a given stack.<jupyter_code>from flyqma import Experiment

experiment_path = '../data/replication_data/tutorial'
experiment = Experiment(experiment_path)

stack_id = 'stack'
stack = experiment.load_stack(stack_id, full=True)<jupyter_output><empty_output><jupyter_text>Define a region of interest Fly-QMA includes a simple tool for manually curating which regions of each layer are included in the analysis. The tool uses a simple GUI implemented via matplotlib. Use the GUI to manually draw a set of vertices enclosing the region of interest. Once running, the GUI offers both click and key commands:  Click on an image to add a vertex. Key commands:> **T**: remove last added point> **Y**: remove all points in layer> **E**: exclude entire layer> **W**: save selection to file> **Q**: disconnect and exit GUI Launch the GUI, then use the ```GUI.connect``` method to activate event handling so that any clicks and keystrokes are recognized. Unfortunately the GUI also requires us to switch our visualization protocol using the ```%matplotlib notebook``` magic. The ```W``` key will automatically save any progress to file within the corresponding stack directory.<jupyter_code>from flyqma.selection import GUI
%matplotlib notebook

gui = GUI(stack)
gui.connect()<jupyter_output><empty_output><jupyter_text>If a region of interest has previously been saved to file using the ```W``` key, we can re-load the same GUI using the ```GUI.load``` method. This also requires us to re-connect event handling.<jupyter_code>from flyqma.selection import GUI
%matplotlib notebook

gui = GUI.load(stack)
gui.connect()<jupyter_output><empty_output><jupyter_text>Any saved regions of interest will now be automatically incorporated into the stack's fluorescence measurements whenever they are loaded from file. The boolean ```selected``` column of the measurement dataframe indicates whether or not each nucleus resides within a region of interest. The same attribute will automatically be used to curate which data are used for bleedthrough correction and annotation.  Annotate clones within an image stackWe are now ready to annotate the image stack using Fly-QMA's unsupervised clustering algorithm. Annotation is performed by training a bivariate Gaussian mixture model to recognize patches of cells with distinct patterns of clonal marker fluorescence.  Train annotatorWe can train the mixture model for an individual layer, or we can use the ```Stack.train_annotator``` method to train a model across all measurements within an entire image stack. See the Fly-QMA documentation for details regarding the various keyword arguments used to control the annotation procedure. We can use the ```save=True``` argument to indicate whether we would like the annotation results to be saved to file.<jupyter_code>attribute = 'ch2' # this is the clonal marker fluorescence channel in this example

stack.train_annotator(attribute=attribute, save=True, max_num_components=8)<jupyter_output><empty_output><jupyter_text>Once saved, the corresponding annotation will be applied whenever the corresponding stack is loaded from file in the future, so long as the ```use_cache=True``` keyword argmument is used. Alternatively, we can explicitly apply the annotation to a given layer using the ```Layer.apply_annotation``` method.<jupyter_code># instantiate layer (without process=True)
layer_id = 0
layer = stack[layer_id]

# annotate layer
layer.apply_annotation()<jupyter_output><empty_output><jupyter_text>Inspect annotation resultsWe can review the annotator performance by using the ```Layer.build_attribute_mask``` method to overlay the labels on a given layer's clonal marker image. The ```interior_only``` argument determines whether clone boundaries are masked for ease of visualization, and the ```selection_only``` argument determines whether nuclei outside the region of interest are neglected.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

channel_id = 2
border_thickness = 1
interior_only = True
selection_only = False

# visualize clonal marker fluorescence
fig, axes = plt.subplots(ncols=2, figsize=(10, 5))
for ax in axes:
    fig = layer.get_channel(channel_id).show(segments=False, cmap='r', ax=ax)

# overlay assigned genotypes
layer.mark_boundaries(max_edges=border_thickness)
mask = layer.build_attribute_mask('genotype', interior_only=interior_only, selection_only=selection_only)
_ = axes[1].imshow(mask)<jupyter_output><empty_output><jupyter_text>Correct for fluorescence bleed-through First we can visualize the control reporter expression as we did before.<jupyter_code>control_reporter_channel = 1
control_reporter = layer.get_channel(channel=control_reporter_channel)
fig = control_reporter.show(segments=False, cmap='g', figsize=(5, 5))<jupyter_output><empty_output><jupyter_text>Bleedthrough correction must be applied on a layer-by-layer basis because each layer contains a different background region. We can fit a linear model to a layer's background bleedthrough trend using the ```LayerCorrection``` object. The ```xvar``` argument specifies the fluorescence channel from with bleedthrough originates, while the ```yvar``` argument specifies which channel is subject to bleedthrough. See the Fly-QMA documentation for additional detail regarding the remaining keyword arguments that control background pixel extraction and bleedthrough coefficient determination.<jupyter_code>from flyqma.bleedthrough.correction import LayerCorrection

# perform background correction
correction = LayerCorrection(layer, xvar=2, yvar=1, niters=20)<jupyter_output><empty_output><jupyter_text>Once the model is fit, we can visualize the background determination using the ```LayerCorrection.show_background_extraction``` method. We can also visualize the bleedthrough model fit using the ```LayerCorrection.show_correction``` method.<jupyter_code># visualize background extraction and model fit
_ = correction.show_background_extraction(invert=True, figsize=(2, 2))
_ = correction.show_correction(selected_only=False, figsize=(4, 1.25))<jupyter_output><empty_output><jupyter_text>Finally, if we save the bleedthrough model using the ```LayerCorrection.save``` method, the measurements will be automatically corrected whenever the layer is loaded from file in the future.<jupyter_code>correction.save()<jupyter_output><empty_output><jupyter_text>Analyze the processed dataOnce each stack has been segmented, corrected, and annotated, we can go ahead and aggregate all of the processed measurements from a collection of related image stacks.<jupyter_code>from flyqma import Experiment

experiment_path = '../data/replication_data/tutorial'
experiment = Experiment(experiment_path)

measurements = experiment.aggregate_measurements(selected_only=False, exclude_boundary=False)<jupyter_output><empty_output><jupyter_text>We can then use Fly-QMA's ```statistics.CloneComparison``` tool to evaluate whether fluorescence levels differ between clones. The ```basis``` keyword argument defines which channel is used in the comparison.<jupyter_code>from flyqma.analysis.statistics import CloneComparison

label_a, label_b = 0, 1

comparison = CloneComparison(measurements, label_a, label_b, basis='ch2', label='genotype')
comparison.plot()

# run statistical comparison
pvalue = comparison.compare(test='MW')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07144517628826036, 'fraction_numerical': 0.0037194885703215807, 'mean_word_length': 6.627659574468085, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 98, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '203145', 'license': 'MIT', 'n_tokens_mistral': 3493, 'n_tokens_neox': 3269, 'n_words': 1436, 'path': 'tutorial.ipynb', 'repo_name': 'sbernasek/flyqma', 'size': '12905'}