"""Collect rasters of responses for all LGN opto experiments as a function of stimulus type (movies, gratings), run state, and opto state. Then, calculate various measures and save results to DataFrames. Some variables (DataFrames and mseu string lists) can be optionally saved to and loaded from .pickle files to avoid long and unnecessary recalculation times. After connecting to server with `djd`, run this script at the IPython/Jupyter command line; the `-i` runs scripts in IPython's namespace instead of an empty one: `run -i main.py` `pvmvis` (default), `ntsrmvis` or `negntsrmvis` can be provided as an `--exptype=` argument to this script, to load and later plot PV-Cre (most of the paper), Ntsr1-Cre (Fig1S4), or Ntsr negative (Fig1S5) data, respectively. This script can also be run without a server connection, if you choose to load from pickles. """ import sys import argparse from numpy import pi from numpy.random import choice import scipy.stats from scipy.stats import (linregress, ttest_rel, ttest_ind, ttest_1samp, wilcoxon, ks_2samp) from scipy.stats.mstats import gmean from scipy.optimize import least_squares #import statsmodels.formula.api as smf import matplotlib.pyplot as plt from matplotlib.offsetbox import AnchoredText from matplotlib.ticker import ScalarFormatter, FormatStrFormatter import seaborn as sns import pandas as pd #from pycircstat.tests import kuiper from collections import OrderedDict as odict import expio from djd.util import (seq2sql, intround, split_tranges, wrap_raster, issubset, key2mseustr, mseustr2key, key2msestr) from djd.plot import wintitle, simpletraster, cf, saveall, simple_plot_rfs from djd.signal import (raster2psth, raster2freqcomp, sparseness, reliability, snr_baden2016, get_psth_peaks_gac, get_trialmat, pairwisecorr, rfs_for_cdata) from djd.model import fit_MUA_RFs from djd.stats import vector_OSI, percentile_ci, rayleigh_test, niell_DSI from djd.model import sum_of_gaussians from djd.eye import i_clean from util import (get_exps, fitmodel, desat, get_max_snr, ms2msstr, msu2msustr, mseustr2msustr, mseustrs2msustrs, mseustr2mstr, mseustrs2mstrs, mseustr2msestr, mseustrs2msestrs, findmse, load, load_all, save, export2csv, axes_disable_scientific, linear_loss, residual_rsquared) #import rf # deprecated ## Define constants: EXPTYPE = 'pvmvis' # pvmvis, ntsrmvis, negntsrmvis, pgnpvmvis if __name__ == '__main__': # parse command-line arguments: parser = argparse.ArgumentParser(description='Calculate script for natfeedback paper') parser.add_argument("--exptype", help="pvmvis, ntsrmvis, negntsrmvis, pgnpvmvis") args = parser.parse_args() EXPTYPE = args.exptype or EXPTYPE print('*** EXPTYPE: %r' % EXPTYPE) PAPERPATH = '/home/mspacek/blab/natfb/paper' VERBOSE = False # only used for gac at the moment RATETHRESH = 0.01 # mean firing rate threshold for unit inclusion, across all trial types, Hz SNRTHRESH = 0.015 #RSQTHRESH = 0.2 # p-value thresh for significance of each point's deviance from diagonal in meanrate, burst # ratio, and reliability scatter plots: SCATTERPTHRESH = 0.05 SGNF2ALPHA = {True: 1, False: 0.4} # visual drive unit inclusion criteria for gratings, ignores st8 and opto: VDNCRIT = 1 # min num points in tuning curve that are at least VDZCRIT sems away from spont rate VDZCRIT = 2.5 # 1.95, 2.5, 3.2 ## TODO: these need to be optimized for mouse instead of cat: PSTHBASELINEMEDIANX = 2 # PSTH median multiplier to estimate baseline PSTHPEAKMINTHRESH = 3 # PSTH peak detection threshold relative to estimated baseline, Hz STIMTYPES = 'mvi', 'grt' stimtype2axislabel = {'grt':'grating', 'mvi':'movie'} # exclude 'clr' kind, since 'nat' is a superset of 'clr' anyway #MVIKINDS = 'nat', 'pnk', 'shf' MVIKINDS = 'nat', ALLSTIMTYPES = MVIKINDS + ('grt',) STIMTYPESPLUSALL = list(STIMTYPES) + ['mvi+grt'] # strings STIMTYPESPLUSALLI = list(STIMTYPES) + [slice(None)] # for use as indices # iterate over run states in this order: ST8S = ['run', 'sit'] # excludes 'none' ALLST8S = ['none', 'run', 'sit'] ST8COMBOS = [['none'], ['run', 'sit'], ['run', 'sit']] OPTOCOMBOS = [[False, True], [False], [True]] st82crit = odict((('run', '(run_speed > 1).mean() > 0.5'), ('sit', '(run_speed < 0.25).mean() > 0.5'), ('none', 'none'))) rungreen = '#' + ''.join([ format(val, '02x') for val in [66, 168, 117] ]) sitorange = '#' + ''.join([ format(val, '02x') for val in [235, 115, 9] ]) st82clr = {'run':rungreen, 'sit':sitorange, 'none':'black'} st82alpha = {None:1, 'run':1, 'sit':0.5} # sitting is like suppression, so desaturate colour st82clrmap = {'run':plt.cm.Greens, 'sit':plt.cm.Oranges, 'none':plt.cm.Greys} # define example mseus and their colors: green = '#007f00' # deep green violet = '#9f3fff' # pure violet (7f00ff) is a little too dark darkblue = '#000064' # control condition optoblue = '#2a7fff' # light shade of blue for opto condition black = '#000000' asterisk = 5, 2, 0 # 5 sided, asterisk style, 0 deg rotation, has rounded corners though DEFSZ = 50 # default scatter plot point size (open points) exmpli2clr = {1:violet, 2:green, 3:optoblue} # example neurons 1, 2, 3 # use different marker shapes to help distinguish example points: exmpli2mrk = {1:'x', 2:'+', 3:asterisk} exmpli2sz = {1:60, 2:70, 3:60} # marker size exmpli2lw = {1:2, 2:2, 3:2} # (marker) line width # use large closed circles for all example points: #exmpli2mrk = {1:'.', 2:'.', 3:'.'} #EXMPLSZ = 70 # example scatter plot point size (closed points) #exmpli2sz = {1:EXMPLSZ, 2:EXMPLSZ, 3:EXMPLSZ} fig1exmpli = 1 # fig1 uses example neuron 1 fig5exmpli = 1 # fig5 uses example neuron 1 mvimseu2exmpli = {'PVCre_2018_0003_s03_e03_u51':1, # example neuron 1 #'PVCre_2018_0003_s03_e03_u52', 'PVCre_2017_0008_s12_e06_u56':2, # example neuron 2 #'PVCre_2017_0008_s09_e07_u22', } grtmseu2exmpli = {'PVCre_2018_0003_s03_e02_u51':1, # example neuron 1 #'PVCre_2018_0003_s03_e02_u52', #'PVCre_2018_0003_s02_e02_u11', #'PVCre_2018_0003_s02_e02_u38': #'PVCre_2018_0001_s02_e06_u64': #'PVCre_2018_0001_s02_e06_u61': #'PVCre_2018_0001_s02_e06_u64': #'PVCre_2018_0001_s05_e02_u16':3, # old example neuron 3 'PVCre_2018_0003_s05_e02_u164':3, # new example neuron 3 (2019-05-31) } mvimsu2exmpli = { mseustr2msustr(k):v for k, v in mvimseu2exmpli.items() } grtmsu2exmpli = { mseustr2msustr(k):v for k, v in grtmseu2exmpli.items() } msu2exmpli = {**mvimsu2exmpli, **grtmsu2exmpli} # merge the two burstclr = 'red' OPTOS = False, True opto2tuni = {False:0, True:1} # for indexing into tun_params field in Tuning table opto2fb = {None:'', False:'feedback', True:'suppression'} opto2clr = {None:darkblue, False:darkblue, True:optoblue} opto2alpha = {None:1, False:1, True:0.5} # with opto, feedback is removed, so desaturate colour modmeasures = ['meanrate', 'meanrate02', 'meanrate35', 'blankmeanrate', 'blankcondmeanrate', 'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio', 'spars', 'rel', 'meanpkw', 'snr'] modmeasuresnoblankcond = modmeasures.copy() modmeasuresnoblankcond.remove('blankcondmeanrate') modmeasuresnoblankcond.remove('blankcondmeanburstratio') modmeasuresnoblank = modmeasuresnoblankcond.copy() modmeasuresnoblank.remove('blankmeanrate') modmeasuresnoblank.remove('blankmeanburstratio') measure2axislabel = {'meanrate':'FR', 'meanrate02':'FR', 'meanrate35':'FR', 'blankmeanrate':'FR', 'blankcondmeanrate':'FR', 'meanburstratio':'burst ratio', 'blankmeanburstratio':'burst ratio', 'blankcondmeanburstratio':'burst ratio', 'spars':'sparseness', 'rel':'reliability', 'meanpkw':'mean peak width', 'snr':'SNR'} measure2axisunits = {'meanrate':' (spk/s)', 'meanrate02':' (spk/s)', 'meanrate35':' (spk/s)', 'meanpkw':' (s)'} short2longaxislabel = {'FR':'Firing rate'} # restrict FiringPattern to our usual burst definition from Lu et al, 1992: BURSTCRITERION = ('(fp_dtsilent BETWEEN 0.099 and 0.101) AND ' '(fp_dtburst BETWEEN 0.0039 AND 0.0041)') fmodes = ['', 'nb', 'b', 'nr'] # (narrow) firing modes fmode2txt = {'':'all', 'nb':'nonburst', 'b':'burst', 'nr':'nonrand'} # Tuning table keys used for gratings: ivs_order = 'grat_orientation, opto1' tun_model = 'sum_of_gaussians' model = 'threshlin' # linear or threshlin relaverage = 'mean' # reliability averaging method: mean or median DEFAULTFIGURESIZE = 2.8, 2.8 # normal DUMBBELLFIGSIZE = 3.3, 4.4 # dumbbell plots RASTERSCALEX = 0.75 OFFSETS = -0.5, 0.75 # wide raster time offsets, sec OFFSETS02 = 0, -3 # limits to spikes from t=0 to 2 s, assuming 5 sec movies OFFSETS35 = 3, 0 # limits to spikes from t=3 to 5 s NTRIALSMVIGRT = 120 # limit movies to first 120 trials for more direct comparison with gratings PSTHHEIGHT = 3 binw, tres, kernel = 0.02, 0.001, 'gauss' # for calculating PSTHs, sec psthplottres = 0.010 # sec ssx = intround(psthplottres / tres) # subsample factor to get the desired plot tres MINNTRIALSAMPLETHRESH = 10 # ON/OFF/transient classification: PDT = 0.050 # propagation delay time for movie onset (s) WINDT = 0.20 # duration separating start of pre and PDT, and PDT and end of post (s) EXCLTR = 0, 2*PDT # time range around PDT to treat as transient response assert EXCLTR[0] <= PDT <= EXCLTR[1] ## Query the user: msg = ('Constants defined. What next?\n' '0: calculate from scratch, including fits\n' '1: calculate from scratch, but load fits from pickle\n' '2: load from pickles and exit\n' '3: load from pickles and assume connected\n' '>> ') response = int(input(msg)) # load fits from pickle? save heaps of CPU time: LOADFITSFROMPICKLE = False if response == 0 else True if response in [2, 3]: subfolder = EXPTYPE if EXPTYPE != 'pvmvis' else '' name2val = load_all(subfolder=subfolder) locals().update(name2val) if response == 2: sys.exit() ## Collect series, experiments, units from database: ''' # all single movie opto experiments in LGN that have been sorted: nats = ((e & 'e_name LIKE "MAS_%%"' & 'e_name NOT LIKE "%%+%%"' & 'e_optowl = 470') & (s & 's_depth > 2500' & 's_depth < 3500') # upper limit excludes TRN series & u) # all pink noise opto experiments in LGN that have been sorted: pnks = ((e & 'e_name LIKE "MAS_%%PNK%%"' & 'e_optowl = 470') & (s & 's_depth > 2500' & 's_depth < 3500') & u) ''' """ PVCre mice: * PVCre_2018_0002_s04 - good RFs, OK tuning, excellent opto ** PVCre_2019_0001_s05 - excellent RFs, OK tuning, good opto - make sure to exclude e01 and e02 due to seizure between e02 and e03 ** PVCre_2019_0001_s06 - excellent RFs and tuning, OK opto, started by Davide - sorted: single unit RFs are quite bad though, excluding - this mouse's pupil pointed very high under visible light illumination, in the dark with dilated pupil, it pointed forward normally. Another good reason to exclude - PVCre_2019_0001_s07 - excellent RFs and tuning, very weak opto * PVCre_2019_0001_s08 - excellent RFs and tuning, OK opto. stained with DiI * PVCre_2019_0002_s05 - excellent RFs and tuning, good opto * PVCre_2019_0002_s06 - excellent RFs, OK tuning, excellent opto * PVCre_2019_0002_s07 - excellent RFs, good tuning and opto ** PVCre_2019_0002_s08 - excellent RFs, tuning and opto! stained with DiI - sorted, single unit RFs are good As of 2018-12-13, both PVCre_2018_0010_s02 and s05 are also excluded via updated 's_region' values. Both were intended to be in dLGN, but may have somehow missed it, based on weird RFs Ntsr1cre mice: * Ntsr1Cre_2019_0002_s03: sorted, mov rasters look wrong (bad sorting? asleep mouse?), 100 trials/cond, no running, but opto effects seem consistent - check itracking movies, resort this one - Ntsr1Cre_2019_0002_s04: no mov * Ntsr1Cre_2019_0002_s05: sorted but only 3 units, only 2 have mov PSTH peaks, ~ consistent - Ntsr1Cre_2019_0003: no mov - Ntsr1Cre_2019_0004 has big cortical hole: - Ntsr1Cre_2019_0004_s02: nice mov PSTHs, but no opto effect at all for mov or grt - Ntsr1Cre_2019_0004_s03: ?? ** Ntsr1Cre_2019_0007_s04 - good - Ntsr1Cre_2019_0007_s05 - suspicious about if it's really LGN, SU results look poor ** Ntsr1Cre_2019_0008_s05 - good - do older mice just for their grating data, do shell/core and other separation """ mvis = get_exps(locals(), EXPTYPE) mviseries = s & mvis mviunits = up & mvis mvispikes = spk & mvis if EXPTYPE == 'negntsrmvis': # force inclusion of a V1 movie recording mvieye = EyeIxtract() & [{'m':'Ntsr1Cre_2020_0001', 's':2, 'e':11}, # V1 {'m':'Ntsr1Cre_2020_0001', 's':4, 'e':10}] # dLGN else: mvieye = EyeIxtract() & mvis mvirun = rt & mvis # fetch unique animals included in data set: #np.unique(mviseries.fetch('m')) ## TODO: need to check which units are plausibly in dLGN, not all units on the shank will be! ## - this is mostly done ahead of time during spike sorting, MUA envelope RF analysis suggests ## some single units have maxchans that are outside the chan range of significant RFs of ## the population, though there are some caveats to that - maybe maxchan is too stringent, ## or maybe it's possible for one or two single units to have visual responses even when ## the population doesn't show much response to sparse noise # all grating oritun opto experiments in LGN that have been sorted that occurred during the # same series as movie experiments: grts = ((e & 'e_name LIKE "%%oritun%%"' & 'e_optowl IN (465, 470)') & mviseries & u) grtseries = mviseries & grts # use the same units for gratings as in movies: #grtunits = mviunits & grts # use units active during a grating experiment, not necessarily the same units as in movies: grtunits = up & grts grtspikes = spk & grts grttun = tun & grts spons = ((e & 'e_name LIKE "%%spont%%"' & 'e_optowl IN (465, 470)') & mviseries & u) # get all possible movie mseus, regardless of rate: mvimseus = [] for mse in mvis.fetch(dj.key): uids = (up & mse).fetch('u') for uid in uids: mseu = mse.copy() mseu['u'] = uid mvimseus.append(mseu) # get all possible grating mseus, regardless of rate: grtmseus = [] for mse in grts.fetch(dj.key): uids = (up & mse).fetch('u') for uid in uids: mseu = mse.copy() mseu['u'] = uid grtmseus.append(mseu) ## TODO: grtmseu filtering by visual responsiveness should be done here, but would currently ## result in index errors in the big grtresp loop below, so it's done piecemeal in the loop: #grtmseus = (grts * up).fetch(dj.key) #grtmseus = (grts * up) & (vd.Units() & 'n_crit=1' & 'z_crit > 2.5') # get all possible spontaneous mseus, regardless of rate: sponmseus = [] for mse in spons.fetch(dj.key): uids = (up & mse).fetch('u') for uid in uids: mseu = mse.copy() mseu['u'] = uid sponmseus.append(mseu) msefmt = '{m}_s{s:02}_e{e:02}' msufmt = '{m}_s{s:02}_u{u:02}' mseufmt = '{m}_s{s:02}_e{e:02}_u{u:02}' mvimseustrs = [ mseufmt.format(**mseu) for mseu in mvimseus ] grtmseustrs = [ mseufmt.format(**mseu) for mseu in grtmseus ] sponmseustrs = [ mseufmt.format(**mseu) for mseu in sponmseus ] mvimsustrs = [ msufmt.format(**unit) for unit in mviunits ] grtmsustrs = [ msufmt.format(**unit) for unit in grtunits ] mvigrtmsustrs = list(np.union1d(mvimsustrs, grtmsustrs)) # superset of the two lists if response == 3: sys.exit() ## Calculate: ## collect movie rasters, meanrates, burst info, PSTHs, sparseness, reliability and PSTH ## peak info by mseu, movie kind, run state, and opto condition: levels = [mvimseustrs, MVIKINDS, ALLST8S, OPTOS] names = ['mseu', 'kind', 'st8', 'opto'] mi = pd.MultiIndex.from_product(levels, names=names) columns = ['dt', 'optotrange', 'trialis', 'raster', 'braster', 'nbraster', 'wraster', 'wbraster', 'wnbraster', # rasters 'rates', 'rate02s', 'rate35s', 'blankrates', # single trial FRs 'meanrate', 'meanrate02', 'meanrate35', 'blankmeanrate', # trial-averaged FRs 'burstis', 'wburstis', 'burstratios', 'blankburstratios', # single trial BRs 'meanburstratio', 'blankmeanburstratio', # trial-averaged BRs 'psth', 'bpsth', 'nbpsth', 'wpsth', 'wbpsth', 'wnbpsth', 't', 'wt', 'bins', 'wbins', 'spars', 'rhos', 'rel', 'snr', 'pkts', 'pkws', 'meanpkw'] mviresp = pd.DataFrame(index=mi, columns=columns) ntotpeaks = 0 # total number of peaks detected for mvi in mvis: msestr = msefmt.format(**mvi) print(msestr) exp = mvi['e'] moviecond = mov & mvi kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping trials = trial & mvi alltrialis = np.sort(trials.fetch('trial_id')) # trialis corresponding to rows in tranges #rates = ra & mvi uids = (up & mvi).fetch('u') # all trials, in trial ID order: rasters, tranges, _, opto, _ = (spk & mvi).get_rasters() wrasters, *_ = (spk & mvi).get_rasters(offsets=OFFSETS) # w = wide, same returned tranges rasters02, *_ = (spk & mvi).get_rasters(offsets=OFFSETS02) # from 0 to 2 s, same tranges rasters35, *_ = (spk & mvi).get_rasters(offsets=OFFSETS35) # from 3 to 5 s, same tranges wtranges = tranges.copy() + OFFSETS #tranges02 = tranges.copy() + OFFSETS02 #tranges35 = tranges.copy() + OFFSETS35 # make sure trial duration is very consistent: dts = tranges.ptp(axis=1) dt = dts.mean() # mean trial duration dtstd = dts.std() assert dtstd / dt < 0.005 # make sure opto tranges are consistent relative to trial tranges: optotrialstimt0s = tranges[opto == True][:, 0] evtranges = (evtimes & mvi).fetch1('ev_tranges') alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1)) optotrange = alloptotranges.mean(axis=0) optotrangestd = alloptotranges.std(axis=0) assert (optotrangestd / optotrange < 0.005).all() # get tranges and rasters for just the period before stimulus onset during which opto # can be on, i.e. from optotrange[0] (-0.272 s) to 0 s before each stimulus onset: blankdt = abs(optotrange[0]) # averages to 0.283 s across all experiments print('mvi blankdt:', blankdt) blanktranges = np.zeros_like(tranges) blanktranges[:, 0] = tranges[:, 0] - blankdt blanktranges[:, 1] = tranges[:, 0] blankrasters, _, _, _, _ = (spk & mvi).get_rasters(tranges=blanktranges) mviopto = evcond & mvi & 'ev_chan="opto1"' mvioptostimis, mvioptovals = mviopto.fetch('stim_id', 'ec_val') # calculate burst ratios for all units: firingpatterns = fp & mvi & BURSTCRITERION brus = firingpatterns.burst_ratio(tranges=tranges) wbrus = firingpatterns.burst_ratio(tranges=wtranges) blankbrus = firingpatterns.burst_ratio(tranges=blanktranges) for kind, kindstimis in kind2stimis.items(): print(kind) if not kindstimis or kind not in MVIKINDS: # exclude empty stimis and 'clr' kind continue kindtrials = (trials & 'stim_id IN %s' % seq2sql(kindstimis)) kindtrialis = kindtrials.fetch('trial_id') trialiis = np.isin(alltrialis, kindtrialis) # bool with same shape as alltrialis # calculate narrow and wide PSTH bins for all units: tmin, tmax = 0, dt bins = split_tranges([(tmin, tmax)], binw, tres) # narrow midbins = bins.mean(axis=1) # narrow wbins = split_tranges([(tmin+OFFSETS[0], tmax+OFFSETS[1])], binw, tres) # wide wmidbins = wbins.mean(axis=1) # wide for st8, st8crit in st82crit.items(): if st8crit == 'none': # 'none' is not in the State.Trial table st8trialis = alltrialis else: st8trialis = (st8t & mvi & {'st8_crit':st8crit}).fetch('trial_id') st8trialis = np.sort(st8trialis) for opto in OPTOS: optois = mvioptovals == opto # limit stimis by opto: optostimis = mvioptostimis[optois] # also limit stimis by movie kind: stimis = np.intersect1d(kindstimis, optostimis) if len(stimis) == 0: continue optotrials = trials & 'stim_id IN %s' % seq2sql(optostimis) optotrialis = optotrials.fetch('trial_id') # intersect kind, st8 and opto trialis: trialis = np.intersect1d(kindtrialis, st8trialis) trialis = np.intersect1d(trialis, optotrialis) trialiis = np.isin(alltrialis, trialis) # bool with same shape as alltrialis if len(trialis) == 0: print('%s %s %s has no trials, skipping' % (kind, st8, opto)) continue for uid in uids: if uid not in rasters: continue mseu = {'m':mvi['m'], 's':mvi['s'], 'e':mvi['e'], 'u':uid} mseustr = mseufmt.format(**mseu) raster = rasters[uid][trialiis] # subset of raster trials wraster = wrasters[uid][trialiis] # subset of wraster trials blankraster = blankrasters[uid][trialiis] # subset of blankraster trials raster02 = rasters02[uid][trialiis][:NTRIALSMVIGRT] raster35 = rasters35[uid][trialiis][:NTRIALSMVIGRT] rates = np.array([ len(trialspikes)/dt for trialspikes in raster ]) rate02s = np.full(rates.shape, np.nan) # init with full shape with NaN rate35s = np.full(rates.shape, np.nan) blankrates = np.full(rates.shape, np.nan) # fill only the first n trials/blanks: rate02s[:len(raster02)] = np.array([ len(trialspikes)/2 for trialspikes in raster02 ]) rate35s[:len(raster35)] = np.array([ len(trialspikes)/2 for trialspikes in raster35 ]) blankrates[:len(blankraster)] = np.array([ len(bspikes)/blankdt for bspikes in blankraster ]) meanrate = len(np.hstack(raster)) / dt / len(raster) #assert np.isclose(meanrate, rates.mean()) # sanity check meanrate02 = len(np.hstack(raster02)) / 2 / len(raster02) meanrate35 = len(np.hstack(raster35)) / 2 / len(raster35) blankmeanrate = len(np.hstack(blankraster)) / blankdt / len(blankraster) if meanrate < RATETHRESH: print('%s %s %s %s does not meet RATETHRESH, skipping' % (mseustr, kind, st8, opto)) continue mviresprow = mviresp.loc[mseustr, kind, st8, opto] # technically, this is the wrong approach, and could result in writing # to a copy (an independent Series), which isn't the intention (see https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy) # but it seems to work because the mviresp DataFrame is being filled one # row at a time, and wasn't initialized explicitly by setting everything # to e.g. np.nan in a single vector operation mviresprow['dt'] = dt mviresprow['optotrange'] = optotrange mviresprow['trialis'] = trialis mviresprow['raster'] = raster # narrow raster mviresprow['wraster'] = wraster # wide raster, including OFFSETS mviresprow['rates'] = rates # based on narrow raster mviresprow['rate02s'] = rate02s mviresprow['rate35s'] = rate35s mviresprow['blankrates'] = blankrates mviresprow['meanrate'] = meanrate # based on narrow raster mviresprow['meanrate02'] = meanrate02 mviresprow['meanrate35'] = meanrate35 mviresprow['blankmeanrate'] = blankmeanrate # fill in the corresponding trial burst ratios: for bru, wbru, blankbru in zip(brus, wbrus, blankbrus): if bru['u'] == uid: assert bru['m'] == mseu['m'] assert bru['s'] == mseu['s'] assert bru['e'] == mseu['e'] burstis = np.array(bru['burst_index'], # narrow, ragged array dtype=object)[trialiis] wburstis = np.array(wbru['burst_index'], # wide dtype=object)[trialiis] #blankburstis = np.array(blankbru['burst_index'], # dtype=object)[trialiis] burstratios = np.array(bru['burst_ratio'])[trialiis] # narrow blankburstratios = np.array(blankbru['burst_ratio'])[trialiis] mviresprow['burstis'] = burstis # narrow mviresprow['wburstis'] = wburstis # wide, for plotting mviresprow['burstratios'] = burstratios # narrow mviresprow['blankburstratios'] = blankburstratios mviresprow['meanburstratio'] = np.nanmean(burstratios) mviresprow['blankmeanburstratio'] = np.nanmean(blankburstratios) break # out of bru loop # collect narrow burst and non-burst spike rasters: braster, nbraster = [], [] for row, bis in zip(raster, burstis): # iterate over trials in raster if len(bis) == 0: bis = [] # none nbis = slice(None) # all else: allis = np.arange(len(row)) # all spikes in this trial nbis = np.setdiff1d(allis, bis) braster.append(row[bis]) nbraster.append(row[nbis]) # collect wide burst and non-burst spike rasters: wbraster, wnbraster = [], [] for row, bis in zip(wraster, wburstis): # iterate over trials in wraster if len(bis) == 0: bis = [] # none nbis = slice(None) # all else: allis = np.arange(len(row)) # all spikes in this trial nbis = np.setdiff1d(allis, bis) wbraster.append(row[bis]) wnbraster.append(row[nbis]) braster = np.array(braster, dtype=object) # ragged array nbraster = np.array(nbraster, dtype=object) wbraster = np.array(wbraster, dtype=object) wnbraster = np.array(wnbraster, dtype=object) mviresprow['braster'] = braster # narrow burst raster mviresprow['nbraster'] = nbraster # narrow non-burst raster mviresprow['wbraster'] = wbraster # wide burst raster mviresprow['wnbraster'] = wnbraster # wide non-burst raster # calculate narrow PSTHs: psth = raster2psth(raster, bins, binw, tres, kernel) bpsth = raster2psth(braster, bins, binw, tres, kernel) nbpsth = raster2psth(nbraster, bins, binw, tres, kernel) # calculate wide PSTHs: wpsth = raster2psth(wraster, wbins, binw, tres, kernel) wbpsth = raster2psth(wbraster, wbins, binw, tres, kernel) wnbpsth = raster2psth(nbraster, wbins, binw, tres, kernel) mviresprow['psth'] = psth # narrow PSTH mviresprow['bpsth'] = bpsth # narrow burst PSTH mviresprow['nbpsth'] = nbpsth # narrow non-burst PSTH mviresprow['wpsth'] = wpsth # wide PSTH mviresprow['wbpsth'] = wbpsth # wide burst PSTH mviresprow['wnbpsth'] = wnbpsth # wide non-burst PSTH mviresprow['t'] = midbins # narrow PSTH timepoints mviresprow['wt'] = wmidbins # wide PSTH timepoints mviresprow['bins'] = bins # narrow PSTH bins mviresprow['wbins'] = wbins # wide PSTH bins # calculate PSTH sparseness on narrow PSTH: spars = sparseness(psth) mviresprow['spars'] = spars # collect single trial signals from narrow PSTH, calc response reliability: trialsignals = [] for trialspiketimes in raster: trialsignal = raster2psth([trialspiketimes], bins, binw, tres, kernel) trialsignals.append(trialsignal) trialsignals = np.vstack(trialsignals) rel, rhos = reliability(trialsignals, average=relaverage) if not np.isnan(rel): assert rel > -0.05 # can come out weakly negative mviresprow['rhos'] = rhos mviresprow['rel'] = max(rel, 0) # clamp weakly negative values to 0 mviresprow['snr'] = snr_baden2016(trialsignals) # detect peaks in narrow PSTH: baseline = PSTHBASELINEMEDIANX * np.median(psth) thresh = baseline + PSTHPEAKMINTHRESH # peak detection threshold spikets = np.sort(np.concatenate(raster)) # flatten raster spike times peakis, lis, ris = get_psth_peaks_gac(spikets, midbins, psth, thresh, verbose=VERBOSE) if VERBOSE: print() # new line npeaks = len(peakis) if npeaks > 0: mviresprow['pkts'] = peakis * tres # PSTH peak times (s) pkws = (ris - lis) * tres # PSTH peak widths (s) mviresprow['pkws'] = pkws mviresprow['meanpkw'] = pkws.mean() # mean PSTH peak width (s) ## TODO: use the peakis themselves for some other kind of measure, ## maybe plot distributions of peaks relative to stimulus start, or ## interpeak intervals, e.g. # use absence of PSTH peaks to decide which units/conditions to exclude? # this would obviously only work for movies, but that's what mviresp # is all about... #if npeaks == 0: # nrnids[state].append(nid) # save nonresponsive nids by state # continue # this PSTH has no peaks, skip all subsequent measures ntotpeaks += npeaks print('Detected %d movie PSTH peaks' % ntotpeaks) ## collect grating meanrates, rasters, tuning, and burst info by mseu, ## run state, and opto condition: levels = [grtmseustrs, ALLST8S, OPTOS] names = ['mseu', 'st8', 'opto'] mi = pd.MultiIndex.from_product(levels, names=names) columns = ['dt', 'optotrange', 'trialis', 'raster', 'braster', 'nbraster', 'wraster', # rasters 'rates', 'blankrates', 'blankcondrates', # single trial FRs 'meanrate', 'blankmeanrate', 'blankcondmeanrate', # trial-averaged FRs 'burstis', 'wburstis', 'burstratios', # single trial BRs 'blankburstratios', 'blankcondburstratios', # single trial BRs 'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio', # trial-averaged BRs 'stimis', 'oris', 'tfreq', 'tun', 'tunrsq'] grtresp = pd.DataFrame(index=mi, columns=columns) for grt in grts: msestr = msefmt.format(**grt) print(msestr) exp = grt['e'] alltrials = trial & grt # all trials, including blanks # all stimis and trialis, in stimis order: allstimis, alltrialis = alltrials.fetch('stim_id', 'trial_id') # map triali to stimi for all trials: triali2stimi = { triali:stimi for stimi, triali in zip(allstimis, alltrialis) } alltrialis = np.sort(alltrialis) # trialis corresponding to rows in tranges trials = alltrials & (grat & 'grat_contrast != 0.0') # filter out blank trials blankcondtrials = alltrials & (grat & 'grat_contrast = 0.0') # blank trials only uids = (up & grt).fetch('u') # all trials, in trial ID (temporal) order: rasters, tranges, _, opto, _ = (spk & grt).get_rasters() wrasters, *_ = (spk & grt).get_rasters(offsets=OFFSETS) # w = wide, same returned tranges wtranges = tranges.copy() + OFFSETS # make sure opto tranges are consistent relative to trial tranges: optotrialstimt0s = tranges[opto == True][:, 0] evtranges = (evtimes & grt).fetch1('ev_tranges') alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1)) optotrange = alloptotranges.mean(axis=0) optotrangestd = alloptotranges.std(axis=0) assert (optotrangestd / optotrange < 0.005).all() # get tranges and rasters for just the period before stimulus onset during which opto # can be on, i.e. from optotrange[0] to 0 s before each stimulus onset: blankdt = abs(optotrange[0]) # averages to 0.285 s across all experiments, as in movies print('grt blankdt:', blankdt) blanktranges = np.zeros_like(tranges) blanktranges[:, 0] = tranges[:, 0] - blankdt blanktranges[:, 1] = tranges[:, 0] blankrasters, _, _, _, _ = (spk & grt).get_rasters(tranges=blanktranges) grtopto = evcond & grt & 'ev_chan="opto1"' grtoptostimis, grtoptovals = grtopto.fetch('stim_id', 'ec_val') # calculate burst ratios for all units: firingpatterns = fp & grt & BURSTCRITERION brus = firingpatterns.burst_ratio(tranges=tranges) wbrus = firingpatterns.burst_ratio(tranges=wtranges) blankbrus = firingpatterns.burst_ratio(tranges=blanktranges) # make sure trial duration is very consistent: dts = tranges.ptp(axis=1) dt = dts.mean() # mean trial duration, for all trials (non-blank and blank) dtstd = dts.std() assert dtstd / dt < 0.005 # make sure opto tranges are consistent relative to trial tranges: optotrialstimt0s = tranges[opto == True][:, 0] evtranges = (evtimes & grt).fetch1('ev_tranges') alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1)) optotrange = alloptotranges.mean(axis=0) optotrangestd = alloptotranges.std(axis=0) assert (optotrangestd / optotrange < 0.005).all() # make sure temporal freq is very consistent: tfreqs = (grat & grt).fetch('grat_temp_freq') tfreq = tfreqs.mean() tfreqstd = tfreqs.std() assert tfreqstd / tfreq < 0.005 for st8, st8crit in st82crit.items(): if st8crit == 'none': # 'none' is not in the State.Trial table st8trialis = alltrialis else: st8trialis = (st8t & grt & {'st8_crit':st8crit}).fetch('trial_id') st8trialis = np.sort(st8trialis) for opto in OPTOS: optois = grtoptovals == opto # limit stimis by opto: optostimis = grtoptostimis[optois] optotrials = trials & 'stim_id IN %s' % seq2sql(optostimis) blankcondoptotrials = blankcondtrials & 'stim_id IN %s' % seq2sql(optostimis) optotrialis = optotrials.fetch('trial_id') blankcondoptotrialis = blankcondoptotrials.fetch('trial_id') # intersect st8 and opto trialis: trialis = np.intersect1d(st8trialis, optotrialis) blankcondtrialis = np.intersect1d(st8trialis, blankcondoptotrialis) if len(trialis) == 0: print('%s %s has no trials, skipping' % (st8, opto)) continue stimis = np.array([ triali2stimi[triali] for triali in trialis ]) # triali order oristimis, oris = (grat & grt).fetch('stim_id', 'grat_orientation') oris = np.hstack([ oris[oristimis == stimi] for stimi in stimis ]) # triali order stimsortis = stimis.argsort(kind='stable') trialis = trialis[stimsortis] # sorted by stimi stimis = stimis[stimsortis] # sorted by stimi oris = oris[stimsortis] # sorted by stimi for uid in uids: if uid not in rasters: continue mseu = {'m':grt['m'], 's':grt['s'], 'e':grt['e'], 'u':uid} mseustr = mseufmt.format(**mseu) ## TODO: this should really be done once, where grtmseus are defined, ## but I don't want to break things right now (2019-07-15) if not (vd.Units() & mseu & ('n_crit=%d' % VDNCRIT) & ('z_crit > %f' % VDZCRIT)): # not a visually responsive mseu (this query ignores st8 and opto) print('%s is not visually responsive, skipping' % mseustr) continue raster = rasters[uid][trialis] # subset of trials, sorted by stimi wraster = wrasters[uid][trialis] # subset of trials, sorted by stimi blankraster = blankrasters[uid][trialis] # subset of trials, sorted by stimi blankcondraster = rasters[uid][blankcondtrialis] # subset of blank condition trials rates = np.array([ len(trialspikes)/dt for trialspikes in raster ]) blankrates = np.array([ len(bspikes)/blankdt for bspikes in blankraster ]) assert len(blankrates) == len(rates) # should be same shape as rates blankcondrates = np.array([ len(bcspikes)/dt for bcspikes in blankcondraster ]) meanrate = len(np.hstack(raster)) / dt / len(raster) #assert np.isclose(meanrate, rates.mean()) # sanity check blankmeanrate = len(np.hstack(blankraster)) / blankdt / len(blankraster) if len(blankcondraster) == 0: blankcondmeanrate = np.nan else: blankcondmeanrate = len(np.hstack(blankcondraster)) / dt / len(blankcondraster) if meanrate < RATETHRESH: print('%s %s %s does not meet RATETHRESH, skipping' % (mseustr, st8, opto)) continue grtresprow = grtresp.loc[mseustr, st8, opto] grtresprow['dt'] = dt grtresprow['optotrange'] = optotrange # overload 'trialis' name in df - includes non-blank trials and blank condition # trials, in that order: grtresprow['trialis'] = np.concatenate([trialis, blankcondtrialis]) grtresprow['raster'] = raster # sorted by stimi grtresprow['wraster'] = wraster grtresprow['rates'] = rates grtresprow['blankrates'] = blankrates grtresprow['blankcondrates'] = blankcondrates grtresprow['meanrate'] = meanrate grtresprow['blankmeanrate'] = blankmeanrate grtresprow['blankcondmeanrate'] = blankcondmeanrate grtresprow['stimis'] = stimis # one entry per row in raster grtresprow['oris'] = oris # one entry per row in raster grtresprow['rates'] = rates grtresprow['tfreq'] = tfreq tunrow = (tun & mseu & {'ivs_order':ivs_order, 'tun_model':tun_model, 'st8_crit':st8crit}) if tunrow: # not empty tunparams = tunrow.fetch1('tun_pars') # 2x5 array of tuning params tunrsq = tunrow.fetch1('tun_rsq') # len 2 array of R2 values tunspon = tunrow.fetch1('tun_spon_mean') # len 2 array of spon activity oti = opto2tuni[opto] # opto Tuning index params = list(tunparams[oti]) + [(tunspon[oti])] grtresprow['tun'] = params # len 6 array: 5 tuning params + spon """[dp, rp, rn, r0, sigma, spon]: oripref, gauss ampl @ oripref, gauss ampl @ antipref, gauss offset, gauss width, spon level so index 1 wrt 3 and 5 is a measure of ori tuning strength""" grtresprow['tunrsq'] = tunrsq[oti] else: print('WARNING: empty Tuning row for %s, %s' % (mseu, st8crit)) # fill in the corresponding trial burst ratios: for bru, wbru, blankbru in zip(brus, wbrus, blankbrus): if bru['u'] == uid: assert bru['m'] == mseu['m'] assert bru['s'] == mseu['s'] assert bru['e'] == mseu['e'] burstis = np.array(bru['burst_index'], # narrow, ragged array dtype=object)[trialis] wburstis = np.array(wbru['burst_index'], # wide dtype=object)[trialis] burstratios = np.array(bru['burst_ratio'])[trialis] # narrow blankburstratios = np.array(blankbru['burst_ratio'])[trialis] blankcondburstratios = np.array(bru['burst_ratio'])[blankcondtrialis] grtresprow['burstis'] = burstis # narrow grtresprow['wburstis'] = wburstis # wide, for plotting grtresprow['burstratios'] = burstratios # narrow grtresprow['blankburstratios'] = blankburstratios # narrow grtresprow['blankcondburstratios'] = blankcondburstratios # narrow # nanmean can raise "Mean of empty slice" warnings if argument # is empty or is all nan, but that's OK, returns a nan in both cases: grtresprow['meanburstratio'] = np.nanmean(burstratios) grtresprow['blankmeanburstratio'] = np.nanmean(blankburstratios) if len(blankcondburstratios) == 0: grtresprow['blankcondmeanburstratio'] = np.nan else: grtresprow['blankcondmeanburstratio'] = np.nanmean(blankcondburstratios) break # out of bru loop # collect narrow burst and non-burst spike rasters: braster, nbraster = [], [] for row, bis in zip(raster, burstis): # iterate over trials in raster if len(bis) == 0: bis = [] # none nbis = slice(None) # all else: allis = np.arange(len(row)) # all spikes in this trial nbis = np.setdiff1d(allis, bis) braster.append(row[bis]) nbraster.append(row[nbis]) braster = np.array(braster, dtype=object) # ragged array nbraster = np.array(nbraster, dtype=object) grtresprow['braster'] = braster # narrow burst raster grtresprow['nbraster'] = nbraster # narrow non-burst raster ## Generate grating ori tuning dataframe for saving to a .pickle, and for convenient ## offline tuning plots: grttundf = grttun.fetch(format='frame') grttuninfo = (TuningInfo() & grttun).fetch(format='frame') # contains stimulus oris # convert separate m, s, e, u columns in grttundf to a single mseu column in a new temporary # df without a MultiIndex: grttunresprows = [] for rowindex, rowvalues in grttundf.iterrows(): mid, sid, eid, ivs_order, tun_model, st8_type, st8_crit, uid = rowindex mean, sem, spon_mean, pars, rsq = (rowvalues['tun_mean'], rowvalues['tun_sem'], rowvalues['tun_spon_mean'], rowvalues['tun_pars'], rowvalues['tun_rsq']) mseustr = key2mseustr({'m':mid, 's':sid, 'e':eid, 'u':uid}) ori = grttuninfo.loc[mid, sid, eid, ivs_order]['ti_axes'][:, 0] newrow = {'mseu':mseustr, 'ivs_order':ivs_order, 'tun_model':tun_model, 'st8_type':st8_type, 'st8_crit':st8_crit, 'ori':ori, 'mean':mean, 'sem':sem, 'spon_mean':spon_mean, 'pars':pars, 'rsq':rsq} grttunresprows.append(newrow) newindex = ['mseu', 'ivs_order', 'tun_model', 'st8_type', 'st8_crit'] grttunresp = pd.DataFrame(grttunresprows).set_index(newindex) ## collect spontaneous meanrates and meanburstratios by mseu and opto cond, ## ignore run st8 for now: levels = [sponmseustrs, OPTOS] names = ['mseu', 'opto'] mi = pd.MultiIndex.from_product(levels, names=names) columns = ['meanrate', 'meanburstratio'] sponresp = pd.DataFrame(index=mi, columns=columns) for spon in spons: msestr = msefmt.format(**spon) print(msestr) #exp = spon['e'] optotranges = (evtimes & spon).fetch1('ev_tranges') # ensure that intervals between end and start of consecutive opto pulses are at least # as long as the pulse length: optodts = np.diff(optotranges, axis=1) # opto pulse durations (sec) maxoptodt = optodts.max() # sec isitranges = np.vstack([optotranges[:-1, 0], optotranges[1:, 1]]).T # shape: npulses-1, 2 assert np.diff(isitranges, axis=1).min() > maxoptodt # sec # use one opto pulse duration before each pulse as the control trange for that pulse: ctrltranges = optotranges - optodts ctrldts = np.diff(ctrltranges, axis=1) # control trange durations (sec) spkspon = spk & spon fpspon = fp & spon & BURSTCRITERION optorasters, _, _, optovals, _ = spkspon.get_rasters(optotranges, offsets=[0, 0]) ctrlrasters, _, _, ctrlvals, _ = spkspon.get_rasters(ctrltranges, offsets=[0, 0]) optobrus = fpspon.burst_ratio(tranges=optotranges) ctrlbrus = fpspon.burst_ratio(tranges=ctrltranges) assert optovals.all() == True assert ctrlvals.all() == False assert list(optorasters) == list(ctrlrasters) # should have exact same set of units uids = list(optorasters) opto2rasters = {False:ctrlrasters, True:optorasters} opto2brus = {False:ctrlbrus, True:optobrus} opto2totaldt = {False:ctrldts.sum(), True:optodts.sum()} for opto in OPTOS: rasters = opto2rasters[opto] brus = opto2brus[opto] totaldt = opto2totaldt[opto] # sec for uid, bru in zip(uids, brus): trialspikets = np.concatenate(rasters[uid]) # concatenate spikes from all "trials" mseu = {'m':spon['m'], 's':spon['s'], 'e':spon['e'], 'u':uid} mseustr = mseufmt.format(**mseu) sponresprow = sponresp.loc[mseustr, opto] nspikes = len(trialspikets) meanrate = nspikes / totaldt sponresprow['meanrate'] = meanrate assert bru['u'] == uid # sanity check burstratios = bru['burst_ratio'] # one per "trial" meanburstratio = np.nanmean(burstratios) sponresprow['meanburstratio'] = meanburstratio ## find best movie and grating response parameters per unit, across experiments: mvilevels = [mvimsustrs, MVIKINDS, ALLST8S, OPTOS] grtlevels = [grtmsustrs, ALLST8S, OPTOS] mvinames = ['msu', 'kind', 'st8', 'opto'] grtnames = ['msu', 'st8', 'opto'] mvimi = pd.MultiIndex.from_product(mvilevels, names=mvinames) grtmi = pd.MultiIndex.from_product(grtlevels, names=grtnames) bestmvicolumns = modmeasuresnoblankcond bestgrtcolumns = ['meanrate', 'blankmeanrate', 'meanburstratio', 'blankmeanburstratio'] bestmviresp = pd.DataFrame(index=mvimi, columns=bestmvicolumns) bestgrtresp = pd.DataFrame(index=grtmi, columns=bestgrtcolumns) stimtype2resp = {'mvi':mviresp, 'grt':grtresp} stimtype2bestresp = {'mvi':bestmviresp, 'grt':bestgrtresp} for stimtype in STIMTYPES: resp, bestresp = stimtype2resp[stimtype], stimtype2bestresp[stimtype] columns = bestresp.columns for index, row in resp.iterrows(): mseu = index[0] msu = mseustr2msustr(mseu) bestindex = msu, *index[1:] for column in columns: newval = row[column] oldbestval = bestresp[column][bestindex] # the biggest value isn't always the best, e.g. lower meanpkw values are better: if column in ['meanpkw']: if pd.isna(oldbestval) or (newval < oldbestval): bestresp[column][bestindex] = newval # overwrite else: if pd.isna(oldbestval) or (newval > oldbestval): bestresp[column][bestindex] = newval # overwrite ## collect movie model fits by kind, state and opto: #levels = [mvimseustrs, MVIKINDS, ALLST8S] levels = [mvimseustrs, ['nat'], ['none']] names = ['mseu', 'kind', 'st8'] mi = pd.MultiIndex.from_product(levels, names=names) fitcolumns = ['m', 'th', 'rsq', 'nbm', 'nbth', 'nbrsq', 'bm', 'bth', 'brsq'] sampfitcolumns = ['ms', 'ths', 'rsqs', 'nbms', 'nbths', 'nbrsqs', 'bms', 'bths', 'brsqs', 'nrms', 'nrths', 'nrrsqs'] fits = pd.DataFrame(index=mi, columns=fitcolumns) sampfits = pd.DataFrame(index=mi, columns=sampfitcolumns) nsamples = 1000 np.random.seed(0) # fix random seed for identical results from np.random.choice() on each run if LOADFITSFROMPICKLE: print('Loading fits') fits = load('fits') sampfits = load('sampfits') nsamples = len(sampfits[sampfits['ms'].notnull()]['ms'].iloc[0]) # convoluted, but works mseustrfititerator = [] else: print('Calculating fits') mseustrfititerator = mvimseustrs for mseustr in mseustrfititerator: for kind in ['nat']: #MVIKINDS: for st8 in ['none']: #ALLST8S: print(mseustr, kind, st8) psth = mviresp.loc[mseustr, kind, st8]['psth'] if psth.isna().any(): # test both ctrl and opto continue # init sampfits lists: sampfits.loc[mseustr, kind, st8] = [], [], [], [], [], [], [], [], [], [], [], [] # iterate over (narrow) firing modes (all, non-burst, burst, non-rand): for fmode in fmodes: # first do single fit of all relevant trials, fit and test data are # identical (no trial resampling): if fmode != 'nr': # don't bother filling nr column in the fits array psth = mviresp.loc[mseustr, kind, st8][fmode+'psth'] ctrlpsth, optopsth = psth[False], psth[True] mm, b, rsq = fitmodel(ctrlpsth, optopsth, ctrlpsth, optopsth, model=model) th = -b / mm # threshold, i.e., x intercept fits.loc[mseustr, kind, st8][fmode+'m'] = mm fits.loc[mseustr, kind, st8][fmode+'th'] = th fits.loc[mseustr, kind, st8][fmode+'rsq'] = rsq # now sample trials multiple times, separate into fit and test data: if fmode != 'nr': ctrlraster = mviresp.loc[mseustr, kind, st8, False][fmode+'raster'] optoraster = mviresp.loc[mseustr, kind, st8, True][fmode+'raster'] else: # fmode == 'nr' # find number of burst spikes per trial and condition, # remove that same number randomly from each trial in each condition: ctrlraster, optoraster = [], [] ctrlbraster = mviresp.loc[mseustr, kind, st8, False]['braster'] optobraster = mviresp.loc[mseustr, kind, st8, True]['braster'] # find number of burst spikes per trial in both conditions: bctrlcounts = [ len(ctrlbtrial) for ctrlbtrial in ctrlbraster ] boptocounts = [ len(optobtrial) for optobtrial in optobraster ] # build nr rasters from full rasters, trial by trial: ctrlfullraster = mviresp.loc[mseustr, kind, st8, False]['raster'] optofullraster = mviresp.loc[mseustr, kind, st8, True]['raster'] for bctrlcount, ctrltrial in zip(bctrlcounts, ctrlfullraster): nspikes = len(ctrltrial) - bctrlcount # num spikes to sample resampctrltrial = choice(ctrltrial, size=nspikes, replace=False) ctrlraster.append(resampctrltrial) for boptocount, optotrial in zip(boptocounts, optofullraster): nspikes = len(optotrial) - boptocount # num spikes to sample resampoptotrial = choice(optotrial, size=nspikes, replace=False) optoraster.append(resampoptotrial) ctrlraster, optoraster = np.array(ctrlraster), np.array(optoraster) nctrltrials = len(ctrlraster) # typically 200 trials when st8 == 'none' noptotrials = len(optoraster) # typically 200 trials when st8 == 'none' for ntrials in [nctrltrials, noptotrials]: if ntrials < MINNTRIALSAMPLETHRESH: print('Not enough trials to sample:', mseustr, kind, st8, fmode) continue # don't bother sampling, leave entry in sampfits empty ctrltrialis = np.arange(nctrltrials) optotrialis = np.arange(noptotrials) ctrlsamplesize = intround(nctrltrials / 2) # half for fitting, half for testing optosamplesize = intround(noptotrials / 2) # half for fitting, half for testing # probably doesn't matter if use ctrl or opto bins, should be the same: bins = mviresp.loc[mseustr, kind, st8, False]['bins'] for samplei in range(nsamples): # randomly sample half the ctrl and opto trials, without replacement: ctrlfitis = np.sort(choice(ctrltrialis, size=ctrlsamplesize, replace=False)) optofitis = np.sort(choice(optotrialis, size=optosamplesize, replace=False)) ctrltestis = np.setdiff1d(ctrltrialis, ctrlfitis) # get the complement optotestis = np.setdiff1d(optotrialis, optofitis) # get the complement ctrlfitraster, ctrltestraster = ctrlraster[ctrlfitis], ctrlraster[ctrltestis] optofitraster, optotestraster = optoraster[optofitis], optoraster[optotestis] # calculate fit and test opto PSTHs, subsampled in time: ctrlfitpsth = raster2psth(ctrlfitraster, bins, binw, tres, kernel)[::ssx] optofitpsth = raster2psth(optofitraster, bins, binw, tres, kernel)[::ssx] ctrltestpsth = raster2psth(ctrltestraster, bins, binw, tres, kernel)[::ssx] optotestpsth = raster2psth(optotestraster, bins, binw, tres, kernel)[::ssx] #if np.isnan(ctrlfitpsth).any() or np.isnan(optofitpsth).any(): # continue #assert not (np.isnan(ctrltestpsth).any() or np.isnan(optotestpsth).any()) mm, b, rsq = fitmodel(ctrlfitpsth, optofitpsth, ctrltestpsth, optotestpsth, model=model) #if np.isnan(rsq) or rsq < RSQTHRESH: # poor fit # continue # skip this sample th = -b / mm # threshold, i.e., x intercept sampfits.loc[mseustr, kind, st8][fmode+'ms'].append(mm) sampfits.loc[mseustr, kind, st8][fmode+'ths'].append(th) sampfits.loc[mseustr, kind, st8][fmode+'rsqs'].append(rsq) ## FMI (feedback modulation index: (feedback - suppression) / (feedback + suppression)) and ## RMI (run modulation index: (run - sit) / (run + sit)) for each mseu, ## for nat mvi and grt experiments, for all measures. # set up mviFMI DataFrame: levels = [mvimseustrs, ALLST8S] names = ['mseu', 'st8'] mi = pd.MultiIndex.from_product(levels, names=names) mviFMI = pd.DataFrame(index=mi, columns=modmeasures) # set up grtFMI DataFrame: levels = [grtmseustrs, ALLST8S] names = ['mseu', 'st8'] mi = pd.MultiIndex.from_product(levels, names=names) grtFMI = pd.DataFrame(index=mi, columns=modmeasures) # lots of unused columns # set up mviRMI DataFrame: levels = [mvimseustrs, OPTOS] names = ['mseu', 'opto'] mi = pd.MultiIndex.from_product(levels, names=names) mviRMI = pd.DataFrame(index=mi, columns=modmeasures) # set up grtRMI DataFrame: levels = [grtmseustrs, OPTOS] names = ['mseu', 'opto'] mi = pd.MultiIndex.from_product(levels, names=names) grtRMI = pd.DataFrame(index=mi, columns=modmeasures) stimtype2FMI = {'mvi': mviFMI, 'grt': grtFMI} stimtype2RMI = {'mvi': mviRMI, 'grt': grtRMI} for stimtype in STIMTYPES: if stimtype == 'mvi': resp = mviresp.xs('nat', level='kind') # keep only nat movie measures else: resp = grtresp fmidf = stimtype2FMI[stimtype] rmidf = stimtype2RMI[stimtype] for st8 in ALLST8S: for measure in modmeasures: if measure not in resp: continue measurevals = resp[measure][(slice(None), st8)] # mseustr, opto feedback = measurevals[(slice(None), False)] # mseustr suppress = measurevals[(slice(None), True)] # mseustr assert (feedback.dropna() >= 0).all() # negative values violate -1 <= FMI <= 1 assert (suppress.dropna() >= 0).all() # filter out cases where measureval in both conditions is 0: keepis = (feedback != 0) | (suppress != 0) if stimtype == 'mvi': # also filter by SNRTHRESH: maxsnrs = resp['snr'][(slice(None), st8)].max(level='mseu') # max SNR per mseu keepis = keepis & (maxsnrs >= SNRTHRESH) feedback, suppress = feedback[keepis], suppress[keepis] fmis = (feedback - suppress) / (feedback + suppress) # indexed by mseustr for mseustr, fmi in zip(fmis.index.values, fmis): fmidf.loc[mseustr, st8][measure] = fmi # fill in the appropriate rows for opto in OPTOS: for measure in modmeasures: if measure not in resp: continue measurevals = resp[measure][(slice(None), slice(None), opto)] # mseustr, st8 runvals = measurevals[(slice(None), 'run')] # mseustr sitvals = measurevals[(slice(None), 'sit')] # mseustr assert (runvals.dropna() >= 0).all() # negative values violate -1 <= RMI <= 1 assert (sitvals.dropna() >= 0).all() # filter out cases where measureval in both conditions is 0: keepis = (runvals != 0) | (sitvals != 0) if stimtype == 'mvi': # also filter by SNRTHRESH: maxsnrs = resp['snr'][(slice(None), slice(None), opto)].max(level='mseu') keepis = keepis & (maxsnrs >= SNRTHRESH) runvals, sitvals = runvals[keepis], sitvals[keepis] rmis = (runvals - sitvals) / (runvals + sitvals) # indexed by mseustr for mseustr, rmi in zip(rmis.index.values, rmis): rmidf.loc[mseustr, opto][measure] = rmi # fill in the appropriate rows ## maxFMI and maxRMI for each unit (msu), for nat mvi and grt experiments, for all measures. ## For each stimtype+st8+measure+unit combination, take the abs max across experiments. Also ## store the argmax (i.e., the mseu). Among other things, this is for comparison ## of movie and grating responses. # set up maxFMI DataFrame: levels = [mvigrtmsustrs, ALLST8S, STIMTYPES] names = ['msu', 'st8', 'stimtype'] mi = pd.MultiIndex.from_product(levels, names=names) columns = ['mseu'] + modmeasures maxFMI = pd.DataFrame(index=mi, columns=columns) # set up maxRMI DataFrame: levels = [mvigrtmsustrs, OPTOS, STIMTYPES] names = ['msu', 'opto', 'stimtype'] mi = pd.MultiIndex.from_product(levels, names=names) columns = ['mseu'] + modmeasures maxRMI = pd.DataFrame(index=mi, columns=columns) for stimtype in STIMTYPES: fmidf = stimtype2FMI[stimtype] rmidf = stimtype2RMI[stimtype] for st8 in ALLST8S: for measure in modmeasures: mseustrs = fmidf.index.levels[0].values for mseustr in mseustrs: fmi = fmidf.loc[mseustr, st8][measure] if pd.isna(fmi): continue msustr = mseustr2msustr(mseustr) oldfmi = maxFMI.loc[msustr, st8, stimtype][measure] # check existing val if pd.isna(oldfmi) or (abs(fmi) > abs(oldfmi)): maxFMI.loc[msustr, st8, stimtype][measure] = fmi # overwrite maxFMI.loc[msustr, st8, stimtype]['mseu'] = mseustr # overwrite for opto in OPTOS: for measure in modmeasures: mseustrs = rmidf.index.levels[0].values for mseustr in mseustrs: rmi = rmidf.loc[mseustr, opto][measure] if pd.isna(rmi): continue msustr = mseustr2msustr(mseustr) oldrmi = maxRMI.loc[msustr, opto, stimtype][measure] # check existing val if pd.isna(oldrmi) or (abs(rmi) > abs(oldrmi)): maxRMI.loc[msustr, opto, stimtype][measure] = rmi # overwrite maxRMI.loc[msustr, opto, stimtype]['mseu'] = mseustr # overwrite ## calculate eye position stdev and pupil area during all 3 stimtypes, ## as a function of run state and opto conditions, store in ipos_st8 and ipos_opto DataFrames: # Requires huxley file system to be mounted. Author: Davide Crombie # parameters for cleaning up the pupil displacement data: min0t = 1 # only interpolate over stretches missing data shorter than min0t (s) min1t = 1 # only keep stretches of continuous data longer than min1t (s) min1prop = .8 # only keep stretches of data with at least 80% non-NaN st8_rows = [] # to be filled with dictionaries containing relevant data opto_rows = [] print('Calculating ipos_st8 and ipos_opto DataFrames...') for stimtype, stims in zip(['mvi', 'grt', 'spon'], [mvis, grts, spons]): print('...for stimtype %r' % stimtype) stimeye = EyeIxtract() & stims for key in stimeye.get_keys(): # these are mse keys if stimtype == 'mvi': moviecond = mov & key kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping # stim_id values considered in this analysis: kindstimis = tuple(kind2stimis['nat']) if len(kindstimis) == 0: continue # wrong kind of movie, e.g. PNK or SHF msestr = "%s_s%02d_e%02d" % (key['m'], key['s'], key['e']) print(msestr) # fetch and clean the eye position data: irow = EyeIxtract() & key ts = irow.fetch1('i_frame_times') eyedt = np.diff(ts).mean() # eye sampling period # x & y pupil positions relative to baseline (median) in deg. visual angle, # after subtraction of corneal reflection position: xpos, ypos = irow.fetch_position(cr_sub=True, deg=True) # in deg area = irow.fetch_area(mm=False, normalize=True) # in pix # position timecourses after elimintation of unclean stretches, smoothing, # and interpolation of remaining NaN values within clean stretches: min0len, min1len = min0t/eyedt, min1t/eyedt # convert from s to n timepoints xpos_clean = i_clean(xpos, min0len=min0len, min1len=min1len, min1prop=min1prop) ypos_clean = i_clean(ypos, min0len=min0len, min1len=min1len, min1prop=min1prop) area_clean = i_clean(area, min0len=min0len, min1len=min1len, min1prop=min1prop) conditions = {} # get keys for trial table restriction based on st8: conditions['run'] = (st8t & key & {'st8_crit':st82crit['run']}).get_keys() conditions['sit'] = (st8t & key & {'st8_crit':st82crit['sit']}).get_keys() if len(conditions['run']) == 0: # perhaps run data wasn't acquired for this exp assert len(conditions['sit']) == 0 # sanity check del conditions['run'] # del empty list del conditions['sit'] # del empty list # get keys for trial table restriction based on opto condition: if stimtype == 'mvi': if len(kindstimis) > 1: # only include opto conditions correspond to nat movie conditions[False] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=0' & 'stim_id IN {}'.format(kindstimis)).get_keys() conditions[True] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=1' & 'stim_id IN {}'.format(kindstimis)).get_keys() else: conditions[False] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=0').get_keys() conditions[True] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=1').get_keys() # compute eye position variance for all conditions: for condition in conditions: # iterate over all existing keys in conditions row = {} row['mse'] = msestr row['stimtype'] = stimtype # start and stop times for trials of the current state only: ## TODO: spons don't have trials, just periods of opto and non-opto trialis, t0s, t1s = (trial & conditions[condition]).fetch('trial_id', 'trial_on_time', 'trial_off_time') if len(t0s) == 0 or len(t1s) == 0: continue # sort pupil displacement data into trial matrices: xpos_trialmat, _ = get_trialmat(xpos_clean, ts, t0s, t1s) ypos_trialmat, _ = get_trialmat(ypos_clean, ts, t0s, t1s) area_trialmat, area_trialts = get_trialmat(area_clean, ts, t0s, t1s) area_trialmean = np.nanmean(area_trialmat, axis=1) # mean pupil area of each trial # save trial IDs: row['trialis'] = trialis # save mean pupil area for each trial: row['area_trialmean'] = area_trialmean # save pupil area signal for each trial: row['area_trialmat'] = area_trialmat row['area_trialts'] = area_trialts # save mean pupil area across trials: row['mean_area'] = np.nanmean(area_trialmean) # save eye position variability across trials: row['std_xpos_cross'] = np.nanstd(xpos_trialmat, axis=0).mean() row['std_ypos_cross'] = np.nanstd(ypos_trialmat, axis=0).mean() # save eye position variability within trials: row['std_xpos_within'] = np.nanmean(np.nanstd(xpos_trialmat, axis=1)) row['std_ypos_within'] = np.nanmean(np.nanstd(ypos_trialmat, axis=1)) # add dictionary to appropriate list: if condition in ['run', 'sit']: row['st8'] = condition st8_rows.append(row) elif condition in [False, True]: row['opto'] = condition opto_rows.append(row) # convert list of dicts to DataFrame: ipos_st8 = pd.DataFrame(st8_rows) if len(st8_rows) > 0: ipos_st8.set_index(['mse', 'stimtype', 'st8'], inplace=True) # convert to MultiIndex ipos_opto = pd.DataFrame(opto_rows) if len(opto_rows) > 0: ipos_opto.set_index(['mse', 'stimtype', 'opto'], inplace=True) ## collect unit spatial position, store in upos DataFrame: rows = [ {'msu':msu2msustr(msu), 'x':x, 'y':y} for (msu, x, y) in zip(*mviunits.fetch(dj.key, 'u_xpos', 'u_ypos')) ] upos = pd.DataFrame(rows) upos = upos.set_index(['msu']) ## collect runspeed traces, store in runspeed DataFrame: runspeedrows = [] # to be filled with dictionaries containing relevant data for mvi in mvis.fetch(dj.key): msestr = "%s_s%02d_e%02d" % (mvi['m'], mvi['s'], mvi['e']) print(msestr) # fetch run data for all trials in this mse, each is an array, sometimes some trials # are missing from the run data: runtrialis, runt, runspeed = (rt & mvi).fetch('trial_id', 'run_t', 'run_speed') if len(runspeed) == 0: continue # no run data for this mse moviecond = mov & mvi kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping #trials = trial & mvi #alltrialis = np.sort(trials.fetch('trial_id')) # trialis corresponding to rows in tranges for kind in ['nat']: #MVIKINDS: kindstimis = kind2stimis[kind] if len(kindstimis) == 0: continue # wrong kind of movie # get stimis from evcond table for this mvi and kind, and restrict to opto1 chan: mviopto = evcond & mvi & 'stim_id IN %s' % seq2sql(kindstimis) & {'ev_chan':'opto1'} mvioptostimis, mvioptovals = mviopto.fetch('stim_id', 'ec_val') assert len(mvioptostimis) == 2 assert len(mvioptovals) == 2 mvioptovals = np.bool_(mvioptovals) # convert from [0, 1] to [False, True] opto2stimi = {False:int(mvioptostimis[mvioptovals == False]), True :int(mvioptostimis[mvioptovals == True])} for opto in OPTOS: row = {} row['mse'] = msestr row['kind'] = kind row['opto'] = opto # extract run trials for the stimi corresponding to this kind and opto: stimi = opto2stimi[opto] stimtrialis = np.sort((trial & mvi & {'stim_id':stimi}).fetch('trial_id')) # make sure that stimtrialis is a subset of runtrialis: assert issubset(stimtrialis, runtrialis) # now it's safe to ask where stimtrialis are found in runtrialis: trialiis = runtrialis.searchsorted(stimtrialis) #print(opto, len(stimtrialis), len(runtrialis)) row['trialis'] = runtrialis[trialiis] row['t'] = runt[trialiis] row['speed'] = runspeed[trialiis] runspeedrows.append(row) runspeed = pd.DataFrame(runspeedrows) if len(rows) > 0: runspeed = runspeed.set_index(['mse', 'kind', 'opto']) # convert to MultiIndex ## classify units in various ways, store classification in celltype DataFrame: # use normal Index, not MultiIndex, since there's only a single index column, otherwise # assigning values with `df.loc[indexval][colname] = foo` won't work index = pd.Index(mvigrtmsustrs, name='msu') columns = ['sbczscore', 'sbc', 'dsi', 'depth', 'normdepth', 'onoff', 'trans', 'mvionsetpsth', 't'] celltype = pd.DataFrame(index=index, columns=columns) print('Classifying SbC') # NOTE: seems that the zscores stored in the Condition.ZScores table are only # for control conditions, i.e. opto == False, which is good SBCZSCORETHRESH = -3 # 0.0, -1.96, -2.58, -3 # classify by SbC, keep the most negative sbczscore of each msu across grating experiments: msustr2maxzsmseu = {} # map msustr to most negative sbczscore mseu key for grt in grts: zscores = zs & grt uids = (up & grt).fetch('u') for uid in uids: msu = {'m':grt['m'], 's':grt['s'], 'u':uid} mseu = {'m':grt['m'], 's':grt['s'], 'e':grt['e'], 'u':uid} msustr = msufmt.format(**msu) mseustr = mseufmt.format(**mseu) unitzscore = zscores & {'u':uid} if len(unitzscore) == 0: print('No Condition.Zscores entry for', mseustr) continue sbckeys = unitzscore.SbC() assert len(sbckeys) == 1 sbczscore = sbckeys[0]['sbc_zscore'] # median, across ori conditions if np.isnan(sbczscore): print('WARNING: SbC() returned NaN for', mseustr) meanrate = grtresp.loc[mseustr, 'none', False]['meanrate'] if np.isnan(meanrate): continue # skip this non visually responsive mseu, based on z_crit and n_crit oldsbczscore = celltype.loc[msustr]['sbczscore'] if np.isnan(oldsbczscore) or sbczscore < oldsbczscore: # keep most -ve sbczscore celltype.loc[msustr]['sbczscore'] = sbczscore # update msustr2maxzsmseu[msustr] = mseu # classify units as SbC == [False, True] according to SBCZSCORETHRESH: for msustr in mvigrtmsustrs: sbczscore = celltype.loc[msustr]['sbczscore'] if np.isnan(sbczscore): sbc = np.nan # this is the default anyway elif sbczscore <= SBCZSCORETHRESH: sbc = True else: sbc = False celltype.loc[msustr]['sbc'] = sbc # plot rasters, for testing SBCZSCORETHRESH, and zscores in general: try: mseu = msustr2maxzsmseu[msustr] except KeyError: # not a single exp with an sbczscore found for this msu print('Did not collect an sbczscore for', msustr) continue ''' # plot the non-blank, opto == False trials of this msu: offsets = [-1, 1] rasters, tranges, allstimis, opto, ukeys = (spk & mseu).get_rasters(offsets=offsets) assert len(rasters) == 1 uid, raster = rasters.popitem() stimis = allstimis[opto == False] raster = raster[opto == False] blankstimis = (grat & mseu & 'stim_id IN %s' % seq2sql(stimis) & 'grat_contrast = 0.0').fetch('stim_id') assert len(blankstimis) == 1 blankstimi = blankstimis[0] # should be 24 raster = raster[stimis != blankstimi] # exclude blank trials title = "%s, zscore=%.3f, SbC=%s" % (msustr, sbczscore, sbc) ax = simpletraster(raster, dt=5, offsets=offsets, scatter=True, title=title) ax.set_title(title) ''' # classify by direction selectivity, keep biggest DSI of each msu across grating experiments: print('Classifying DSI') resp = grtresp.xs('none', level='st8') # keep only none run state resp = resp.xs(False, level='opto') # keep only control opto state for mseustr, row in resp.iterrows(): print(mseustr) tunparams = row['tun'] if np.isnan(tunparams).all(): continue # skip this mseustr dp, rp, rn, r0, sigma, spon = tunparams dsi = niell_DSI(rp, rn, r0) msustr = mseustr2msustr(mseustr) celltyperow = celltype.loc[msustr] olddsi = celltyperow['dsi'] if np.isnan(olddsi) or dsi > olddsi: celltyperow['dsi'] = dsi # update # estimate dLGN depth of each unit using sparse noise MUA RFs to help delineate # at what channel dLGN starts and ends: print('Estimating dLGN depth') # first plot MUA RFs for each sparse noise experiment in each series, visually inspect them, # and manually enter estimated start and end channel of dLGN for each series into # msstr2dLGNchanrange: for series in mviseries: snkeys = (Series.Experiment() & series & 'e_name LIKE "%%Noise%%"').fetch(dj.key) for snkey in snkeys: print('Sparse noise exp:', snkey) rfs, axes, chans = rfs_for_cdata(snkey) # plot MUA RFs for inspection simple_plot_rfs(snkey, rfs, axes, chans, interpolation=None, nrows=4) # Now save all plots to disk in MUA_RFs folder: saveall(format='pdf') saveall(format='png') # Next, manually map MOUSE_SERIES string to tuple of uppermost and lowermost clear dLGN # channel (end inclusive) based on qualitative estimation of start and end of MUA RF # progression of each series, potentially considering multiple sparse noise experiments # per series: ## TODO: switch to using LGNChans table and its ms2maxchanrange() method instead of ## duplicating it here: msstr2dLGNchanrange = {'PVCre_2017_0006_s03': (29, 8), 'PVCre_2017_0008_s09': (25, 7), 'PVCre_2017_0008_s12': (24, 2), 'PVCre_2017_0015_s03': (28, 5), 'PVCre_2017_0015_s07': (19, 1), 'PVCre_2018_0001_s02': (11, 1), 'PVCre_2018_0001_s05': (25, 6), 'PVCre_2018_0003_s02': (20, 3), 'PVCre_2018_0003_s03': (14, 2), 'PVCre_2018_0003_s05': (20, 9), 'PVCre_2019_0002_s08': (21, 1), 'Ntsr1Cre_2020_0001_s04': (31, 5), } # get dLGN depth of each unit in each series, add to celltype df: for series in mviseries: msstr = ms2msstr(series) hichan, lochan = msstr2dLGNchanrange[msstr] # hi and lo chans (in terms of depth) # could use coords entered into ProbeModel table, but easier to use SiteLoc dict # in probes.expio: #probemodel, chans, coords = (pm & (pr & series)).fetch1('pm', 'pm_chans', 'pm_coords') probemodel = (pr & series).fetch1('pm') # string probe = expio.probes.getprobe(probemodel) # y coord of upper and lowermost dLGN chan, use upper as reference and upper-lower as depth # normalization for all units in this series: y0 = probe.SiteLoc[hichan][1] y1 = probe.SiteLoc[lochan][1] uids, ys = (up & series).fetch('u', 'u_ypos') depths = ys - y0 dLGNspan = y1 - y0 normdepths = depths / dLGNspan for uid, depth, normdepth in zip(uids, depths, normdepths): msustr = msstr + '_u%02d' % uid if depth < 0: print('WARNING: %s has depth %.1f < 0, skipping entry in celltype df' % (msustr, depth)) continue if normdepth > 1.5: print('WARNING: %s has normdepth %.1f > 1.5, skipping entry in celltype df' % (msustr, normdepth)) continue celltype.loc[msustr, 'depth'] = depth celltype.loc[msustr, 'normdepth'] = normdepth ''' # classify by layer (shell/core), pretty much deprecated, don't trust the code in rf.py # that does all of this... # NOTE: several (13/80 as of 2020-02-11) msu end up with no layer classification because # they fall outside of the detected range of chans with MUA envelope RFs of sufficient score print('Classifying shell/core') print('NOTE: something is wrong with this code, hangs on first series at full CPU') for series in mviseries: try: unitlayers = rf.get_layer(series, snexpid=-1) # list of dicts except RuntimeError as err: print(err) continue # no .envl.dat file? for k in unitlayers: msu = {'m':k['m'], 's':k['s'], 'u':k['u']} msustr = msufmt.format(**msu) celltype.loc[msustr]['layer'] = k['layer'] ''' # classify as ON/OFF/transient according to movie onset response. Movies are generally # darker than mid grey, so a decrease in PSTH after onset is an ON response. # After examining lots of raster plots, simplest and most informative may be # +/- 0.2 sec on either side of movie onset, after taking into account 0.05 s propagation # delay, i.e. from -0.15 to 0 as the spont (pre) rate, and from 0.05 to 0.25 as the evoked # (post) rate. Taking just a short period of time around the movie onset helps reduce the # impact of movie motion on the result print('Classifying ON/OFF/transient') resp = mviresp.xs('nat', level='kind') # keep only nat movie kind resp = resp.xs('none', level='st8') # keep only none run state resp = resp.xs(False, level='opto') # keep only control opto state for mseustr, row in resp.iterrows(): print(mseustr) wpsth = row['wpsth'] wt = row['wt'] if np.isnan(wt).any(): continue # skip this mseustr prei0, prei1 = wt.searchsorted([PDT-WINDT, EXCLTR[0]]) posti0, posti1 = wt.searchsorted([EXCLTR[1], PDT+WINDT]) mvionsetpsth = wpsth[prei0:posti1] t = wt[prei0:posti1] prepsth = wpsth[prei0:prei1] postpsth = wpsth[posti0:posti1] insidepsth = wpsth[prei1:posti0] outsidepsth = np.concatenate([prepsth, postpsth]) pre = prepsth.mean() post = postpsth.mean() ins = insidepsth.mean() out = outsidepsth.mean() if pre == 0.0 or post == 0.0 or ins == 0.0 or out == 0.0: print('No pre/post/inside/outside firing rate, skipping') continue onoff = (pre - post)/(pre + post) trans = (ins - out)/(ins + out) Z = max(prepsth.max(), postpsth.max()) # normalization factor, ignore transient period msustr = mseustr2msustr(mseustr) celltyperow = celltype.loc[msustr] oldonoff = celltyperow['onoff'] oldtrans = celltyperow['trans'] if np.isnan(oldonoff) or (abs(onoff) > abs(oldonoff)): celltyperow['onoff'] = onoff # update celltyperow['mvionsetpsth'] = mvionsetpsth / Z # normalize for plotting, update celltyperow['t'] = t # update if np.isnan(oldtrans) or (abs(trans) > abs(oldtrans)): celltyperow['trans'] = trans # update ## code for avoiding annoying dataframe pointer vs copy issues: ''' # overwriting/resetting values of an existing column: celltype['sbc'][:] = np.nan # I think the explicit slice makes the difference # code for adding a new column to an existing dataframe: celltype = celltype.copy() celltype['dsi'] = np.nan celltype['dsi'] = celltype['dsi'].astype(object) celltype = celltype.copy() celltype['onoff'] = np.nan celltype['trans'] = np.nan celltype['mvionsetpsth'] = np.nan celltype['t'] = np.nan celltype['onoff'] = celltype['onoff'].astype(object) celltype['trans'] = celltype['trans'].astype(object) celltype['mvionsetpsth'] = celltype['mvionsetpsth'].astype(object) celltype['t'] = celltype['t'].astype(object) ''' ## collect MUA envelope RF positions for each msu, to use in calculating how well centered ## the screen was on the RF of each unit: index = pd.Index(mvigrtmsustrs, name='msu') columns = ['x0', 'y0'] cellscreenpos = pd.DataFrame(index=index, columns=columns) print('Getting RF params from MUA envelope') snexpid = -1 # use last one in each series MUAENVLRSQTHRESH = 0.2 # consider only decent MUA envelope RF fits series2rffitdf = {} # collect all RF fit params for all series for series in mviseries: erows = Series.Experiment() & series & 'e_name LIKE "%%Noise%%"' expvals = erows.fetch(dj.key, 'e_name', 'e_displayazim', 'e_displayelev') snkeys, enames, azims, elevs = expvals snkey, ename, azim, elev = (snkeys[snexpid], enames[snexpid], azims[snexpid], elevs[snexpid]) print('Sparse noise exp:', snkey, ename) rfs, axes, chans = rfs_for_cdata(snkey) #simple_plot_rfs(snkey, rfs, axes, chans, interpolation=None, nrows=4) fitdf = fit_MUA_RFs(snkey, rfs, chans) msestr = key2msestr(snkey) msstr = mseustr2mstr(msestr) series2rffitdf[msstr] = fitdf # save chans = fitdf.index.levels[0].to_numpy() # all the chans in the df chan2screenpos = {} for chan in chans: chanrow = fitdf.loc[chan, 'mix'] if chanrow['rsq'] < MUAENVLRSQTHRESH: continue # remove azim and elev correction to get RF position wrt screen center: chan2screenpos[chan] = chanrow['x0']-azim, chanrow['y0']-elev urows = up & series uids, maxchans = urows.fetch('u', 'u_maxchan') for uid, maxchan in zip(uids, maxchans): if maxchan in chan2screenpos: screenpos = chan2screenpos[maxchan] msu = {'m':series['m'], 's':series['s'], 'u':uid} msustr = msufmt.format(**msu) cellscreenpos.loc[msustr] = screenpos ## TODO: somehow convert series2rffitdf dict to DataFrame multiindexed by msstr, chan, rftype: #muarfs = pd.DataFrame() ''' # test reliability measure using sparse trials with randomly placed spikes: ntrials, nt, nspikes = 200, 1000, 50 trials = np.zeros((ntrials, nt)) for trial in trials: for spikei in range(nspikes): i = np.random.randint(0, nt) trial[i] = 1 print('reliability:', reliability(trials)[0]) ''' ''' # how to get pairs of ON and OFF values from mviresp instead of resp: # multiindex row slicing is confusing, need the axis=0 for some reason, # see https://pandas.pydata.org/pandas-docs/stable/advanced.html#using-slicers li = mviresp.loc(axis=0) sparsons = np.float64(li[:, st8, True]['spars'].values) sparsoffs = np.float64(li[:, st8, False]['spars'].values) nanis = np.isnan(sparsons) if not (np.isnan(sparsoffs) == nanis).all(): raise ValueError("ON and OFF conditions have different nan indices") sparsons = sparsons[~nanis] sparsoffs = sparsoffs[~nanis] ''' ## TODO: now that mseu are the row indicies instead of msu, they're all unique, and ## the group kwarg isn't doing any grouping, I think ## TODO: I think the MLM stats don't work because currently the various different MVIKINDS ## and ST8S are disabled, only nat and none are being used. So comment this out for now: ''' print('ALL MVI MixedLM') mvirespr = mviresp.reset_index() # make a copy, convert mi to columns mvirespr = mvirespr.dropna() # drop rows with nans, but not really necessary mvirespr = mvirespr[mvirespr['st8'] != 'none'] # exclude 'none' state for col in ['meanrate', 'spars', 'rel', 'meanburstratio']: # each of these columns are object arrs for some reason, instead of float, which # confuses statsmodels and prevents printout of summary. See: # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation mvirespr[col] = mvirespr[col].astype(np.float64) # convert each column to float formulas = ('meanrate ~ kind * opto * st8', 'spars ~ kind * opto * st8', 'rel ~ kind * opto * st8', 'meanburstratio ~ kind * opto * st8') for formula in formulas: # group kwarg specifies how data should be treated, prevents measures from different units # from being lumped together: mlm = smf.mixedlm(formula=formula, data=mvirespr, groups=mvirespr['mseu']) results = mlm.fit() print(formula) print('ngroups:', results.mlm.n_groups) paramp = pd.concat([results.params, results.pvalues], axis=1, keys=['param', 'p']) # Dataframe from 2 Series print(paramp) #print(results.summary()) print() print('----------------------------------------') print() print('NAT only MixedLM') mvirespr = mviresp.reset_index() # make a copy, convert mi to columns natrespr = mvirespr[mvirespr['kind'] == 'nat'] natrespr = natrespr.dropna() # drop rows with nans, but not really necessary natrespr = natrespr[natrespr['st8'] != 'none'] # exclude 'none' state for col in ['meanrate', 'spars', 'rel', 'meanburstratio']: # each of these columns are object arrs for some reason, instead of float, which # confuses statsmodels and prevents printout of summary. See: # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation natrespr[col] = natrespr[col].astype(np.float64) # convert each column to float formulas = ('meanrate ~ opto * st8', 'spars ~ opto * st8', 'rel ~ opto * st8', 'meanburstratio ~ opto * st8') for formula in formulas: # group kwarg specifies how data should be treated, prevents measures from different units # from being lumped together: mlm = smf.mixedlm(formula=formula, data=natrespr, groups=natrespr['mseu']) results = mlm.fit() print(formula) print('ngroups:', results.mlm.n_groups) paramp = pd.concat([results.params, results.pvalues], axis=1, keys=['param', 'p']) # Dataframe from 2 Series print(paramp) #print(results.summary()) print() print('----------------------------------------') print() print('PNK only MixedLM') mvirespr = mviresp.reset_index() # make a copy, convert mi to columns pnkrespr = mvirespr[mvirespr['kind'] == 'pnk'] pnkrespr = pnkrespr.dropna() # drop rows with nans, but not really necessary pnkrespr = pnkrespr[pnkrespr['st8'] != 'none'] # exclude 'none' state for col in ['meanrate', 'spars', 'rel', 'meanburstratio']: # each of these columns are object arrs for some reason, instead of float, which # confuses statsmodels and prevents printout of summary. See: # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation pnkrespr[col] = pnkrespr[col].astype(np.float64) # convert each column to float formulas = ('meanrate ~ opto * st8', 'spars ~ opto * st8', 'rel ~ opto * st8', 'meanburstratio ~ opto * st8') for formula in formulas: # group kwarg specifies how data should be treated, prevents measures from different units # from being lumped together: mlm = smf.mixedlm(formula=formula, data=pnkrespr, groups=pnkrespr['mseu']) results = mlm.fit() print(formula) print('ngroups:', results.mlm.n_groups) paramp = pd.concat([results.params, results.pvalues], axis=1, keys=['param', 'p']) # Dataframe from 2 Series print(paramp) #print(results.summary()) print() print('----------------------------------------') print() print('SHF only MixedLM') mvirespr = mviresp.reset_index() # make a copy, convert mi to columns shfrespr = mvirespr[mvirespr['kind'] == 'shf'] shfrespr = shfrespr.dropna() # drop rows with nans, but not really necessary shfrespr = shfrespr[shfrespr['st8'] != 'none'] # exclude 'none' state for col in ['meanrate', 'spars', 'rel', 'meanburstratio']: # each of these columns are object arrs for some reason, instead of float, which # confuses statsmodels and prevents printout of summary. See: # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation shfrespr[col] = shfrespr[col].astype(np.float64) # convert each column to float formulas = ('meanrate ~ opto * st8', 'spars ~ opto * st8', 'rel ~ opto * st8', 'meanburstratio ~ opto * st8') for formula in formulas: # group kwarg specifies how data should be treated, prevents measures from different units # from being lumped together: mlm = smf.mixedlm(formula=formula, data=shfrespr, groups=shfrespr['mseu']) results = mlm.fit() print(formula) print('ngroups:', results.mlm.n_groups) paramp = pd.concat([results.params, results.pvalues], axis=1, keys=['param', 'p']) # Dataframe from 2 Series print(paramp) #print(results.summary()) print() print('----------------------------------------') print() print('Grating MixedLM') grtrespr = grtresp.reset_index() # make a copy, convert mi to columns grtrespr = grtrespr.dropna() # drop rows with nans, but not really necessary grtrespr = grtrespr[grtrespr['st8'] != 'none'] # exclude 'none' state for col in ['meanrate', 'meanburstratio']: # each of these columns are object arrs for some reason, instead of float, which # confuses statsmodels and prevents printout of summary. See: # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation grtrespr[col] = grtrespr[col].astype(np.float64) # convert each column to float formulas = ('meanrate ~ opto * st8', 'meanburstratio ~ opto * st8') for formula in formulas: # group kwarg specifies how data should be treated, prevents measures from different units # from being lumped together: mlm = smf.mixedlm(formula=formula, data=grtrespr, groups=grtrespr['mseu']) results = mlm.fit() print(formula) print('ngroups:', results.mlm.n_groups) paramp = pd.concat([results.params, results.pvalues], axis=1, keys=['param', 'p']) # Dataframe from 2 Series print(paramp) #print(results.summary()) print() print('----------------------------------------') print() '''