diff --git a/README.md b/README.md index 77dcd364..3635783e 100644 --- a/README.md +++ b/README.md @@ -91,7 +91,7 @@ The software has been heavily tested on Windows 10 and Ubuntu 18.04, and less we 2. Run `conda env create --name suite2p` 3. To activate this new environment, run `conda activate suite2p` (you will have to activate every time you want to run suite2p) 4. Install the local version of suite2p into this environment in develop mode with the command `pip install -e .` -5. Run tests: `python setup.py test` or `pytest tests\`, this will automatically download the test data into your `suite2p` folder. The test data is split into two parts: test inputs and expected test outputs which will be downloaded in `data/test_inputs` and `data/test_outputs` respectively. The .zip files for these two parts can be downloaded from these links: [test_inputs](https://www.suite2p.org/static/test_data/test_inputs.zip) and [test_outputs](https://www.suite2p.org/static/test_data/test_outputs.zip). +5. Run tests: `python setup.py test` or `pytest -vs`, this will automatically download the test data into your `suite2p` folder. The test data is split into two parts: test inputs and expected test outputs which will be downloaded in `data/test_inputs` and `data/test_outputs` respectively. The .zip files for these two parts can be downloaded from these links: [test_inputs](https://www.suite2p.org/static/test_data/test_inputs.zip) and [test_outputs](https://www.suite2p.org/static/test_data/test_outputs.zip). ## Examples diff --git a/conftest.py b/conftest.py index 68f95599..6109c22b 100644 --- a/conftest.py +++ b/conftest.py @@ -81,6 +81,7 @@ def download_url_to_file(url, dst, progress=True): dst = os.path.expanduser(dst) dst_dir = os.path.dirname(dst) f = tempfile.NamedTemporaryFile(delete=False, dir=dst_dir) + print(f"\nDownloading: {url}") try: with tqdm(total=file_size, disable=not progress, unit='B', unit_scale=True, unit_divisor=1024) as pbar: diff --git a/suite2p/io/nwb.py b/suite2p/io/nwb.py index 409f782a..c29b552b 100644 --- a/suite2p/io/nwb.py +++ b/suite2p/io/nwb.py @@ -1,35 +1,38 @@ import datetime +import gc import os -from natsort import natsorted +import time +from pathlib import Path import numpy as np -import time import scipy -import gc +from natsort import natsorted +from .. import run_s2p from ..detection.stats import roi_stats from . import utils from .. import run_s2p, default_ops try: - from pynwb import NWBFile + from pynwb import NWBHDF5IO, NWBFile from pynwb.base import Images from pynwb.image import GrayscaleImage - from pynwb.device import Device - from pynwb.ophys import OpticalChannel - from pynwb.ophys import TwoPhotonSeries - from pynwb.ophys import ImageSegmentation - from pynwb.ophys import RoiResponseSeries - from pynwb.ophys import Fluorescence - from pynwb import NWBHDF5IO + from pynwb.ophys import ( + Fluorescence, + ImageSegmentation, + OpticalChannel, + RoiResponseSeries, + TwoPhotonSeries, + ) + NWB = True except ModuleNotFoundError: NWB = False def nwb_to_binary(ops): - """ convert nwb file to binary (experimental) - + """convert nwb file to binary (experimental) + converts single plane single channel nwb file to binary for suite2p processing Parameters @@ -48,61 +51,62 @@ def nwb_to_binary(ops): 'frames_per_folder', 'nframes', 'meanImg', 'meanImg_chan2' """ - + # force 1 plane 1 chan for now - ops['nplanes'] = 1 - ops['nchannels'] = 1 + ops["nplanes"] = 1 + ops["nchannels"] = 1 # initialize ops with reg_file and raw_file paths, etc ops = utils.init_ops(ops)[0] - t0=time.time() - nplanes = ops['nplanes'] - nchannels = ops['nchannels'] - - batch_size = ops['batch_size'] - batch_size = int(nplanes*nchannels*np.ceil(batch_size/(nplanes*nchannels))) + t0 = time.time() + nplanes = ops["nplanes"] + nchannels = ops["nchannels"] + + batch_size = ops["batch_size"] + batch_size = int(nplanes * nchannels * np.ceil(batch_size / (nplanes * nchannels))) # open reg_file (and when available reg_file_chan2) - if 'keep_movie_raw' in ops and ops['keep_movie_raw']: - reg_file = open(ops['raw_file'], 'wb') - if nchannels>1: - reg_file_chan2 = open(ops['raw_file_chan2'], 'wb') + if "keep_movie_raw" in ops and ops["keep_movie_raw"]: + reg_file = open(ops["raw_file"], "wb") + if nchannels > 1: + reg_file_chan2 = open(ops["raw_file_chan2"], "wb") else: - reg_file = open(ops['reg_file'], 'wb') - if nchannels>1: - reg_file_chan2 = open(ops['reg_file_chan2'], 'wb') + reg_file = open(ops["reg_file"], "wb") + if nchannels > 1: + reg_file_chan2 = open(ops["reg_file_chan2"], "wb") nwb_driver = None - if ops.get('nwb_driver') and isinstance(nwb_driver, str): - nwb_driver = ops['nwb_driver'] + if ops.get("nwb_driver") and isinstance(nwb_driver, str): + nwb_driver = ops["nwb_driver"] - - with NWBHDF5IO(ops['nwb_file'], 'r', driver=nwb_driver) as fio: + with NWBHDF5IO(ops["nwb_file"], "r", driver=nwb_driver) as fio: nwbfile = fio.read() # get TwoPhotonSeries - if not ops.get('nwb_series'): + if not ops.get("nwb_series"): TwoPhotonSeries_names = [] for v in nwbfile.acquisition.values(): if isinstance(v, TwoPhotonSeries): TwoPhotonSeries_names.append(v.name) - if len(TwoPhotonSeries_names)==0: - raise ValueError('no TwoPhotonSeries in NWB file') - elif len(TwoPhotonSeries_names)>1: - raise Warning('more than one TwoPhotonSeries in NWB file, choosing first one') - ops['nwb_series'] = TwoPhotonSeries_names[0] - - series = nwbfile.acquisition[ops['nwb_series']] - series_shape = nwbfile.acquisition[ops['nwb_series']].data.shape - ops['nframes'] = series_shape[0] - ops['frames_per_file'] = np.array([ops['nframes']]) - ops['frames_per_folder'] = np.array([ops['nframes']]) - ops['meanImg'] = np.zeros((series_shape[1], series_shape[2]), np.float32) - for ik in np.arange(0, ops['nframes'], batch_size): - ikend = min(ik+batch_size, ops['nframes']) - im = series.data[ik : ikend] - + if len(TwoPhotonSeries_names) == 0: + raise ValueError("no TwoPhotonSeries in NWB file") + elif len(TwoPhotonSeries_names) > 1: + raise Warning( + "more than one TwoPhotonSeries in NWB file, choosing first one" + ) + ops["nwb_series"] = TwoPhotonSeries_names[0] + + series = nwbfile.acquisition[ops["nwb_series"]] + series_shape = nwbfile.acquisition[ops["nwb_series"]].data.shape + ops["nframes"] = series_shape[0] + ops["frames_per_file"] = np.array([ops["nframes"]]) + ops["frames_per_folder"] = np.array([ops["nframes"]]) + ops["meanImg"] = np.zeros((series_shape[1], series_shape[2]), np.float32) + for ik in np.arange(0, ops["nframes"], batch_size): + ikend = min(ik + batch_size, ops["nframes"]) + im = series.data[ik:ikend] + # check if uint16 if im.dtype.type == np.uint16: im = (im // 2).astype(np.int16) @@ -112,49 +116,58 @@ def nwb_to_binary(ops): im = im.astype(np.int16) reg_file.write(bytearray(im)) - ops['meanImg'] += im.astype(np.float32).sum(axis=0) + ops["meanImg"] += im.astype(np.float32).sum(axis=0) - if ikend%(batch_size*4)==0: - print('%d frames of binary, time %0.2f sec.'%(ikend,time.time()-t0)) + if ikend % (batch_size * 4) == 0: + print( + "%d frames of binary, time %0.2f sec." % (ikend, time.time() - t0) + ) gc.collect() # write ops files - do_registration = ops['do_registration'] - ops['Ly'],ops['Lx'] = ops['meanImg'].shape - ops['yrange'] = np.array([0,ops['Ly']]) - ops['xrange'] = np.array([0,ops['Lx']]) - ops['meanImg'] /= ops['nframes'] - if nchannels>1: - ops['meanImg_chan2'] /= ops['nframes'] + ops["Ly"], ops["Lx"] = ops["meanImg"].shape + ops["yrange"] = np.array([0, ops["Ly"]]) + ops["xrange"] = np.array([0, ops["Lx"]]) + ops["meanImg"] /= ops["nframes"] + if nchannels > 1: + ops["meanImg_chan2"] /= ops["nframes"] # close all binary files and write ops files - np.save(ops['ops_path'], ops) + np.save(ops["ops_path"], ops) reg_file.close() - if nchannels>1: + if nchannels > 1: reg_file_chan2.close() return ops def read_nwb(fpath): - """ read NWB file for use in the GUI """ - with NWBHDF5IO(fpath, 'r') as fio: + """read NWB file for use in the GUI""" + with NWBHDF5IO(fpath, "r") as fio: nwbfile = fio.read() - + # ROIs try: - rois = nwbfile.processing['ophys']['ImageSegmentation']['PlaneSegmentation']['pixel_mask'] + rois = nwbfile.processing["ophys"]["ImageSegmentation"][ + "PlaneSegmentation" + ]["pixel_mask"] multiplane = False - except: - rois = nwbfile.processing['ophys']['ImageSegmentation']['PlaneSegmentation']['voxel_mask'] + except Exception: + rois = nwbfile.processing["ophys"]["ImageSegmentation"][ + "PlaneSegmentation" + ]["voxel_mask"] multiplane = True stat = [] for n in range(len(rois)): if isinstance(rois[0], np.ndarray): stat.append( { - "ypix": np.array([rois[n][i][0].astype("int") for i in range(len(rois[n]))]), - "xpix": np.array([rois[n][i][1].astype("int") for i in range(len(rois[n]))]), - "lam": np.array([rois[n][i][-1] for i in range(len(rois[n]))]) + "ypix": np.array( + [rois[n][i][0].astype("int") for i in range(len(rois[n]))] + ), + "xpix": np.array( + [rois[n][i][1].astype("int") for i in range(len(rois[n]))] + ), + "lam": np.array([rois[n][i][-1] for i in range(len(rois[n]))]), } ) else: @@ -162,15 +175,17 @@ def read_nwb(fpath): { "ypix": rois[n]["x"].astype("int"), "xpix": rois[n]["y"].astype("int"), - "lam": rois[n]["weight"] + "lam": rois[n]["weight"], } ) if multiplane: stat[-1]['iplane'] = int(rois[n][0][-2]) ops = default_ops() - + if multiplane: - nplanes = np.max(np.array([stat[n]['iplane'] for n in range(len(stat))]))+1 + nplanes = ( + np.max(np.array([stat[n]["iplane"] for n in range(len(stat))])) + 1 + ) else: nplanes = 1 stat = np.array(stat) @@ -182,227 +197,343 @@ def read_nwb(fpath): bg_strs = ['meanImg', 'Vcorr', 'max_proj', 'meanImg_chan2'] ops['nchannels'] = 1 for bstr in bg_strs: - if bstr in nwbfile.processing['ophys']['Backgrounds_%d'%iplane].images: - ops[bstr] = np.array(nwbfile.processing['ophys']['Backgrounds_%d'%iplane][bstr].data) - if bstr=='meanImg_chan2': - ops['nchannels'] = 2 - ops['Ly'], ops['Lx'] = ops[bg_strs[0]].shape - ops['yrange'] = [0, ops['Ly']] - ops['xrange'] = [0, ops['Lx']] - ops['tau'] = 1.0 - ops['fs'] = nwbfile.acquisition['TwoPhotonSeries'].rate + if ( + bstr + in nwbfile.processing["ophys"]["Backgrounds_%d" % iplane].images + ): + ops[bstr] = np.array( + nwbfile.processing["ophys"]["Backgrounds_%d" % iplane][ + bstr + ].data + ) + if bstr == "meanImg_chan2": + ops["nchannels"] = 2 + ops["Ly"], ops["Lx"] = ops[bg_strs[0]].shape + ops["yrange"] = [0, ops["Ly"]] + ops["xrange"] = [0, ops["Lx"]] + ops["tau"] = 1.0 + ops["fs"] = nwbfile.acquisition["TwoPhotonSeries"].rate ops1.append(ops.copy()) stat = roi_stats(stat, ops['Ly'], ops['Lx'], ops['aspect'], ops['diameter']) - + # fluorescence - F = np.array(nwbfile.processing['ophys']['Fluorescence']['Fluorescence'].data) - Fneu = np.array(nwbfile.processing['ophys']['Neuropil']['Neuropil'].data) - spks = np.array(nwbfile.processing['ophys']['Deconvolved']['Deconvolved'].data) - dF = F - ops['neucoeff'] * Fneu + ophys = nwbfile.processing["ophys"] + + def get_fluo(name: str) -> np.ndarray: + """Extract Fluorescence data.""" + roi_response_series = ophys[name].roi_response_series + if name in roi_response_series.keys(): + fluo = ophys[name][name].data[:] + elif "plane0" in roi_response_series.keys(): + for key, value in roi_response_series.items(): + if key == "plane0": + fluo = value.data[:] + else: + fluo = np.concatenate((fluo, value.data[:]), axis=1) + fluo = np.transpose(fluo) + else: + raise AttributeError(f"Can't find {name} container in {fpath}") + return fluo + + F = get_fluo("Fluorescence") + Fneu = get_fluo("Neuropil") + spks = get_fluo("Deconvolved") + dF = F - ops["neucoeff"] * Fneu for n in range(len(stat)): - stat[n]['skew'] = scipy.stats.skew(dF[n]) + stat[n]["skew"] = scipy.stats.skew(dF[n]) # cell probabilities - iscell = [nwbfile.processing['ophys']['ImageSegmentation']['PlaneSegmentation']['iscell'][n] - for n in range(len(stat))] + iscell = [ + ophys["ImageSegmentation"]["PlaneSegmentation"]["iscell"][n] + for n in range(len(stat)) + ] iscell = np.array(iscell) - probcell = iscell[:,1] - iscell = iscell[:,0].astype('bool') + probcell = iscell[:, 1] + iscell = iscell[:, 0].astype("bool") redcell = np.zeros_like(iscell) probredcell = np.zeros_like(probcell) - + if multiplane: ops = ops1[0].copy() - Lx = ops['Lx'] - Ly = ops['Ly'] - nX = np.ceil(np.sqrt(ops['Ly'] * ops['Lx'] * len(ops1))/ops['Lx']) + Lx = ops["Lx"] + Ly = ops["Ly"] + nX = np.ceil(np.sqrt(ops["Ly"] * ops["Lx"] * len(ops1)) / ops["Lx"]) nX = int(nX) - nY = int(np.ceil(len(ops1)/nX)) for j in range(len(ops1)): - ops1[j]['dx'] = (j%nX) * Lx - ops1[j]['dy'] = int(j/nX) * Ly - - LY = int(np.amax(np.array([ops['Ly']+ops['dy'] for ops in ops1]))) - LX = int(np.amax(np.array([ops['Lx']+ops['dx'] for ops in ops1]))) + ops1[j]["dx"] = (j % nX) * Lx + ops1[j]["dy"] = int(j / nX) * Ly + + LY = int(np.amax(np.array([ops["Ly"] + ops["dy"] for ops in ops1]))) + LX = int(np.amax(np.array([ops["Lx"] + ops["dx"] for ops in ops1]))) meanImg = np.zeros((LY, LX)) max_proj = np.zeros((LY, LX)) - if ops['nchannels']>1: + if ops["nchannels"] > 1: meanImg_chan2 = np.zeros((LY, LX)) Vcorr = np.zeros((LY, LX)) - for k,ops in enumerate(ops1): - xrange = np.arange(ops['dx'],ops['dx']+ops['Lx']) - yrange = np.arange(ops['dy'],ops['dy']+ops['Ly']) - meanImg[np.ix_(yrange, xrange)] = ops['meanImg'] - Vcorr[np.ix_(yrange, xrange)] = ops['Vcorr'] - max_proj[np.ix_(yrange, xrange)] = ops['max_proj'] - if ops['nchannels']>1: - if 'meanImg_chan2' in ops: - meanImg_chan2[np.ix_(yrange, xrange)] = ops['meanImg_chan2'] - for j in np.nonzero(np.array([stat[n]['iplane']==k for n in range(len(stat))]))[0]: - stat[j]['xpix'] += ops['dx'] - stat[j]['ypix'] += ops['dy'] - stat[j]['med'][0] += ops['dy'] - stat[j]['med'][1] += ops['dx'] - ops['Vcorr'] = Vcorr - ops['max_proj'] = max_proj - ops['meanImg'] = meanImg - if 'meanImg_chan2' in ops: - ops['meanImg_chan2'] = meanImg_chan2 - ops['Ly'], ops['Lx'] = LY, LX - ops['yrange'] = [0, LY] - ops['xrange'] = [0, LX] + for k, ops in enumerate(ops1): + xrange = np.arange(ops["dx"], ops["dx"] + ops["Lx"]) + yrange = np.arange(ops["dy"], ops["dy"] + ops["Ly"]) + meanImg[np.ix_(yrange, xrange)] = ops["meanImg"] + Vcorr[np.ix_(yrange, xrange)] = ops["Vcorr"] + max_proj[np.ix_(yrange, xrange)] = ops["max_proj"] + if ops["nchannels"] > 1: + if "meanImg_chan2" in ops: + meanImg_chan2[np.ix_(yrange, xrange)] = ops["meanImg_chan2"] + for j in np.nonzero( + np.array([stat[n]["iplane"] == k for n in range(len(stat))]) + )[0]: + stat[j]["xpix"] += ops["dx"] + stat[j]["ypix"] += ops["dy"] + stat[j]["med"][0] += ops["dy"] + stat[j]["med"][1] += ops["dx"] + ops["Vcorr"] = Vcorr + ops["max_proj"] = max_proj + ops["meanImg"] = meanImg + if "meanImg_chan2" in ops: + ops["meanImg_chan2"] = meanImg_chan2 + ops["Ly"], ops["Lx"] = LY, LX + ops["yrange"] = [0, LY] + ops["xrange"] = [0, LX] return stat, ops, F, Fneu, spks, iscell, probcell, redcell, probredcell def save_nwb(save_folder): - """ convert folder with plane folders to NWB format """ - - plane_folders = natsorted([ f.path for f in os.scandir(save_folder) if f.is_dir() and f.name[:5]=='plane']) - ops1 = [np.load(os.path.join(f, 'ops.npy'), allow_pickle=True).item() for f in plane_folders] - - if NWB and not ops1[0]['mesoscan']: - if len(ops1)>1: + """convert folder with plane folders to NWB format""" + + plane_folders = natsorted( + [ + Path(f.path) + for f in os.scandir(save_folder) + if f.is_dir() and f.name[:5] == "plane" + ] + ) + ops1 = [ + np.load(f.joinpath("ops.npy"), allow_pickle=True).item() for f in plane_folders + ] + nchannels = min([ops["nchannels"] for ops in ops1]) + + if NWB and not ops1[0]["mesoscan"]: + if len(ops1) > 1: multiplane = True else: multiplane = False ops = ops1[0] + if "date_proc" in ops: + session_start_time = ops["date_proc"] + if not session_start_time.tzinfo: + session_start_time = session_start_time.astimezone() + else: + session_start_time = datetime.datetime.now().astimezone() - ### INITIALIZE NWB FILE + # INITIALIZE NWB FILE nwbfile = NWBFile( - session_description='suite2p_proc', - identifier=str(ops['data_path'][0]), - session_start_time=(ops['date_proc'] if 'date_proc' in ops - else datetime.datetime.now()) + session_description="suite2p_proc", + identifier=str(ops["data_path"][0]), + session_start_time=session_start_time, ) print(nwbfile) device = nwbfile.create_device( - name='Microscope', - description='My two-photon microscope', - manufacturer='The best microscope manufacturer' + name="Microscope", + description="My two-photon microscope", + manufacturer="The best microscope manufacturer", ) optical_channel = OpticalChannel( - name='OpticalChannel', - description='an optical channel', - emission_lambda=500. + name="OpticalChannel", + description="an optical channel", + emission_lambda=500.0, ) imaging_plane = nwbfile.create_imaging_plane( - name='ImagingPlane', + name="ImagingPlane", optical_channel=optical_channel, - imaging_rate=ops['fs'], - description='standard', + imaging_rate=ops["fs"], + description="standard", device=device, - excitation_lambda=600., - indicator='GCaMP', - location='V1', - grid_spacing=([2.0,2.0,30.0] if multiplane else [2.0,2.0]), - grid_spacing_unit='microns' + excitation_lambda=600.0, + indicator="GCaMP", + location="V1", + grid_spacing=([2.0, 2.0, 30.0] if multiplane else [2.0, 2.0]), + grid_spacing_unit="microns", ) # link to external data image_series = TwoPhotonSeries( - name='TwoPhotonSeries', - dimension=[ops['Ly'], ops['Lx']], - external_file=(ops['filelist'] if 'filelist' in ops else ['']), + name="TwoPhotonSeries", + dimension=[ops["Ly"], ops["Lx"]], + external_file=(ops["filelist"] if "filelist" in ops else [""]), imaging_plane=imaging_plane, - starting_frame=[0], - format='external', - starting_time=0.0, - rate=ops['fs'] * ops['nplanes'] + starting_frame=[0], + format="external", + starting_time=0.0, + rate=ops["fs"] * ops["nplanes"], ) nwbfile.add_acquisition(image_series) # processing img_seg = ImageSegmentation() ps = img_seg.create_plane_segmentation( - name='PlaneSegmentation', - description='suite2p output', + name="PlaneSegmentation", + description="suite2p output", imaging_plane=imaging_plane, - reference_images=image_series + reference_images=image_series, ) ophys_module = nwbfile.create_processing_module( - name='ophys', - description='optical physiology processed data' + name="ophys", description="optical physiology processed data" ) ophys_module.add(img_seg) - - file_strs = ['F.npy', 'Fneu.npy', 'spks.npy'] - traces = [] - ncells_all = 0 - Nfr = np.array([ops['nframes'] for ops in ops1]).max() + + file_strs = ["F.npy", "Fneu.npy", "spks.npy"] + file_strs_chan2 = ["F_chan2.npy", "Fneu_chan2.npy"] + traces, traces_chan2 = [], [] + ncells = np.zeros(len(ops1), dtype=np.int_) + Nfr = np.array([ops["nframes"] for ops in ops1]).max() for iplane, ops in enumerate(ops1): - if iplane==0: - iscell = np.load(os.path.join(ops['save_path'], 'iscell.npy')) + if iplane == 0: + iscell = np.load(os.path.join(ops["save_path"], "iscell.npy")) for fstr in file_strs: - traces.append(np.load(os.path.join(ops['save_path'], fstr))) + traces.append(np.load(os.path.join(ops["save_path"], fstr))) + if nchannels > 1: + for fstr in file_strs_chan2: + traces_chan2.append( + np.load(plane_folders[iplane].joinpath(fstr)) + ) + PlaneCellsIdx = iplane * np.ones(len(iscell)) else: - iscell = np.append(iscell, np.load(os.path.join(ops['save_path'], 'iscell.npy')), axis=0) - for i,fstr in enumerate(file_strs): - trace = np.load(os.path.join(ops['save_path'], fstr)) + iscell = np.append( + iscell, + np.load(os.path.join(ops["save_path"], "iscell.npy")), + axis=0, + ) + for i, fstr in enumerate(file_strs): + trace = np.load(os.path.join(ops["save_path"], fstr)) if trace.shape[1] < Nfr: - fcat = np.zeros((trace.shape[0],Nfr-trace.shape[1]), 'float32') - trace = np.concatenate((trace, fcat), axis=1) - traces[i] = np.append(traces[i], trace, axis=0) - - stat = np.load(os.path.join(ops['save_path'], 'stat.npy'), allow_pickle=True) - ncells = len(stat) - for n in range(ncells): + fcat = np.zeros( + (trace.shape[0], Nfr - trace.shape[1]), "float32" + ) + trace = np.concatenate((trace, fcat), axis=1) + traces[i] = np.append(traces[i], trace, axis=0) + if nchannels > 1: + for i, fstr in enumerate(file_strs_chan2): + traces_chan2[i] = np.append( + traces_chan2[i], + np.load(plane_folders[iplane].joinpath(fstr)), + axis=0, + ) + PlaneCellsIdx = np.append( + PlaneCellsIdx, iplane * np.ones(len(iscell) - len(PlaneCellsIdx)) + ) + + stat = np.load( + os.path.join(ops["save_path"], "stat.npy"), allow_pickle=True + ) + ncells[iplane] = len(stat) + for n in range(ncells[iplane]): if multiplane: - pixel_mask = np.array([stat[n]['ypix'], stat[n]['xpix'], - iplane*np.ones(stat[n]['npix']), - stat[n]['lam']]) + pixel_mask = np.array( + [ + stat[n]["ypix"], + stat[n]["xpix"], + iplane * np.ones(stat[n]["npix"]), + stat[n]["lam"], + ] + ) ps.add_roi(voxel_mask=pixel_mask.T) else: - pixel_mask = np.array([stat[n]['ypix'], stat[n]['xpix'], - stat[n]['lam']]) + pixel_mask = np.array( + [stat[n]["ypix"], stat[n]["xpix"], stat[n]["lam"]] + ) ps.add_roi(pixel_mask=pixel_mask.T) - ncells_all+=ncells - ps.add_column('iscell', 'two columns - iscell & probcell', iscell) - - rt_region = ps.create_roi_table_region( - region=list(np.arange(0, ncells_all)), - description='all ROIs' - ) + ps.add_column("iscell", "two columns - iscell & probcell", iscell) + + rt_region = [] + for iplane, ops in enumerate(ops1): + if iplane == 0: + rt_region.append( + ps.create_roi_table_region( + region=list( + np.arange(0, ncells[iplane]), + ), + description=f"ROIs for plane{int(iplane)}", + ) + ) + else: + rt_region.append( + ps.create_roi_table_region( + region=list( + np.arange( + np.sum(ncells[:iplane]), + ncells[iplane] + np.sum(ncells[:iplane]), + ) + ), + description=f"ROIs for plane{int(iplane)}", + ) + ) # FLUORESCENCE (all are required) - file_strs = ['F.npy', 'Fneu.npy', 'spks.npy'] - name_strs = ['Fluorescence', 'Neuropil', 'Deconvolved'] - - for i, (fstr,nstr) in enumerate(zip(file_strs, name_strs)): - roi_resp_series = RoiResponseSeries( - name=nstr, - data=traces[i], - rois=rt_region, - unit='lumens', - rate=ops['fs'] - ) - fl = Fluorescence(roi_response_series=roi_resp_series, name=nstr) + name_strs = ["Fluorescence", "Neuropil", "Deconvolved"] + name_strs_chan2 = ["Fluorescence_chan2", "Neuropil_chan2"] + + for i, (fstr, nstr) in enumerate(zip(file_strs, name_strs)): + for iplane, ops in enumerate(ops1): + roi_resp_series = RoiResponseSeries( + name=f"plane{int(iplane)}", + data=np.transpose(traces[i][PlaneCellsIdx == iplane]), + rois=rt_region[iplane], + unit="lumens", + rate=ops["fs"], + ) + if iplane == 0: + fl = Fluorescence(roi_response_series=roi_resp_series, name=nstr) + else: + fl.add_roi_response_series(roi_response_series=roi_resp_series) ophys_module.add(fl) + if nchannels > 1: + for i, (fstr, nstr) in enumerate(zip(file_strs_chan2, name_strs_chan2)): + for iplane, ops in enumerate(ops1): + roi_resp_series = RoiResponseSeries( + name=f"plane{int(iplane)}", + data=np.transpose(traces_chan2[i][PlaneCellsIdx == iplane]), + rois=rt_region[iplane], + unit="lumens", + rate=ops["fs"], + ) + + if iplane == 0: + fl = Fluorescence( + roi_response_series=roi_resp_series, name=nstr + ) + else: + fl.add_roi_response_series(roi_response_series=roi_resp_series) + + ophys_module.add(fl) + # BACKGROUNDS # (meanImg, Vcorr and max_proj are REQUIRED) - bg_strs = ['meanImg', 'Vcorr', 'max_proj', 'meanImg_chan2'] - nplanes = ops['nplanes'] + bg_strs = ["meanImg", "Vcorr", "max_proj", "meanImg_chan2"] + nplanes = ops["nplanes"] for iplane in range(nplanes): - images = Images('Backgrounds_%d'%iplane) + images = Images("Backgrounds_%d" % iplane) for bstr in bg_strs: if bstr in ops: - if bstr=='Vcorr' or bstr=='max_proj': - img = np.zeros((ops['Ly'], ops['Lx']), np.float32) - img[ops['yrange'][0]:ops['yrange'][-1], - ops['xrange'][0]:ops['xrange'][-1]] = ops[bstr] + if bstr == "Vcorr" or bstr == "max_proj": + img = np.zeros((ops["Ly"], ops["Lx"]), np.float32) + img[ + ops["yrange"][0] : ops["yrange"][-1], + ops["xrange"][0] : ops["xrange"][-1], + ] = ops[bstr] else: img = ops[bstr] images.add_image(GrayscaleImage(name=bstr, data=img)) - + ophys_module.add(images) - with NWBHDF5IO(os.path.join(save_folder, 'ophys.nwb'), 'w') as fio: + with NWBHDF5IO(os.path.join(save_folder, "ophys.nwb"), "w") as fio: fio.write(nwbfile) else: print('pip install pynwb OR don"t use mesoscope recording') diff --git a/suite2p/run_s2p.py b/suite2p/run_s2p.py index ff22e7bd..0e48b3d1 100644 --- a/suite2p/run_s2p.py +++ b/suite2p/run_s2p.py @@ -226,7 +226,7 @@ def run_plane(ops, ops_path=None, stat=None): """ ops = {**default_ops(), **ops} - ops['date_proc'] = datetime.now() + ops['date_proc'] = datetime.now().astimezone() # for running on server or on moved files, specify ops_path if ops_path is not None: diff --git a/tests/test_io.py b/tests/test_io.py index be52de6c..afde72ac 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -1,6 +1,8 @@ """ Tests for the Suite2p IO module """ +import pathlib +import platform import re from pathlib import Path @@ -8,18 +10,19 @@ import pytest from natsort import natsorted from pynwb import NWBHDF5IO - from suite2p import io -from suite2p.io.nwb import save_nwb +from suite2p.io.nwb import read_nwb, save_nwb from suite2p.io.utils import get_suite2p_path @pytest.fixture() def binfile1500(test_ops): - test_ops['tiff_list'] = ['input_1500.tif'] + test_ops["tiff_list"] = ["input_1500.tif"] op = io.tiff_to_binary(test_ops) - bin_filename = str(Path(op['save_path0']).joinpath('suite2p/plane0/data.bin')) - with io.BinaryFile(Ly=op['Ly'], Lx=op['Lx'], read_filename=bin_filename) as bin_file: + bin_filename = str(Path(op["save_path0"]).joinpath("suite2p/plane0/data.bin")) + with io.BinaryFile( + Ly=op["Ly"], Lx=op["Lx"], read_filename=bin_filename + ) as bin_file: yield bin_file @@ -30,9 +33,14 @@ def replace_ops_save_path_with_local_path(request): by its local path version """ - # Get the `data_folder` variable from the running test name - data_folder = re.search(r"\[(.*)\]", request.node.name).group(1) - save_folder = Path("data").joinpath("test_inputs", data_folder, "suite2p") + # Workaround to load pickled NPY files on Windows containing + # `PosixPath` objects + if platform.system() == 'Windows': + pathlib.PosixPath = pathlib.WindowsPath + + # Get the `data_folder` variable from the running test name + data_folder = re.search(r"\[(.*?)(-.*?)?\]", request.node.name).group(1) + save_folder = Path("data").joinpath("test_outputs", data_folder, "suite2p") save_path = {} plane_folders = [ @@ -48,17 +56,27 @@ def replace_ops_save_path_with_local_path(request): ops1.item(0)["save_path"] = str(plane_dir.absolute()) np.save(plane_dir.joinpath("ops.npy"), ops1) - # Concatenate iscell arrays from the NumPy files - if plane_idx == 0: - iscell_npy = np.load(plane_dir.joinpath("iscell.npy"), allow_pickle=True) - else: - iscell_npy = np.append( - iscell_npy, - np.load(plane_dir.joinpath("iscell.npy"), allow_pickle=True), - axis=0, - ) - - yield save_folder, iscell_npy + def concat_npy(name: str) -> np.ndarray: + """Concatenate arrays from NUmPy files.""" + for plane_idx, plane_dir in enumerate(plane_folders): + if plane_idx == 0: + array_npy = np.load( + plane_dir.joinpath(f"{name}.npy"), allow_pickle=True + ) + else: + array_npy = np.append( + array_npy, + np.load(plane_dir.joinpath(f"{name}.npy"), allow_pickle=True), + axis=0, + ) + return array_npy + + iscell_npy = concat_npy("iscell") + F_npy = concat_npy("F") + Fneu_npy = concat_npy("Fneu") + spks_npy = concat_npy("spks") + + yield save_folder, iscell_npy, F_npy, Fneu_npy, spks_npy # Teardown the fixture for plane_dir in plane_folders: @@ -67,15 +85,21 @@ def replace_ops_save_path_with_local_path(request): ops1.item(0)["save_path"] = save_path[plane_dir] np.save(plane_dir.joinpath("ops.npy"), ops1) + def test_h5_to_binary_produces_nonnegative_output_data(test_ops): - test_ops['h5py'] = Path(test_ops['data_path'][0]).joinpath('input.h5') - test_ops['nplanes'] = 3 - test_ops['nchannels'] = 2 - test_ops['data_path'] = [] + test_ops["h5py"] = Path(test_ops["data_path"][0]).joinpath("input.h5") + test_ops["nplanes"] = 3 + test_ops["nchannels"] = 2 + test_ops["data_path"] = [] op = io.h5py_to_binary(test_ops) - output_data = io.BinaryFile(read_filename=Path(op['save_path0'], 'suite2p/plane0/data.bin'), Ly=op['Ly'], Lx=op['Lx']).data + output_data = io.BinaryFile( + read_filename=Path(op["save_path0"], "suite2p/plane0/data.bin"), + Ly=op["Ly"], + Lx=op["Lx"], + ).data assert np.all(output_data >= 0) + def test_that_bin_movie_without_badframes_results_in_a_same_size_array(binfile1500): mov = binfile1500.bin_movie(bin_size=1) assert mov.shape == (1500, binfile1500.Ly, binfile1500.Lx) @@ -88,7 +112,9 @@ def test_that_bin_movie_with_badframes_results_in_a_smaller_array(binfile1500): mov = binfile1500.bin_movie(bin_size=1, bad_frames=bad_frames, reject_threshold=0) assert len(mov) < binfile1500.n_frames, "bin_movie didn't produce a smaller array." - assert len(mov) == len(bad_frames) - sum(bad_frames), "bin_movie didn't produce the right size array." + assert len(mov) == len(bad_frames) - sum( + bad_frames + ), "bin_movie didn't produce the right size array." def test_that_binaryfile_data_is_repeatable(binfile1500): @@ -104,39 +130,73 @@ def test_that_binaryfile_data_is_repeatable(binfile1500): @pytest.mark.parametrize( "data_folder", [ - ("1plane1chan"), ("1plane1chan1500"), - # ("1plane2chan"), - # ("1plane2chan-scanimage"), - # TODO: Make the test work with the commented folders above - # `np.load("ops.npy")` with `allow_pickle=True` currently fails with: - # NotImplementedError: cannot instantiate 'WindowsPath' on your system - ("2plane2chan"), ("2plane2chan1500"), + ("bruker"), ], ) -def temp_test_save_nwb(replace_ops_save_path_with_local_path, data_folder): - save_folder, iscell_npy = replace_ops_save_path_with_local_path - +def test_nwb_round_trip(replace_ops_save_path_with_local_path, data_folder): + + # Get expected data already saved as NumPy files + ( + save_folder, + expected_iscell, + expected_F, + expected_Fneu, + expected_spks, + ) = replace_ops_save_path_with_local_path + + # Save as NWB file save_nwb(save_folder) - with NWBHDF5IO(str(save_folder.joinpath("ophys.nwb")), "r") as io: + # Check (some of) the structure of the NWB file saved + nwb_path = save_folder.joinpath("ophys.nwb") + assert nwb_path.exists() + with NWBHDF5IO(str(nwb_path), "r") as io: read_nwbfile = io.read() assert read_nwbfile.processing - assert read_nwbfile.processing["ophys"].data_interfaces["Deconvolved"] - assert read_nwbfile.processing["ophys"].data_interfaces["Fluorescence"] - assert read_nwbfile.processing["ophys"].data_interfaces["Neuropil"] + ophys = read_nwbfile.processing["ophys"] + assert ophys.data_interfaces["Deconvolved"] + Fluorescence = ophys.data_interfaces["Fluorescence"] + assert Fluorescence + if "2plane" in data_folder: + assert Fluorescence["plane0"] + assert Fluorescence["plane1"] + assert ophys.data_interfaces["Neuropil"] + if "2chan" in data_folder: + Fluorescence_chan2 = ophys.data_interfaces["Fluorescence_chan2"] + assert Fluorescence_chan2 + assert ophys.data_interfaces["Neuropil_chan2"] iscell_nwb = ( - read_nwbfile.processing["ophys"] - .data_interfaces["ImageSegmentation"] + ophys.data_interfaces["ImageSegmentation"] .plane_segmentations["PlaneSegmentation"] .columns[2] .data[:] ) - np.testing.assert_array_equal(iscell_nwb, iscell_npy) + np.testing.assert_array_equal(iscell_nwb, expected_iscell) + + # Extract Suite2p info from NWB file + stat, ops, F, Fneu, spks, iscell, probcell, redcell, probredcell = read_nwb( + nwb_path + ) + + # Check we get the same data back as the original data + # after saving to NWB + reading + np.testing.assert_array_equal(F, expected_F) + np.testing.assert_array_equal(Fneu, expected_Fneu) + np.testing.assert_array_equal(spks, expected_spks) + np.testing.assert_array_equal( + np.transpose(np.array([iscell, probcell])), expected_iscell + ) + # TODO: assert round trip for `stat` and `ops` + # Probably need to recreate the data files as some fields are missing in the dict + # expected_stat = np.load(save_folder.joinpath("plane0", "stat.npy"), allow_pickle=True) + # expected_ops = np.load(save_folder.joinpath("plane0", "ops.npy"), allow_pickle=True) + # np.testing.assert_equal(stat, expected_stat) + # np.testing.assert_equal(ops, expected_ops) # Remove NWB file - save_folder.joinpath("ophys.nwb").unlink() + nwb_path.unlink() @pytest.mark.parametrize(