OSCR

Human adherent cortical organoids in a multi-well format.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

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

Jupyter notebook · 565 lines · 27 KB · Apache-2.0

  1. # %%
  2. import pandas as pd
  3. from matplotlib import pyplot as plt
  4. from matplotlib import gridspec as gridspec
  5. from matplotlib.patches import Rectangle
  6. from matplotlib.lines import Line2D
  7. from pathlib import Path
  8. import numpy as np
  10. plt.style.use('seaborn-darkgrid')
  12. ## input folder (seen relatively from this script location, use ../ to go back 1 folder)
  13. input_folder = 'data/adjusted for doubles/'
  14. ## output folder (seen relatively from this script location, use ../ to go back 1 folder)
  15. output_folder = "test_output_folder"
  17. ## spike calling method
  18. # use percentage based threshold of the data distribution instead of a global/local method, if on global/local will not be used! --- default: False
  19. percentage_based_threshold = False
  20. # percentage of data distribution for simple method --- default: 0.99 (min 0.00 & max 1.00)
  21. percentage_threshold = 0.95
  22. # global standard deviation (stddev) based on data, threshold can be set to n times the stddev --- default: 3 (recommended: multiplier not lower than 2.5, needs optimalisation)
  23. global_standard_deviation_threshold_multiplier = 3
  24. # under noise standard deviation based on data and global stddev multiplier, threshold can be set to n times the stddev --- default: 5 (needs optimalisation)
  25. under_noise_standard_deviation_threshold_multiplier = 5
  27. ## spike calling settings
  28. # how many seconds to wait after a spike peak has been to defined before another spike can start to be defined --- default: 3 (needs optimalisation)
  29. spike_offset_time = 3
  30. # time offset in seconds to detect overlap in spikes from other video (as measurements are not on the exact same time scale) --- default: 0.2 (needs optimalisation)
  31. nb_time_offset = 1
  33. ## network burst detection settings
  34. # percentage of how many traces are required to contribute to a Network Burst (NB) to have it count for plotting and statistics --- default: 0.3 (30%, depends on data, n_traces per group)
  35. n_traces_for_NB_per_group_percentage_threshold = 0.3
  36. # absolute amount of traces that are required to contribute to a Network Burst (NB) to have it count for plotting and statistics --- default: 2 (minimally 1, needs optimalisation)
  37. n_traces_for_NB_per_group_threshold = 2
  39. ## plotting settings
  40. # quality control plots, used to check if using stddev_threshold made sense --- default: False
  41. to_plot_distributions = False
  42. # plot equal trace lengths for neatness or false to show each individual trace length --- default: False (Not needed now every group gets its own subplot)
  43. show_equal_trace_length = True
  44. # font size used for any element found on plot --- default: 10
  45. font_size = 24
  46. # to adjust legend positions based on x and y axes respectively, x bigger = to the right, y bigger = upward etc --- default (1.05, .895)
  47. legend_pos_xy = (1.05, .895)
  49. # %%
  50. # create output folder if it doesn't exist
  51. path = Path(output_folder)
  52. path.mkdir(parents=True, exist_ok=True)
  54. # %%
  55. data = Path(input_folder).glob('*.xlsx')
  56. dataframes = []
  58. for xlsx in data:
  59. df = pd.read_excel(xlsx, index_col = 0, engine = 'openpyxl')
  60. # renaming df index for clarity and brevity
  61. index = []
  62. raw_smooth = 'smooth'
  63. for name in df.index.values.tolist():
  64. if 'raw' in name:
  65. raw_smooth = 'raw'
  66. # keep well.video.neuronID
  67. name = name.split(' ')[1]
  68. index.append('{}.{}'.format(name, raw_smooth))
  69. df.index = index
  70. dataframes.append(df)
  72. # merge all dataframes from separate .xlsx files
  73. df_merged = pd.DataFrame()
  74. for df in dataframes:
  75. df_merged = df_merged.append(df)
  76. df = df_merged.T
  78. # DEVNOTE: removed all raw columns, for now not working with raw data!
  79. for name in df.columns:
  80. if 'raw' in name:
  81. del df[name]
  83. # %%
  84. def define_spikes(df, global_sd_multiplier = 3, under_noise_sd_multiplier = 3, spike_offset = 0.5, percentage_based_threshold = True, percentage = 0.99):
  85. start = []
  86. peak = []
  87. # end = [] # DEPRECATED, peak are new ends
  88. index = []
  89. stats_index = []
  90. stats_n_spikes = []
  91. stats_spikes_time = []
  92. stats_seconds_measured = []
  93. threshold_data = []
  94. threshold_indices = []
  95. for (col, data) in df.iteritems():
  96. # create numpy data to enable masking
  97. npdata = np.array(df[col]).astype(np.double)
  98. # create mask to obtain pure trace data
  99. mask = np.isfinite(npdata)
  101. # use simple percentage based on data distribution method or mediocre sophisticated global/local threshold
  102. if percentage_based_threshold:
  103. sorted_data = data[mask].sort_values()
  104. threshold_index = int(percentage*len(data[mask]))
  105. threshold = sorted_data.iloc[threshold_index]
  106. else:
  107. # define 'global' standard deviation - global meaning for whole trace data
  108. global_stddev = np.std(df[col][mask])
  109. # define noise threshold
  110. noise_threshold = global_stddev * global_sd_multiplier
  111. # data under noise threshold, used to calculate under noise stddev
  112. under_noise_threshold = data[mask] < noise_threshold
  113. # define under noise data stddev from all data under noise threshold
  114. under_noise_stddev = np.std(data[mask][under_noise_threshold])
  115. # define threshold
  116. threshold = under_noise_stddev * under_noise_sd_multiplier
  118. # save threshold for plotting
  119. threshold_data.append(threshold)
  120. # trace name
  121. threshold_indices.append(col)
  123. # add a starting point when start_possible and peak_possible and value > threshold
  124. # add peak when trace starting to decreaes, peak_possible and not start_possible
  125. # reset peak_possible and start_possible to true when enough time passed (spike_offset) after peak detected
  126. count = 0
  127. old_value = 0
  128. peak_value = 0
  129. peak_index = 0
  130. peak_possible = True
  131. start_possible = True
  132. for (idx, value) in data[mask].iteritems():
  133. increasing = value > old_value
  134. if value > threshold:
  135. # print(idx, value, peak_index - idx)
  137. if increasing == True and peak_possible == True:
  138. peak_value = value
  139. peak_index = idx
  140. if start_possible == True:
  141. # print('first over threshold')
  142. start.append(idx)
  143. start_possible = False
  144. if increasing == False and peak_possible == True and start_possible == False:
  145. # print('peak found')
  146. # peak found, add as spike
  147. peak_possible = False
  148. peak.append(peak_index)
  149. index.append(col)
  150. count += 1
  151. if increasing == False and (idx - peak_index) > spike_offset and start_possible == False and peak_possible == False:
  152. # print('reset for time')
  153. peak_possible = True
  154. start_possible = True
  155. if value < threshold and (idx - peak_index) > spike_offset and start_possible == False and peak_possible == False:
  156. # print('reset for threshold')
  157. peak_possible = True
  158. start_possible = True
  159. # set to check if trace increasing or decreasing
  160. old_value = value
  162. # extract stats - where idx = total measure time per video
  163. stats_index.append(col)
  164. stats_n_spikes.append(count)
  165. stats_spikes_time.append((count/idx) * 60)
  166. stats_seconds_measured.append(idx)
  168. # # sometimes a spike did not recover ('end') as measuring was cut off, append closing time to ends
  169. if len(start) > len(peak):
  170. peak.append(idx)
  171. index.append(col)
  173. threshold_data = {'threshold': threshold_data}
  174. df_threshold = pd.DataFrame(data = threshold_data, index = threshold_indices)
  176. # initiate stats data to create data frame
  177. stats_data = {'n_spikes': stats_n_spikes, 'spikes_per_minute': stats_spikes_time, 'seconds_measured': stats_seconds_measured}
  178. # create spike stats dataframe
  179. df_stats = pd.DataFrame(data = stats_data, index = stats_index)
  181. data = {'start': start, 'end': peak}
  182. df_spikes = pd.DataFrame(data = data, index = index)
  184. return df_spikes, df_stats, df_threshold
  186. df_spikes, df_stats_spikes, df_threshold = define_spikes(df, global_sd_multiplier = global_standard_deviation_threshold_multiplier, under_noise_sd_multiplier = under_noise_standard_deviation_threshold_multiplier, spike_offset = spike_offset_time, percentage_based_threshold = percentage_based_threshold, percentage = percentage_threshold)
  189. # %%
  190. def define_network_bursts(df_spikes, nb_time_offset, n_traces_for_NB_per_group_percentage_threshold, n_traces_for_NB_per_group_threshold):
  191. # initiate information loop and statistics dictionaries
  192. nb_info = {}
  193. nb_group_stats = {}
  194. nb_trace_stats = {}
  195. nb_trace_interval_mean = {}
  196. nb_trace_interval_sd = {}
  197. nb_group_interval_data = {}
  198. nb_group_interval_mean = {}
  199. nb_group_interval_sd = {}
  200. for ind in df_spikes.index:
  201. well = int(ind.split('.')[0])
  202. video = int(ind.split('.')[1])
  203. trace = int(ind.split('.')[2])
  204. group = '{}.{}'.format(well, video)
  205. if group not in nb_info.keys():
  206. nb_info[group] = 0
  207. nb_group_stats[group] = 0
  208. nb_group_interval_data[group] = []
  209. nb_group_interval_mean[group] = 0
  210. nb_group_interval_sd[group] = 0
  211. old_trace = 0
  212. if trace != old_trace:
  213. nb_info[group] += 1
  214. nb_trace_stats[ind] = 0
  215. nb_trace_interval_mean[ind] = 0
  216. nb_trace_interval_sd[ind] = 0
  217. old_trace = trace
  219. # initiate all NB hits dictionary for plotting NBs later
  220. all_hits = {'name': [], 'start': []}
  221. # loop for each group: g.g.x
  222. for group in nb_info.keys():
  223. # get group data in df
  224. df_group = df_spikes[df_spikes.index.str.startswith(group)]
  225. # while df_group is not empty, take first value and query
  226. while not df_group.empty:
  227. val = df_group['start'][0]
  228. # # check if spike values of all traces within group lie within range based on NB time offset
  229. nb_values = df_group.query('@val - @nb_time_offset < start < @val + @nb_time_offset')
  230. df_group = df_group[~df_group['start'].between(val-nb_time_offset, val+nb_time_offset, inclusive=False)]
  231. if len(nb_values.index) > len(nb_values.index.unique()):
  232. print(1, nb_values)
  233. duplicate = [x for x in nb_values.index if list(nb_values.index).count(x) > 1][0]
  234. ind = list(nb_values.index).index(duplicate)
  235. rename = list(nb_values.index)
  236. rename[ind] = 'delete_row'
  237. nb_values.index = rename
  238. nb_values = nb_values.drop(labels=['delete_row'])
  239. n_traces = 0
  240. for nb_ind in nb_values.index:
  241. n_traces += 1
  242. if n_traces/nb_info[group] >= n_traces_for_NB_per_group_percentage_threshold and n_traces > n_traces_for_NB_per_group_threshold:
  243. nb_group_stats[group] += 1
  244. for nb_ind, nb_start in zip(nb_values.index, nb_values['start']):
  245. nb_trace_stats[nb_ind] += 1
  246. all_hits['name'].append(nb_ind)
  247. all_hits['start'].append(nb_start)
  249. # create NB stats data frame with per trace and per group values
  250. df_stats_network_bursts = pd.DataFrame(data = list(nb_trace_stats.values()), index = list(nb_trace_stats.keys()), columns = ['n_trace_NBs'])
  251. n_group_NBs = []
  252. for name in df_stats_network_bursts.index:
  253. well = int(name.split('.')[0])
  254. video = int(name.split('.')[1])
  255. index = '{}.{}'.format(well, video)
  256. n_group_NBs.append(nb_group_stats[index])
  257. df_stats_network_bursts['n_group_NBs'] = n_group_NBs
  259. # create hits dataframe for plotting
  260. df_hits = pd.DataFrame(data = all_hits['start'], index = all_hits['name'], columns = ['start'])
  263. # calculate trace interval stats
  264. for i in df_hits.index.unique():
  265. group = '{}.{}'.format(int(i.split('.')[0]), int(i.split('.')[1]))
  266. n_trace_NBs = df_stats_network_bursts.loc[i, 'n_trace_NBs']
  267. if n_trace_NBs > 1:
  268. sorted_trace_values = sorted(df_hits.loc[i, 'start'])
  269. # if i == "1.3.2.smooth":
  270. # print(i, sorted_trace_values)
  271. # calculate mean
  272. for v, v2 in zip(sorted_trace_values[1:], sorted_trace_values[:-1]):
  273. nb_trace_interval_mean[i] += v-v2 # sum interval by taking difference of NBtx-NBtx+1
  274. nb_group_interval_data[group].append(v-v2)
  275. nb_trace_interval_mean[i] /= (n_trace_NBs-1) # amount of NB intervals
  276. # calculate standard deviation
  277. for v, v2 in zip(sorted_trace_values[1:], sorted_trace_values[:-1]):
  278. nb_trace_interval_sd[i] += (v-v2 - nb_trace_interval_mean[i])**2 # sum intervals (differences) - respective mean by taking NBtx-NBtx+1
  279. nb_trace_interval_sd[i] /= (n_trace_NBs-1) # amount of NB intervals
  280. nb_trace_interval_sd[i] = nb_trace_interval_sd[i]**(1/2) # take square root
  281. # calculate group interval stats
  282. n_group_NBs = 0
  283. for group in nb_group_interval_data.keys():
  284. # calculate mean
  285. for v in nb_group_interval_data[group]:
  286. nb_group_interval_mean[group] += v
  287. n_group_NBs += 1
  288. if n_group_NBs > 0:
  289. nb_group_interval_mean[group] /= n_group_NBs
  290. # calculate standard deviation
  291. for v in nb_group_interval_data[group]:
  292. nb_group_interval_sd[group] += (v - nb_group_interval_mean[group])**2
  293. if n_group_NBs > 0:
  294. nb_group_interval_sd[group] /= n_group_NBs
  295. nb_group_interval_sd[group] = nb_group_interval_sd[group]**(1/2)
  296. n_group_NBs = 0
  297. # add trace and group meand-sd-coefficient of variation summary statistics to dataframe
  298. df_stats_network_bursts['trace_NB_interval_mean'] = nb_trace_interval_mean.values()
  299. df_stats_network_bursts['trace_NB_interval_sd'] = nb_trace_interval_sd.values()
  300. df_stats_network_bursts['trace_NB_interval_variation'] = (df_stats_network_bursts['trace_NB_interval_sd'] / df_stats_network_bursts['trace_NB_interval_mean']) * 100
  301. nb_group_interval_mean_values = []
  302. nb_group_interval_sd_values = []
  303. for k in df_stats_network_bursts.index:
  304. well = int(k.split('.')[0])
  305. video = int(k.split('.')[1])
  306. group = '{}.{}'.format(well, video)
  307. nb_group_interval_mean_values.append(nb_group_interval_mean[group])
  308. nb_group_interval_sd_values.append(nb_group_interval_sd[group])
  309. df_stats_network_bursts['group_NB_interval_mean'] = nb_group_interval_mean_values
  310. df_stats_network_bursts['group_NB_interval_sd'] = nb_group_interval_sd_values
  311. df_stats_network_bursts['group_NB_interval_variation'] = (df_stats_network_bursts['group_NB_interval_sd'] / df_stats_network_bursts['group_NB_interval_mean']) * 100
  313. # print(df_stats_network_bursts)
  315. return df_hits, df_stats_network_bursts
  317. df_network_bursts, df_stats_network_bursts = define_network_bursts(df_spikes, nb_time_offset = nb_time_offset, n_traces_for_NB_per_group_percentage_threshold = n_traces_for_NB_per_group_percentage_threshold, n_traces_for_NB_per_group_threshold = n_traces_for_NB_per_group_threshold)
  319. # %%
  320. def plot_distribution(df, type, well, output_folder):
  321. custom_colors = [(0,0,0), (146,0,0), (73,0,146), (219,109,0), (0,109,219), (182,109,255), (182,219,255), (146,0,0), (109,182,255),
  322. (36,255,36), (255,182,219), (0,73,73), (255,255,109), (146,73,0), (255,109,182), (58,58,0), (0,146,146)]
  324. # calculate how many traces within condition (correct well and type)
  325. nrows = 0
  326. for col in df.columns:
  327. if well == int(col.split('.')[0]) and type == col.split('.')[3]:
  328. nrows += 1
  330. # instantiate plot
  331. fig, axes = plt.subplots(nrows = nrows, ncols = 1, sharex=True, figsize=(25, 20))
  332. plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.05)
  333. plt.suptitle('Distribution - All videos, well: {}, type: {}'.format(well, type), fontsize=24)
  334. # xlabel
  335. fig.text(0.5, 0.04, 'Timepoints (s)', ha='center', fontsize=24)
  336. # ylabel
  337. fig.text(0.04, 0.5, 'Occurence (%)', va='center', rotation='vertical', fontsize=24)
  339. nrow = 0
  340. for col in df.columns:
  341. if well == int(col.split('.')[0]) and type == col.split('.')[3]:
  342. # get color from color blind friendly color list (custom)
  343. color = list(map(lambda x: x/255, custom_colors[int(col.split('.')[1])]))
  345. # get xdata
  346. xdata = np.array(df[col]).astype(np.double)
  347. xmask = np.isfinite(xdata)
  349. # plot data
  350. axes[nrow].hist(xdata[xmask], bins = 10, label = col, color=color)
  352. # bump index for plotting in correct matplotlib plot figure axes
  353. nrow += 1
  355. # add legend
  356. fig.legend(prop={'size': 20})
  357. # save plot
  358. fig.savefig('{}/Distribution-Well_{}-Type_{}'.format(output_folder, well, type))
  361. # %%
  362. def plot_traces(df, type, well, spikes, network_bursts, thresholds, show_equal_trace_length, font_size, legend_pos_xy, output_folder):
  363. custom_colors = [(0,0,0), (146,0,0), (73,0,146), (219,109,0), (0,109,219), (182,109,255), (182,219,255), (146,0,0), (109,182,255),
  364. (36,255,36), (255,182,219), (0,73,73), (255,255,109), (146,73,0), (255,109,182), (58,58,0), (0,146,146)]
  366. plt.rc('font', size=font_size) # controls default text sizes
  367. plt.rc('axes', titlesize=font_size) # fontsize of the axes title
  368. plt.rc('axes', labelsize=font_size) # fontsize of the x and y labels
  369. plt.rc('xtick', labelsize=font_size) # fontsize of the tick labels
  370. plt.rc('ytick', labelsize=font_size) # fontsize of the tick labels
  371. plt.rc('legend', fontsize=font_size+50) # legend fontsize
  372. plt.rc('figure', titlesize=font_size) # fontsize of the figure title
  374. # calculate how many traces within condition (correct well and type)
  375. nrows = 0
  376. last_valid_trace_indices = []
  377. for col in df.columns:
  378. if well == int(col.split('.')[0]) and type == col.split('.')[3]:
  379. nrows += 1
  380. # from col grab last value, check with index at what final time that is, plot all traces to minimal final time
  381. last_valid_trace_indices.append(df[col].last_valid_index())
  382. minimal_last_trace_index = min(last_valid_trace_indices)
  384. # for specific well, get amount of groups based on final trace of well
  385. group_trace = {}
  386. group_trace_data = {}
  387. for col in df.columns:
  388. wel = int(col.split('.')[0])
  389. if well == wel:
  390. group = int(col.split('.')[1])
  391. trace = int(col.split('.')[2])
  392. group_trace[group] = trace
  393. group_trace_data[str(group) + str(trace)] = df[col]
  395. n_groups = group
  397. # inititate plot
  398. fig = plt.figure(figsize=((25, 20)))
  399. # initiate grid for amount of groups to subplot iteratively
  400. outer = gridspec.GridSpec(n_groups, 1, wspace=0.2, hspace=0.2)
  401. # plot title
  402. plt.suptitle('All videos, well: {}, type: {}'.format(well, type), fontsize=font_size+6)
  403. # xlabel
  404. fig.text(0.5, 0.04, 'Time (s)', ha='center', fontsize=font_size)
  405. # ylabel
  406. ylabel = u'Calcium trace (ΔF/F)' if 'raw' in type else u'Inferred calcium trace (ΔF/F)'
  407. fig.text(0.04, 0.5, ylabel, va='center', rotation='vertical', fontsize=font_size)
  409. # get minimally shared x-axis data
  410. xdata = np.array(df.index).astype(np.double)
  412. legend_lines = []
  413. legend_names = []
  414. # initialize plot grid layout
  415. for i in range(n_groups):
  416. # initiate inner grid (traces) in outer grid (groups)
  417. inner = gridspec.GridSpecFromSubplotSpec(group_trace[i+1], 1,
  418. subplot_spec=outer[i], wspace=0.1, hspace=0.1)
  419. # get color from color blind friendly color list (custom)
  420. color = list(map(lambda x: x/255, custom_colors[i+1]))
  422. legend_lines.append(Line2D([0], [0], color=color, lw=4))
  423. legend_names.append("Cluster " + str(i+1))
  424. for j in range(group_trace[i+1]):
  426. try:
  427. # create dataframe column identifier
  428. tuple = (str(well), str(i+1), str(j+1), type)
  429. col = '.'.join(tuple)
  431. # get fluorescence spiking data
  432. ydata = np.array(group_trace_data[str(i+1)+str(j+1)]).astype(np.double)
  433. # get mask to filter None/nan data (to plot continuous lines)
  434. ymask = np.isfinite(ydata)
  435. # initiate axis in inner grid
  436. ax = plt.Subplot(fig, inner[j])
  437. fig.add_subplot(ax)
  439. # remove x-axis ticks when not last trace of group for visibility
  440. if j+1 != group_trace[i+1]:
  441. ax.set_xticks([])
  442. # plot data
  443. ax.plot(xdata[ymask], ydata[ymask], linestyle='-', marker='', color=color, alpha = .9)
  446. # network bursts: define with lighter colored retangles
  447. for (ind, start) in zip(network_bursts.index, network_bursts.start):
  448. width = 4
  449. if col == ind:
  450. ax.add_patch(Rectangle(xy = (start-width, 0), width = width*2, height = max(ydata[ymask])/1.2, facecolor='grey', alpha=.5))
  452. # spikes: define with lighter colored retangles
  453. for (ind, start, end) in zip(spikes.index, spikes.start, spikes.end):
  454. if col == ind:
  455. # ax.add_patch(Rectangle(xy = (start, 0), width = end-start + 1, height = max(ydata[ymask]), facecolor='black', alpha=0.3))
  456. ax.plot(start, max(ydata[ymask]), marker="*", color='black', alpha=0.3)
  458. # plot used spike threshold
  459. ax.axhline(y=thresholds['threshold'][col], xmin=0.00, xmax=1, color='black', linestyle='dotted', alpha = 0.6)
  461. # set y-axes ticks every 25% of data (max/4)
  462. ax.set_yticks(np.arange(0, max(ydata[ymask]), round(max(ydata[ymask]) / 2, -1)))
  464. if show_equal_trace_length:
  465. # set x-axes limit to minimal_last_trace_index for neat visualization (all traces show equal length of time)
  466. ax.set_xlim([0, minimal_last_trace_index])
  468. except KeyError:
  469. pass
  471. # add 'artists' (class) to legend
  472. threshold_artist = plt.Line2D((0,1),(0,0), color='k', linestyle = 'dotted')
  473. spike_artist = plt.Line2D((0,1),(0,0), color='grey', marker="*", linestyle='')
  474. nb_artist = Rectangle(xy=(0,1),width=1, height = 1, color='black', alpha=.5)
  475. legend_lines.append(threshold_artist)
  476. legend_lines.append(spike_artist)
  477. legend_lines.append(nb_artist)
  478. legend_names.append("Spike threshold")
  479. legend_names.append("Spike")
  480. legend_names.append("Network burst")
  481. fig.legend(legend_lines, legend_names, bbox_to_anchor=[legend_pos_xy[0], legend_pos_xy[1]], loc='upper right', prop={'size': font_size})
  483. # save plot
  484. fig.savefig('{}/Well_{}-Type_{}'.format(output_folder, well, type), bbox_inches='tight')
  487. # %%
  488. # take final column (trace) and extract it's well to get total amount of wells
  489. n_wells = int(df.columns[-1].split('.')[0])
  490. data_types = ['smooth'] # DEVNOTE: for now removed 'raw'
  491. # create separate figures for wells and raw/smooth data types
  492. for n in range(1, n_wells + 1):
  493. for type in data_types:
  494. # TODO reinput defaults: define_spikes = to_define_spikes, define_network_bursts = to_define_network_bursts
  495. plot_traces(df = df, type = type, well = n, spikes = df_spikes, network_bursts = df_network_bursts, thresholds = df_threshold, show_equal_trace_length = show_equal_trace_length, font_size = font_size, legend_pos_xy = legend_pos_xy, output_folder = output_folder)
  496. if to_plot_distributions:
  497. plot_distribution(df = df, type = type, well = n, output_folder = output_folder)
  499. # %%
  500. pd_stats = pd.concat([df_stats_spikes, df_stats_network_bursts], axis = 1)
  502. # calculate amount of spikes per group
  503. group_dict = {}
  504. for ind in pd_stats.index:
  505. well = int(ind.split('.')[0])
  506. video = int(ind.split('.')[1])
  507. key = '{}.{}'.format(well, video)
  508. if key not in group_dict.keys():
  509. group_dict[key] = 0
  510. value = pd_stats.loc[ind, 'n_spikes']
  511. group_dict[key] += value
  512. n_spikes_per_group = []
  513. for ind in pd_stats.index:
  514. well = int(ind.split('.')[0])
  515. video = int(ind.split('.')[1])
  516. key = '{}.{}'.format(well, video)
  517. n_spikes_per_group.append(group_dict[key])
  518. pd_stats['n_group_spikes'] = n_spikes_per_group
  520. # rename n_spikes to n_trace_spikes specifically for clarity
  521. renamed = list(pd_stats.columns)
  522. print(renamed)
  523. renamed[0] = 'n_trace_spikes'
  524. renamed[1] = 'trace_spikes/min'
  525. pd_stats.columns = renamed
  527. print(pd_stats.columns)
  529. # perform summary statistics calculations
  530. pd_stats['group_spikes/min'] = pd_stats['n_group_spikes'] / (pd_stats['seconds_measured'] / 60)
  531. pd_stats['trace_NBs/min'] = pd_stats['n_trace_NBs'] / (pd_stats['seconds_measured'] / 60)
  532. pd_stats['group_NBs/min'] = pd_stats['n_group_NBs'] / (pd_stats['seconds_measured'] / 60)
  533. pd_stats['trace_spike/trace_NBs'] = pd_stats['n_trace_NBs'] / pd_stats['n_trace_spikes']
  534. pd_stats['group_spike/trace_NBs'] = pd_stats['n_trace_NBs'] / pd_stats['n_group_spikes']
  535. pd_stats['trace_spike/group_NBs'] = pd_stats['n_group_NBs'] / pd_stats['n_trace_spikes']
  536. pd_stats['group_spike/group_NBs'] = pd_stats['n_group_NBs'] / pd_stats['n_group_spikes']
  538. # order dataframe for clarity
  539. pd_stats = pd_stats[['seconds_measured', 'n_trace_spikes', 'n_group_spikes', 'n_trace_NBs', 'n_group_NBs',
  540. 'trace_spikes/min', 'group_spikes/min', 'trace_NBs/min', 'group_NBs/min',
  541. 'trace_spike/trace_NBs', 'group_spike/trace_NBs', 'trace_spike/group_NBs', 'group_spike/group_NBs',
  542. 'trace_NB_interval_mean', 'trace_NB_interval_sd', 'trace_NB_interval_variation',
  543. 'group_NB_interval_mean', 'group_NB_interval_sd', 'group_NB_interval_variation']]
  545. pd_stats.to_excel('{}/stats.xlsx'.format(output_folder))
  547. ### conclusions
  548. # If nb_time_offset (too) big, 'NBs' are found based on multiple spikes counting from a single trace (should be impossible!)
  549. # Some spikes seem to be very close to the NB but are not plotted as such: example in 1.2.1 4th-last visible double spike peak, because first gets called as spike and spike_offset_time big enough not to call second one as spike as well, second spike actually seems to belong to the NB based on the other traces!
  550. ## similar examples found in orange traces
  551. # visualize spikes and NBs in one go!
  552. # need to explain the summary stats and need to discuss which to keep and where to document such that it is clear for users!
  554. # %%
  556. # Stats documentatie: kopjes in excel moeten toch alles uitleggen :)
  557. # protocol CNMF-E .xlsx --> input .xlsx
  558. ## selecteer 'well x video x' Excel cell
  559. ## toetsenboord combinatie: shift + down_arrow, voor alle traces
  560. ## toetsenboord combinatie: shift + control + right_arrow, dan kopiëren (control + c)
  561. ## plakken in nieuwe Excel file (control + V), opslaan als bestand met logische naam voor die groep etc
  562. ## alle files met wells.videos.X die dan unieke identifiers hebben (X.X.X C/C_raw) kunnen dan in één input folder, deze ingeven bovenaan script
  563. ### als je te veel traces hebt om dit voor te doen kunnen we kijken of ik toch de CNMF-E data op een andere manier als zelfde geheel kan importeren!
  565. # stappenplan commandos installatie, overdracht code

calcium_imaging.ipynb at commit 33e5e71, under Apache-2.0 · at the source

Overview

Authors: Mark van der Kroeg1, Sakshi Bansal1, Maurits A Unkel1, Hilde Smeenk1, Steven A Kushner1,2,3, Femke MS de Vrij1,4
ORCID iDs: Mark van der Kroeg, Sakshi Bansal, Maurits A Unkel, Hilde Smeenk, Steven A Kushner, Femke MS de Vrij
  1. Department of Psychiatry, Erasmus MC Rotterdam Netherlands
  2. Stavros Niarchos Foundation (SNF) Center for Precision Psychiatry & Mental Health, Columbia University New York United States
  3. Department of Psychiatry, Columbia University Irving Medical Center New York United States
  4. ENCORE Expertise Center for Neurodevelopmental Disorders, Erasmus MC Rotterdam Netherlands
Institutions: Erasmus MC (Netherlands); Columbia University (United States); Columbia University Irving Medical Center (United States)
Journal: eLife, volume 13, article RP98340
Dates: published online 5 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.98340 · PMID 42083454 · PMCID PMC13143278 · OpenAlex W4400884419
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism)
Methods: Statistics, Evoked potentials, Single-unit activity, calcium imaging
Keywords: brain organoids, stem cells, disease modeling, neural differentiation, Human
MeSH: Cell Culture Techniques*, Cerebral Cortex*, Induced Pluripotent Stem Cells*, Neural Stem Cells*, Organoids*, Astrocytes, Cell Differentiation, Humans, Neurons (* major topic)
Journal subjects: Neuroscience, Stem Cells and Regenerative Medicine
Topic: Pluripotent Stem Cells Research (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Nederlandse Organisatie voor Wetenschappelijk Onderzoek (024.003.001); Hersenstichting ((F2012(1)-39); ZonMw (114025201)
Citations: cited by 2 papers (Europe PMC); 55 references in the paper
Research resources: 2012 from hiPSC line IPSC0028 RRID:CVCL_EE38, RRID:CVCL_Y803

Abstract

In the growing diversity of human induced pluripotent stem cell (iPSC)-derived models of brain development, we present here a novel method that exhibits 3D cortical layer formation in a reproducible topography of minimal dimensions. The resulting adherent cortical organoids (ACOs) develop by self-organization after seeding frontal cortex-patterned iPSC-derived neural progenitor cells in 384-well plates during 8 weeks of differentiation. The organoids have stereotypical dimensions of 3 × 3 × 0.2 mm, contain multiple subtypes of neurons, astrocytes, and oligodendrocyte lineage cells, and are amenable to extended culture for at least 10 months. Longitudinal imaging revealed morphologically mature dendritic spines, axonal myelination, and robust neuronal activity. Moreover, ACOs compare favorably to existing free-floating brain organoid models on the basis of robust reproducibility in obtaining topographically standardized radial cortical structures and circumventing internal necrosis. Adherent human cortical organoids hold considerable potential for high-throughput drug discovery applications, neurotoxicological screening, and mechanistic pathophysiological studies of brain disorders.

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

Repository

Its files are read in the Code ↔ Paper reader above.

deVrijLab/calcium_imaging

License: Apache-2.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 33e5e71975987dbf27496de4ab92d60a888ac97f, 16 March 2026
Languages: Jupyter (1), Python (1)
Size: 26 files, 2 scripts
Software Heritage: not archived
Found in: the text, “Calcium imaging”
Holds: license file, 1 notebook
Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
3 files
At the source: github.com/deVrijLab/calcium_imaging

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • 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

No dataset and no data link were found in the paper.

Data availability

Raw calcium imaging videos, ROIs, tracing, and analysis are deposited as an openly available dataset on DataverseNL: https://doi.org/10.34894/5E8AHT. The code for the calcium imaging analyses presented in this paper is openly accessible at https://github.com/deVrijLab/calcium_imaging (copy archived at Unkel, 2026).

The following dataset was generated:

van der KroegM BansalS UnkelM KushnerSA de VrijFMS 2026Calcium imaging adherent cortical organoids Van der Kroeg et al 2026DataverseNL10.34894/5E8AHTPMC1314327842083454

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

Versions

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

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 6 authors, 5 keywords, 9 MeSH terms, 3 funders, 54 references, 2 RRIDs.

Cite

This paper

van der Kroeg, M., Bansal, S., Unkel, M. A., Smeenk, H., Kushner, S. A., & de Vrij, F. M. (2026). Human adherent cortical organoids in a multi-well format. eLife, 13, RP98340. https://doi.org/10.7554/elife.98340

BibTeX

@article{vanderkroeg2026human,
author = {van der Kroeg, Mark and Bansal, Sakshi and Unkel, Maurits A and Smeenk, Hilde and Kushner, Steven A and de Vrij, Femke MS},
title = {{Human adherent cortical organoids in a multi-well format}},
journal = {eLife},
year = {2026},
month = may,
volume = {13},
pages = {RP98340},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.98340},
url = {https://doi.org/10.7554/elife.98340},
pmid = {42083454},
pmcid = {PMC13143278}
}

RIS

TY - JOUR
AU - van der Kroeg, Mark
AU - Bansal, Sakshi
AU - Unkel, Maurits A
AU - Smeenk, Hilde
AU - Kushner, Steven A
AU - de Vrij, Femke MS
TI - Human adherent cortical organoids in a multi-well format
T2 - eLife
J2 - eLife
PY - 2026
DA - 2026/05/05
VL - 13
SP - RP98340
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.98340
UR - https://doi.org/10.7554/elife.98340
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.98340",
"type": "article-journal",
"title": "Human adherent cortical organoids in a multi-well format",
"container-title": "eLife",
"author": [
{
"family": "van der Kroeg",
"given": "Mark"
},
{
"family": "Bansal",
"given": "Sakshi"
},
{
"family": "Unkel",
"given": "Maurits A"
},
{
"family": "Smeenk",
"given": "Hilde"
},
{
"family": "Kushner",
"given": "Steven A"
},
{
"family": "de Vrij",
"given": "Femke MS"
}
],
"container-title-short": "eLife",
"volume": "13",
"page": "RP98340",
"DOI": "10.7554/elife.98340",
"PMID": "42083454",
"PMCID": "PMC13143278",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.98340",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
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.1371/journal.pbio.3003757 [code]
Cell type-agnostic transcriptomic signatures enable uniform comparisons of neural maturation.
Journal: PLoS biology
In common: pandas, Matplotlib, NumPy, 8 references
[2] doi:10.1126/sciadv.aec5080 [code]
Label-free biochemical imaging and time point analysis of neural organoids via deep learning-enhanced Raman microspectroscopy.
Journal: Science advances
In common: Matplotlib, NumPy, 5 references
[3] doi:10.1016/j.cell.2026.07.023
Metabolic atlas of early human cortex reveals glycolytic remodeling and pentose phosphate pathway control of cell fate transitions.
Journal: Cell
In common: 6 references
[4] doi:10.1002/advs.202519893 [code]
NeuroSuite for Long-Term Functional and Structural Studies of Air-Liquid Interface Cerebral Organoids.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: pandas, Matplotlib, NumPy, 4 references
[5] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: 6 references
[6] doi:10.1126/sciadv.adu3955 [code]
Defective EV-mediated transport of SHH alters neural fate specification in EPM1 epilepsy.
Journal: Science advances
In common: pandas, Matplotlib, NumPy, 4 references
[7] doi:10.1038/s41593-026-02316-x [code]
Single-cell multi-omic atlas and morphogen screening informs midbrain and hindbrain organoid engineering.
Journal: Nature neuroscience
In common: pandas, Matplotlib, NumPy, 4 references
[8] doi:10.1038/s41586-026-10877-x [code]
Human brain organoids record the passage of time over multiple years.
Journal: Nature
In common: 5 references
[9] doi:10.1016/j.crmeth.2026.101425
Human cerebral organoids with microglia and vasculature model glioma stem cell interactions and radiotherapy response.
Journal: Cell reports methods
In common: 5 references
[10] doi:10.1038/s41467-026-74320-5 [code]
Spatial architecture of autism pathogenesis reveals mosaic structural disarray during early development.
Journal: Nature communications
In common: pandas, Matplotlib, NumPy, 3 references

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.