added code for all the plots in plot_photometry. also added the preprocessing needed for the new plots within the convert_photometry.py file. There is still some things that need to be rewritten given the new nwb converted data. Basically anything that calls the previous sampledata dataframe and visinds

Author: j2figuer

src/jdb_to_nwb/convert_photometry.py

@@ -10,6 +10,9 @@
 from scipy.sparse import diags, eye, csc_matrix
 from scipy.sparse.linalg import spsolve
 from sklearn.linear_model import Lasso
+from scipy.stats import linregress
+
+from plot_photometry import .
 
 # Some of these imports are unused for now but will be used for photometry metadata
 from ndx_fiber_photometry import (
@@ -430,22 +433,105 @@ def process_ppd_photometry(nwbfile: NWBFile, ppd_file_path):
     relative_raw_signal = raw_green / raw_405   
 
     sampling_rate = ppd_data['sampling_rate']
+    time_seconds = ppd_data['time']/1000
     visits = ppd_data['pulse_inds_1'][1:]
 
+    pulse_times_in_mins = [time / 60000 for time in visits]
+    pulse_times_in_seconds = np.array(pulse_times_in_mins) * 60
+
+    xvals = np.arange(0,len(raw_green))/sampling_rate/60 
+
+    plot_raw_ratiometric_and_565_photometry_signals(visits, xvals, pulse_times_in_mins, raw_green, raw_red, raw_405, relative_raw_signal, sampling_rate)
+
+    # Calculate the correlation between the signals. This I am not too sure if we want to keep or not. Best to discuss with Josh? 
+    slope_405x470, intercep_405x470, r_value_405x470 = linregress(x=raw_405, y=raw_green)
+    plot_405_470_correlation(raw_405, raw_green,slope_405x470, intercep_405x470, r_value_405x470)
+
+    slope_405x565, intercep_405x565, r_value_405x565 = linregress(x=raw_405, y=raw_red)
+    plot_405_565_correlation(raw_405, raw_red, slope_405x565, intercep_405x565, r_value_405x565)
+
+    slope_470x565, intercep_470x565, r_value_470x565 = linregress(x=raw_green, y=raw_red)
+    plot_470_565_correlation(raw_green, raw_red, slope_470x565, intercep_470x565, r_value_470x565)
+
     # low pass at 10Hz to remove high frequency noise
     print('Filtering data...')
     b,a = butter(2, 10, btype='low', fs=sampling_rate)
     green_denoised = filtfilt(b,a, raw_green)
     red_denoised = filtfilt(b,a, raw_red)
     ratio_denoised = filtfilt(b,a, relative_raw_signal)
     denoised_405 = filtfilt(b,a, raw_405)
+
     # high pass at 0.001Hz which removes the drift due to bleaching, but will also remove any physiological variation in the signal on very slow timescales.
     b,a = butter(2, 0.001, btype='high', fs=sampling_rate)
     green_highpass = filtfilt(b,a, green_denoised, padtype='even')
     red_highpass = filtfilt(b,a, red_denoised, padtype='even')
     ratio_highpass = filtfilt(b,a, ratio_denoised, padtype='even')
     highpass_405 = filtfilt(b,a, denoised_405, padtype='even')
 
+    # Plot the filtered signals of interest against each other 
+    slope_filtered, intercept_filtered, r_value_filtered = linregress(x=red_highpass, y=ratio_highpass)
+    plot_ratio_565_correlation(ratio_highpass, red_highpass, slope_filtered, intercept_filtered, r_value_filtered)
+
+    plot_interactive_filtered_signals(time_seconds, ratio_highpass, green_highpass, highpass_405, red_highpass)
+
+
+    time_window_in_sec = input('Enter the time window in seconds for the photometry signal alignment: ') # >:)
+    
+    samples_window = time_window_in_sec * sampling_rate
+    total_rewarded_trials = (rwd == 1).sum()
+    print("Total rewarded trials: "+str(total_rewarded_trials))
+
+    total_omitted_trials = (rwd == 0).sum()
+    print("Total omitted trials: "+str(total_omitted_trials))
+
+    # Initialize lists for aligned traces
+    aligned_ratio = []
+    aligned_green = []
+    aligned_405 = []
+    aligned_red = []
+
+#################################################################################
+    # QUESTION FOR STEPH:
+    # I need this code to run some of the plotting from this point forward.
+    # Is there a way to get this information from convert_behvaior?
+    # For now, I'll keep everything the way it is using "visinds" and "port_entry_indeces"
+    visinds = sampledata.loc[sampledata.port.notnull()].index.values
+    port_entry_indices = visinds 
+#################################################################################
+
+    # Align windows
+    for idx in visits:
+        if idx - samples_window >= 0 and idx + samples_window < len(time_seconds):  # Ensure indices are in bounds
+            aligned_ratio.append(ratio_highpass[idx - samples_window:idx + samples_window])
+            aligned_green.append(green_highpass[idx - samples_window:idx + samples_window])
+            aligned_405.append(highpass_405[idx - samples_window:idx + samples_window])
+            aligned_red.append(red_highpass[idx - samples_window:idx + samples_window])
+
+    # Convert to numpy arrays for easier averaging
+    aligned_ratio = np.array(aligned_ratio)
+    aligned_green = np.array(aligned_green)
+    aligned_405 = np.array(aligned_405)
+    aligned_red = np.array(aligned_red)
+
+    # Calculate means and SEMs
+    mean_ratio = np.mean(aligned_ratio, axis=0)
+    sem_ratio = np.std(aligned_ratio, axis=0) / np.sqrt(len(aligned_ratio))
+
+    mean_green = np.mean(aligned_green, axis=0)
+    sem_green = np.std(aligned_green, axis=0) / np.sqrt(len(aligned_green))
+
+    mean_405 = np.mean(aligned_405, axis=0)
+    sem_405 = np.std(aligned_405, axis=0) / np.sqrt(len(aligned_405))
+
+    mean_red = np.mean(aligned_red, axis=0)
+    sem_red = np.std(aligned_red, axis=0) / np.sqrt(len(aligned_red))
+
+    # Generate time axis for aligned data
+    aligned_time = np.linspace(-time_window_in_sec, time_window_in_sec, 2 * samples_window)
+
+
+    plot_signals_aligned_port_entry(port_entry_indices, aligned_time, mean_ratio, sem_ratio, mean_green, sem_green, mean_405, sem_405, mean_red, sem_red, total_rewarded_trials, total_omitted_trials)
+
     # Z-score of each signal to normalize the data
     print('Z-scoring data...')
     green_zscored = np.divide(np.subtract(green_highpass,green_highpass.mean()),green_highpass.std())
@@ -457,6 +543,110 @@ def process_ppd_photometry(nwbfile: NWBFile, ppd_file_path):
     ratio_zscored = np.divide(np.subtract(ratio_highpass,ratio_highpass.mean()),ratio_highpass.std())
     print('Done processing photometry data!')
 
+    plot_normalized_signals(xvals, pulse_times_in_mins, green_zscored, zscored_405, red_zscored, ratio_zscored)
+
+    ratio_traces = []
+    red_traces = [] 
+    time_window = int(time_window_in_sec)
+
+    time_window_xvals = np.arange(-time_window*86,time_window*86+1)/86
+
+    # Loop through each index in visinds
+    for i in visinds:
+        # Extract a time window of data around the current index
+        start_idx = i - 86 * time_window
+        end_idx = i + 86 * time_window
+        ratio_trace = sampledata.loc[start_idx:end_idx, "ratio_z_scored"].values
+        red_trace = sampledata.loc[start_idx:end_idx, "red_z_scored"].values
+
+        # Only use traces of the correct length (matching xvals)
+        if len(ratio_trace) == len(time_window_xvals):
+            ratio_traces.append(ratio_trace)
+
+        if len(red_trace) == len(time_window_xvals):
+            red_traces.append(red_trace)
+
+    # Calculate mean and SEM 
+    if ratio_traces: 
+        ratio_traces_mean = np.mean(ratio_traces, axis=0)
+        ratio_traces_sem = np.std(ratio_traces, axis=0) / np.sqrt(len(ratio_traces))
+    else:  # If no traces, create empty arrays
+        ratio_traces_rwd_mean = np.zeros_like(time_window_xvals)
+        ratio_traces_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM 
+    if red_traces: 
+        red_traces_mean = np.mean(red_traces, axis=0)
+        red_traces_sem = np.std(red_traces, axis=0) / np.sqrt(len(red_traces))
+    else:  # If no traces, create empty arrays
+        red_traces_mean = np.zeros_like(time_window_xvals)
+        red_traces_sem = np.zeros_like(time_window_xvals)
+
+    plot_ratio_and_565_signals_aligned_port_entry(rat, date, time_window, time_window_xvals, ratio_traces_mean, ratio_traces_sem, red_traces_mean, red_traces_sem, visinds) # 'rat' and 'date' are only used for the title of the plot
+
+    ratio_rwd_traces = []
+    ratio_om_traces = []
+    red_rwd_traces = []
+    red_om_traces = []
+
+    # Loop through each index in visinds
+    for i in visinds:
+        # Extract a time window of data around the current index
+        start_idx = i - 86 * time_window_in_sec
+        end_idx = i + 86 * time_window_in_sec
+        ratio_trace = sampledata.loc[start_idx:end_idx, "ratio_z_scored"].values
+        red_trace = sampledata.loc[start_idx:end_idx, "red_z_scored"].values
+
+        # Only use traces of the correct length (matching xvals)
+        if len(ratio_trace) == len(time_window_xvals):
+            # Check if the trial is rewarded or omitted
+            if sampledata.loc[i, 'rwd'] == 1:
+                ratio_rwd_traces.append(ratio_trace)  # Add to rewarded traces
+            else:
+                ratio_om_traces.append(ratio_trace)  # Add to omitted traces
+
+            # Only use traces of the correct length (matching xvals)
+        if len(red_trace) == len(time_window_xvals):
+            # Check if the trial is rewarded or omitted
+            if sampledata.loc[i, 'rwd'] == 1:
+                red_rwd_traces.append(red_trace)  # Add to rewarded traces
+            else:
+                red_om_traces.append(red_trace)  # Add to omitted traces
+
+    # Calculate mean and SEM for rewarded traces
+    if ratio_rwd_traces:  # If there are any rewarded traces
+        ratio_rwd_mean = np.mean(ratio_rwd_traces, axis=0)
+        ratio_rwd_sem = np.std(ratio_rwd_traces, axis=0) / np.sqrt(len(ratio_rwd_traces))
+    else:  # If no rewarded traces, create empty arrays
+        ratio_rwd_mean = np.zeros_like(time_window_xvals)
+        ratio_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for rewarded traces
+    if red_rwd_traces:  # If there are any rewarded traces
+        red_rwd_mean = np.mean(red_rwd_traces, axis=0)
+        red_rwd_sem = np.std(red_rwd_traces, axis=0) / np.sqrt(len(red_rwd_traces))
+    else:  # If no rewarded traces, create empty arrays
+        red_rwd_mean = np.zeros_like(time_window_xvals)
+        red_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for omitted traces
+    if ratio_om_traces:  # If there are any omitted traces
+        ratio_om_mean = np.mean(ratio_om_traces, axis=0)
+        ratio_om_sem = np.std(ratio_om_traces, axis=0) / np.sqrt(len(ratio_om_traces))
+    else:  # If no omitted traces, create empty arrays
+        ratio_om_mean = np.zeros_like(time_window_xvals)
+        ratio_om_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for omitted traces
+    if red_om_traces:  # If there are any omitted traces
+        red_om_mean = np.mean(red_om_traces, axis=0)
+        red_om_sem = np.std(red_om_traces, axis=0) / np.sqrt(len(red_om_traces))
+    else:  # If no omitted traces, create empty arrays
+        red_om_mean = np.zeros_like(time_window_xvals)
+        red_om_sem = np.zeros_like(time_window_xvals)
+
+    plot_ratio_and_565_signals_aligned_port_entry_separated_by_rewarded_or_omitted(rat, date, time_window_in_sec, time_window_xvals, ratio_rwd_mean, ratio_rwd_sem, ratio_om_mean, ratio_om_sem, red_rwd_mean, red_rwd_sem, red_om_mean, red_om_sem, total_rewarded_trials, total_omitted_trials)
+    
     # Add actual photometry data to the NWB
     print("Adding photometry signals to NWB ...")