OSCR

Projection targeting with phototagging to study the structure and function of retinal ganglion cells.

Code ↔ Paper

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

The 3 matches
  1. [1] § STAR★Methods › Quantification and statistical analysis › ReaChR cell identification ↔ Phototagging_methods_code/phototagging_functions/Reachr_classification_AMR.m, lines 54–101 · score 0.66 · ReaChR signal, PC1, strengths, variance, deviations, histogram
  2. [2] § STAR★Methods › Method details › Visual stimulation ↔ src/yass/rf/run.py, lines 135–246 · score 0.58 · frame rate, white noise, MATLAB, stimuli, temporal
  3. [3] § STAR★Methods › Quantification and statistical analysis › Receptive field estimation ↔ src/yass/rf/run.py, lines 135–246 · score 0.55 · Gaussian fit, pixels, triggered, frames, RF, STA

Paper

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

The paper is loaded when this pane is shown.

The authors' code

Python · 542 lines · 20 KB · Apache-2.0 · 2 matches

  1. import h5py
  2. import parmap
  3. import numpy as np
  4. import scipy.io as sio
  5. import os
  6. from tqdm import tqdm
  7. import scipy
  8. import scipy.io
  9. import scipy.optimize as opt
  10. from scipy.spatial.distance import cdist
  11. from scipy.ndimage import gaussian_filter
  12. import networkx as nx
  13. from pkg_resources import resource_filename
  15. from yass import read_config
  16. from yass.template import upsample_resample, shift_chans
  17. from yass.rf.sta_fit import get_fit_on_sta
  18. from yass.rf.util import get_rf, get_circle_plotting_data, classifiy_contours
  20. def run():
  21. """RF computation
  22. """
  24. CONFIG = read_config()
  26. stim_movie_file = os.path.join(CONFIG.data.root_folder, CONFIG.data.stimulus)
  27. triggers_fname = os.path.join(CONFIG.data.root_folder, CONFIG.data.triggers)
  28. spike_train_fname = os.path.join(CONFIG.path_to_output_directory,
  29. 'spike_train.npy')
  30. saving_dir = os.path.join(CONFIG.path_to_output_directory, 'rf')
  32. rf = RF(stim_movie_file, triggers_fname, spike_train_fname, saving_dir)
  33. rf.calculate_STA()
  34. rf.detect_multi_rf()
  35. rf.classification()
  38. class RF(object):
  39. def __init__(self, saving_dir, stim_movie_file,triggers_fname,
  40. spike_train_fname, soft_assignment_fname=None,
  41. fname_classification_boundary=None, matlab_bin='matlab'):
  43. # default parameter
  44. self.n_color_channels = 3
  45. self.sp_frame_rate = 20000
  46. #self.data_sample_len = 36000000 # len of white noise data (this script doesn't look at natural scenes)
  48. self.load_spike_train(spike_train_fname)
  49. if soft_assignment_fname is not None:
  50. self.soft_assignment = np.load(soft_assignment_fname)
  51. else:
  52. self.soft_assignment = np.ones(self.sps.shape[0])
  54. self.stim_movie_file = stim_movie_file
  55. self.triggers_fname = triggers_fname
  57. self.save_dir = saving_dir
  58. if not os.path.exists(self.save_dir):
  59. os.makedirs(self.save_dir)
  60. if self.save_dir[-1] != '/':
  61. self.save_dir += '/'
  63. self.matlab_bin = matlab_bin
  65. self.fname_classification_boundary = fname_classification_boundary
  67. print("spike train:\t{}".format(self.sps.shape))
  68. print("Number of units:\t{}".format(self.Ncells))
  70. def load_stimulus_trigger(self, stim_movie_file, triggers_fname):
  72. print('Loading Stimulus...')
  74. # Load stim file
  75. h5_temp = h5py.File(stim_movie_file, 'r')
  76. self.WN_stim = h5_temp['movie'][:]
  77. h5_temp.close()
  78. self.stim_size = self.WN_stim.shape[2:4]
  79. self.WN_stim = self.WN_stim.reshape((-1, self.n_color_channels,
  80. self.stim_size[0]*self.stim_size[1]))
  82. ## Load triggers
  83. if triggers_fname.split('.')[-1] == 'trig':
  84. with open(triggers_fname, 'rb'):
  85. self.WN_trigger_times = np.fromfile(triggers_fname, dtype='int16')
  86. elif triggers_fname.split('.')[-1] == 'mat':
  87. self.WN_trigger_times = sio.loadmat(triggers_fname)
  88. self.WN_trigger_times = self.WN_trigger_times['triggers'].flatten().astype('float')
  90. print("stim movie:\t{}".format(self.WN_stim.shape))
  92. np.save(self.save_dir+'stim_size.npy',self.stim_size)
  94. def calculate_frame_times(self):
  96. frame_per_pulse = 100
  98. ## Find pulses and calculate frame times
  99. # Get first locations of pulses in seconds
  100. pulses = np.where(np.diff(self.WN_trigger_times)==-2048)[0]+1 # find where pulse starts (diff+1)
  101. pulses_seconds = pulses / float(self.sp_frame_rate) # divide by 20k Hz to get seconds
  103. self.frame_times = np.interp(
  104. np.arange(0,frame_per_pulse * pulses_seconds.shape[0]),
  105. np.arange(0,frame_per_pulse * pulses_seconds.shape[0], frame_per_pulse),
  106. pulses_seconds)
  109. def load_spike_train(self, spike_train_fname):
  111. print('Loading Spike Train...')
  113. ## Load spikes
  114. sps_file_ext = os.path.splitext(spike_train_fname)[1]
  115. if sps_file_ext == '.mat':
  116. #sps = sio.loadmat(spike_train_fname)['spike_train'].astype('int32')
  118. # for single columnd data
  119. sps_temp = sio.loadmat(spike_train_fname)['spike_train'].astype('int32')
  120. unique_ids = np.unique(sps_temp[:,1])
  121. unique_ids = unique_ids[unique_ids>0] - 1
  122. self.sps = np.zeros(sps_temp.shape, 'int32')
  123. self.sps[:, 0] = sps_temp[:, 0]
  125. for i, k in enumerate(unique_ids):
  126. idx = sps_temp[:, 1] == (k+1)
  127. self.sps[idx, 1] = i
  129. elif sps_file_ext == '.npy':
  130. self.sps = np.load(spike_train_fname)
  132. # Get number of cells/units
  133. self.Ncells = int(np.max(self.sps[:,1])+1)
  135. def calculate_STA(self):
  137. self.load_stimulus_trigger(self.stim_movie_file, self.triggers_fname)
  138. self.calculate_frame_times()
  140. tmp_dir = os.path.join(self.save_dir, 'tmp')
  141. if not os.path.exists(tmp_dir):
  142. os.makedirs(tmp_dir)
  144. tmp_dir_sta = os.path.join(tmp_dir, 'sta')
  145. if not os.path.exists(tmp_dir_sta):
  146. os.makedirs(tmp_dir_sta)
  148. tmp_dir_rgc = os.path.join(tmp_dir, 'rgc')
  149. if not os.path.exists(tmp_dir_rgc):
  150. os.makedirs(tmp_dir_rgc)
  152. print('Calculating STA...')
  154. ############################################
  155. ## Get full STAs and spatial/temporal STA ##
  156. ############################################
  158. STA_temporal_length = 30 # how many bins/frames to include in STA
  159. Ncells = self.Ncells
  160. stim_size = self.stim_size
  161. n_color_channels = self.n_color_channels
  162. n_pixels = stim_size[0]*stim_size[1]
  164. unique_ids = np.unique(self.sps[:,1])
  166. args_in = []
  167. for i_cell in np.arange(Ncells):
  168. fname = os.path.join(tmp_dir_sta, 'unit_'+str(i_cell)+'.mat')
  169. if not os.path.exists(fname):
  171. ##################################
  172. ### Get spikes in stimulus bins ##
  173. ##################################
  175. # Get spike times of this cell in seconds
  176. idx_ = np.where(self.sps[:,1]==i_cell)[0]
  177. these_sps = self.sps[idx_, 0]
  178. #spikes before 36000000 are white noise spikes, divide by frame rate to get seconds
  179. these_sps = these_sps / float(self.sp_frame_rate)
  181. weight = self.soft_assignment[idx_]
  183. ## Line up spikes with frames
  184. binned_spikes = weighted_histogram(these_sps, weight, self.frame_times)
  185. which_spikes = np.where(binned_spikes>0)[0]
  186. which_spikes = which_spikes[which_spikes>STA_temporal_length]
  188. args_in.append([
  189. self.WN_stim,
  190. binned_spikes,
  191. which_spikes,
  192. STA_temporal_length,
  193. stim_size,
  194. fname
  195. ])
  197. if False:
  198. n_processors = 6
  199. parmap.map(sta_calculation_parallel,
  200. args_in,
  201. processes=n_processors,
  202. pm_pbar=True)
  203. else:
  204. for unit in tqdm(range(len(args_in))):
  205. sta_calculation_parallel(args_in[unit])
  207. sta_array = np.zeros((Ncells, stim_size[0], stim_size[1],
  208. n_color_channels, STA_temporal_length))
  209. for unit in range(Ncells):
  210. fname = os.path.join(tmp_dir_sta, 'unit_'+str(unit)+'.mat')
  211. sta = sio.loadmat(fname)['temp_stas']
  212. sta_array[unit] = sta
  214. ## run matlab code
  215. #print('running matlab code')
  216. #rf_matlab_loc = resource_filename('yass', 'rf/rf_matlab')
  217. #command = '{} -nodisplay -r \"cd(\'{}\'); fit_sta_liam_parallel(\'{}\', \'{}\'); exit\"'.format(
  218. # self.matlab_bin, rf_matlab_loc, tmp_dir_sta, tmp_dir_rgc)
  219. #print(command)
  220. #os.system(command)
  221. #print('done running matlab code')
  223. #STA_spatial = np.zeros((self.Ncells, stim_size[0], stim_size[1], n_color_channels))
  224. #STA_temporal = np.zeros((self.Ncells, STA_temporal_length, n_color_channels))
  225. #gaussian_fits = np.zeros((self.Ncells, 5))
  226. #for unit in unique_ids:
  227. # fname = os.path.join(tmp_dir_rgc, 'rgc_{}.mat'.format(unit))
  228. # try:
  229. # data = scipy.io.loadmat(fname)
  230. # if 'temp_rf' in data.keys():
  231. # STA_spatial[unit] = data['temp_rf']
  232. # STA_temporal[unit] = data['fit_tc']
  233. # gaussian_fits[unit] = data['temp_fit_params']['fit_params'][0][0][0][:5]
  234. # except:
  235. # print('unit {} corrupted'.format(unit))
  237. STA_spatial, STA_temporal, gaussian_fits = get_fit_on_sta(sta_array)
  239. # hack for now
  240. STA_spatial = np.tile(STA_spatial[:, :, :, None],
  241. (1, 1, 1, n_color_channels))
  242. STA_temporal = STA_temporal.transpose(0, 2, 1)
  244. np.save(os.path.join(self.save_dir, 'STA_spatial.npy'), STA_spatial)
  245. np.save(os.path.join(self.save_dir, 'STA_temporal.npy'), STA_temporal)
  246. np.save(os.path.join(self.save_dir, 'gaussian_fits.npy'), gaussian_fits)
  248. def detect_multi_rf(self):
  250. STA_spatial = np.load(self.save_dir+'STA_spatial.npy')
  251. n_units = STA_spatial.shape[0]
  253. n_rfs = np.zeros(n_units)
  254. for j in range(n_units):
  255. rf = STA_spatial[j][:, :, 1]
  256. n_rfs[j] = len(get_rf(rf - np.mean(rf), 2))
  258. # yass
  259. idx = np.where(n_rfs==1)[0]
  260. np.save(self.save_dir+'idx_single_rf.npy', idx)
  262. idx = np.where(n_rfs==0)[0]
  263. np.save(self.save_dir+'idx_no_rf.npy', idx)
  265. idx = np.where(n_rfs > 1)[0]
  266. np.save(self.save_dir+'idx_multi_rf.npy', idx)
  268. def load_data_for_classification(self, load_contours=False):
  270. # load data
  271. sta_spatial = np.load(os.path.join(self.save_dir, 'STA_spatial.npy'))
  272. sta_spatial[np.isnan(sta_spatial)] = 0
  273. sta_temporal = np.load(os.path.join(self.save_dir, 'STA_temporal.npy'))
  274. sta_temporal[np.isnan(sta_temporal)] = 0
  275. gaussian_fits = np.load(os.path.join(self.save_dir, 'gaussian_fits.npy'))
  277. n_units = sta_temporal.shape[0]
  279. spike_train = self.sps
  280. unique_ids, n_spikes = np.unique(spike_train[:,1], return_counts=True)
  281. firing_rates = np.zeros(n_units)
  282. firing_rates[unique_ids] = n_spikes/(np.ptp(spike_train[:,0])/self.sp_frame_rate)
  284. max_loc = np.abs(sta_temporal[:,:,1]).argmax(1)
  285. sign = np.sign(sta_temporal[np.arange(n_units), max_loc])
  287. peak_val = np.zeros((n_units, 3))
  288. for j in range(n_units):
  289. sta_ = (sta_spatial[j].reshape(-1, 3))*sign[j][None]
  290. peak_val[j] = sta_[np.max(sta_, 1).argmax()]
  291. peak_val = peak_val*sign
  292. green_val = peak_val[:, 1]
  294. gaussian_sd = gaussian_fits[:, 3:5]
  296. if load_contours:
  297. contours = np.zeros((n_units, 64, 2))
  298. for j in range(n_units):
  299. xy = get_circle_plotting_data(j, gaussian_fits)
  300. contours[j] = xy.T
  301. else:
  302. contours = None
  304. return gaussian_sd, green_val, firing_rates, contours
  306. def classification(self, fname_classification_boundary=None):
  308. gaussian_sd, green_val, f_rates, _ = self.load_data_for_classification()
  310. idx_single = np.load(os.path.join(self.save_dir, 'idx_single_rf.npy'))
  312. if fname_classification_boundary is None:
  313. fname_classification_boundary = self.fname_classification_boundary
  315. temp = np.load(fname_classification_boundary)
  316. sd_mean_noise_th = temp['sd_mean_noise_th']
  317. sd_ratio_noise_th = temp['sd_ratio_noise_th']
  318. green_noise_th = temp['green_noise_th']
  319. midget_on_th = temp['midget_on_th']
  320. midget_off_th = temp['midget_off_th']
  321. large_on_th = temp['large_on_th']
  322. large_off_th = temp['large_off_th']
  323. sbc_fr_th = temp['sbc_fr_th']
  325. labels_single, cell_types = classifiy_contours(
  326. gaussian_sd[idx_single],
  327. green_val[idx_single],
  328. f_rates[idx_single],
  329. sd_mean_noise_th,
  330. sd_ratio_noise_th,
  331. green_noise_th,
  332. midget_on_th,
  333. midget_off_th,
  334. large_on_th,
  335. large_off_th,
  336. sbc_fr_th)
  338. labels = np.ones(self.Ncells, 'int32')*-1
  339. labels[idx_single] = labels_single
  341. np.save(os.path.join(self.save_dir, 'labels.npy'), labels)
  342. np.save(os.path.join(self.save_dir, 'cell_types.npy'), cell_types)
  344. def twoD_Gaussian(self, xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
  345. ## Define 2D Gaussian that we'll fit to spatial STAs
  346. (x, y) = xdata_tuple
  347. xo = float(xo)
  348. yo = float(yo)
  349. a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
  350. b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
  351. c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
  352. g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)+c*((y-yo)**2)))
  353. return g.ravel()
  355. def fit_gaussian(self):
  357. if not os.path.exists(self.save_dir+'Gaussian_params.npy'):
  359. print('Fitting Gaussian on STA...')
  361. stim_size = self.stim_size
  363. STA_spatial = np.load(self.save_dir+'STA_spatial.npy')
  365. ## Fit Gaussian to STA
  366. use_green_only = True
  367. if use_green_only:
  368. this_STA_spatial = STA_spatial[:,1]
  369. else:
  370. this_STA_spatial = STA_spatial_colorcat
  372. Gaussian_params = np.zeros((self.Ncells,7))
  373. Gaussian_params[:]=np.nan
  374. nonconverged_Gaussian_cells = np.empty(0,) # keep track of cells where fitting procedure doesn't converge
  376. # Loop over cells
  377. for i_cell in tqdm(range(self.Ncells)):
  379. # Get STA for this cell
  380. this_STA = this_STA_spatial[i_cell].reshape((-1,))
  382. # Create x and y indices for grid for Gaussian fit
  383. x = np.arange(0, stim_size[1], 1)
  384. y = np.arange(0, stim_size[0], 1)
  385. x, y = np.meshgrid(x, y)
  387. # Get initial guess for Gaussian parameters (helps with fitting)
  388. init_amp = this_STA[np.argmax(np.abs(this_STA))] # get amplitude guess from most extreme (max or min) amplitude of this_STA
  389. init_x,init_y = np.unravel_index(np.argmax(np.abs(this_STA)),(stim_size[0],stim_size[1])) # guess center of Gaussian as indices of most extreme (max or min) amplitude
  390. initial_guess = (init_amp,init_y,init_x,2,2,0,0)
  392. # Try to fit, if it doesn't converge, log that cell
  393. try:
  394. popt, pcov = opt.curve_fit(self.twoD_Gaussian, (x, y), this_STA, p0=initial_guess)
  395. Gaussian_params[i_cell] = popt
  396. Gaussian_params[i_cell,3:5] = np.abs(popt[3:5]) # sometimes sds are negative (in Gaussian def above, they're always squared)
  397. except:
  398. nonconverged_Gaussian_cells = np.append(nonconverged_Gaussian_cells,i_cell)
  400. np.save(self.save_dir+'Gaussian_params.npy',Gaussian_params)
  403. def get_denoiser(STA):
  404. max_timecourse = np.max(np.abs(STA[:,20:]), axis=1)
  405. max_timecourse_mean = np.mean(max_timecourse, axis=2)
  406. max_timecourse_std = np.std(max_timecourse, axis=2)
  407. good_ones = max_timecourse > (max_timecourse_mean + 4*max_timecourse_std)[:,:,np.newaxis]
  409. denoiser = np.zeros((3, 3, STA.shape[1]))
  410. for color in range(3):
  411. unit_id, pixel_id = np.where(good_ones[:,color])
  412. good_timecourse = STA[unit_id, :, color, pixel_id]
  413. [U,S,V] = np.linalg.svd(good_timecourse.T)
  414. denoiser[color] = U[:,:3].T
  416. return denoiser
  418. def denoise_STA(STA):
  420. denoiser = get_denoiser(STA)
  422. n_units, _, _, n_pixels = STA.shape
  423. n_colors, n_filter, n_time = denoiser.shape
  425. STA_denoised = np.zeros(STA.shape)
  426. for color in range(n_colors):
  427. deno = np.matmul(denoiser[color].T, denoiser[color])
  428. STA_temp = STA[:,:,color].transpose(0,2,1).reshape(-1, n_time)
  429. STA_denoised[:,:,color] = np.matmul(STA_temp, deno).reshape(n_units, n_pixels, n_time).transpose(0,2,1)
  431. return STA_denoised
  433. def get_temp_spat_filters(STA):
  434. [U,S,V] = np.linalg.svd(STA)
  436. temp_filter = U[:, 0]
  437. spatial_filter = V[0]
  438. temp_sign = np.sign(temp_filter[np.abs(temp_filter).argmax()])
  439. spat_sign = np.sign(spatial_filter[np.abs(spatial_filter).argmax()])
  440. sign = temp_sign*spat_sign
  441. if temp_sign != sign:
  442. temp_filter *= -1.0
  443. if spat_sign != sign:
  444. spatial_filter *= -1.0
  446. return temp_filter, spatial_filter
  449. def sta_calculation_parallel(arg_in):
  451. WN_stim = arg_in[0]
  452. binned_spikes = arg_in[1]
  453. which_spikes = arg_in[2]
  454. STA_temporal_length = arg_in[3]
  455. stim_size = arg_in[4]
  456. fname = arg_in[5]
  458. ####################
  459. ### Calculate STA ##
  460. ####################
  462. ## Swap out fastest version here
  463. _, n_color_channels, n_pixels = WN_stim.shape
  464. STA = np.zeros((STA_temporal_length, n_color_channels, n_pixels))
  465. for i in range(which_spikes.shape[0]):
  466. bin_number = which_spikes[i]
  467. STA += binned_spikes[bin_number]*WN_stim[bin_number-(STA_temporal_length-1):bin_number+1]
  469. if which_spikes.shape[0] == 0:
  470. STA += 0.5
  472. # full sta
  473. if np.sum(binned_spikes[STA_temporal_length:])>0:
  474. STA = STA/np.sum(binned_spikes[STA_temporal_length:])
  475. STA = STA.reshape(STA_temporal_length, n_color_channels,
  476. stim_size[0], stim_size[1])
  477. STA = STA.transpose(2,3,1,0)
  479. scipy.io.savemat(fname, mdict={'temp_stas': STA})
  481. def align_tc(tc, ref):
  482. n_units, n_timepoints, n_channels = tc.shape
  484. max_channels = np.abs(tc).max(1).argmax(1)
  486. main_tc = np.zeros((n_units, n_timepoints))
  487. for j in range(n_units):
  488. main_tc[j] = np.abs(tc[j][:,max_channels[j]])
  491. best_shifts = align_get_shifts_tc(main_tc, ref, upsample_factor=1)
  492. shifted_tc = shift_chans(tc, best_shifts)
  494. return shifted_tc
  496. def align_get_shifts_tc(wf, ref, upsample_factor = 5, nshifts = 21):
  498. ''' Align all waveforms on a single channel
  500. wf = selected waveform matrix (# spikes, # samples)
  501. max_channel: is the last channel provided in wf
  503. Returns: superresolution shifts required to align all waveforms
  504. - used downstream for linear interpolation alignment
  505. '''
  506. # convert nshifts from timesamples to #of times in upsample_factor
  507. nshifts = (nshifts*upsample_factor)
  508. if nshifts%2==0:
  509. nshifts+=1
  511. # or loop over every channel and parallelize each channel:
  512. #wf_up = []
  513. wf_up = upsample_resample(wf, upsample_factor)
  515. wf_start = 15*upsample_factor
  516. wf_trunc = wf_up[:,wf_start:]
  517. wlen_trunc = wf_trunc.shape[1]
  519. # align to last chanenl which is largest amplitude channel appended
  520. ref_upsampled = upsample_resample(ref[np.newaxis], upsample_factor)[0]
  521. ref_shifted = np.zeros([wf_trunc.shape[1], nshifts])
  523. for i,s in enumerate(range(-int((nshifts-1)/2), int((nshifts-1)/2+1))):
  524. ref_shifted[:,i] = np.roll(ref_upsampled, -s)[wf_start:]
  526. bs_indices = np.matmul(wf_trunc[:,np.newaxis], ref_shifted).squeeze(1).argmax(1)
  527. best_shifts = (np.arange(-int((nshifts-1)/2), int((nshifts-1)/2+1)))[bs_indices]
  529. return best_shifts/np.float32(upsample_factor)
  531. def weighted_histogram(data, weights, bin_range):
  532. bin_counts = np.zeros(len(bin_range)-1)
  533. j = 0
  534. ii = 0
  535. data = data[data < bin_range.max()]
  536. while ii < len(data):
  537. if data[ii] < bin_range[j+1]:
  538. bin_counts[j] += weights[ii]
  539. ii += 1
  540. else:
  541. j += 1
  542. return bin_counts

run.py at commit b18d13a, under Apache-2.0 · at the source

Overview

Authors: Martin O. Bohlen1, Andra M. Rudzite2, Tierney B. Daw1,3, Genevieve M. Kuczewski1, Ergi Spiro4, Cassie Hammond1, Darienne R. Rogers1, Alejandro Gallego-Ortega5, Michael B. Manookin6, Suva Roy2,7, Kimberly Ritola8, Marc A. Sommer1,2,4, Greg D. Field2,5
ORCID iDs: Martin O. Bohlen, Greg D. Field
  1. Department of Biomedical Engineering, Pratt School of Engineering, Duke University, Durham, NC 27708, USA
  2. Department of Neurobiology, Duke University School of Medicine, Durham, NC 27708, USA
  3. Department of Neurobiology, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  4. Center for Cognitive Neuroscience, Duke University, Durham, NC 27708, USA
  5. Jules Stein Eye Institute, Department of Ophthalmology, University California, Los Angeles, Los Angeles, CA 90095, USA
  6. Department of Ophthalmology, University of Washington, Seattle, WA 98109, USA
  7. John A. Moran Eye Center, Department of Ophthalmology and Visual Sciences, University of Utah, Salt Lake City, UT 84132, USA
  8. Department of Pharmacology, University of North Carolina, Chapel Hill, NC 27514, USA
Institutions: Duke University (United States); Duke Medical Center (United States); University of California, Los Angeles (United States); University of Washington (United States); University of Utah (United States); University of North Carolina at Chapel Hill (United States)
Journal: Cell reports methods, volume 6, issue 3, article 101308
Dates: received 2 July 2025; accepted 7 January 2026; published online 5 March 2026; in print March 2026
Type: Research article · Language: English
License: CC BY-NC
Identifiers: DOI 10.1016/j.crmeth.2026.101308 · PMID 41791371 · PMCID PMC13030990 · OpenAlex W7133641340
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: extracellular electrophysiology (units, LFP) (modality), rat (organism), systems (subfield)
Methods: Smoothing, state filtering, decompositions, Machine learning, Evoked potentials, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: optogenetics, AAV, multi-electrode array, superior colliculus, rat, optotagging, retina
MeSH: Optogenetics*, Retinal Ganglion Cells*, Animals, Rats, Superior Colliculi, Visual Pathways (* major topic)
Topic: Retinal Development and Disorders (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Duke Institute for Brain Sciences Incubator Award; NIH (R01 EY034004, K99 EY032119, R01 EY027323, P30 EY005722, P30 EY000331); Research at Duke University
Citations: not cited yet (Europe PMC); 99 references in the paper
Research resources: Goat Anti-ChAT RRID:AB_2079751, Goat Anti-GFP RRID:AB_218182, Normal Donkey Serum RRID:AB_2337258, Donkey Anti-Goat AF 488 RRID:AB_2534102, Donkey Anti-Chicken IgY Biotinylated RRID:AB_2535477, Donkey Anti-Goat AF 555 RRID:AB_2762839, Donkey Anti-Chicken AF 488 RRID:AB_2921070

Abstract

Understanding the structure-function relationships across neurons is challenging, particularly when circuits are composed of dozens of distinct cell types. We refined an approach, called “projection targeting with phototagging”, that allows simultaneous elucidation of the projections, morphology, and visual response properties of diverse retinal ganglion cell (RGC) types in the mammalian retina. The approach combines retrograde virally mediated phototagging of RGCs, microscopy, and large-scale multi-electrode array (MEA) measurements. Importantly, the approach does not rely on transgenic animals and thus is potentially generalizable across species. We validated this approach in rats by targeting retinal projections to the superior colliculus (SC). We showed that multiple RGC types project to the SC and that these results in rats align well with prior findings from transgenic mouse studies.

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

Repositories

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

gdfield/Phototagging_methods

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 1128b36d99c73a7c5c51debc9589ea94b56eac62, 2 September 2025
Languages: MATLAB (19)
Size: 19 files, 19 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (2 files)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
19 files
At the source: github.com/gdfield/Phototagging_methods

Zenodo 17795220

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the text, “Key resources table”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (2 files)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
19 files
At the source: zenodo.org/records/17795220

paninski-lab/yass

License: Apache-2.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: b18d13a69946c1fee28fbc1f67215d3a89d892af, 29 September 2022
Languages: Python (174), Shell (5), C++ (3), CUDA (3), Jupyter (2), C (2)
Size: 269 files, 189 scripts
Software Heritage: archived
Found in: the text, “Key resources table”
Holds: README, license file, environment (requirements.txt, setup.cfg, setup.py, src/gpu_bspline_interp/setup.py, src/gpu_rowshift/setup.py), tests, continuous integration, documentation, 2 notebooks
Not found: CITATION.cff
Tools: NumPy (125 files), SciPy (49 files), PyTorch (27 files), Matplotlib (25 files), scikit-learn (15 files), NetworkX (9 files), h5py (2 files), statsmodels (1 file), TensorFlow (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
192 files
At the source: github.com/paninski-lab/yass

The paper's code and data availability statement is in the Data section.

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

Data and code availability

• Traced retinal ganglion cell morphologies will be deposited at Neuromorpho.org. Accession information is in the key resources table. • Example MEA data will be deposited at DANDI Archive upon publication. Accession information is in the key resources table. Additional data are available upon request from the lead contact. • Example code for identifying ReaChR-positive neurons and plotting EIs over micrographs is available at https://github.com/gdfield/Phototagging_methods (archived at Zenodo). Accession information is in the key resources table. Additional code is available upon request from the lead contact. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 7 keywords, 6 MeSH terms, 3 funders, 95 references, 7 RRIDs.

Cite

This paper

Bohlen, M. O., Rudzite, A. M., Daw, T. B., Kuczewski, G. M., Spiro, E., Hammond, C., Rogers, D. R., Gallego-Ortega, A., Manookin, M. B., Roy, S., Ritola, K., Sommer, M. A., & Field, G. D. (2026). Projection targeting with phototagging to study the structure and function of retinal ganglion cells. Cell reports methods, 6(3), 101308. https://doi.org/10.1016/j.crmeth.2026.101308

BibTeX

@article{bohlen2026projection,
author = {Bohlen, Martin O. and Rudzite, Andra M. and Daw, Tierney B. and Kuczewski, Genevieve M. and Spiro, Ergi and Hammond, Cassie and Rogers, Darienne R. and Gallego-Ortega, Alejandro and Manookin, Michael B. and Roy, Suva and Ritola, Kimberly and Sommer, Marc A. and Field, Greg D.},
title = {{Projection targeting with phototagging to study the structure and function of retinal ganglion cells}},
journal = {Cell reports methods},
year = {2026},
month = mar,
volume = {6},
number = {3},
pages = {101308},
publisher = {Elsevier},
issn = {2667-2375},
doi = {10.1016/j.crmeth.2026.101308},
url = {https://doi.org/10.1016/j.crmeth.2026.101308},
pmid = {41791371},
pmcid = {PMC13030990}
}

RIS

TY - JOUR
AU - Bohlen, Martin O.
AU - Rudzite, Andra M.
AU - Daw, Tierney B.
AU - Kuczewski, Genevieve M.
AU - Spiro, Ergi
AU - Hammond, Cassie
AU - Rogers, Darienne R.
AU - Gallego-Ortega, Alejandro
AU - Manookin, Michael B.
AU - Roy, Suva
AU - Ritola, Kimberly
AU - Sommer, Marc A.
AU - Field, Greg D.
TI - Projection targeting with phototagging to study the structure and function of retinal ganglion cells
T2 - Cell reports methods
J2 - Cell Rep Methods
PY - 2026
DA - 2026/03/05
VL - 6
IS - 3
SP - 101308
SN - 2667-2375
PB - Elsevier
DO - 10.1016/j.crmeth.2026.101308
UR - https://doi.org/10.1016/j.crmeth.2026.101308
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.crmeth.2026.101308",
"type": "article-journal",
"title": "Projection targeting with phototagging to study the structure and function of retinal ganglion cells",
"container-title": "Cell reports methods",
"author": [
{
"family": "Bohlen",
"given": "Martin O."
},
{
"family": "Rudzite",
"given": "Andra M."
},
{
"family": "Daw",
"given": "Tierney B."
},
{
"family": "Kuczewski",
"given": "Genevieve M."
},
{
"family": "Spiro",
"given": "Ergi"
},
{
"family": "Hammond",
"given": "Cassie"
},
{
"family": "Rogers",
"given": "Darienne R."
},
{
"family": "Gallego-Ortega",
"given": "Alejandro"
},
{
"family": "Manookin",
"given": "Michael B."
},
{
"family": "Roy",
"given": "Suva"
},
{
"family": "Ritola",
"given": "Kimberly"
},
{
"family": "Sommer",
"given": "Marc A."
},
{
"family": "Field",
"given": "Greg D."
}
],
"container-title-short": "Cell Rep Methods",
"volume": "6",
"issue": "3",
"page": "101308",
"DOI": "10.1016/j.crmeth.2026.101308",
"PMID": "41791371",
"PMCID": "PMC13030990",
"ISSN": "2667-2375",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.crmeth.2026.101308",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
5
]
]
}
}

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: TensorFlow, NetworkX, h5py, 5 other tools, extracellular electrophysiology (units, LFP), 2 references
[2] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: NetworkX, h5py, statsmodels, 6 other tools, systems, 1 reference
[3] doi:10.3389/fnsys.2026.1822122 [code]
Convergence-divergence circuits for multimodal integration of innate and learned opponent valences.
Journal: Frontiers in systems neuroscience
In common: TensorFlow, NetworkX, h5py, 5 other tools, systems, 1 reference
[4] doi:10.1371/journal.pbio.3003831 [code]
Disinhibitory signaling enables flexible coding of top-down information in cortical networks.
Journal: PLoS biology
In common: TensorFlow, NetworkX, h5py, 6 other tools, systems
[5] doi:10.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: NetworkX, h5py, statsmodels, 6 other tools, systems
[6] doi:10.1093/nar/gkag706 [code]
scDifformer: diffusion-based post-training for virtual cell modeling across large-scale single-cell data.
Journal: Nucleic acids research
In common: TensorFlow, NetworkX, h5py, 6 other tools
[7] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: TensorFlow, NetworkX, statsmodels, 6 other tools
[8] doi:10.1016/j.isci.2026.117088 [code]
Spatial biases in visual feature representation of mouse dorsal lateral geniculate nucleus boutons.
Journal: iScience
In common: Statistics and Machine Learning Toolbox, systems, 5 references
[9] doi:10.1038/s41467-026-72152-x [code]
Centralized brain networks controlling antennal grooming coordination.
Journal: Nature communications
In common: TensorFlow, NetworkX, h5py, 5 other tools, systems
[10] doi:10.1038/s41467-026-75455-1 [code]
Shared latent representations of speech production for cross-patient speech decoding.
Journal: Nature communications
In common: TensorFlow, h5py, statsmodels, 6 other tools

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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