<jupyter_start><jupyter_text>Live Tutorial 1b – Single-particle tracking in Python - Part II---------- Qbio Summer School 2021--------------```Instructor: <NAME>Author: <NAME>Contact Info: <EMAIL>Copyright (c) 2021 Dr. <NAME>. Dr. <NAME>, Will RaymondColorado State University.Licensed under MIT License.```   Abstract This notebook provides a pipeline to track single-molecule RNA in a cell.  At the end of the tutorial, the student is expected to acquire the computational skills to implement the following list of objectives independently. List of objectives1.   To use the scientific libraries needed to segment and track single molecules inside live cells.2.   To load the experimental data.3.   To perform a cell segmentation using [Cellpose](https://github.com/MouseLand/cellpose).4.   To perform particle tracking using [Trackpy](http://soft-matter.github.io/trackpy/v0.5.0/).<jupyter_code># Loading libraries
import matplotlib.pyplot as plt # Library used for plotting
import urllib.request # importing library to download data
import numpy as np # library for array manipulation
import pandas as pd # Library to work with data frames
import random # Library to generate random numbers
import ipywidgets as widgets # Importing library
from ipywidgets import interact, interactive, HBox, Layout, VBox #  importing modules and functions.
# Installing and updating libraries
%%capture
!pip uninstall scikit-image -y
!pip install -U scikit-image
!pip install trackpy 
import skimage # Library for image manipulation
from skimage.io import imread # sublibrary from skimage
import trackpy as tp # Library for particle tracking
# Import cellpose
%%capture
!pip install cellpose
from cellpose import models
from cellpose import plot<jupyter_output><empty_output><jupyter_text>Dowloading data Downloading images from Github<jupyter_code>urls = ['https://github.com/MunskyGroup/uqbio2021/raw/main/module_1/images/simulated_cell.tif'] # loading simulated data
print('Downloading file...')
urllib.request.urlretrieve(urls[0], './simulated_cell.tif')
figName = './simulated_cell.tif'
# Importing the video as a numpy array
video = imread(figName)  #Numpy array [T,Y,X,C]
# dimensions in the video
number_timepoints, y_dim, x_dim, number_channels = video.shape[0], video.shape[1], video.shape[2], video.shape[3] # obtaining the dimensions size<jupyter_output><empty_output><jupyter_text>Visualization<jupyter_code>#plt.figure(figsize=(7,7))
#plt.imshow(video[0,:,:,0],cmap='gray')
#plt.show()
#@title
from matplotlib import animation
from IPython.display import HTML
def make_movie(image_tensor, cmaps = ['Reds_r','Greens_r','Blues_r'], channels = True):
  '''
  returns an ipython.html obj to display an image tensor
  '''
  if channels:  # if channels = True we will make a 3 color / 3 plot movie
    fig,axes = plt.subplots(1,3,dpi=120,tight_layout=True,)
  else: # otherwise just one subplot is needed
    fig,axes = plt.subplots(1,1,dpi=120,tight_layout=True,) #tight layout sets the axis's apart
  #fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)
  if channels: #if channels set up inital images here
    im1 = axes[0].imshow(image_tensor[0,:,:,0], cmap=cmaps[0])  #get blank axes.image objects to maniuplate later
    im2 = axes[1].imshow(image_tensor[0,:,:,1], cmap=cmaps[1])
    im3 = axes[2].imshow(image_tensor[0,:,:,2], cmap=cmaps[2])
    for k in range(3):
      axes[k].set_title('Channel %i'%k)
  else:
    im1 = axes.imshow(image_tensor[0,:,:,:], cmap=cmaps[0]) #if merge just make one
    axes.set_title('Merge')
  plt.close()

  #here set up the function that is called by FuncAnimation, it returns the artist objects and updates the frame by i
  def movieFrame(i):
    ## Extract all the image channels
    #Define  frames
    if channels:  #get the channels if 3 color
      channel_1_img = image_tensor[i,:,:,0]
      channel_2_img = image_tensor[i,:,:,1]
      channel_3_img = image_tensor[i,:,:,2]
      images = [channel_1_img, channel_2_img, channel_3_img]  # concatenate channels
      axes_handles = [im1,im2,im3]  #Axes handles from the imshow
    else:
      channel_img = image_tensor[i,:,:,:]   
      images = [channel_img,]  # concatenate channels
      axes_handles = [im1,]  #Axes handles from the imshow

    #update the axes imshow objects with the new data per k channels
    for k,image_n in enumerate(images):  #for loop that returns an object and an index
      axes_handles[k].set_array(images[k])  # 
    return axes_handles  #return the object
  
  #here call the animation function from matplotlib and then return the final HTML video object to the user
  anim = animation.FuncAnimation(fig, movieFrame, frames=image_tensor.shape[0], interval=100, blit=True)
  return HTML(anim.to_html5_video())
make_movie(video)<jupyter_output><empty_output><jupyter_text>Cell segmentation<jupyter_code># RUN CELLPOSE
imgs_2D = video[0,:,:,0] # [T,Y,X,C]
use_GPU = models.use_gpu()

# DEFINE CELLPOSE MODEL
model = models.Cellpose(gpu=use_GPU, model_type='cyto') # model_type='cyto' or model_type='nuclei'

# Running the models
masks, flows, styles, diams = model.eval(imgs_2D, diameter=200, flow_threshold=None, channels=[0,0])
# Cellpose returns an array with numbers 0 to i, where i!=0 represent the mask for each cell in the image.
plt.figure(figsize=(5,5))
plt.imshow(masks,cmap='Greys')
plt.show()
# We select the specific mask
selected_masks= masks==3 # Change this value to the specific cell that you need to select

# Plotting
plt.figure(figsize=(5,5))
plt.imshow(selected_masks)
plt.show()
# Removing the background from the video.
video_removed_mask = np.einsum('ijkl,jk->ijkl', video, selected_masks) # Tensor Matrix multiplication
# Plotting the masked image.
plt.figure(figsize=(5,5))
plt.imshow(video_removed_mask[0,:,:,0],cmap='gray')
plt.show()<jupyter_output><empty_output><jupyter_text>Particle tracking using Trackpy The particle tracking is performed using [trackpy](http://http://soft-matter.github.io/trackpy/v0.5.0/)The library documentation can be accessed in the following [link](https://buildmedia.readthedocs.org/media/pdf/trackpy/v0.2.3/trackpy.pdf).  Manual selection. (Initial step)  Visualizing the distribution for the particles' intensity<jupyter_code># This section will generate an histograme with the intensity of the detected particles in the video.

particle_size = 5 # according to the documentation must be an even number 3,5,7,9 etc.
minimal_intensity_for_selection = 0 # minimal intensity to detect a particle.

# "f" is a pandas data freame that contains the infomation about the detected spots
f = tp.locate(video_removed_mask[0,:,:,0], particle_size, minmass=minimal_intensity_for_selection) 

plt.rcParams["figure.figsize"] = (5,5)
fig, ax = plt.subplots()
ax.hist(f['mass'], bins=100, color = "orangered", ec="orangered")
ax.set(xlabel='mass', ylabel='count');
ax.set_ylim([0,200])
plt.show()<jupyter_output><empty_output><jupyter_text>Select an intensity threshold value and particle size In the following widget, the user is asked to manually input the particle size and the Intensity threshold, these numbers will be used for during the rest of the code to calculate intensities. Notice that tracking is only performed on channel 0.<jupyter_code># To start visualization move the time slider.
plt.rcParams["figure.figsize"] = (10,10) # if movie is too big, change size to (7,7)
def figure_viewer_tr(time,mass_text, drop_size):
    ch = 0  
    f = tp.locate(video_removed_mask[time,:,:,ch],drop_size, minmass=mass_text,maxsize=7,percentile=60) # "f" is a pandas data freame that contains the infomation about the detected spots
    tp.annotate(f,video_removed_mask[time,:,:,ch]);  # tp.anotate is a trackpy function that displays the image with the detected spots  
values_size=[3,5,7,9] # Notice value must be an EVEN number.
interactive_plot_tr = interactive(figure_viewer_tr,mass_text = widgets.IntText(value=1500,min=10,description='min Intensity'),drop_size = widgets.Dropdown(options=values_size,value=7,description='Particle Size'),time = widgets.IntSlider(min=0,max=video.shape[0]-1,step=1,value=0,description='Time'))
controls = HBox(interactive_plot_tr.children[:-1], layout = Layout(flex_flow='row wrap'))
output = interactive_plot_tr.children[-1]
display(VBox([controls, output]))<jupyter_output><empty_output><jupyter_text>Saving the threshold parameters that will be used to define the spot selection criteria<jupyter_code># This section saves the parameters adjusted in the previous widget in two variables that will be use for the rest of the code.
selected_intensity = interactive_plot_tr.kwargs_widgets[1].value
selected_size = interactive_plot_tr.kwargs_widgets[2].value<jupyter_output><empty_output><jupyter_text>Detecting the spots in all frames<jupyter_code># "f" is a pandas data freame that contains the infomation about the detected spots.
# tp.batch is a trackpy function that detects spots for multiple frames in a video.
f = tp.batch(video_removed_mask[:,:,:,0],selected_size, minmass=selected_intensity,percentile=70)<jupyter_output>Frame 49: 33 features<jupyter_text>Linking all detected spots accross all frames<jupyter_code># tp.link_df is a trackpy function that links spots detected in multiple frames, this generates the spots trajectories in time.
t = tp.link_df(f,5, memory=2) # tp.link_df(data_frame, min_distance_particle_moves, min_time_particle_vanishes).<jupyter_output>Frame 49: 39 trajectories present.<jupyter_text>Eliminating spurious trajectories.<jupyter_code>#trackpy.filtering.filter_stubs(tracks, threshold=100)
t1 = tp.filter_stubs(t, 10)  # selecting trajectories that appear in at least 10 frames.
# Compare the number of particles in the unfiltered and filtered data.
print('Before:', t['particle'].nunique())
print('After:', t1['particle'].nunique())<jupyter_output>Before: 99
After: 54<jupyter_text>Particle trajectories are stored in a Pandas data frame [Pandas data Library ](https://pandas.pydata.org) and a [link](https://pandas.pydata.org/docs/user_guide/10min.html) with a 10min guide to learn pandas.<jupyter_code># Pandas data frame
t1 
# Showing information for a selected_particle
selected_particle = 1
t1.loc[t1['particle']==selected_particle ]
# Showing information for a selected_frame
selected_frame = 10
t1.loc[t1['frame']==10]
# Extracting the intensity values for selected_particle
selected_particle = 1
t1.loc[t1['particle']==selected_particle].mass.values[:]
# Extracting the intensity values for particles with intensity higher than a threshold and in frame 1
selected_particle = 1
t1.loc[(t1['mass']>=3000) & (t1['frame']==1)  ]<jupyter_output><empty_output><jupyter_text>Exporting tracks to CSV.<jupyter_code># Optional section that saves the particles trajectories and intensities as a CSV file
t1.to_csv(r'tracking_particles.csv', index = False)<jupyter_output><empty_output><jupyter_text>Plotting trajectories in x-y.<jupyter_code># Plotting 
plt.figure(figsize=(10,10))
for i in range(0,t1['particle'].nunique() ):
  x_val = t1.loc[t1['particle']==i ].x.values[:]
  y_val = t1.loc[t1['particle']==i ].y.values[:]
  plt.plot(x_val,y_val, '-',linewidth = 3 )
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.title('Particle movement 2D')
plt.show()<jupyter_output><empty_output><jupyter_text>Plotting trajectories versus intensity using trackpy Notice that the particle number is not consecutive, for that reason the indexing is as follows:t1['particle'].unique()   Where .unique() is a method that gives a numpy array with the particle names.<jupyter_code># Deffining the number of detected spots.
n_particles = t1['particle'].nunique()

# Plotting
plt.figure(figsize=(20,5))

# plotting intensities in au
for id in range(0,n_particles):
    plt.plot(t1.loc[t1['particle']==t1['particle'].unique()[id]].frame, t1.loc[t1['particle']==t1['particle'].unique()[id]].mass , color ='orangered')
plt.xlabel('Time (sec)')
plt.ylabel('Intensity a.u.')
plt.title('Intensity vs Time')
plt.show()<jupyter_output><empty_output><jupyter_text>References Allan, <NAME>., <NAME>, <NAME>., <NAME>, <NAME>., & Verweij, <NAME>. (2021, April 13). soft-matter/trackpy: Trackpy v0.5.0 (Version v0.5.0). Zenodo. http://doi.org/10.5281/zenodo.4682814 <NAME>., <NAME>., <NAME>. and <NAME>., 2021.Cellpose: a generalist algorithm for cellular segmentation. Nature Methods, 18(1), pp.100-106. https://github.com/MouseLand/cellpose<jupyter_code><jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup