From 146370dcede2f65c7c4c463e4dc720a7900ab855 Mon Sep 17 00:00:00 2001 From: dhruvm9 Date: Mon, 18 Sep 2023 16:37:44 -0400 Subject: [PATCH] makeRing --- Analysis/makeRing.py | 209 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 209 insertions(+) create mode 100644 Analysis/makeRing.py diff --git a/Analysis/makeRing.py b/Analysis/makeRing.py new file mode 100644 index 0000000..9513255 --- /dev/null +++ b/Analysis/makeRing.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon May 8 09:21:03 2023 + +@author: Dhruv +""" + +import numpy as np +import pandas as pd +import scipy.io +import pynapple as nap +import os, sys +import matplotlib.pyplot as plt +from matplotlib.colors import hsv_to_rgb +import scipy.signal +from sklearn.manifold import Isomap + +#%% + +# data_directory = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Data/AdrianPoSub/###AllPoSub' +data_directory = '/media/adrien/LaCie/PoSub-UPState/Data/###AllPoSub' +# datasets = np.genfromtxt(os.path.join(data_directory,'dataset_Hor_DM.list'), delimiter = '\n', dtype = str, comments = '#') +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_test.list'), delimiter = '\n', dtype = str, comments = '#') + +rwpath = '/media/adrien/LaCie/PoSub-UPState/Project/Data' + +alltruex = [] +alltruey = [] + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + rawpath = os.path.join(rwpath,s) + + data = nap.load_session(rawpath, 'neurosuite') + data.load_neurosuite_xml(rawpath) + spikes = data.spikes + epochs = data.epochs + + # ############################################################################################### + # # LOAD MAT FILES + # ############################################################################################### + filepath = os.path.join(path, 'Analysis') + listdir = os.listdir(filepath) + file = [f for f in listdir if 'BehavEpochs' in f] + behepochs = scipy.io.loadmat(os.path.join(filepath,file[0])) + + file = [f for f in listdir if 'CellDepth' in f] + celldepth = scipy.io.loadmat(os.path.join(filepath,file[0])) + depth = celldepth['cellDep'] + + file = [f for f in listdir if 'MeanFR' in f] + mfr = scipy.io.loadmat(os.path.join(filepath,file[0])) + r_wake = mfr['rateS'] + + file = [f for f in listdir if 'Layers' in f] + lyr = scipy.io.loadmat(os.path.join(filepath,file[0])) + layer = lyr['l'] + + file = [f for f in listdir if 'Velocity' in f] + vl = scipy.io.loadmat(os.path.join(filepath,file[0]), simplify_cells = True) + vel = nap.Tsd(t = vl['vel']['t'], d = vl['vel']['data'] ) + + file = [f for f in listdir if 'CellTypes' in f] + celltype = scipy.io.loadmat(os.path.join(filepath,file[0])) + + pyr = [] + interneuron = [] + hd = [] + + for i in range(len(spikes)): + if celltype['ex'][i] == 1 and celltype['gd'][i] == 1: + pyr.append(i) + + for i in range(len(spikes)): + if celltype['fs'][i] == 1 and celltype['gd'][i] == 1: + interneuron.append(i) + + for i in range(len(spikes)): + if celltype['hd'][i] == 1 and celltype['gd'][i] == 1: + hd.append(i) + + +#%% + + # ############################################################################################### + # # LOAD UP AND DOWN STATE, NEW SWS AND NEW WAKE EPOCHS + # ############################################################################################### + + + file = os.path.join(rawpath, name +'.evt.py.dow') + if os.path.exists(file): + tmp = np.genfromtxt(file)[:,0] + tmp = tmp.reshape(len(tmp)//2,2)/1000 + down_ep = nap.IntervalSet(start = tmp[:,0], end = tmp[:,1], time_units = 's') + + file = os.path.join(rawpath, name +'.evt.py.upp') + if os.path.exists(file): + tmp = np.genfromtxt(file)[:,0] + tmp = tmp.reshape(len(tmp)//2,2)/1000 + up_ep = nap.IntervalSet(start = tmp[:,0], end = tmp[:,1], time_units = 's') + + file = os.path.join(rawpath, name +'.DM.new_sws.evt') + if os.path.exists(file): + tmp = np.genfromtxt(file)[:,0] + tmp = tmp.reshape(len(tmp)//2,2)/1000 + new_sws_ep = nap.IntervalSet(start = tmp[:,0], end = tmp[:,1], time_units = 's') + + file = os.path.join(rawpath, name +'.DM.new_wake.evt') + if os.path.exists(file): + tmp = np.genfromtxt(file)[:,0] + tmp = tmp.reshape(len(tmp)//2,2)/1000 + new_wake_ep = nap.IntervalSet(start = tmp[:,0], end = tmp[:,1], time_units = 's') + +#%% + + filepath = os.path.join(path, 'Analysis') + data = pd.read_csv(filepath + '/Tracking_data.csv', header = None) + position = pd.DataFrame(index = data[0].values, data = data[[1,2,3]].values, columns=['x', 'y', 'ang']) + position = position.loc[~position.index.duplicated(keep='first')] + position['ang'] = position['ang'] *(np.pi/180) #convert degrees to radian + position['ang'] = (position['ang'] + 2*np.pi) % (2*np.pi) #convert [-pi, pi] to [0, 2pi] + position = nap.TsdFrame(position) + position = position.restrict(epochs['wake']) + + angle = position['ang'].loc[vel.index.values] + +#%% Wake binning + + wake_bin = 0.2 #0.4 #400ms binwidth + + ep = epochs['wake'] + + dx = position['x'].bin_average(wake_bin,ep) + dy = position['y'].bin_average(wake_bin,ep) + ang = angle.bin_average(wake_bin, ep) + v = vel.bin_average(wake_bin,ep) + + v = v.threshold(2) + + ang = ang.loc[v.index.values] + + wake_count = spikes[hd].count(wake_bin, ep) + wake_count = wake_count.loc[v.index.values] + + wake_count = wake_count.as_dataframe() + wake_rate = np.sqrt(wake_count/wake_bin) + wake_rate = wake_rate.rolling(window = 50, win_type = 'gaussian', center = True, min_periods = 1, axis = 0).mean(std = 3) + + ang = ang.dropna() + ang = ang.rename('ang') + + wake_rate = wake_rate.loc[ang.index.values] + + wake_rate = pd.concat([wake_rate, pd.DataFrame(ang)], axis = 1) + wake_rate = wake_rate.sample(frac = 0.5).sort_index() + +#%% Sleep binning + + sleep_dt = 0.01 #0.015 #0.025 #0.015 + sleep_binwidth = 0.05 #0.1 #0.03 + + tokeep = np.where((up_ep['end'] - up_ep['start']) > 0.5) + + du = nap.IntervalSet(start = up_ep.iloc[tokeep]['start'] - 0.25, end = up_ep.iloc[tokeep]['start'] + 0.25) + longdu = nap.IntervalSet(start = up_ep.iloc[tokeep]['start'] - 0.25, end = up_ep.iloc[tokeep]['start'] + 0.5) + + num_overlapping_bins = int(sleep_binwidth/sleep_dt) + + sleep_count = spikes[hd].count(sleep_dt,longdu) + sleep_count = sleep_count.as_dataframe() + sleep_rate = np.sqrt(sleep_count/sleep_dt) + sleep_rate = sleep_rate.rolling(window = num_overlapping_bins, win_type = 'gaussian', center = True, min_periods = 1, axis = 0).mean(std = 3) + + fit_rate = nap.TsdFrame(sleep_rate).restrict(du) + + rate = np.vstack([wake_rate.loc[:, wake_rate.columns != 'ang'].values, fit_rate.as_dataframe().sample(frac = 0.01).values]) #Take 1000 random values + + fit = Isomap(n_components = 2, n_neighbors = 200).fit(rate) + p_wake = fit.transform(wake_rate.loc[:, wake_rate.columns != 'ang']) + p_sleep = fit.transform(sleep_rate) + + + H = wake_rate['ang'].values/(2*np.pi) + HSV = np.vstack((H, np.ones_like(H), np.ones_like(H))).T + RGB = hsv_to_rgb(HSV) + + truex = p_sleep[np.where(sleep_rate.sum(axis=1)==0)][0][0] + truey = p_sleep[np.where(sleep_rate.sum(axis=1)==0)][0][1] + + alltruex.append(truex) + alltruey.append(truey) + + p_wake = p_wake - [truex, truey] + p_sleep = p_sleep - [truex, truey] + + projection = nap.TsdFrame(t = sleep_rate.index.values, d = p_sleep, columns = ['x', 'y']) + projection.to_pickle(rawpath + '/' + s + '_projection.pkl') + + np.save(rawpath + '/' + s + '_pwake.npy',p_wake) + np.save(rawpath + '/' + s + '_RGB.npy', RGB) + +#%% + +np.save(rwpath + '/alltruex.npy' , alltruex) +np.save(rwpath + '/alltruey.npy' , alltruey) +