catalystneuro / leifer_lab_to_nwb

Conversion scripts for the Leifer lab. Includes the publication Neural signal propagation atlas of Caenorhabditis elegans (Nature, 2023).
BSD 3-Clause "New" or "Revised" License
0 stars 0 forks source link

Details of point spread function #5

Closed CodyCBakerPhD closed 1 month ago

CodyCBakerPhD commented 3 months ago

Extended data Fig 2 shows the point-spread function of the temporal focusing but does not report the mean and standard deviation of a fit; haven't found this in the data share yet either

emosb commented 2 months ago

So depending on how you are viewing the figure, there should be a hyperlink at the bottom (below the caption) where you can download the source data for the extended figures (for this I am accessing online via https://www.nature.com/articles/s41586-023-06683-4#Sec34). There are multiple sheets in the workbook with data corresponding to the different panels of the figure. For Extended Figure 2a there is a few rows of data with the values for the x, y, and z spread shown. The mean and std are not directly reported there but can be easily computed. I can also email you the spreadsheet if that is easier on your end.

CodyCBakerPhD commented 1 month ago

Wonderful, thank you!

CodyCBakerPhD commented 3 days ago

For reference:

import numpy
import pandas
import scipy

# Manually download the .xlsx, import it into Google Sheets, clear the column names and other tabs, and export to .tsv

file_path = "C:/Users/theac/Downloads/point_spread_function.tsv"
table = pandas.read_table(filepath_or_buffer=file_path, header=None)

def estimate_gaussian_parameters(domain: numpy.ndarray, intensity: numpy.ndarray) -> tuple[int, int]:
    absolute_intensity = numpy.absolute(intensity)
    probability = absolute_intensity / scipy.integrate.trapezoid(y=absolute_intensity, x=domain)
    mean = scipy.integrate.trapezoid(y=probability*domain, x=domain)
    var = scipy.integrate.trapezoid(y=probability*(probability-domain)**2, x=domain)
    std = numpy.sqrt(var)
    return mean, std

x_domain = numpy.array(table.iloc[0])
x_intensity = numpy.array(table.iloc[1])  # 'normalized' by what appears to be the peak intensity; perhaps shifted too?
x_mean, x_std = estimate_gaussian_parameters(domain=x_domain, intensity=x_intensity)

# remove nans from array
x_domain = x_domain[~numpy.isnan(x_intensity)]

y_domain = numpy.array(table.iloc[2])
y_intensity = numpy.array(table.iloc[3])
nan_indices = numpy.isnan(y_domain)
y_domain = y_domain[~nan_indices]
y_intensity = y_intensity[~nan_indices]
y_mean, y_std = estimate_gaussian_parameters(domain=y_domain, intensity=y_intensity)

z_domain = numpy.array(table.iloc[4])
z_intensity = numpy.array(table.iloc[5])
nan_indices = numpy.isnan(z_domain)
z_domain = z_domain[~nan_indices]
z_intensity = z_intensity[~nan_indices]
z_mean, z_std = estimate_gaussian_parameters(domain=z_domain, intensity=z_intensity)