main.py 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. """Collect rasters of responses for all LGN opto experiments as a function of stimulus type
  2. (movies, gratings), run state, and opto state. Then, calculate various measures and save
  3. results to DataFrames. Some variables (DataFrames and mseu string lists) can be optionally
  4. saved to and loaded from .pickle files to avoid long and unnecessary recalculation times.
  5. Run this script at the IPython/Jupyter command line; the `-i` runs scripts in IPython's
  6. namespace instead of an empty one:
  7. `run -i main.py`
  8. `pvmvis` (default), `ntsrmvis` or `negntsrmvis` can be provided as an `--exptype=` argument
  9. to this script, to load and later plot PV-Cre (most of the paper), Ntsr1-Cre (Fig1S4),
  10. or Ntsr negative (Fig1S5) data, respectively.
  11. """
  12. import os
  13. import sys
  14. import argparse
  15. from numpy import pi
  16. from numpy.random import choice
  17. import scipy.stats
  18. from scipy.stats import (linregress, ttest_rel, ttest_ind, ttest_1samp, wilcoxon,
  19. ks_2samp)
  20. from scipy.stats.mstats import gmean
  21. from scipy.optimize import least_squares
  22. import matplotlib.pyplot as plt
  23. from matplotlib.offsetbox import AnchoredText
  24. from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
  25. import seaborn as sns
  26. import pandas as pd
  27. from util import (load, load_all, save, export2csv,
  28. desat, axes_disable_scientific,
  29. ms2msstr, msu2msustr, mseustr2msustr, mseustrs2msustrs, mseustr2mstr,
  30. mseustrs2mstrs, mseustr2msestr, mseustrs2msestrs,
  31. findmse, fitmodel, residual_rsquared, linear_loss, get_max_snr,
  32. intround, split_tranges, wrap_raster,
  33. cf, saveall, wintitle, simpletraster,
  34. raster2psth, raster2freqcomp, sparseness, reliability, snr_baden2016,
  35. get_psth_peaks_gac, vector_OSI, percentile_ci, sum_of_gaussians)
  36. ## Define constants:
  37. EXPTYPE = 'pvmvis' # pvmvis, ntsrmvis, negntsrmvis
  38. if __name__ == '__main__':
  39. # parse command-line arguments:
  40. parser = argparse.ArgumentParser(description='Main script for natfeedback paper')
  41. parser.add_argument("--exptype", help="pvmvis, ntsrmvis, negntsrmvis")
  42. args = parser.parse_args()
  43. EXPTYPE = args.exptype or EXPTYPE
  44. print('*** EXPTYPE: %r' % EXPTYPE)
  45. RATETHRESH = 0.01 # mean firing rate threshold for unit inclusion, across all trial types, Hz
  46. SNRTHRESH = 0.015
  47. SCATTERPTHRESH = 0.05
  48. SGNF2ALPHA = {True: 1, False: 0.4}
  49. # visual drive unit inclusion criteria for gratings, ignores st8 and opto:
  50. VDNCRIT = 1 # min num points in tuning curve that are at least VDZCRIT sems away from spont rate
  51. VDZCRIT = 2.5
  52. PSTHBASELINEMEDIANX = 2 # PSTH median multiplier to estimate baseline
  53. PSTHPEAKMINTHRESH = 3 # PSTH peak detection threshold relative to estimated baseline, Hz
  54. STIMTYPES = 'mvi', 'grt'
  55. stimtype2axislabel = {'grt':'grating', 'mvi':'movie'}
  56. MVIKINDS = 'nat',
  57. ALLSTIMTYPES = MVIKINDS + ('grt',)
  58. STIMTYPESPLUSALL = list(STIMTYPES) + ['mvi+grt'] # strings
  59. STIMTYPESPLUSALLI = list(STIMTYPES) + [slice(None)] # for use as indices
  60. # iterate over run states in this order:
  61. ST8S = ['run', 'sit'] # excludes 'none'
  62. ALLST8S = ['none', 'run', 'sit']
  63. ST8COMBOS = [['none'], ['run', 'sit'], ['run', 'sit']]
  64. OPTOCOMBOS = [[False, True], [False], [True]]
  65. st82crit = {'run':'(run_speed > 1).mean() > 0.5',
  66. 'sit':'(run_speed < 0.25).mean() > 0.5',
  67. 'none':'none'}
  68. rungreen = '#' + ''.join([ format(val, '02x') for val in [66, 168, 117] ])
  69. sitorange = '#' + ''.join([ format(val, '02x') for val in [235, 115, 9] ])
  70. st82clr = {'run':rungreen, 'sit':sitorange, 'none':'black'}
  71. st82alpha = {None:1, 'run':1, 'sit':0.5} # sitting is like suppression, so desaturate colour
  72. st82clrmap = {'run':plt.cm.Greens, 'sit':plt.cm.Oranges, 'none':plt.cm.Greys}
  73. # define example mseus and their colors:
  74. green = '#007f00' # deep green
  75. violet = '#9f3fff' # pure violet (7f00ff) is a little too dark
  76. darkblue = '#000064' # control condition
  77. optoblue = '#2a7fff' # light shade of blue for opto condition
  78. black = '#000000'
  79. asterisk = 5, 2, 0 # 5 sided, asterisk style, 0 deg rotation, has rounded corners though
  80. DEFSZ = 50 # default scatter plot point size (open points)
  81. exmpli2clr = {1:violet, 2:green, 3:optoblue} # example neurons 1, 2, 3
  82. # use different marker shapes to help distinguish example points:
  83. exmpli2mrk = {1:'x', 2:'+', 3:asterisk}
  84. exmpli2sz = {1:60, 2:70, 3:60} # marker size
  85. exmpli2lw = {1:2, 2:2, 3:2} # (marker) line width
  86. fig1exmpli = 1 # fig1 uses example neuron 1
  87. fig5exmpli = 1 # fig5 uses example neuron 1
  88. # designate movie example units:
  89. pvmvimseu2exmpli = {'PVCre_2018_0003_s03_e03_u51':1, # example neuron 1
  90. 'PVCre_2017_0008_s12_e06_u56':2, # example neuron 2
  91. }
  92. ntsrmvimseu2exmpli = {'Ntsr1Cre_2019_0008_s03_e04_u18':1,
  93. }
  94. exptype2mviexmplis = {'pvmvis':pvmvimseu2exmpli, 'ntsrmvis':ntsrmvimseu2exmpli}
  95. mvimseu2exmpli = exptype2mviexmplis[EXPTYPE]
  96. # designate grating example units:
  97. pvgrtmseu2exmpli = {'PVCre_2018_0003_s03_e02_u51':1, # example neuron 1
  98. 'PVCre_2018_0003_s05_e02_u164':3, # example neuron 3
  99. }
  100. ntsrgrtmseu2exmpli = {'Ntsr1Cre_2019_0008_s06_e02_u46':1,
  101. }
  102. exptype2grtexmplis = {'pvmvis':pvgrtmseu2exmpli, 'ntsrmvis':ntsrgrtmseu2exmpli}
  103. grtmseu2exmpli = exptype2grtexmplis[EXPTYPE]
  104. mvimsu2exmpli = { mseustr2msustr(k):v for k, v in mvimseu2exmpli.items() }
  105. grtmsu2exmpli = { mseustr2msustr(k):v for k, v in grtmseu2exmpli.items() }
  106. msu2exmpli = {**mvimsu2exmpli, **grtmsu2exmpli} # merge the two
  107. burstclr = 'red'
  108. OPTOS = False, True
  109. opto2tuni = {False:0, True:1} # for indexing into tun_params field in Tuning table
  110. opto2fb = {None:'', False:'feedback', True:'suppression'}
  111. opto2clr = {None:darkblue, False:darkblue, True:optoblue}
  112. opto2alpha = {None:1, False:1, True:0.5} # with opto, feedback is removed, so desaturate colour
  113. modmeasures = ['meanrate', 'meanrate02', 'meanrate35', 'blankmeanrate', 'blankcondmeanrate',
  114. 'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio',
  115. 'spars', 'rel', 'meanpkw', 'snr']
  116. modmeasuresnoblankcond = modmeasures.copy()
  117. modmeasuresnoblankcond.remove('blankcondmeanrate')
  118. modmeasuresnoblankcond.remove('blankcondmeanburstratio')
  119. modmeasuresnoblank = modmeasuresnoblankcond.copy()
  120. modmeasuresnoblank.remove('blankmeanrate')
  121. modmeasuresnoblank.remove('blankmeanburstratio')
  122. measure2axislabel = {'meanrate':'FR', 'meanrate02':'FR', 'meanrate35':'FR',
  123. 'blankmeanrate':'FR', 'blankcondmeanrate':'FR',
  124. 'meanburstratio':'burst ratio', 'blankmeanburstratio':'burst ratio',
  125. 'blankcondmeanburstratio':'burst ratio',
  126. 'spars':'sparseness', 'rel':'reliability',
  127. 'meanpkw':'mean peak width', 'snr':'SNR'}
  128. measure2axisunits = {'meanrate':' (spk/s)', 'meanrate02':' (spk/s)', 'meanrate35':' (spk/s)',
  129. 'meanpkw':' (s)'}
  130. short2longaxislabel = {'FR':'Firing rate'}
  131. fmodes = ['', 'nb', 'b', 'nr'] # (narrow) firing modes
  132. fmode2txt = {'':'all', 'nb':'nonburst', 'b':'burst', 'nr':'nonrand'}
  133. # Tuning table keys used for gratings:
  134. ivs_order = 'grat_orientation, opto1'
  135. tun_model = 'sum_of_gaussians'
  136. model = 'threshlin' # linear or threshlin
  137. relaverage = 'mean' # reliability averaging method: mean or median
  138. DEFAULTFIGURESIZE = 2.8, 2.8 # normal
  139. DUMBBELLFIGSIZE = 3.3, 4.4 # dumbbell plots
  140. RASTERSCALEX = 0.75
  141. OFFSETS = -0.5, 0.75 # wide raster time offsets, sec
  142. OFFSETS02 = 0, -3 # limits to spikes from t=0 to 2 s, assuming 5 sec movies
  143. OFFSETS35 = 3, 0 # limits to spikes from t=3 to 5 s
  144. NTRIALSMVIGRT = 120 # limit movies to first 120 trials for more direct comparison with gratings
  145. PSTHHEIGHT = 3
  146. binw, tres, kernel = 0.02, 0.001, 'gauss' # for calculating PSTHs, sec
  147. psthplottres = 0.010 # sec
  148. ssx = intround(psthplottres / tres) # subsample factor to get the desired plot tres
  149. MINNTRIALSAMPLETHRESH = 10
  150. # ON/OFF/transient classification:
  151. PDT = 0.050 # propagation delay time for movie onset (s)
  152. WINDT = 0.20 # duration separating start of pre and PDT, and PDT and end of post (s)
  153. EXCLTR = 0, 2*PDT # time range around PDT to treat as transient response
  154. assert EXCLTR[0] <= PDT <= EXCLTR[1]
  155. # load data from pickles:
  156. subfolder = EXPTYPE if EXPTYPE != 'pvmvis' else ''
  157. name2val = load_all(subfolder=subfolder)
  158. locals().update(name2val)