In [ ]:

    

# load the custom modules
%load_ext autoreload
%autoreload 2

import os, sys, inspect
import matplotlib.pylab as plt
import numpy as np

# Add source code directory to path to enable module import
curr_frame = inspect.getfile(inspect.currentframe())
curr_dir = os.path.dirname(os.path.abspath(curr_frame))
parent_dir = os.path.dirname(curr_dir)
module_dir = os.path.join(parent_dir, 'src')
os.sys.path.insert(0, module_dir)

import parse_trackmate as pt
import diffusion as dif


    

In [ ]:

    

# Load data
import glob

# Path to raw TrackMate data
data_dir = '../data/raw/IRE1_foci_diffusion'

data_files = sorted(glob.glob(os.path.join(data_dir,'*.xml')))
stripped_names = [os.path.basename(f) for f in data_files]
file_names = [os.path.splitext(fn)[0] for fn in stripped_names]

# Read data from the provided source files
raw_data = [pt.read_trackmate_file(f) for f in data_files]


    

In [ ]:

    

# Parse the data into tracks, creating lists of spots that
# make up each individual track 

# Store all track data
processed_tracks_by_movie = []
spots_by_movie = []

# Process each trace
for file, trace in zip(data_files, raw_data):
    # Parse the data into tracks, creating lists of spots that
    # make up each individual track 
    raw_tracks = pt.get_raw_tracks(trace)
    spots = pt.list_spots(trace)
    processed_tracks = pt.build_spot_lists_by_track(raw_tracks, spots)
    track_coords, track_edges, track_spot_ids = processed_tracks
    
    # Get intensity values for each frame of each track
    all_track_intensities = pt.pull_property_by_track('TOTAL_INTENSITY',
                                                   track_spot_ids, spots)
   
    # Check that spot radii are constant (otherwise comparison of 
    # mean intensities is meaningless)
    all_spot_radii = pt.pull_property_by_track('RADIUS',
                                                track_spot_ids, spots)
    radii_norm = np.concatenate(all_spot_radii - all_spot_radii[0][0])
    diff_radii = np.count_nonzero(radii_norm)
    if diff_radii > 0: print('WARNING: spot sizes are not uniform')

    processed_tracks_by_movie.append(processed_tracks)
    spots_by_movie.append(spots)


    

In [ ]:

    

# Clean up garbage (the dom files can take up lots of RAM)
import gc

del raw_tracks
del raw_data
gc.collect()


    

In [ ]:

    

# Plot all MSDs (one plot per movie)

plot_figs = False
save_figs = True
save_dir = '../reports/figures/IRE1_foci_tracking/MSD-all'

# Input imaging parameters:
frame_interval = 5 # in seconds

# Fitting parameters
fit_dc = True # Fit diffusion constant
diff_dim = 2 # 1D, 2D, or 3D diffusion
dc_fit_nframes = 10 # number of frames to fit to get diffusion constant

# interactive plotting
plt.close('all')
if plot_figs: 
    %matplotlib
    plt.ion()
else:
    plt.ioff()
    
# create save figure dirs
if save_figs:
    save_dir_svg, save_dir_png = dif.make_fig_save_dirs(save_dir)

# Calculate and plot a separate MSD graph for each movie

for curr_tracks, file_name in zip(processed_tracks_by_movie, file_names):
    
    track_coords, track_edges, track_spot_ids = curr_tracks
    
    # perform all MSD calculations
    t_dsq, msd_data = dif.calc_msd(track_coords, frame_interval)
    time, mean_dsq, std_dsq, sterr_dsq = msd_data
    
    if fit_dc:
        fit_params = dif.fit_diffusion_const(msd_data, diff_dim, dc_fit_nframes)
    else: 
        fit_params = None

    # Plot results
    fig, axarr = dif.plot_msd(t_dsq, msd_data, fit_params=fit_params)
    
    fig.canvas.set_window_title(file_name)
    plt.suptitle(file_name)
    if save_figs:
        dif.save_fig(fig, file_name, save_dir_svg, save_dir_png)
    if not plot_figs:
        del fig
        del axarr
        gc.collect()

if plot_figs:
    plt.show()
else:
    plt.close('all')
    print('done!')


    

In [ ]:

    

# Plot MSD sorted by spot intensity

plot_figs = False
save_figs = True
save_dir = '../reports/figures/IRE1_foci_tracking/MSD-by-Int'

# Input imaging parameters:
frame_interval = 5 # in seconds
num_intensity_bins = 20 # number of separate MSD plots to make

# Fitting parameters
fit_dc = True # Fit diffusion constant
fits = [] # storing fit results for each bin
diff_dim = 2 # 1D, 2D, or 3D diffusion
dc_fit_nframes = 10 # number of frames to fit to get diffusion constant
    
# create save figure dirs
if save_figs:
    save_dir_svg, save_dir_png = dif.make_fig_save_dirs(save_dir)

if plot_figs:
    %matplotlib
    
# Pool all tracks together
all_track_c = [] # coordinates (list of x,y,t numpy arrays)
all_track_int = [] # intensities

for curr_tracks, spots in zip(processed_tracks_by_movie, spots_by_movie):
    
    track_coords, track_edges, track_spot_ids = curr_tracks
    
    # Get intensity values for each frame of each track
    track_intensities = pt.pull_property_by_track('TOTAL_INTENSITY',
                                                   track_spot_ids, spots)
    all_track_c = all_track_c + track_coords
    all_track_int = all_track_int + track_intensities

# Get mean and starting intensities
mean_ints = [np.mean(x) for x in all_track_int]
start_ints = [x[0] for x in all_track_int]

# Bin all tracks by intensities
ints = np.array(mean_ints, dtype=float) # or replace with start_ints
ints_s = sorted(ints)
bin_size = int(len(ints) / (num_intensity_bins))
bin_cutoffs = ints_s[0:len(ints):bin_size]
bin_cutoffs[0] = min(ints_s)-1
bin_cutoffs[-1] = max(ints_s)

binned_coords = []
binned_intensities = []
for x in range(num_intensity_bins):
    lower = bin_cutoffs[x]
    upper = bin_cutoffs[x+1]
    inds = (ints > lower) & (ints <= upper)
    coords_bin = [x for i,x in enumerate(all_track_c) if inds[i]]
    ints_bin = [x for i,x in enumerate(ints) if inds[i]]
    binned_coords.append(coords_bin)
    binned_intensities.append(ints_bin)

for i, curr_coords in enumerate(binned_coords):
    # perform all MSD calculations
    t_dsq, msd_data = dif.calc_msd(curr_coords, frame_interval)
    time, mean_dsq, std_dsq, sterr_dsq = msd_data
    
    if fit_dc:
        fit_params = dif.fit_diffusion_const(msd_data, diff_dim, dc_fit_nframes)
        fits.append(fit_params)
    else: 
        fit_params = None
    
    # Plot results
    title = 'Intensities: {0:.2f} to {1:.2f}'.format(bin_cutoffs[i], 
                                                     bin_cutoffs[i+1])
    filename = 'Intensities_{0:.2f}_to_{1:.2f}'.format(bin_cutoffs[i], 
                                                     bin_cutoffs[i+1])
    if plot_figs or save_figs:
        fig, axarr = dif.plot_msd(t_dsq, msd_data, fit_params)

        fig.canvas.set_window_title(title)
        plt.suptitle(title)
    if save_figs:
        dif.save_fig(fig, filename, save_dir_svg, save_dir_png)
    
    if not plot_figs:
        del fig
        del axarr
        gc.collect()
    
if plot_figs:
    plt.show()
else:
    plt.close('all')
    print('done!')


    

In [ ]:

    

# Plot diffusion constants vs. cluster size
%matplotlib
save_figs = True
save_dir = '../reports/figures/IRE1_foci_tracking/MSD-by-Int-summary'
filename = 'Summary'

if save_figs:
    save_dir_svg, save_dir_png = dif.make_fig_save_dirs(save_dir)

if fit_dc:
    bin_ints = [np.mean(x) for x in binned_intensities]
    dcs = [x['dc'] for x in fits]
    stdevs = [x['stdev'] for x in fits]
    fig = plt.figure()
    plt.errorbar(bin_ints, dcs, yerr=stdevs, fmt='o')
    plt.ylim([0, max(dcs)*1.2])
    plt.xlabel('Cluster intensity (A.U.)')
    plt.ylabel(r'Diffusion constant ($\mu m^2/s$)')
    dif.save_fig(fig, filename, save_dir_svg, save_dir_png)


    

In [ ]:

    

# Plot cumulative distribution function of any desired property

from scipy.stats import truncnorm, norm

prop = 'TOTAL_INTENSITY'
#prop = 'MEAN_INTENSITY'

save_figs = True
plot_pdf = True
save_dir = '../reports/figures/IRE1_foci_tracking/MSD-by-Int-summary'
filename = 'CDF'
filename_pdf = 'PDF'

if save_figs:
    save_dir_svg, save_dir_png = dif.make_fig_save_dirs(save_dir)

all_prop = []
for curr_tracks, spots in zip(processed_tracks_by_movie, spots_by_movie):
    track_coords, track_edges, track_spot_ids = curr_tracks
    # Get intensity values for each frame of each track
    track_prop = pt.pull_property_by_track(prop, track_spot_ids, spots)
    all_prop = all_prop + track_prop
#print(spots[0]) # uncomment to see what options are available

# Get mean and starting intensities
data_us = [np.mean(x) for x in all_prop]
data = sorted(data_us)

plt.rcParams["figure.figsize"] = [10,8]
fig = plt.figure()
x = np.array(data) / float(max(data))
x = x - min(x)
print(max(x))
x_sim = np.linspace(0,1,100)
y = np.arange(len(data))
print(len(range(len(data))))
#x = np.cumsum(data)

# simulate truncated gaussian distribution for comparison
sim_std = np.std(x)
sim_mean = np.mean(x)
sim_max = max(x)
cdf_ceil = len(x)
a, b = (0 - sim_mean) / sim_std, (sim_max - sim_mean) / sim_std
sim_dist = cdf_ceil*truncnorm.cdf(x_sim,a,b, loc=sim_mean, scale=sim_std)
#sim_dist = cdf_ceil*norm.cdf(x_sim, loc=sim_mean, scale=sim_std)
mean, var = truncnorm.stats(a, b, loc=sim_mean, scale=sim_std, moments='mv')

print(mean,var, np.mean(x), np.var(x))

plt.step(x,y, linewidth=2.0, label='IRE1 foci')
line2 = plt.plot(x_sim, sim_dist, 'r--', label='Simulated normal distribution')
plt.setp(line2, linewidth=2.0)
plt.title('CDF for ' + prop)
plt.xlabel('Normalized ' + prop)
plt.ylabel('Cumulative counts')
plt. legend(loc='lower right', shadow=True)

if save_figs:
    dif.save_fig(fig, filename, save_dir_svg, save_dir_png)

if plot_pdf:
    data_norm = x - min(x)
    bottom_lim = (max(x)-min(x)) * -0.01
    
    fig2 = plt.figure()
    ax2 = plt.axes()
    
    num_bins = 50
    counts, bins = np.histogram(x, bins=num_bins)
    bins = bins[:-1] + (bins[1] - bins[0])/2
    probs = counts/float(counts.sum())
    plt.bar(bins, probs, 1.0/num_bins, label='IRE1 foci')
    
    #plt.hist(x, bins=50, normed=True, label='IRE1 foci')
    #pdf_norm = max(probs) # np.sum(data_norm)
    pdf_norm = 1.0 / num_bins
    sim_pdf = pdf_norm * truncnorm.pdf(x_sim,a,b, loc=sim_mean, scale=sim_std)
    #sim_pdf = pdf_norm * norm.pdf(x_sim, loc=sim_mean, scale=sim_std)
    plt.plot(x_sim, sim_pdf,'r--', linewidth=2.0, label='Simulated normal distribution')
    plt.title('PDF for ' + prop)
    plt.xlabel('Normalized ' + prop)
    plt.ylabel('Probability')
    ax2.set_xlim(left=0)
    #ax2.set_ylim(bottom=bottom_lim)
    plt.legend(loc='upper center', shadow=True)
    
    if save_figs:
        dif.save_fig(fig2, filename_pdf, save_dir_svg, save_dir_png)

plt.show()


    

In [ ]:

    

# Plot trajectories
%matplotlib

# Various plotting tools
f2, axarr = plt.subplots(1,2)
axarr[0].set_title('All tracks')
axarr[0].set_xlabel('Position (μm)')
axarr[0].set_ylabel('Position (μm)')
axarr[1].set_xlabel('Frame')
axarr[1].set_ylabel('Intensity (A.U.)')
axarr[1].set_title('All intensities')
for i, track in enumerate(track_coords):
    axarr[0].plot(track[:,0], track[:,1])
    axarr[1].plot(track[:,2], all_track_intensities[i])
plt.show()
#plt.plot(track_coords[testind][:,2], all_track_intensities[testind])