Merge pull request #230 from LiuDaveLiu/master

Update oralfacial_analyses and plotting

MATLAB/pipeline/+tracking/LickPortTracking3DBot.m

@@ -0,0 +1,48 @@
+%{
+# 3D trajectories of the lickport
+-> tracking.Tracking
+---
+lickport_x    :   longblob
+lickport_y    :   longblob
+lickport_z    :   longblob
+%}
+
+classdef LickPortTracking3DBot < dj.Computed
+    properties
+        keySource = (experiment.getSchema().v.Session & 'rig = "RRig-MTL"')...
+            & (tracking.getSchema().v.Tracking & 'tracking_device = "Camera 3"')...
+            & (tracking.getSchema().v.Tracking & 'tracking_device = "Camera 4"');
+    end
+    methods(Access=protected)        
+        function makeTuples(self, key)            
+            [bot_lickport_x, bot_lickport_y] = fetchn(tracking.getSchema().v.TrackingLickPortTracking & 'tracking_device = "Camera 4"' & key, 'lickport_x', 'lickport_y', 'ORDER BY trial');
+            [sid_lickport_x, sid_lickport_y] = fetchn(tracking.getSchema().v.TrackingLickPortTracking & 'tracking_device = "Camera 3"' & key, 'lickport_x', 'lickport_y', 'ORDER BY trial');
+
+            if length(bot_lickport_x) ~= length(sid_lickport_x)
+                disp('Mismatch in tracking bottom and side trials')
+                return
+            end
+
+            load('Calib_Results_stereo.mat')
+            key.tracking_device='Camera 4';
+            key_insert=repmat(key,1,length(bot_lickport_x));
+
+            counter=0;
+            for i = 1:length(bot_lickport_x) % loop over trials
+                x_left_1=[cell2mat(bot_lickport_x(i))'; cell2mat(bot_lickport_y(i))']; x_right_1=[cell2mat(sid_lickport_x(i))'; cell2mat(sid_lickport_y(i))'];
+
+                if size(x_left_1,2)==size(x_right_1,2)
+                    [Xc_1_left,~] = utils.stereo_triangulation(x_left_1,x_right_1,om,T,fc_left,cc_left,kc_left,alpha_c_left,fc_right,cc_right,kc_right,alpha_c_right);
+                    key_insert(i-counter).lickport_x = Xc_1_left(1,:);
+                    key_insert(i-counter).lickport_y = Xc_1_left(2,:);
+                    key_insert(i-counter).lickport_z = Xc_1_left(3,:);
+                    key_insert(i-counter).trial=i;
+                else
+                    key_insert=key_insert([1:i-counter-1 i-counter+1:end]);
+                    counter=counter+1;
+                end
+            end
+            insert(self,key_insert);
+        end
+    end
+end 

MATLAB/pipeline/populate.m

@@ -5,5 +5,6 @@
 init
 
 %% run the populate routines
-parpopulate(tracking.TongueTracking3DBot)
+% parpopulate(tracking.LickPortTracking3DBot)
+% parpopulate(tracking.TongueTracking3DBot)
 parpopulate(tracking.JawTracking3DSid)

pipeline/mtl_analysis/helper_functions.py

@@ -3,6 +3,7 @@
 from scipy import optimize
 
 import matplotlib.pyplot as plt
+plt.rcParams['font.size'] = 48
 import numpy as np
 
 # ======== Define some useful variables ==============
@@ -199,7 +200,7 @@ def get_trial_all(trial_tracks, trial_breathing, trial_whisking):
 
 def plot_tracking(session_key, unit_key,
                   tracking_feature='jaw_y', camera_key=_side_cam,
-                  trial_offset=0, trial_limit=10, xlim=(0, 5), axs=None):
+                  trial_offset=0, trial_limit=3, xlim=(0, 5), axs=None):
     """
     Plot jaw movement per trial, time-locked to trial-onset, with spike times overlay
     :param session_key: session where the trials are from
@@ -251,8 +252,8 @@ def get_trial_track(trial_tracks):
 
     fig = None
     if axs is None:
-        fig, axs = plt.subplots(3, 3, figsize=(16, 16))
-    assert len(axs) == 9
+        fig, axs = plt.subplots(3, 3, figsize=(84, 12))
+    assert len(axs.flatten()) == 9
 
     h_spacing = 150
     for trial_tracks, ax, ax_name, spk_color in zip(
@@ -265,10 +266,10 @@ def get_trial_track(trial_tracks):
             ax.set_yticks([])
             continue
         for tr_id, (trk_feat, tongue_out_bool, spike_times, tvec) in enumerate(get_trial_track(trial_tracks)):
-            ax.plot(tvec, trk_feat + tr_id * h_spacing, '.k', markersize=1)
-            ax.plot(tvec[tongue_out_bool], trk_feat[tongue_out_bool] + tr_id * h_spacing, '.', color='lime', markersize=2)
+            ax.plot(tvec, trk_feat + tr_id * h_spacing, '.k', markersize=8)
+            ax.plot(tvec[tongue_out_bool], trk_feat[tongue_out_bool] + tr_id * h_spacing, '.', color='lime', markersize=8)
             ax.plot(spike_times, np.full_like(spike_times, trk_feat[tongue_out_bool].mean() + h_spacing/3) + tr_id * h_spacing,
-                        color=spk_color, marker='$I$', linestyle='None', markersize=6)
+                        color=spk_color, marker='$I$', linestyle='None', markersize=16)
             ax.set_title(ax_name)
             ax.axvline(x=0, linestyle='--', color='k')
 
@@ -285,7 +286,7 @@ def get_trial_track(trial_tracks):
 
     return fig
 
-def plot_breathing(session_key, unit_key, trial_offset=0, trial_limit=10, xlim=(0, 5), axs=None):
+def plot_breathing(session_key, unit_key, trial_offset=0, trial_limit=3, xlim=(0, 5), axs=None):
     """
     Plot breathing per trial, time-locked to trial-onset, with spike times overlay
     :param session_key: session where the trials are from
@@ -318,7 +319,7 @@ def get_trial_breathing(trial_breathing):
 
     fig = None
     if axs is None:
-        fig, ax = plt.subplots(1, 1, figsize=(16, 16))
+        fig, ax = plt.subplots(1, 1, figsize=(48, 4))
 
     h_spacing = 3000
 
@@ -519,4 +520,50 @@ def rose_plot(ax, angles, bins=16, density=None, offset=0, lab_unit="degrees", s
     if lab_unit == "radians":
         label = ['$0$', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$',
                   r'$\pi$', r'$5\pi/4$', r'$3\pi/2$', r'$7\pi/4$']
-        ax.set_xticklabels(label)
+        ax.set_xticklabels(label)
+
+def plot_glmfit(unit_key, num_time=2000):
+
+    taus = np.arange(-5,6)
+    t_vec=np.arange(num_time)*0.017
+    test_y, predict_y, test_x, r2_t, weights = (oralfacial_analysis.GLMFit() & unit_key).fetch('test_y','predict_y','test_x','r2_t','weights')
+    max_r2_idx=np.where(r2_t[0]==np.max(r2_t[0]))
+    weights_r=np.round(weights[0],decimals=2)
+
+    fig, ax = plt.subplots(1, 1, figsize=(24, 16))
+
+    ax.plot(t_vec, test_x[0][:num_time, 1]+6,color='k',label='jaw ' +str(weights_r[max_r2_idx[0][0],2]))
+    ax.plot(t_vec, test_x[0][:num_time, 4]+6,color='g',label='tongue '+str(weights_r[max_r2_idx[0][0],5]))
+    ax.plot(t_vec, test_x[0][:num_time, 6]+12,color='b',label='breathing ' +str(weights_r[max_r2_idx[0][0],7]))
+    ax.plot(t_vec, test_x[0][:num_time, 7]+14,color='r',label='whisking '+str(weights_r[max_r2_idx[0][0],8]))
+    test_y_roll=np.roll(test_y[0],taus[max_r2_idx[0][0]])
+    ax.plot(t_vec,test_y_roll[:num_time],color='k',label='test')
+    ax.plot(t_vec,predict_y[0][max_r2_idx[0][0]][:num_time],color='r',label='prediction')
+    ax.legend(loc='upper left')
+    ax.set_xlabel('s')
+    ax.set_title('r2 = '+ str(np.round(r2_t[0][max_r2_idx[0][0]],decimals=3)))
+
+    # cosmetic
+    ax.set_xlim((t_vec[0],t_vec[-1]))
+    ax.set_yticks([])
+    ax.spines['left'].set_visible(False)
+    ax.spines['right'].set_visible(False)
+    ax.spines['top'].set_visible(False)
+
+    return fig
+
+def plot_dir_tuning(unit_key):
+
+    dir_tuning=(oralfacial_analysis.DirectionTuning() & unit_key).fetch('direction_tuning')
+    fr_mat=[[dir_tuning[0][2],dir_tuning[0][5],dir_tuning[0][8]], [dir_tuning[0][1],dir_tuning[0][4],dir_tuning[0][7]], [dir_tuning[0][0],dir_tuning[0][3],dir_tuning[0][6]]]
+
+    fig, ax = plt.subplots(1, 1, figsize=(3, 3))
+    neg=ax.imshow(fr_mat,vmin = 0, cmap='Reds', interpolation='none')
+    fig.colorbar(neg,ax=ax, location='right', anchor=(0, 0.3), shrink=0.7)
+    ax.set_xticks([])
+    ax.set_yticks([])
+    ax.spines['left'].set_visible(False)
+    ax.spines['right'].set_visible(False)
+    ax.spines['top'].set_visible(False)
+
+    return fig

pipeline/oralfacial_analysis.py

@@ -206,7 +206,7 @@ def make(self, key):
 
         fs=(tracking.TrackingDevice & 'tracking_device="Camera 4"').fetch1('sampling_rate')
 
-        amp, phase=behavior_plot.compute_insta_phase_amp(session_traces_w, float(fs), freq_band=(5, 20))
+        amp, phase=behavior_plot.compute_insta_phase_amp(session_traces_w, float(fs), freq_band=(3, 25))
         phase = phase + np.pi
 
         # compute phase and MI
@@ -302,6 +302,12 @@ def make(self, key):
         session_traces_t_l_t[np.where((session_traces_s_l_t > tongue_thr) & (session_traces_b_l_t > tongue_thr))] = 1
         session_traces_t_l_t[np.where((session_traces_s_l_t <= tongue_thr) | (session_traces_b_l_t <= tongue_thr))] = 0
         session_traces_t_l_t = np.hstack(session_traces_t_l_t)
+
+        session_traces_s_l_f = np.vstack(session_traces_s_l_o)
+        session_traces_b_l_f = np.vstack(session_traces_b_l_o)
+        session_traces_t_l_f = session_traces_b_l_f
+        session_traces_t_l_f[np.where((session_traces_s_l_f > tongue_thr) & (session_traces_b_l_f > tongue_thr))] = 1
+        session_traces_t_l_f[np.where((session_traces_s_l_f <= tongue_thr) | (session_traces_b_l_f <= tongue_thr))] = 0
 
         # from 3D calibration
         traces_s = v_tracking.JawTracking3DSid & key & [{'trial': tr} for tr in trial_key_o]
@@ -311,9 +317,18 @@ def make(self, key):
         session_traces_s_y_o = stats.zscore(np.vstack(session_traces_s_y_o),axis=None)
         session_traces_s_x_o = stats.zscore(np.vstack(session_traces_s_x_o),axis=None)
         session_traces_s_z_o = stats.zscore(np.vstack(session_traces_s_z_o),axis=None)
-        session_traces_b_y_o = stats.zscore(np.vstack(session_traces_b_y_o),axis=None) 
-        session_traces_b_x_o = stats.zscore(np.vstack(session_traces_b_x_o),axis=None)
-        session_traces_b_z_o = stats.zscore(np.vstack(session_traces_b_z_o),axis=None)
+        session_traces_b_y_o = np.vstack(session_traces_b_y_o)
+        traces_y_mean=np.mean(session_traces_b_y_o[np.where(session_traces_t_l_f == 1)])
+        traces_y_std=np.std(session_traces_b_y_o[np.where(session_traces_t_l_f == 1)])
+        session_traces_b_y_o = (session_traces_b_y_o - traces_y_mean)/traces_y_std
+        session_traces_b_x_o = np.vstack(session_traces_b_x_o)
+        traces_x_mean=np.mean(session_traces_b_x_o[np.where(session_traces_t_l_f == 1)])
+        traces_x_std=np.std(session_traces_b_x_o[np.where(session_traces_t_l_f == 1)])
+        session_traces_b_x_o = (session_traces_b_x_o - traces_x_mean)/traces_x_std
+        session_traces_b_z_o = np.vstack(session_traces_b_z_o)
+        traces_z_mean=np.mean(session_traces_b_z_o[np.where(session_traces_t_l_f == 1)])
+        traces_z_std=np.std(session_traces_b_z_o[np.where(session_traces_t_l_f == 1)])
+        session_traces_b_z_o = (session_traces_b_z_o - traces_z_mean)/traces_z_std
 
         session_traces_s_y = session_traces_s_y_o[trial_key-1]
         session_traces_s_x = session_traces_s_x_o[trial_key-1]
@@ -538,4 +553,98 @@ def make(self, key):
         proc = process.run(video_files_l, proc=roi_data)
 
         self.insert1({**key, 'mot_svd': proc['motSVD'][1][:, :3]})
-
+
+@schema
+class ContactLick(dj.Computed):
+    definition = """
+    -> tracking.Tracking
+    ---
+    contact_times: mediumblob
+    """
+
+    key_source = experiment.Session & v_tracking.TongueTracking3DBot & v_tracking.LickPortTracking3DBot & ephys.Unit & 'rig = "RRig-MTL"'
+
+    def make(self, key):
+        ts = 0.0034
+        radius=1
+        ton_thr = 0.95
+
+        bot_ton_x, bot_ton_y, bot_ton_z,trials = (v_tracking.TongueTracking3DBot & key).fetch('tongue_x','tongue_y','tongue_z','trial',order_by = 'trial')
+        bot_tongue_l = (v_tracking.Tracking.TongueTracking & key & 'tracking_device = "Camera 4"' & [{'trial': tr} for tr in trials]).fetch('tongue_likelihood', order_by = 'trial')
+        sid_tongue_l = (v_tracking.Tracking.TongueTracking & key & 'tracking_device = "Camera 3"' & [{'trial': tr} for tr in trials]).fetch('tongue_likelihood', order_by = 'trial')
+        bot_lic_x, bot_lic_y, bot_lic_z = (v_tracking.LickPortTracking3DBot & key).fetch('lickport_x','lickport_y','lickport_z', order_by = 'trial')
+
+        bot_tongue_l = np.vstack(bot_tongue_l)
+        sid_tongue_l = np.vstack(sid_tongue_l)
+        likelihood = bot_tongue_l
+        likelihood[np.where((sid_tongue_l > ton_thr) & (bot_tongue_l > ton_thr))] = 1
+        likelihood[np.where((sid_tongue_l <= ton_thr) | (bot_tongue_l <= ton_thr))] = 0
+
+        bot_ton_x=np.vstack(bot_ton_x)
+        bot_ton_y=np.vstack(bot_ton_y)
+        bot_ton_z=np.vstack(bot_ton_z)
+        bot_lic_x=np.vstack(bot_lic_x)
+        bot_lic_y=np.vstack(bot_lic_y)
+        bot_lic_z=np.vstack(bot_lic_z)
+
+        trial_contact = []
+
+        for i in np.arange(np.size(bot_ton_x,axis=0)):
+            lickSpan=np.where(likelihood[i,:]==1)[0]
+            lickBreak=np.diff(lickSpan)
+            lickS=np.concatenate(([0], np.where(lickBreak>1)[0]+1))
+
+            contacts = []
+            if len(lickS)>1:
+                lickS1=lickSpan[lickS]
+                lickE1=np.concatenate((lickSpan[lickS[1:]-1], [lickSpan[-1]]))
+                lick_x_med=np.median(bot_lic_x[i,350:])
+                lick_y_med=np.median(bot_lic_y[i,350:])
+                lick_z_med=np.median(bot_lic_z[i,350:])       
+
+                for j in np.arange(len(lickS1)):
+                    xp=bot_ton_x[i,lickS1[j]:lickE1[j]]
+                    yp=bot_ton_y[i,lickS1[j]:lickE1[j]]
+                    zp=bot_ton_z[i,lickS1[j]:lickE1[j]]
+                    inside=np.where(((xp-lick_x_med)**2 + (yp-lick_y_med)**2 + (zp-lick_z_med)**2) < radius**2)
+                    if lickE1[j]-lickS1[j]>10 and lickE1[j]-lickS1[j]<35  and np.size(inside)>0:
+                        contacts.append(lickS1[j]*ts)              
+            trial_contact.append({**key, 'trial': trials[i], 'tracking_device': 'Camera 4', 'contact_times': contacts})
+
+        self.insert(trial_contact, ignore_extra_fields=True)
+
+@schema
+class DirectionTuning(dj.Computed):
+    definition = """
+    -> ephys.Unit
+    ---
+    direction_tuning: mediumblob
+    
+    """
+
+    key_source = experiment.Session & v_oralfacial_analysis.ContactLick & ephys.Unit & 'rig = "RRig-MTL"'
+
+    def make(self, key):
+        good_units=ephys.Unit * ephys.ClusterMetric * ephys.UnitStat & key & 'presence_ratio > 0.9' & 'amplitude_cutoff < 0.15' & 'avg_firing_rate > 0.2' & 'isi_violation < 10' & 'unit_amp > 150'
+        unit_keys=good_units.fetch('KEY')
+
+        contact_times, trials,water_port=(v_oralfacial_analysis.ContactLick * experiment.MultiTargetLickingSessionBlock.BlockTrial * experiment.MultiTargetLickingSessionBlock.WaterPort & key).fetch('contact_times','trial','water_port', order_by = 'trial')
+
+        unit_dir=[]
+        for unit_key in unit_keys: # loop for each neuron
+            all_spikes=(ephys.Unit.TrialSpikes & unit_key & [{'trial': tr} for tr in trials]).fetch('spike_times', order_by='trial')
+            direction_spk=np.zeros(9)
+            direction_lick=np.zeros(9)
+            for i in np.arange(len(trials)):
+                tr_fr=np.zeros(len(contact_times[i]))
+                dir_idx=int(water_port[i][-1])-1
+                for j in np.arange(len(tr_fr)):
+                    tr_fr[j], _ = np.histogram(all_spikes[i], bins=1, range=(contact_times[i][j]-.05, contact_times[i][j]+.1))
+                direction_spk[dir_idx]=direction_spk[dir_idx]+sum(tr_fr)
+                direction_lick[dir_idx]=direction_lick[dir_idx]+len(tr_fr)
+
+            direction_tun=direction_spk/direction_lick
+
+            unit_dir.append({**unit_key, 'direction_tuning': direction_tun})
+
+        self.insert(unit_dir, ignore_extra_fields=True)