Low-dimensional population dynamics in the brainstem gate REM sleep.

Code ↔ Paper

17 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 17 matches

[1] § Methods › Identification of laser-modulated units ↔ PySleep/sleepy.py, lines 2946–3001 · score 0.71 · interval preceding laser, laser pulse, laser stimulation, artifacts, train, onset
[2] § Methods › Identification of laser-modulated units ↔ PySpike/spyke.py, lines 722–771 · score 0.67 · spontaneous spikes, laser pulse, spike waveforms, correlation, mice
[3] § Methods › Dimensionality reduction using PCA ↔ neuropyx.py, lines 1610–1681 · score 0.66 · excluded long wake, wake episodes, longer, PCA, Firing rates, matrix
[4] § Results › Stereotypic NREM→REM transitions in state space ↔ neuropyx.py, lines 5096–5160 · score 0.64 · REM offset, state space, REM onset, ellipse, subspaces, trajectory
[5] § Methods › Definition of subspaces ↔ neuropyx.py, lines 5546–5677 · score 0.63 · ellipses capturing, state space, eigenvalue, subspaces, eigenvector, angles
[6] § Methods › Cross-correlation between single units, σ power and PC2 ↔ Photometry/pyphi.py, lines 1332–1481 · score 0.63 · temporal resolution, cross correlation, EEG spectrogram, lags, consecutive, overlap
[7] § Methods › Cross-correlation between single units, σ power and PC2 ↔ PySpike/spyke.py, lines 5273–5328 · score 0.61 · temporal resolution, cross correlation, EEG spectrogram, consecutive, overlap, spike
[8] § Methods › Neuron activity during infraslow cycles and inter-REM ↔ Photometry/pyphi.py, lines 3580–3688 · score 0.58 · fast Fourier transforms, band, microarousals, windows, signal, box
[9] § Methods › Spike sorting ↔ neuropyx.py, lines 750–825 · score 0.57 · Gaussian kernel, standard deviation, encoded, firing rates, smoothed, binned
[10] § Results › Low-dimensional population activity in midbrain and pons during sleep ↔ Photometry/pyphi.py, lines 1587–1645 · score 0.56 · confidence intervals, Bonferroni correction, EMG amplitude, EEG spectrogram, denote, error
[11] § Methods › Cross-correlation between single units, σ power and PC2 ↔ neuropyx.py, lines 750–825 · score 0.55 · Gaussian kernel, standard deviation, firing rates, episode, smoothed, bins
[12] § Methods › Habituation to head fixation and polysomnographic recordings ↔ Photometry/pyphi.py, lines 3580–3688 · score 0.55 · fast Fourier transforms, Brain states, REM sleep, window, signals, scored
[13] § Methods › Functional connectivity analysis ↔ PySpike/spyke.py, lines 5273–5328 · score 0.53 · negative peak, cross correlated, resolution, consecutive, Spike, windows
[14] § Methods › Analysis of coupling frequency ↔ Photometry/pyphi.py, lines 1587–1645 · score 0.53 · mouse identity, Bonferroni corrected, CI
[15] § Methods › Electrode tract reconstruction ↔ basic_analysis_howto.ipynb, lines 69–106 · score 0.53 · MB, PRNc, PRNr, MRN, CS, RPO
[16] § Methods › Spike sorting ↔ PySpike/spyke.py, lines 1657–1770 · score 0.51 · spike trains binned, Clusters, encoded, waveform, downsampled
[17] § Methods › Electrophysiological recordings ↔ PySleep/sleepy.py, lines 2946–3001 · score 0.51 · pulse trains, laser stimulation, REM sleep, interval, wake, mice

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 5,681 lines · 195 KB · MIT · 5 matches

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
neuropyx.py
"""
import sys
sys.path.append('/Users/tortugar/My Drive/Penn/Programming/PySleep')
import sleepy
import numpy as np
import scipy
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import os
import csv
import pingouin as pg
import seaborn as sns
import scipy.io as so
import re
import scipy.stats as stats
from scipy import linalg
import matplotlib as mpl
from sklearn.decomposition import PCA
import math
import h5py
import shutil
# debugger
import pdb
def brstate_class(np_path, sleep_path, sleep_rec, mouse, tend=-1, tstart=0, pzscore=True,
class_mode='', pearson=True, pnorm_spec=True, single_mice=True,
ma_thr=10, ma_rem_exception=False, box_filt=[],
pplot=True, config_file='mouse_config.txt'):
"""
calculate average firing rate during each brain state and then
perform statistics for units to classify them into REM-max, Wake-max, or NREM-max.
For each ROI anova is performed, followed by Tukey-test
:param np_path: folder where firing rates are located
:param sleep_path: folder where EEG data and sleep annotation is located
:param sleep_rec: name of sleep recording
Note: easiest way to get np_path, sleep_path, and sleep_rec:
paths = neuropyx.load_config(config_file)[mouse]
ppath, name = os.path.split(paths['SL_PATH'])
np_path = paths['NP_PATH']
:param mouse: mouse name
:param pzscore: if True, z-score DF/F traces
:param class_mode: class_mode == 'basic': classify ROIs into
REM-max, Wake-max and NREM-max ROIs
class_mode == 'rem': further separate REM-max ROIs
into REM > Wake > NREM (R>W>N) and REM > NREM > Wake (R>N>W) ROIs
REM-max neurons where Wake and NREM is not significantly different
are classified as 'R>N=W'
Neurons that are signficantly modulated by the brain state
but not part of any of these different classes are labeled Z
Neurons that are not significantly modulated by the brain state
are labeled X
NOTE: A given unit can only be part of one subclass.
R-Off comes before W-max.
:param single_mice: boolean, if True use separate colors for single mice in
summary plots
:return df_class: pd.DataFrame
with columns
['ID', 'R', 'W', 'N', 'F-anova', 'P-anova', 'P-tukey', 'Type',
'Depth', 'Quality', 'mouse', 'brain_region']
:rem df_pearson: pd.DataFrame
with columns
['unit', 'r', 'p', 'sig', 'state','depth','Type','Quality']
"""
units, cell_info, M, kcut = load_mouse(mouse, config_file)
# import file that has information on depth
with open(os.path.join(np_path,'cluster_info.TSV')) as inf:
reader = csv.reader(inf, delimiter="\t")
cell_info=list(reader)
cell_list=[]
for x in cell_info[1:]:
cell_list.append([x[0],x[3],x[6]])
# import sleep annotation file and also do MA functions if chosen
sr = sleepy.get_snr(sleep_path,sleep_rec)
nbin = int(np.round(sr)*2.5)
sdt = nbin * (1.0/sr)
# flatten out MAs #########################################################
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*sdt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
###########################################################################
# cut out kcuts: ###############
#print('Applying kcut')
tidx = kcut_idx(M, units, kcut)
if tidx[-1] >= units.shape[0]:
tidx = tidx[0:-1]
M = M[tidx]
units = units.iloc[tidx,:]
################################
istart = int(np.round(tstart/sdt))
if tend == -1:
iend = M.shape[0]
else:
iend = int(np.round(tend/sdt))
if iend >= len(units):
M = M[istart:len(units)]
else:
M = M[istart:iend]
# calculate sigma power
if pearson:
band=[10,15]
state_map = {1:'REM', 2:'Wake', 3:'NREM'}
ddir = os.path.join(sleep_path, sleep_rec)
P = so.loadmat(os.path.join(ddir, 'sp_%s.mat' % sleep_rec), squeeze_me=True)
SP = P['SP']
freq = P['freq']
ifreq = np.where((freq >= band[0]) & (freq <= band[1]))[0]
df = freq[1] - freq[0]
if len(box_filt) > 0:
filt = np.ones(box_filt)
filt = np.divide(filt, filt.sum())
SP = scipy.signal.convolve2d(SP, filt, boundary='symm', mode='same')
if pnorm_spec:
sp_mean = SP.mean(axis=1)
SP = np.divide(SP, np.tile(sp_mean, (SP.shape[1], 1)).T)
pow_band = SP[ifreq, :].mean(axis=0)
else:
pow_band = SP[ifreq,:].sum(axis=0)*df
pow_band = pow_band[istart:iend]
# make a nested dict that has each of the units and inside the fr values for all the different sleep bins
units_stateval = {}
for unit in units:
units_stateval[unit] = {1:[], 2:[], 3:[],'depth':[]}
for unit in units:
depth_idx=units.columns.get_loc(unit)
if pzscore:
values = np.array((units[unit] - units[unit].mean())/units[unit].std(ddof=0))
else:
values=np.array(units[unit])
for state in [1,2,3]:
seq = np.where(M==state)[0]
units_stateval[unit][state]=values[seq].tolist()
units_stateval[unit]['depth']=cell_list[depth_idx][2]
columns = ['ID', 'R', 'W', 'N', 'F-anova', 'P-anova', 'P-tukey', 'Type','Depth','Quality']
data = []
data_p = []
for unit in units_stateval:
stateval = units_stateval[unit]
val = np.concatenate([stateval[1], stateval[2], stateval[3]],axis=0)
state = ['R']*len(stateval[1]) + ['W']*len(stateval[2]) + ['N']*len(stateval[3])
depth=float(units_stateval[unit]['depth'])
if 'good' in unit:
unit_quality='good'
elif 'noise' in unit:
continue
else:
unit_quality='mua'
continue
d = {'state':state, 'val':val}
df = pd.DataFrame(d)
res = pg.anova(data=df, dv='val', between='state')
try:
res2 = pg.pairwise_tukey(data=df, dv='val', between='state')
except:
print("Unit %s did not have enough data points for tukey test" % unit)
print(df)
def _get_mean(s):
return df[df['state']==s]['val'].mean()
rmean = _get_mean('R')
wmean = _get_mean('W')
nmean = _get_mean('N')
if class_mode == 'basic':
roi_type = 'X'
# REM-max
if (rmean > wmean) and (rmean > nmean):
cond1 = res2[(res2['A'] == 'N') & (res2['B'] == 'R')]
cond2 = res2[(res2['A'] == 'R') & (res2['B'] == 'W')]
if cond1['p-tukey'].iloc[0] < 0.05 and cond2['p-tukey'].iloc[0] < 0.05 and res['p-unc'].iloc[0] < 0.05:
roi_type = 'R-max'
#REM-Off (R-Off)
elif (rmean < wmean) and (rmean < nmean):
cond1 = res2[(res2['A'] == 'R') & (res2['B'] == 'W')]
cond2 = res2[(res2['A'] == 'R') & (res2['B'] == 'N')]
if cond1['p-tukey'].iloc[0] < 0.05 and cond2['p-tukey'].iloc[0] < 0.05 and res['p-unc'].iloc[0] < 0.05:
roi_type = 'R-Off'
# W-max
elif (wmean > nmean) and (wmean > rmean):
cond1 = res2[(res2['A'] == 'N') & (res2['B'] == 'W')]
cond2 = res2[(res2['A'] == 'R') & (res2['B'] == 'W')]
if cond1['p-tukey'].iloc[0] < 0.05 and cond2['p-tukey'].iloc[0] < 0.05 and res['p-unc'].iloc[0] < 0.05:
roi_type = 'W-max'
# N-max
elif (nmean > wmean) and (nmean > rmean):
cond1 = res2[(res2['A'] == 'N') & (res2['B'] == 'W')]
cond2 = res2[(res2['A'] == 'N') & (res2['B'] == 'R')]
if cond1['p-tukey'].iloc[0] < 0.05 and cond2['p-tukey'].iloc[0] < 0.05 and res['p-unc'].iloc[0] < 0.05:
roi_type = 'N-max'
else:
roi_type = 'X'
tmp = [unit, rmean, wmean, nmean, res.F.iloc[0], res['p-unc'].iloc[0], res2['p-tukey'].iloc[0], roi_type,depth,unit_quality]
data.append(tmp)
# REM mode:
else:
roi_type = 'X'
if res['p-unc'].iloc[0] < 0.05:
p_nr = res2[(res2['A'] == 'N') & (res2['B'] == 'R')]['p-tukey'].iloc[0]
p_rw = res2[(res2['A'] == 'R') & (res2['B'] == 'W')]['p-tukey'].iloc[0]
p_nw = res2[(res2['A'] == 'N') & (res2['B'] == 'W')]['p-tukey'].iloc[0]
# R>W>N
if (rmean > wmean) and (rmean > nmean) and (wmean > nmean) and p_nr < 0.05 and p_rw<0.05 and p_nw < 0.05:
roi_type = 'R>W>N'
# R>N>W
elif (rmean > wmean) and (rmean > nmean) and (nmean > wmean) and p_nr < 0.05 and p_rw<0.05 and p_nw < 0.05:
roi_type = 'R>N>W'
# NEW:R>N=W #####################################################
# The remaining REM-max units: R>N and R>W, but N and W are not significantly different
# I'm calling these units R>N=W
elif (rmean > wmean) and (rmean > nmean) and p_nr < 0.05 and p_rw<0.05 and p_nw >= 0.05:
roi_type = 'R>N=W'
# END[NEW:R>N=W] ###############################################
# Rem-off
elif (rmean < wmean) and (rmean < nmean) and p_nr < 0.05 and p_rw<0.05:
roi_type = 'R-Off'
# W-max
elif (wmean > nmean) and (wmean > rmean) and p_nw < 0.05 and p_rw < 0.05:
roi_type = 'W-max'
# N-max
elif (nmean > wmean) and (nmean > rmean) and p_nw < 0.05 and p_nr < 0.05:
roi_type = 'N-max'
# NEW:Z ######################################################
# significantly modulated by brainstate (according to ANOVA),
# but not part of any of these subclasses
else:
roi_type = 'Z'
else:
roi_type = 'X'
tmp = [unit, rmean, wmean, nmean, res.F.iloc[0], res['p-unc'].iloc[0], res2['p-tukey'].iloc[0], roi_type,depth,unit_quality]
data.append(tmp)
if pearson:
for s in [1,2,3]:
idx = np.where(M==s)[0]
r,p = scipy.stats.pearsonr(np.array(units[unit])[idx], pow_band[idx])
if p < 0.05:
sig = 'yes'
else:
sig = 'no'
pearson_temp=[unit, r, p, sig, state_map[s],depth,roi_type,unit_quality]
data_p.append(pearson_temp)
df_class = pd.DataFrame(data, columns=columns)
# NEW 3/11/26: replaced 'unit' with 'ID'
df_pearson = pd.DataFrame(data_p,columns=['ID', 'r', 'p', 'sig', 'state','depth','Type','Quality'])
if pplot:
# plotting for unit type
mice = [mouse]
j = 0
mdict = {}
for m in mice:
mdict[m] = j
j+=1
clrs = sns.color_palette("husl", len(mice))
types = df_class['Type'].unique()
types = [i for i in types if not (i=='X')]
types.sort()
j = 1
plt.figure()
for typ in types:
mouse_shown = {m:0 for m in mice}
plt.subplot(int('1%d%d' % (len(types), j)))
df = df_class[df_class['Type']==typ][['R', 'N', 'W']]
sns.barplot(data=df[['R', 'N', 'W']], color='gray')
for index, row in df.iterrows():
if single_mice:
if mouse_shown[m] > 0:
plt.plot(['R', 'N', 'W'], row[['R', 'N', 'W']], color=clrs[mdict[m]])
else:
plt.plot(['R', 'N', 'W'], row[['R', 'N', 'W']], color=clrs[mdict[m]], label=m)
mouse_shown[m] += 1
else:
plt.plot(['R', 'N', 'W'], row[['R', 'N', 'W']], color='black')
sns.despine()
plt.title(typ)
plt.legend()
if j == 1:
if not pzscore:
plt.ylabel('DF/F (%)')
else:
plt.ylabel('Firing Rate (z-scored)')
j += 1
# plot swarm plot with depth vs type
df_class_good=df_class.loc[df_class['Quality']=='good']
plt.figure()
plt.title('Unit type sorted by depth ')
sns.swarmplot(data=df_class_good, x='Type', y='Depth',palette="husl")
if pearson:
#plotting pearson plot
test_df=df_pearson.loc[df_pearson['sig']=='yes']
test_df=test_df.loc[test_df['Quality']=='good']
plt.figure()
sns.swarmplot(data=test_df, x='Type', y='r', hue='state',palette="husl")
plt.title('correllation of firing rate to'+ ' '+ 'eeg band:'+ str(band))
plt.figure()
sns.swarmplot(data=test_df, x='Type', y='depth', hue='state',palette="husl")
plt.title('correllation of firing rate to'+ ' '+ 'eeg band:'+ str(band))
return df_class, df_pearson
def load_config(config):
"""
Load a config file specifying the file locations for each mouse of
the sleep recording and the neuropixel recording;
Syntax:
MOUSE: Mouse_name
SL_PATH: Sleep_recording_folder
NP_PATH: Neuropixels_data_folder
MOUSE: Mouse_name
SL_PATH: Sleep_recording_folder
NP_PATH: Neuropixels_data_folder
KCUT: t0-t1;t2-t3
NO_REGION: A-B-C
EXCLUDE: unit_ID1,unit_ID2,...,unit_ID2
NOTE:
Different mice are separated by new lines
After a newline the first entry must be MOUSE.
KCUT: is optional and allows to specify a time frame to be discarded in the
recording.
NO_REGION: brain_regions to exclude
Parameters
----------
config : str
text file (including path)
Returns
-------
recordings : dict
dictionary: mouse_ID --> SL_PATH, NP_PATH
"""
fid = open(config, 'r')
lines = fid.readlines()
mouse = 'X'
recordings = dict()
newline = True
for l in lines:
if re.match(r'^\s*#', l):
continue
if len(l)==0 or re.match(r'^\s*$', l):
newline = True
continue
a = re.split(r'\s+', l)
# cut away the ':' or ' :'
field = re.split(r'\s*:', a[0])[0]
value = a[1]
if not newline:
recordings[mouse][field] = value
else: # newline == True
mouse = a[1]
if not mouse in recordings:
recordings[mouse] = {}
newline = False
for m in recordings:
if 'KCUT' in recordings[m]:
a = recordings[m]['KCUT']
a = re.split(';', a)
kcuts = []
for b in a:
c = re.split('-', b)
c = [s.strip() for s in c]
k1 = float(c[0])
if re.match(r'[\d\.]+', c[1]):
k2 = float(c[1])
else:
k2 = c[1]
kcuts.append( [k1, k2] )
recordings[m]['KCUT'] = kcuts
for m in recordings:
if 'EXCLUDE' in recordings[m]:
a = recordings[m]['EXCLUDE']
k = re.split(r',\s*', a)
recordings[m]['EXCLUDE'] = k
return recordings
def load_mouse(mouse_id, config_file):
"""
Load neuropixels recordings as described in mouse_config.txt
Parameters
----------
mouse_id : TYPE
DESCRIPTION.
config_file : TYPE
DESCRIPTION.
Returns
-------
units : pd.DataFrame
unit DataFrame: The columns correspond to the firing rates;
the column name is the ID of the unit, the rows correspond to single time points.
cell_info : pd.DataFrame
DESCRIPTION.
"""
dt = 2.5
recs = load_config(config_file)
sl_path = recs[mouse_id]['SL_PATH']
(sl_path, sl_name) = os.path.split(sl_path)
np_path = recs[mouse_id]['NP_PATH']
# Get sleep annotation:
M = sleepy.load_stateidx(sl_path, sl_name)[0]
# load kcuts, i.e. regions at beginning or end to discard from the recording --
if 'KCUT' in recs[mouse_id]:
kcut = recs[mouse_id]['KCUT']
for k in kcut:
if k[1] == '$':
k[1] = len(M)*dt
else:
kcut = ()
traind_file = ''
if os.path.isfile(os.path.join(np_path, 'traind.csv')):
traind_file = 'traind.csv'
elif os.path.isfile(os.path.join(np_path, 'spike_train.csv')):
traind_file = 'spike_train.csv'
else:
traind_file = 'traind_lfp.csv'
units = pd.read_csv(os.path.join(np_path, traind_file))
if os.path.isfile(os.path.join(np_path,'channel_locations.json')):
regions=pd.read_json(os.path.join(np_path,'channel_locations.json')).T
regions=regions.iloc[0:-1]
regions['ch']=regions.index.str.split('_').str[-1].astype('int64')
cell_info = pd.read_csv(os.path.join(np_path,'cluster_info.TSV'),delimiter="\t")
cell_info['group']=cell_info['group'].fillna(cell_info['KSLabel'])
cl_id = ''
if cell_info.columns.isin(['cluster_id']).sum():
cl_id = 'cluster_id'
else:
cl_id = 'id'
cell_info['ID'] = cell_info[cl_id].astype(str) +'_'+cell_info['group'].astype(str)
if os.path.isfile(os.path.join(np_path,'channel_locations.json')):
cell_info=cell_info.merge(regions,on='ch')
# NEW 07/15/23
if len(M) > units.shape[0]:
M = M[0:units.shape[0]]
return units, cell_info, M, kcut
def exclude_units(units, mouse, config_file):
"""
Exclude units listed in mouse_config.txt under 'EXCLUDE:'
Parameters
----------
units : pd.pandas
each colums hold the firing rate of a unit.
mouse : string
the mouse name.
config_file : string
string of the name of the config file as loaded by &load_config().
Returns
-------
None.
"""
paths = load_config(config_file)[mouse]
if 'EXCLUDE' in paths:
ex_units = paths['EXCLUDE']
units.drop(columns=ex_units, inplace=True)
def fr_corr_state(units, M, idx1=[], idx2=[], win=60, state=3, ma_thr=10, mode='cross', pzscore=True, pplot=True, dt=2.5):
"""
Perform cross-correlation between firing rates for a given brain state.
Calculate the correlation for each pair of the provided neurons unitIDs (in idx1 and idx2)
Parameters
----------
units : pd.DataFrame
Each column corresponds to one unit. The column name is the unitID.
So to get all unit ID, get all column names
M : np.array
hypnogram.
idx1 : list, optional
List of unitIDs for neuron1. If empty, use all neurons (unitIDs)
The default is [].
idx2 : list, optional
List of unitIDs for neuron1. If empty, use all neurons (unitIDs)
The default is [].
win : float, optional
DESCRIPTION. The default is 120.
state : TYPE, optional
DESCRIPTION. The default is 3.
ma_thr : TYPE, optional
DESCRIPTION. The default is 10.
mode : TYPE, optional
DESCRIPTION. The default is 'cross'.
pzscore : TYPE, optional
DESCRIPTION. The default is True.
pplot : TYPE, optional
DESCRIPTION. The default is True.
dt : TYPE, optional
DESCRIPTION. The default is 2.5.
Returns
-------
df : TYPE
DESCRIPTION.
"""
# If no IDs are provided for unit1, use all units as "unit1"
if len(idx1) == 0:
idx1 = units.columns.unique()
# If no IDs are provided for unit2, use all units as "unit2"
if len(idx2) == 0:
idx2 = units.columns.unique()
# get all firing rates and cast them to np.array
arr1 = units[idx1]
arr2 = units[idx2]
fr1 = np.array(arr1)
fr2 = np.array(arr2)
n1 = fr1.shape[1]
n2 = fr2.shape[1]
# z-scoring
if pzscore:
for i in range(n1):
fr1[:,i] = (fr1[:,i]-fr1[:,i].mean()) / fr1[:,i].std()
for i in range(n2):
fr2[:,i] = (fr2[:,i]-fr2[:,i].mean()) / fr2[:,i].std()
iwin = int(win/dt)
t = np.arange(-iwin, iwin+1) * dt
seq = sleepy.get_sequences(np.where(M==2)[0])
if ma_thr > 0:
seq = sleepy.get_sequences(np.where(M == 2)[0])
for s in seq:
if len(s) * dt < ma_thr:
M[s] = 3
seq = sleepy.get_sequences(np.where(M==state)[0])
seq = [s for s in seq if len(s)*dt >= 2*win]
data = []
for i in range(n1):
for j in range(n2):
CC = []
for s in seq:
m = len(s)
fr1_cut = fr1[s,i]
fr2_cut = fr2[s,j]
fr1_cut = fr1_cut - fr1_cut.mean()
fr2_cut = fr2_cut - fr2_cut.mean()
norm = np.nanstd(fr1_cut) * np.nanstd(fr2_cut)
# for used normalization, see: https://en.wikipedia.org/wiki/Cross-correlation
if norm > 0:
xx = (1/m) * scipy.signal.correlate(fr1_cut, fr2_cut) / norm
ii = np.arange(len(xx) / 2 - iwin, len(xx) / 2 + iwin + 1)
ii = [int(i) for i in ii]
ii = np.concatenate((np.arange(m-iwin-1, m), np.arange(m, m+iwin, dtype='int')))
# note: point ii[iwin] is the "0", so xx[ii[iwin]] corresponds to the 0-lag correlation point
CC.append(xx[ii])
if len(CC) > 0:
CC = np.array(CC).mean(axis=0)
#pdb.set_trace()
m = CC.shape[0]
un1 = list(idx1)[i]
un2 = list(idx2)[j]
label = r'%s~%s' % (un1, un2)
data += zip(t, CC, [un1]*m, [un2]*m, [label]*m)
df = pd.DataFrame(data=data, columns=['time', 'cc', 'unit1', 'unit2', 'label'])
return df
def sort_xcorr_byregion(type1, type2, corr_frame, unit_info):
"""
Go through all brain regions contained in DataFrame unit_info['brain_region'],
For each pair of brain regions take all neurons of type1 in area A and all neurons
of type2 in area B and determine the average cross-correlation for this type/area pair.
Parameters
----------
type1 : TYPE
type2 : TYPE
DESCRIPTION.
corr_frame : pd.DataFrame with columns ['time', 'cc', 'unit1', 'unit2', 'label']
DataFrame as returned by the function neuropyx.fr_corr_state()
The DataFrame contains for pairs of neurons (indicated by the IDs in columns 'unit1' and 'unit2')
the cross-correlation
unit_info : pd.DataFrame
Necessary columns: 'ID' (IDs of units), 'brain_region' (brain region of unit 'ID'),
'Type' (neuron subclass as computed by neuropyx.brstate_fr())
Returns
-------
None.
"""
regions = unit_info.brain_region.unique()
nregions = len(regions)
data = []
for i in range(nregions):
for j in range(0, nregions):
region1 = regions[i]
region2 = regions[j]
ids1 = unit_info[(unit_info['brain_region']==region1) & (unit_info['Type']==type1)]['ID']
ids2 = unit_info[(unit_info['brain_region']==region2) & (unit_info['Type']==type2)]['ID']
print(ids2)
if len(ids1) > 0 and len(ids2) > 0:
df = corr_frame[(corr_frame.unit1.isin(ids1)) & (corr_frame.unit2.isin(ids2))]
dfm = df.groupby(['time']).mean()
time = np.array(dfm.index)
cc = dfm['cc']
m = len(cc)
label = '%s~%s' % (region1, region2)
data += zip(time, cc, [region1]*m, [region2]*m, [label]*m)
df = pd.DataFrame(data=data, columns=['time', 'cc', 'region1', 'region2', 'label'])
return df
def fr_transitions(units, M, unit_info, transitions, pre, post, si_threshold, sj_threshold,
ma_thr=10, ma_rem_exception=False, sdt=2.5, pzscore=True, sf=0, ma_mode=False,
attributes=[], kcuts=[],
pspec=False, fmax=20, spe_filt=[], mouse='', config_file='mouse_config.txt'):
"""
Note: If you like to also calculate the average EEG specotrogram for each transition and unit,
set $psec=True, set the $mouse name and select the right $config_file
Parameters
----------
units : pd.DataFrame
Each columns, one unit.
M : np.array
hypnogram; 1 - REM, 2 - Wake, 3 - NREM.
unit_info : pd.DataFrame
Additional unit information, such as brain_region. Units are referenced
using the same IDs as in @units
transitions : list of tuples (2 element lists)
Specific the transitions to be calculated.
REM - 1, Wake - 2, NREM - 3, MA - 4;
So a NREM to REM transition is specified as [3,1]
pre : float
time before transition.
post : float
time after transition.
si_threshold : list with 3 floats
For example, if we're looking at a NREM to REM transition, for how long
should the mouse be in NREM before the actual transition. Specify for
NREM, Wake, REM the minimum preceding bout duration
sj_threshold : list with 3 floats
For example, if we're looking at a NREM to REM transition, for how long
should the mouse be in NREM before the actual transition. Specify for
NREM, Wake, REM the minimum following bout duration.
ma_thr : TYPE, float
Wake sequences <= $ma_thr s are interpreted as NREM (3). The default is 10.
ma_rem_exception : bool, optional
If True, then the MA rule does not apply for wake episodes directly following REM.
The default is False.
sdt : float, optional
Time binning for firing rates and hypnogram. The default is 2.5.
pzscore : bool, optional
If True, z-score firing rates. The default is True.
sf : float, optional
Standard deviation for Gaussian kernel to smooth firing rates. The default is 0.
ma_mode : bool, optional
If True, treat MAs as their own brain state.
attributes : list of strings, optional
Allows you to transfer columns in DataFrame @unit_info to the returned DataFrame @df. The default is [].
kcuts : list of tuples or lists with two elements.
Discard the time interval ranging from kcuts[i][0] to kcuts[i][1] seconds
pspect : bool
if True, also claculate EEG spectrogram
mouse : str,
If $pspec == True, needs to be set to an existing
mouse name
config_file : str
Name of mouse recording configuration file
fmax : float
Maximum frequency for EEG spectrogram
Returns
-------
df : pd.DataFrame
with columns ['ID', 'time', 'fr', 'trans'], i.e.
unit IDs, time point relative to transition, firing rate value, transition type
encoded as NR, NW etc.
mx_transspe : dict: Transition type --> np.array
The dict holds for each calculated transition (encoded as NR, NW, etc.)
and each unit (axis=2) the average spectrogram during that transition.
"""
dt = 2.5
states = {1:'R', 2:'W', 3:'N', 4:'M'}
# cut out kcuts: ###############
tidx = kcut_idx(M, units, kcuts)
M = M[tidx]
units = units.iloc[tidx,:]
################################
ipre = int(np.round(pre/sdt))
ipost = int(np.round(post/sdt))
m = ipre + ipost
t = np.arange(-ipre, ipost) * sdt
unitIDs = units.columns.unique()
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*sdt) <= ma_thr:
if ma_rem_exception:
if (s[0]>=1) and (M[s[0] - 1] != 1):
if ma_mode:
M[s] = 4
else:
M[s] = 3
else:
if ma_mode:
M[s] = 4
else:
M[s] = 3
# NEW: load spectrogram
# load spectrogram and normalize
if pspec:
path = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(path)
P = so.loadmat(os.path.join(ppath, name, 'sp_%s.mat' % name), squeeze_me=True)
SP = P['SP']
freq = P['freq']
ifreq = np.where(freq <= fmax)[0]
if len(spe_filt) > 0:
filt = np.ones(spe_filt)
filt = np.divide(filt, filt.sum())
SP = scipy.signal.convolve2d(SP, filt, boundary='symm', mode='same')
sp_mean = SP.mean(axis=1)
SP = np.divide(SP, np.tile(sp_mean, (SP.shape[1], 1)).T)
SP[:,tidx]
unit_transspe = {}
mx_transspe = {}
for (si,sj) in transitions:
# string label for type of transition:
sid = states[si] + states[sj]
unit_transspe[sid] = {r:[] for r in unitIDs}
mx_transspe[sid] = np.zeros((len(ifreq), len(t), len(unitIDs)))
##########################################################################
data = []
for unit in unitIDs:
unit_annotation = unit_info[unit_info.ID == unit]
if unit_annotation.shape[0] > 0:
attr = unit_annotation[attributes].values.tolist()[0]
else:
attr = ['X', 'X']
fr = np.array(units[unit])
if sf > 0:
fr = sleepy.smooth_data(fr, sf)
if pzscore:
fr = (fr-fr.mean()) / fr.std()
for (si,sj) in transitions:
# string label for type of transition:
sid = states[si] + states[sj]
seq = sleepy.get_sequences(np.where(M==si)[0])
for s in seq:
ti = s[-1]
# check if next state is sj; only then continue
if ti < len(M)-1 and M[ti+1] == sj:
# go into future
p = ti+1
while p<len(M)-1 and M[p] == sj:
p += 1
p -= 1
sj_idx = list(range(ti+1, p+1))
# so the indices of state si are seq
# the indices of state sj are sj_idx
if ipre <= ti < len(M)-ipost and len(s)*sdt >= si_threshold[si-1] and len(sj_idx)*sdt >= sj_threshold[sj-1]:
rem_ID = '%d' % (s[0])
act = fr[ti-ipre+1:ti+ipost+1]
# Note: ti+1 is the first time point of the "post" state
# i = 10, ipre = 2, ipost = 2
# 8,9,10
# np.arange(8,12) = 8,9,10,11,12
if pspec:
spe = SP[ifreq, ti-ipre+1:ti+ipost+1]
unit_transspe[sid][unit].append(spe)
dur_post = len(sj_idx)*dt
fr_post = np.mean(fr[sj_idx])
new_data = zip([unit]*m, t, act, [sid]*m, [rem_ID]*m, [dur_post]*m, [fr_post]*m)
new_data = [list(x) + attr for x in list(new_data)]
data += new_data
df = pd.DataFrame(data=data, columns=['ID', 'time', 'fr', 'trans', 'remID', 'dur_post', 'fr_post'] + attributes)
if pspec:
for (si,sj) in transitions:
for i,unit in enumerate(unitIDs):
# string label for type of transition:
sid = states[si] + states[sj]
tmp = np.array(unit_transspe[sid][unit])
mx_transspe[sid][:,:,i] = np.nanmean(tmp, axis=0)
if not pspec:
return df, []
else:
return df, mx_transspe
def fr_transitions_stats(df_trans, base_int, unit_avg=True, dt=2.5, time_mode='midpoint'):
"""
Parameters
----------
df_trans : TYPE
DESCRIPTION.
base_int : TYPE
DESCRIPTION.
unit_avg : TYPE, optional
DESCRIPTION. The default is True.
dt : TYPE, optional
DESCRIPTION. The default is 2.5.
time_mode : str, optional
options: 'midpoint' or 'endpoint'. The default is 'midpoint'.
Returns
-------
df_stats : TYPE
DESCRIPTION.
"""
# test if df has column ms_id
if not 'ms_id' in df_trans.columns:
mice = list(df_trans['mouse'])
ids = list(df_trans['ID'])
ms_ids = [m + '-' + i for m,i in zip(mice, ids)]
df_trans['ms_id'] = ms_ids
ids = list(df_trans.ms_id.unique())
dfm_trans = df_trans[['ms_id', 'time', 'fr', 'trans']].groupby(['ms_id', 'time', 'trans',]).mean().reset_index()
t = dfm_trans.time.unique()
# number of bins per time bin
ibin = int(base_int / dt)
pre = t[0]
post = t[-1]
nbin = int(np.floor((abs(pre)+post)/base_int))
trans_dict = {}
for tr in df_trans.trans.unique():
trans_mx = np.zeros((len(ids), len(t)))
for i,ID in enumerate(ids):
fr = dfm_trans.loc[(dfm_trans.ms_id == ID) & (dfm_trans.trans==tr), 'fr']
trans_mx[i,:] = fr
trans_dict[tr] = trans_mx
tinit = t[0]
# Statistics: When does activity becomes significantly different from baseline?
ibin = int(np.round(base_int / dt))
nbin = int(np.floor((abs(pre)+post)/base_int))
data = []
for tr in trans_dict:
trans = trans_dict[tr]
base = trans[:,0:ibin].mean(axis=1)
for i in range(1,nbin):
test_vals = trans[:, i * ibin:(i + 1) * ibin].mean(axis=1)
ttest_res = stats.ttest_rel(base, trans[:,i*ibin:(i+1)*ibin].mean(axis=1))
tval = ttest_res.statistic
pval = ttest_res.pvalue
dof = len(base) - 1
# Compute Cohen's d
diff = base - test_vals
cohend = diff.mean() / diff.std(ddof=1)
sig = 'no'
if pval < (0.05 / (nbin-1)):
sig = 'yes'
if time_mode == 'midpoint':
tpoint = i*(ibin*dt)+tinit + ibin*dt/2
else:
tpoint = i*(ibin*dt)+tinit + ibin*dt
tpoint = float('%.2f'%tpoint)
pval = pval * (nbin-1)
if pval > 1:
pval = 1
data.append([tpoint, pval, sig, tr, tval, dof, cohend])
df_stats = pd.DataFrame(data = data, columns = ['time', 'p-value', 'sig', 'trans', 'tval', 'dof', 'cohend'])
return df_stats
def pc_transitions(PC, M, transitions, pre, post, si_threshold, sj_threshold,
ma_thr=10, ma_rem_exception=False, sdt=2.5, ma_mode=False,
kcuts=[], allowed_idx=[], pzscore_pc=False):
"""
Calculate timecourse of PCs in population activity relative to brain state
transitions.
Parameters
----------
PC : np.array
Number of PCx x number of time bins
PCs; each row corresponds to one PC.
M : np.array
hynpogram.
transitions : TYPE
DESCRIPTION.
pre : TYPE
DESCRIPTION.
post : TYPE
DESCRIPTION.
si_threshold : TYPE
DESCRIPTION.
sj_threshold : TYPE
DESCRIPTION.
ma_thr : TYPE, optional
DESCRIPTION. The default is 10.
ma_rem_exception : bool, optional
If True, don't touch wake following REM sleep. The default is False.
sdt : float, optional
time bin duration in seconds of one brain state. The default is 2.5.
ma_mode : bool, optional
If True, then specifically analyze transitions from and to MAs.
Note that when considering transitions from and to NREM, MAs are
considered as NREM sleep.
kcuts : TYPE, optional
DESCRIPTION. The default is [].
allowed_idx : TYPE, optional
DESCRIPTION. The default is [].
pzscore_pc : bool
If True, z-score PCs across entire recording
Returns
-------
df : pd.DataFrame
with columns ['event', 'pc', 'time', 'fr', 'trans'].
"""
states = {1:'R', 2:'W', 3:'N', 4:'M'}
# cut out kcuts: ###############
tidx = kcut_idx(M, PC, kcuts)
M = M[tidx]
################################
ipre = int(np.round(pre/sdt))
ipost = int(np.round(post/sdt))
m = ipre + ipost
t = np.arange(-ipre, ipost) * sdt
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*sdt) <= ma_thr:
if ma_rem_exception:
if (s[0]>=1) and (M[s[0] - 1] != 1):
if ma_mode:
M[s] = 4
else:
M[s] = 3
else:
if ma_mode:
M[s] = 4
else:
M[s] = 3
Mrepr = M.copy()
M[M==4] = 3
if not ma_mode:
#just forget about MAs:
Mrepr = M
if len(allowed_idx) == 0:
allowed_idx = range(0, len(M))
# zscore pcs: #############################################################
if pzscore_pc:
for i in range(PC.shape[0]):
PC[i,:] = (PC[i,:]- PC[i,:].mean()) / PC[i,:].std()
###########################################################################
data = []
ev = 0
for count,fr in enumerate(PC):
label = 'pc%d' % int(count+1)
for (si,sj) in transitions:
# string label for type of transition:
sid = states[si] + states[sj]
if si == 3 and not ma_mode:
seq = sleepy.get_sequences(np.where(M==si)[0])
elif si==3 and ma_mode:
seq = sleepy.get_sequences(np.where(Mrepr==si)[0])
else:
seq = sleepy.get_sequences(np.where(Mrepr==si)[0])
for s in seq:
# ti is the last bin in the current sequence
ti = s[-1]
if si == 3 and ma_mode:
p = s[0]
p = p-1
while p >0 and M[p] == 3:
p = p-1
p = p+1
s = np.arange(p, ti+1)
# check if next state is sj; only then continue
if ti < len(M)-1 and Mrepr[ti+1] == sj:
# go into future
p = ti+1
if sj == 3:
# if sj == 3, we're treating MAs as NREM
while p<len(M)-1 and M[p] == sj:
p += 1
else:
while p<len(M)-1 and Mrepr[p] == sj:
p += 1
p -= 1
sj_idx = list(range(ti+1, p+1))
# so the indices of state si are seq
# the indices of state sj are sj_idx
if ti in allowed_idx and ipre <= ti < len(M)-ipost and len(s)*sdt >= si_threshold[si-1] and len(sj_idx)*sdt >= sj_threshold[sj-1]:
act = fr[ti-ipre+1:ti+ipost+1]
# Note: ti+1 is the first time point of the "post" state
# i = 10, ipre = 2, ipost = 2
# 8,9,10
# np.arange(8,12) = 8,9,10,11,12
data += zip([s[0]]*m, [label]*m, t, act, [sid]*m)
ev += 1
df = pd.DataFrame(data=data, columns=['event', 'pc', 'time', 'fr', 'trans'])
return df
def pc_transitions_laser(mouse, PC, M, transitions, pre, post, si_threshold, sj_threshold,
ma_thr=10, ma_rem_exception=False, sdt=2.5, ma_mode=False, rnd_laser=False,
kcuts=[], allowed_idx=[], pzscore_pc=True, config_file='', laser_dur=-1):
"""
Compare spontaneous and laser-induced brain state transitions
Parameters
----------
PC : np.array
Number of PCx x number of time bins
PCs; each row corresponds to one PC.
M : np.array
hynpogram.
transitions : TYPE
DESCRIPTION.
pre : TYPE
DESCRIPTION.
post : TYPE
DESCRIPTION.
si_threshold : TYPE
DESCRIPTION.
sj_threshold : TYPE
DESCRIPTION.
ma_thr : TYPE, optional
DESCRIPTION. The default is 10.
ma_rem_exception : bool, optional
If True, don't touch wake following REM sleep. The default is False.
sdt : float, optional
time bin duration in seconds of one brain state. The default is 2.5.
ma_mode : bool, optional
If True, then specifically analysis transitions from and to MAs.
Note that when considering transitions from and to NREM, MAs are
considered as NREM sleep.
kcuts : TYPE, optional
DESCRIPTION. The default is [].
allowed_idx : TYPE, optional
DESCRIPTION. The default is [].
Returns
-------
df : TYPE
DESCRIPTION.
"""
dt = 2.5
states = {1:'R', 2:'W', 3:'N', 4:'M'}
nhypno = M.shape[0]
ndim = PC.shape[0]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
ipre = int(np.round(pre/sdt))
ipost = int(np.round(post/sdt))
# m = ipre + ipost
t = np.arange(-ipre, ipost) * sdt
m = len(t)
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*sdt) <= ma_thr:
if ma_rem_exception:
if (s[0]>=1) and (M[s[0] - 1] != 1):
if ma_mode:
M[s] = 4
else:
M[s] = 3
else:
if ma_mode:
M[s] = 4
else:
M[s] = 3
Mrepr = M.copy()
M[M==4] = 3
if not ma_mode:
#just forget about MAs:
Mrepr = M
if len(allowed_idx) == 0:
allowed_idx = range(0, len(M))
#######################################################################
# get laser start and end index after excluding kcuts: ################
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*dt)
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
if laser_dur == -1:
idxe = [int(i/nbin) for i in idxe]
b = [i for i in idxe if i in tidx]
a = [i for i in idxs if i in tidx]
b = np.array(b)
a = np.array(a)
laser_dur = (np.mean(b-a) + 1) * dt
print('Laser duration: %f' % laser_dur)
else:
idxe = [(int(i + int(laser_dur/dt))) for i in idxs]
# randomize laser #########################################################
dur = int(laser_dur/dt)
if rnd_laser:
idxs_rnd = []
idxe_rnd = []
tmp = np.random.randint(dur, idxs[0]-dur)
idxs_rnd.append(tmp)
idxe_rnd.append(tmp+dur)
for (a,b) in zip(idxe[0:-1], idxs[1:]):
if a+2*dur < b-dur:
tmp = np.random.randint(a+2*dur,b-dur)
idxs_rnd.append(tmp)
idxe_rnd.append(tmp+dur)
idxs = idxs_rnd
idxe = idxe_rnd
###########################################################################
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
laser_idx = np.array(laser_idx)
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser_idx = laser_idx[laser_idx < nlsr]
laser[laser_idx] = 1
laser = laser[tidx]
# get again indices after kcut
laser_idx = np.where(laser == 1)[0]
idxs = [s[0] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
idxe = [s[-1] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
#######################################################################
# zscore pcs:
if pzscore_pc:
for i in range(PC.shape[0]):
PC[i,:] = (PC[i,:]- PC[i,:].mean()) / PC[i,:].std()
data = []
ev = 0
for count,fr in enumerate(PC):
label = 'pc%d' % int(count+1)
for (si,sj) in transitions:
# string label for type of transition:
sid = states[si] + states[sj]
if si == 3 and not ma_mode:
seq = sleepy.get_sequences(np.where(M==si)[0])
elif si==3 and ma_mode:
seq = sleepy.get_sequences(np.where(Mrepr==si)[0])
else:
seq = sleepy.get_sequences(np.where(Mrepr==si)[0])
for s in seq:
# ti is the last bin in the current sequence
ti = s[-1]
if si == 3 and ma_mode:
p = s[0]
p = p-1
while p > 0 and M[p] == 3:
p = p-1
p = p+1
s = np.arange(p, ti+1)
# check if next state is sj; only then continue
if ti < len(M)-1 and Mrepr[ti+1] == sj:
# go into future
p = ti+1
if sj == 3:
# if sj == 3, we're treating MAs as NREM
while p<len(M)-1 and M[p] == sj:
p += 1
else:
while p<len(M)-1 and Mrepr[p] == sj:
p += 1
p -= 1
sj_idx = list(range(ti+1, p+1))
# so the indices of state si are seq
# the indices of state sj are sj_idx
if ti in allowed_idx and ipre <= ti < len(M)-ipost and len(s)*sdt >= si_threshold[si-1] and len(sj_idx)*sdt >= sj_threshold[sj-1]:
act = fr[ti-ipre+1:ti+ipost+1]
# Note: ti+1 is the first time point of the "post" state
# i = 10, ipre = 2, ipost = 2
# 8,9,10
# np.arange(8,12) = 8,9,10,11,12
delay = -1
laser_on = 'no'
if ti+1 in laser_idx:
laser_on = 'yes'
a = ti
while laser[a] == 1:
a = a-1
a = a+1
delay = (ti - a + 1) * dt
laser_cut = laser[ti-ipre+1:ti+ipost+1]
data += zip([s[0]]*m, [label]*m, t, act, [sid]*m, [laser_on]*m, laser_cut, [delay]*m)
ev += 1
df = pd.DataFrame(data=data, columns=['event', 'pc', 'time', 'fr', 'trans', 'laser_on', 'laser', 'delay'])
return df
def fr_svd(units, nsmooth=0, pzscore=False):
"""
perform SVD on firing rates of a population of neurons
Parameters
----------
units : pd.DataFrame
each column is a unit; the column names are the unitIDs.
ndim : int, optional
DESCRIPTION. The default is 3.
nsmooth : float, optional
If > 0, som. The default is 0.
Returns
-------
PC : np.array
Each row vector corresponds to one principal component.
dimensions: ndim x timepoints
V : np.array
Variance captured by the principal components
"""
# first transform pd.DataFrame units to np.array:
# Matrix arrangement:
# rows - neurons; columns - time points
# NOTE: the rows (neurons) are the variables (dimensions),
# the time points are the samples (trials)
# through PCA we want to keep the number of samples, but we want to reduce
# the dimensions. Again, for our matrix R, we have the arrangement:
# samples (=time)
# variables x
#
# In our case, that's
# timepoints
# units x
unitIDs = [unit for unit in units.columns if re.split('_', unit)[1] == 'good']
nsample = units.shape[0] # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
i = 0
for unit in unitIDs:
R[i,:] = sleepy.smooth_data(np.array(units[unit]),nsmooth)
if pzscore:
R[i,:] = (R[i,:] - R[i,:].mean()) / R[i,:].std()
i += 1
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
for i in range(nvar):
R[i,:] = R[i,:] - R[i,:].mean()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R.T / np.sqrt(nsample-1)
U,S,Vh = scipy.linalg.svd(Y)
return U, S, Vh
def pr_components(units, ndim=3, nsmooth=0, pzscore=False):
"""
Calculate principal components of simultaneously recorded firing rates.
Parameters
----------
units : pd.DataFrame
each column is a unit; the column names are the unitIDs.
ndim : int, optional
DESCRIPTION. The default is 3.
nsmooth : float, optional
If > 0, som. The default is 0.
Returns
-------
PC : np.array
Each row vector corresponds to one principal component.
dimensions: ndim x timepoints
V : np.array
Variance captured by the principal components
"""
# first transform pd.DataFrame units to np.array:
# Matrix arrangement:
# rows - neurons; columns - time points
# NOTE: the rows (neurons) are the variables (dimensions),
# the time points are the samples (trials)
# through PCA we want to keep the number of samples, but we want to reduce
# the dimensions. Again, for our matrix R, we have the arrangement:
# samples (=time)
# variables x
#
# In our case, that's
# timepoints
# units x
unitIDs = [unit for unit in units.columns if re.split('_', unit)[1] == 'good']
nsample = units.shape[0] # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
i = 0
for unit in unitIDs:
R[i,:] = sleepy.smooth_data(np.array(units[unit]),nsmooth)
if pzscore:
R[i,:] = (R[i,:] - R[i,:].mean()) / R[i,:].std()
i += 1
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
for i in range(nvar):
R[i,:] = R[i,:] - R[i,:].mean()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R.T / np.sqrt(nsample-1)
U,S,Vh = scipy.linalg.svd(Y)
# project data @R onto principal direction @Vh
# Vh: nvar x nvar; each row is an eigenvector of the covariance matrix.
# So using Vh * R, we project the data into the covariance space.
PC = np.dot(Vh, R)[0:ndim,:]
V = S**2
# Side note; the same result you can get also with:
# U * S ... in more detail:
# SM = np.zeros((nsample, nsample))
# for i in range(nvar):
# SM[i,i] = S[i]
# PC = np.dot(U, SM) * np.sqrt(nsample-1)
###########################################################################
# plots
plt.figure()
plt.subplot(211)
plt.plot(V, '.')
plt.subplot(212)
plt.plot(PC.T)
var = S**2
var_total = np.sum(S**2)
# Calculate the cumulative variance explained, i.e. how much
# of the variance is captured by the first i principal components.
p = []
for i in range(1,len(var)+1):
s = np.sum(var[0:i])
p.append(s/var_total)
plt.figure()
plt.plot(p, '.')
return PC, V
def sleep_components(units, M, wake_dur=600, wake_break=60, ndim=3, nsmooth=0,
pzscore=False, ddof=0, detrend=False, pplot=True, pc_sign=[], scale_pc=False,
kcuts=[], mode='w', mouse='', collapse=False,
ylim=[], config_file='mouse_config.txt'):
"""
Calculate principal components, excluding long wake periods in the recording.
Long wake periods are defined by the parameters $wake_dur (duration of wake episodes)
and $wake_break (wake episodes separated by less than $wake_break seconds are fused).
Parameters
----------
units : pd.DataFrame
each column corresponds to a unit; each row is a time bin
M : np.array
hypnogram.
wake_dur : float, optional
Exclude wake periods that are longer than $wake_dur seconds. The default is 600.
wake_break : float, optional
Two wake periods that are separated by less than $wake_break seconds
are merged to one period. The default is 60.
ndim : int, optional
reduce data (matrix of firing rate vector) to $ndim dimensions using PCA. The default is 3.
nsmooth : float, optional
Smooth firing rate vector. The default is 0.
pzscore : boolean, optional
If True, zscore firing rates. The default is False.
detrend : boolean, optional
If True, detrend each firing rate vector.
pc_sign : list with $ndim elements, either 1 or -1.
The sign of PCs is ambiguous, so if preferred multiply, PC i with pc_sign[i]
scale_pc: if True, scale PCs by 1 / np.sqrt(nsample - 1); nsample is the number of time bins used
for PC calculation; Vh * S /sqrt(nsample-1) is conventionally used as PCs
ddof: degrees of freedom for z-score calculation;
kcuts : list of tuples or lists with two elements.
Discard the time interval ranging from kcuts[i][0] to kcuts[i][1] seconds
mode : string with characters 'r' and/or 'w'
if 'r' in mode, remove all REM indices for PC computation
if 'w' in mode, remove all Wake indices for PC computation
collapse: boolean
if True, plot each PC in its own axis
ylim: empty list, or tuple
if empty list, don't fix ylims, otherwise set plt.ylim(ylim) for each PC
config_file: str
mouse configuration file as loaded by &load_config()
mouse: str
To plot laser, set $mouse to mouse name
Returns
-------
PC : np.array
The $ndim principal components. Note although the PCs have been calculated
without long wake periods, the returned PCs do include all wake periods.
V : np.array
Eigenvalues of the covariance matrix = Variance associated with each PC
Vh : np.array
each row in Vh is an eigenvector of the covariance matrix
idx : np.array
Indices of time bins used for PC calculation.
NOTE that @idx are the indices obtained
AFTER cutting out the KCUT intervals!
"""
dt = 2.5
nhypno = np.min((len(M), units.shape[0]))
Morig = M.copy()
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# KCUT ####################################################################
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
print(len(tidx))
###########################################################################
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nsample = nhypno # number of time points
nvar = len(unitIDs) # number of units
# Note the number of samples stays the same; while the number is variables
# (or dimensions) is reduced! We want to keep the same number of time points,
# but have only a few 'modes'.
R = np.zeros((nvar, nsample))
#@tidx are the indices we're further considering.
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]),nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
R[i,:] = (tmp - tmp.mean()) / tmp.std(ddof=ddof)
else:
R[i,:] = tmp
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
#for i in range(nvar):
# R[i,:] = R[i,:] - R[i,:].mean()
# find long wake blocks:
widx = sleepy.get_sequences(np.where(M==2)[0], ibreak=int(wake_break/dt))
# all REM sequences
ridx = sleepy.get_sequences(np.where(M==1)[0])
nidx = sleepy.get_sequences(np.where(M==3)[0])
tmp = []
for w in widx:
if len(w) * dt > wake_dur:
tmp += list(w)
widx = tmp
tmp = []
for r in ridx:
tmp += list(r)
ridx = tmp
tmp = []
for r in nidx:
tmp += list(r)
nidx = tmp
nhypno = np.min((len(M), units.shape[0]))
idx = np.arange(0, np.min((len(M), units.shape[0])))
if 'w' in mode:
idx = np.setdiff1d(idx, widx)
else:
widx = []
if 'r' in mode:
idx = np.setdiff1d(idx, ridx)
if 'n' in mode:
idx = np.setdiff1d(idx, nidx)
else:
ridx = []
### NEW: 5/14/25:
nsample = len(idx)
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R[:,idx].T / np.sqrt(nsample-1)
# make sure the the columns of Y are mean-zero
for i in range(nvar):
Y[:,i] = Y[:,i] - Y[:,i].mean()
# redundant (see few lines below $mu); removed the line on 5/14/25
# R[i,:] = R[i,:] - R[i,idx].mean()
# SVD
# Note that each column is a neuron,
# and each row is a time point
U,S,Vh = scipy.linalg.svd(Y)
# each row in Vh is an eigenvector of the COV matrix;
# to get the PCs, we project R onto the eigenvectors:
# NEW 05/14/25!!!
#PC = np.dot(Vh, R)[0:ndim,:]
mu = R[:, idx].mean(axis=1, keepdims=True)
sigma = R[:, idx].std(axis=1, ddof=1, keepdims=True) if pzscore else 1.0
R0 = (R - mu) / sigma
if scale_pc:
PC = (Vh[:ndim] @ R0) / np.sqrt(nsample - 1)
else:
PC = np.dot(Vh, R0)[0:ndim,:]
V = S**2
if len(pc_sign) > 0:
i = 0
for s in pc_sign:
PC[i,:] = PC[i,:] * s
i += 1
if pplot:
add_laser = False
t = np.arange(0, nhypno) * dt
plt.figure()
tmp = widx+ridx
tmp.sort()
widx = sleepy.get_sequences(np.array(tmp))
axes_exc = plt.axes([0.2, 0.9, 0.7, 0.05])
for w in widx:
if len(w) > 1:
if w[-1] < len(M):
plt.plot([t[w[0]], t[w[-1]]], [1, 1], 'k', lw=2)
plt.ylim([0, 2])
plt.xlim((t[0], t[-1]))
# if laser exists also add laser here --
if mouse:
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*2.5)
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
idxe = [int(i/nbin) for i in idxe]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser[laser_idx] = 1
laser = laser[tidx]
lsr_seq = sleepy.get_sequences(np.where(laser==1)[0])
for w in lsr_seq:
if w[-1] < len(M):
plt.plot([t[w[0]], t[w[-1]]], [1.5, 1.5], 'b', lw=2)
plt.xlim((t[0], t[-1]))
add_laser = True
sleepy._despine_axes(axes_exc)
axes_brs = plt.axes([0.2, 0.85, 0.7, 0.05], sharex=axes_exc)
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
sleepy._despine_axes(axes_brs)
if collapse:
plt.axes([0.2, 0.2, 0.7, 0.6], sharex=axes_brs)
plt.plot(t, PC[:,:].T)
plt.xlim((t[0], t[-1]))
sns.despine()
plt.xlabel('Time (s)')
plt.ylabel('PC')
else:
clrs = sns.color_palette("husl", ndim)
d = (0.6 / ndim) * 0.3
ny = (0.6 / ndim)-d
for i in range(ndim):
ax = plt.axes([0.2, 0.2+i*(ny+d), 0.7, ny], sharex=axes_exc)
ax.plot(t, PC[ndim-1-i,:], color=clrs[i])
plt.xlim([t[0], t[-1]])
sleepy.box_off(ax)
plt.ylabel('PC%d' % (ndim-i))
if add_laser:
ylims = ax.get_ylim()
for w in lsr_seq:
if w[-1] < len(M):
if 1 in M[w]:
color = 'b'
else:
color = 'r'
#if M[w[0]] == 3:
#plt.plot(t[w[0]], PC[ndim-1-i,w[0]-1:w[0]+1].mean(), '.', color=color, lw=2)
plt.plot(t[w[0]], PC[ndim-1-i,w[0]], '.', color=color, lw=2)
laser_tend = (w[-1] - w[0]) * dt
dy = ylims[1] - ylims[0]
ax.add_patch(patches.Rectangle(
(t[w[0]], ylims[0]), laser_tend, dy,
facecolor=[0.6, 0.6, 1], edgecolor=[0.6, 0.6, 1]))
if i > 0:
ax.spines["bottom"].set_visible(False)
ax.axes.get_xaxis().set_visible(False)
else:
plt.xlabel('Time (s)')
if len(ylim) > 0:
plt.ylim(ylim)
var = S**2
var_total = np.sum(S**2)
# Calculate the cumulative variance explained, i.e. how much
# of the variance is captured by the first i principal components.
p = []
for i in range(1,len(var)+1):
s = np.sum(var[0:i])
p.append(s/var_total)
plt.figure(figsize=(4,4))
plt.plot(p, '.', color='gray')
plt.xlabel(r'$\mathrm{PC_i}$')
plt.ylabel('Cum. variance')
plt.subplots_adjust(bottom=0.2, left=0.2)
plt.ylim([0, 1.1])
sns.despine()
if len(kcuts) > 0:
# find long wake blocks:
widx = sleepy.get_sequences(np.where(Morig==2)[0], ibreak=int(wake_break/dt))
tmp = []
for w in widx:
if len(w) * dt > wake_dur:
tmp += list(w)
widx = tmp
nhypno = np.min((len(Morig), units.shape[0]))
idx_total = np.arange(0, nhypno)
idx_total = np.setdiff1d(idx_total, kidx)
idx_total = np.setdiff1d(idx_total, widx)
return PC, V, Vh, idx
def sleep_components_laser(units, M, wake_dur=600, wake_break=60, ndim=3, nsmooth=0,
pzscore=False, detrend=False, pplot=True, pc_sign=[],
kcuts=[], mode='w', mouse='', collapse=False,
ylim=[], config_file='mouse_config.txt'):
"""
Similar to &sleep_components(), but allows for discarding laser intervals
for PC estimation.
Calculate principal components, excluding long wake periods in the recording.
Long wake periods are defined by the parameters $wake_dur (duration of wake episodes)
and $wake_break (wake episodes separated by less than $wake_break seconds are fused).
Note: To make sure that the laser is read, set $mouse to the mouse name
Parameters
----------
units : pd.DataFrame
each column corresponds to a unit; each row is a time bin
M : np.array
hypnogram.
wake_dur : float, optional
Exclude wake periods that are longer than $wake_dur seconds. The default is 600.
wake_break : float, optional
Two wake periods that are separated by less than $wake_break seconds
are merged to one period. The default is 60.
ndim : int, optional
reduce data (matrix of firing rate vector) to $ndim dimensions using PCA. The default is 3.
nsmooth : float, optional
Smooth firing rate vector. The default is 0.
pzscore : boolean, optional
If True, zscore firing rates. The default is False.
detrend : boolean, optional
If True, detrend each firing rate vector.
pc_sign : list with $ndim elements, either 1 or -1.
The sign of PCs is ambiguous, so if preferred multiply, PC i with pc_sign[i]
kcuts : list of tuples or lists with two elements.
Discard the time interval ranging from kcuts[i][0] to kcuts[i][1] seconds
mode : string with characters 'r' and/or 'w'
if 'r' in mode, remove all REM indices for PC computation
if 'w' in mode, remove all Wake indices for PC computation
collapse: boolean
if True, plot each PC in its own axis
ylim: empty list, or tuple
if empty list, don't fix ylims, otherwise set plt.ylim(ylim) for each PC
config_file: str
mouse configuration file as loaded by &load_config()
Returns
-------
PC : np.array
The $ndim principal components. Note although the PCs have been calculated
without long wake periods, the returned PCs do include all wake periods.
V : np.array
Eigenvalues of the covariance matrix = Variance associated with each PC
Vh : np.array
each row in Vh is an eigenvector of the covariance matrix
idx : np.array
Indices of time bins used for PC calculation.
NOTE that @idx are the indices obtained
AFTER cutting out the KCUT intervals!
"""
dt = 2.5
nhypno = np.min((len(M), units.shape[0]))
Morig = M.copy()
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# KCUT ####################################################################
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
print(len(tidx))
###########################################################################
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nsample = nhypno # number of time points
nvar = len(unitIDs) # number of units
# Note the number of samples stays the same; while the number is variables
# (or dimensions) is reduced! We want to keep the same number of time points,
# but have only a few 'modes'.
R = np.zeros((nvar, nsample))
#@tidx are the indices we're further considering.
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]),nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
R[i,:] = (tmp - tmp.mean()) / tmp.std()
else:
R[i,:] = tmp
# find long wake blocks:
widx = sleepy.get_sequences(np.where(M==2)[0], ibreak=int(wake_break/dt))
# all REM sequences
ridx = sleepy.get_sequences(np.where(M==1)[0])
nidx = sleepy.get_sequences(np.where(M==3)[0])
tmp = []
for w in widx:
if len(w) * dt > wake_dur:
tmp += list(w)
widx = tmp
tmp = []
for r in ridx:
tmp += list(r)
ridx = tmp
tmp = []
for r in nidx:
tmp += list(r)
nidx = tmp
nhypno = np.min((len(M), units.shape[0]))
idx = np.arange(0, np.min((len(M), units.shape[0])))
if 'w' in mode:
idx = np.setdiff1d(idx, widx)
else:
widx = []
if 'r' in mode:
idx = np.setdiff1d(idx, ridx)
if 'n' in mode:
idx = np.setdiff1d(idx, nidx)
else:
ridx = []
# collect laser information ###############################################
if mouse:
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*2.5)
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
idxe = [int(i/nbin) for i in idxe]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser[laser_idx] = 1
laser = laser[tidx]
laser_idx = np.where(laser == 1)[0]
print('Removing laser indices from mouse %s' % mouse)
idx = np.setdiff1d(idx, laser_idx)
###########################################################################
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R[:,idx].T / np.sqrt(nsample-1)
# make sure the the columns of Y are mean-zero
for i in range(nvar):
Y[:,i] = Y[:,i] - Y[:,i].mean()
R[i,:] = R[i,:] - R[i,idx].mean()
# SVD
# Note that each column is a neuron,
# and each row is a time point
U,S,Vh = scipy.linalg.svd(Y)
# each row in Vh is an eigenvector of the COV matrix;
# to get the PCs, we project R onto the eigenvectors:
PC = np.dot(Vh, R)[0:ndim,:]
V = S**2
if len(pc_sign) > 0:
i = 0
for s in pc_sign:
PC[i,:] = PC[i,:] * s
i += 1
if pplot:
add_laser = False
t = np.arange(0, nhypno) * dt
plt.figure()
tmp = widx+ridx
tmp.sort()
widx = sleepy.get_sequences(np.array(tmp))
axes_exc = plt.axes([0.2, 0.9, 0.7, 0.05])
for w in widx:
if len(w) > 1:
if w[-1] < len(M):
plt.plot([t[w[0]], t[w[-1]]], [1, 1], 'k', lw=2)
plt.ylim([0, 2])
plt.xlim((t[0], t[-1]))
# if laser exists also add laser here --
if mouse:
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*2.5)
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
idxe = [int(i/nbin) for i in idxe]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser[laser_idx] = 1
laser = laser[tidx]
lsr_seq = sleepy.get_sequences(np.where(laser==1)[0])
for w in lsr_seq:
if w[-1] < len(M):
plt.plot([t[w[0]], t[w[-1]]], [1.5, 1.5], 'b', lw=2)
plt.xlim((t[0], t[-1]))
add_laser = True
sleepy._despine_axes(axes_exc)
axes_brs = plt.axes([0.2, 0.85, 0.7, 0.05], sharex=axes_exc)
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
sleepy._despine_axes(axes_brs)
if collapse:
plt.axes([0.2, 0.2, 0.7, 0.6], sharex=axes_brs)
plt.plot(t, PC[:,:].T)
plt.xlim((t[0], t[-1]))
sns.despine()
plt.xlabel('Time (s)')
plt.ylabel('PC')
else:
clrs = sns.color_palette("husl", ndim)
d = (0.6 / ndim) * 0.3
ny = (0.6 / ndim)-d
for i in range(ndim):
ax = plt.axes([0.2, 0.2+i*(ny+d), 0.7, ny], sharex=axes_exc)
ax.plot(t, PC[ndim-1-i,:], color=clrs[i])
plt.xlim([t[0], t[-1]])
sleepy.box_off(ax)
plt.ylabel('PC%d' % (ndim-i))
if add_laser:
ylims = ax.get_ylim()
for w in lsr_seq:
if w[-1] < len(M):
if 1 in M[w]:
color = 'b'
else:
color = 'r'
#if M[w[0]] == 3:
#plt.plot(t[w[0]], PC[ndim-1-i,w[0]-1:w[0]+1].mean(), '.', color=color, lw=2)
plt.plot(t[w[0]], PC[ndim-1-i,w[0]], '.', color=color, lw=2)
laser_tend = (w[-1] - w[0]) * dt
dy = ylims[1] - ylims[0]
ax.add_patch(patches.Rectangle(
(t[w[0]], ylims[0]), laser_tend, dy,
facecolor=[0.6, 0.6, 1], edgecolor=[0.6, 0.6, 1]))
if i > 0:
ax.spines["bottom"].set_visible(False)
ax.axes.get_xaxis().set_visible(False)
else:
plt.xlabel('Time (s)')
if len(ylim) > 0:
plt.ylim(ylim)
var = S**2
var_total = np.sum(S**2)
# Calculate the cumulative variance explained, i.e. how much
# of the variance is captured by the first i principal components.
p = []
for i in range(1,len(var)+1):
s = np.sum(var[0:i])
p.append(s/var_total)
plt.figure(figsize=(4,4))
plt.plot(p, '.', color='gray')
plt.xlabel(r'$\mathrm{PC_i}$')
plt.ylabel('Cum. variance')
plt.subplots_adjust(bottom=0.2, left=0.2)
plt.ylim([0, 1.1])
sns.despine()
if len(kcuts) > 0:
# find long wake blocks:
widx = sleepy.get_sequences(np.where(Morig==2)[0], ibreak=int(wake_break/dt))
tmp = []
for w in widx:
if len(w) * dt > wake_dur:
tmp += list(w)
widx = tmp
nhypno = np.min((len(Morig), units.shape[0]))
idx_total = np.arange(0, nhypno)
idx_total = np.setdiff1d(idx_total, kidx)
idx_total = np.setdiff1d(idx_total, widx)
return PC, V, Vh, idx
def sleep_components_fine(ids, wake_dur=600, wake_break=60, ndim=3, nsmooth=0,
pzscore=False, detrend=False, pplot=True, pc_sign=[], ndown=250,
kcuts=[], mode='w', ppath='', mouse='', collapse=False,
ylim=[], config_file='mouse_config.txt'):
dt = 2.5
NDOWN = ndown
NUP = int(dt / (0.001 * NDOWN))
fine_scale = True
if len(config_file) == 0:
config_file = 'mouse_config.txt'
path = load_config(config_file)[mouse]['SL_PATH']
ppath, file = os.path.split(path)
M = sleepy.load_stateidx(ppath, file)[0]
###########################################################################
recs = load_config(config_file)[mouse]
if 'TR_PATH' in recs:
tr_path = load_config(config_file)[mouse]['TR_PATH']
else:
tr_path = load_config(config_file)[mouse]['NP_PATH']
units = np.load(os.path.join(tr_path,'1k_train.npz'))
unitIDs = [unit for unit in list(units.keys()) if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
# re-scale dt
dt = dt / NUP
M = upsample_mx(M, NUP)
nhypno = int(np.min((len(M), units[unitIDs[0]].shape[0]/NDOWN)))
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
nsample = nhypno # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
#@tidx are the indices we're further considereing.
print('Starting downsampling, smoothing and z-scoring')
fr_file = os.path.join(tr_path, 'fr_fine_ndown%d.mat' % NDOWN)
if not os.path.isfile(fr_file):
for i,unit in enumerate(unitIDs):
tmp = sleepy.downsample_vec(np.array(units[unit]), NDOWN)
R[i,:] = tmp[tidx]
so.savemat(fr_file, {'R':R, 'ndown':NDOWN})
else:
# NOTE: the save matrix are the responses AFTER KCUT!
R = so.loadmat(fr_file, squeeze_me=True)['R']
for i,unit in enumerate(unitIDs):
tmp = R[i,:]
tmp = sleepy.smooth_data(tmp, nsmooth)
if pzscore:
R[i,:] = (tmp[:] - tmp[:].mean()) / tmp[:].std()
else:
R[i,:] = tmp[:]
# find long wake blocks:
widx = sleepy.get_sequences(np.where(M==2)[0], ibreak=int(wake_break/dt))
# all REM sequences
ridx = sleepy.get_sequences(np.where(M==1)[0])
nidx = sleepy.get_sequences(np.where(M==3)[0])
tmp = []
for w in widx:
if len(w) * dt > wake_dur:
tmp += list(w)
widx = tmp
tmp = []
for r in ridx:
tmp += list(r)
ridx = tmp
tmp = []
for r in nidx:
tmp += list(r)
nidx = tmp
nhypno = np.min((len(M), R.shape[1]))
idx = np.arange(0, nhypno)
if 'w' in mode:
idx = np.setdiff1d(idx, widx)
else:
widx = []
if 'r' in mode:
idx = np.setdiff1d(idx, ridx)
if 'n' in mode:
idx = np.setdiff1d(idx, nidx)
else:
ridx = []
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R[:,idx].T / np.sqrt(nsample-1)
# make sure the the columns of Y are mean-zero
for i in range(nvar):
Y[:,i] = Y[:,i] - Y[:,i].mean()
R[i,:] = R[i,:] - R[i,idx].mean()
# SVD
# Note that each column is a neuron,
# and each row is a time point
U,S,Vh = scipy.linalg.svd(Y)
# each row in Vh is an eigenvector of the COV matrix;
# to get the PCs, we project R onto the eigenvectors:
PC = np.dot(Vh, R)[0:ndim,:]
V = S**2
if len(pc_sign) > 0:
i = 0
for s in pc_sign:
PC[i,:] = PC[i,:] * s
i += 1
## add figure
if pplot:
t = np.arange(0, nhypno) * dt
plt.figure()
axes_brs = plt.axes([0.2, 0.85, 0.7, 0.05])
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
sleepy._despine_axes(axes_brs)
plt.axes([0.2, 0.2, 0.7, 0.6], sharex=axes_brs)
plt.plot(t, PC[:,:].T)
plt.xlim((t[0], t[-1]))
sns.despine()
plt.xlabel('Time (s)')
plt.ylabel('PC')
return PC, V
def kcut_idx(M, X, kcuts, dt=2.5):
"""
Using the values defined in the mouse_config.txt file (field KCUT:), determine
the indices used for further calculations.
@param: np.array, brainstate sequence
@param X: np.pandas or np.array, array or DataFrame with time axis using same binning as @M
@param kcuts: list of tuples, areas at the beginning or end that should be discarded
@return tidx: np.array, list of indices in @M used for further calculation, i.e. indices
that are NOT within the ranges defined in @kcuts.
"""
#n = np.max(X.shape)
#nhypno = np.min((len(M), n))
nhypno = len(M)
#M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
return tidx
def kcut_idx2(M, kcuts, X = [], dt=2.5):
"""
Using the values defined in the mouse_config.txt file (field KCUT:), determine
the indices used for further calculations.
@param: np.array, brainstate sequence
@param X: np.pandas or np.array, array or DataFrame with time axis using same binning as @M
@param kcuts: list of tuples, areas at the beginning or end that should be discarded
@return tidx: np.array, list of indices in @M used for further calculation, i.e. indices
that are NOT within the ranges defined in @kcuts.
"""
nhypno = len(M)
if len(X) > 0:
n = np.max(X.shape)
nhypno = np.min((len(M), n))
tidx = np.arange(0, nhypno)
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
return tidx
def align_pcsign(PC, mouse, rem_pca=0, kcuts=[], align_pc3=True, config_file='', pnorm_spec=True):
"""
Automatically determine the sign of the given PCs.
Fix PC1 (normally increased during REM) though activity during REM sleep.
Fix PC2 (normally positively correlated with sigma power) through its
correlation with the sigma power.
Fix PC3...PCn the same way as PC1.
Parameters
----------
PC : np.array
Matrix with PCs; each row corresponds to one PC.
mouse : str
mouse name.
rem_pca : int, optional
0 or 1. If 0, assume that PC1 (PC[0,:]) is the "REM-PC".
kcuts : list of tuples
Areas at beginning or end to remove from recording. The default is [].
align_pc3 : bool, optional
If True, adjust sign of PC3, ... PCn using the same strategy as for PC1
config_file : str, optional
File name of mouse config file. The default is ''.
pnorm_spec : bool, optional
If True, normalize EEG spectrogram. The default is True.
Returns
-------
pc_sign : list of length PC.shape[0]
1 or -1 depending on whether the orientation of the PC should be changed or not.
"""
ndim = PC.shape[0]
pc_sign = [1 for i in range(ndim)]
if len(config_file) == 0:
config_file = 'mouse_config.txt'
path = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(path)
M = sleepy.load_stateidx(ppath, name)[0]
sigma = [10,15]
P = so.loadmat(os.path.join(ppath, name, 'sp_%s.mat'%name), squeeze_me=True)
SP = P['SP']
freq = P['freq']
dfreq = freq[1]-freq[0]
isigma = np.where((freq>=sigma[0])&(freq<=sigma[1]))[0]
# cut out kcuts: ###############
tidx = kcut_idx(M, PC, kcuts)
if tidx[-1] > PC.shape[1]-1:
tidx = tidx[0:-1]
M = M[tidx]
#PC = PC[:,tidx]
SP = SP[:,tidx]
################################
if SP.shape[1] != PC.shape[1]:
print('Check kcut!!!')
print('Shape of SP = %d; shape of PC = %d' % (SP.shape[1], PC.shape[1]))
if rem_pca == 0:
sigma_pca = 1
else:
sigma_pca = 0
# fix PC1
mmin = np.min((len(M), PC.shape[1]))
M = M[0:mmin]
PC = PC[:,0:mmin]
pc1 = PC[rem_pca,:]
rem_idx = np.where(M==1)[0]
a = pc1[rem_idx].mean()
if a < 0:
pc_sign[rem_pca] = -1
# fix PC3...PCn the same way
if align_pc3:
if PC.shape[0] >= 3:
for j in range(2, PC.shape[0]):
pc3 = PC[j,:]
rem_idx = np.where(M==1)[0]
a = pc3[rem_idx].mean()
if a < 0:
pc_sign[j] = -1
# fix PC2 through correlation with sigma power
if pnorm_spec:
sp_mean = SP.mean(axis=1)
SP = np.divide(SP, np.tile(sp_mean, (SP.shape[1], 1)).T)
sigma_pow = SP[isigma,:].mean(axis=0)
else:
sigma_pow = SP[isigma, :].sum(axis=0)*dfreq
CC, t = state_correlation(PC[sigma_pca,:], sigma_pow, M, win=60, pplot=False)
CC = np.nanmean(CC, axis=0)
i = np.argmax(np.abs(CC))
if CC[i] < 0:
pc_sign[sigma_pca] = -1
return pc_sign
def plot_pcs_withsigma(PC, M, mouse, ndim=2, pc_sign=[], kcuts=[],
ma_thr=10, ma_rem_exception=False,
tstart=0, tend=-1,
dt=2.5, sigma=[10,15], fmax=20, vm=[],
box_filt=[], pnorm_spec=True,
reverse_pcs=False, tlegend=120, zoomin=[], r_mu=[10,200],
pplot=True, config_file=''):
"""
Plot principal components along with sigma power.
Parameters
----------
PC : TYPE
DESCRIPTION.
M : TYPE
DESCRIPTION.
mouse : TYPE
DESCRIPTION.
ndim : TYPE, optional
DESCRIPTION. The default is 2.
pc_sign : TYPE, optional
DESCRIPTION. The default is [].
kcuts : TYPE, optional
DESCRIPTION. The default is [].
ma_thr : TYPE, optional
DESCRIPTION. The default is 10.
ma_rem_exception : TYPE, optional
DESCRIPTION. The default is False.
tstart : TYPE, optional
DESCRIPTION. The default is 0.
tend : TYPE, optional
DESCRIPTION. The default is -1.
dt : TYPE, optional
DESCRIPTION. The default is 2.5.
sigma : TYPE, optional
DESCRIPTION. The default is [10,15].
fmax : TYPE, optional
DESCRIPTION. The default is 20.
vm : list, optional
List with two elements defining vmin and vmax for the EEG spectrogram colormap.
If [], matplotlib will automatically set the color range.
box_filt : list, optional
Dimensions of box filter to smooth EEG spectrogram. If [], no filteringThe default is [].
pnorm_spec : bool, optional
If True, normalize spectrogram.
reverse_pcs : TYPE, optional
DESCRIPTION. The default is False.
tlegend : TYPE, optional
DESCRIPTION. The default is 120.
zoomin : TYPE, optional
DESCRIPTION. The default is [].
pplot : TYPE, optional
DESCRIPTION. The default is True.
config_file : TYPE, optional
DESCRIPTION. The default is ''.
Returns
-------
PC_orig : TYPE
DESCRIPTION.
M : TYPE
DESCRIPTION.
sigma_pow : TYPE
DESCRIPTION.
"""
if len(r_mu) == 0:
no_emg = True
else:
no_emg = False
if len(config_file) == 0:
config_file = 'mouse_config.txt'
path = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(path)
M = sleepy.load_stateidx(ppath, name)[0]
# load EEG spectrogram
P = so.loadmat(os.path.join(ppath, name, 'sp_%s.mat'%name), squeeze_me=True)
SP = P['SP']
# load EMG spectrogram
P = so.loadmat(os.path.join(ppath, name, 'msp_%s.mat'%name), squeeze_me=True)
SPM = P['mSP']
# cut out kcuts: ###############
tidx = kcut_idx(M, PC, kcuts)
tidx = tidx[0:PC.shape[1]]
M = M[tidx]
#PC = PC[:,tidx]
SP = SP[:,tidx]
if not(len(M) == SP.shape[1] == PC.shape[1]):
print('Something went wrong with KCUT')
print('returning')
return
################################
# flatten out MAs #########################################################
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
###########################################################################
freq = P['freq']
dfreq = freq[1]-freq[0]
isigma = np.where((freq>=sigma[0])&(freq<=sigma[1]))[0]
if len(box_filt) > 0:
filt = np.ones(box_filt)
filt = np.divide(filt, filt.sum())
SP = scipy.signal.convolve2d(SP, filt, boundary='symm', mode='same')
if pnorm_spec:
sp_mean = SP.mean(axis=1)
SP = np.divide(SP, np.tile(sp_mean, (SP.shape[1], 1)).T)
sigma_pow = SP[isigma,:].mean(axis=0)
else:
sigma_pow = SP[isigma, :].sum(axis=0)*dfreq
if not no_emg:
#EMG amplitude
i_mu = np.where((freq >= r_mu[0]) & (freq <= r_mu[1]))[0]
p_mu = np.sqrt(SPM[i_mu, :].sum(axis=0) * dfreq)
istart = int(tstart/dt)
if tend == -1:
iend = len(M)
else:
iend = int(tend/dt)
# cut out time interval istart - iend:
M = M[istart:iend]
PC = PC[:,istart:iend]
PC_orig = PC.copy()
PC_orig = PC_orig[istart:iend]
sigma_pow = sigma_pow[istart:iend]
p_mu = p_mu[istart:iend]
t = np.arange(0, len(M))*dt
if pplot:
plt.figure()
if no_emg:
### define all axes:
axes_brs = plt.axes([0.1, 0.9, 0.8, 0.05])
axes_spec = plt.axes([0.1, 0.72, 0.8, 0.15], sharex=axes_brs)
# axes for colorbar
axes_cbar = plt.axes([0.9, 0.72, 0.05, 0.15])
axes_sig = plt.axes([0.1, 0.57, 0.8, 0.1], sharex=axes_brs)
axes_pcs = plt.axes([0.1,0.05,0.8,0.5], sharex=axes_brs)
axes_legend = plt.axes([0.1,0.04,0.8,0.04], sharex=axes_brs)
else:
# -----------------------------------------------------------------
# Define a vertical stack: (rel. heights add up to 1.00)
# -----------------------------------------------------------------
lay = {
"axes_brs" : 0.05, # hypnogram / brs bar
"axes_spec" : 0.15, # spectrogram
"axes_emg" : 0.10, # NEW: EMG, same height as sig
"axes_sig" : 0.10, # significance trace
"axes_pcs" : 0.52, # PC/trace panel
"axes_legend": 0.04 # legend bar
}
pad = 0.01 # small vertical padding between rows
# cumulative bottom edges
y_bottom = 1.0
axes_dict = {}
for name, h in lay.items():
y_bottom -= h
ax = plt.axes([0.10, y_bottom, 0.80, h], sharex=axes_dict.get("axes_brs"))
axes_dict[name] = ax
y_bottom -= pad
# separate colour bar next to spectrogram
axes_dict["axes_cbar"] = plt.axes([0.91, # x‑pos
1.0 - lay["axes_brs"] - pad - lay["axes_spec"],
0.04, # width
lay["axes_spec"]])
# unpack for backward compatibility
axes_brs = axes_dict["axes_brs"]
axes_spec = axes_dict["axes_spec"]
axes_emg = axes_dict["axes_emg"]
axes_sig = axes_dict["axes_sig"]
axes_pcs = axes_dict["axes_pcs"]
axes_legend = axes_dict["axes_legend"]
axes_cbar = axes_dict["axes_cbar"]
# show brainstate
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
sleepy._despine_axes(axes_brs)
# show EEG spectrogram
# calculate median for choosing right saturation for heatmap
med = np.median(SP.max(axis=0))
if len(vm) == 0:
vm = [0, med*2.0]
ifreq = np.where(freq <= fmax)[0]
im = axes_spec.pcolorfast(t, freq[ifreq], SP[ifreq, istart:iend], cmap='jet', vmin=vm[0], vmax=vm[1])
axes_spec.axis('tight')
axes_spec.set_xticklabels([])
axes_spec.set_xticks([])
axes_spec.spines["bottom"].set_visible(False)
axes_spec.set_ylabel('Freq (Hz)')
sleepy.box_off(axes_spec)
axes_spec.set_xlim([t[0], t[-1]])
# colorbar for EEG spectrogram
cb = plt.colorbar(im, ax=axes_cbar, pad=0.0, aspect=10.0)
if pnorm_spec:
cb.set_label('Norm. power')
else:
cb.set_label('Power ($\mathrm{\mu}$V$^2$s)')
#if len(cb_ticks) > 0:
# cb.set_ticks(cb_ticks)
axes_cbar.set_alpha(0.0)
sleepy._despine_axes(axes_cbar)
#im = axes_spec.pcolorfast(t, freq[ifreq], SP[ifreq, istart:iend], cmap='jet', vmin=vm[0], vmax=vm[1])
if not(no_emg):
axes_emg.plot(t, p_mu, color='k')
axes_emg.set_xlim([t[0], t[-1]])
axes_emg.set_ylabel('EMG ampl. $\mathrm{\mu}$V')
axes_emg.spines["top"].set_visible(False)
axes_emg.spines["right"].set_visible(False)
axes_emg.spines["bottom"].set_visible(False)
axes_emg.axes.get_xaxis().set_visible(False)
# show sigmapower
axes_sig.plot(t, sigma_pow, color='gray')
axes_sig.set_xlim([t[0], t[-1]])
axes_sig.set_ylabel('$\mathrm{\sigma}$ power')
axes_sig.spines["top"].set_visible(False)
axes_sig.spines["right"].set_visible(False)
axes_sig.spines["bottom"].set_visible(False)
axes_sig.axes.get_xaxis().set_visible(False)
a = np.percentile(sigma_pow, 99.5)
plt.ylim([0, a+a*0.1])
# show PCs
for i in range(ndim):
PC[i,:] = PC[i,:] - np.min(PC[i,:])
pc_max = []
for i in range(ndim):
p = np.max(PC[i,:])
pc_max.append(p)
mmax = np.max(np.array(pc_max))
pos = [0]
for i in range(1,ndim):
if reverse_pcs:
PC[i,:] = PC[i,:] + i*mmax
pos.append(i*mmax)
else:
PC[i,:] = PC[i,:] - i*mmax
pos.append(-i*mmax)
# axes for PCs
# colors for PCs
cmap = sns.color_palette("husl", ndim)
for i in range(ndim):
axes_pcs.plot(t, PC[i,:], c=cmap[i])
axes_pcs.text(t[-1], pos[i], 'PC%d' % (i+1), fontsize=14, color=cmap[i])
plt.xlim((t[0], t[-1]))
axes_pcs.spines["left"].set_visible(False)
axes_pcs.spines["right"].set_visible(False)
axes_pcs.axes.get_yaxis().set_visible(False)
sleepy._despine_axes(axes_pcs)
if len(zoomin) > 0:
zoomin = [int(z/dt) for z in zoomin]
for z in zoomin:
pos = np.array(pos)
plt.plot([t[z], t[z]], [mmax, -(ndim-1)*mmax], 'k--')
# axes for time legend
axes_legend.plot([0, tlegend], [1, 1], lw=2, color='k')
axes_legend.set_ylim([-1, 1])
axes_legend.text(0, -2, '%d s' % tlegend, verticalalignment='bottom', horizontalalignment='left')
plt.xlim((t[0], t[-1]))
sleepy._despine_axes(axes_legend)
return PC_orig, M, sigma_pow
def sleep_subspaces(units, mouse, ndim=2, nsmooth=0, ma_thr=10, ma_rem_exception=False,
pzscore=False, detrend=False, pplot=True, coords=[], traj_mode='full', trig_state=1,
local_rotation=True, pspec=True,
kcuts=[], proj_3d=False, config_file=''):
"""
Calculate (and plot) PCA separately for each state (REM, Wake, NREM).
Parameters
----------
units : np.DataFrame or []
If [], use 1 ms spike trains to calculate firing rates.
mouse : TYPE
DESCRIPTION.
ndim : int, optional
Number of PCs to keep for dimensionality reduction.
The default is 2.
nsmooth : TYPE, optional
DESCRIPTION. The default is 0.
ma_thr : float, optional
Wake sequences <= $ma_thr s are interpreted as NREM (3). The default is 10.
ma_rem_exception : bool, optional
If True, then the MA rule does not apply for wake episodes directly following REM.
The default is False.
pzscore : TYPE, optional
DESCRIPTION. The default is False.
pplot : TYPE, optional
DESCRIPTION. The default is True.
coords : list, optional
Specific the two PCs that should be shown on the x and y-axis.
PC1 corresponds to "0".
The default is [].
traj_mode : string, optional
If 'full', show complete trajectory.
If 'trig', only show trajectories for the specific state transitions (-> trig_state)
The default is 'full'.
trig_state : int, optional
1,2, or 3. If traj_mode == 'trig', only show the the transitions to state $trig_state
instead of the full trajectories across the whole recording session
The default is 1.
local_rotation : bool, optional
If True, rotate the ndim dimensional space (by performing another PCA) and
use the resulting first two dimensions to get the PCs.
If True, the parameter @coords has no effect.
The default is True.
pspec : TYPE, optional
DESCRIPTION. The default is True.
kcuts : TYPE, optional
DESCRIPTION. The default is [].
proj_3d : bool, optional
If True, plot 3D subspaces
config_file : str
file name of mouse configuration file, as loaded by &load_config()
Returns
-------
pc_dict : TYPE
DESCRIPTION.
vh_dict : TYPE
DESCRIPTION.
idx_dict : TYPE
DESCRIPTION.
"""
dt = 2.5
NDOWN = 500
NUP = int(dt / (0.001 * NDOWN))
fine_scale = False
if len(config_file) == 0:
config_file = 'mouse_config.txt'
path = load_config(config_file)[mouse]['SL_PATH']
ppath, file = os.path.split(path)
M = sleepy.load_stateidx(ppath, file)[0]
###########################################################################
if len(units) == 0:
fine_scale = True
tr_path = load_config(config_file)[mouse]['TR_PATH']
units = np.load(os.path.join(tr_path,'1k_train.npz'))
unitIDs = [unit for unit in list(units.keys()) if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
dt = dt / NUP
M = upsample_mx(M, NUP)
nhypno = int(np.min((len(M), units[unitIDs[0]].shape[0]/NDOWN)))
else:
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nhypno = np.min((len(M), units.shape[0]))
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# OLD:
# if len(kcut) > 0:
# kidx = np.arange(int(kcut[0]/dt), int(kcut[-1]/dt))
# tidx = np.setdiff1d(tidx, kidx)
# M = M[tidx]
# nhypno = len(tidx)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
#print(len(tidx))
nsample = nhypno # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
#@tidx are the indices we're further considereing.
print('Starting downsampling, smoothing and z-scoring')
if fine_scale:
fr_file = os.path.join(tr_path, 'fr_fine_ndown%d.mat' % NDOWN)
if not os.path.isfile(fr_file):
for i,unit in enumerate(unitIDs):
tmp = sleepy.downsample_vec(np.array(units[unit]), NDOWN)
R[i,:] = tmp[tidx]
so.savemat(fr_file, {'R':R, 'ndown':NDOWN})
else:
R = so.loadmat(fr_file, squeeze_me=True)['R']
for i,unit in enumerate(unitIDs):
tmp = R[i,:]
tmp = sleepy.smooth_data(tmp, nsmooth)
if pzscore:
R[i,:] = (tmp[tidx] - tmp[tidx].mean()) / tmp[tidx].std()
else:
R[i,:] = tmp[tidx]
if not fine_scale:
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]), nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
#R[i,:] = (tmp[tidx] - tmp[tidx].mean()) / tmp[tidx].std()
R[i,:] = (tmp - tmp.mean()) / tmp.std()
else:
#R[i,:] = tmp[tidx]
R[i,:] = tmp
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 4
# all REM sequences
ridx = sleepy.get_sequences(np.where(M==1)[0])
# find long wake blocks:
widx = sleepy.get_sequences(np.where(M==2)[0])
# all NREM sequences
nidx = sleepy.get_sequences(np.where(M>=3)[0])
midx = sleepy.get_sequences(np.where(M==4)[0])
ridx = [list(a) for a in ridx]
widx = [list(a) for a in widx]
nidx = [list(a) for a in nidx]
midx = [list(a) for a in midx]
ridx = sum(ridx, [])
widx = sum(widx, [])
nidx = sum(nidx, [])
midx = sum(midx, [])
idx_dict = {'REM':ridx, 'Wake':widx, 'NREM':nidx, 'MA':midx}
pc_dict = {}
vh_dict = {}
labels = ['REM', 'Wake', 'NREM', 'MA']
for idx,label in zip([ridx, widx, nidx], labels):
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R[:,idx].T / np.sqrt(nsample-1)
# make sure the the columns of Y are mean-zero
for i in range(nvar):
Y[:,i] = Y[:,i] - Y[:,i].mean()
# Note the PCs are orthogonal to each other, but
# only for the segments for which we calculate PCA
# and only if we center the firing rates across these
# segments!
R[i,:] = R[i,:] - R[i,idx].mean()
#R[i,:] = R[i,:] - R[i,:].mean()
# SVD
U,S,Vh = scipy.linalg.svd(Y)
# each row in Vh is an eigenvector of the COV matrix:
PC = np.dot(Vh, R)[0:ndim,:]
if local_rotation:
if not proj_3d:
PC = pca(PC.copy().T, dims=2)[0].T
coords = [0,1]
else:
PC = pca(PC.copy().T, dims=3)[0].T
coords = [0,1,2]
pc_dict[label] = PC
vh_dict[label] = Vh[0:ndim,:]
if coords == []:
coords = list(range(ndim))
if pplot:
sleepy.set_fontsize(12)
plt.figure(figsize=(12,10))
t = np.arange(0, len(M))*dt
if pspec:
# load spectrogram
tmp = so.loadmat(os.path.join(ppath, file, 'sp_%s.mat'%file), squeeze_me=True)
SP = tmp['SP']
freq = tmp['freq']
ifreq = np.where(freq < 30)[0]
axes_spec = plt.axes([0.4, 0.85, 0.55, 0.08])
axes_spec.pcolorfast(t, freq[ifreq], SP[ifreq,:], vmin=0, vmax=2000, cmap='jet')
sleepy._despine_axes(axes_spec)
axes_brs = plt.axes([0.4, 0.8, 0.55, 0.025], sharex=axes_spec)
else:
axes_brs = plt.axes([0.4, 0.8, 0.55, 0.025])
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
axes_brs.axes.get_xaxis().set_visible(False)
axes_brs.axes.get_yaxis().set_visible(False)
axes_brs.spines["top"].set_visible(False)
axes_brs.spines["right"].set_visible(False)
axes_brs.spines["bottom"].set_visible(False)
axes_brs.spines["left"].set_visible(False)
ax1 = plt.axes([0.4, 0.6, 0.55, 0.2], sharex=axes_brs)
sns.despine()
ax2 = plt.axes([0.4, 0.35, 0.55, 0.2], sharex=ax1)
sns.despine()
ax3 = plt.axes([0.4, 0.1, 0.55, 0.2], sharex=ax1)
sns.despine()
axes = [ax1, ax2, ax3]
i = 0
for ax, label in zip(axes, pc_dict):
ax.plot(t, pc_dict[label][coords].T)
ax.set_xlim([t[0], t[-1]])
if i < 2:
#ax.set_xticklabels([])
pass
else:
ax.set_xlabel('Time (s)')
i += 1
if not proj_3d:
ax1 = plt.axes([0.1, 0.6, 0.2, 0.2])
sns.despine()
ax2 = plt.axes([0.1, 0.35, 0.2, 0.2])
sns.despine()
ax3 = plt.axes([0.1, 0.1, 0.2, 0.2])
sns.despine()
else:
ax1 = plt.axes([0.1, 0.6, 0.2, 0.2], projection='3d')
sns.despine()
ax2 = plt.axes([0.1, 0.35, 0.2, 0.2], projection='3d')
sns.despine()
ax3 = plt.axes([0.1, 0.1, 0.2, 0.2], projection='3d')
sns.despine()
axes = [ax1, ax2, ax3]
clrs = [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8], [1, 0.2, 0.2]]
if traj_mode == 'full':
for ax, label in zip(axes, pc_dict):
PC = pc_dict[label]
p = M[0]
k = 0
kold = k
while k < len(M)-1:
while M[k] == p and k < len(M)-1:
k+=1
if not proj_3d:
ax.plot(PC[coords[0],kold:k], PC[coords[1], kold:k], color=clrs[int(p)], lw=0.5)
else:
ax.plot3D(PC[coords[0],kold:k], PC[coords[1], kold:k], PC[coords[2], kold:k], color=clrs[int(p)], lw=0.5)
p = M[k]
kold = k-1
if local_rotation:
ax.set_xlabel("PC1'")
ax.set_ylabel("PC2'")
else:
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
else:
for ax, label in zip(axes, pc_dict):
PC = pc_dict[label]
plot_trajectories(PC, M, 300, 0, dt=dt, min_dur=20, istate=trig_state, pre_state=2,
state_num=[], ma_thr=20, kcuts=(), coords=coords, ax=ax, lw=0.5)
if local_rotation:
ax.set_xlabel("PC1'")
ax.set_ylabel("PC2'")
else:
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
return pc_dict, vh_dict, idx_dict
def subspace_mixing_mx(units, mouse, ndim=2, nsmooth=0, ma_thr=10, ma_rem_exception=False,
pzscore=False, detrend=False, pplot=True, coords=[], traj_mode='full', trig_state=1,
pspec=True, kcuts=[], proj_3d=False, config_file=''):
"""
Determine coordindates systems for REM, NREM, and Wake subspace.
Then, project the firing rate vectors for REM, NREM, and Wake into each subspace
and determine the variance of the projected data.
See also function &sleep_subspaces()
Calculate (and plot) PCA separately for each state (REM, Wake, NREM).
Parameters
----------
units : TYPE
DESCRIPTION.
mouse : TYPE
DESCRIPTION.
ndim : int, optional
Number of PCs to keep for dimensionality reduction.
The default is 2.
nsmooth : float, optional
Smooth firing rates using sleepy.smooth_data(data, nsmooth). The default is 0.
ma_thr : float, optional
Wake sequences <= $ma_thr s are interpreted as NREM (3). The default is 10.
ma_rem_exception : bool, optional
If True, then the MA rule does not apply for wake episodes directly following REM.
The default is False.
pzscore : TYPE, optional
DESCRIPTION. The default is False.
pplot : TYPE, optional
DESCRIPTION. The default is True.
coords : list, optional
Specific the two PCs that should be shown on the x and y-axis.
PC1 corresponds to "0".
The default is [].
traj_mode : string, optional
If 'full', show complete trajectory.
If 'trig', only show trajectories for the specific state transitions (-> trig_state)
The default is 'full'.
trig_state : int, optional
1,2, or 3. If traj_mode == 'trig', only show the the transitions to state $trig_state
instead of the full trajectories across the whole recording session
The default is 1.
pspec : TYPE, optional
DESCRIPTION. The default is True.
kcuts : TYPE, optional
DESCRIPTION. The default is [].
proj_3d : bool, optional
If True, plot 3D subspaces
config_file : str
file name of mouse configuration file, as loaded by &load_config()
Returns
-------
pc_dict : dict
dict: state --> PCs.
vh_dict : TYPE
DESCRIPTION.
idx_dict : TYPE
DESCRIPTION.
"""
dt = 2.5
if len(config_file) == 0:
config_file = 'mouse_config.txt'
path = load_config(config_file)[mouse]['SL_PATH']
ppath, file = os.path.split(path)
M = sleepy.load_stateidx(ppath, file)[0]
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nhypno = np.min((len(M), units.shape[0]))
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
#@tidx are the indices we're further considereing.
###########################################################################
nsample = nhypno # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
# R: each row is one neuron with firing rates for len(tidx) time points
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]), nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
R[i,:] = (tmp - tmp.mean()) / tmp.std()
else:
R[i,:] = tmp
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
# all REM sequences
ridx = sleepy.get_sequences(np.where(M==1)[0])
# find long wake blocks:
widx = sleepy.get_sequences(np.where(M==2)[0])
# all NREM sequences
nidx = sleepy.get_sequences(np.where(M>=3)[0])
midx = sleepy.get_sequences(np.where(M==4)[0])
ridx = [list(a) for a in ridx]
widx = [list(a) for a in widx]
nidx = [list(a) for a in nidx]
midx = [list(a) for a in midx]
ridx = sum(ridx, [])
widx = sum(widx, [])
nidx = sum(nidx, [])
midx = sum(midx, [])
idx_dict = {'REM':ridx, 'Wake':widx, 'NREM':nidx, 'MA':midx}
pc_dict = {}
vh_dict = {}
cc = np.sqrt(nsample-1)
cc = 1
labels = ['REM', 'Wake', 'NREM']
for idx,label in zip([ridx, widx, nidx], labels):
Rsub = R.copy()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R[:,idx].T / cc
# make sure the the columns of Y are mean-zero
for i in range(nvar):
Y[:,i] = Y[:,i] - Y[:,i].mean()
# Note the PCs are orthogonal to each other, but
# only for the segments for which we calculate PCA
# and only if we center the firing rates across these
# segments!
Rsub[i,:] = Rsub[i,:] - Rsub[i,idx].mean()
# SVD
U,S,Vh = scipy.linalg.svd(Y)
# each row in Vh is an eigenvector of the COV matrix:
PC = np.dot(Vh, Rsub)[0:ndim,:]
pc_dict[label] = PC
vh_dict[label] = Vh[0:ndim,:] * cc
MX = np.zeros((len(labels), len(labels)))
# Build mixing matrix
for i in range(len(labels)):
PCsub = pc_dict[labels[i]]
# For example take NREM subspace and project
# NREM, Wake and REM firing rates into this space, and
# then calculate the total variance of NREM, Wake, and REM
# neurons within this space
for j in range(0, len(labels)):
state_idx = idx_dict[labels[j]]
A = PCsub[:,state_idx]
MX[i,j] = np.trace(np.cov(A))
if coords == []:
coords = list(range(ndim))
if pplot:
sleepy.set_fontsize(12)
plt.figure(figsize=(12,10))
t = np.arange(0, len(M))*dt
if pspec:
# load spectrogram
tmp = so.loadmat(os.path.join(ppath, file, 'sp_%s.mat'%file), squeeze_me=True)
SP = tmp['SP']
freq = tmp['freq']
ifreq = np.where(freq < 30)[0]
axes_spec = plt.axes([0.4, 0.85, 0.55, 0.08])
axes_spec.pcolorfast(t, freq[ifreq], SP[ifreq,:], vmin=0, vmax=2000, cmap='jet')
sleepy._despine_axes(axes_spec)
axes_brs = plt.axes([0.4, 0.8, 0.55, 0.025], sharex=axes_spec)
else:
axes_brs = plt.axes([0.4, 0.8, 0.55, 0.025])
cmap = plt.cm.jet
my_map = cmap.from_list('brs', [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8]], 4)
tmp = axes_brs.pcolorfast(t, [0, 1], np.array([M]), vmin=0, vmax=3)
tmp.set_cmap(my_map)
axes_brs.axis('tight')
axes_brs.axes.get_xaxis().set_visible(False)
axes_brs.axes.get_yaxis().set_visible(False)
axes_brs.spines["top"].set_visible(False)
axes_brs.spines["right"].set_visible(False)
axes_brs.spines["bottom"].set_visible(False)
axes_brs.spines["left"].set_visible(False)
ax1 = plt.axes([0.4, 0.6, 0.55, 0.2], sharex=axes_brs)
sns.despine()
ax2 = plt.axes([0.4, 0.35, 0.55, 0.2], sharex=ax1)
sns.despine()
ax3 = plt.axes([0.4, 0.1, 0.55, 0.2], sharex=ax1)
sns.despine()
axes = [ax1, ax2, ax3]
i = 0
for ax, label in zip(axes, pc_dict):
ax.plot(t, pc_dict[label][coords,:].T)
ax.set_xlim([t[0], t[-1]])
if i < 2:
#ax.set_xticklabels([])
pass
else:
ax.set_xlabel('Time (s)')
i += 1
if not proj_3d:
ax1 = plt.axes([0.1, 0.6, 0.2, 0.2])
sns.despine()
ax2 = plt.axes([0.1, 0.35, 0.2, 0.2])
sns.despine()
ax3 = plt.axes([0.1, 0.1, 0.2, 0.2])
sns.despine()
else:
ax1 = plt.axes([0.1, 0.6, 0.2, 0.2], projection='3d')
sns.despine()
ax2 = plt.axes([0.1, 0.35, 0.2, 0.2], projection='3d')
sns.despine()
ax3 = plt.axes([0.1, 0.1, 0.2, 0.2], projection='3d')
sns.despine()
axes = [ax1, ax2, ax3]
clrs = [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8], [1, 0.2, 0.2]]
if traj_mode == 'full':
for ax, label in zip(axes, pc_dict):
PC = pc_dict[label]
p = M[0]
k = 0
kold = k
while k < len(M)-1:
while M[k] == p and k < len(M)-1:
k+=1
if not proj_3d:
ax.plot(PC[coords[0],kold:k], PC[coords[1], kold:k], color=clrs[int(p)], lw=0.5)
else:
ax.plot3D(PC[coords[0],kold:k], PC[coords[1], kold:k], PC[coords[2], kold:k], color=clrs[int(p)], lw=0.5)
p = M[k]
kold = k-1
else:
for ax, label in zip(axes, pc_dict):
PC = pc_dict[label]
plot_trajectories(PC, M, 300, 0, dt=dt, min_dur=20, istate=trig_state, pre_state=2,
state_num=[], ma_thr=20, kcuts=(), coords=coords, ax=ax, lw=0.5)
return pc_dict, vh_dict, idx_dict, MX
def pc_reconstruction(units, cell_info, ndim=3, nsmooth=0, pnorm=True, pzscore=False,
pc_sign=[], dt=2.5):
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nsample = units.shape[0] # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
# OLD:
i = 0
for unit in unitIDs:
R[i,:] = sleepy.smooth_data(np.array(units[unit]),nsmooth)
if pzscore:
R[i,:] = (R[i,:] - R[i,:].mean()) / R[i,:].std()
i += 1
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
for i in range(nvar):
R[i,:] = R[i,:] - R[i,:].mean()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
Y = R.T / np.sqrt(nsample-1)
U,S,Vh = scipy.linalg.svd(Y)
PC = np.dot(Vh, R)[0:ndim,:]
plt.figure()
plt.plot(PC.T)
# reconstruct neural responses with the first $ndim eigenvectors
SM = np.zeros((nsample, ndim))
for i in range(ndim):
SM[i,i] = S[i]
# That's the reconstruction:
Yhat = (np.dot(U[:,0:ndim], np.dot(SM[0:ndim,:], Vh[0:ndim,:]))) * np.sqrt(nsample-1)
# Alternative:
# neurons x dim * [dim x neurons * neurons x time]
# Yhat2 = np.dot(Vh[0:ndim,:].T, np.dot(Vh[0:ndim,:], R)).T
C = Vh[0:ndim,:].T
if pnorm:
for i in range(C.shape[0]):
C[i,:] = C[i,:] / np.sqrt(np.sum(C[i,:]**2))
if len(pc_sign) > 0:
i = 0
for s in pc_sign:
C[:,i] = C[:,i] * s
i += 1
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
#df['brain_region'] = list(cell_info[cell_info.ID.isin(unitIDs)]['brain_region'])
df['brain_region'] = [cell_info[cell_info.ID == i].brain_region.iloc[0] for i in unitIDs]
plt.figure()
plt.subplot(211)
sns.histplot(data=df, x='brain_region', y='c1')
plt.subplot(212)
sns.histplot(data=df, x='brain_region', y='c2')
return C, unitIDs, df, Yhat
def pc_reconstruction2(units, cell_info, time_idx=[], ndim=3, nsmooth=0, detrend=False, pnorm=True, pzscore=False,
pc_sign=[], scale_pc=False, pearson_scaling=False,
dim_reconstr=[], dt=2.5, kcuts=[], pearson=False, pplot=True, sign_plot=True):
"""
Use whole time axis for smoothing and z-scoring.
Calculate SVD only using time points in @time_idx and reconstruct firing rates
only for timepoints in @time_idx.
When calculating the coefficients for each PC this function
takes only time points in time_idx into account!
Note on SVD:
Assume Y is a matrix with each column corresponding to the firing rate vector
of one recorded unit (Y ~ time bins x units).
Assume
U,S,Vh = scipy.linalg.svd(Y)
is the SVD of matrix Y
Then, PC = U * S are the PCs;
PC[i,:] is the i-th PC
Vh[i,:] are the coefficients of each neurons for PCi
The firing rates of unit fr_i can be reconstructed using,
fr_i = PC1 * c[0,i] + PC2 * c2[1,i] + ...
which we rewrite as
fr_i = PC1 * c1_i + PC2 * c2_i + ...
Parameters
----------
units : pd.DataFrame
Firing rates (time_bins x num_units). Each column is a unit with unit ID as column name.
cell_info : pd.DataFrame
Cell metadata. Must contain columns:
- 'ID': Unit identifiers (matching units.columns)
- 'brain_region': Brain region assignment for each unit
time_idx : list or np.ndarray, optional
Time indices to use for SVD computation. If [], uses all available timepoints.
Default is [].
ndim : int, optional
Number of principal components to extract and use for reconstruction. Default is 3.
nsmooth : float, optional
Gaussian kernel smoothing factor (σ) applied to firing rates before analysis. Default is 0.
detrend : bool, optional
If True, detrend firing rates before analysis. Default is False.
pnorm : bool, optional
If True, normalize PC coefficients (L2 norm). Default is True.
pzscore : bool, optional
If True, z-score firing rates after smoothing. Default is False.
pc_sign : list or np.ndarray, optional
Sign correction vector: PC_i is multiplied by pc_sign[i-1]. Must have length ndim.
Use to ensure consistent PC orientation across analyses. Default is [].
scale_pc : bool, optional
If True, scale PC coefficients by PC variance. Default is False.
pearson_scaling : bool, optional
If True, apply Pearson correlation-based scaling. Default is False.
dim_reconstr : list, optional
Specific dimensions to reconstruct. If [], reconstructs all ndim dimensions. Default is [].
dt : float, optional
Sampling period (in seconds) for converting between time and sample indices. Default is 2.5.
kcuts : list of tuples, optional
Timepoints to exclude. Each tuple (start, end) specifies interval in seconds to discard.
Default is [].
pearson : bool, optional
If True, calculate Pearson correlation between each PC and firing rates, returning
r and p values in the output DataFrame. Default is False.
pplot : bool, optional
If True, display plots of principal components. Default is True.
sign_plot : bool, optional
If True, plots show PCs after multiplying by pc_sign. Default is True.
Returns
-------
C : np.ndarray with shape (num_units, ndim)
PC coefficients for each unit. C[i, j] is the coefficient of unit i for PC j.
df : pd.DataFrame
Coefficient summary table with columns:
- 'c1', 'c2', ..., 'cn': PC coefficients for each unit
- 'ID': Unit identifier
- 'brain_region': Brain region assignment (from cell_info)
- 'r1', 'r2', ..., 'rn': Pearson r values (if pearson=True)
- 'p1', 'p2', ..., 'pn': Pearson p-values (if pearson=True)
units : pd.DataFrame
Original firing rates (after smoothing/detrending if applicable).
units_hat : pd.DataFrame
Reconstructed firing rates: all_vars @ C.T. Shape (num_timepoints, num_units).
"""
nhypno = units.shape[0]
tidx = np.arange(0, nhypno)
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
nhypno = len(tidx)
###########################################################################
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nsample = units.shape[0] # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nhypno)) # dimensions: number of units x time bins
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]),nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
R[i,:] = (tmp - tmp.mean()) / tmp.std()
else:
R[i,:] = tmp
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
R = R[:,time_idx]
for i in range(nvar):
R[i,:] = R[i,:] - R[i,:].mean()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
# NEW: 05/14/25!!!
nsample = R.shape[1]
Y = R.T / np.sqrt(nsample-1)
U,S,Vh = scipy.linalg.svd(Y)
PC = np.dot(Vh, R)[0:ndim,:]
if scale_pc:
PC = PC / np.sqrt(nsample - 1)
## double check code using chat-gpt:
# n_units = len(unitIDs)
# R = np.empty((n_units, len(tidx)))
# for i, unit in enumerate(unitIDs):
# x = sleepy.smooth_data(np.asarray(units[unit]), nsmooth)[tidx]
# if detrend:
# x = scipy.signal.detrend(x)
# R[i] = (x - x.mean()) / x.std(ddof=1) if pzscore else x
# # --- restrict to analysis window ------------------------------------------
# R = R[:, time_idx] # if time_idx == tidx you can drop one of them
# nsample = R.shape[1]
# # --- mean-centre (skipped automatically if pzscore) ------------------------
# R -= R.mean(axis=1, keepdims=True)
# # --- PCA via SVD -----------------------------------------------------------
# Y = R.T / np.sqrt(nsample - 1)
# U, S, Vh = scipy.linalg.svd(Y, full_matrices=False)
# scores = (S[:ndim, None] * U[:, :ndim].T) # PC time series
# loadings = Vh[:ndim] # neuron weights
# eigvals = S[:ndim]**2 # variance of each PC
# expl_var = eigvals / eigvals.sum()
if len(pc_sign) == []:
pc_sign = np.ones((ndim,))
if sign_plot:
plt.figure()
for i in range(ndim):
if i < len(pc_sign):
plt.plot(PC[i,:]*pc_sign[i])
else:
plt.plot(PC[i,:])
# reconstruct neural responses with the first $ndim eigenvectors
SM = np.zeros((nsample, ndim))
for i in range(ndim):
SM[i,i] = S[i]
if len(dim_reconstr) == 0:
dim_reconstr = range(0, ndim)
# That's the reconstruction:
Yhat = (np.dot(U[:,dim_reconstr], np.dot(SM[dim_reconstr,:], Vh[dim_reconstr,:]))) * np.sqrt(nsample-1)
# Alternative:
# neurons x dim * [dim x neurons * neurons x time]
# Yhat2 = np.dot(Vh[0:ndim,:].T, np.dot(Vh[0:ndim,:], R)).T
# 11/08/23 added the factor '* np.sqrt(nsample-1)'
C = Vh[0:ndim,:].T #* np.sqrt(nsample-1)
if scale_pc:
C = C * np.sqrt(nsample-1)
if pnorm:
for i in range(C.shape[0]):
C[i,:] = C[i,:] / np.sqrt(np.sum(C[i,:]**2))
if pearson_scaling:
C = C * S[:ndim]
if len(pc_sign) > 0:
i = 0
for s in pc_sign:
C[:,i] = C[:,i] * s
i += 1
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
if pplot:
df['brain_region'] = [cell_info[cell_info.ID == i].brain_region.iloc[0] for i in unitIDs]
plt.figure()
plt.subplot(211)
sns.histplot(data=df, x='brain_region', y='c1')
plt.subplot(212)
sns.histplot(data=df, x='brain_region', y='c2')
#convert Yhat to pd.DataFrame with rows of Yhat as columns
units_hat = pd.DataFrame(data=Yhat, columns=unitIDs)
units = pd.DataFrame(data=Y*np.sqrt(nsample-1), columns=unitIDs)
if pearson:
for j in range(ndim):
p = []
cc = []
for i,ID in enumerate(unitIDs):
r = R[i,:]
res = scipy.stats.pearsonr(r, PC[j,:])
p.append(res.pvalue)
cc.append(res.statistic)
df['r' + str(j+1)] = cc
df['p' + str(j+1)] = p
return C, df, units, units_hat
def partialpc_reconstruction(units, cell_info, time_idx, neuron_ids=[], ndim=3, nsmooth=0,
pnorm=True, pzscore=False, pc_sign=[], fit_timeidx=True, pplot=True):
"""
NOTE: In its current implementation, the function fits the PC coefficients
to the complete time axis, although the principal axes have been fit using
a (potentially smaller) time range.
Parameters
----------
units : TYPE
DESCRIPTION.
cell_info : TYPE
DESCRIPTION.
time_idx : TYPE
DESCRIPTION.
neuron_ids : TYPE, optional
DESCRIPTION. The default is [].
ndim : TYPE, optional
DESCRIPTION. The default is 3.
nsmooth : TYPE, optional
DESCRIPTION. The default is 0.
pnorm : TYPE, optional
DESCRIPTION. The default is True.
pzscore : TYPE, optional
DESCRIPTION. The default is False.
pc_sign : TYPE, optional
DESCRIPTION. The default is [].
fit_timeidx: boolean, optional
If True, take for fit of firing rates using the PCs only time points in time_idx into account
Returns
-------
C : TYPE
DESCRIPTION.
PC : TYPE
DESCRIPTION.
df : TYPE
DESCRIPTION.
units: pd.DataFrame
smoothed and centered firing rates
units_hat: pd.DataFrame
fitted firing rates using $ndim PCs
"""
unitIDs = [unit for unit in units.columns if '_' in unit]
unitIDs = [unit for unit in unitIDs if re.split('_', unit)[1] == 'good']
nsample = units.shape[0] # number of time points
nvar = len(unitIDs) # number of units
R = np.zeros((nvar, nsample))
#all_units = False
if neuron_ids == []:
neuron_ids = unitIDs
#all_units = True
i = 0
neuron_idx = []
for unit in unitIDs:
R[i,:] = sleepy.smooth_data(np.array(units[unit]),nsmooth)
if pzscore:
R[i,:] = (R[i,:] - R[i,:].mean()) / R[i,:].std()
if unit in neuron_ids:
neuron_idx.append(i)
i += 1
# first, for each varible (dimension), we first need to remove the mean:
# mean-zero rows:
for i in range(nvar):
R[i,:] = R[i,:] - R[i,:].mean()
# divide by sqrt(nsample - 1) to make SVD equivalent to PCA
R2 = R[neuron_idx,:]
Y = R2[:,time_idx].T / np.sqrt(len(time_idx)-1)
# make sure the the columns of Y are mean-zero
for i in range(Y.shape[1]):
Y[:,i] = Y[:,i] - Y[:,i].mean()
U,S,Vh = scipy.linalg.svd(Y)
# that's the projections of all neuron_idx neurons (rows) into the PCA space,
# giving us the principal components
PC = np.dot(Vh, R2)[0:ndim,:]
# Now, for ALL neurons and time points the to optimal coefficients to
# reconstruct the original firing rates using the $ndim PCs
A = np.ones((PC.shape[1], ndim+1))
A[:,0:ndim] = PC.T
Rhat = np.zeros(R.shape)
C = np.zeros((nvar, ndim+1))
for i in range(nvar):
r = R[i,:]
if fit_timeidx:
w = np.linalg.lstsq(A[time_idx], r[time_idx])[0]
else:
w = np.linalg.lstsq(A, r)[0]
C[i,:] = w
# reconstruction of units using ndeim PCs
rhat = np.dot(A, w)
Rhat[i,:] = rhat
C = C[:,0:ndim]
#C = Vh[0:ndim,:].T
if pnorm:
for i in range(C.shape[0]):
C[i,:] = C[i,:] / np.sqrt(np.sum(C[i,:]**2))
if len(pc_sign) > 0:
for i,s in enumerate(pc_sign):
C[:,i] = C[:,i] * s
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
labels = ['c'+str(i) for i in range(1, ndim+1)]
df = pd.DataFrame(data=C, columns=labels)
df['ID'] = unitIDs
#df['brain_region'] = list(cell_info[cell_info.ID.isin(unitIDs)]['brain_region'])
df['brain_region'] = [cell_info[cell_info.ID == i].brain_region.iloc[0] for i in unitIDs]
### FIGURE ################################################################
if pplot:
plt.figure()
plt.subplot(211)
sns.histplot(data=df, x='brain_region', y='c1')
plt.subplot(212)
sns.histplot(data=df, x='brain_region', y='c2')
# reconstruct neural responses with the first $ndim eigenvectors
#SM = np.zeros((nsample, ndim))
#for i in range(ndim):
# SM[i,i] = S[i]
# That's the reconstruction:
#Yhat = (np.dot(U[:,0:ndim], np.dot(SM[0:ndim,:], Vh[0:ndim,:]))) * np.sqrt(nsample-1)
# UNDER CONSTRUCTION
units_hat = pd.DataFrame(data=Rhat.T, columns=unitIDs)
units = pd.DataFrame(data=R.T, columns=unitIDs)
return C, PC, df, units, units_hat
def optimal_direction(dfc, pref_direction, thr=0, pplot=True, ax='', brain_region=False):
"""
Parameters
----------
dfc : pd.DataFrame
DataFrame with coefficients, c1, ... c_n (loadings) for PCs.
Columns: 'c1', ... 'c_n', 'ID', 'brain_region'
pref_direction : TYPE
DESCRIPTION.
thr : TYPE
DESCRIPTION.
Returns
-------
df2 : TYPE
DESCRIPTION.
"""
pref_direction = pref_direction / scipy.linalg.norm(pref_direction)
cols = ['c'+str(i+1) for i in range(pref_direction.shape[0])]
data = []
for index, row in dfc.iterrows():
vec = np.array(row[cols])
# normalize each vector
vec = vec / scipy.linalg.norm(vec)
d = np.dot(vec, pref_direction)
if brain_region:
data += [[d, row['ID'], row['brain_region']]]
else:
data += [[d, row['ID']]]
if brain_region:
df2 = pd.DataFrame(data=data, columns=['proj', 'ID', 'brain_region'])
else:
df2 = pd.DataFrame(data=data, columns=['proj', 'ID'])
if thr > 0:
dfs = dfc[(df2.proj > thr)]
else:
dfs = dfc[(df2.proj < thr)]
if pplot:
if ax == '':
plt.figure()
ax = plt.axes([0.2, 0.2, 0.7, 0.7])
#ax.axis('equal')
ax.axhline(0, linestyle='--', color='k') # horizontal lines
ax.axvline(0, linestyle='--', color='k') # vertical lines
if thr != 0:
sns.scatterplot(data=dfc, x='c1', y='c2', color='gray')
sns.scatterplot(data=dfs, x='c1', y='c2', color='red')
else:
sns.scatterplot(data=dfc, x='c1', y='c2', hue='brain_region')
plt.grid(False)
sns.despine()
return df2
def laser_triggered_pcs(PC, pre, post, M, mouse, kcuts=[], min_laser=20, pzscore_pc=False, local_pzscore=True,
pplot=True, ci=None, refractory_rule=False, ma_thr=10, ma_rem_exception=False, rnd_laser=False, seed=1,
config_file='mouse_config.txt', laser_dur=-1, start_mode='floor'):
"""
Calculated the time course of the provided PCs relative to the laser onset.
Parameters
----------
PC : np.array
Each row corresponds to one PC.
pre : float
Time before laser onset.
post : float
Time after laser onset.
M : np.array
Hypnogram.
mouse : str
Mouse name.
kcuts : list of tuples, optional
DESCRIPTION. The default is [].
min_laser : float, optional
Minimum duration of laser train. If laser duration < $min_laser,
disregard the laser trial.
pzscore_pc : bool, optional
If true, z-score the PCs (across entire recording)
local_pzscore : bool, optional
DESCRIPTION. The default is True.
pplot : bool, optional
If True, plot figures summarizing results.
ci : float or None, optional
Confidence interval for plots.
refractory_rule : TYPE, optional
DESCRIPTION. The default is False.
config_file : str, optional
Mouse recordings configuration file.
laser_dur : float
Duration of laser; if -1, calculate laser duration using laser_*.mat file
start_mode : str
'floor': The first bins that touches the laser, is the laser onset;
'before': Time point 0 is the bin that just does not touch the laser.
Returns
-------
df : pandas.DataFrame
with columns ['mouse', 'time', 'val', 'valz',
'pc', 'lsr_start', 'start_state', 'start_state_int',
'state', 'success', 'refractory']
'state': brain state sequence
'start_state_int': How long does it take till the brain state at laser onset
switches to a different state
'valz': PC values during laser trial, for each trial the PC vector
is z-scored.
'rem_delay': Delay from laser to REM onset; if there's no REM the value is set to -1
"""
if rnd_laser:
np.random.seed(seed)
dt = 2.5
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
nhypno = M.shape[0]
ndim = PC.shape[0]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
# flatten out MAs #########################################################
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>0) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
###########################################################################
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*dt)
dt = nbin * (1.0/sr)
ipre = int(pre/dt)
ipost = int(post/dt)
t = np.arange(-ipre, ipost)*dt
nt = len(t)
#######################################################################
# get laser start and end index after excluding kcuts: ################
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
idxs, idxe = sleepy.laser_start_end(lsr)
if start_mode == 'floor':
idxs = [int(i/nbin) for i in idxs]
elif start_mode == 'before':
[int(i/nbin)-1 for i in idxs]
else:
idxs = [int(i/nbin)+1 for i in idxs]
if laser_dur == -1:
idxe = [int(i/nbin) for i in idxe]
b = [i for i in idxe if i in tidx]
a = [i for i in idxs if i in tidx]
b = np.array(b)
a = np.array(a)
laser_dur = (np.mean(b-a) + 1) * dt
print('Laser duration: %f' % laser_dur)
else:
idxe = [(int(i + int(laser_dur/dt))) for i in idxs]
# randomize laser #####################################################
dur = int(laser_dur/dt)
if rnd_laser:
idxs_rnd = []
idxe_rnd = []
tmp = np.random.randint(dur, idxs[0]-dur)
idxs_rnd.append(tmp)
idxe_rnd.append(tmp+dur)
for (a,b) in zip(idxe[0:-1], idxs[1:]):
if a+2*dur < b-dur:
tmp = np.random.randint(a+2*dur,b-dur)
idxs_rnd.append(tmp)
idxe_rnd.append(tmp+dur)
idxs = idxs_rnd
idxe = idxe_rnd
#######################################################################
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
laser_idx = np.array(laser_idx)
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser_idx = laser_idx[laser_idx < nlsr]
laser[laser_idx] = 1
laser = laser[tidx]
idxs = [s[0] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
idxe = [s[-1] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
#######################################################################
label = []
for p in range(1, ndim+1):
l = 'pc' + str(p)
label.extend([l]*nt)
# zscore pcs:
if pzscore_pc:
for i in range(PC.shape[0]):
PC[i,:] = (PC[i,:]- PC[i,:].mean()) / PC[i,:].std()
data = []
ev = 0
for i,j in zip(idxs, idxe):
if i >= ipre and i+ipost < nhypno:
if (j-i)*dt < min_laser:
continue
idx = np.arange(i-ipre, i+ipost).astype('int')
pc_cut = PC[:,idx]
pc_cut_z = pc_cut.copy()
for k in range(ndim):
pc_cut_z[k,:] = (pc_cut[k,:] - pc_cut[k,:].mean()) / pc_cut[k,:].std()
m_lsr = M[i:j+1]
# repeat m_cut ndim-times
m_cut = np.tile(M[idx], (ndim,))
vec = np.reshape(pc_cut, (ndim*nt,))
vecz = np.reshape(pc_cut_z, (ndim*nt,))
tvec = np.tile(t, (ndim,))
lsr_rem = 'no'
if 1 in m_lsr:
lsr_rem='yes'
start_state = M[i]
rem_delay = -1
if lsr_rem == 'yes':
rem_delay = np.where(m_lsr == 1)[0][0] * dt
# duration of interval of brainstate start_state till the
# brain_state switches:
start_state_int = len(m_lsr)*dt
a = np.where(m_lsr != start_state)[0]
if len(a) > 0:
start_state_int = len(a) * dt
l = i-1
while M[l] != 1 and l>0:
l = l-1
refr = 'no'
if M[l] == 1 and l != i-1:
v = l
while M[v] == 1 and v > 0:
v = v-1
v = v+1
dur_rem_pre = (l-v+1)*dt
inrem = len(np.where(M[l+1:i] == 3)[0]) * dt
if inrem <= dur_rem_pre * 2:
refr = 'yes'
data += zip([mouse]*nt*ndim, [ev]*nt*ndim, tvec, vec, vecz,
label, [i]*nt*ndim, [start_state]*nt*ndim, [start_state_int]*nt*ndim,
m_cut, [lsr_rem]*nt*ndim, [rem_delay]*nt*ndim, [refr]*nt*ndim)
ev += 1
df = pd.DataFrame(data=data, columns=['mouse', 'ev', 'time', 'val', 'valz',
'pc', 'lsr_start', 'start_state', 'start_state_int',
'state', 'success', 'rem_delay', 'refractory'])
if pplot:
plt.figure()
if local_pzscore:
sns.lineplot(data=df, x='time', y='valz', hue='pc', palette='husl')
else:
sns.lineplot(data=df, x='time', y='val', hue='pc', palette='husl')
plt.xlim([t[0], t[-1]])
sns.despine()
pcs = df.pc.unique()
data = {p:[] for p in pcs}
data_start = []
for p in pcs:
for si in df.lsr_start.unique():
a = np.array(df[(df.lsr_start==si) & (df.pc == p)]['val'])
m_cut = np.array(df[(df.lsr_start==si) & (df.pc == p)]['state'])
t = np.array(df[(df.lsr_start==si) & (df.pc == p)]['time'])
l = np.array(df[(df.lsr_start==si) & (df.pc == p)]['success'])[0]
if p == 'pc1':
if l == 'yes':
rem_start = np.where((t >= 0) & (m_cut == 1) )[0][0]
data_start.append(t[rem_start])
else:
data_start.append(-1)
data[p].append(a)
for p in pcs:
data[p] = np.array(data[p])
f, axes = plt.subplots(nrows=ndim, ncols=1, sharex='all')
for ax,p in zip(axes, pcs):
mx = data[p].copy()
if local_pzscore:
for i in range(mx.shape[0]):
mx[i,:] = (mx[i,:] - mx[i,:].mean()) / mx[i,:].std()
ax.pcolorfast(t, range(0, mx.shape[0]+1), mx, cmap='jet')
for ii in range(mx.shape[0]):
tt = data_start[ii]
if tt >= 0:
ax.plot([tt,tt], [ii,ii+1], color='black', lw=3)
ax.set_ylabel('')
return df
def laser_triggered_pcs_fine(PC, pre, post, M, mouse, kcuts=[], min_laser=20, pzscore_pc=False, local_pzscore=True, ndown=250,
pplot=True, ci=None, refractory_rule=False, ma_thr=10, ma_rem_exception=False, rnd_laser=False, seed=1,
config_file='mouse_config.txt', laser_dur=-1):
"""
Calculated the time course of the provided PCs relative to the laser onset using finer timescale.
Parameters
----------
PC : np.array
Each row corresponds to one PC.
pre : float
Time before laser onset.
post : float
Time after laser onset.
M : np.array
Hypnogram.
mouse : str
Mouse name.
kcuts : list of tuples, optional
DESCRIPTION. The default is [].
min_laser : float, optional
Minimum duration of laser train. If laser duration < $min_laser,
disregard the laser trial.
pzscore_pc : bool, optional
If true, z-score the PCs (across entire recording)
local_pzscore : bool, optional
DESCRIPTION. The default is True.
pplot : bool, optional
If True, plot figures summarizing results.
ci : float or None, optional
Confidence interval for plots.
refractory_rule : TYPE, optional
DESCRIPTION. The default is False.
config_file : str, optional
Mouse recordings configuration file.
laser_dur : float
Duration of laser; if -1, calculate laser duration using laser_*.mat file
Returns
-------
df : pandas.DataFrame
with columns ['mouse', 'time', 'val', 'valz',
'pc', 'lsr_start', 'start_state', 'start_state_int',
'state', 'success', 'refractory']
'state': brain state sequence
'start_state_int': How long does it take till the brain state at laser onset
switches to a different state
'valz': PC values during laser trial, for each trial the PC vector
is z-scored.
'rem_delay': Delay from laser to REM onset; if there's no REM the value is set to -1
"""
dt = 2.5
NDOWN = ndown
NUP = int(dt / (0.001 * NDOWN))
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
M = sleepy.load_stateidx(ppath, name)[0]
# flatten out MAs #########################################################
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
###########################################################################
dt = dt / NUP
M = upsample_mx(M, NUP)
nhypno = len(M)
M = M[0:nhypno]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
sr = sleepy.get_snr(ppath, name)
nbin = NDOWN
ipre = int(pre/dt)
ipost = int(post/dt)
t = np.arange(-ipre, ipost)*dt
nt = len(t)
#######################################################################
# get laser start and end index after excluding kcuts: ################
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
if laser_dur == -1:
idxe = [int(i/nbin) for i in idxe]
b = [i for i in idxe if i in tidx]
a = [i for i in idxs if i in tidx]
b = np.array(b)
a = np.array(a)
laser_dur = (np.mean(b-a) + 1) * dt
print('Laser duration: %f' % laser_dur)
else:
idxe = [(int(i + int(laser_dur/dt))) for i in idxs]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
laser_idx = np.array(laser_idx)
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser_idx = laser_idx[laser_idx < nlsr]
laser[laser_idx] = 1
laser = laser[tidx]
idxs = [s[0] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
idxe = [s[-1] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
#######################################################################
ndim = np.min(PC.shape)
label = []
for p in range(1, ndim+1):
l = 'pc' + str(p)
label.extend([l]*nt)
# zscore pcs:
if pzscore_pc:
for i in range(PC.shape[0]):
PC[i,:] = (PC[i,:]- PC[i,:].mean()) / PC[i,:].std()
data = []
ev = 0
for i,j in zip(idxs, idxe):
if i > ipre and i+ipost < nhypno:
if (j-i)*dt < min_laser:
continue
idx = np.arange(i-ipre, i+ipost).astype('int')
pc_cut = PC[:,idx]
pc_cut_z = pc_cut.copy()
for k in range(ndim):
pc_cut_z[k,:] = (pc_cut[k,:] - pc_cut[k,:].mean()) / pc_cut[k,:].std()
m_lsr = M[i:j+1]
# repeat m_cut ndim-times
m_cut = np.tile(M[idx], (ndim,))
vec = np.reshape(pc_cut, (ndim*nt,))
vecz = np.reshape(pc_cut_z, (ndim*nt,))
tvec = np.tile(t, (ndim,))
lsr_rem = 'no'
if 1 in m_lsr:
lsr_rem='yes'
start_state = M[i]
rem_delay = -1
if lsr_rem == 'yes':
rem_delay = np.where(m_lsr == 1)[0][0] * dt
# duration of interval of brainstate start_state till the
# brain_state switches:
start_state_int = len(m_lsr)*dt
a = np.where(m_lsr != start_state)[0]
if len(a) > 0:
start_state_int = len(a) * dt
l = i-1
while M[l] != 1 and l>0:
l = l-1
refr = 'no'
if M[l] == 1 and l != i-1:
v = l
while M[v] == 1 and v > 0:
v = v-1
v = v+1
dur_rem_pre = (l-v+1)*dt
inrem = len(np.where(M[l+1:i] == 3)[0]) * dt
if inrem <= dur_rem_pre * 2:
refr = 'yes'
data += zip([mouse]*nt*ndim, [ev]*nt*ndim, tvec, vec, vecz,
label, [i]*nt*ndim, [start_state]*nt*ndim, [start_state_int]*nt*ndim,
m_cut, [lsr_rem]*nt*ndim, [rem_delay]*nt*ndim, [refr]*nt*ndim)
ev += 1
df = pd.DataFrame(data=data, columns=['mouse', 'ev', 'time', 'val', 'valz',
'pc', 'lsr_start', 'start_state', 'start_state_int',
'state', 'success', 'rem_delay', 'refractory'])
if pplot:
plt.figure()
if local_pzscore:
sns.lineplot(data=df, x='time', y='valz', hue='pc', palette='husl')
else:
sns.lineplot(data=df, x='time', y='val', hue='pc', palette='husl')
plt.xlim([t[0], t[-1]])
sns.despine()
pcs = df.pc.unique()
data = {p:[] for p in pcs}
data_start = []
for p in pcs:
for si in df.lsr_start.unique():
a = np.array(df[(df.lsr_start==si) & (df.pc == p)]['val'])
m_cut = np.array(df[(df.lsr_start==si) & (df.pc == p)]['state'])
t = np.array(df[(df.lsr_start==si) & (df.pc == p)]['time'])
l = np.array(df[(df.lsr_start==si) & (df.pc == p)]['success'])[0]
if p == 'pc1':
if l == 'yes':
rem_start = np.where((t >= 0) & (m_cut == 1) )[0][0]
data_start.append(t[rem_start])
else:
data_start.append(-1)
data[p].append(a)
for p in pcs:
data[p] = np.array(data[p])
f, axes = plt.subplots(nrows=ndim, ncols=1, sharex='all')
for ax,p in zip(axes, pcs):
mx = data[p].copy()
if local_pzscore:
for i in range(mx.shape[0]):
mx[i,:] = (mx[i,:] - mx[i,:].mean()) / mx[i,:].std()
ax.pcolorfast(t, range(0, mx.shape[0]+1), mx, cmap='jet')
for ii in range(mx.shape[0]):
tt = data_start[ii]
if tt >= 0:
ax.plot([tt,tt], [ii,ii+1], color='black', lw=3)
ax.set_ylabel('')
return df
# --------------------------------------------------------------
# replacement for laser_triggered_pcs()
# randomises *where* a laser trial is assumed to occur
# but only keeps trials whose PC2 value (row index 1) exceeds pc_thr
# --------------------------------------------------------------
def random_pc2_triggered_pcs(
PC, # (n_dim × T) principal components
pre, post, # flank lengths in seconds
M, # hypnogram (1=REM,2=Wake,3=NREM ...)
mouse,
n_trials = 40, # how many pseudo-laser trials
pc_thr = 2.0, # threshold for PC2
dt = 2.5, # bin size of M and PC (s)
kcuts = None, # list of (t0,t1) tuples to exclude
pzscore_pc = False,
local_pzscore= True,
pplot = True,
config_file = 'mouse_config.txt'):
"""
Pick `n_trials` random time points where PC2 > pc_thr and
build the same trial-aligned PC dataframe as laser_triggered_pcs().
"""
# <- your helpers
# ---------- basic sizes ----------
ndim, T = PC.shape
ipre = int(round(pre / dt))
ipost = int(round(post / dt))
nt = ipre + ipost
tvec = np.arange(-ipre, ipost) * dt
# ---------- mask out k-cuts (if any) ------------------------
kmask = np.ones(T, dtype=bool)
if kcuts:
for t0, t1 in kcuts:
kmask[int(t0 / dt) : int(t1 / dt)] = False
# ---------- admissible centres --------------------------------
ok = np.ones(T, dtype=bool)
ok[:ipre] = False # need full window on the left
ok[-ipost:] = False # ... and on the right
ok &= kmask
# PC2 threshold
ok &= PC[1, :] > pc_thr
centres = np.where(ok)[0]
if len(centres) < n_trials:
raise RuntimeError(f"Only {len(centres)} time points fulfil PC2>{pc_thr}")
rng = np.random.default_rng(1)
sel = rng.choice(centres, size=n_trials, replace=False)
# ---------- assemble dataframe (same columns as before) -------
rows = []
for ev, c in enumerate(sel):
idx = np.arange(c-ipre, c+ipost)
start_state = M[c]
m_window = M[idx]
# per-trial z-score
win = PC[:, idx].copy()
win_z = (win - win.mean(1, keepdims=True)) / win.std(1, keepdims=True)
for k, lab in enumerate([f'pc{i+1}' for i in range(ndim)]):
rows.append(pd.DataFrame({
'mouse' : mouse,
'ev' : ev,
'time' : tvec,
'val' : win[k],
'valz' : win_z[k] if local_pzscore else win[k],
'pc' : lab,
'lsr_start' : c,
'start_state' : start_state,
'start_state_int': (np.where(m_window != start_state)[0][0] * dt
if np.any(m_window != start_state) else len(m_window)*dt),
'state' : m_window,
'success' : 'rnd', # no real “success”
'rem_delay' : -1,
'refractory' : 'na'
}))
df = pd.concat(rows, ignore_index=True)
# ---------- optional global z-score ---------------------------
if pzscore_pc:
df['val'] = (df['val'] - df['val'].mean()) / df['val'].std()
df['valz'] = (df['valz'] - df['valz'].mean()) / df['valz'].std()
# ---------- optional quick QC plot ---------------------------
if pplot:
sns.lineplot(data=df, x='time', y='valz' if local_pzscore else 'val',
hue='pc', palette='husl', errorbar='se')
plt.axvline(0, ls='--', c='k', lw=.8)
plt.xlabel('Time (s)'); plt.ylabel('PC value (z-score)')
sns.despine(); plt.tight_layout()
return df
def laser_triggered_frs(units, pre, post, mouse, kcuts=[], ma_thr=10, ma_rem_exception=True,
pzscore=True, nsmooth=0, detrend=True, min_laser=20,
config_file='mouse_config.txt'):
"""
Parameters
----------
units : TYPE
DESCRIPTION.
pre : TYPE
DESCRIPTION.
post : TYPE
DESCRIPTION.
mouse : TYPE
DESCRIPTION.
kcuts : TYPE, optional
DESCRIPTION. The default is [].
ma_thr : TYPE, optional
DESCRIPTION. The default is 10.
ma_rem_exception : TYPE, optional
DESCRIPTION. The default is True.
pzscore : TYPE, optional
DESCRIPTION. The default is True.
nsmooth : TYPE, optional
DESCRIPTION. The default is 0.
detrend : TYPE, optional
DESCRIPTION. The default is True.
min_laser : TYPE, optional
DESCRIPTION. The default is 20.
config_file : TYPE, optional
DESCRIPTION. The default is 'mouse_config.txt'.
Returns
-------
df : TYPE
DESCRIPTION.
"""
dt = 2.5
ddir = load_config(config_file)[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
M = sleepy.load_stateidx(ppath, name)[0]
nhypno = M.shape[0]-1
ndim = units.shape[1]
tidx = np.arange(0, nhypno)
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
# flatten out MAs #########################################################
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>0) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
###########################################################################
if os.path.isfile(os.path.join(ddir, 'laser_%s.mat' % name)):
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*dt)
dt = nbin * (1.0/sr)
ipre = int(pre/dt)
ipost = int(post/dt)
t = np.arange(-ipre, ipost)*dt
nt = len(t)
#######################################################################
# get laser start and end index after excluding kcuts: ################
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
idxe = [int(i/nbin) for i in idxe]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser[laser_idx] = 1
laser = laser[tidx]
idxs = [s[0] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
idxe = [s[-1] for s in sleepy.get_sequences(np.where(laser == 1)[0])]
#######################################################################
unitIDs = [unit for unit in units.columns if '_' in unit]
nvar = len(unitIDs)
R = np.zeros((nvar, nhypno)) # dimensions: number of units x time bins
for i,unit in enumerate(unitIDs):
tmp = sleepy.smooth_data(np.array(units[unit]),nsmooth)
tmp = tmp[tidx]
if detrend:
tmp = scipy.signal.detrend(tmp)
if pzscore:
R[i,:] = (tmp - tmp.mean()) / tmp.std()
else:
R[i,:] = tmp
label = []
for p in range(0, nvar):
l = unitIDs[p]
label.extend([l]*nt)
ms_id = [mouse + '_' + i for i in label]
data = []
for i,j in zip(idxs, idxe):
if i > ipre and i+ipost < nhypno:
if (j-i)*dt < min_laser:
continue
idx = np.arange(i-ipre, i+ipost).astype('int')
r_cut = R[:,idx]
m_lsr = M[i:j+1]
# repeat m_cut ndim-times
m_cut = np.tile(M[idx], (nvar,))
vec = np.reshape(r_cut, (nvar*nt,))
tvec = np.tile(t, (ndim,))
lsr_rem = 'no'
if 1 in m_lsr:
lsr_rem='yes'
start_state = M[i]
# duration of interval of brainstate start_state till the
# brain state switches:
start_state_int = len(m_lsr)*dt
a = np.where(m_lsr != start_state)[0]
if len(a) > 0:
start_state_int = len(a) * dt
data += zip([mouse]*nt*nvar, tvec, vec,
label, ms_id, [i]*nt*nvar, [start_state]*nt*nvar, [start_state_int]*nt*nvar,
m_cut, [lsr_rem]*nt*nvar)
df = pd.DataFrame(data=data, columns=['mouse', 'time', 'fr', 'ID',
'ms_id', 'lsr_start', 'start_state', 'start_state_int',
'state', 'success'])
return df
def plot_trajectories(PC, M, pre, post, istate=1, dt=2.5, ma_thr=10, ma_rem_exception=False,
kcuts=[],min_dur=0, pre_state=0, state_num=[], ax='', lw=1, coords=[0,1], mouse=''):
tidx = np.arange(0, len(M))
# NEW 07/01/22:
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
nhypno = np.min((len(M), PC.shape[1]))
M = M[0:nhypno]
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
if ax == '':
plt.figure()
ax = plt.axes([0.15, 0.15, 0.7, 0.7])
ax_3d = False
if ax == '3D':
ax_3d = True
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
clrs = [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8], [1, 0.2, 0.2]]
ipre = int(pre/dt)
ipost = int(post/dt)
if istate in [1,2,3]:
if istate > 0:
seq = sleepy.get_sequences(np.where(M==istate)[0])
else:
seq = [np.arange(0, len(M))]
if len(state_num) > 0:
seq = [seq[i] for i in state_num]
for s in seq:
si = s[0]
sj = s[-1]
if si-ipre > 0 and sj+ipost < nhypno and len(s)*dt > min_dur and M[si-1]==pre_state:
#print(len(s)*dt)
mcut = M[si-ipre:sj+ipost]
k = si-ipre
p = mcut[0]
kold = k
while k < sj+ipost:
while M[k] == p and k < sj+ipost:
k+=1
if ax_3d:
ax.plot(PC[coords[0],kold:k+1], PC[coords[1], kold:k+1], PC[coords[2], kold:k+1], color=clrs[int(p)], lw=lw)
else:
ax.plot(PC[coords[0],kold:k+1], PC[coords[1], kold:k+1], color=clrs[int(p)], lw=lw)
p = M[k]
kold = k
plt.xlabel('PC1')
plt.ylabel('PC2')
if istate == 'laser':
ddir = load_config('mouse_config.txt')[mouse]['SL_PATH']
ppath, name = os.path.split(ddir)
sr = sleepy.get_snr(ppath, name)
nbin = int(np.round(sr)*2.5)
dt = nbin * (1.0/sr)
#######################################################################
# get laser start and end index after excluding kcuts: ################
lsr = so.loadmat(os.path.join(ddir, 'laser_%s.mat' % name), squeeze_me=True)['laser']
idxs, idxe = sleepy.laser_start_end(lsr)
idxs = [int(i/nbin) for i in idxs]
idxe = [int(i/nbin) for i in idxe]
laser_idx = []
for (si,sj) in zip(idxs, idxe):
laser_idx += list(range(si,sj+1))
nlsr = int(np.floor(lsr.shape[0]/nbin))
laser = np.zeros((nlsr,))
laser[laser_idx] = 1
laser = laser[tidx]
idxs = np.array([s[0] for s in sleepy.get_sequences(np.where(laser == 1)[0])])
idxe = np.array([s[-1] for s in sleepy.get_sequences(np.where(laser == 1)[0])])
#######################################################################
for (si,sj) in zip(idxs[state_num], idxe[state_num]):
if (sj-si) * dt < 20:
continue
mcut = M[si:sj+1]
if 1 in mcut:
ax.plot(PC[coords[0],si-ipre:si], PC[coords[1], si-ipre:si], color='cornflowerblue', lw=lw)
ax.plot(PC[coords[0],si-1:sj], PC[coords[1], si-1:sj], color='blue', lw=lw)
ax.plot(PC[coords[0],si-1], PC[coords[1], si-1], color='orange', marker='o', lw=lw)
else:
ax.plot(PC[coords[0],si-ipre:si], PC[coords[1], si-ipre:si], color='pink', lw=lw)
ax.plot(PC[coords[0],si-1:sj], PC[coords[1], si-1:sj], color='red', lw=lw)
ax.plot(PC[coords[0],si-1], PC[coords[1], si-1], color='green', marker='o', lw=lw)
def pc_state_space(PC, M, ma_thr=10, ma_rem_exception=False, kcuts=[], dt=2.5, ax='', nrem2wake=False, nrem2wake_step=4,
pscatter=True, local_coord=False, outline_std=True, rem_onset=False, rem_offset=False, rem_offset_only_nrem=False, show_avgtraj=False,
pre_win=30, post_win=0, rem_min_dur=0, break_out=True, break_in=False, prefr=False, scale=1.645, pzscore_pc=False):
"""
Plot for each time point the population activity within the 2D state space spanned
by PC[0,:] and PC[1,:]
Parameters
----------
PC : np.array
each row are the PC coefficients or scores
M : np.array
brain state annotation.
1 - REM, 2 - Wake, 3 - NREM
ma_thr : float, optional
Microarousal threshold. The default is 10.
ma_rem_exception : bool, optional
If True, don't set a wake episodes after REM that is shorter then $ma_thr to NREM.
The default is False.
kcuts : list of tuples, optional
Each tuple describes the start and end of an interval to be discarded. The default is [].
dt : float, optional
Time bin of firing rates and hypnogram. The default is 2.5.
ax : figure axis handle, optional
If $as is provided, use it to draw all plots on it.
Otherwise generate new figure.
The default is ''.
nrem2wake : Bool, optional
if True, show NREM->Wake transitions
nrem2wake_step: int, optional
Show every nrem2wake_step-th NREM->Wake transition; otherwise
it gets too clusttered
pscatter : bool, optional
If True, draw each (dimensionally reduced) population vector in a scatter plot.
The default is True.
local_coord : bool, optional
If True, draw into the NREM ellipse a local coordinate system. The default is False.
outline_std : bool
If True, draw outline of one std using an ellipse;
if False, use sns.kdeplot to draw outline of data spread
rem_onset : bool
If True, color-code the REM onset and the preceding $pre_win seconds
rem_offset : bool
If True, plot REM->Wake->NREM transitions
rem_offset_only_nrem: bool
If True, plot only REM->Wake->NREM transitions where NREM occurs within
the $post_win interval following the REM offset.
show_avgtraj : bool
If True, show average across all trajectories instead of each individual
trajectory
rem_mindur: float
Only REM episodes >= rem_mindur are considered
break_out: bool
if True, draw a dot where a NREM to REM trajectory
leaves the NREM subspace, defined by $scale
prefr: bool,
If True, also draw outline of refractory period within state space
scale: float
$scale == 1 means that the drawn ellipse outlines one standard deviation
for each subspace.
$scale == 1.645 outlines the area (of the fitted Gaussian)
that comprises 90% of the data distribution.
$scale == 1.96 outlines 95% of the distribution
pzscore_pc: bool
If True, zscore PCs.
Returns
-------
ax : plt.axes
return axes of current figure.
df_breakout : pd.DataFrame
with columns ['angle':first_angles, 'pc1':c1, 'pc2':c2, 'pc1_org':porg1, 'pc2_org':porg2]
Describes along each NREM->REM trajectory the first point and angle leaving the
NREM subspace
df_breakin : pd.DataFrame
Describes for each REM->Wake->NREM transition trajectory the first point within NREM.
df_ampl: pd.DataFrame
Describes for each NREM->REM ($nrem2wake=False) or NREM->Wake ($nrem2wake=True) the maximum
PC1 and PC2 values of the preceding trajectory of during $pre_win.
"""
bs_map = {'rem':[0, 1, 1], 'wake':[0.6, 0, 1], 'nrem':[0.8, 0.8, 0.8]}
tidx = np.arange(0, len(M))
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
nhypno = np.min((len(M), PC.shape[1]))
M = M[0:nhypno]
# zscore PCs:
if pzscore_pc:
for i in range(PC.shape[0]):
PC[i,:] = (PC[i,:] - PC[i,:].mean()) / PC[i,:].std()
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
if ax == '':
plt.figure()
ax = plt.axes([0.2, 0.15, 0.7, 0.7])
clrs = [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.6, 0.6, 0.6], [1, 0.2, 0.2]]
state_idx = {}
for s in [1,2,3]:
idx = np.where(M==s)[0]
state_idx[s] = idx
# get all indices for REM, Wake, and NREM
for s in [1,2,3]:
idx = state_idx[s]
C = PC[0:2,idx].T
if outline_std:
mean = C.mean(axis=0)
covar = np.cov(C.T)
v, w = linalg.eigh(covar)
# columns of w are the eigenvectors
# v are the eigenvalues in ascending order
# 2 * scale * std (the 2 is to transform the radius to diameter)
v = 2.0 * scale * np.sqrt(v)
u = w[0] / linalg.norm(w[0])
# Plot an ellipse to show the Gaussian component
angle = np.arctan(u[1] / u[0])
angle = 180.0 * angle / np.pi # convert to degrees
ell = mpl.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=clrs[s], lw=2)
ell.set_clip_box(ax.bbox)
ell.set_alpha(0.3)
ax.add_artist(ell)
else:
sns.kdeplot(x=C[:, 0], y=C[:, 1], ax=ax, color=clrs[s], fill=True, alpha=0.8, levels=[0.25, 0.5, 0.75, 1])
# Add refractory period to state space
if prefr:
refr_color = 'maroon'
df_refr, refr_vec, _ = add_refr(M)
refr_idx = np.where(refr_vec == 1)[0]
nr_idx = np.where(M==3)[0]
idx = np.intersect1d(refr_idx, nr_idx)
C = PC[0:2,idx].T
if outline_std:
mean = C.mean(axis=0)
covar = np.cov(C.T)
v, w = linalg.eigh(covar)
# columns of w are the eigenvectors
# v are the eigenvalues in ascending order
# 2 * scale * std (the 2 is to transform the radius to diameter)
v = 2.0 * scale * np.sqrt(v)
u = w[0] / linalg.norm(w[0])
# Plot an ellipse to show the Gaussian component
angle = np.arctan(u[1] / u[0])
angle = 180.0 * angle / np.pi # convert to degrees
ell = mpl.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=refr_color, lw=2)
ell.set_clip_box(ax.bbox)
ell.set_alpha(0.3)
ax.add_artist(ell)
else:
sns.kdeplot(x=C[:, 0], y=C[:,1], ax=ax, color=refr_color, fill=True, alpha=0.8, levels=[0.25, 0.5, 0.75, 1])
if local_coord:
idx = state_idx[3]
C = PC[0:2,idx].T
pca = PCA(n_components=2)
pca.fit(C)
mean_x, mean_y = np.mean(C, axis=0)
pc1, pc2 = pca.components_
if pc1[1] < 0:
pc1 = -pc1
# Calculate the standard deviations along the principal components
std_x, std_y = np.sqrt(pca.explained_variance_)
# Scale the principal components by the standard deviations
pc1_scaled = pc1 * std_x * 1 * scale
pc2_scaled = pc2 * std_y * 1 * scale
origin = [mean_x], [mean_y]
ax.quiver(*origin, pc1_scaled[0], pc1_scaled[1], angles='xy', scale_units='xy', scale=1, color='black', label='PC1 (Scaled)')
ax.quiver(*origin, pc2_scaled[0], pc2_scaled[1], angles='xy', scale_units='xy', scale=1, color='black', label='PC2 (Scaled)')
pc1_min = PC[0,:].min()
pc1_max = PC[0,:].max()
pc2_min = PC[1,:].min()
pc2_max = PC[1,:].max()
d1 = pc1_max - pc1_min
d2 = pc2_max - pc2_min
ax.set_xlim(pc1_min - 0.1*d1, pc1_max + 0.1*d1)
ax.set_ylim(pc2_min - 0.1*d2, pc2_max + 0.1*d2)
if pscatter:
for s in [1,2,3]:
idx = state_idx[s]
ax.scatter(PC[0,idx], PC[1,idx], color=clrs[s], s=1, alpha=1)
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
sns.despine()
# just the onset of REM as single dot:
ipre_win = int(pre_win/dt)
ipost_win = int(post_win/dt)
rem_start = [s[0] for s in sleepy.get_sequences(np.where(M==1)[0]) if len(s)*dt >= rem_min_dur and s[0]*dt >= pre_win and s[0]+ipost_win < len(M)]
if nrem2wake:
rem_start = [s[0] for s in sleepy.get_sequences(np.where(M==2)[0]) if len(s)*dt >= rem_min_dur and s[0]*dt >= pre_win and M[s[0]-1]==3]
rem_start = rem_start[1::nrem2wake_step]
# show trajectories for REM onset
if rem_onset:
if not show_avgtraj:
for r in rem_start[:]:
if not nrem2wake:
plt.plot(PC[0,r], PC[1,r], '*', color=bs_map['rem'], markersize=10, zorder=3)
else:
plt.plot(PC[0,r], PC[1,r], '*', color=bs_map['wake'], markersize=10, zorder=3)
for r in rem_start:
# NEW - 04/17/24:
pc1, pc2 = PC[0,r-ipre_win:r+ipost_win+1], PC[1,r-ipre_win:r+ipost_win+1]
sm = _jet_plot(pc1, pc2, ax, lw=2, cmap='magma')
else:
tmp1, tmp2 = [], []
for r in rem_start:
pc1, pc2 = PC[0,r-ipre_win:r+ipost_win+1], PC[1,r-ipre_win:r+ipost_win+1]
tmp1.append(pc1)
tmp2.append(pc2)
pc1 = np.array(tmp1).mean(axis=0)
pc2 = np.array(tmp2).mean(axis=0)
sm = _jet_plot(pc1, pc2, ax, lw=2, cmap='magma')
cbar = plt.colorbar(sm, ax=ax, orientation='vertical', pad=0.05, shrink=0.6) # pad adjusts the distance between the plot and colorbar
sm.set_clim(-pre_win, post_win)
cbar.set_ticks([-pre_win, post_win])
cbar.set_label("Time (s)")
if rem_offset:
ipre_win = int(pre_win/dt)
ipost_win = int(post_win/dt)
rem_end_all = [s[-1]+1 for s in sleepy.get_sequences(np.where(M==1)[0]) if len(s)*dt >= rem_min_dur and s[-1]+ipost_win < len(M) and s[-1]>ipre_win]
tmp = []
# for each REM offset search for the end of the following wake episode
for r in rem_end_all:
i = r
while i < len(M) and M[i] != 3:
i += 1
dur = (i - r) * dt
if dur <= post_win:
tmp.append(r)
rem_end_nrem = tmp
if rem_offset_only_nrem:
rem_end = rem_end_nrem
else:
rem_end = rem_end_all
for r in rem_end:
pc1, pc2 = PC[0,r-ipre_win:r+ipost_win+1], PC[1,r-ipre_win:r+ipost_win+1]
sm = _jet_plot(pc1, pc2, ax, lw=2, cmap='magma')
plt.plot(PC[0,r+1], PC[1,r+1], '*', color=bs_map['wake'], markersize=10, zorder=3)
cbar = plt.colorbar(sm, ax=ax, orientation='vertical', pad=0.05, shrink=0.6) # pad adjusts the distance between the plot and colorbar
sm.set_clim(-pre_win, post_win)
cbar.set_ticks([-pre_win, post_win])
cbar.set_label("Time (s)")
df_breakout = []
data_ampl = []
df_ampl = []
if break_out:
seq = sleepy.get_sequences(np.where(M==1)[0])
idx = state_idx[3]
C = PC[0:2,idx].T
meanc = C.mean(axis=0)
covar = np.cov(C.T)
pca = PCA(n_components=2)
pca.fit(C)
# get the eigenvectors of the covariance matrix:
pc1, pc2 = pca.components_
if pc1[1] < 0:
pc1 = -pc1
w = np.zeros((2,2))
w[:,0] = pc1
w[:,1] = pc2
first = []
if not nrem2wake:
# NREM -> REM
for r in rem_start:
ifirst = last_subspace_point(r, PC[0:2,:], meanc, covar, scale=scale)
first.append(ifirst)
plt.plot(PC[0,ifirst], PC[1,ifirst], 'ro')
else:
# NREM -> Wake
for r in rem_start:
if not is_in_ellipse(r, PC[0:2,:], meanc, covar, scale=scale):
ifirst = last_subspace_point(r, PC[0:2,:], meanc, covar, scale=scale)
first.append(ifirst)
#plt.plot(PC[0,ifirst], PC[1,ifirst], 'ro')
for r in rem_start:
tmp1 = PC[0,r-ipre_win:r+1]
tmp2 = PC[1,r-ipre_win:r+1]
# calculation current REM duration
state = M[r]
ii = r
while ii < len(M) and M[ii] == state:
ii = ii+1
dur = (ii-r) * dt
# END #################################
#data_ampl += [[tmp1.max(), tmp2.max(), dur]]
data_ampl += [[tmp1.max(), tmp2.max(), dur]]
df_ampl = pd.DataFrame(data=data_ampl, columns=['pc1', 'pc2', 'dur_post'])
# Alternative way of calculating the eigenvectors of the
# covariance matrix:
# v, w = linalg.eigh(covar)
# ii = np.argsort(v)[::-1]
# v = v[ii]
# w = w[:,ii]
first_angles = []
c1, c2 = [], []
porg1, porg2 = [], []
porg1_rel, porg2_rel = [], []
for ifirst in first:
# take the first point outside the NREM subspace
p = PC[0:2,ifirst]
# center the point
pctr = p - meanc
# project it onto the eigenvectors
a = np.dot(pctr, w)
x, y = a[0], a[1]
# calculate the angle in radians
theta = math.atan2(x, y)
# Convert the angle to degrees;
# NOTE: 0 deg. corresponds to 3; 90 deg. corresponds to 12
angle_degrees = math.degrees(theta)
first_angles.append(angle_degrees)
c1.append(a[0])
c2.append(a[1])
porg1_rel.append(pctr[0])
porg2_rel.append(pctr[1])
porg1.append(p[0])
porg2.append(p[1])
df_breakout = pd.DataFrame({'angle':first_angles, 'pc1':c1, 'pc2':c2,
'pc1_org':porg1, 'pc2_org':porg2,
'pc1_rel':porg1_rel, 'pc2_rel':porg2_rel})
df_breakin = []
if break_in:
rem_end = rem_end_nrem
seq = sleepy.get_sequences(np.where(M==1)[0])
idx = state_idx[3]
C = PC[0:2,idx].T
meanc = C.mean(axis=0)
covar = np.cov(C.T)
pca = PCA(n_components=2)
pca.fit(C)
# get the eigenvectors of the covariance matrix:
pc1, pc2 = pca.components_
if pc1[1] < 0:
pc1 = -pc1
w = np.zeros((2,2))
w[:,0] = pc1
w[:,1] = pc2
first = []
for r in rem_end:
ifirst = first_subspace_point(r, PC[0:2,:], meanc, covar, scale=scale)
first.append(ifirst)
plt.plot(PC[0,ifirst], PC[1,ifirst], 'bo')
first_angles = []
c1, c2 = [], []
porg1, porg2 = [], []
porg1_rel, porg2_rel = [], []
for ifirst in first:
# take the first point outside the NREM subspace
p = PC[0:2,ifirst]
# center the point
pctr = p - meanc
# project it onto the eigenvectors
a = np.dot(pctr, w)
x, y = a[0], a[1]
# calculate the angle in radians
theta = math.atan2(x, y)
# Convert the angle to degrees;
# NOTE: 0 deg. corresponds to 3; 90 deg. corresponds to 12
angle_degrees = math.degrees(theta)
first_angles.append(angle_degrees)
c1.append(a[0])
c2.append(a[1])
porg1_rel.append(pctr[0])
porg2_rel.append(pctr[1])
porg1.append(p[0])
porg2.append(p[1])
df_breakin = pd.DataFrame({'angle':first_angles, 'pc1':c1, 'pc2':c2,
'pc1_org':porg1, 'pc2_org':porg2,
'pc1_rel':porg1_rel, 'pc2_rel':porg2_rel})
return ax, df_breakout, df_breakin, df_ampl
def state_space_geometry(PC, M, ma_thr=10, ma_rem_exception=False, kcuts=[], dt=2.5, ax='',
outline_std=True, prefr=True, show_nrem=True, scale=1.645):
"""
(1) Distance between different subspaces
(2) Refractory and permissive state space; draw ellipses capturing the distribution
of the refractory and permissive period. Note that these periods only include
NREM sleep
Returns
-------
df_geom: pd.DataFrame
with columns ['pc1', 'pc2', 'area', 'state'] that
describe for each $state the coordindates of the mean of its subspace
spanned by 'pc1' and 'pc2' and the 'area' of this subspace.
df_distr: pd.DataFrame
with columns ['pc1', 'pc2', 'state']
All PC1 and PC2 values within refractory and permissive state space
"""
state_map = {1:'REM', 2:'Wake', 3:'NREM'}
# KCUTS
tidx = np.arange(0, len(M))
# get the indices (in brainstate time) that we're going to completely discard:
if len(kcuts) > 0:
kidx = []
for kcut in kcuts:
a = int(kcut[0]/dt)
b = int(kcut[-1]/dt)
if b > len(M):
b = len(M)
kidx += list(np.arange(a, b))
tidx = np.setdiff1d(tidx, kidx)
M = M[tidx]
nhypno = len(tidx)
###########################################################################
nhypno = np.min((len(M), PC.shape[1]))
M = M[0:nhypno]
# flatten out MAs
if ma_thr>0:
seq = sleepy.get_sequences(np.where(M==2)[0])
for s in seq:
if np.round(len(s)*dt) <= ma_thr:
if ma_rem_exception:
if (s[0]>1) and (M[s[0] - 1] != 1):
M[s] = 3
else:
M[s] = 3
if ax == '':
plt.figure()
ax = plt.axes([0.2, 0.15, 0.7, 0.7])
clrs = [[0, 0, 0], [0, 1, 1], [0.6, 0, 1], [0.8, 0.8, 0.8], [1, 0.2, 0.2]]
state_idx = {}
for s in [1,2,3]:
idx = np.where(M==s)[0]
state_idx[s] = idx
# get all indices for REM, Wake, and NREM
data_geom = []
if not show_nrem:
shown_states = [1,2]
else:
shown_states = [1,2,3]
for s in shown_states:
idx = state_idx[s]
C = PC[0:2,idx].T
if outline_std:
mean = C.mean(axis=0)
covar = np.cov(C.T)
v, w = linalg.eigh(covar)
# columns of w are the eigenvectors
# v are the eigenvalues in ascending order
# 2 * scale * std (the 2 is to transform the radius to diameter)
v = 2.0 * scale * np.sqrt(v)
u = w[0] / linalg.norm(w[0])
# Plot an ellipse to show the Gaussian component
angle = np.arctan(u[1] / u[0])
angle = 180.0 * angle / np.pi # convert to degrees
ell = mpl.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=clrs[s], lw=2, fill=True)
#else:
# ell = mpl.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=clrs[s], lw=2)
ell.set_clip_box(ax.bbox)
ell.set_alpha(0.3)
ax.add_artist(ell)
area = np.pi * v[0] * v[1]
data_geom += [list(C.mean(axis=0)) + [area] + [state_map[s]]]
else:
sns.kdeplot(x=C[:, 0], y=C[:, 1], ax=ax, color=clrs[s], fill=True, alpha=0.8, levels=[0.25, 0.5, 0.75, 1])
# Add refractory period to state space
if prefr:
refr_color = 'maroon'
perm_color = 'dodgerblue'
rp_clrs = [perm_color, refr_color]
df_refr, refr_vec, _ = add_refr(M)
refr_idx = np.where(refr_vec == 1)[0]
perm_idx = np.where(refr_vec == 2)[0]
nr_idx = np.where(M==3)[0]
nr_idx_refr = np.intersect1d(refr_idx, nr_idx)
nr_idx_perm = np.intersect1d(perm_idx, nr_idx)
Crefr = PC[0:2,nr_idx_refr].T
Cperm = PC[0:2,nr_idx_perm].T
data_distr = []
data_distr += zip(Crefr[:,0], Crefr[:,1], ['refr']*Crefr.shape[0])
data_distr += zip(Cperm[:,0], Cperm[:,1], ['perm']*Cperm.shape[0])
df_distr = pd.DataFrame(data=data_distr, columns=['pc1', 'pc2', 'state'])
if outline_std:
for C,clr,label in zip([Cperm, Crefr], rp_clrs, ['perm', 'refr']):
mean = C.mean(axis=0)
covar = np.cov(C.T)
v, w = linalg.eigh(covar)
# columns of w are the eigenvectors
# v are the eigenvalues in ascending order

neuropyx.py at commit ded0273, under MIT · at the source

Overview

Authors: David E. Lozano¹, Jiso Hong¹, Xi Jin¹, Joseph A. Stucynski¹, Christian K. Machens², Shinjae Chung¹, Franz Weber¹

ORCID iDs: Christian K. Machens, Shinjae Chung, Franz Weber

Institutions: University of Pennsylvania (United States); Champalimaud Foundation (Portugal)

Journal: Nature neuroscience, volume 29, issue 7, pages 1625-1637

Dates: received 8 October 2024; accepted 22 April 2026; published online 25 May 2026; in print 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41593-026-02314-z · PMID 42185574 · PMCID PMC13270127 · OpenAlex W7162341888

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), mouse (organism), systems (subfield)

Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, Single-unit activity, calcium imaging, Physiology & signal measures

Keywords: Circadian rhythms and sleep, Neural circuits

MeSH: Brain Stem*, Neurons*, Sleep, REM*, Animals, Electroencephalography, Male, Mice, Mice, Inbred C57BL (* major topic)

Topic: Sleep and Wakefulness Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: HHS | NIH | National Heart, Lung, and Blood Institute (NHLBI) (R01HL149133, R01HL172283); NCCIH (R01AT013188); NINDS (T32NS105607)

Citations: cited by 1 paper (Europe PMC); 75 references in the paper

Research resources: RRID:Addgene_20298, RRID:Addgene_27056, we injected AAVrg-Ef1a-mCherry-IRES-Cre RRID:Addgene_55632

Abstract

Rapid-eye-movement (REM) sleep is generated in the brainstem, but the brainstem population dynamics that drive transitions to REM sleep remain largely unknown. Here, combining mouse Neuropixels recordings and dimensionality reduction, we found that population activity in the midbrain and pons is dominated by two components, one of which captures strong infraslow fluctuations in neural activity. During transitions from non-REM (NREM) to REM sleep, the population activity followed a stereotypic trajectory that was preceded by an increase in the infraslow component. Our analysis revealed—across all brainstem areas—subpopulations of REM sleep-activated and REM sleep-inhibited neurons with opposing infraslow dynamics and diverging ramping activity between REM sleep episodes, reinforced through antagonistic functional connections. Activation of REM sleep-promoting medullary neurons rapidly enhanced the infraslow component, whose strength gated the ability of upstream circuits to induce REM sleep. Collectively, our results identify a population-level mechanism for gating REM sleep, suggesting that NREM-to-REM sleep transitions are coordinated by low-dimensional, antagonistic brainstem dynamics.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 17 matches between paragraphs and lines of code.

tortugar/Lab

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: bcb8dae1594e64a511545e34f6050e2a417c1f45, 19 March 2025

Languages: Python (33)

Size: 35 files, 33 scripts

Software Heritage: not archived

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (29 files), SciPy (25 files), Matplotlib (19 files), h5py (11 files), pandas (10 files), seaborn (9 files), OpenCV (5 files), Pingouin (2 files), scikit-learn (2 files), statsmodels (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

33 files

Fear/fear.py, Python, 605 lines
Fear/fear_processing.py, Python, 306 lines
Fear/sev.py, Python, 90 lines
MiniscImaging/Utility.py, Python, 361 lines
MiniscImaging/alignsession.py, Python, 523 lines
MiniscImaging/alignsession_script.py, Python, 265 lines
MiniscImaging/data_processing_cl.py, Python, 414 lines
MiniscImaging/imaging.py, Python, 5,179 lines
MiniscImaging/imaging_gui.py, Python, 712 lines
MiniscImaging/roipoly.py, Python, 147 lines
MiniscImaging/sleepy.py, Python, 5,364 lines
Photometry/fib_processing.py, Python, 330 lines
Photometry/pyphi.py, Python, 4,242 lines, 5 matches
Photometry/sev.py, Python, 88 lines
PyScripts/PhasicEvents/spont_rem.py, Python, 142 lines
PyScripts/graph_conditions_stats.py, Python, 1,868 lines
PyScripts/spectralfield_analysis.py, Python, 64 lines
PySleep/data_processing.py, Python, 529 lines
PySleep/sleep_annotation_qt.py, Python, 1,229 lines
PySleep/sleepy.py, Python, 5,349 lines, 2 matches
PySpike/spike_view.py, Python, 1,355 lines
PySpike/spyke.py, Python, 5,619 lines, 4 matches
TDTFear/extract_videoframes.py, Python, 72 lines
TDTFear/run_con.py, Python, 176 lines
TDTPhotometry/fib_processing.py, Python, 327 lines
TDTPhotometry/run_fibpho_twomice.py, Python, 199 lines
TDTPy/run_spikes_recording.py, Python, 179 lines
VideoProcessing/crop_movie.py, Python, 63 lines
VideoProcessing/video_processor_stack.py, Python, 1,096 lines
VideoProcessing/vypro.py, Python, 1,810 lines
VirusExpression/detect_contour.py, Python, 392 lines
VirusExpression/draw_contours.py, Python, 100 lines
VirusExpression/test_contours.py, Python, 27 lines

At the source: github.com/tortugar/Lab/

tortugar/Npx

License: MIT

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: ded027368c37f19c837baf02b81802007886ba46, 16 April 2026

Languages: Jupyter (1), Python (1)

Size: 6 files, 2 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, 1 notebook

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (2 files), NumPy (2 files), seaborn (2 files), h5py (1 file), pandas (1 file), Pingouin (1 file), scikit-learn (1 file), SciPy (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

4 files

basic_analysis_howto.ipynb, Jupyter, 362 lines, 1 match
neuropyx.py, Python, 5,681 lines, 5 matches
LICENSE, License, 21 lines
README.md, Text, 70 lines

At the source: github.com/tortugar/Npx/

Code availability

Code for sleep annotation, analysis of sleep data and EEG/EMG recordings is available on GitHub via https://github.com/tortugar/Lab/. Code for analysis of Neuropixels data is available at https://github.com/tortugar/Npx/.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
35 scripts, each with its path and the digest of its content;
17 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

zenodo:19457135, at Zenodo; found in the references
zenodo:19462601, at Zenodo; found in “Data availability”

Data availability

Neuropixels datasets generated as part of this study are available on Zenodo via 10.5281/zenodo.19462601 (ref. 75). Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 2 keywords, 8 MeSH terms, 3 funders, 75 references, 3 RRIDs.

Cite

This paper

Lozano, D. E., Hong, J., Jin, X., Stucynski, J. A., Machens, C. K., Chung, S., & Weber, F. (2026). Low-dimensional population dynamics in the brainstem gate REM sleep. Nature neuroscience, 29(7), 1625-1637. https://doi.org/10.1038/s41593-026-02314-z

BibTeX

@article{lozano2026low, author = {Lozano, David E. and Hong, Jiso and Jin, Xi and Stucynski, Joseph A. and Machens, Christian K. and Chung, Shinjae and Weber, Franz}, title = {{Low-dimensional population dynamics in the brainstem gate REM sleep}}, journal = {Nature neuroscience}, year = {2026}, month = may, volume = {29}, number = {7}, pages = {1625--1637}, publisher = {Nature Portfolio}, issn = {1097-6256}, doi = {10.1038/s41593-026-02314-z}, url = {https://doi.org/10.1038/s41593-026-02314-z}, pmid = {42185574}, pmcid = {PMC13270127} }

RIS

TY - JOUR AU - Lozano, David E. AU - Hong, Jiso AU - Jin, Xi AU - Stucynski, Joseph A. AU - Machens, Christian K. AU - Chung, Shinjae AU - Weber, Franz TI - Low-dimensional population dynamics in the brainstem gate REM sleep T2 - Nature neuroscience J2 - Nat Neurosci PY - 2026 DA - 2026/05/25 VL - 29 IS - 7 SP - 1625 EP - 1637 SN - 1097-6256 PB - Nature Portfolio DO - 10.1038/s41593-026-02314-z UR - https://doi.org/10.1038/s41593-026-02314-z LA - en ER -

CSL-JSON

{ "id": "10.1038/s41593-026-02314-z", "type": "article-journal", "title": "Low-dimensional population dynamics in the brainstem gate REM sleep", "container-title": "Nature neuroscience", "author": [ { "family": "Lozano", "given": "David E." }, { "family": "Hong", "given": "Jiso" }, { "family": "Jin", "given": "Xi" }, { "family": "Stucynski", "given": "Joseph A." }, { "family": "Machens", "given": "Christian K." }, { "family": "Chung", "given": "Shinjae" }, { "family": "Weber", "given": "Franz" } ], "container-title-short": "Nat Neurosci", "volume": "29", "issue": "7", "page": "1625-1637", "DOI": "10.1038/s41593-026-02314-z", "PMID": "42185574", "PMCID": "PMC13270127", "ISSN": "1097-6256", "publisher": "Nature Portfolio", "URL": "https://doi.org/10.1038/s41593-026-02314-z", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 25 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

Low-dimensional population dynamics in the brainstem gate REM sleep.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

tortugar/Lab

tortugar/Npx

Code availability

Tracing map

Data

Datasets cited

Data availability

Versions

Version 2, 28 September 2026

Version 1, 28 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity