+- Python Forum (https://python-forum.io)
+-- Forum: Python Coding (https://python-forum.io/forum-7.html)
+--- Forum: Data Science (https://python-forum.io/forum-44.html)
+--- Thread: Fiber Photometry analysis (/thread-18303.html)
Fiber Photometry analysis - Gius_ - May-12-2019
Hi all,
I have a problem extracting events from a fiber photometry file using python.
I have 2 different files: the first one contains raw data (col 1: time milliseconds, col 2: GCaMP raw data) and the second file contains events of opto stimulation (col 1: time milliseconds and col 2 a series of 0 and 1 where length of 1 corresponding to the opto stimulation and length of zeros corresponding to no-stim). I'm able to extract the raw data, correct for bleaching and plot the Calcium signaling from the first file. However, I'm having troubles extracting the data around the opto stimulation from the second file and plot the signal around (+/- 5 seconds) the stimulation event.
here is the code I'm using for data extraction and plot:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import medfilt
from scipy.optimize import curve_fit
import seaborn as sns
from scipy.signal import medfilt, butter, filtfilt
data = pd.read_csv('GG1_fip_20Hz2019-05-11T11_38_00.csv', header=None, delimiter=',')
stim_data = pd.read_csv('GG1_20Hz_stim2019-05-11T11_37_58.csv', header=None, delimiter=',')
# extract values from FIP file
GCaMP = np.array(data[2])
row_time_fip = np.array(data[0])
time_fip = np.array(data[0] / 1000) # millisecond
normalized_time_sec = np.array(time_fip - min(time_fip))
normalized_time_min = np.array(normalized_time_sec / 60)
# extract values from stim file
stim_all = np.array(stim_data)
stim_onset = np.array(stim_data[1])
time_event = np.array(stim_data[0]/1000)
time_event_raw = np.array(stim_data[0])
normalized_time_event_sec = np.array(time_event - min(time_event))
normalized_time_event_min = np.array(normalized_time_sec / 60)
# denoise GCaMP signal (remove artifacts --> kernel should be odd number
GCaMP_denoised = (medfilt(GCaMP, kernel_size=51))
# method 1 bleaching correction = polyfit
coefs_GCaMP = np.polyfit(time_fip, GCaMP_denoised, deg=4)
GCaMP_polyfit = np.polyval(coefs_GCaMP, time_fip)
plt.subplot(2,1,1)
plt.plot(normalized_time_min, GCaMP_denoised, 'g', label='GCaMP')
plt.plot(normalized_time_min, GCaMP_polyfit,'k', linewidth=1.5)
plt.subplot(2,1,2)
GCaMP_ps = GCaMP_denoised - GCaMP_polyfit
plt.plot(normalized_time_min, GCaMP_ps, 'g', label='GCaMP')
plt.subplots_adjust(top=0.95, bottom=0.1, left=0.12, right=0.95, hspace=0.6, wspace=0.35)
plt.xlabel('Time (min)')
# method 2 bleaching correction = highpass slow frequency
b,a = butter(2, 0.001, btype='high')
GCaMP_highpass = filtfilt(b,a, GCaMP_denoised, padtype='even')
plt.plot(normalized_time_min, GCaMP_highpass)
plt.show()raw data: https://drive.google.com/open?id=1oIkV1H3AVk25SwcCxr4D8Oxnp2QM9OMb
https://drive.google.com/open?id=1UkxSaum4xv5ITXXJiE_lzhWL23wHepHt