
from collections.abc import Iterable
import os

from datetime import datetime, timedelta
from dateutil import tz
from hdmf.backends.hdf5.h5_utils import H5DataIO
from hdmf.container import Container
from hdmf.data_utils import DataChunkIterator
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pynwb import load_namespaces, get_class, register_class, NWBFile, TimeSeries, NWBHDF5IO
from pynwb.file import MultiContainerInterface, NWBContainer, Device, Subject
from pynwb.ophys import ImageSeries, OnePhotonSeries, OpticalChannel, ImageSegmentation, PlaneSegmentation, Fluorescence, DfOverF, CorrectedImageStack, MotionCorrection, RoiResponseSeries, ImagingPlane
from pynwb.core import NWBDataInterface, NWBData
from pynwb.epoch import TimeIntervals
from pynwb.behavior import SpatialSeries, Position
from pynwb.image import ImageSeries
import scipy.io as sio
import skimage.io as skio
from tifffile import TiffFile

# ndx_mulitchannel_volume is the novel NWB extension for multichannel optophysiology in C. elegans
from ndx_multichannel_volume import CElegansSubject, OpticalChannelReferences, OpticalChannelPlus, ImagingVolume, VolumeSegmentation, MultiChannelVolume, MultiChannelVolumeSeries, SegmentationLabels


#Change to location of NWB_tutorial folder on your local system
datapath = '/Users/danielysprague/foco_lab/data/NWB_tutorial/'


nwbfile = NWBFile(
    session_description = 'Add a description for the experiment/session. Can be just long form text',
    #Can use any identity marker that is specific to an individual trial. We use date-time to specify trials
    identifier = '20230322-21-41-10',
    #Specify date and time of trial. Datetime entries are in order Year, Month, Day, Hour, Minute, Second. Not all entries are necessary
    session_start_time = datetime(2023, 3, 22, 21, 41, 10, tzinfo=tz.gettz("US/Pacific")),
    lab = 'FOCO lab',
    institution = 'UCSF',
    related_publications = ''
)

nwbfile #in the latest version of PyNWB, you can output the PyNWB object in Jupyter notebook to easily visualize all of the fields contained within the object


nwbfile.subject = CElegansSubject(
    #This is the same as the NWBFile identifier for us, but does not have to be. It should just identify the subject for this trial uniquely.
    subject_id = '20230322-21-41-10',
    #Age is optional but should be specified in ISO 8601 duration format similarly to what is shown here for growth_stage_time
    #age = pd.Timedelta(hours=2, minutes=30).isoformat(),
    #Date of birth is a required field but if you do not know or if it's not relevant, you can just use the current date or the date of the experiment
    date_of_birth = datetime(2023, 3, 20, tzinfo=tz.gettz("US/Pacific")), 
    #Specify growth stage of worm - should be one of two-fold, three-fold, L1-L4, YA, OA, dauer, post-dauer L4, post-dauer YA, post-dauer OA
    growth_stage = 'YA',
    #Optional: specify time the worm has spent in current growth stage
    growth_stage_time=pd.Timedelta(hours=2, minutes=30).isoformat(),
    #Specify temperature at which animal was cultivated
    cultivation_temp = 20.,
    description = "free form text description, can include whatever you want here",
    #Currently using the ontobee species link until NWB adds support for C. elegans
    species  =  "http://purl.obolibrary.org/obo/NCBITaxon_6239",
    #Currently just using O for other until support added for other gender specifications
    sex = "O", 
    strain = "FC128"
)

nwbfile.subject


device = nwbfile.create_device(
    name = "Spinning disk confocal",
    description = "Leica DMi8 Inverted Microscope with Yokogawa CSU-W1 SoRA, 40x WI objective 1.1 NA",
    manufacturer = "Leica, Yokagawa"
)

device


#Channels is a list of tuples where each tuple contains the fluorophore used, the specific emission filter used, and a short description
#structured as "excitation wavelength - emission filter center point- width of emission filter in nm"
#Make sure this list is in the same order as the channels in your data
channels = [("mTagBFP2", "Chroma ET 460/50", "405-460-50m"), ("CyOFP1", "Chroma ET 605/70","488-605-70m"), ("CyOFP1-high filter", "Chroma ET 700/75","488-700-75m"), ("GFP-GCaMP", "Chroma ET 525/50","488-525-50m"), ("mNeptune 2.5", "Chroma ET 700/75", "561-700-75m"), ("Tag RFP-T", "Chroma ET 605/70", "561-605-70m"), ("mNeptune 2.5-far red", "Chroma ET 700/75", "639-700-75m")]


OptChannels = []
OptChanRefData = []
#The loop below takes the list of channels and converts it into a list of OpticalChannelPlus objects which hold the metadata
#for the optical channels used in the experiment
for fluor, des, wave in channels:
    excite = float(wave.split('-')[0])
    emiss_mid = float(wave.split('-')[1])
    emiss_range = float(wave.split('-')[2][:-1])
    OptChan = OpticalChannelPlus(
        name = fluor,
        description = des,
        excitation_lambda = excite,
        excitation_range = [excite-1.5, excite+1.5],
        emission_range = [emiss_mid-emiss_range/2, emiss_mid+emiss_range/2],
        emission_lambda = emiss_mid
    )

    OptChannels.append(OptChan)
    OptChanRefData.append(wave)


#This object just contains references to the order of channels because OptChannels does not preserve ordering by itself
OpticalChannelRefs = OpticalChannelReferences(
    name = 'OpticalChannelRefs',
    channels = OptChanRefData
)

ImagingVol = ImagingVolume(
    name= 'NeuroPALImVol',
    #Add connections to the OptChannels and OpticalChannelRefs objects
    optical_channel_plus = OptChannels,
    order_optical_channels = OpticalChannelRefs,
    #Free form description of what is being imaged in this volume
    description = 'NeuroPAL image of C. elegans brain',
    #Reference the device created earlier that was used to acquire this data
    device = device,
    #Specifies where in the C. elegans body the image is being taken of
    location = "Head",
    #Specifies the voxel spacing in x, y, z respectively. The values specified should be how many micrometers of physical
    #distance are covered by a single pixel in each dimension
    grid_spacing = [0.3208, 0.3208, 0.75],
    grid_spacing_unit = 'micrometers',
    #Origin coords, origin coords unit, and reference frames are carry over fields from other model organisms where you 
    #are likely only looking at a small portion of the brain. These fields are unfortunately required but feel free to put
    #whatever feels right here
    origin_coords = [0, 0, 0],
    origin_coords_unit = "micrometers",
    reference_frame = "Worm head"
)

nwbfile.add_imaging_plane(ImagingVol) #add this ImagingVol to the nwbfile 

ImagingVol


OpticalChannelRefs


raw_file = datapath+ 'NP_im_raw.tif'

#Data is in order C, Z, Y, X when reading in 
data = skio.imread(raw_file)
print(data.shape)
data = np.transpose(data) #Tranpose data to proper order
print(data.shape)

RGBW_channels = [0,1,4,6]

Image = MultiChannelVolume(
    name = 'NeuroPALImageRaw',
    #This is the same OpticalChannelRefs used in the associated Imaging Volume
    order_optical_channels = OpticalChannelRefs,
    description = 'free form description of image',
    #Specifies which channels in the image are associated with the RGBW channels - should be a list of channel indices as shown above
    RGBW_channels = RGBW_channels,
    #This is the raw data numpy array that we loaded above
    data = H5DataIO(data=data, compression=True),
    #This is a reference to the Imaging Volume object we defined previously
    imaging_volume = ImagingVol
)

Image

(7, 48, 240, 1000)
(1000, 240, 48, 7)


processed_file = datapath+'NP_im_proc.tif'

proc_data = np.transpose(skio.imread(processed_file),[2,1,0,3])

print(proc_data.shape)

proc_channels = [channels[i] for i in RGBW_channels]

ProcOptChannels = []
ProcOptChanRefData = []
for fluor, des, wave in proc_channels:
    excite = float(wave.split('-')[0])
    emiss_mid = float(wave.split('-')[1])
    emiss_range = float(wave.split('-')[2][:-1])
    OptChan = OpticalChannelPlus(
        name = fluor,
        description = des,
        excitation_lambda = excite,
        excitation_range = [excite-1.5, excite+1.5],
        emission_range = [emiss_mid-emiss_range/2, emiss_mid+emiss_range/2],
        emission_lambda = emiss_mid
    )

    ProcOptChannels.append(OptChan)
    ProcOptChanRefData.append(wave)


ProcOpticalChannelRefs = OpticalChannelReferences(
    name = 'OpticalChannelRefs',
    channels = ProcOptChanRefData
)

ProcImagingVol = ImagingVolume(
    name= 'ProcessedImVol',
    optical_channel_plus = ProcOptChannels,
    order_optical_channels = ProcOpticalChannelRefs,
    description = 'NeuroPAL image of C. elegans brain',
    device = device,
    location = "Head",
    grid_spacing = [0.3208, 0.3208, 0.75],
    grid_spacing_unit = 'micrometers',
    origin_coords = [0, 0, 0],
    origin_coords_unit = "micrometers",
    reference_frame = "Worm head"
)

ProcImage = MultiChannelVolume(
    name = 'ProcessedImage',
    order_optical_channels = ProcOpticalChannelRefs,
    description = 'Can describe the specific pre-processing steps taken here',
    RGBW_channels = [0,1,2,3],
    data = proc_data,
    imaging_volume = ProcImagingVol
)

nwbfile.add_imaging_plane(ProcImagingVol)

ProcImage

(1000, 240, 48, 4)


blob_file =  datapath + 'blobs.csv'
blobs = pd.read_csv(blob_file)

IDs = blobs['ID']
labels = IDs.replace(np.nan,'',regex=True)
labels = list(np.asarray(labels)) 
#Need list of ID labels to add as extra column to PlaneSegmentation

vs = PlaneSegmentation(
    name = 'NeuroPALNeurons',
    description = 'Neuron centers for multichannel volumetric image. Weight set at 1 for all voxels. Labels refers to cell ID of segmented neurons',
    #Reference the same ImagingVolume that your image was taken with
    imaging_plane = ImagingVol,
)

#Uncomment the code below to use image_mask instead of voxel_mask
'''
for i, row in blobs.iterrows():
    #image_mask should have the same size as the x, y, z dimensions of the raw image
    image_mask = np.zeros(data.shape[:-1])
    x = row['X']
    y = row['Y']
    z = row['Z']
    ID = row['ID']

    #we define the ROI here as a 3x3x3 region surrounding the neuron center
    image_mask[x-1:x+1, y-1:y+1, z-1:z+1] = 1
    vs.add_roi(image_mask=image_mask)

'''

#Uncomment the code below to use voxel_mask instead of image_mask
for i, row in blobs.iterrows():
    voxel_mask = []
    x = row['X']
    y = row['Y']
    z = row['Z']

    voxel_mask.append([np.uint(x),np.uint(y),np.uint(z),1]) #If using a segmentation mask rather than neuron center, voxel_mask should be a list of all voxels associated with the ROI

    vs.add_roi(voxel_mask=voxel_mask)

vs.add_column(
     name = 'ID_labels',
     description = 'ROI ID labels',
     data = labels,
     index=True, 
)

NeuroPALImSeg = ImageSegmentation(
    name = 'NeuroPALSegmentation',
)

NeuroPALImSeg.add_plane_segmentation(vs)

NeuroPALImSeg


NeuroPALImSeg = MultiChannelVolume(
    name = 'NeuroPALSegmentation',
    imaging_volume = ImagingVol, #can just use the same as the associated raw image
    order_optical_channels = OpticalChannelRefs, #can just use the same as the associated raw image
    description = 'Indexed mask of neuron ROIs',
    RGBW_channels = RGBW_channels,
    data = H5DataIO(data=data, compression=True), 
    #data here should be an array with the same shape as the raw image with values as described in the markdown above
)

Seglabels = SegmentationLabels(
    name = 'NeuroPALLabels',
    labels = labels, #should be array of text labels with the same length as the number of ROIs
    description = 'NeuroPAL segmentation labels',
    #Use one of the below to link either the ImageSegmentation, MultiChannelVolume, or MultiChannelVolumeSeries object where
    #associated ROIs are defined
    #ImageSegmentation = NeuroPALImSeg, 
    MCVSegmentation = NeuroPALImSeg,
    #MCVSeriesSegmentation = NeuroPALImSeg,
)

Seglabels

['',
 '',
 '',
 'IL1VL',
 '',
 '',
 'IL2DR',
 'IL2VR',
 'IL2R',
 'URADR',
 '',
 'IL1R',
 'IL1VR',
 'IL1DL',
 'URYVR',
 'IL1DR',
 'OLLR',
 '',
 'URBR',
 'OLQVL',
 'URYDR',
 '',
 '',
 'OLQVR',
 'OLQDL',
 '',
 'RMER',
 '',
 'RMED',
 'OLQDR',
 '',
 'RMEV',
 'RMDVR',
 '',
 '',
 'RMDVL',
 '',
 'AVAR',
 '',
 '',
 'AVAL',
 '',
 '',
 '',
 'ALA',
 '',
 'RMDR',
 '',
 'SMDVR',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 '',
 'AVJR',
 '',
 '',
 'AVDR',
 '',
 '',
 '',
 '',
 '',
 '',
 'FLPR',
 '',
 'ADER',
 '',
 'AQR',
 'ADAR']


from tifffile import TiffFile

tif = TiffFile(datapath+'Calc_series_raw.tiff')

print(len(tif.pages))

page = tif.pages[0]

print(page.shape)
print(page.dtype)
print(page.axes)

series = tif.series[0]
print(series.shape)
print(series.axes)

print(len(tif.pages))

tif.close()

18000
(240, 1000)
uint16
YX
(18000, 240, 1000)
QYX
18000


#define a data generator function that will yield a single data entry, in our case we are iterating over time points and creating a Z stack of images for each time point
def iter_calc_tiff(filename, numZ):

    #TiffFile object allows you to access metadata for the tif file and selectively load individual pages/series
    tif = TiffFile(filename)

    #In this dataset, one page is one XY plane and every 12 pages comprises one Z stack for an individual time point
    pages = len(tif.pages)
    timepoints = int(pages/numZ)

    pageshape = tif.pages[0].shape

    #We iterate through all of the timepoints and yield each timepoint back to the DataChunkIterator
    for i in range(timepoints):
        tpoint = np.zeros((pageshape[1],pageshape[0], numZ))
        for j in range(numZ):
            image = np.transpose(tif.pages[i*numZ+j].asarray())
            tpoint[:,:,j] = image

        #Make sure array ends up as the correct dtype coming out of this function (the dtype that your data was collected as)
        yield tpoint.astype('uint16')

    tif.close()

    return

#The DataChunkIterator wraps the data generator function and will stitch together the chunks as it iteratively reads over the full file
data = DataChunkIterator(
    data= iter_calc_tiff(datapath+'Calc_series_raw.tiff', 12),
    #this will be the max shape of the final image. Can leave blank or set as the size of your full data if you know that ahead of time
    maxshape = None,
    buffer_size = 10,
)

CalcChannels = [("GFP-GCaMP", "Chroma ET 525/50","488-525-50m")]

CalcOptChannels = []
CalcOptChanRefData = []
#The loop below takes the list of channels and converts it into a list of OpticalChannelPlus objects which hold the metadata
#for the optical channels used in the experiment
for fluor, des, wave in CalcChannels:
    excite = float(wave.split('-')[0])
    emiss_mid = float(wave.split('-')[1])
    emiss_range = float(wave.split('-')[2][:-1])
    OptChan = OpticalChannelPlus(
        name = fluor,
        description = des,
        excitation_lambda = excite,
        excitation_range = [excite-1.5, excite+1.5],
        emission_range = [emiss_mid-emiss_range/2, emiss_mid+emiss_range/2],
        emission_lambda = emiss_mid
    )

    CalcOptChannels.append(OptChan)
    CalcOptChanRefData.append(wave)


#This object just contains references to the order of channels because OptChannels does not preserve ordering
CalcOptChanRefs = OpticalChannelReferences(
    name = 'OpticalChannelRefs',
    channels = CalcOptChanRefData
)

CalcImagingVolume = ImagingVolume(
    name = 'CalciumImVol',
    description = 'Imaging volume used to acquire calcium imaging data',
    optical_channel_plus = CalcOptChannels,
    order_optical_channels = CalcOptChanRefs,
    device = device,
    location = 'Worm head',
    grid_spacing = [0.3208, 0.3208, 2.5],
    grid_spacing_unit = 'um',
    reference_frame = 'Worm head'
)

wrapped_data = H5DataIO(data=data, compression="gzip", compression_opts=4)

calcium_image_series = MultiChannelVolumeSeries(
    name="CalciumImageSeries",
    description = "Raw GCaMP series images",
    comments = "Include information about what each channel represents here: e.g. GFP-GCaMP channel represents GCaMP signal, mNeptune represents reference signal",
    data=wrapped_data,
    device = device,
    unit="Voxel gray counts", 
    scan_line_rate = 2995.,
    dimension = [1000,240,12],
    resolution = 1., #smallest meaningful difference (in specified unit) between values in data: i.e. level of precision
    rate = 1.04, #sampling rate in hz
    imaging_volume = CalcImagingVolume,
)


nwbfile.add_imaging_plane(CalcImagingVolume)

calcium_image_series


gce_quant = pd.read_csv(datapath+'gce_quant.csv')

print(len(gce_quant['blob_ix'].unique())) #Count the number of unique blobs in this file
print(len(gce_quant['T'].unique())) #Count the number of unique time points in this file

quant = gce_quant[['X', 'Y', 'Z', 'gce_quant', 'ID', 'T', 'blob_ix']] #Reorder columns to order we want

blobquant = None
for idx in gce_quant['blob_ix'].unique():
    blob = quant[quant['blob_ix']==idx]
    blobarr = np.asarray(blob[['X','Y','Z','gce_quant','ID']]) 
    blobarr = blobarr[np.newaxis, :, :]
    if blobquant is None:
        blobquant=blobarr

    else:
        blobquant = np.vstack((blobquant, blobarr))

print(blobquant.shape) #Now dimensions are blob_ix, time, and data columns (X, Y, Z, gce_quant, ID). We are now ready to add this data to NWB objects.

volsegs = []

for t in range(blobquant.shape[1]):
    volseg = PlaneSegmentation(
        name = 'Seg_tpoint_'+str(t),
        description = 'Neuron segmentation for time point ' +str(t) + ' in calcium image series',
        imaging_plane = CalcImagingVolume,
        reference_images = calcium_image_series,
    )

    for i in range(blobquant.shape[0]):
        voxel_mask = blobquant[i,t,0:3] #X, Y, Z columns
        if np.any(np.isnan(voxel_mask)):
            voxel_mask = np.asarray([0,0,0,0]) #if blob does not exist at time point (nan values in row) we replace values with 0 and set weight to 0
        else:
            voxel_mask = np.hstack((voxel_mask, 1)) #add weight of one to each blob
        voxel_mask = voxel_mask[np.newaxis,:] #add empty new axis to make shape compatible

        volseg.add_roi(voxel_mask=voxel_mask)

    volsegs.append(volseg)

ImSeg = ImageSegmentation(
    name = 'CalciumSeriesSegmentation', #use if tracking neurons across frames (correspondence between segmentations)
    #name = 'CalciumSeriesSegmentationUntracked' #use if not trakcing across frames (ie raw segmentation in each frame)
    plane_segmentations = volsegs
)

gce_data = np.transpose(blobquant[:,:,3]) #Take only gce quantification column and transpose so time is in the first dimension

rt_region = volsegs[0].create_roi_table_region(
    description = 'All segmented neurons associated with calcium image series',
    region = list(np.arange(blobquant.shape[0]))
)

#If you have raw fluorescence values rather than DFoF use the Fluorescence object instead of the DfOverF object to save your RoiResponseSeries
SignalRoiResponse = RoiResponseSeries(
    #See https://pynwb.readthedocs.io/en/stable/pynwb.ophys.html#pynwb.ophys.RoiResponseSeries for additional key word argument options
    name = 'SignalCalciumImResponseSeries',
    description = 'DF/F activity for calcium imaging data',
    data = gce_data, #first dimension should represent time and second dimension should represent ROIs
    rois = rt_region,
    unit = 'Percentage', #the unit of measurement for the data input here
    resolution = 0.01,
    rate = 1.04
)

SignalFluor = DfOverF( #Change to Fluorescence if using raw fluorescence
    name = 'SignalDFoF', #Change name to SignalRawFluor if using raw fluorescence, rename reference and processed object accordingly
    roi_response_series = SignalRoiResponse
)

'''

RefRoiResponse = RefRoiResponseSeries(
    #See https://pynwb.readthedocs.io/en/stable/pynwb.ophys.html#pynwb.ophys.RoiResponseSeries for additional key word argument options
    name = 'ReferenceCalciumImResponseSeries',
    description = 'Fluorescence for reference channel in calcium imaging',
    data = ref_data, #first dimension should represent time and second dimension should represent ROIs
    rois = rt_region, 
    unit = '', #the unit of measurement for the data input here
    rate = 4.0
)

RefFluor = DfOverF(
    name = 'ReferenceDFoF',
    roi_response_series = RefRoiResponse
)

ProcRoiResponse = RoiResponseSeries(
    #See https://pynwb.readthedocs.io/en/stable/pynwb.ophys.html#pynwb.ophys.RoiResponseSeries for additional key word argument options
    name = 'ProcessedCalciumImResponseSeries',
    description = 'processed calcium fluorescence activity - describe processing steps taken to obtain',
    data = proc_gce_data, #first dimension should represent time and second dimension should represent ROIs
    rois = rt_region,
    unit = '', #the unit of measurement for the data input here
    rate = 4.0
)

ProcFluor = DfoF(
    name = 'ProcessedDFoF',
    roi_response_series = RoiResponse
)

'''

ImSeg

229
1500
(229, 1500, 5)


CalciumSegSeries = MultiChannelVolumeSeries(
    name="CalciumSeriesSegmentation",
    description = "Series of indexed masks associated with calcium segmentation",
    comments = "Include here whether ROIs are tracked across frames or any other comments",
    data=data, # data here should be series of indexed masks
    #Elements below can be kept the same as the CalciumImageSeries defined above
    device = device,
    unit="Voxel gray counts", 
    scan_line_rate = 2995.,
    dimension = [1000,240,12],
    resolution = 1., #smallest meaningful difference (in specified unit) between values in data: i.e. level of precision
    rate = 1.04, #sampling rate in hz
    imaging_volume = CalcImagingVolume,
)

FirstFrameSeg = PlaneSegmentation(
        name = 'Seg_tpoint_'+str(t),
        description = 'Neuron segmentation for time point ' +str(t) + ' in calcium image series',
        imaging_plane = CalcImagingVolume,
        reference_images = calcium_image_series,
)

for i in range(blobquant.shape[0]):
    voxel_mask = blobquant[i,0,0:3] #X, Y, Z columns for time point 0
    if np.any(np.isnan(voxel_mask)):
        voxel_mask = np.asarray([0,0,0,0]) #if blob does not exist at time point (nan values in row) we replace values with 0 and set weight to 0
    else:
        voxel_mask = np.hstack((voxel_mask, 1)) #add weight of one to each blob
    voxel_mask = voxel_mask[np.newaxis,:] #add empty new axis to make shape compatible

    FirstFrameSeg.add_roi(voxel_mask=voxel_mask)

rt_region = FirstFrameSeg.create_roi_table_region( #add ROIResponseSeries as defined above using this rt_region
    description = 'All segmented neurons associated with calcium image series',
    region = list(np.arange(blobquant.shape[0]))
)


#we add our raw NeuroPAL image to the acquisition module of the base NWB file
nwbfile.add_acquisition(Image)
nwbfile.add_acquisition(calcium_image_series)

#we create a processing module for our neuroPAL data 
neuroPAL_module = nwbfile.create_processing_module(
    name = 'NeuroPAL',
    description = 'NeuroPAL image metadata and segmentation'
)

neuroPAL_module.add(NeuroPALImSeg)
#neuroPAL_module.add(Seglabels) #optional, include if defining labels in separate SegmentationLabels object
neuroPAL_module.add(OpticalChannelRefs)

#we create a processing module for the pre-processed neuroPAL image
processed_im_module = nwbfile.create_processing_module(
    name = 'ProcessedImage',
    description = 'Data and metadata associated with the processed neuroPAL image'
)

processed_im_module.add(ProcImage)
processed_im_module.add(ProcOpticalChannelRefs)

ophys = nwbfile.create_processing_module(
    name = 'CalciumActivity',
    description = 'Calcium time series metadata, segmentation, and fluorescence data'
)

ophys.add(ImSeg)
#ophys.add(CalciumSegSeries) # comment out above line and uncomment this one if using indexed mask approach
#ophys.add(FirstFrameSeg) # uncomment if using indexed mask approach
ophys.add(SignalFluor)
ophys.add(CalcOptChanRefs)
#ophys.add(RefFluor)
#ophys.add(ProcFluor)


#each NWB file should be named with a unique identifier 
identifier = '20230322-21-41-10'

#specify the file path you want to save this NWB file to 
io = NWBHDF5IO(datapath+'example.nwb', mode='w')
io.write(nwbfile)
io.close()

/Applications/anaconda3/envs/NWB-dev/lib/python3.11/site-packages/hdmf/build/objectmapper.py:259: DtypeConversionWarning: Spec 'ImagingVolume/origin_coords': Value with data type int64 is being converted to data type float64 (min specification: float32).
  warnings.warn(full_warning_msg, DtypeConversionWarning)
/Applications/anaconda3/envs/NWB-dev/lib/python3.11/site-packages/hdmf/build/objectmapper.py:668: MissingRequiredBuildWarning: ImagingVolume 'NeuroPALImVol' is missing required value for attribute 'optical_channel'.
  warnings.warn(msg, MissingRequiredBuildWarning)
/Applications/anaconda3/envs/NWB-dev/lib/python3.11/site-packages/hdmf/build/objectmapper.py:668: MissingRequiredBuildWarning: ImagingVolume 'CalciumImVol' is missing required value for attribute 'optical_channel'.
  warnings.warn(msg, MissingRequiredBuildWarning)
/Applications/anaconda3/envs/NWB-dev/lib/python3.11/site-packages/hdmf/build/objectmapper.py:668: MissingRequiredBuildWarning: ImagingVolume 'ProcessedImVol' is missing required value for attribute 'optical_channel'.
  warnings.warn(msg, MissingRequiredBuildWarning)


with NWBHDF5IO(datapath+'example.nwb', mode='r', load_namespaces=True) as io:
    read_nwbfile = io.read()
    print(read_nwbfile) # you can print the nwbfile after reading to see the contents of the file
    print(read_nwbfile.processing['CalciumActivity'])
    print(read_nwbfile.processing['NeuroPAL'])
    print(read_nwbfile.processing['ProcessedImage'])
    subject = read_nwbfile.subject #get the metadata about the experiment subject
    growth_stage = subject.growth_stage
    image = read_nwbfile.acquisition['NeuroPALImageRaw'].data[:] #get the neuroPAL image as a np array
    channels = read_nwbfile.acquisition['NeuroPALImageRaw'].RGBW_channels[:] #get which channels of the image correspond to which RGBW pseudocolors
    im_vol = read_nwbfile.acquisition['NeuroPALImageRaw'].imaging_volume #get the metadata associated with the imaging acquisition
    seg = read_nwbfile.processing['NeuroPAL']['NeuroPALSegmentation']['NeuroPALNeurons'].voxel_mask[:] #get the locations of neuron centers
    labels = read_nwbfile.processing['NeuroPAL']['NeuroPALSegmentation']['NeuroPALNeurons']['ID_labels'][:]
    optchans = im_vol.optical_channel_plus[:] #get information about all of the optical channels used in acquisition
    chan_refs = read_nwbfile.processing['NeuroPAL']['OpticalChannelRefs'].channels[:] #get the order of the optical channels in the image
    proc_image = read_nwbfile.processing['ProcessedImage']['ProcessedImage'].data[:] #get the pre-processed image
    calcium_frames = read_nwbfile.acquisition['CalciumImageSeries'].data[0:15, :,:,:] #load the first 15 frames of the calcium images
    fluor = read_nwbfile.processing['CalciumActivity']['SignalDFoF']['SignalCalciumImResponseSeries'].data[:]
    #calc_seg = read_nwbfile.processing['CalciumActivity']['CalciumSeriesSegmentation']['Seg_tpoint_0'].voxel_mask[:]

/Applications/anaconda3/envs/NWB-dev/lib/python3.11/site-packages/hdmf/spec/namespace.py:531: UserWarning: Ignoring cached namespace 'ndx-multichannel-volume' version 0.1.5 because version 0.1.7 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."

root pynwb.file.NWBFile at 0x12094531344
Fields:
  acquisition: {
    CalciumImageSeries <class 'ndx_multichannel_volume.ndx_multichannel_volume.MultiChannelVolumeSeries'>,
    NeuroPALImageRaw <class 'ndx_multichannel_volume.ndx_multichannel_volume.MultiChannelVolume'>
  }
  devices: {
    Spinning disk confocal <class 'pynwb.device.Device'>
  }
  file_create_date: [datetime.datetime(2023, 8, 18, 13, 3, 32, 838112, tzinfo=tzoffset(None, -25200))]
  identifier: 20230322-21-41-10
  imaging_planes: {
    CalciumImVol <class 'ndx_multichannel_volume.ndx_multichannel_volume.ImagingVolume'>,
    NeuroPALImVol <class 'ndx_multichannel_volume.ndx_multichannel_volume.ImagingVolume'>,
    ProcessedImagingVolume <class 'ndx_multichannel_volume.ndx_multichannel_volume.ImagingVolume'>
  }
  institution: UCSF
  lab: FOCO lab
  processing: {
    CalciumActivity <class 'pynwb.base.ProcessingModule'>,
    NeuroPAL <class 'pynwb.base.ProcessingModule'>,
    ProcessedImage <class 'pynwb.base.ProcessingModule'>
  }
  related_publications: ['']
  session_description: Add a description for the experiment/session. Can be just long form text
  session_start_time: 2023-03-22 21:41:10-07:00
  subject: subject abc.CElegansSubject at 0x12096049808
Fields:
  age__reference: birth
  cultivation_temp: 20.0
  date_of_birth: 2023-03-20 00:00:00-07:00
  description: free form text description, can include whatever you want here
  growth_stage: YA
  growth_stage_time: P0DT2H30M0S
  sex: O
  species: http://purl.obolibrary.org/obo/NCBITaxon_6239
  strain: FC128
  subject_id: 20230322-21-41-10

  timestamps_reference_time: 2023-03-22 21:41:10-07:00

CalciumActivity pynwb.base.ProcessingModule at 0x12094802064
Fields:
  data_interfaces: {
    CalciumSeriesSegmentation <class 'pynwb.ophys.ImageSegmentation'>,
    OpticalChannelRefs <class 'abc.OpticalChannelReferences'>,
    SignalDFoF <class 'pynwb.ophys.DfOverF'>
  }
  description: Calcium time series metadata, segmentation, and fluorescence data

NeuroPAL pynwb.base.ProcessingModule at 0x12094814608
Fields:
  data_interfaces: {
    NeuroPALSegmentation <class 'pynwb.ophys.ImageSegmentation'>,
    OpticalChannelRefs <class 'abc.OpticalChannelReferences'>
  }
  description: NeuroPAL image metadata and segmentation

ProcessedImage pynwb.base.ProcessingModule at 0x12096062672
Fields:
  data_interfaces: {
    OpticalChannelRefs <class 'abc.OpticalChannelReferences'>,
    ProcessedImage <class 'ndx_multichannel_volume.ndx_multichannel_volume.MultiChannelVolume'>
  }
  description: Data and metadata associated with the processed neuroPAL image


labels = ["".join(label) for label in labels]
print(labels)

['', '', '', 'IL1VL', '', '', 'IL2DR', 'IL2VR', 'IL2R', 'URADR', '', 'IL1R', 'IL1VR', 'IL1DL', 'URYVR', 'IL1DR', 'OLLR', '', 'URBR', 'OLQVL', 'URYDR', '', '', 'OLQVR', 'OLQDL', '', 'RMER', '', 'RMED', 'OLQDR', '', 'RMEV', 'RMDVR', '', '', 'RMDVL', '', 'AVAR', '', '', 'AVAL', '', '', '', 'ALA', '', 'RMDR', '', 'SMDVR', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', 'AVJR', '', '', 'AVDR', '', '', '', '', '', '', 'FLPR', '', 'ADER', '', 'AQR', 'ADAR']


blobs = pd.DataFrame.from_records(seg, columns = ['X', 'Y', 'Z', 'weight', 'ID'])
blobs = blobs.drop(['weight'], axis=1)
blobs = blobs.replace('nan', np.nan, regex=True) 

print(proc_image.shape)

RGB = proc_image[:,:,:,:-1]

print(RGB.shape)

Zmax = np.max(RGB, axis=2)
Ymax = np.max(RGB, axis=1)

plt.figure()

plt.imshow(np.transpose(Zmax, [1,0,2]))
plt.scatter(blobs['x'], blobs['y'], s=5)
plt.xlim((0, Zmax.shape[0]))
plt.ylim((0, Zmax.shape[1]))
plt.gca().set_aspect('equal')

plt.show()

plt.figure()

plt.imshow(np.transpose(Ymax, [1,0,2]))
plt.scatter(blobs['x'], blobs['z'], s=5)
plt.xlim((0, Ymax.shape[0]))
plt.ylim((0, Ymax.shape[1]))
plt.gca().set_aspect('equal')

plt.show()

(1000, 240, 48, 4)
(1000, 240, 48, 3)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).


#TODO: implement example of editing existing NWB file

NWB conversion tutorial¶

Data used in this tutorial¶

Imports and setup¶

Setting up the NWB file and adding basic metadata¶

root (NWBFile)

subject (CElegansSubject)

Spinning disk confocal (Device)

Adding NeuroPAL images¶

NeuroPALImVol (ImagingVolume)

OpticalChannelRefs (OpticalChannelReferences)

NeuroPALImageRaw (MultiChannelVolume)

ProcessedImage (MultiChannelVolume)

Adding segmentation and ID results¶

NeuroPALSegmentation (ImageSegmentation)

Adding calcium imaging data¶

CalciumImageSeries (MultiChannelVolumeSeries)

CalciumSeriesSegmentation (ImageSegmentation)

Writing the NWB file¶

Reading from existing NWB files¶

Editing existing NWB file¶