"""Figure 1 and 1S2 plots, and first and last 2 sec scatter plots for 4S1, use run -i fig1.py""" mi = pd.MultiIndex.from_product([mvimseustrs, OPTOS], names=['mseu', 'opto']) # include single trial and trial-averaged columns: columns = ['trialis', 'rates', 'rate02s', 'rate35s', 'burstratios', 'blankrates', 'blankburstratios'] + modmeasuresnoblankcond fig1 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8 # plot movie rasters for each unit by kind, state and opto condition: showbursts = True SCATTER = True # True: low quality, but fast; False: high quality, but slow trialiis = None #np.arange(120) # plot only these trials, None plots all trials for mseustr in []:#mvimseu2exmpli: #mvimseustrs: for kind in ['nat']: #MVIKINDS for st8 in ['none']: for opto in OPTOS: dt = mviresp.loc[mseustr, kind, st8, opto]['dt'] # stimulus duration (s) wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1] optotrange = mviresp.loc[mseustr, kind, st8, opto]['optotrange'] raster = mviresp.loc[mseustr, kind, st8, opto]['wraster'] if np.any(pd.isna(raster)): # no raster for this condition continue title = '%s movie %s %s %s raster' % (mseustr, kind, st8, opto) if trialiis is not None: raster = raster[trialiis] title += ' %dtrials' % len(trialiis) burstis = None if showbursts: wburstis = mviresp.loc[mseustr, kind, st8, opto]['wburstis'] a = simpletraster(raster=raster, alpha=0.5, s=1, dt=dt, offsets=OFFSETS, scatter=SCATTER, title=title, inchespersec=1.5*RASTERSCALEX, inchespertrial=1/25*RASTERSCALEX, xaxis=True, burstis=wburstis) if opto: # plot horizontal line signifying opto ON period, just above axes: ton, toff = optotrange a.hlines(-5, ton, toff, colors=optoblue, lw=4, clip_on=False, in_layout=False) # plot horizontal line signifying stimulus period, just above axes: a.hlines(-2, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False) # plot vertical coloured line just right of axes: #vlinexpos = dt + OFFSETS[1] + wdt*0.01 #y0, y1 = a.get_ylim() #clr = st82clr[st8] #alpha = opto2alpha[opto] #a.vlines(vlinexpos, y0, y1, colors=clr, lw=4, alpha=alpha, clip_on=False, # in_layout=False) a.set_xlim(tmin, tmax) # plot movie PSTHs for each unit by kind, overplot state and opto combinations: for mseustr in []: #mvimseu2exmpli:#['Ntsr1Cre_2020_0001_s04_e10_u12']: #mvimseustrs for kind in MVIKINDS: wt = mviresp.loc[mseustr, kind]['wt'] # DataFrame wpsth = mviresp.loc[mseustr, kind]['wpsth'] # DataFrame wbpsth = mviresp.loc[mseustr, kind]['wbpsth'] # DataFrame #t = mviresp.loc[mseustr, kind]['t'] # DataFrame psth = mviresp.loc[mseustr, kind]['psth'] # DataFrame if np.isnan(np.hstack(wpsth)).any(): # no PSTH for at least one condition for this mseu continue optorho, optorho_p = scipy.stats.pearsonr(psth['none', False], psth['none', True]) runsitrho, runsitrho_p = scipy.stats.pearsonr(psth['run', False], psth['sit', False]) skewoff = scipy.stats.skew(psth['none', False]) skewon = scipy.stats.skew(psth['none', True]) skewrun = scipy.stats.skew(psth['run', False]) skewsit = scipy.stats.skew(psth['sit', False]) _, optoKS_p = ks_2samp(psth['none', False], psth['none', True]) _, runsitKS_p = ks_2samp(psth['run', False], psth['sit', False]) print('%s example PSTH:\n' 'optorho=%g optorho_p=%g skewoff=%g skewon=%g optoKS_p=%g\n' 'runsitrho=%g runsitrho_p=%g skewrun=%g skewsit=%g runsitKS_p=%g' % (mseustr, optorho, optorho_p, skewoff, skewon, optoKS_p, runsitrho, runsitrho_p, skewrun, skewsit, runsitKS_p)) for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS): st8combostr = ' '.join(st8combo) optocombostr = ' '.join([str(opto) for opto in optocombo]) dt = mviresp.loc[mseustr, kind, 'none']['dt'].mean() # mean stimulus duration (s) wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1] figsize = 0.872 + wdt*1.5*RASTERSCALEX, PSTHHEIGHT f, a = plt.subplots(figsize=figsize) title = '%s %s %s %s PSTH' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) # subsample as in PSTH scatter plots, and plot all 4 st8/opto combinations: for st8 in st8combo: color = st82clr[st8] for opto in optocombo: if st8combo == ['run', 'sit'] and optocombo == [True]: optoalpha = 1 else: optoalpha = opto2alpha[opto] if optocombo == [False]: # desaturate bursts by st8, not opto burstalpha = st82alpha[st8] else: burstalpha = opto2alpha[opto] c = desat(color, optoalpha) # do manual alpha mixing bc = desat(burstclr, burstalpha) # do manual alpha mixing fb = opto2fb[opto].title() if st8 == 'none': label = fb else: label = ', '.join([fb, st8.title()]) a.plot(wt[st8, opto][::ssx], wpsth[st8, opto][::ssx], '-', color=c, label=label) a.plot(wt[st8, opto][::ssx], wbpsth[st8, opto][::ssx], '-', color=bc, label=label+', Burst') l = a.legend(frameon=False) l.set_draggable(True) a.set_xlabel('Time (s)') a.set_ylabel('Firing rate (spk/s)') a.set_xlim(tmin, tmax) # plot horizontal line signifying stimulus period, just below data: ymin, ymax = a.get_ylim() a.hlines(-ymax*0.025, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False) a.set_ylim(ymin, ymax) # restore y limits from before plotting the hline # plot PDFs, CDFs and QQ of PSTHs for each unit by kind, overplot state and opto combinations: percentiles = np.arange(0, 100, 2) figsize = DEFAULTFIGURESIZE for mseustr in mvimseu2exmpli: #mvimseustrs if EXPTYPE != 'pvmvis': break # skip these plots for kind in MVIKINDS: psths = mviresp.loc[mseustr, kind]['psth'] if np.isnan(np.hstack(psths)).any(): # no PSTH for at least one condition for this mseu continue for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS): st8combostr = ' '.join(st8combo) optocombostr = ' '.join([str(opto) for opto in optocombo]) pdff, pdfa = plt.subplots(figsize=figsize) # PDF title = '%s %s %s %s PSTH PDF' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) hpdff, hpdfa = plt.subplots(figsize=(figsize[0], 3)) # horizontal, match PSTH height title = '%s %s %s %s horizontal PSTH PDF' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) npdff, npdfa = plt.subplots(figsize=figsize) # normalized PDF title = '%s %s %s %s PSTH nPDF' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) cdff, cdfa = plt.subplots(figsize=figsize) title = '%s %s %s %s PSTH CDF' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) qqf, qqa = plt.subplots(figsize=figsize) title = '%s %s %s %s PSTH QQ' % (mseustr, kind, st8combostr, optocombostr) wintitle(title) # plot PDFs, CDFs and QQ of all st8/opto combinations: for st8 in st8combo: c = st82clr[st8] qq = {} for opto in optocombo: psth = psths[st8, opto][::ssx] npsth = psth / psth.max() # normalize uvals = np.unique(npsth) qq[opto] = np.percentile(psth, percentiles) fb = opto2fb[opto].title() if st8 == 'none': label = fb else: label = ', '.join([fb, st8.title()]) alpha = opto2alpha[opto] nbins = intround(np.sqrt(len(psth))) pdfa.hist(psth, bins=nbins, density=True, histtype='step', color=c, alpha=alpha, label=label) hpdfa.hist(psth, bins=nbins, density=True, histtype='step', color=c, alpha=alpha, label=label, orientation='horizontal') npdfa.hist(npsth, bins=nbins, density=True, histtype='step', color=c, alpha=alpha, label=label) cdfa.hist(npsth, bins=uvals, density=True, cumulative=True, histtype='step', color=c, alpha=alpha, label=label) if len(optocombo) == 2: qqa.scatter(qq[False], qq[True], marker='.', c=c) #l = a.legend(frameon=False) #l.set_draggable(True) pdfa.set_xlabel('Firing rate (spk/s)') pdfa.set_ylabel('Probability density') #pdfa.set_yticks([0, 0.5, 1]) pdff.show() hpdfa.set_xlabel('Probability density') hpdfa.set_ylabel('Firing rate (spk/s)') hpdff.show() npdfa.set_xlabel('Normalized rate') npdfa.set_ylabel('Probability density') #pdfa.set_yticks([0, 0.5, 1]) npdff.show() cdfa.set_xlabel('Normalized rate') cdfa.set_ylabel('Cumulative probability') cdfa.set_yticks([0, 0.5, 1]) cdff.show() if len(optocombo) == 2: qqmax = np.ceil(np.array(qqa.dataLim).max()) qqxyline = [0, qqmax], [0, qqmax] qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1) qqa.set_xlabel('Feedback rate (spk/s)') qqa.set_ylabel('Suppression rate (spk/s)') qqa.set_xlim([0, qqmax]) qqa.set_ylim([0, qqmax]) qqa.set_aspect('equal') #qqa.set_xticks([0, 0.5, 1]) #qqa.set_yticks([0, 0.5, 1]) qqf.show() # plot one PSTH PDF, CDF and QQ, for each combination of kind, state and opto, # each pooled across all mseu: percentiles = np.arange(0, 100, 2) figsize = DEFAULTFIGURESIZE for kind in MVIKINDS: if EXPTYPE != 'pvmvis': break # skip these plots for st8combo in ST8COMBOS: pdff, pdfa = plt.subplots(figsize=figsize) title = 'mean %s %s PSTH PDF' % (kind, ' '.join(st8combo)) wintitle(title) cdff, cdfa = plt.subplots(figsize=figsize) title = 'mean %s %s PSTH CDF' % (kind, ' '.join(st8combo)) wintitle(title) qqf, qqa = plt.subplots(figsize=figsize) title = 'mean %s %s PSTH QQ' % (kind, ' '.join(st8combo)) wintitle(title) for st8 in st8combo: c = st82clr[st8] qq = {} for opto in OPTOS: psths, npsths = [], [] for mseustr in mvimseustrs: #mvimseu2exmpli psth = mviresp.loc[mseustr, kind, st8, opto]['psth'] if np.isnan(psth).any(): # no PSTH for this mseu continue #psth = psth[::ssx] # no need to subsample, only calculating a few CDFs psths.append(psth) npsths.append(psth / psth.max()) # normalize psths, npsths = np.concatenate(psths), np.concatenate(npsths) uvals = np.unique(npsths) qq[opto] = np.percentile(psths, percentiles) fb = opto2fb[opto].title() if st8 == 'none': label = fb else: label = ', '.join([fb, st8.title()]) alpha = opto2alpha[opto] pdfa.hist(npsths, bins=50, density=True, histtype='step', color=c, alpha=alpha, label=label) cdfa.hist(npsths, uvals, density=True, cumulative=True, histtype='step', color=c, alpha=alpha, label=label) qqa.scatter(qq[False], qq[True], marker='.', c=c) #l = a.legend(frameon=False) #l.set_draggable(True) pdfa.set_xlabel('Normalized rate') pdfa.set_ylabel('Probability density') #pdfa.set_yticks([0, 0.5, 1]) pdff.show() cdfa.set_xlabel('Normalized rate') cdfa.set_ylabel('Cumulative probability') cdfa.set_yticks([0, 0.5, 1]) cdff.show() qqmax = np.ceil(np.array(qqa.dataLim).max()) qqxyline = [0, qqmax], [0, qqmax] qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1) qqa.set_xlabel('Feedback rate (spk/s)') qqa.set_ylabel('Suppression rate (spk/s)') #qqa.set_aspect('equal', 'datalim') qqa.set_xlim([0, qqmax]) qqa.set_ylim([0, qqmax]) #qqa.set_xticks([0, 0.5, 1]) #print('%s xticks: %r' % (title, qqa.get_xticks())) #qqa.set_yticks(qqa.get_xticks()) qqf.show() # scatter plot movie meanrates: figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak to make space for log units logmin, logmax = -1, 2 logticks = np.array([-1, 0, 1, 2]) #log0min = logmin + 0.05 for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto meanrate movie %s %s' % (kind, st8)) rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: meanrate = mviresp.loc[mseustr, kind, st8]['meanrate'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if meanrate.isna().any(): # missing one or both meanrates continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue rons.append(meanrate[True]) roffs.append(meanrate[False]) trialis = mviresp.loc[mseustr, kind, st8]['trialis'] rates = mviresp.loc[mseustr, kind, st8]['rates'] _, pval = ttest_ind(rates[False], rates[True], equal_var=False) sgnfs.append(pval < SCATTERPTHRESH) # bool fig1.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save fig1.loc[mseustr]['trialis'] = trialis # save, for joining DataFrames trial-wise fig1.loc[mseustr]['rates'] = rates # save trial-wise values if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment rons = np.asarray(rons) roffs = np.asarray(roffs) sgnfs = np.asarray(sgnfs) # plot y=x line: xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) insgnfis, = np.where(sgnfs == False) sgnfis, = np.where(sgnfs == True) normlinsgnfis = np.intersect1d(normlis, insgnfis) normlsgnfis = np.intersect1d(normlis, sgnfis) if len(normlinsgnfis) > 0: # plot normal (non-example) insignificant points: c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) if len(normlsgnfis) > 0: # plot normal (non-example) significant points: c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(rons[exmpli], roffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xlabel('Suppression FR (spk/s)') a.set_ylabel('Feedback FR (spk/s)') a.set_xscale('log') a.set_yscale('log') a.set_xlim(10**logmin, 10**logmax) a.set_ylim(10**logmin, 10**logmax) a.set_xticks(10.0**logticks) a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks a.minorticks_off() axes_disable_scientific(a) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(rons, roffs) # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False)) #mu = rons.mean(), roffs.mean() #txt = '$\mathregular{\mu=%.1f, %.1f}$' % mu #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False)) # plot mean firing rate PDF: f, a = plt.subplots(figsize=figsize) wintitle('opto meanrate PDF movie %s %s' % (kind, st8)) bins = np.logspace(-1, 2, num=25, base=10) # end inclusive a.hist(roffs, bins=bins, histtype='step', density=False, color='k', label='Feedback', clip_on=False) a.hist(rons, bins=bins, histtype='step', density=False, color=optoblue, label='V1 Suppression', clip_on=False) a.set_xlabel('Mean firing rate (spk/s)') a.set_ylabel('Cell count') a.set_xscale('log') a.set_xlim(10**logmin, 10**logmax) a.set_xticks(10.0**logticks) a.minorticks_off() axes_disable_scientific(a) #l = a.legend(frameon=False) #l.set_draggable(True) # scatter plot movie burst ratio: figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units logmin, logmax = -3.55, 0 if EXPTYPE == 'ntsrmvis': logmin, logmax = -3.83, 0 logticks = np.array([-3, -2, -1, 0]) for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto burst ratio movie %s %s' % (kind, st8)) brons, broffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: br = mviresp.loc[mseustr, kind, st8]['meanburstratio'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if br.isna().any(): # missing for at least one opto condition continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue brons.append(br[True]) broffs.append(br[False]) burstratios = mviresp.loc[mseustr, kind, st8]['burstratios'] _, pval = ttest_ind(burstratios[False], burstratios[True], equal_var=False) sgnfs.append(pval < SCATTERPTHRESH) # bool fig1.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save fig1.loc[mseustr]['burstratios'] = burstratios # save trial-wise values if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment brons, broffs = np.asarray(brons), np.asarray(broffs) sgnfs = np.asarray(sgnfs) # plot y=x line: xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) insgnfis, = np.where(sgnfs == False) sgnfis, = np.where(sgnfs == True) normlinsgnfis = np.intersect1d(normlis, insgnfis) normlsgnfis = np.intersect1d(normlis, sgnfis) if len(normlinsgnfis) > 0: # plot normal (non-example) insignificant points: c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing a.scatter(brons[normlinsgnfis], broffs[normlinsgnfis], clip_on=True, marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails if len(normlsgnfis) > 0: # plot normal (non-example) significant points: c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing a.scatter(brons[normlsgnfis], broffs[normlsgnfis], clip_on=True, marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(brons[exmpli], broffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xlabel('Suppression BR') # keep it short to maximize space for axes a.set_ylabel('Feedback BR') a.set_xscale('log') a.set_yscale('log') a.set_xlim(10**logmin, 10**logmax) a.set_ylim(10**logmin, 10**logmax) a.set_xticks(10.0**logticks) a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks a.minorticks_off() axes_disable_scientific(a) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(brons, broffs) # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False)) # scatter plot movie PSTH sparseness: figsize = DEFAULTFIGURESIZE linmin, linmax, linstep = 0, 1, 0.5 ticks = np.arange(intround(linmin), intround(linmax)+linstep, linstep) for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto sparseness %s %s' % (kind, st8)) sparsons, sparsoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: spars = mviresp.loc[mseustr, kind, st8]['spars'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if spars.isna().any(): # missing for at least one opto condition continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue sparsons.append(spars[True]) sparsoffs.append(spars[False]) fig1.loc[mseustr, 'spars'] = spars[False], spars[True] # save if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment sparsons, sparsoffs = np.asarray(sparsons), np.asarray(sparsoffs) # plot y=x line: xyline = [linmin, linmax], [linmin, linmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) # plot normal (non-example) points: a.scatter(sparsons[normlis], sparsoffs[normlis], clip_on=False, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(sparsons[exmpli], sparsoffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xlabel('Suppression spars') # keep it short to maximize space for axes a.set_ylabel('Feedback spars') a.set_xlim(linmin, linmax) a.set_ylim(linmin, linmax) a.set_xticks(ticks) a.set_yticks(ticks) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(sparsons, sparsoffs) # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False)) # scatter plot movie PSTH reliability: figsize = DEFAULTFIGURESIZE[0]*0.98, DEFAULTFIGURESIZE[1] # tweak to match others linmin, linmax = 0, 0.52 ticks = 0, 0.5 for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto reliability %s %s' % (kind, st8)) relons, reloffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: rel = mviresp.loc[mseustr, kind, st8]['rel'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if rel.isna().any(): # missing for at least one opto condition continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue relons.append(rel[True]) reloffs.append(rel[False]) # collect rhos from all trial pairs, reduce to single rho per trial, use # to estimate significance of individual units: rhos = mviresp.loc[mseustr, kind, st8]['rhos'] fullrhos, meanrhos = {}, {} for opto in OPTOS: npairs = len(rhos[opto]) # from quadratic formula solution to npairs = ntrials*(ntrials-1) / 2 ntrials = (1 + np.sqrt(1 + 8*npairs)) / 2 assert np.isclose(ntrials, int(ntrials)) # make sure this float is int ntrials = int(ntrials) uti = np.triu_indices(ntrials, k=1) # upper indices, excludes diagonal lti = np.tril_indices(ntrials, k=-1) # lower indices, excludes diagonal # empty full corr matrix, nans will remain on diagonal and will be ignored: fullrhos[opto] = np.tile(np.nan, (ntrials, ntrials)) fullrhos[opto][uti] = rhos[opto] # fill the upper triangle # copy upper to lower triangle, # see https://stackoverflow.com/a/42209263/2020363: fullrhos[opto][lti] = fullrhos[opto].T[lti] meanrhos[opto] = np.nanmean(fullrhos[opto], axis=1) # along either axis is fine _, pval = ttest_ind(meanrhos[False], meanrhos[True], equal_var=False) sgnfs.append(pval < SCATTERPTHRESH) # bool fig1.loc[mseustr, 'rel'] = rel[False], rel[True] # save if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment relons, reloffs = np.asarray(relons), np.asarray(reloffs) sgnfs = np.asarray(sgnfs) # plot y=x line: xyline = [linmin, linmax], [linmin, linmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) insgnfis, = np.where(sgnfs == False) sgnfis, = np.where(sgnfs == True) normlinsgnfis = np.intersect1d(normlis, insgnfis) normlsgnfis = np.intersect1d(normlis, sgnfis) if len(normlinsgnfis) > 0: # plot normal (non-example) insignificant points: c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing a.scatter(relons[normlinsgnfis], reloffs[normlinsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) if len(normlsgnfis) > 0: # plot normal (non-example) significant points: c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing a.scatter(relons[normlsgnfis], reloffs[normlsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(relons[exmpli], reloffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xlabel('Suppression rel') # keep it short to maximize space for axes a.set_ylabel('Feedback rel') a.set_xlim(linmin, linmax) a.set_ylim(linmin, linmax) a.set_xticks(ticks) a.set_yticks(ticks) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(relons, reloffs, nan_policy='omit') # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False)) # plot movie SNR PDFs, show SNRTHRESH: figsize = DEFAULTFIGURESIZE bins = np.arange(0, 0.5, 0.005) #bins = np.logspace(-2, -0.2, 50, base=10) for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: for opto in [None, False, True]: if opto is not None: snrs = mviresp['snr'][(slice(None), kind, st8, opto)].dropna() # one val per mseu else: snrs = mviresp['snr'][(slice(None), kind, st8)].dropna() # two vals per mseu? snrs = np.float64(snrs) # convert from object array f, a = plt.subplots(figsize=figsize) titlestr = 'SNR PDF %s %s' % (kind, st8) if opto is not None: titlestr += ' %s' % opto2fb[opto] wintitle(titlestr, f) c = st82clr[st8] a.hist(snrs, bins=bins, density=True, color=c, alpha=opto2alpha[opto], label=opto2fb[opto]) a.axvline(x=SNRTHRESH, ls='-', marker='', c='gray')#, zorder=-np.inf) txt = 'n=%d' % len(snrs) a.add_artist(AnchoredText(txt, loc='upper right', frameon=False, prop=dict(color=c, alpha=opto2alpha[opto]))) #a.set_xscale('log') #l = a.legend(frameon=False) #l.set_draggable(True) a.set_xlabel('SNR') if opto is not None: a.set_xlabel('%s SNR' % opto2fb[opto].title()) a.set_ylabel('Probability density') # scatter plot movie SNR: figsize = DEFAULTFIGURESIZE linmin, linmax = 0, 0.6 ticks = 0, 0.5 for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto SNR %s %s' % (kind, st8)) snrons, snroffs, exmplis, exmplmseustrs, normlis = [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: snr = mviresp.loc[mseustr, kind, st8]['snr'] if snr.isna().any(): # missing for at least one opto condition continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue snrons.append(snr[True]) snroffs.append(snr[False]) fig1.loc[mseustr, 'snr'] = snr[False], snr[True] # save if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment snrons, snroffs = np.asarray(snrons), np.asarray(snroffs) # plot y=x line: xyline = [linmin, linmax], [linmin, linmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) # plot normal (non-example) points: a.scatter(snrons[normlis], snroffs[normlis], clip_on=False, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(snrons[exmpli], snroffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xlabel('Suppression SNR') a.set_ylabel('Feedback SNR') a.set_xlim(linmin, linmax) a.set_ylim(linmin, linmax) a.set_xticks(ticks) a.set_yticks(ticks) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(snrons, snroffs, nan_policy='omit') # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False)) # plot movie PSTH peak width distributions: normalizepeakwidths = True figsize = DEFAULTFIGURESIZE for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) titlestr = ('opto peak width PDF %s %s PSTHBASELINEMEDIANX=%d PSTHPEAKMINTHRESH=%d' % (kind, st8, PSTHBASELINEMEDIANX, PSTHPEAKMINTHRESH)) if normalizepeakwidths: titlestr += " normed" wintitle(titlestr, f) pkwons, pkwoffs = [], [] for mseustr in mvimseustrs: pkws = mviresp.loc[mseustr, kind, st8]['pkws'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if pkws.isna().any(): continue # skip this mseu if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue if normalizepeakwidths: maxpkw = np.concatenate(pkws.values).max() # max peak width for this mseu pkws = pkws / maxpkw # normalize all peak widths by maxpkw pkwons.append(pkws[True]) pkwoffs.append(pkws[False]) c = st82clr[st8] pkwons, pkwoffs = np.concatenate(pkwons), np.concatenate(pkwoffs) for opto, pkws, loc in zip(OPTOS, [pkwoffs, pkwons], ['upper right', 'upper left']): a.hist(pkws, bins=50, density=True, histtype='step', color=c, alpha=opto2alpha[opto], label=opto2fb[opto]) txt = 'n=%d' % len(pkws) a.add_artist(AnchoredText(txt, loc=loc, frameon=False, prop=dict(color=c, alpha=opto2alpha[opto]))) a.set_xscale('log') #a.set_yscale('log') axes_disable_scientific(a) #l = a.legend(frameon=False) #l.set_draggable(True) a.set_xlabel('PSTH peak width (s)') a.set_ylabel('Probability density') #t, p = ttest_rel(pkwoffs, pkwons, nan_policy='omit') # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False)) # scatter plot movie PSTH mean peak widths: figsize = DEFAULTFIGURESIZE[0]*1.05, DEFAULTFIGURESIZE[1] # tweak to make space for log units linmin, linmax = 0.025, 0.2 ticks = 0.05, 0.1, 0.2 for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: f, a = plt.subplots(figsize=figsize) wintitle('opto mean peak width %s %s' % (kind, st8)) meanpkwons, meanpkwoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: meanpkw = mviresp.loc[mseustr, kind, st8]['meanpkw'] snr = mviresp.loc[mseustr, kind, st8]['snr'] if meanpkw.isna().any(): # missing for at least one opto condition continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR continue meanpkwons.append(meanpkw[True]) meanpkwoffs.append(meanpkw[False]) fig1.loc[mseustr, 'meanpkw'] = meanpkw[False], meanpkw[True] # save if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment meanpkwons, meanpkwoffs = np.asarray(meanpkwons), np.asarray(meanpkwoffs) # plot y=x line: xyline = [linmin, linmax], [linmin, linmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) # plot normal (non-example) points: a.scatter(meanpkwons[normlis], meanpkwoffs[normlis], clip_on=False, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(meanpkwons[exmpli], meanpkwoffs[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) a.set_xscale('log', basex=2) a.set_yscale('log', basey=2) a.set_xlabel('Suppression pkw (s)') # keep it short to maximize space for axes a.set_ylabel('Feedback pkw (s)') a.set_xlim(linmin, linmax) a.set_ylim(linmin, linmax) a.set_xticks(ticks) a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks a.minorticks_off() axes_disable_scientific(a) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) a.xaxis.set_major_formatter(ScalarFormatter()) a.yaxis.set_major_formatter(ScalarFormatter()) a.xaxis.set_major_formatter(FormatStrFormatter("%.2g")) a.yaxis.set_major_formatter(FormatStrFormatter("%.2g")) #t, p = ttest_rel(meanpkwoffs, meanpkwons, nan_policy='omit') # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False)) """Figure 1S2 FMI plots""" # scatter plot FMI of all movie measures vs meanrate FMI: figsize = DEFAULTFIGURESIZE linmin, linmax, linstep = -1, 1, 1 ticks = np.arange(linmin, linmax+linstep, linstep) for kind in ['nat']:#MVIKINDS: for st8 in ['none']:#ALLST8S: for measure in modmeasuresnoblank: if measure.startswith('meanrate'): continue # don't bother plotting meanrate* against meanrate axislabel = measure2axislabel[measure] rfmis, msrfmis = [], [] # rate FMIs, measure FMIs exmplis, exmplmseustrs, normlis = [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: rfmi = mviFMI.loc[mseustr, st8]['meanrate'] msrfmi = mviFMI.loc[mseustr, st8][measure] if pd.isna(rfmi) or pd.isna(msrfmi): # missing at least one FMI value continue rfmis.append(rfmi) msrfmis.append(msrfmi) if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment rfmis, msrfmis = np.asarray(rfmis), np.asarray(msrfmis) ## scatter plot FMIs: f, a = plt.subplots(figsize=figsize) wintitle('opto FMI %s vs FMI rate %s %s' % (measure, kind, st8)) # plot x=0 and y=0 lines: a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf) a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf) # plot normal (non-example) points: a.scatter(rfmis[normlis], msrfmis[normlis], clip_on=False, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(rfmis[exmpli], msrfmis[exmpli], clip_on=False, marker=marker, c=c, s=sz, lw=lw) # get fname of appropriate LMM .cvs file: if measure == 'meanburstratio': # use Steffen's LMM linregress fit, for fig1S2 fname = os.path.join('stats', 'figure_1_S2c_coefs.csv') elif measure == 'spars': # use Steffen's LMM linregress fit, for fig1S2 fname = os.path.join('stats', 'figure_1_S2d_coefs.csv') elif measure == 'rel': # use Steffen's LMM linregress fit, for fig1S2 fname = os.path.join('stats', 'figure_1_S2e_coefs.csv') elif measure == 'snr': # use Steffen's LMM linregress fit, for fig1S2 fname = os.path.join('stats', 'figure_1_S2f_coefs.csv') elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig1S2 fname = os.path.join('stats', 'figure_1_S2g_coefs.csv') else: print("WARNING: No LMM stats for %s measure" % measure) fname = None # fetch LMM linregress fit params from .csv: if fname: try: df = pd.read_csv(fname) foundregression = True except FileNotFoundError: print('Missing file: %s' % fname) foundregression = False if foundregression: mm = df['slope'][0] b = df['intercept'][0] x = np.array([rfmis.min(), rfmis.max()]) y = mm * x + b a.plot(x, y, '-', color='red') # plot linregress fit a.set_xlabel('FR FMI') if axislabel.islower(): axislabel = axislabel.capitalize() ylabel = '%s FMI' % axislabel a.set_ylabel(ylabel) a.set_xlim(linmin, linmax) a.set_ylim(linmin, linmax) a.set_xticks(ticks) a.set_yticks(ticks) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) """Some figure 4S1 plots""" # scatter plot meanrates of first 2 and last 2 s of movies, limited to first 120 trials # per condition, as per standard grating experiments: figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units logmin, logmax = -2, 2 logticks = np.array([-2, 0, 2]) #log0min = logmin + 0.05 kind = 'nat' st8 = 'none' for meanratestr, rangestr in zip(['meanrate02', 'meanrate35'], ['0-2', '3-5']): f, a = plt.subplots(figsize=figsize) wintitle('opto meanrate movie t=%s s %d trials %s %s' % (rangestr, NTRIALSMVIGRT, kind, st8)) rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], [] keptmseui = 0 # manually init and increment instead of using enumerate() for mseustr in mvimseustrs: meanrate = mviresp.loc[mseustr, kind, st8][meanratestr] snr = mviresp.loc[mseustr, kind, st8]['snr'] if meanrate.isna().any(): # missing one or both meanrates continue if (snr < SNRTHRESH).all(): # neither condition has decent SNR on the full range continue rons.append(meanrate[True]) roffs.append(meanrate[False]) ratesstr = meanratestr.split('mean')[1]+'s' rates = mviresp.loc[mseustr, kind, st8][ratesstr] _, pval = ttest_ind(rates[False], rates[True], equal_var=False) sgnfs.append(pval < SCATTERPTHRESH) # bool fig1.loc[mseustr, meanratestr] = meanrate[False], meanrate[True] # save fig1.loc[mseustr][ratesstr] = rates # save trial-wise values if mvimseu2exmpli.get(mseustr) == fig1exmpli: exmplis.append(keptmseui) exmplmseustrs.append(mseustr) else: normlis.append(keptmseui) keptmseui += 1 # manually increment rons = np.asarray(rons) roffs = np.asarray(roffs) sgnfs = np.asarray(sgnfs) # plot y=x line: xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax] a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1) insgnfis, = np.where(sgnfs == False) sgnfis, = np.where(sgnfs == True) normlinsgnfis = np.intersect1d(normlis, insgnfis) normlsgnfis = np.intersect1d(normlis, sgnfis) # plot normal (non-example) insignificant points: c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) # plot normal (non-example) significant points: c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False, marker='.', c='None', edgecolor=c, s=DEFSZ) # plot normal (non-example) points: #a.scatter(rons[normlis], roffs[normlis], clip_on=False, # marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ) # plot example points, one at a time: for exmpli, mseustr in zip(exmplis, exmplmseustrs): marker = exmpli2mrk[mvimseu2exmpli[mseustr]] c = exmpli2clr[mvimseu2exmpli[mseustr]] sz = exmpli2sz[mvimseu2exmpli[mseustr]] lw = exmpli2lw[mvimseu2exmpli[mseustr]] a.scatter(rons[exmpli], roffs[exmpli], marker=marker, clip_on=False, c=c, s=sz, lw=lw) a.set_xlabel('Suppression FR (spk/s)') a.set_ylabel('Feedback FR (spk/s)') a.set_xscale('log') a.set_yscale('log') a.set_xlim(10**logmin, 10**logmax) a.set_ylim(10**logmin, 10**logmax) a.set_xticks(10.0**logticks) a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks a.minorticks_off() axes_disable_scientific(a) a.set_aspect('equal') a.spines['left'].set_position(('outward', 4)) a.spines['bottom'].set_position(('outward', 4)) #t, p = ttest_rel(rons, roffs) # paired t-test #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False)) #mu = rons.mean(), roffs.mean() #txt = '$\mathregular{\mu=%.1f, %.1f}$' % mu #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))